# Federated Domain Adaptation for ASR with Full Self-Supervision

###### Abstract

Cross-device federated learning (FL) protects user privacy by collaboratively
training a model on user devices, therefore eliminating the need for
collecting, storing, and manually labeling user data. While important topics
such as the FL training algorithm, non-IID-ness, and Differential Privacy have
been well studied in the literature, this paper focuses on two challenges of
practical importance for improving on-device ASR: the lack of ground-truth
transcriptions and the scarcity of compute resource and network bandwidth on
edge devices. First, we propose a FL system for on-device ASR domain
adaptation with full self-supervision, which uses self-labeling together with
data augmentation and filtering techniques. The system can improve a strong
Emformer-Transducer based ASR model pretrained on out-of-domain data, using
in-domain audio without any ground-truth transcriptions. Second, to reduce the
training cost, we propose a self-restricted RNN Transducer (SR-RNN-T) loss, a
variant of alignment-restricted RNN-T that uses Viterbi alignments from self-
supervision. To further reduce the compute and network cost, we systematically
explore adapting only a subset of weights in the Emformer-Transducer. Our best
training recipe achieves a $12.9\%$ relative WER reduction over the strong
out-of-domain baseline, which equals $70\%$ of the reduction achievable with
full human supervision and centralized training.

Index Terms: ASR, Adaptation, Federated Learning 111This paper was submitted
to Interspeech 2022.

##  1 Introduction

With the increasing adoption of deep learning in artificial intelligence
applications, data privacy is of growing concern [1, 2]. Automatic speech
recognition (ASR) is a particularly sensitive use case due to the personalized
nature of speech [3, 4]. Aiming to protect user privacy, cross-device
federated learning (FL) has been proposed and applied to a variety of tasks
[5, 6, 7, 8]. FL is a distributed training paradigm that allows a loose
federation of trainers to collectively improve a shared model [9, 10, 11, 12].

Cross-device FL for ASR faces several unique challenges. These include the
lack of ground truth transcriptions; the high costs of computation and network
communication; the non independent and identical distribution of training data
(non-IID-ness); and the difficulty of providing privacy guarantees. Several
recent works have considered cross-device FL for ASR applications [13, 14, 15,
16, 17, 18]. In particular, weighted model averaging [13, 14] and federated
variation noise [16] have been proposed to address the challenge of training
on non-IID data, while differential privacy has been proposed to provide
formal privacy guarantees [19, 20]. In this paper, however, we address the
other two practical challenges for on-device FL: the lack of ground truth
transcriptions, and the scarcity of compute resources and network bandwidth.

Figure 1: A FL system for on-device ASR domain adaptation with self-
supervision. Each device runs a local trainer alongside the ASR service.

A common assumption in existing literature on FL for ASR is the availability
of ground-truth transcriptions [13, 14, 15, 16, 17, 18]. However, this
assumption is not practical. Users of on-device ASR have neither the incentive
to manually transcribe their audio, nor the desire to frequently edit
inaccurate automatic transcriptions. Assuming no labeled audio data, [11]
proposes to use a large pretrained teacher model to label target domain audio
on cloud servers; other works, not focused on ASR, have studied semi-
supervised FL for keyword spotting, audio classification [21], and image
classification [22, 23]. However, to the best of our knowledge, there has been
limited literature that explored on-device adaptation for ASR with self-
supervision. To this end, we propose an FL system for on-device model domain
adaptation with self-supervision that runs alongside an ASR application.

On a high level, user audio and respective automatic transcriptions (pseudo
labels) from an on-device ASR model are stored temporarily on the device in an
encrypted queue. Then, a trainer that runs alongside the ASR process consumes
minibatches of examples from the queue and performs local model updates. Model
updates are periodically sent to a central FL server for aggregation. We
develop data augmentation and filtering techniques to improve robustness of
the self-supervised learning. Our offline experiments on two diverse FL
adaptation scenarios confirm the effectiveness of the proposed system.

Another important challenge in cross-device FL is the cost of training and
synchronization with the central FL server [10], since edge devices typically
have limited compute, battery, and network resources [9]. In this work, we
present a combination of techniques to improve an Emformer-Transducer based
ASR model [24, 25], which is trained with the RNN-T loss [26]. First, to
reduce the compute cost during on-device adaptation, we propose a self-
restricted (SR)-RNN-T loss, which extends the alignment-restricted (AR)-RNN-T
[27] by replacing the external alignment with a Viterbi alignment (similar to
[28]) from the model itself. Secondly, to further reduce the network cost
during synchronization, we consider only adapting a subset of model weights.
Inspired by a recent work on speaker adaptation [29], we study which subset of
weights in the ASR model would give the best domain-adaptation performance. We
find that adapting the key/value matrices in the self-attention module can
match (or surpass) full model adaptation, even though those matrices only
account for $5.7\%$ model weights. A related approach called partial variation
training adapts different subset of weights on different devices [30].

Combining all techniques, we achieve 12.9% relative WER reduction for on-
device adaptation — 70% of the reduction achievable with human supervision and
centralized training.

##  2 Methodology

Input: pretrained model $w_{0}$, # of rounds $T$, # of devices $N$, # of
updates per round $K$

Output: finetuned model $w_{T}$

1

2

3Procedure  _CentralFLServer(__)_ :

4 for  _t $\in$ 1 …T_ do

5 for  _i $\in$ 1 …N_ do

6
${\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\texttt{SendModel}}(w_{t-1},i,t)$

7
$\Delta_{t}^{(i)}\leftarrow{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\texttt{ReceiveUpdate}}(i,t)$

8 $\Delta_{t}\leftarrow\textstyle\frac{1}{N}\sum_{i=1}^{N}\Delta_{t}^{(i)}$

9 $w_{t}\leftarrow w_{t-1}+\beta(w_{t-1}-w_{t-2})+\Delta_{t}$

10

11

12Procedure  _Device(_$i$_)_:

13 $Q_{i}\leftarrow\\{\\}$

14 $\texttt{ASRService}(Q_{i})$

15 $\texttt{Trainer}(i,Q_{i})$

16

17Procedure  _ASRService(_$Q_{i}$_)_:

18 loop  __

19 $\bm{x}\leftarrow\texttt{GetRequest}()$

20 $\hat{\bm{y}}\leftarrow\texttt{EncodeDecode}(w_{0},\bm{x})$

21 $\texttt{PushToQueue}(Q_{i},(\bm{x},\hat{\bm{y}}))$

22

23

24Procedure  _Trainer(_$i$ , $Q_{i}$_)_:

25 for  _t $\in$ 1 …T_ do

26
$w_{t-1}\leftarrow{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\texttt{ReceiveModel}}(i,t)$

27 $w_{t,0}^{(i)}\leftarrow w_{t-1}$

28 for  _k $\in$ 1 …K_ do

29 $B_{k}\leftarrow\texttt{PopBatchFromQueue}(Q_{i})$

30 $\bar{B}_{k}\leftarrow\texttt{FilterAndAugment}(B_{k})$

31 $w_{t,k}^{(i)}\leftarrow\texttt{ModelUpdate}(w_{t,k-1}^{(i)},\bar{B}_{k})$

32 $\Delta_{t}^{(i)}\leftarrow w_{t,K}^{(i)}-w_{t,0}^{(i)}$

33
${\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\texttt{SendUpdate}}(\Delta_{t}^{(i)},i,t)$

34

Algorithm 1 FL system for ASR domain adaptation.

###  2.1 Federated Learning Alongside the ASR Process

This section presents the design of the proposed FL system. An initial ASR
model is pretrained on large amounts of supervised data from a broad range of
sources, using ground truth transcriptions or high-quality machine generated
transcriptions from a large offline teacher model. Then, the FL system adapts
the pretrained ASR model to target domain data using pseudo labels generated
by the pretrained model itself, via a hierarchical SGD method known in the FL
community as Federated Averaging with Momentum, abbreviated as FedAvgM [31],
although the same method was first developed in the speech community under the
name blockwise model-update filtering, or BMUF [32].

Figure 1 illustrates the cross-device FL system, comprising a central FL
server and many ASR devices, implementing Algorithm 1. Each device $i$
executes an ASR process that runs the pretrained model $w_{0}$, which is fixed
throughout the FL process. A device-specific trainer manages a second copy of
the ASR model $w_{t,k}^{(i)}$ that is updated via FL. Each time the device
receives an ASR request, the pretrained model is used to transcribe audio
$\bm{x}$ to generate transcription $\hat{\bm{y}}$, which we refer as pseudo
labels (this differs from [11] which uses a large offline teacher model for
this purpose). The $(\bm{x},\hat{\bm{y}})$ pairs are added to a queue for
later consumption by the trainer. In each local model update, the trainer
samples a minibatch from the queue and computes model gradients with the SR-
RNN-T loss. We find that it is important for adaptation performance to use
confidence filtering and data augmentation. Finally, after $K$ updates, the
difference $\Delta^{(i)}$ is sent back to the central FL server for
aggregation.

###  2.2 Efficient RNN-T Adaptation with Self-Alignment

Figure 2: Illustration of the proposed SR-RNN-T loss. In contrast with AR-
RNN-T, which relies on an external alignment (red), SR-RNN-T uses the Viterbi
forced-alignment (cyan).

We use an Emformer-Transducer [26, 33] as our ASR model, pretrained with AR-
RNN-T loss and centralized distributed data-parallel (DDP) training. Let
$\bm{x}$ denote the audio features over $T$ time steps and $\bm{y}$ denote the
reference text over $U$ target tokens. First, the encoder maps the audio
features $\bm{x}$ into input embeddings $\bm{h}$, and the predictor maps the
reference tokens $\bm{y}$ into output embeddings $\bm{g}$. Then, a joiner
network combines every pair of $h_{t},g_{u}\in\bm{h},\bm{g}$ to compute the
probability distribution $\textrm{Pr}(\cdot|h_{t},g_{u})$ over all possible
output tokens including blank. During training, the RNN-T loss optimizes the
sum of the probabilities over all possible audio-text alignments. During
inference, a beam search algorithm is used to output the most probable
sequence, which we refer to as $\hat{\bm{y}}$.

In FL scenario, $\bm{y}$ is not available. Instead, we use the most probable
sequence $\hat{\bm{y}}$ from inference as target labels. A straight-forward
implementation of the RNN-T loss would directly use $\bm{x}$ and
$\hat{\bm{y}}$ and optimize over all possible alignments. However, this
approach is computationally intensive. As a potential alternative, [27] showed
that the AR-RNN-T loss can greatly reduce memory and compute by restricting
the forward-backward algorithm to a band around a ground truth alignment.
However, AR-RNN-T relies on external ground-truth alignments generated offline
with a separate hybrid model. Motivated by some recent studies [28], we
propose to use Viterbi forced-alignment from the model itself in place of an
external alignment. Particularly, the forced-alignments are computed on the
fly over the RNN-T grid using the backtracing (see Fig. 2). We refer the
resulting loss function as self-restricted (SR)-RNN-T loss. In our design of
the FL system, the Viterbi forced-alignment is further approximated with the
most probable alignment encountered during beam search decoding after pruning.

###  2.3 Data Filtering and Augmentation

One main concern with self-supervision is the quality of the pseudo labels
$\hat{\bm{y}}$. Incorrect pseudo labels carry the risk of reinforcing errors.
We find an effective mitigation to be data filtering, where all utterances
with $\log{\textrm{Pr}(\hat{\bm{y}})}<\theta$ are discarded. The threshold
$\theta$ is a tunable hyper-parameter, which is empirically determined on a
small development dataset. We find that data filtering is especially important
if adaptation domain data is noisy (see section 3.4).

On the other hand, when the target domain audios are clean, we find the semi-
supervised self-training can be further boosted with data augmentation (speed
perturbation and additive noise) while keeping the pseudo labels unchanged
(see section 3.3).

###  2.4 Finetuning Only a Subset of Model Weights

To reduce the cost of communication (and computation), we consider only
adapting and communicating a subset of the model weights. Inspired by a recent
study [29], we systemically studied which subset of weights in the Emformer-
Transducer model provides the best trade-off between domain-adaptation
performance vs. communication cost. We explore adaptation of encoder,
predictor and joiner separately. Within the encoder, we carefully analyze
adapting only attention or attention key/value matrices. We find adapting the
key/value matrices is very effective for the target domain, and freezing the
rest of weights serves as a regularization technique, which prevents
catastrophic forgetting on the source domain.

##  3 Experiments

Table 1: Different experiment setups for simulating diverse set of FL
experiments. The 1st row in each setup denotes the source domain data used for
model pretraining, and 2nd row denotes the target domain data used for FL
adaptation.

Setup | Dataset | Train (hours) | Eval (utterances)
---|---|---|---
Device | Mixture | 1.5M | 3.6K
Edge | 9.7K | 28.6K
Video | Clean | 17.4K | 1.2K
Noisy | 1.2K | 1.5K

###  3.1 Data Setup

To simulate realistic federated learning in practical scenarios, we consider
two different experiment setups. Each setup takes a converged model pretrained
on a source domain and performs FL adaptation to a target domain. The adapted
model is tested on the evaluation data from both the source and target
domains. Table 1 gives more information on the data used in each setup.

Device: This setup closely resembles a cross-device FL use case. For
pretraining, we use a mixture of data sources, including public videos, audios
collected by paid third-party vendors on portal, and synthetic audios
generated by text-to-speech. This dataset contains 1.5 Million hours of
audios, which are transcribed either through human annotators or a large
teacher model (see [34] for details). For adaptation, we use assistant and
dictation audios recorded on edge devices, which are collected by external
vendors under clean acoustic conditions. The source domain evaluation data is
a separately prepared conversation dataset, while the target domain evaluation
data is a held out subset of the on-device recordings.

Video: For pretraining, we use 17.4K hours of long form videos (averaged 45
sec) with relative clean audios. For adaptation, we use short form videos
(averaged 12 sec) that might be recorded in noisy environments. Both datasets
are de-identified before transcription by human annotators. A subset of the
audios from both the clean and noisy videos are held out for evaluation.

###  3.2 Modeling

Our ASR model is an Emformer-Transducer [33]. Features are 80 mel-spaced log
filter-bank channels with a sliding window of 25ms at 10ms frame shift,
regularized with SpecAugment [35]. A time reduction layer stacks features with
a stride of 6, and the resulting 480-dimensional vectors form the input to the
encoder. The predictor is a 2-layer LSTM model with a hidden dimension of 256.
Both the encoder and the predictor outputs are projected to 768-dimension
tensors and combined in the RNN-T joiner, which consists of a ReLU followed by
a softmax layer with a sentence-piece vocabulary size of 4096.

For each task, we first pretrain the ASR model on the source domain until
convergence, using an AR-RNN-T loss with buffer size $b_{l}=5,b_{r}=25$ and an
Adam optimizer with learning rate $10^{-3}$ and tri-stage learning rate
scheduling. We then finetune the ASR model to the target domain using 32
parallel workers. Each worker uses an Adam optimizer with a fixed learning
rate of $10^{-4}$ for local model updates. The local models on different
workers are synchronized every 20 local updates using FedAvgM/BMUF with a
block momentum of 0.8.

###  3.3 Device Adaptation

Figure 3:  Convergence of the Device adaptation experiments. The combination
of FedAvgM/BMUF, pseudo-labels, SR-RNN-T and KeyValue adaptation not only
achieves a low WER on the target domain, but also prevents catastrophic
forgetting on the source domain.

The Device adaptation setup adapts from the source domain “Mixture” to the
target domain “Edge”. The baseline pretrained model achieves a lower WER on
the target domain (4.96) than on the source domain (6.07), since the latter
captures diverse acoustic conditions while the former is a commissioned data
collection that is mostly clean.

Table 2: Summary of Device adaptation experiments. Model | Target | Source
---|---|---
Pretrained | 4.96 | 6.07
1 (DDP, True Labels, AR-RNN-T) | 4.05 | 6.12
2 (1 \+ FedAvgM/BMUF) | 4.09 | 6.16
3 (2 \+ Pseudo Labels) | 4.35 | 6.17
4 (3 \+ SR-RNN-T) | 4.32 | 6.08
5 (4 \+ KeyValue Adaptation) | 4.32 | 5.95

Lower bound: We first consider supervised adaptation with centralized DDP
training, using the ground-truth labels and the AR-RNN-T loss. The adapted
model’s WER on the source and target domains across different checkpoints are
plotted in figure 3 1. In particular, the best checkpoint (see table 2 1)
reduces the WER $4.96\rightarrow 4.05$ on the target domain, which serves as a
lower bound for semi-supervised FL adaptations. We also observe a slight WER
increase $6.07\rightarrow 6.12$ on the source domain.

FedAvgM/BMUF: Next, we look at the effect of using the FL training algorithm,
by replacing DDP with FedAvgM/BMUF while keeping the labels and the loss
function the same (2). We observe a marginal WER increase on both the target
domain (to $4.09$) and the source domain (to $6.16$).

Table 3: The effect of data augmentation. When target domain is clean, adaptation is less effective without data augmentation. Model | Target | Source
---|---|---
3 | 4.35 | 6.17
3 \- Data Augment | 4.48 | 6.44
3 \- Data/Spec Augment | 4.89 | 6.67

Pseudo labels: Now, we introduce pseudo labels to determine the effects of
self-supervision (3). Although not as effective as the true labels, self-
supervision still reduces the pretrained model’s WER on the target domain
$4.96\rightarrow 4.35$. We note the data filtering techniques discussed in
section 2.3 greatly contributes to the success of self-supervision. Without
filtering, adaptation with the pseudo labels becomes much less effective,
causing a WER increase on the target domain $4.35\rightarrow 4.48$ (see table
3). If we also remove SpecAugment, the target domain WER further increases to
$4.89$. On the other hand, confidence filtering makes little difference since
the target domain WER of the pretrained model is already low.

SR-RNN-T: To reduce on-device computation cost, the proposed SR-RNN-T loss
restricts the forward/backward to a limited window of 2 frames around the most
probable alignment ($b_{l}=b_{r}=2$) according to the respective current
model. The SR-RNN-T trained model is marginally more accurate than the AR-
RNN-T on both target ($4.32$) and source domain ($6.08$).

Table 4: Adapting a subset of model weights. Finetuning only the key/value weights gives the best adaptation performance. Adaptation | Params (M) | Target | Source
---|---|---|---
All (4) | 72.4 (100%) | 4.32 | 6.08
Encoder | 67.4 (93%) | 4.31 | 6.15
Attention | 8.2 (11%) | 4.35 | 6.16
KeyValue (5) | 4.1 (5.7%) | 4.32 | 5.95
Predictor | 1.7 (2.5%) | 4.68 | 6.19
Joiner | 3.1 (4.4%) | 4.70 | 6.20
Bias | 0.2 (0.3%) | 4.35 | 5.75

KeyValue adaptation: Lastly, we experiment with finetuning only subsets of the
model weights. The results are summarized in table 4. Finetuning only the
key/value matrices of the self-attention module not only achieves a low target
domain WER ($4.32$), but also prevents catastrophic forgetting on the source
domain as the number of updates increases (see figure 3). Most importantly,
key/value matrices only account for $5.7\%$ of all weights. Hence,
communication cost is greatly reduced.

###  3.4 Video Adaptation

Figure 4: Convergence of Video adaptation experiments.

The Video adaptation setup adapts from the source domain “Clean” to the target
domain “Noisy”. The adaptation results are summarized in figures 4 and 5, and
they confirms most of our conclusions in section 3.3. In particular, FL
adaptation of key/value matrices with pseudo labels (5) effectively reduces
the WER on the target domain ($14.91\rightarrow 13.36$). Moreover, we have two
additional learnings.

Table 5: Summary of Video adaptation experiments. Model | Target | Source
---|---|---
Pretrained | 14.91 | 10.30
1 (DDP, True Labels, AR-RNN-T) | 12.71 | 10.26
2 (1 \+ FedAvgM/BMUF) | 13.09 | 10.36
3 (2 \+ Pseudo Labels) | 13.89 | 10.33
4 (3 \+ SR-RNN-T) | 13.74 | 10.08
5 (4 \+ KeyValue Adaptation) | 13.36 | 9.43
6 (3 \- Filtering) | 14.72 | 11.09

Filtering: Since the audios in the target domain is noisy, the pseudo labels
from the pretrained model have a significant error rate. We find it necessary
to apply confidence filtering using the decoding scores to remove low-
confidence examples. The optimal threshold parameter $\theta=-0.2$ removes
$27\%$ of training examples. As shown in table 5 row 6, semi-supervised
adaptation is much less effective when confidence filtering is removed, which
negates most of the gains ($13.89\rightarrow 14.72$).

KeyValue adaptation: Adapting only the key/value weights outperforms full
model adaptation by large margins on both target ($13.36$ vs. $13.74$) and
source domain ($9.43$ vs. $10.08$). Moreover, the adapted model even
outperforms the pretrained model on the source domain ($9.43$ vs. $10.30$).
This is likely because the target domain contains more diverse audio.

##  4 Conclusions

This paper addressed two practical challenges for domain adaptation of an
Emformer-Transducer ASR model on edge devices using Federated Learning. First,
we show that the combination of self-restricted RNN-T loss with data
augmentation and filtering techniques can improve an on-device ASR model on
target domain data with full self-supervision. Second, we explored techniques
reducing the FL compute and network communication costs by adapting only a
subset of weights in the model. We will further explore effects of non-IID-
ness and differential privacy on our proposed techniques as part of future
work.

## References

  * [1] F. Mireshghallah, M. Taram, P. Vepakomma, A. Singh, R. Raskar, and H. Esmaeilzadeh, “Privacy in deep learning: A survey,” _arXiv preprint arXiv:2004.12254_ , 2020.
  * [2] T. Ha, T. K. Dang, H. Le, and T. A. Truong, “Security and privacy issues in deep learning: a brief review,” _SN Computer Science_ , vol. 1, no. 5, pp. 1–15, 2020.
  * [3] B. Schuller and A. Batliner, _Computational paralinguistics: emotion, affect and personality in speech and language processing_. John Wiley & Sons, 2013.
  * [4] D. Cherkassky, “The voice privacy problem,” Kardome, Blog, 2021. [Online]. Available: https://www.kardome.com/blog-posts/voice-privacy-concerns
  * [5] T. Yang, G. Andrew, H. Eichner, H. Sun, W. Li, N. Kong, D. Ramage, and F. Beaufays, “Applied federated learning: Improving google keyboard query suggestions,” _arXiv preprint arXiv:1812.02903_ , 2018.
  * [6] A. Hard, K. Rao, R. Mathews, S. Ramaswamy, F. Beaufays, S. Augenstein, H. Eichner, C. Kiddon, and D. Ramage, “Federated learning for mobile keyboard prediction,” _arXiv preprint arXiv:1811.03604_ , 2018.
  * [7] D. Leroy, A. Coucke, T. Lavril, T. Gisselbrecht, and J. Dureau, “Federated learning for keyword spotting,” in _ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2019, pp. 6341–6345.
  * [8] S. Ramaswamy, R. Mathews, K. Rao, and F. Beaufays, “Federated learning for emoji prediction in a mobile keyboard,” _arXiv preprint arXiv:1906.04329_ , 2019.
  * [9] J. Konečnỳ, H. B. McMahan, F. X. Yu, P. Richtárik, A. T. Suresh, and D. Bacon, “Federated learning: Strategies for improving communication efficiency,” _arXiv preprint arXiv:1610.05492_ , 2016.
  * [10] B. McMahan, E. Moore, D. Ramage, S. Hampson, and B. A. y Arcas, “Communication-efficient learning of deep networks from decentralized data,” in _Artificial intelligence and statistics_. PMLR, 2017, pp. 1273–1282.
  * [11] K. Nandury, A. Mohan, and F. Weber, “Cross-silo federated training in the cloud with diversity scaling and semi-supervised learning,” in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 3085–3089.
  * [12] Q. Li, Z. Wen, Z. Wu, S. Hu, N. Wang, Y. Li, X. Liu, and B. He, “A survey on federated learning systems: vision, hype and reality for data privacy and protection,” _IEEE Transactions on Knowledge and Data Engineering_ , 2021\.
  * [13] D. Dimitriadis, K. Kumatani, R. Gmyr, Y. Gaur, and S. E. Eskimez, “A federated approach in training acoustic models.” in _Interspeech_ , 2020, pp. 981–985.
  * [14] Y. Gao, T. Parcollet, S. Zaiem, J. Fernandez-Marques, P. P. de Gusmao, D. J. Beutel, and N. D. Lane, “End-to-end speech recognition from federated acoustic models,” _arXiv preprint arXiv:2104.14297_ , 2021.
  * [15] W. Yu, J. Freiwald, S. Tewes, F. Huennemeyer, and D. Kolossa, “Federated learning in asr: Not as easy as you think,” in _Speech Communication; 14th ITG Conference_. VDE, 2021, pp. 1–5.
  * [16] D. Guliani, F. Beaufays, and G. Motta, “Training speech recognition models with federated learning: A quality/cost framework,” in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 3080–3084.
  * [17] X. Cui, S. Lu, and B. Kingsbury, “Federated acoustic modeling for automatic speech recognition,” in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 6748–6752.
  * [18] M. Paulik, M. Seigel, H. Mason, D. Telaar, J. Kluivers, R. van Dalen, C. W. Lau, L. Carlson, F. Granqvist, C. Vandevelde _et al._ , “Federated evaluation and tuning for on-device personalization: System design & applications,” _arXiv preprint arXiv:2102.08503_ , 2021.
  * [19] C. Dwork, A. Roth _et al._ , “The algorithmic foundations of differential privacy.” _Found. Trends Theor. Comput. Sci._ , vol. 9, no. 3-4, pp. 211–407, 2014.
  * [20] M. Abadi, A. Chu, I. Goodfellow, H. B. McMahan, I. Mironov, K. Talwar, and L. Zhang, “Deep learning with differential privacy,” in _Proceedings of the 2016 ACM SIGSAC conference on computer and communications security_ , 2016, pp. 308–318.
  * [21] V. Tsouvalas, A. Saeed, and T. Ozcelebi, “Federated self-training for semi-supervised audio recognition,” _arXiv preprint arXiv:2107.06877_ , 2021\.
  * [22] W. Jeong, J. Yoon, E. Yang, and S. J. Hwang, “Federated semi-supervised learning with inter-client consistency,” _arXiv preprint arXiv:2006.12097_ , 2020.
  * [23] C. He, Z. Yang, E. Mushtaq, S. Lee, M. Soltanolkotabi, and S. Avestimehr, “Ssfl: Tackling label deficiency in federated learning via personalized self-supervision,” _arXiv preprint arXiv:2110.02470_ , 2021.
  * [24] Y. He, T. N. Sainath, R. Prabhavalkar, I. McGraw, R. Alvarez, D. Zhao, D. Rybach, A. Kannan, Y. Wu, R. Pang _et al._ , “Streaming end-to-end speech recognition for mobile devices,” in _ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2019, pp. 6381–6385.
  * [25] K. Rao, H. Sak, and R. Prabhavalkar, “Exploring architectures, data and units for streaming end-to-end speech recognition with rnn-transducer,” in _2017 IEEE Automatic Speech Recognition and Understanding Workshop (ASRU)_. IEEE, 2017, pp. 193–199.
  * [26] A. Graves, “Sequence transduction with recurrent neural networks,” _arXiv preprint arXiv:1211.3711_ , 2012.
  * [27] J. Mahadeokar, Y. Shangguan, D. Le, G. Keren, H. Su, T. Le, C.-F. Yeh, C. Fuegen, and M. L. Seltzer, “Alignment restricted streaming recurrent neural network transducer,” in _2021 IEEE Spoken Language Technology Workshop (SLT)_. IEEE, 2021, pp. 52–59.
  * [28] J. Kim, H. Lu, A. Tripathi, Q. Zhang, and H. Sak, “Reducing streaming asr model delay with self alignment,” _arXiv preprint arXiv:2105.05005_ , 2021\.
  * [29] Y. Huang, G. Ye, J. Li, and Y. Gong, “Rapid speaker adaptation for conformer transducer: Attention and bias are all you need,” _Proc. Interspeech 2021_ , pp. 1309–1313, 2021.
  * [30] T.-J. Yang, D. Guliani, F. Beaufays, and G. Motta, “Partial variable training for efficient on-device federated learning,” _arXiv preprint arXiv:2110.05607_ , 2021.
  * [31] T.-M. H. Hsu, H. Qi, and M. Brown, “Measuring the effects of non-identical data distribution for federated visual classification,” _arXiv preprint arXiv:1909.06335_ , 2019.
  * [32] K. Chen and Q. Huo, “Scalable training of deep learning machines by incremental block training with intra-block parallel optimization and blockwise model-update filtering,” in _2016 ieee international conference on acoustics, speech and signal processing (icassp)_. IEEE, 2016, pp. 5880–5884.
  * [33] Y. Shi, Y. Wang, C. Wu, C.-F. Yeh, J. Chan, F. Zhang, D. Le, and M. Seltzer, “Emformer: Efficient memory transformer based acoustic model for low latency streaming speech recognition,” in _ICASSP 2021-2021 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2021, pp. 6783–6787.
  * [34] A. Xiao, W. Zheng, G. Keren, D. Le, F. Zhang, C. Fuegen, O. Kalinli, Y. Saraf, and A. Mohamed, “Scaling asr improves zero and few shot learning,” _arXiv preprint arXiv:2111.05948_ , 2021.
  * [35] D. S. Park, W. Chan, Y. Zhang, C.-C. Chiu, B. Zoph, E. D. Cubuk, and Q. V. Le, “Specaugment: A simple data augmentation method for automatic speech recognition,” _arXiv preprint arXiv:1904.08779_ , 2019.


# A Graph-Based Modeling Framework for Tracing
Hydrological Pollutant Transport in Surface Waters

David L. Cole†, Gerardo J. Ruiz-Mercado‡, and Victor M. Zavala† Corresponding
Author<EMAIL_ADDRESS>
(†Department of Chemical and Biological Engineering,
University of Wisconsin-Madison, Madison, WI 53706
‡Office of Research and Development,
U.S. Environmental Protection Agency, Cincinnati, OH 45268, USA,
and Chemical Engineering Graduate Program,
Universidad del Atlántico, Puerto Colombia 080007, Colombia )

###### Abstract

Anthropogenic pollution of hydrological systems affects diverse communities
and ecosystems around the world. Data analytics and modeling tools play a key
role in fighting this challenge, as they can help identify key sources as well
as trace transport and quantify impact within complex hydrological systems.
Several tools exist for simulating and tracing pollutant transport throughout
surface waters using detailed physical models; these tools are powerful, but
can be computationally intensive, require significant amounts of data to be
developed, and require expert knowledge for their use (ultimately limiting
application scope). In this work, we present a graph modeling framework—which
we call HydroGraphs—for understanding pollutant transport and fate across
waterbodies, rivers, and watersheds. This framework uses a simplified
representation of hydrological systems that can be constructed based purely on
open-source data (National Hydrography Dataset and Watershed Boundary
Dataset). The graph representation provides an flexible intuitive approach for
capturing connectivity and for identifying upstream pollutant sources and for
tracing downstream impacts within small and large hydrological systems.
Moreover, the graph representation can facilitate the use of advanced
algorithms and tools of graph theory, topology, optimization, and machine
learning to aid data analytics and decision-making. We demonstrate the
capabilities of our framework by using case studies in the State of Wisconsin;
here, we aim to identify upstream nutrient pollutant sources that arise from
agricultural practices and trace downstream impacts to waterbodies, rivers,
and streams. Our tool ultimately seeks to help stakeholders design effective
pollution prevention/mitigation practices and evaluate how surface waters
respond to such practices.

Keywords: graph theory, hydrology, connectivity, pollutants, nutrients,
watersheds, lakes, rivers.

Synopsis Statement We present a graph-based modeling framework for
hydrological systems that facilitates pollutant tracing. The graph
representation provides a simple, intuitive, and flexible way to rapidly
analyze pollutant transport in complex systems.

##  1 Introduction

Anthropogenic pollution in hydrological systems has significant impacts on
communities and ecosystems around the globe. These pollutants arise from
diverse sources and come in many forms; common pollutants include nutrients
such as nitrogen- and phosphorus-based fertilizers 1, emerging contaminants
(ECs) 2, microplastics 3, heavy metals 4, and microbes 5. These contaminants
can move through surface waters (e.g., lakes, rivers, streams) and
groundwaters by following complex network pathways, leading to impacts near
the contaminant release as well as to downstream impacts that span thousands
of miles. Studies have shown that many of these pollutants can be toxic to
humans and animals, and many of these impacts are still not fully understood.
For instance, ECs include chemicals such as pharmaceuticals, personal care
products, and per- and poly-fluoroalkyl substances (PFAS) 6. Risk assessment
studies on some ECs suggest that they can impact the immune system or cause
cancer, while other ECs have little published health information 7, 8. Other
pollutants, such as heavy metals or microbes, can likewise lead to cancer or
waterborne diseases, such as dysentery or diarrhea 9, 10.

Pollutants can also cause significant environmental and economic problems; for
instance, nutrient pollution is a major driver of harmful algae blooms (HABs)
in both freshwater and saline waterbodies 11, 12, 13. Nutrient pollution and
HABs can destroy marine wildlife through anoxia, poisoning, and other
mechanisms 14, 15, 16, can cause human health impact, and they can decrease
property values and hurt recreational and fishing operations on the order of
billions of US dollars 17, 18. Many pollutants—such as ECs, microplastics, and
heavy metals—have also been shown to bioaccumulate in wildlife, leading to
long-term effects 19, 2, 20. Furthermore, the fate of some pollutants such as
plastics is not fully understood; many of these contaminants have uncertain
environmental and human health impacts and can degrade into smaller compounds
with unknown properties 21, 2. These pollutants also tend to travel (via
hydrological systems) over long distances and find their way to oceans 22, 16.
Moreover, such pollutants may disproportionately impact vulnerable
communities, such as those in rural areas and developing countries 23.

To better understand and combat the environmental, economic, and health
impacts of hydrological pollutants, there is a need for tools that provide
intuitive and easy-to-use tools that can help navigate complexity and answer
questions of interest. For instance, for a given pollutant release, it is
important to understand which parts of a hydrological system will likely be
impacted or how far a pollutant can travel. Similarly, if a contaminant is
discovered in a given river or lake, we might be interested in identifying
what are possible upstream sources from which it originated. These types of
questions are challenging to address because hydrological systems involve
large and highly interconnected networks. For instance, the conterminous
United States contains more than 85,000 lakes 24 and over one million
kilometers of rivers and streams 25. Moreover, interconnections between
waterbodies span multiple spatial scales and are often non-intuitive. For
instance, the Mississippi river is connected to waterbodies in Pennsylvania
but not to waterbodies in Michigan. As a result, pollutant transport can also
be complex, as many pollutants can originate from point and non-point sources
1, 12, 26, can come from far upstream 27, can involve significant
spatial/temporal scales (i.e., legacy pollutants) 28, 29, 30, and can be
dependent on many factors such as weather, topology, soil type, or land cover
31, 32, 33.

Diverse tools exist for understanding pollutant transport using detailed
physical models 34, 35, 36, 6, 37, 38; these tools are powerful, but are
computationally intensive and may require significant expertise and data to be
used (ultimately limiting application scope). A modeling approach that can
help trace pollutant pathways in a more simplified manner consists of
representing hydrological systems as graphs (networks). Graphs are
mathematical representations (models) that are comprised of sets of nodes and
edges; nodes represent different objects (e.g., lakes), and edges are links
placed between nodes (e.g., rivers and streams connecting lakes). A wide range
of applications of graph representations have been explored in science and
engineering (from cosmology to social networks and infrastructure networks).
The success of graph modeling tools in such applications has been due to the
availability of diverse analysis/visualization tools and of underlying
algorithms that enable scalable analysis. In addition, graphs provide an
intuitive and flexible approach for analyzing connectivity, which can be an
important factor in understanding pollutant impacts and transport within
hydrological systems 39, 40, 41.

Graph models have been used in different studies to model hydrological
systems. Recently, King and co-workers 24 used graphs for analyzing
connectivity of lakes in the US; however, this work did not target
hydrological pollutant tracing and the graph representation used did not
capture rivers and streams as nodes within the graph. Other works have used
graphs to represent river systems, but these did not include explicit
connectivity to waterbodies or analyzed pollutant transport to waterbodies 42,
43, 44, 45, 46, 47. In addition, there are other datasets that capture
connectivity of some hydrological systems, such as the Watershed Index Online
(WSIO; captures connectivity of HUC12 watersheds) 48 or the Wisconsin
Department of Natural Resources (DNR) 24K Hydrography Geodatabase 49, 50.
These sources do not explicitly represent the waterbody-river system as a
graph, but they provide important data that can be incorporated into a graph
model of rivers and waterbodies. As will be shown later, a graph model can
provide an easy way to incorporate pollutant data, and it can be applied to
areas not covered by WSIO or the Wisconsin DNR 24K Hydrography Geodatabase.

In this work, we present a graph modeling framework—which we call
HydroGraphs—for capturing watershed-river-waterbody connectivity in
hydrological systems. The proposed framework was implemented in Python and was
developed with the goal of analyzing pollutant transport (Figure 1). A unique
aspect of our tool is that it allows the user to trace pollutant transport in
surface waters from a given watershed or pollutant source through downstream
waterbodies and watersheds. Here, we provide a detailed description of the
methods used for building a graph representation starting from public and
open-source databases such as the National Hydrography Dataset (NHDPlusV2) 51
and Watershed Boundary Dataset (WBD) 52 by using GeoPandas 53 and NetworkX 54.

We demonstrate the capabilities of HydroGraphs by providing case studies in
the State of Wisconsin, where our resulting graph contains over 45,000 nodes,
representing more than 3,000 km2 of waterbodies and 79,000 km of rivers and
streams. We also provide case studies showing how the graph framework can be
applied to analyze upstream sources and downstream impacts of anthropogenic
pollution. In these studies, we analyze nutrient pollution in Wisconsin, a
challenge that has lasted for decades and that originates from intensive
agricultural practices and other anthropogenic sources. We emphasize with
these studies how HydroGraphs can easily incorporate data for point and non-
point pollutant sources as well as impact data; in particular, we use data for
hundreds of concentrated animal feeding operations (CAFOs) as point sources
within the graph and include over 93,000 km2 of agricultural land as non-point
sources of pollution and we show how to link such source data to impact data
(chlorophyll-a concentration in lakes to monitor the onset of HABs). With
this, we aim to show how HydroGraphs can be a valuable tool for researchers
and decision-makers to conduct quick assessments of pollutant impacts on the
environment. Moreover, we discuss how the framework can be used in conjunction
with supply chain optimization models to understand how changes in
agricultural practices or in infrastructure can increase (or decrease) the
quality of hydrological systems.

Figure 1: An overview of HydroGraphs, a graph-based modeling framework for
incorporating GIS data and pollutant data to trace anthropogenic pollution
transport through surface waters. The framework uses open data from the
National Hydrography Dataset (NHDPlusV2) 51 and the Watershed Boundary Dataset
52, and uses GeoPandas 53 and NetworkX 54 for building and analyzing the
resulting graph representation. Total phosphorus and chlorophyll-a data shown
in the top right is from the Wisconsin Department of Natural Resources for
Lake Winnebago, WI55.

##  2 Graph Representation of Hydrological Systems

The focus of our work was to create a graph modeling tool, HydroGraphs, to
link upstream pollutant sources and downstream pollutant impacts.
Specifically, we aim to develop a tool that help us answer questions such as:
What waterbodies will be affected by a pollutant release in a specific
watershed? What upstream pollutant releases may be impacting a given
waterbody? What waterbodies may be $``$storing$"$ pollutants along a given
pathway? To answer such questions, it is necessary to model the connectivity
between the objects of interest (e.g., pollutant sources or waterbodies) and
creating simple and intuitive ways of analyzing these interconnections. Graphs
provide a natural mathematical representation to achieve these goals. In this
section, we provide an overview on graph representations, outline the data
needed to build such graphs, and outline the specific steps to express the
graph connectivity. Our framework can be used to capture diverse hydrological
systems in the United States (or in the world, provided that data is available
in the required format). In the next section, we illustrate these general
capabilities by building graph representations to capture interconnectivity of
surface waters in the State of Wisconsin and we explore pollutant tracing
applications.

An illustration of the methodology followed by our framework is provided in
Figure 2. We start with river segments (line features) and waterbodies
(polygon features) within watersheds (panel a). From existing data in the
NHDPlusV2, we create a graph of the river network (panel b) which does not
include any waterbodies. We then determine the connectivity of the waterbodies
and add them to the river graph (panel c).

Figure 2: Visualization of methods used for building a graph from geospatial
data. We start with river segments and waterbodies from the NHDPlusV2 dataset
(panel a) and form a directed graph of the river systems (panel b). We then
identify lake connectivity to the river graph and add waterbody nodes to the
directed graph (Panel c). Watersheds shown are the HUC10 watersheds
0709000205, 0709000206, and 0709000207, primarily in Dane County, Wisconsin.

###  2.1 Graph Theory Overview

Graphs are modeling abstractions that are comprised of a set of nodes and
edges. Nodes are used to represent diverse objects/elements of a system while
edges are used to model connectivity between nodes. As a graph is a
mathematical model, there are many different ways in which nodes and edges can
be defined and a given selection is often driven by the insights needed from
the model and from the data available. To build the watershed-river-waterbody
system of interest as a graph, we chose to represent river segments and
waterbodies as nodes. Edges are then placed between river segments and
waterbodies that flow into one another through surface waters (e.g., by rivers
or streams).

Nodes and edges of graphs can contain attributes (data) that are useful in
manipulating and visualizing a graph. For example, an attribute that we used
in our model is the watershed in which the node resides (the node encodes
spatial/geographical context); this attribute can be used to filter out nodes
that lie on a specific watershed. Our representation uses a directed (rather
than undirected) graph. In undirected graphs, edges only capture connectivity
(with no notion of directionality); while, in directed graphs, edges capture
directionality. In our context, we are interested in tracing nutrient
pollution, and we thus need to capture flow directionality.

A key benefit of using graph representations is that there are a wide range of
theory and computational techniques for analyzing large-scale graphs; for
instance, one can use algorithms to identify a set of nodes that is connected
to a given node by using pathway analysis. Moreover, it is possible to
visualize, aggregate, and partition graphs to gain insights into the
connectivity and properties of a graph. In addition, it is possible to compute
statistics of a given graph object, such as the number of nodes in a given
pathway, the fractal dimension of a graph (e.g., a measure of complexity), or
the node degree distribution (e.g., number of connections of a node). Another
key benefit of using graph representations is that there are a wide range of
open-source tools that can be leveraged for building and visualizing graphs.

###  2.2 Data Overview

The graph representation was constructed directly using GIS data from the
NHDPlusV2 and the WBD datasets. While the data within these datasets are for
the United States, similar methods could be applied to other geographical
areas provided that the data is in similar formats as what is discussed below.
We used GeoPandas 53 in Python to work with the GIS data, and we use NetworkX
54 for building and managing the graph. The code for our framework is
available at https://github.com/zavalab/JuliaBox/tree/master/HydroGraphs.

In our representation there are three main types of objects used to build the
graph (see Figure 2a): rivers (NHDFlowline from NHDPlusV2; represented by line
features), waterbodies (NHDWaterbody from NHDPlusV2; represented by polygon
features), and watersheds (from WBD; represented by polygon features). Each of
these lines or polygons has a geographic location (e.g., edges of the polygons
were represented by specific geographic coordinates) and has a unique
identifier. The NHDPlusV2 dataset uses a unique common identifier (COMID) for
every individual river segment or waterbody, while the WBD uses a unique
hydrologic unit code (HUC) for individual watersheds. In addition, the WBD has
different hierarchical levels, with 8-digit, 10-digit, and 12-digit codes
depending on the size of the watersheds, where the higher digit codes are
partitions of the lower digit codes (e.g., each HUC10 watershed is made up of
multiple, smaller HUC12 watersheds). We will use these identifiers in talking
about their corresponding objects; for example, we will form nodes out of each
river segment or waterbody, and we will identify these nodes by their
corresponding COMID.

###  2.3 Expressing Connectivity

One of the key technical challenges in building the graph representation is
identifying the connectivity between the rivers, waterbodies, and watersheds.
Conveniently, the NHDPlusV2 dataset provides connectivity between river
segments by giving directed pairs of river segments identified by their COMID
(i.e., these are given as pairs of $``$from$"$ COMIDs and their corresponding
$``$to$"$ COMIDs). This is essentially a list of directed edges of a graph and
could be used to build a graph of the river system. In this case, each river
segment corresponds to a node, where the node is identified by the river
segment’s COMID. The graph formed by this list of directed edges (which we
will refer to as the ”river graph”; see Figure 2b) was a basis for building
our overall graph. Note, however, that the river graph does not include
specific nodes that correspond to waterbodies ($``$waterbody nodes$"$ will be
added in a later process). The NHDPlusV2 data overlaps river segments with
waterbody polygons, and many polygons may overlap with several river segments.
Furthermore, some river segments may overlap with multiple waterbodies, making
it difficult to identify which waterbodies flow into one another. Our methods
for identifying which river segments overlapped with waterbodies and for
adding waterbodies to the given river segment graph are outlined in this
subsection.

The first step in adding the waterbodies to the graph was identifying the
connectivity of waterbodies and rivers. We first built empty lists for every
river segment and waterbody polygon. These lists would contain the COMIDs of
all intersecting rivers (for waterbodies) or all intersecting waterbodies (for
river segments). We then tested every river segment against every waterbody to
see if the river segment intersected the waterbody polygon. This was performed
within two loops using the GeoPandas’ intersects function. If a river segment
intersected a waterbody, the river segment COMID was added to the waterbody’s
list of intersecting rivers, and the waterbody COMID was added to the river
segment’s list of intersecting lakes. These lists could then be used to add
the respective lakes to the graph.

Figure 3: Visualization of methods used for representing the river-waterbody
system as a graph. The red lines represent river segments and are identified
by numbers while the blue polygons represent waterbodies and are identified
with letters. River segment 1 only intersects waterbody A, so it is replaced
by a single node corresponding to waterbody A. River segment 4 intersects two
waterbodies, B and C, so it is replaced by two nodes within the graph
corresponding to waterbodies B and C.

After identifying all intersections between rivers and waterbodies, we added
waterbodies to the river graph by replacing river COMIDs with waterbody COMIDs
(i.e., replacing river nodes with waterbody nodes) and adding edges between
the resulting waterbody nodes and other river nodes (See Figure 2c). We first
looped through every river segment; if the list of intersecting waterbodies
for a given river segment contained only one waterbody, we replaced the river
segment COMID in the river graph’s edge list with the intersecting waterbody
COMID. This added several waterbodies to the river graph. After completing
this loop, we looped over every river segment again. If the river segment
intersected multiple waterbodies, then the river segment was replaced by all
waterbodies which it intersected. This meant that more nodes were added to the
graph. A visualization of this process is given in Figure 3. Note that this
results in waterbodies that intersect the same river segment not necessarily
being directly connected (e.g., waterbody B in Figure 3 should connect to
waterbody C, but it does not). This was a simplification that had to be made
because individual river segments do not have an inherent flow direction to
them, meaning that we cannot tell which waterbody is intersected first by the
river segment. We believe that this is a reasonable simplification because the
general connectivity of the full graph is still maintained; moreover, this
highlights how the graph representation used is inherently limited by the
availability of data. The affected waterbodies still have the same upstream
and downstream connections with the minor exception of not being connected to
waterbody(ies) that intersect their same line segment. For example, waterbody
A still flows into waterbodies B and C in Figure 3 even though waterbodies B
and C are not directly connected to each other. In other words, waterbody A’s
downstream graph (and river node 5’s upstream graph) includes the same set of
nodes as it would if this simplification was not applied. For the State of
Wisconsin, this latter simplification applied to less than 10% of the
waterbodies in the graph. This simplification is very localized (average river
segment length is $<$ 2 km, and the impacts are only on waterbodies connected
by the same river segment), and thus we believe this will produce minimal
error. Further details about this method can be found in the supporting
information. The resulting graph obtained for Wisconsin is presented in Figure
4.

Figure 4: Wisconsin hydrological system (a), representation as a directed
graph (b), and aggregated form of the directed graph (c). The full directed
graph contains over 45,000 nodes and 47,000 edges, representing more than
3,000 km2 of waterbodies and 79,000 km of rivers or streams.

We make a few additional remarks about our data and methods. First, we chose
to exclude a few waterbodies from our graph. We excluded Lake Michigan and
Lake Superior from the graph to make the visualizations simpler. Ultimately,
for the Wisconsin area explored in the next section, everything that does not
flow into the Mississippi River flows into these two lakes (i.e., all nodes in
the 04 HUC2 watershed are connected to the Great Lakes), so the connectivity
to these lakes is established by nodes being within the 04 HUC2 watershed.
Representing either of these lakes by a single node makes the visualizations
more difficult to follow, and thus were undesirable for our needs. However,
these could be included in our methods simply by not removing the COMIDs for
these two lakes from the original dataset.

In addition, while swamps and marshes are included in the NHDWaterbody
dataset, we removed these objects to simplify the analysis. Swamps and marshes
impact nutrient transport differently than some other waterbodies, so we did
not want to include them in the same category with lakes and reservoirs.
However, we note that swamps and marshes could be influential in nutrient
transport (in fact they can be used to control nutrient pollution) 56, 57, 58,
59. Thus, these will be a subject of future research, but they are outside the
scope of this study. Wetlands like swamps and marshes can introduce
significant complexity because the pollutant transport can be dependent on
soil type and vegetative processes (as is the case for nutrient pollution 57,
59), and they could significantly impact the time scales of the pollutant
transport. This study does not focus on the temporal aspect of hydrological
pollutant transport, but rather on the pathways and fate of the pollutants in
hydrological systems.

We also note that much of the connectivity we give here could be elucidated
from the LAGOS-US NETWORKS v1 data set 24. The LAGOS-US NETWORKS v1 dataset
includes information on the lakes that are connected to the river graph edge
list (see nets_flow_medres.csv60) by indicating whether an edge of the river
graph also goes to or from a lake in the dataset. However, we chose to build
our connectivity list from scratch because we wanted to include several
waterbodies that were not included in the LAGOS-US NETWORKS v1 dataset. For
example, they omit lakes that are $<1$ hectare in size, and they do not
include reservoirs (NHDWaterbodies attribute FTYPE equal to $``$Reservoir$"$).
Both of these sets of waterbodies could be areas that anthropogenic pollutants
accumulate and could be a focus of pollutant studies (see for example 61, 62).
Excluding these waterbodies could thus lead to incorrect results when seeking
to identify areas of pollutant impacts.

We recognize that the above methods are only focused on waterbodies that are
connected through surface waters. The resulting graph outlined above does not
include every waterbody in a geographic area because many waterbodies are
isolated and not connected by surface waters to other objects in the graph.
Furthermore, this graph is specific to surface water and does not include
transport through other means such as groundwater, which can be important
factors 63, 64, 65. Including other transport mechanisms greatly impacts the
complexity and will be explored in future work.

###  2.4 Aggregating River Nodes

The above methods for building this graph result in several intermediate river
nodes upstream or downstream of waterbodies. In many cases, it may be
desirable to aggregate nodes to simplify the model representation, either to
make the visualizations simpler or to reduce the number of nodes involved when
analyzing the graph with varying degrees of spatial resolution. Our framework
provides capabilities for automating aggregation; details of this aggregation
procedure are included in the supporting information.

##  3 Wisconsin Case Studies

We highlight how HydroGraphs can be used to identify pollutant sources and
their potential destinations; we do this by developing some specific case
studies. Case studies focus on nutrient pollution in Wisconsin, a challenge
that has existed for decades due in part to the large amount of agricultural
land and CAFOs throughout the state that result in nitrogen (N) and phosphorus
(P) flowing into nearby waterways. Nutrient losses from these sources
frequently lead to HABs, which can have negative health, economic, and
environmental impacts for local communities.

The first case study compares a couple of lakes in Wisconsin with differing
total phosphorus (TP) concentrations and looks at their upstream graphs and
likely P sources that contribute to these differences. The second case study
looks more generally at several hundred lakes for which we have TP and
chlorophyll-a data and compares connectivity attributes of the graph between
polluted and clean lakes. The final case study looks at how our framework can
be used to identify impacts that a potential pollutant source could have. This
is done by inspecting the nodes in the graph that are downstream of the
source.

The case studies presented herein are intended as examples of ways that this
graph framework could be applied. While these case studies do incorporate real
data, they are not rigorous studies intended to give exact causation or to
make policy recommendations. Rather, their purpose is to present how the graph
could be applied by experts and researchers in this field. Further, they
highlight how HydroGraphs could be used in helping decision-makers approach
complex problems involving pollutants in hydrological systems. For the code to
replicate these case studies, see
https://github.com/zavalab/JuliaBox/tree/master/HydroGraphs.

###  3.1 Case Study I: Identifying Upstream Sources

In this case study, we look at how this graph can enable identifying upstream
influences to a given waterbody. We focus in this case on P pollution in
waterbodies, but the principles in this case study could easily be applied to
other pollutants, such as ECs, microplastics, or heavy metals. Here, we build
the upstream graphs for two lakes in Wisconsin—Lake Altoona and Mohawksin
Lake—and identify possible pollutant sources that contribute to these lakes.
We choose these lakes because we have TP concentration measurements from the
Wisconsin DNR for each lake 66, 67. Based on the measured TP concentrations
and reported lake perception, Lake Altoona has poorer water quality (average
measured TP of 103 mg/m3) and has noticeably worse problems with algae than
Mohawksin Lake (average measured TP of 40 mg/m3). Data for these lakes and
their reported perception measurements are available in the supporting
information.

There are several upstream factors that can impact pollutant transport to
waterbodies. We look at three specific factors that may influence the TP
concentrations within the waterbodies. The first factor is the waterbodies
that are upstream to a given waterbody. Upstream waterbodies could accumulate
pollutants and could impact how much of a pollutant reaches a downstream
waterbody and when it reaches it 39, 28, 68. These waterbodies are naturally a
part of the upstream graph from a specific waterbody. The second factor we
consider are CAFOs; these can be a source of pollutants in Wisconsin because
of their large production of manure, and many sources show that these can be
significant contributors of P to waterbodies 69, 70, 71. For our analysis
here, we use data from Hu and co-workers72 to identify locations for more than
200 CAFOs. We add the CAFOs from this dataset to the directed graph outlined
above by adding a directed edge from the location of the CAFO to the closest
node in the same HUC12 watershed as the given CAFO. The other source we
consider is agricultural land; agricultural land is a significant non-point
source of N/P and can be an indicator of surface water pollutant
concentrations 1, 73, 74, 75. This source can be closely related to CAFOs
because the manure from CAFOs is often applied as a nutrient source for crops.
We use shapefiles from 76 to identify agricultural land, and we only look at
agricultural land that shares a HUC12 watershed with the given waterbody or
its upstream nodes.

Because we are using a directed graph representation, we can easily form the
upstream graphs for both Altoona Lake and Lake Mohawksin. This allows us to
frame potential hypotheses of why Altoona Lake has much higher TP
concentrations than Lake Mohawksin. These upstream graphs for these lakes can
be seen in Figures 5 and 6. Panel a) shows all upstream nodes (including
CAFOs) overlayed on the HUC12 watersheds and where these watersheds lie in
Wisconsin. Panel b) includes the agricultural land polygons. These figures
also include the upstream waterbody polygons to give an idea of the size of
the waterbodies to which waterbody nodes correspond.

From Figures 5 and 6, it is clear that both Altoona Lake and Lake Mohawksin
have large upstream graphs that include several waterbodies. Altoona Lake
includes 34 upstream waterbodies that cover an area of more than 7 km2, while
Lake Mohawksin includes 265 upstream waterbodies covering more than 345 km2.
Altoona Lake also has upstream connections to two CAFOs while Mohawksin Lake
has no CAFOs in its upstream watersheds. The figures also show that there is
significantly more agricultural land in Altoona Lake’s upstream graph. The
agricultural land makes up 65.4% of Altoona Lake upstream watersheds while it
makes up only 30.8% of Lake Mohawksin’s upstream watersheds. The addition of
CAFOs and the higher fraction of agricultural land can thus be likely
contributors to the high TP concentrations of Altoona Lake.

Figure 5: Upstream graph of Altoona Lake in Wisconsin. Panel a) shows all
upstream river and waterbody nodes along with CAFOS connected to the closest
node within their HUC12 watershed. Panel b) includes agricultural land
polygons in green. Figure 6: Upstream graph of Lake Mohawksin in Wisconsin.
Panel a) shows all upstream river and waterbody nodes. Panel b) includes
agricultural land polygons in green.

The above analysis gives some examples of how HydroGraphs could be used for
analysis. It is very possible that the high amount of agricultural land
(agricultural land fraction is twice as high for Lake Altoona as Mohawksin
Lake) and the upstream CAFOs (two CAFOs for Lake Altoona compared to zero
CAFOs for Mohawksin Lake) are contributors to the high TP concentrations
within Altoona Lake. Further, the high number of upstream waterbodies for Lake
Mohawksin may influence the transfer of pollutants to Lake Mohawksin (e.g.,
through the accumulation of nutrient pollution in upstream waterbody
sediments) 39, 77, 78. Formulating these systems as a graph enables the above
visualizations and simplifies analysis. It makes it easier to identify
upstream point sources that could contribute to pollutant concentrations in
waterbodies. In addition, just as CAFOs were added to the graph, other
pollutant sources (such as wastewater treatment plants 79, 80 or landfills 72)
could also be added. This could be useful in the event that a pollutant, such
as an EC, is discovered in a given lake or stream. Building the upstream graph
would allow researchers and decision makers to identify where this pollutant
may be coming from, and to identify other upstream waterbodies or rivers that
may need to be tested to see if they are likewise contaminated. As seen from
the CAFOs in Figure 5, some of these upstream pollutant sources could be far
upstream but are more easily identifiable by building the graph.

###  3.2 Case Study II: Graph Connectivity Metrics

In this case study, we study upstream graph metrics for waterbodies for which
we have TP and chlorophyll-a (a measure that relates to the level of algae in
the lake) data from the DNR 55. The DNR provides water quality data collected
by volunteers for hundreds of waterbodies in Wisconsin. We compiled their data
for more than 700 unique waterbodies. Data for each waterbody could vary in
terms of frequency of measurements and type of measurements taken. To ensure
that the lakes considered had significant data, we only studied waterbodies
for which there were at least 50 measurements for both TP and chlorophyll-a
which resulted in a set of 241 waterbodies. We chose the cutoff of 50 to
ensure that we had several data of both TP and chlorophyll-a while maintaining
a reasonable subset of lakes to analyze (e.g., choosing a higher cutoff such
as 100 resulted in too few lakes, while choosing a smaller cutoff like 10
could result in too little data for the lakes). For more details on how this
data was compiled and what was included, see the Supporting Information.

Based on this data for TP and chlorophyll-a, we studied five waterbodies that
had the highest and lowest levels of TP and chlorophyll-a within our graph.
The times at which lakes were sampled varied between waterbodies; as such, we
averaged all reported measurements of TP and chlorophyll-a. This is a primary
reason why we required that lakes have at least 50 data points for both TP and
chlorophyll-a as we assumed that this would reduce or eliminate any
differences between waterbodies caused by different sampling times. We then
looked at the five waterbodies that were in our graph that had the highest
average TP measurement and an average chlorophyll-a measurement of at least 30
mg/m3, and we looked at the five waterbodies with the lowest average TP
measurement and an average chlorophyll-a measurement of no more than 5 mg/m3.
The cutoff values for chlorophyll-a were chosen to ensure that the given
waterbodies were at the high or low extremes of both TP and chlorophyll-a
within our dataset (for further analysis of and explanation for these cutoff
values, please see the Supporting Information). After determining these five
most and least polluted waterbodies for which we have data, we built the
upstream graphs using similar methods as discussed in Case Study I, where we
added CAFOs and agricultural land to the visualizations. The results can be
seen in Figure 7.

Figure 7: Comparison of upstream graphs for five polluted waterbodies and five
clean waterbodies, as determined by data from the Wisconsin DNR 55. ”TP avg”
and ”chla avg” are the average TP and chlorophyll-a measurement, and ”ag land”
is the agricultural land fraction in the HUC12 watersheds comprising the
graph.

By creating the graphs in Figure 7, we are able to observe some trends between
the polluted and clean waterbodies. The polluted waterbodies all had
agricultural land fractions higher than the clean waterbodies, and no polluted
waterbody had an agricultural land fraction lower than 80%. Both Lake Tomah
and Puckaway Lake (polluted waterbodies) had upstream CAFOs while no clean
waterbody had an upstream CAFO. The polluted waterbodies also had three
waterbodies with significant upstream connections. In addition, there were two
waterbodies in the polluted waterbodies (Carstens Lake and Little Green Lake)
and one in the clean waterbodies (Maiden Lake) that had no upstream
connections. This is because they are in the graph, but only have downstream
connections.

The analysis of Figure 7 can also be expanded to include more waterbodies. To
look at how well these metrics apply to other waterbodies, we took a bigger
subset of waterbodies by taking as ”clean” all waterbodies with TP
concentrations of $<$ 15 mg/m3 and chlorophyll-a concentrations of $<$ 5 mg/m3
which totaled 60 waterbodies. For the ”polluted”, we used all waterbodies with
TP concentrations of $>$ 60 mg/m3 and chlorophyll-a concentrations of $>$ 15
mg/m3 which totaled 18 waterbodies (for further analysis of and explanation
for these cutoff values, please see the Supporting Information). We then
looked at how many of these waterbodies are connected to CAFOs upstream (or,
in the case of waterbodies without an upstream graph, if there is a CAFO in
their HUC12 watershed), what is their agricultural land fraction, how many
waterbodies are not in the graph (i.e., how many are not connected to a river
or stream in the graph), and how many upstream nodes these waterbodies have
that were in the graph. The results of this analysis are shown in Table 1. In
addition, we also found that the average agricultural land fraction for the
polluted waterbodies was 73.1% while the average for the clean waterbodies was
41.7%.

Table 1: Metrics for Polluted and Clean waterbodies using HydroGraphs. Metrics indicate the number of waterbodies that are connected upstream to a CAFO (CAFO); the number of waterbodies that are in the graph (In Graph); the number of waterbodies that could be considered headwaters as they are in the graph but have no upstream nodes (Headwater); the number of waterbodies that have an agricultural land fraction of more than 50% in their watershed ($>$ 50% ag frac); and the number of waterbodies that have at least 10 nodes in their upstream graph (10$+$ nodes). |  | Total | CAFO | In Graph | Headwater |  $>$ 50% ag frac | 10$+$ nodes
---|---|---|---|---|---|---|---
Polluted Waterbodies | Number | 18 | 6 | 15 | 2 | 17 | 9
Fraction |  | 0.33 | 0.83 | 0.11 | 0.94 | 0.50
Clean Waterbodies | Number | 60 | 4 | 24 | 10 | 28 | 3
Fraction |  | 0.07 | 0.40 | 0.17 | 0.47 | 0.05

Overall, this analysis suggests that some graph metrics may be feasible
indicators that help identify polluted lakes. One third of the polluted
waterbodies were connected to an upstream CAFO, and the polluted waterbodies
were generally in watersheds with much higher agricultural land fractions.
Further, the polluted waterbodies often exhibited much higher upstream
connectivity than the clean waterbodies. These metrics are all relatively easy
to compute using the graph representation of the hydrological system, and the
results could potentially be extrapolated to other lakes for which we do not
have data compiled. Building these systems as a graph provides new tools for
pollutant fate and transport. We would like to highlight that the results
found in this study are only speculative and do not aim to provide a final
recommendation on the origins of Lake nutrient pollution (which can be the
result of many factors).

###  3.3 Case Study III: Downstream Impacts

For this final case study, we consider a hypothetical example of placing a
potential pollutant source (e.g., a new CAFO or wastewater treatment plant)
and identifying its downstream impacts. In determining where to place this new
pollutant source, we need to consider the potential downstream impacts this
could have. We look at a couple of potential locations for its placement
within the same area. Interestingly, the locations are just within five
kilometers of each other, but they have drastically different downstream
destinations (because they lie in different watersheds). These pollutant
sources are added to the graph by placing a directed edge from the pollutant
source to the nearest node within the pollutant source’s HUC12 watershed. The
downstream, aggregated graphs for these two pollutant sources are shown in
Figure 8.

Figure 8: Downstream aggregated graphs for two potential pollutant sources.
Source locations are separated by less than 5 km

Location 1 (Figure 8a) results in the pollutant passing through 15 waterbodies
and going into the Mississippi River (the Western border of Wisconsin), while
location 2 (Figure 8b) shows the pollutant also passing through 15 waterbodies
and moving into Lake Michigan. This analysis can be useful in identifying
waterbodies potentially impacted by new pollutant releases. In the event of a
contaminant release, this graph can help identify waterbodies that may be
impacted by such release and this can help develop mitigation/response
strategies. Furthermore, the graph can also be useful for decision-makers in
identifying where to place potential pollutant sources (such as building a new
wastewater treatment plant) to minimize that sources impacts on the
environment and communities. The locations shown in Figure 8 travel through
completely different waterbodies, and these waterbodies may be at different
stages of eutrophication and may be of varying importance/priority (e.g., if
some waterbodies serve as drinking water for a locality, it may be more
important than another waterbody). Thus, building the downstream graph enables
researchers or decision-makers to see potential impacts of a newly introduced
pollutant source.

###  3.4 Conclusions and Future Work

HydroGraphs provides a framework for analyzing hydrological pollution pathways
to identify upstream sources and downstream impacts. The above case studies
have shown how HydroGraphs can be used with point and non-point pollution
sources to identify upstream sources and link attributes within the graph to
pollutant data. Further, it can also help identify potential downstream
impacts from a given pollutant source. While the case studies in this paper
focused on anthropogenic nutrient pollution, similar methods could be applied
for other pollutants such as ECs or microplastics. Point or non-point sources
of these contaminants could be added to the graph following a similar analysis
as that done with CAFOs and agricultural land. Building these hydrological
systems as a graph ultimately provides simple visualization and rapid
analysis.

There are two areas we would like to address in the future using HydroGraphs.
First, there is additional data that we can incorporate into the graph. For
example, the NHDPlusV2 dataset includes attributes such as average stream
flowrates that could be added to the graph as edge weights. As waterbodies
often have multiple streams flowing into them, these flowrates could show
which upstream sources have a stronger impact on the waterbody (e.g., the
upstream sources connected through the larger stream may have a stronger
impact). Second, we would like to incorporate HydroGraphs into additional
decision-making models, such as into supply chain optimization. For example,
Tominac et al. 81 included environmental policy-makers as stakeholders within
their supply chain model. One of the challenges in doing so is quantifying the
environmental or social impacts within the supply chain. HydroGraphs could
provide a tool for quantifying these impacts, such as quantifying the number
of lakes that would be impacted by the introduction of a new pollutant source.
This would allow for these decision-making models to highlight the
environmental, economic, and social impacts of many pollutant sources.

## Supporting Information

Additional methodological details on graph construction and aggregation,
details of the DNR data presented herein, and an overview of functionality for
working with the graph representation are provided in the SI.

## Acknowledgments

We acknowledge support from the U.S. EPA (contract number EP-18-C-000016). We
thank Eric Booth for helpful feedback on an early version of this manuscript.

The views expressed in this article are those of the authors and do not
necessarily reflect the views or policies of the U.S. Environmental Protection
Agency. Mention of trade names, products, or services does not convey, and
should not be interpreted as conveying, official U.S. EPA approval,
endorsement, or recommendation.

## References

  * Carpenter et al. 1998 Carpenter, S. R.; Caraco, N. F.; Correll, D. L.; Howarth, R. W.; Sharpley, A. N.; Smith, V. H. Nonpoint pollution of surface waters with phosphorus and nitrogen. _Ecological applications_ 1998, _8_ , 559–568
  * Wilkinson et al. 2017 Wilkinson, J.; Hooda, P. S.; Barker, J.; Barton, S.; Swinden, J. Occurrence, fate and transformation of emerging contaminants in water: An overarching review of the field. _Environmental Pollution_ 2017, _231_ , 954–970
  * Haddout et al. 2022 Haddout, S.; Gimiliani, G.; Priya, K.; Hoguane, A.; Casila, J. C. C.; Ljubenkov, I. Microplastics in surface waters and sediments in the sebou estuary and Atlantic Coast, Morocco. _Analytical Letters_ 2022, _55_ , 256–268
  * Ciazela et al. 2018 Ciazela, J.; Siepak, M.; Wojtowicz, P. Tracking heavy metal contamination in a complex river-oxbow lake system: Middle Odra Valley, Germany/Poland. _Science of the Total Environment_ 2018, _616_ , 996–1006
  * Nawab et al. 2016 Nawab, J.; Khan, S.; Ali, S.; Sher, H.; Rahman, Z.; Khan, K.; Tang, J.; Ahmad, A. Health risk assessment of heavy metals and bacterial contamination in drinking water sources: a case study of Malakand Agency, Pakistan. _Environmental monitoring and assessment_ 2016, _188_ , 1–12
  * Tong et al. 2022 Tong, X.; Mohapatra, S.; Zhang, J.; Tran, N. H.; You, L.; He, Y.; Gin, K. Y.-H. Source, fate, transport and modelling of selected emerging contaminants in the aquatic environment: Current status and future perspectives. _Water Research_ 2022, 118418
  * Bonato et al. 2020 Bonato, M.; Corrà, F.; Bellio, M.; Guidolin, L.; Tallandini, L.; Irato, P.; Santovito, G. PFAS environmental pollution and antioxidant responses: an overview of the impact on human field. _International journal of environmental research and public health_ 2020, _17_ , 8020
  * Cousins et al. 2020 Cousins, I. T.; DeWitt, J. C.; Glüge, J.; Goldenman, G.; Herzke, D.; Lohmann, R.; Miller, M.; Ng, C. A.; Scheringer, M.; Vierke, L.; Wang, Z. Strategies for grouping per-and polyfluoroalkyl substances (PFAS) to protect human and environmental health. _Environmental Science: Processes & Impacts_ 2020, _22_ , 1444–1460
  * Khan et al. 2013 Khan, S.; Shahnaz, M.; Jehan, N.; Rehman, S.; Shah, M. T.; Din, I. Drinking water quality and human health risk in Charsadda District, Pakistan. _Journal of cleaner production_ 2013, _60_ , 93–101
  * Lim et al. 2008 Lim, H.-S.; Lee, J.-S.; Chon, H.-T.; Sager, M. Heavy metal contamination and health risk assessment in the vicinity of the abandoned Songcheon Au–Ag mine in Korea. _Journal of Geochemical Exploration_ 2008, _96_ , 223–230
  * Committee on the Causes and Management of Coastal Eutrophication; Ocean Studies Board; Water Science and Technology Board; Commission on Geosciences, Environment, and Resources; National Research Council 2000 Committee on the Causes and Management of Coastal Eutrophication; Ocean Studies Board; Water Science and Technology Board; Commission on Geosciences, Environment, and Resources; National Research Council, _Clean coastal waters: understanding and reducing the effects of nutrient pollution_ ; National Academies Press, 2000
  * Nie et al. 2018 Nie, J.; Feng, H.; Witherell, B. B.; Alebus, M.; Mahajan, M. D.; Zhang, W.; Yu, L. Causes, assessment, and treatment of nutrient (N and P) pollution in rivers, estuaries, and coastal waters. _Current pollution reports_ 2018, _4_ , 154–161
  * Shortle and Horan 2017 Shortle, J.; Horan, R. D. Nutrient pollution: A wicked challenge for economic instruments. _Water Economics and Policy_ 2017, _3_ , 1650033
  * Bauman et al. 2010 Bauman, A. G.; Burt, J. A.; Feary, D. A.; Marquis, E.; Usseglio, P. Tropical harmful algal blooms: An emerging threat to coral reef communities? _Marine pollution bulletin_ 2010, _60_ , 2117–2122
  * Bruslé 1994 Bruslé, J. The impact of harmful algal blooms on finfish. Mortality, pathology and toxicology. _Repères Océans_ 1994,
  * Rabotyagov et al. 2020 Rabotyagov, S. S.; Kling, C. L.; Gassman, P. W.; Rabalais, N. N.; Turner, R. E. The economics of dead zones: Causes, impacts, policy challenges, and a model of the Gulf of Mexico hypoxic zone. _Review of Environmental Economics and Policy_ 2020,
  * Dodds et al. 2009 Dodds, W. K.; Bouska, W. W.; Eitzmann, J. L.; Pilger, T. J.; Pitts, K. L.; Riley, A. J.; Schloesser, J. T.; Thornbrugh, D. J. Eutrophication of US freshwaters: Analysis of potential economic damages. _Environmental Science and Technology_ 2009, _43_ , 12–19
  * Sampat et al. 2021 Sampat, A. M.; Hicks, A.; Ruiz-Mercado, G. J.; Zavala, V. M. Valuing economic impact reductions of nutrient pollution from livestock waste. _Resources, Conservation and Recycling_ 2021, _164_ , 105199
  * Copat et al. 2012 Copat, C.; Bella, F.; Castaing, M.; Fallico, R.; Sciacca, S.; Ferrante, M. Heavy metals concentrations in fish from Sicily (Mediterranean Sea) and evaluation of possible health risks to consumers. _Bulletin of Environmental Contamination and Toxicology_ 2012, _88_ , 78–83
  * Zhang et al. 2019 Zhang, S.; Ding, J.; Razanajatovo, R. M.; Jiang, H.; Zou, H.; Zhu, W. Interactive effects of polystyrene microplastics and roxithromycin on bioaccumulation and biochemical status in the freshwater fish red tilapia (Oreochromis niloticus). _Science of the total environment_ 2019, _648_ , 1431–1439
  * Gogoi et al. 2018 Gogoi, A.; Mazumder, P.; Tyagi, V. K.; Chaminda, G. T.; An, A. K.; Kumar, M. Occurrence and fate of emerging contaminants in water environment: A review. _Groundwater for Sustainable Development_ 2018, _6_ , 169–180
  * Ho et al. 2019 Ho, J. C.; Michalak, A. M.; Pahlevan, N. Widespread global increase in intense lake phytoplankton blooms since the 1980s. _Nature_ 2019, _574_ , 667–670
  * Ashbolt 2004 Ashbolt, N. J. Microbial contamination of drinking water and disease outcomes in developing regions. _Toxicology_ 2004, _198_ , 229–238
  * King et al. 2021 King, K. B.; Wang, Q.; Rodriguez, L. K.; Cheruvelil, K. S. Lake networks and connectivity metrics for the conterminous US (LAGOS-US NETWORKS v1). _Limnology and Oceanography Letters_ 2021, _6_ , 293–307
  * U.S. Environmental Protection Agency 2020 U.S. Environmental Protection Agency, _National rivers and streams assessment 2013-2014: a collaborative survey_ ; 2020; EPA 841-R-19-001
  * Xue et al. 2022 Xue, J.; Wang, Q.; Zhang, M. A review of non-point source water pollution modeling for the urban–rural transitional areas of China: Research status and prospect. _Science of The Total Environment_ 2022, 154146
  * Saul et al. 2019 Saul, B. C.; Hudgens, M. G.; Mallin, M. A. Downstream effects of upstream causes. _Journal of the American Statistical Association_ 2019, _114_ , 1493–1504
  * Motew et al. 2017 Motew, M.; Chen, X.; Booth, E. G.; Carpenter, S. R.; Pinkas, P.; Zipper, S. C.; Loheide, S. P.; Donner, S. D.; Tsuruta, K.; Vadas, P. A.; Kucharik, C. J. The influence of legacy P on lake water quality in a Midwestern agricultural watershed. _Ecosystems_ 2017, _20_ , 1468–1482
  * Li et al. 2018 Li, A.; Guo, J.; Li, Z.; Lin, T.; Zhou, S.; He, H.; Ranansinghe, P.; Sturchio, N. C.; Rockne, K. J.; Giesy, J. P. Legacy polychlorinated organic pollutants in the sediment of the Great Lakes. _Journal of Great Lakes Research_ 2018, _44_ , 682–692
  * Sharpley et al. 2013 Sharpley, A.; Jarvie, H. P.; Buda, A.; May, L.; Spears, B.; Kleinman, P. Phosphorus legacy: Overcoming the effects of past management practices to mitigate future water quality impairment. _Journal of environmental quality_ 2013, _42_ , 1308–1326
  * Sharpley et al. 1993 Sharpley, A. N.; Daniel, T.; Edwards, D. Phosphorus movement in the landscape. _Journal of Production Agriculture_ 1993, _6_ , 492–500
  * Van Es et al. 2004 Van Es, H.; Schindelbeck, R.; Jokela, W. Effect of manure application timing, crop, and soil type on phosphorus leaching. _Journal of Environmental Quality_ 2004, _33_ , 1070–1080
  * Zhu et al. 2021 Zhu, L.; Jiang, C.; Panthi, S.; Allard, S. M.; Sapkota, A. R.; Sapkota, A. Impact of high precipitation and temperature events on the distribution of emerging contaminants in surface water in the Mid-Atlantic, United States. _Science of The Total Environment_ 2021, _755_ , 142552
  * Costa et al. 2021 Costa, C. M. d. S. B.; Leite, I. R.; Almeida, A. K.; de Almeida, I. K. Choosing an appropriate water quality model—a review. _Environmental Monitoring and Assessment_ 2021, _193_ , 1–15
  * Lindim et al. 2016 Lindim, C.; Van Gils, J.; Cousins, I. T. A large-scale model for simulating the fate & transport of organic contaminants in river basins. _Chemosphere_ 2016, _144_ , 803–810
  * Mispan et al. 2015 Mispan, M. R.; Abd Rahman, N. F.; Khalid, K.; Haron, S.; Abdul Rasid, M.; Mohd, M. Nutrient Transport Modeling: A Review on models capabilities. _Int. J. Innov. Sci. Eng. Technol_ 2015, _2_ , 908–914
  * Wellen et al. 2015 Wellen, C.; Kamran-Disfani, A.-R.; Arhonditsis, G. B. Evaluation of the current state of distributed watershed nutrient water quality modeling. _Environmental science & technology_ 2015, _49_ , 3278–3290
  * Yuan et al. 2020 Yuan, L.; Sinshaw, T.; Forshay, K. J. Review of watershed-scale water quality and nonpoint source pollution models. _Geosciences_ 2020, _10_ , 25
  * Carpenter and Lathrop 2014 Carpenter, S. R.; Lathrop, R. C. Phosphorus loading, transport and concentrations in a lake chain: a probabilistic model to compare management options. _Aquatic Sciences_ 2014, _76_ , 145–154
  * Cheruvelil et al. 2022 Cheruvelil, K. S.; Webster, K. E.; King, K.; Poisson, A. C.; Wagner, T. Taking a macroscale perspective to improve understanding of shallow lake total phosphorus and chlorophyll a. _Hydrobiologia_ 2022, 1–15
  * Soranno et al. 2015 Soranno, P. A.; Cheruvelil, K. S.; Wagner, T.; Webster, K. E.; Bremigan, M. T. Effects of land use on lake nutrients: The importance of scale, hydrologic connectivity, and region. _PloS one_ 2015, _10_ , e0135454
  * Abed-Elmdoust et al. 2017 Abed-Elmdoust, A.; Singh, A.; Yang, Z.-L. Emergent spectral properties of river network topology: An optimal channel network approach. _Scientific reports_ 2017, _7_ , 1–9
  * Heckmann et al. 2015 Heckmann, T.; Schwanghart, W.; Phillips, J. D. Graph theory—Recent developments of its application in geomorphology. _Geomorphology_ 2015, _243_ , 130–146
  * Schmidt et al. 2020 Schmidt, C.; Kumar, R.; Yang, S.; Büttner, O. Microplastic particle emission from wastewater treatment plant effluents into river networks in Germany: Loads, spatial patterns of concentrations and potential toxicity. _Science of the Total Environment_ 2020, _737_ , 139544
  * Tejedor et al. 2015 Tejedor, A.; Longjas, A.; Zaliapin, I.; Foufoula-Georgiou, E. Delta channel networks: 1. A graph-theoretic approach for studying connectivity and steady state transport on deltaic surfaces. _Water Resources Research_ 2015, _51_ , 3998–4018
  * Tejedor et al. 2015 Tejedor, A.; Longjas, A.; Zaliapin, I.; Foufoula-Georgiou, E. Delta channel networks: 2. Metrics of topologic and dynamic complexity for delta comparison, physical inference, and vulnerability assessment. _Water Resources Research_ 2015, _51_ , 4019–4045
  * Zaliapin et al. 2010 Zaliapin, I.; Foufoula-Georgiou, E.; Ghil, M. Transport on river networks: A dynamic tree approach. _Journal of Geophysical Research: Earth Surface_ 2010, _115_
  * U.S. Environmental Protection Agency 2022 U.S. Environmental Protection Agency, Watershed Index Online. 2022; https://www.epa.gov/wsio, Accessed 09-21-2022
  * U.S. Environmnetal Protection Agency 2014 U.S. Environmnetal Protection Agency, _Wisconsin Integrated Assessment of Watershed Health_ ; 2014; EPA 841-R-14-001
  * Wisconsin Department of Natural Resources 2017 Wisconsin Department of Natural Resources, 24K Hydro Geodatabase. 2017; https://www.arcgis.com/home/item.html?id=cb1c7f75d14f42ee819a46894fd2e771, Accessed 09-21-2022
  * McKay et al. 2012 McKay, L.; Bondelid, T.; Dewald, T.; Johnston, C.; Moore, R.; Rea, A. _NHDPlus Version 2: User Guide_ ; 2012; Accessed on 09-21-2022
  * United States Geological Survey 2022 United States Geological Survey, Watershed Boundary Dataset. 2022; https://apps.nationalmap.gov/downloader/#/, Accessed on 04-20-2022
  * Jordahl et al. 2020 Jordahl, K. et al. geopandas/geopandas: v0.8.1. 2020; https://doi.org/10.5281/zenodo.3946761, Accessed 11-07-2022
  * Hagberg et al. 2008 Hagberg, A. A.; Schult, D. A.; Swart, P. J. Exploring Network Structure, Dynamics, and Function using NetworkX. Proceedings of the 7th Python in Science Conference. Pasadena, CA USA, 2008; pp 11 – 15
  * 55 Wisconsin Department of Natural Resources, Water Quality Data. https://dnr.wi.gov/lakes/waterquality/, Accessed 08-10-2021
  * Dolph et al. 2019 Dolph, C. L.; Boardman, E.; Danesh-Yazdi, M.; Finlay, J. C.; Hansen, A. T.; Baker, A. C.; Dalzell, B. Phosphorus transport in intensively managed watersheds. _Water Resources Research_ 2019, _55_ , 9148–9172
  * Fisher and Acreman 2004 Fisher, J.; Acreman, M. Wetland nutrient removal: a review of the evidence. _Hydrology and Earth system sciences_ 2004, _8_ , 673–685
  * Verhoeven et al. 2006 Verhoeven, J. T.; Arheimer, B.; Yin, C.; Hefting, M. M. Regional and global concerns over wetlands and water quality. _Trends in ecology & evolution_ 2006, _21_ , 96–103
  * Walton et al. 2020 Walton, C. R.; Zak, D.; Audet, J.; Petersen, R. J.; Lange, J.; Oehmke, C.; Wichtmann, W.; Kreyling, J.; Grygoruk, M.; Jabłońska, E., et al. Wetland buffer zones for nitrogen and phosphorus retention: Impacts of soil type, hydrology and vegetation. _Science of the Total Environment_ 2020, _727_ , 138709
  * King et al. 2021 King, K. B.; Wang, Q.; Rodriguez, L. K.; Haite, M.; Danila, L.; Tan, P.-N.; Zhou, J.; Cheruvelil, K. S. LAGOS-US NETWORKS v1.0: Data module of surface water networks characterizing connections among lakes, streams, and rivers in the conterminous U.S. 2021; https://doi.org/10.6073/pasta/98c9f11df55958065985c3e84a4fe995, Accessed 05-02-2022
  * Wang et al. 2018 Wang, X.; Zhang, L.; Zhao, Z.; Cai, Y. Heavy metal pollution in reservoirs in the hilly area of southern China: Distribution, source apportionment and health risk assessment. _Science of the Total Environment_ 2018, _634_ , 158–169
  * Oliver et al. 2019 Oliver, S.; Corburn, J.; Ribeiro, H. Challenges regarding water quality of eutrophic reservoirs in urban landscapes: a mapping literature review. _International Journal of Environmental Research and Public Health_ 2019, _16_ , 40
  * Meinikmann et al. 2015 Meinikmann, K.; Hupfer, M.; Lewandowski, J. Phosphorus in groundwater discharge–A potential source for lake eutrophication. _Journal of Hydrology_ 2015, _524_ , 214–226
  * Valiela et al. 1990 Valiela, I.; Costa, J.; Foreman, K.; Teal, J. M.; Howes, B.; Aubrey, D. Transport of groundwater-borne nutrients from watersheds and their effects on coastal waters. _Biogeochemistry_ 1990, _10_ , 177–197
  * Wang and Baerenklau 2015 Wang, J.; Baerenklau, K. A. How inefficient are nutrient application limits? A dynamic analysis of groundwater nitrate pollution from concentrated animal feeding operations. _Applied Economic Perspectives and Policy_ 2015, _37_ , 130–150
  * 66 Wisconsin Department of Natural Resources, Altoona Lake – Center of Lake Water Quality Data. https://dnr.wi.gov/lakes/waterquality/Station.aspx?id=183082, Accessed 05-16-2022
  * 67 Wisconsin Department of Natural Resources, Lake Mohawksin – Deep Hole Water Quality Data. https://dnr.wi.gov/lakes/waterquality/Station.aspx?id=353089, Accessed 05-16-2022
  * Jones 2010 Jones, N. E. Incorporating lakes within the river discontinuum: longitudinal changes in ecological characteristics in stream–lake networks. _Canadian Journal of Fisheries and Aquatic Sciences_ 2010, _67_ , 1350–1362
  * Burkholder et al. 2007 Burkholder, J.; Libra, B.; Weyer, P.; Heathcote, S.; Kolpin, D.; Thorne, P. S.; Wichman, M. Impacts of waste from concentrated animal feeding operations on water quality. _Environmental health perspectives_ 2007, _115_ , 308–312
  * Long et al. 2018 Long, C. M.; Muenich, R. L.; Kalcic, M. M.; Scavia, D. Use of manure nutrients from concentrated animal feeding operations. _Journal of Great Lakes Research_ 2018, _44_ , 245–252
  * Parry 1998 Parry, R. Agricultural phosphorus and water quality: A US Environmental Protection Agency perspective. _Journal of Environmental Quality_ 1998, _27_ , 258–261
  * Hu et al. 2018 Hu, Y.; Scarborough, M.; Aguirre-Villegas, H.; Larson, R. A.; Noguera, D. R.; Zavala, V. M. A supply chain framework for the analysis of the recovery of biogas and fatty acids from organic waste. _ACS Sustainable Chemistry & Engineering_ 2018, _6_ , 6211–6222
  * Le et al. 2010 Le, C.; Zha, Y.; Li, Y.; Sun, D.; Lu, H.; Yin, B. Eutrophication of lake waters in China: cost, causes, and control. _Environmental management_ 2010, _45_ , 662–668
  * Motew et al. 2019 Motew, M.; Chen, X.; Carpenter, S. R.; Booth, E. G.; Seifert, J.; Qiu, J.; Loheide II, S. P.; Turner, M. G.; Zipper, S. C.; Kucharik, C. J. Comparing the effects of climate and land use on surface water quality using future watershed scenarios. _Science of the Total Environment_ 2019, _693_ , 133484
  * Robertson et al. 2006 Robertson, D. M.; Graczyk, D. J.; Garrison, P. J.; Wang, L.; LaLiberte, G.; Bannerman, R. Nutrient Concentrations and Their Relations to the Biotic Integrity of Wadeable Streams in Wisconsin. 2006, Professional Paper 1722
  * James and Tomer 2020 James, D.; Tomer, M. Agricultural Land Use by Field: Wisconsin 2010-2019. 2020; https://doi.org/10.15482/USDA.ADC/1520625, Ag Data Commons. Accessed 08-26-2021
  * Leavitt et al. 2006 Leavitt, P. R.; Brock, C. S.; Ebel, C.; Patoine, A. Landscape-scale effects of urban nitrogen on a chain of freshwater lakes in central North America. _Limnology and Oceanography_ 2006, _51_ , 2262–2277
  * Soranno et al. 1999 Soranno, P. A.; Webster, K. E.; Riera, J. L.; Kratz, T. K.; Baron, J. S.; Bukaveckas, P. A.; Kling, G. W.; White, D. S.; Caine, N.; Lathrop, R. C., et al. Spatial variation among lakes within landscapes: ecological organization along lake chains. _Ecosystems_ 1999, _2_ , 395–410
  * Brooker et al. 2018 Brooker, M.; Longnecker, K.; Kujawinski, E.; Evert, M.; Mouser, P. Discrete organic phosphorus signatures are evident in pollutant sources within a Lake Erie tributary. _Environmental science & technology_ 2018, _52_ , 6771–6779
  * Makarewicz et al. 2012 Makarewicz, J. C.; Booty, W. G.; Bowen, G. S. Tributary phosphorus loading to Lake Ontario. _Journal of Great Lakes Research_ 2012, _38_ , 14–20
  * Tominac et al. 2020 Tominac, P.; Aguirre-Villegas, H.; Sanford, J.; Larson, R.; Zavala, V. Evaluating landfill diversion strategies for municipal organic waste management using environmental and economic factors. _ACS Sustainable Chemistry & Engineering_ 2020, _9_ , 489–498