# Streamflow Prediction with Uncertainty Quantification for Water Management:
A Constrained Reasoning and Learning Approach

Mohammed Amine Gharsallaoui  Bhupinderjeet Singh  Supriya Savalkar  Aryan
Deshwal  Yan Yan  Ananth Kalyanaraman  Kirti Rajagopalan &Janardhan Rao Doppa
Washington State University, Pullman, WA, USA {m.gharsallaoui,
bhupinderjeet.singh, supriya.savalkar, aryan.deshwal, yan.yan1, ananth,
kirtir<EMAIL_ADDRESS>
###### Abstract

Predicting the spatiotemporal variation in streamflow along with uncertainty
quantification enables decision-making for sustainable management of scarce
water resources. Process-based hydrological models (aka physics-based models)
are based on physical laws, but using simplifying assumptions which can lead
to poor accuracy. Data-driven approaches offer a powerful alternative, but
they require large amount of training data and tend to produce predictions
that are inconsistent with physical laws. This paper studies a constrained
reasoning and learning (CRL) approach where physical laws represented as
logical constraints are integrated as a layer in the deep neural network. To
address small data setting, we develop a theoretically-grounded training
approach to improve the generalization accuracy of deep models. For
uncertainty quantification, we combine the synergistic strengths of Gaussian
processes (GPs) and deep temporal models (i.e., deep models for time-series
forecasting) by passing the learned latent representation as input to a
standard distance-based kernel. Experiments on multiple real-world datasets
demonstrate the effectiveness of both CRL and GP with deep kernel approaches
over strong baseline methods.

##  1 Introduction

Figure 1: Overview of our constrained reasoning and learning approach for
streamflow forecasting to enable water management decisions.

Streamflow is fundamental to meeting societal needs including food and energy
security, and environmental health. In “fully appropriated” Western US
watersheds—where every drop of water is associated with a beneficial use under
a legal water rights system—an accurate representation of the spatio-temporal
variation in streamflow is essential for sustainable management of scarce
water resources. A typical decision context in the snow-dominant Western US is
reservoir operations leading into the spring and summer months when most of
the streamflow is generated (high streamflow months from March through July
which is the focus of our study).

The reservoirs need to be operated to manage competing objectives based on
expectations of upcoming streamflow. Typical objectives include meeting flood
control, irrigation demands for the season, and hydro power generation needs
(see Fig. 1). One concrete example of the conflicting objectives is as
follows. The reservoirs need to be evacuated enough to capture upcoming
streamflow and avoid flood risk. However, this cannot be at the expense of not
being able to fill the reservoir before the irrigation season starts. Having a
good forecast of expected flows in high flow months, and reliable uncertainty
bounds around it can help optimize storage and release decisions to meet all
objectives. Moreover, different decisions will involve utilizing uncertainty
bounds in different ways. When flood risk needs to be minimized, the reservoir
manager will likely utilize flows at the upper bound of the uncertainty
interval in water release decisions. In contrast, while focusing on meeting
irrigators’ needs, the reservoir manager will likely make storage decisions
based on the lower bound of expected streamflow. Under the auspices of a
federal grant, this project team has been interacting with United States
Bureau of Reclamation (USBR)—the federal agency tasked with managing reservoir
operations for all reservoirs in the Western United States where water for
irrigation is a key operational priority along with a multitude of other uses
(see Fig. 1). Therefore, this work has a real potential for impact at the U.S.
national scale.

Process-based hydrological models (aka physics-based models) that take
weather, soil, and other inputs to translate precipitation to streamflow using
biophysical representations of the hydrological process are commonly used to
estimate streamflow. These models have been known to yield best results when
calibrated to individual watersheds but they continue to struggle in regional-
scale modeling across multiple watersheds because hydrological response is
heterogeneous and watershed-specific Kratzert et al. (2019).

Data-driven (aka deep learning) models can capture this heterogeneity and push
the frontiers in learning universal, regional, and local hydrological behavior
Kratzert et al. (2019); Peters-Lidard et al. (2017); Karpatne et al. (2016)
while increasing streamflow prediction accuracy. However, this learning
methodology has to maintain the basic laws governing hydrological processes
(e.g., conservation of mass for water balance) and capture the essence of
important intermediary processes (e.g., evapotranspiration) that translate
precipitation to streamflow. Since regional-scale water management decisions
are made on “seasonal” time scales, we seek to make streamflow prediction for
high-flow months to allow a targeted impact. As a result, we lack large
amounts of training data required by the deep learning methods—e.g., even with
25 years of monthly data, there are only 125 samples.

Furthermore, in watersheds with competing water uses, when streamflow
predictions go wrong, their societal impacts are significant—including losses
for agriculture and the rural economy given less-than-ideal drought
management; and negative impacts on endangered fish species due to
higher/lower flows than expected, e.g., the 2015 drought year resulted in
$\sim$$0.7 billion agricultural revenue loss in the U.S. Washington State
McLain et al. (2015). Therefore, quantifying uncertainty under which water
managers make these decisions is critical for public trust, conflict
resolution, and informed decision making.

Challenges. There are three key technical challenges: 1) Predicting outputs
which are consistent with the physical laws; 2) Training predictive models
from small data which generalize well to unseen data; and 3) Quantifying the
uncertainty in terms of a prediction interval such that the ground-truth
output lies in this interval in most cases.

Contributions. This paper’s key contribution is a principled approach for
streamflow prediction and uncertainty quantification to address the above-
mentioned three challenges.

Prior work has attempted to address the first challenge using physics-guided
training or physics-informed neural networks Karpatne et al. (2017a); Cai et
al. (2021) where the key idea is to penalize the deep model when the
predictions are not consistent with the physical laws. This approach was shown
to improve over the pure data-driven approach on predicting depth-specific
lake water temperatures Jia et al. (2021, 2019); Karpatne et al. (2017b) and
other use-cases Muralidhar et al. (2020). However, this approach and similar
methods Kamruzzaman et al. (2021); Hoedt et al. (2021); Sit et al. (2021);
Xiang and Demir (2022) do not guarantee that every test-time prediction is
consistent with the domain knowledge. To overcome this drawback, as our first
contribution, _we study a new constrained reasoning and learning (CRL)
approach using a semantic probabilistic layer (SPL)_ Ahmed et al. (2022b). SPL
compiles the domain knowledge specified in the form of logical constraints
into a tractable circuit representation Vergari et al. (2020) and adds it as a
replacement for predictive layers in deep networks for efficient end-to-end
training. We empirically show that CRL approach with SPL instantiation uses
the domain knowledge to improve the predictive accuracy over baseline methods.

To address the small data challenge, _we propose a novel importance-weighted
(IW) training approach to improve the generalization accuracy of predictive
models learned over the latent representation from the deep temporal model_.

The key principle is to assign importance weights directly proportional to the
errors on data samples during the training process, and to compute the
weighted gradients using the data samples for gradient-based training. We
provide theoretical analysis to justify the IW training approach. Our
experiments show that IW training further improves the predictive accuracy
over the standard training. Our IW algorithm and its theoretical analysis is
specific to the small data challenge unlike prior work Xu et al. (2021).

For uncertainty quantification (UQ), we study a principled approach by
combining the synergistic strengths of Gaussian processes (GPs) and deep
temporal models. Deep neural networks are known to extrapolate over-
confidently Valentin Jospin et al. (2020). On the other hand, GP models are
gold standard Williams and Rasmussen (2006) for UQ providing principled
Bayesian uncertainty estimates, unlike existing ad-hoc techniques Liu et al.
(2023), but have difficulty scaling to high-dimensional structured inputs
(e.g., time-series). Deep kernel learning (DKL) Wilson et al. (2016) is a
recent promising technique towards handling this drawback while maintaining
the good UQ properties of a GP. Our experiments demonstrate that GP based
approach with deep kernel produces prediction intervals with high coverage
(i.e., ground truth output lies within the interval). We make our code and
datasets available to the AI community at
https://github.com/aminegha/StreamPred to catalyze research in this important
problem space.

##  2 Background and Problem Setup

In streamflow prediction, utilizing graph representations is an efficient
method for illustrating watersheds. Within a defined geographical expanse, a
watershed represents the entirety of land that contributes water to a shared
outlet, which might be a river, lake, or ocean. Within the watershed, a river
network showcases the flow paths and connections of water in that region. A
watershed is subdivided into grids, representing discrete small sections that
facilitate the analysis and modeling of factors such as water movement. In the
graph representation, the _nodes_ represent the grids within the watershed,
and the _edges_ represent the flows between them. Each node has a set of
_environmental_ features (including variables related to precipitation and
temperature). Given this rich graph representation, our objective is to
predict the streamflow at the watershed outlet at fixed time intervals. Since
water management decisions are made seasonally, we aim to make streamflow
predictions for high-flow months (March to July).

Domain Knowledge. The water balance equation (WBE) is a fundamental concept in
hydrology: the amount of water in a system should be balanced between the
inputs (e.g., precipitation), outputs (e.g., evapotranspiration, streamflow),
and changes in storage (e.g., soil moisture or snowpack). It is given by
$P\geq ET+S$, where $P$ is precipitation, $ET$ is evapotranspiration, and $S$
is streamflow. WBE is not an equality because all components of the water
balance are not present here and holds only at an annual aggregated level.

Problem Setup. For a given watershed, we consider a sequence of graphs over
discrete time steps denoted as $G_{1},G_{2},\cdots,G_{T}$, where each graph
$G_{t}$ represents a directed graph with a set of nodes $V$ and a set of edges
$E\subseteq V\times V$. At each time step, the corresponding features matrix
$X_{t}\in\mathbb{R}^{n\times m}$, where $n=|V|$ and $m$ is the number of
features. The graph $G_{t}$ is also characterized by an adjacency matrix
$A\in\mathbb{R}^{n\times n}$, with the property that $A_{t}$ remains invariant
since grid connections are static. At each (monthly) time step $t$, we also
have the target variable $y_{t}$ which is the streamflow at the watershed
_outlet_ node. Given a training data of time-series
$\mathcal{D}=\\{(X_{t},y_{t})\\}_{t=1}^{T}$ and domain knowledge $\mathcal{K}$
in the form of water balance equation, we have two goals:

  1. a)

To learn a predictor $F$ that can make accurate predictions consistent with
the knowledge $\mathcal{K}$ for unseen inputs $X_{t^{\prime}}$ where
$t^{\prime}>T$ (high-flow months only).

  2. b)

To get uncertainty estimates in the form of a prediction interval $[lb,ub]$
for $X_{t^{\prime}}$, where the true output $y_{t^{\prime}}\in[lb,ub]$ with a
high probability.

##  3 Technical Approach

In this section, we outline our overall methodology shown in Algorithm 1 to
address the three technical challenges, namely, predictions consistent with
domain knowledge, small training data, and uncertainty quantification.

###  3.1 Constrained Reasoning and Learning

Recurrent Convolutional Graph Neural Network. We employ a recurrent
convolutional graph convolutional neural network (RCGNN) Seo et al. (2016);
Sit et al. (2021); Xiang and Demir (2022) to solve the streamflow prediction
task. RCGNN is a generalization of classical graph convolutional network that
incorporates the recurrent nature of RNN to capture temporal dependencies in
graph-structured data. Notably, in graph embedding learning, for a given node,
we integrate information from both preceding time steps and adjacent river
segments (i.e., edges).

The RCGNN model establishes transition connections for the derived latent
representation using a recurrent cell configuration. In the scope of this
study, we employ the Gated Recurrent Unit (GRU) cell, known for its prowess in
capturing sequential patterns and dependencies in graph-structured data. For
time step $t$, given a graph represented by its adjacency matrix $A$ and its
input feature matrix $X$, a graph convolutional GRU layer is parameterized
with multiple weight matrices $W$’s, $U$’s and bias terms $b$’s, and starts by
computing the reset gate $r_{t}$ and updates gate $z_{t}$ from time step
$t-1$:

$\displaystyle r_{t}=\sigma(W_{r}\cdot X+U_{r}\cdot(A\cdot H^{(t-1)})+b_{r})$
$\displaystyle z_{t}=\sigma(W_{z}\cdot X+U_{z}\cdot(A\cdot H^{(t-1)})+b_{z}),$

where $H^{(t-1)}$ is the hidden state of the previous time step $t$, while
$\sigma$ denotes the sigmoid activation function. Next, it computes candidate
hidden state ($h_{t}^{\prime}$):

$\displaystyle h_{t}^{\prime}=\tanh(W_{h}\cdot X+U_{h}\cdot(r_{t}\odot(A\cdot
H^{(t-1)}))+b_{h}).$

Here, $\odot$ represents element-wise multiplication. Using $z_{t}$ and
$h_{t}^{\prime}$, we update the current hidden state $H^{(t)}$ as follows:

$\displaystyle H^{(t)}=z_{t}\odot H^{(t-1)}+(1-z_{t})\odot h_{t}^{\prime}.$

Therefore, the learnable weight matrices of the model are specifically denoted
above by $W_{r}\text{, }U_{r}\text{, }W_{z}\text{, }U_{z}\text{, }W_{h}\text{
and }U_{h}$, and trainable bias terms are $b_{r}\text{, }b_{z}\text{ and
}b_{h}$. For the subsequent output layers, two distinct activation functions
are employed alongside two linear layers. For the activation function, we use
LeakyRelu to introduce non-linearity. We build our RCGNN model to learn the
graph embedding according to the above architecture, denoted by $f(X)$.

Semantic Probabilistic Layer. To address the challenge of maintaining domain
knowledge consistency, the constrained reasoning and learning (CRL) approach
instantiates the concept of semantic probabilistic layer (SPL) Ahmed et al.
(2022b). SPL can be replaced with the predictive layer of any deep neural
network. It allows to model domain knowledge as logical constraints over
multiple inter-related output variables and guarantees consistency of
predicted outputs.

Given a structured input $X$, the computation of the probability of a
candidate structured output $Y$ is decomposed as
$p(Y|f(X))=q_{\Theta}(Y|f(X))\cdot c_{\mathcal{K}}(X,Y)/\mathcal{Z}(X)$, where
$f(X)$ represents the feature embedding or latent representation of input $X$;
$q_{\Theta}(Y|f(X))$ parameterized by $\Theta$ allows us to perform
probabilistic reasoning over candidate structured outputs $Y$;
$c_{\mathcal{K}}(X,Y)$ is a constraint reasoner to ensure that predicted
structured outputs $Y$ are consistent with the domain knowledge $\mathcal{K}$
and produces evaluation $1$ (i.e., $c_{\mathcal{K}}(X,Y)$=1) only when $Y$
satisfies the declarative knowledge $\mathcal{K}$; and $\mathcal{Z}(x)$ is a
normalization constant. SPL leverages the beneficial properties of
probabilistic circuits Vergari et al. (2020); Choi et al. (2020) for tractable
learning and reasoning. It encodes $q_{\Theta}(Y|f(X))$ and
$c_{\mathcal{K}}(X,Y)$ as circuits, and considers constrained computational
graphs to achieve tractability. Therefore, SPL allows us to produce consistent
and correct structured outputs for every inference computation task, i.e.,
satisfies the domain knowledge in the form of water balance equation. We refer
the reader to SPL paper Ahmed et al. (2022b) for additional details.

Algorithm 1 Overall Training Algorithm

Input: Training data of graph time-series
$\mathcal{D}=\\{(X_{t},y_{t})\\}_{t=1}^{T}$; domain knowledge $\mathcal{K}$;
Recurrent GNN model $f_{\theta_{fe}}$; Predictive model $F_{\theta_{p}}$; GP
kernel $\kappa$
Output: Trained feature extractor $f_{\theta_{fe}}$, predictive model
$F_{\theta_{p}}$, and GP model with hyper-parameters of $\kappa$

1: Train the recurrent GNN model $f_{\theta_{fe}}$ on $\mathcal{D}$ with SPL
layer for domain knowledge $\mathcal{K}$ to create latent representations
$z_{i}$=$f(X_{i})$ // Section 3.1

2: Train the predictive regression model $F_{\theta_{p}}$ using
$\mathcal{D}_{\text{tr}}=\\{z_{i},y_{i}\\}_{i=1}^{T}$ with $T$ regression
examples using the importance weighted training approach // Section 3.2

3: Train the GP model on $\mathcal{D}_{\text{tr}}=\\{z_{i},y_{i}\\}_{i=1}^{T}$
to optimize the hyper-parameters of the kernel $\kappa$ // Section 3.3

4: return feature extractor, predictor, and GP model

###  3.2 Importance Weighted Training

Given the learned latent representations $\\{z_{i}=f(X_{i})\\}_{i=1}^{T}$ from
the recurrent GNN model, we describe a novel importance weighted training
approach to improve the generalization accuracy of the predictive model $F$.

Notations. Our goal is to train regression models over the training data
$\mathcal{D}_{\text{tr}}=\\{z_{i},y_{i}\\}_{i=1}^{T}$ with $T$ examples.
Suppose each data sample $(z_{i},y_{i})$ is drawn from a target distribution
$\mathcal{P}$ over the space $\mathcal{Z}\times\mathcal{Y}$ such that $z_{i}$
is an input from the input space $\mathcal{Z}$ and $y_{i}\in\mathcal{Y}$ is
the corresponding ground-truth output. Let
$F:\mathcal{Z}\rightarrow\mathcal{Y}$ denote a predictive model trained on the
training set $\mathcal{D}_{\text{tr}}$. A loss function
$\ell:\mathcal{Y}\times\mathcal{Y}\rightarrow\mathbb{R}$ is used to measure
the accuracy of predictions made by $F$. In this paper, we consider the
residual $\ell(F(z),y)=|F(z)-y|$ as loss and assume that there exists an upper
bound of loss as
$L=\max_{(z,y)\in\mathcal{Z}\times\mathcal{Y}}\ell(F(x),y)<\infty$. Typically,
$F$ is expected to reach a small true risk, defined as
$R(F)=\mathbb{E}_{(z,y)\sim\mathcal{P}}[\ell(F(z),y)]$ in the sense of
population. Unfortunately, we cannot measure $R(F)$ due to the availability of
only finite training samples drawn from $\mathcal{P}$.

Hence, the empirical risk
$\widehat{R}(F)=\frac{1}{T}\sum_{i=1}^{T}\ell(F(z_{i}),y_{i})$ on
$\mathcal{D}_{\text{tr}}$ is employed to estimate it, leaving a standard
estimation error bound in $O(1/\sqrt{T})$ (see Lemma 1). The rate of
convergence of generalization error bound is particularly important when the
number of training samples $T$ is small. For example, to ensure the same
$O(10^{-2})$ generalization error bound, $O(1/\sqrt{T})$ requires
$T=O(10^{4})$, while the faster rate $O(1/T)$ only needs $T=10^{2}$,
significantly reducing the requirement of large training sets from modern ML
methods.

High-level Algorithm. Importance-weighted (IW) training is an iterative
gradient-based approach. It assigns importance weight $\omega(z,y)$ to each
sample $(z,y)\in\mathcal{D}_{\text{tr}}$ such that $\omega(z,y)$ is directly
proportional to the loss $\ell(F(z),y)$ based on the current model $F$ and
employs the weighted gradient to update the model parameters of $F$ in each
iteration. Analogous to the empirical risk $\widehat{R}(F)$, we define the IW
empirical risk with weights $\omega(z,y)$ for $F$ as

$\displaystyle\widehat{R}_{\omega}(F)=\frac{1}{T}\sum_{i=1}^{T}\omega(z_{i},y_{i})\cdot\ell(F(z_{i}),y_{i}).$
(1)

$\widehat{R}_{\omega}(F)$ is general enough to include the empirical risk
$R(F)$ as a special case: setting $\omega(z_{i},y_{i})=1$ for all data samples
reduces $\widehat{R}_{\omega}(F)$ to $\widehat{R}(F)$. In what follows, we
show that under some configurations of importance weighting function $\omega$,
$\widehat{R}_{\omega}(F)$ gives a tighter estimation of the true risk $R(F)$.

Theoretical Analysis. The following lemma shows the generalization error bound
for IW empirical risk $R_{\omega}(F)$.

######  Lemma 1.

(Error bound of IW risk, Theorem 1 in Cortes et al. (2010)) Let
$M=\sup_{z,y}\omega(z,y)$ denote the infinity norm of $\omega$ on the domain.
For given $F$ and $\delta>0$, with probability at least $1-\delta$, the
following bound holds:

$R(F)-\widehat{R}_{\omega}(F)\leq\frac{2M\log(1/\delta)}{3T}+\sqrt{\frac{2d_{2}(\mathcal{P}||\frac{\mathcal{P}}{\omega})\log(1/\delta)}{T}},$
(2)

where $d_{2}(\mathcal{P}\left|\right|\mathcal{Q})$ =
$\int_{x}\mathcal{P}(x)\cdot\frac{\mathcal{P}(x)}{\mathcal{Q}(x)}dx$ is the
base-$2$ exponential for Rényi divergence of order $2$ between two
distributions $\mathcal{P}$ and $\mathcal{Q}$ and $T$ is the number of
training samples.

In the above prior result, the last term on the right side
($O(\sqrt{d_{2}/T})$) typically dominates the first term ($O(M/T)$), i.e.,
$O(M/T)\ll O(\sqrt{d_{2}/T})$. Consequently, the choice of $\omega$ mainly
changes the quantity of Rényi divergence $d_{2}$ and thus the dominating term
$O(\sqrt{d_{2}/T})$. In what follows, we theoretically show that a simple
principle to choose group-wise weights $\omega$ can easily shrink $d_{2}$ in
the dominating term compared with the conventional empirical risk.
Specifically, the grouping principle is based on the prediction loss
$\ell(F(z),y)$ of $F$ overall data samples. Given a total number of groups
$K$, denote the group $k$ by $\mathcal{D}_{k}$ such that any data from
$\mathcal{D}_{k}$ suffers smaller loss than that from the next group
$\mathcal{D}_{k+1}$, i.e., $\ell(F(z),y)\leq\ell(F(z^{\prime}),y^{\prime})$
for all $(z,y)\in\mathcal{D}_{k}$ and
$(z^{\prime},y^{\prime})\in\mathcal{D}_{k+1}$ with $k\in\\{1,...,K-1\\}$ (see
discussion in “Practical Algorithm”). Accordingly, the probability of drawing
data from $\mathcal{D}_{k}$ is denoted by
$\mathcal{P}_{k}=\mathbb{P}\\{(z,y)\in\mathcal{D}_{k}\\}$. Specifically, for
the fixed $F$ and the same analysis framework under Lemma 1, let $\widehat{B}$
and $\widehat{B}_{\omega}$ be the generalization error bounds of empirical
risk and IW empirical risk, respectively, i.e.,
$R(F)-\widehat{R}(F)\leq\widehat{B}$ and
$R(F)-\widehat{R}_{\omega}(F)\leq\widehat{B}_{\omega}$. Below we prove that
$\widehat{B}_{\omega}<\widehat{B}$ under a certain group-wise setting for
$\omega$. Before proceeding to the main theorem, we define group-wise weights
as $\\{\omega_{k}\\}_{k=1}^{K}$, where $\omega_{k}$ is assigned to data from
group $\mathcal{D}_{k}$ as the weight shown by $\widehat{R}_{\omega}(F)$ in
(1).

######  Theorem 1.

(IW Improves generalization bound) Assume
$\frac{2M\log(1/\delta)}{3T}\leq\sqrt{\frac{2d_{2}(\mathcal{P}||\frac{\mathcal{P}}{\omega})\log(1/\delta)}{T}}$
in (2) of Lemma 1. Set $\omega(z,y)=\omega_{k}$ for $(z,y)\in\mathcal{D}_{k}$
such that $\mathcal{P}_{k}\cdot\omega_{k}=\frac{k^{a}}{K^{b+1}}$ for any
$k\in\\{1,...,K\\}$. Under the conditions $0<a\leq\min\\{K-1,b-1/\ln(K)\\}$,
and $K\geq 2$, the IW empirical risk $\widehat{R}_{\omega}(F)$ gives a tighter
generalization bound than empirical risk $\widehat{R}(F)$ with
$\widehat{B}_{\omega}=\sqrt{\frac{2\log(1/\delta)}{(a+1)T}}$ and
$\widehat{B}=\sqrt{\frac{2\log(1/\delta)}{T}}$.

The complete proof of Theorem 1 is in the Appendix.

Remark. Theorem 1 shows some mild conditions to achieve a tighter
generalization bound by strategically setting the group-wise weights
$\omega_{k}$ in IW method than empirical risk, i.e.,
$\widehat{B}_{\omega}<\widehat{B}$. Furthermore, it demonstrates how fast the
hyper-parameter $a$ reduces $\widehat{B}_{\omega}$, and thus builds the
connection between $a$ and the tightness of the bound. For example, setting
$a=10^{2}$ makes $\widehat{B}_{\omega}$ $10$ times smaller than $\widehat{B}$.

Practical Algorithm. To devise a practical IW algorithm, we need to
appropriately set the values of $K$ (no. of groups) and $\omega_{k}$ (weights
for each of the $k$ groups). The key idea of setting the group-wise weight is
to make $\omega_{k}\propto k/\mathcal{P}_{k}$ when we group data into
$\\{\mathcal{D}_{k}\\}_{k=1}^{K}$ such that all data in $\mathcal{D}_{k}$
suffers smaller loss than that the data in $\mathcal{D}_{k+1}$. We propose
range-wise weighting, i.e., grouping data into several bins (loss intervals)
according to their loss values. For example, for $K=2$, a threshold $\tau$
partitions all data into two bins, i.e., the group suffering smaller loss
$\mathcal{D}_{1}=\\{(z,y):\ell(z,y)\leq\tau\\}$ and the other one suffering
larger loss $\mathcal{D}_{2}=\\{(z,y):\ell(z,y)>\tau\\}$, with
$\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ denoting the percentage of data in the
two groups. This strategy can be easily generalized to $K>2$ by using multiple
thresholds, such as uniformly grouping data into $K$ ranges with interval
$\frac{L}{K}$ for maximal loss $L$:
$\mathcal{D}_{k}=\\{(z,y):\frac{L(k-1)}{K}\leq\ell(F(z),y)\leq\frac{Lk}{K}\\}$
for $k\in\\{1,\cdots,K\\}$. In practice, typically more data samples suffer
from smaller training loss, i.e., concentrating around the true output, which
implies $\mathcal{D}_{k}$ with smaller group ID $k$ leads to larger
$\mathcal{P}_{k}$. Hence, setting $\omega_{k}\propto k/\mathcal{P}_{k}$
indicates that group $\mathcal{D}_{k}$ with smaller group ID $k$ gets smaller
weight $\omega_{k}$, In contrast, larger weights $\omega_{k}$ are encouraged
for $\mathcal{D}_{k}$ with larger group ID $k$. Training model with this IW
configuration improves the generalization performance of $F$ according to
Theorem 1.

In our implementation, we employed this strategy for assigning IWs with
$K$=10.

###  3.3 Gaussian Process with Deep Kernels for UQ

Gaussian processes (GPs) Williams and Rasmussen (2006) are considered as gold
standard for uncertainty quantification (UQ). GPs are non-parametric models
and allows the practitioner to incorporate domain knowledge in the form of
kernels $\kappa(x,x^{\prime})$ to measure similarity between pairs of inputs
$x$ and $x^{\prime}$ (e.g., RBF kernel).

GPs are Bayesian models whose posterior mean and standard deviation provide
prediction and uncertainty for each input. A major drawback of GPs is that
they have difficulty scaling to high-dimensional structured inputs (e.g.,
time-series of graphs) as in our problem setting. On the other hand, deep
kernel learning (DKL) Wilson et al. (2016) is a recent technique that has
shown to be a promising solution towards handling this drawback. To
synergistically combine the strengths of GPs and DKL, we parameterize the
kernel function of the GP model with a neural network feature extractor that
is passed to a canonical distance-aware kernel such as RBF kernel.
Specifically, we add a GP layer on top of the latent representations learned
by the recurrent GNN model. Given the learned latent representations
$z_{i}=f(X_{i})$ for all $i$=1 to $T$, we create a GP using
$\mathcal{D}_{\text{tr}}=\\{z_{i},y_{i}\\}_{i=1}^{T}$ with $T$ training
examples. We employ the latent representations $z_{i}=f(X_{i})$ to define a
RBF kernel as shown below.

$k(z,z^{\prime})=\sigma_{f}e^{-||z-z^{\prime}||_{2}^{2}/(2\gamma^{2})}.$ (3)

Training GPs refers to estimating the hyper-parameters of the GP kernel. For
example, the RBF kernel in equation (3) has the length-scale $\gamma$ and the
signal variance $\sigma_{f}$ hyperparameters. To learn the hyperparameters of
the kernel, we maximize the marginal likelihood of the observed data using the
training data $\mathcal{D}_{\text{tr}}$. Another advantage of GPs for UQ is
that we can update the model online based on new training examples to improve
the uncertainty estimates for future predictions.

The predictive uncertainty estimate for an input with latent representation
$z$ can be computed in closed form (for the noiseless observations setting) as
given below:

$\displaystyle\sigma^{2}(z)=k(z,z)-k_{*}^{T}\mathbf{K}^{-1}k_{*}$ (4)

where $k_{*}=[k(z,z_{1}),k(z,z_{2})\cdots k(z,z_{T})]$ and $\mathbf{K}$ is the
kernel matrix computed over the latent representation of training inputs. We
can compute a GP driven confidence interval for input $z$ as
[$\mu(z)-2\sigma(z)$, $\mu(z)+2\sigma(z)$], where $\mu(z)$ is the posterior
mean prediction of the GP.

##  4 Experiments and Results

In this section, we describe our empirical results and analysis.

###  4.1 Experimental Setup

Watershed | BR | C | CCR | F | SFC | Y
---|---|---|---|---|---|---
#grids | 57 | 437 | 103 | 100 | 97 | 449
Table 1: Description of the watershed datasets. All watersheds were observed
at a monthly rate from Oct. 1979 through Sep. 2014 (i.e., 420 months).
Temporal features for each grid (graph node) includes incoming precipitation,
max and min temperature, average temperature, wind speed, specific humidity,
solar radiation, max and min relative humidity, and evapotranspiration. The
environmental variables are obtained from a 4-km gridded dataset, regridded to
the VIC-CropSyst model’s resolution of 6-km via linear interpolation.

Watershed Datasets. We collected US Geological Survey’s observed monthly
streamflow at the outlet of six watersheds in the Columbia River basin of the
Pacific Northwest US : Boise River (BR), Clearwater (C), Clearwater Canyon
Ranger (CCR), Flathead (F), South Fork Clearwater (SFC), and Yakima (Y). Our
intent was to primarily focus on ”natural” watersheds with minimal human
influence, and this criterion along with a constraint on the watershed size
(minimum of 1800 sq. km) resulted in the first five watersheds. Additionally,
we included the Yakima River basin as an example of human-influenced watershed
knowing that the model performance will likely be poor. Table 1 provides a
high-level summary. Each watershed has different number of grids (no. of nodes
in graph) and river network corresponds to the adjacency matrix of the graph.
Each graph node has multiple time-varying (monthly) environmental features.

The watersheds encompass a span of water years, commonly used in hydrology to
capture the hydrological cycles, starting in October and ending in September
next year. We employed the following training/validation/testing splits. The
first 20 years of data was used for training and the next 6 years for hyper-
parameter tuning. We employ the last 9 years of data for testing the
generalization of various methods.

Configuration of Predictive Algorithms. We compare our proposed approach with
several baseline methods.

1) Physics-based model (VIC-CropSyst): We use VIC-CropSyst Malek et al.
(2017); Rajagopalan et al. (2018) which comprises two coupled models—the VIC
model, which is a semi-distributed macroscale hydrological model Liang et al.
(1994), and the CropSyst model Stöckle et al. (2003). For this study, the
simulations were performed with naturalized water settings—i.e., without
extractions of water from the streams such as irrigation.

2) Recurrent Convolutional GNN (RCGNN): We train the RCGNN model (Section 3.1)
using PyTorch to optimize the Mean Squared Error (MSE) loss. It is a pure
data-driven model without the knowledge of water balance equation.

3) Physics-guided RCGNN (PG-RCGNN): This method incorporates water balance
equation (WBE) into the training of RCGNN following approach by Karpatne et
al. (2017a). It employs the combined loss
$\text{loss}_{mse}+\lambda_{wbe}\cdot\text{loss}_{wbe}$, where the
$\text{loss}_{wbe}=max(0,\hat{y}+ET-P)$ term penalizes predictions deviating
from the water balance equation.

4) RCGNN with SPL layer (RCGNN-CRL): This approach employs constrained
reasoning and learning (CRL) by adding a semantic probabilistic layer to
incorporate the water balance equation in yearly aggregates. We use the Pylon
software library Ahmed et al. (2022a) for this implementation.

5) RCGNN-CRL with Importance Weighted Training (RCGNN-CRL-IW): This is a
variant of RCGNN-CRL which uses importance weighted training described in
Section 3.2. We employ $K$=10 groups and the specific weighting strategy
detailed as part of the practical algorithm.

Please see the appendix section D for hyper-parameter tuning.

Configuration of UQ Algorithms. We compare our proposed Gaussian process
approach with a dropout baseline.

1) Gaussian process w/ deep kernel learning (GP w/ DKL). We define a GP model
by defining a RBF kernel over the latent representations
$\\{z_{i}$=$f(X_{i})\\}_{i=1}^{T}$ learned by the recurrent GNN model. We find
the hyper-parameters of the RBF kernel by optimizing the negative log-
likelihood of data.

2) Dropout method: We introduce two dropout layers after obtaining the latent
representation to inject stochasticity into inference Hinton et al. (2012). To
obtain output variance $\sigma$, we generate distinct sets of dropped neurons
within each iteration, to capture the intrinsic uncertainty of deep model. The
prediction interval is [$\hat{y}-2\sigma$, $\hat{y}+2\sigma$] where $\sigma$
is calculated over $30$ iterations with a dropout rate of 0.2.

Evaluation Methodology. For comparing the predictive accuracy of different
methods, we employ two different metrics: Normalized Nash–Sutcliffe Efficiency
(NNSE) Nossent and Bauwens (2012) and Mean Absolute Error (MAE). NSE is a
widely recognized metric to assess the precision of hydrological models Nash
and Sutcliffe (1970) and is calculated as:

$\displaystyle
NSE=1-\frac{\sum_{i=1}^{n}(y_{i}-\hat{y}_{i})^{2}}{\sum_{i=1}^{n}(y_{i}-\bar{y})^{2}}$
(5)

Here, $n$ is the number of distinct year-months during testing for the
streamflow prediction, $y_{i}$ is the observed value at month $i$,
$\hat{y}_{i}$ is the predicted value at month $i$, and $\bar{y}$ is the mean
of the observed values. Since the theoretical lower limit of the NSE
$(-\infty)$ can lead to issues in implementations, we normalize the NSE within
the range of [0, 1]. The Normalized NSE is given by $NNSE$=$\frac{1}{2-NSE}$.
NNSE takes into account both accuracy of predictions as well as the spread of
predicted values around ground truth. Hence, a higher NNSE value indicates
better predictive performance. For comparing the performance of UQ, we employ
two metrics: prediction interval (PI) size and marginal coverage (whether
ground truth lies in the PI or not) and . There is an inherent trade-off
between coverage and PI size. We can get large coverage by producing large PIs
and small PIs which achieve low coverage. Ideally, we prefer to achieve some
minimum target coverage (e.g., 90%) and smaller prediction intervals. We
report all results with a single training run.

###  4.2 Results and Discussion

Performance of Streamflow Predictive Models. Table 2 shows the accuracy
results comparing different streamflow prediction models across six watersheds
in the high-flow months. We make the following observations. 1) All data-
driven models perform better than physics-based model VIC-CropSyst across all
watersheds. 2) All methods that combine data and knowledge perform better than
the pure data-driven model RCGNN for all cases. This demonstrates the value of
domain knowledge in the form of the water balance equation 3) Importance-
weighted training improves the predictive accuracy of RCGNN-CRL-IW over the
standard training. This result demonstrates the practical value of our
theoretically-justified IW training approach. 4) Accuracy enhancement,
relative to the physics-based model, varies by watershed, and time of the year
(Fig 2, Fig 5 and Fig 6 in Appendix).

|  VIC- CropSyst |  RCGNN |  PG- RCGNN |  RCGNN- CRL |  RCGNN- CRL-IW
---|---|---|---|---|---
BR |  0.823 |  0.915 |  0.919 |  0.93 |  0.934
C |  0.837 |  0.902 |  0.925 |  0.948 |  0.952
CCR |  0.817 |  0.883 |  0.896 |  0.931 |  0.935
F |  0.684 |  0.885 |  0.902 |  0.917 |  0.936
SFC |  0.741 |  0.834 |  0.866 |  0.91 |  0.92
Y |  0.212 |  0.405 |  0.509 |  0.635 |  0.654
Table 2: Comparison of predictive models across watersheds during high flow
months. Values are reported in terms of the NNSE metric (higher the better).
NNSE values less than 0.5 indicate that the mean is a better predictor than
the model’s output.

Table 3 illustrates the frequency and magnitude of water balance equation
violations by PG-RCGNN and RCGNN-CRL-IW across various datasets during test
years. Reported magnitudes for PG-RCGNN vary between datasets, reflecting the
disparity in violation severity. Only years with constraint violations are
included in the magnitude calculations, for a fair comparison. Conversely,
RCGNN-CRL-IW consistently maintains the water balance equation across all test
years and datasets. This result demonstrates the effectiveness of RCGNN-CRL-IW
in always obeying the water balance equation during inference.

|  BR |  C |  CCR |  F |  SF |  Y
---|---|---|---|---|---|---
PG- RCGNN |  0.44 (51.41) |  0.44 (25.93) |  0.44 (106.1) |  0.67 (52.26) |  0.11 (8.09) |  0.0 (0.0)
RCGNN- CRL-IW |  0.0 (0.0) |  0.0 (0.0) |  0.0 (0.0) |  0.0 (0.0) |  0.0 (0.0) |  0.0 (0.0)
Table 3: Comparison of fraction (magnitude) of constraint violation for PG-
RCGNN and RCGNN-CRL-IW. Fraction is the proportion of test years violating the
constraint, while magnitude is the average yearly constraint violation
measured in millimeters. Figure 2: Relative performance chart comparing all
models for high flow months (a) Boise, (b) Clearwater, (c) Clearwater Canyon
Ranger, (d) Flathead, (e) Southfork Clearwater, and (f) Yakima. For each test
month and year, each model’s MAE was compared against the best performing
model’s MAE for that month and year (difference shown in X-axis). The fraction
of test month-years for which that performance is achieved is shown along the
Y-axis. Closer a curve is aligned with the Y-axis for larger fraction, the
better.

Fig 2 shows the relative performance of each model compared to other models.
RCGNN-CRL-IW (red line) yields the best performance most consistently for the
high flow months. Even, when it is not the best model, its degradation in MAE
relative to the best model is low (i.e., the red line sticks closer to the
Y-axis. On the other hand, the VIC-CropSyst and RCGNN models are best
performing in very few cases, and when they are not the best models, the
degradation in performance is relatively high (the lines are further away from
the Y-axis). This is primarily due to the failure to capture the timing and
magnitude of peak streamflow (Fig 5 ).

Performance of UQ Methods. Table 4 shows the UQ performance of GP w/ deep
kernel learning (DKL) and dropout. The prediction intervals are in streamflow
units. We make the following observations. 1) Dropout produces very small
prediction intervals but the coverage is low (less than 80 percent for all but
one watershed) which is not desirable. 2) GP w/ DKL achieves significantly
higher coverage (more than 91 percent for all watersheds) at the expense of a
slight increase in the prediction interval size over the dropout method.

| Dropout | GP w/ DKL
---|---|---
coverage | PI width | Coverage | PI width
BR | 74.07 | 4.77 | 91.66 | 5.99
C | 76.85 | 4.7 | 93.52 | 6.3
CCR | 79.62 | 5.47 | 92.59 | 6.27
F | 75.92 | 4.92 | 94.44 | 6.28
SFC | 78.7 | 4.59 | 95.73 | 6.08
Y | 81.48 | 4.37 | 94.44 | 5.84
Table 4: Comparison of marginal coverage and prediction interval size for
dropout and GP w/ DKL approaches.

Modeling streamflow in human-influenced watersheds. The VIC-CropSyst model
simulates “natural” flows without human interventions such as reservoir
operations or irrigation. Hence, its performance tends to be poor in human-
influenced watersheds such as the Yakima River basin, where much river water
is diverted for agriculture. The RCGNN-CRL-IW model shows better performance
in such contexts, capturing some human influence nuances. Nonetheless, it
exhibits relatively lower NNSE values, highlighting the challenge of modeling
streamflow in human-influenced watersheds. This underscores the necessity to
explicitly incorporate elements of human influence into the learning approach.

##  5 Roadmap to Deployment

The intended application base for our streamflow prediction methodology are
federal and state/regional agencies that manage reservoirs or engage in other
forms of water resources management and planning. However, to enable such
partnerships in the near term, several improvements and extensions have been
planned. First, the problem formulation needs to be extended to a forecast
setting where expectations of meteorological and other input variables inform
the streamflow predictions. Datasets for expectations of weather are available
in the same format as our current inputs and can be directly integrated. Our
current prediction accuracies provide an upper bound which will diminish when
errors associated with expectations of input variables are added in. To this
end, the provision of uncertainty bounds supported by our framework will be
crucial for the end-user. User interfaces that display the uncertainty metrics
in a way that is most useful for decision-making are needed for deployment. We
have been exploring the current decision-making interfaces used by the agency
to seamlessly integrate. Finally, our stakeholder agency identified
explainability of the model as a key feature of interest in transitioning to
deployment and we can adapt existing work on explainable AI for this purpose.

##  6 Summary and Future Work

This paper studied a constrained reasoning and learning approach for
streamflow prediction by combining domain knowledge with deep temporal models,
and Gaussian process with deep kernel learning approach for uncertainty
quantification. Our experimental results demonstrated the effectiveness of our
overall methodology on diverse real-world watershed datasets and in
guaranteeing that constraints are satisfied for every test-time prediction.
One unexplored but important challenge is to model and reason about the human
influence on water availability to further improve accuracy.

Acknowledgements. This research was supported in part by United States
Department of Agriculture (USDA) NIFA award No. 2021-67021-35344 (AgAID AI
Institute), and by the USDA NIFA Award #1016467 under the Water for
Agriculture program.

## References

  * Ahmed et al. [2022a] Kareem Ahmed, Tao Li, Thy Ton, Quan Guo, Kai-Wei Chang, Parisa Kordjamshidi, Vivek Srikumar, Guy Van den Broeck, and Sameer Singh.  Pylon: A pytorch framework for learning with constraints.  In NeurIPS 2021 Competitions and Demonstrations Track, pages 319–324. PMLR, 2022.
  * Ahmed et al. [2022b] Kareem Ahmed, Stefano Teso, Kai-Wei Chang, Guy Van den Broeck, and Antonio Vergari.  Semantic probabilistic layers for neuro-symbolic learning.  In Advances in Neural Information Processing Systems 35 (NeurIPS), dec 2022.
  * Cai et al. [2021] Shengze Cai, Zhiping Mao, Zhicheng Wang, Minglang Yin, and George Em Karniadakis.  Physics-informed neural networks (PINNs) for fluid mechanics: A review, 2021.
  * Choi et al. [2020] Y Choi, Antonio Vergari, and Guy Van den Broeck.  Probabilistic circuits: A unifying framework for tractable probabilistic models.  UCLA. URL: http://starai. cs. ucla. edu/papers/ProbCirc20. pdf, 2020.
  * Cortes et al. [2010] Corinna Cortes, Yishay Mansour, and Mehryar Mohri.  Learning bounds for importance weighting.  In Advances in neural information processing systems, volume 23, 2010.
  * Hinton et al. [2012] Geoffrey E Hinton, Nitish Srivastava, Alex Krizhevsky, Ilya Sutskever, and Ruslan R Salakhutdinov.  Improving neural networks by preventing co-adaptation of feature detectors.  arXiv preprint arXiv:1207.0580, 2012.
  * Hoedt et al. [2021] Pieter-Jan Hoedt, Frederik Kratzert, Daniel Klotz, Christina Halmich, Markus Holzleitner, Grey S Nearing, Sepp Hochreiter, and Günter Klambauer.  Mc-lstm: Mass-conserving lstm.  In International conference on machine learning, pages 4275–4286. PMLR, 2021.
  * Jia et al. [2019] Xiaowei Jia, Jared Willard, Anuj Karpatne, Jordan Read, Jacob Zwart, Michael Steinbach, and Vipin Kumar.  Physics guided rnns for modeling dynamical systems: A case study in simulating lake temperature profiles.  In Proceedings of the 2019 SIAM international conference on data mining, pages 558–566. SIAM, 2019.
  * Jia et al. [2021] Xiaowei Jia, Jared Willard, Anuj Karpatne, Jordan S Read, Jacob A Zwart, Michael Steinbach, and Vipin Kumar.  Physics-guided machine learning for scientific discovery: An application in simulating lake temperature profiles.  ACM/IMS Transactions on Data Science, 2(3):1–26, 2021.
  * Kamruzzaman et al. [2021] Sarker Md Kamruzzaman, Zhou Lu, Eberhart Aaron, and Hitzler Pascal.  Neuro-symbolic artificial intelligence: Current trends.  arXiv preprint arXiv: 2105.05330, 2021.
  * Karpatne et al. [2016] Anuj Karpatne, Ankush Khandelwal, Xi Chen, Varun Mithal, James Faghmous, and Vipin Kumar.  Global monitoring of inland water dynamics: State-of-the-art, challenges, and opportunities.  Computational sustainability, pages 121–147, 2016.
  * Karpatne et al. [2017a] Anuj Karpatne, Gowtham Atluri, James H Faghmous, Michael Steinbach, Arindam Banerjee, Auroop Ganguly, Shashi Shekhar, Nagiza Samatova, and Vipin Kumar.  Theory-guided data science: A new paradigm for scientific discovery from data.  IEEE Transactions on knowledge and data engineering, 29(10):2318–2331, 2017.
  * Karpatne et al. [2017b] Anuj Karpatne, William Watkins, Jordan Read, and Vipin Kumar.  Physics-guided neural networks (pgnn): An application in lake temperature modeling.  arXiv preprint arXiv:1710.11431, 2, 2017.
  * Kratzert et al. [2019] Frederik Kratzert, Daniel Klotz, Mathew Herrnegger, Alden K Sampson, Sepp Hochreiter, and Grey S Nearing.  Toward improved predictions in ungauged basins: Exploiting the power of machine learning.  Water Resources Research, 55(12):11344–11354, 2019.
  * Liang et al. [1994] Xu Liang, Dennis P Lettenmaier, Eric F Wood, and Stephen J Burges.  A simple hydrologically based model of land surface water and energy fluxes for general circulation models.  Journal of Geophysical Research: Atmospheres, 99(D7):14415–14428, 1994.
  * Liu et al. [2023] Siyan Liu, Dan Lu, Scott L Painter, Natalie A Griffiths, and Eric M Pierce.  Uncertainty quantification of machine learning models to improve streamflow prediction under changing climate and environmental conditions.  Frontiers in Water, 5:1150126, 2023.
  * Malek et al. [2017] Keyvan Malek, Claudio Stöckle, Kiran Chinnayakanahalli, Roger Nelson, Mingliang Liu, Kirti Rajagopalan, Muhammad Barik, and Jennifer C Adam.  Vic–cropsyst-v2: A regional-scale modeling platform to simulate the nexus of climate, hydrology, cropping systems, and human decisions.  Geoscientific Model Development, 10(8):3059–3084, 2017.
  * McLain et al. [2015] Kelly McLain, Jaclyn Hancock, and Margaret Drennan.  2015 drought and agriculture, 2015.  Accessed: .
  * Muralidhar et al. [2020] Nikhil Muralidhar, Jie Bu, Ze Cao, Long He, Naren Ramakrishnan, Danesh Tafti, and Anuj Karpatne.  Phynet: Physics guided neural networks for particle drag force prediction in assembly.  In Proceedings of the 2020 SIAM International Conference on Data Mining, pages 559–567. SIAM, 2020.
  * Nash and Sutcliffe [1970] J Eamonn Nash and Jonh V Sutcliffe.  River flow forecasting through conceptual models part i—a discussion of principles.  Journal of hydrology, 10(3):282–290, 1970.
  * Nossent and Bauwens [2012] Jiri Nossent and Willy Bauwens.  Application of a normalized nash-sutcliffe efficiency to improve the accuracy of the sobol’sensitivity analysis of a hydrological model.  In EGU general assembly conference abstracts, page 237, 2012.
  * Peters-Lidard et al. [2017] Christa D Peters-Lidard, Martyn Clark, Luis Samaniego, Niko EC Verhoest, Tim Van Emmerik, Remko Uijlenhoet, Kevin Achieng, Trenton E Franz, and Ross Woods.  Scaling, similarity, and the fourth paradigm for hydrology.  Hydrology and earth system sciences, 21(7):3701–3713, 2017.
  * Rajagopalan et al. [2018] K Rajagopalan, KJ Chinnayakanahalli, CO Stockle, RL Nelson, CE Kruger, MP Brady, K Malek, ST Dinesh, ME Barber, AF Hamlet, et al.  Impacts of near-term climate change on irrigation demands and crop yields in the columbia river basin.  Water Resources Research, 54(3):2152–2182, 2018.
  * Seo et al. [2016] Youngjoo Seo, Michaël Defferrard, Pierre Vandergheynst, and Xavier Bresson.  Structured sequence modeling with graph convolutional recurrent networks, 2016.
  * Sit et al. [2021] Muhammed Sit, Bekir Demiray, and Ibrahim Demir.  Short-term hourly streamflow prediction with graph convolutional gru networks.  arXiv preprint arXiv:2107.07039, 2021.
  * Stöckle et al. [2003] Claudio O Stöckle, Marcello Donatelli, and Roger Nelson.  Cropsyst, a cropping systems simulation model.  European journal of agronomy, 18(3-4):289–307, 2003.
  * Valentin Jospin et al. [2020] Laurent Valentin Jospin, Wray Buntine, Farid Boussaid, Hamid Laga, and Mohammed Bennamoun.  Hands-on bayesian neural networks–a tutorial for deep learning users.  arXiv e-prints, pages arXiv–2007, 2020.
  * Vergari et al. [2020] Antonio Vergari, Y Choi, Robert Peharz, and Guy Van den Broeck.  Probabilistic circuits: Representations, inference, learning and applications.  AAAI Tutorial, 2020.
  * Williams and Rasmussen [2006] Christopher KI Williams and Carl Edward Rasmussen.  Gaussian processes for machine learning, volume 2.  MIT press Cambridge, MA, 2006.
  * Wilson et al. [2016] Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P Xing.  Deep kernel learning.  In Artificial intelligence and statistics, pages 370–378. PMLR, 2016.
  * Xiang and Demir [2022] Zhongrun Xiang and Ibrahim Demir.  Fully distributed rainfall-runoff modeling using spatial-temporal graph neural network.  2022\.
  * Xu et al. [2021] Da Xu, Yuting Ye, and Chuanwei Ruan.  Understanding the role of importance weighting for deep learning.  In 9th International Conference on Learning Representations, ICLR 2021, Virtual Event, Austria, May 3-7, 2021. OpenReview.net, 2021.

##  Appendix A Appendix: Recurrent Convolutional Graph Neural Network

Figure 3: Conceptual illustration of RCGNN architecture.

Figure 3 depicts the RCGNN model, an architecture that captures insights from
both spatial and temporal dimensions.

##  Appendix B Appendix: Proof for theorem 1

##  Proof for Theorem 1

######  Theorem 2.

(Theorem 1 restated: IW Improves generalization bound) Assume
$\frac{2M\log(1/\delta)}{3T}\leq\sqrt{\frac{2d_{2}(\mathcal{P}||\frac{\mathcal{P}}{\omega})\log(1/\delta)}{T}}$
in (2) of Lemma 1. Set $\omega(z,y)=\omega_{k}$ for $(z,y)\in\mathcal{D}_{k}$
such that $\mathcal{P}_{k}\cdot\omega_{k}=\frac{k^{a}}{K^{b+1}}$ for any
$k\in\\{1,...,K\\}$. Under the conditions $0<a\leq\min\\{K-1,b-1/\ln(K)\\}$,
and $K\geq 2$, the IW empirical risk $\widehat{R}_{\omega}(F)$ gives a tighter
generalization bound than empirical risk $\widehat{R}(F)$ with
$\widehat{B}_{\omega}=\sqrt{\frac{2\log(1/\delta)}{(a+1)T}}$ and
$\widehat{B}=\sqrt{\frac{2\log(1/\delta)}{T}}$.

###### Proof.

(of Theorem 1) Recall the base-$2$ exponential of Rényi divergence of order
$2$ between $\mathcal{P}$ and $\frac{\mathcal{P}}{\omega}$:

$\displaystyle
d_{2}(\mathcal{P}||\frac{\mathcal{P}}{\omega})=\int_{(z,y)}\mathcal{P}(z,y)\cdot\frac{\mathcal{P}(z,y)}{\frac{\mathcal{P}(z,y)}{\omega(z,y)}}d(z,y)=\int_{(z,y)}\mathcal{P}(z,y)\cdot\omega(z,y)d(z,y).$

It it worth noting a special case of $\omega(z,y)=1$ reduces IW empirical risk
$R_{\omega}(F)$ back to empirical risk $R(F)$:

$\displaystyle
d_{2}(\mathcal{P}||\frac{\mathcal{P}}{\omega})=d_{2}(\mathcal{P}||\frac{\mathcal{P}}{1})=\int_{(z,y)}\mathcal{P}(z,y)\cdot
1d(z,y)=1.$

Then, to prove the main result regarding generalization bounds, under the
assumption that
$\frac{2M\log(1/\delta)}{3T}\leq\sqrt{\frac{2d_{2}(\mathcal{P}||\frac{\mathcal{P}}{\omega})\log(1/\delta)}{T}}$
in the same analysis result (2) of Lemma 1, we have the following
inequalities:

$\displaystyle R(F)-\widehat{R}(F)\leq$ $\displaystyle
2\sqrt{\frac{2d_{2}(\mathcal{P}||\mathcal{P})\log(1/\delta)}{T}}=2\sqrt{\frac{2\log(1/\delta)}{T}}=\widehat{B},$
$\displaystyle R(F)-\widehat{R}_{\omega}(F)\leq$ $\displaystyle
2\sqrt{\frac{2d_{2}(\mathcal{P}||\frac{\mathcal{P}}{\omega})\log(1/\delta)}{T}}=\widehat{B}_{\omega}.$
(6)

To show that IW empirical risk $\widehat{R}_{\omega}(F)$ gives a tighter
generalization bound than empirical risk, it suffixes to prove that
$d_{2}(\mathcal{P}||\frac{\mathcal{P}}{\omega})\leq
d_{2}(\mathcal{P}||\mathcal{P})=1$. Now we divide
$d_{2}(\mathcal{P}||\frac{\mathcal{P}}{\omega})$ that involves all possible
data into $K$ groups as follows

$\displaystyle d_{2}(\mathcal{P}||\frac{\mathcal{P}}{\omega})=$
$\displaystyle\int_{(z,y)}\mathcal{P}(z,y)\cdot\Bigg{(}\sum_{k=1}^{K}\mathbb{I}[(z,y)\in\mathcal{D}_{k}]\Bigg{)}\cdot\omega(z,y)d(z,y)$
$\displaystyle=$
$\displaystyle\sum_{k=1}^{K}\int_{(z,y)}\mathcal{P}(z,y)\cdot\mathbb{I}[(z,y)\in\mathcal{D}_{k}]\cdot\omega(z,y)d(z,y)$
$\displaystyle\stackrel{{\scriptstyle(a)}}{{=}}$
$\displaystyle\sum_{k=1}^{K}\int_{(z,y)}\mathcal{P}(z,y)\cdot\mathbb{I}[(z,y)\in\mathcal{D}_{k}]\cdot\omega_{k}d(z,y)=\sum_{k=1}^{K}\mathcal{P}_{k}\cdot\omega_{k},$

where the above equality $(a)$ is due to the setting of group-wise important
weights $\omega(z,y)=\omega_{k}$ for $(z,y)\in\mathcal{D}_{k}$. Then under the
setting of group-wise weights $\omega_{k}$ such that
$\mathcal{P}_{k}\cdot\omega_{k}=\frac{k^{a}}{K^{b}}\cdot\frac{1}{K}$ with
$a>0$ a for all $k\in[K]$, we can further upper bound $d_{2}$ as follows:

$\displaystyle d_{2}(\mathcal{P}||\frac{\mathcal{P}}{\omega})=$
$\displaystyle\sum_{k=1}^{K}\mathcal{P}_{k}\cdot\omega_{k}=\sum_{k=1}^{K}\frac{k^{a}}{K^{b}}\cdot\frac{1}{K}\leq\frac{\int_{k=1}^{K+1}k^{a}dk}{K^{b+1}}\leq\frac{(K+1)^{(a+1)}}{(a+1)K^{b+1}}=\frac{1}{a+1}\cdot\bigg{(}\frac{K+1}{K}\bigg{)}^{a+1}\cdot
K^{a-b}$ $\displaystyle=$
$\displaystyle\frac{1}{a+1}\cdot\bigg{(}1+\frac{1}{K}\bigg{)}^{a+1}\cdot
K^{a-b}\leq\frac{1}{a+1}\cdot\exp\bigg{(}\frac{a+1}{K}\bigg{)}\cdot
K^{a-b}\leq\frac{1}{a+1}\cdot\exp(1)\cdot K^{a-b}\leq\frac{1}{a+1},$ (7)

where the first inequality is due to the left rule of Riemann sum for the
monotonically increasing function $k^{a}$ with $a>0$, the third inequality is
due to $1+x\leq\exp(x)$, the fourth inequality is due to the condition $a\leq
K-1$, and the last inequality is due to the condition $a\leq b-1/\ln(K)$.

As can be seen, in (Proof.), $\widehat{B}_{\omega}<\widehat{B}$ holds:

$\displaystyle\widehat{B}_{\omega}=2\sqrt{\frac{2d_{2}(\mathcal{P}||\frac{\mathcal{P}}{\omega})\log(1/\delta)}{T}}\leq
2\sqrt{\frac{2\log(1/\delta)}{(a+1)T}}<2\sqrt{\frac{2\log(1/\delta)}{T}}=\widehat{B},$

where the inequality is due to (Proof.) and the condition $0<a$. This
immediately shows that IW empirical risk $\widehat{R}_{\omega}(F)$ gives a
tighter generalization bound than empirical risk $\widehat{R}(F)$. ∎

##  Appendix C Appendix: Additional Experimental Results

Figure 4: Additional analysis for both low- and high-flow months: Relative
performance chart comparing all models for (a) Boise, (b) Clearwater, (c)
Clearwater Canyon Ranger, (d) Flathead, (e) Southfork Clearwater, and (f)
Yakima. For each test month and year, each model’s MAE was compared against
the best performing model’s MAE for that month and year (difference shown in
X-axis). The fraction of test month-years for which that performance is
achieved is shown along the Y-axis. Closer a curve is aligned with the Y-axis
for larger fraction, the better.

Figure 4 shows the relative performance of the different models based on their
performance over all the test months and years. From an overall perspective,
the VIC-CropSyst model yields the best performance in the majority of the
month-year cases, but as highlighted in the main paper Fig 2., it fails to
perform well in high-flow months. This is also demonstrated in Figure 5 and
Figure 6. Accurately capturing flows during these high-flow months is crucial
for our application context because the estimated inflows directly affect
reservoir operations and are the primary focus of this paper. Figure 5 shows
that the VIC-CropSyst and RCGNN model fails to capture the peak of the
streamflow that leads to higher error in the flow during these months (Figure
6). In addition, in some watersheds, we observe the shift in the timing of the
peak streamflow (e.g., Flathead) that degrades the performance of the VIC-
CropSyst and RCGNN models. While other data-driven models with included
knowledge about the constraints perform better.

Figure 5: Time-series plot of simulated and observed streamflow for the water year 2011.  Figure 6: Median of monthly error (prediction - actual) for different watersheds. A positive value of error indicates the overestimation of simulated flow by model and a negative value indicates under-estimation |  VIC- CropSyst |  RCGNN |  PG- RCGNN |  RCGNN- CRL |  RCGNN- CRL-IW
---|---|---|---|---|---
BR |  0.839 |  0.842 |  0.8443 |  0.8578 |  0.8651
C |  0.8551 |  0.8536 |  0.8676 |  0.9018 |  0.9101
CCR |  0.8324 |  0.8412 |  0.851 |  0.9058 |  0.9107
F |  0.7133 |  0.7937 |  0.8218 |  0.855 |  0.8867
SFC |  0.7692 |  0.7659 |  0.81 |  0.8578 |  0.8723
Y |  0.2315 |  0.2873 |  0.3798 |  0.5073 |  0.5277
Table 5: Comparison of predictive models across watersheds during all months.
Values are reported in terms of the NNSE metric (higher the better). NNSE
values less than 0.5 indicate that the mean is a better predictor than the
model’s output.

Table 5 underscores the superior performance of RCGNN-CRL-IW over all other
baselines, even when evaluating across all months. These results align with
the findings drawn from comparing model performance specifically during high
flow months, as illustrated in Table 2. This table provides supplementary
details regarding the evaluation of model performance discussed in the Table
2. While the primary focus of the paper is on high flow months, overall
results are also included here for comprehensive understanding. It is
noteworthy that the overall performance, as presented in this appendix,
exhibits slight variations as we go from high flow months to all months.

| PG-RCGNN | RCGNN-CRL | RCGNN-CRL-IW
---|---|---|---
BR | 0.8192 | 0.8242 | 0.8432
C | 0.8471 | 0.8588 | 0.8713
CCR | 0.8288 | 0.8461 | 0.8571
F | 0.7563 | 0.7703 | 0.7835
SF | 0.7622 | 0.7811 | 0.7986
Y | 0.3675 | 0.3851 | 0.4368
Table 6: Comparison of predictive models across watersheds using bi-weekly
data. Values are reported in terms of the NNSE metric (higher the better).
NNSE values less than 0.5 indicate that the mean is a better predictor than
the model’s output.

Table 6 presents the NNSE comparison for the models utilizing bi-weekly data,
highlighting the effectiveness of RCGNN-CRL-IW even within a lower granularity
framework. Additionally, the table reveals a marginal decrease in overall
performance across all models as we go from monthly to bi-weekly datasets.

| PG-RCGNN | RCGNN-CRL | RCGNN-CRL-IW
---|---|---|---
BR | 0.8136 | 0.8204 | 0.8317
C | 0.8337 | 0.8396 | 0.8526
CCR | 0.8053 | 0.8178 | 0.8382
F | 0.7349 | 0.7461 | 0.7512
SF | 0.7472 | 0.7522 | 0.7639
Y | 0.3418 | 0.3698 | 0.3943
Table 7: Comparison of predictive models across watersheds using weekly data.
Values are reported in terms of the NNSE metric (higher the better). NNSE
values less than 0.5 indicate that the mean is a better predictor than the
model’s output.

Table 7 reaffirms the dominance of RCGNN-CRL-IW when utilizing weekly
observations, outperforming other models. Moreover, we observe a consistent
trend of decreasing performance as granularity decreases.

| high flow months | all months
---|---|---
| PG-
---
RCGNN
| PG-
---
RCGNN-IW
| PG-
---
RCGNN
| PG-
---
RCGNN-IW
BR | 0.919 | 0.924 | 0.844 | 0.853
C | 0.925 | 0.931 | 0.868 | 0.879
CCR | 0.896 | 0.905 | 0.851 | 0.875
F | 0.902 | 0.911 | 0.822 | 0.846
SFC | 0.866 | 0.883 | 0.81 | 0.832
Y | 0.509 | 0.542 | 0.38 | 0.407
Table 8: Comparison of PG-RCGNN and PG-RCGNN-IW (PG-RCGNN with importance
weighted training) across watersheds during high flow and overall months.
Values are reported in terms of the NNSE metric (higher the better). NNSE
values less than 0.5 indicate that the mean is a better predictor than the
model’s output.

Table 8 shows that the importance weighted training contributes to the
performance of PG-RCGNN in both high flow months and all months. This further
showcases the generalizability and effectiveness of IW training in improving
performance.

##  Appendix D Hyper-parameter Tuning

We employed the validation data to select the hyper-parameters for all deep
neural network approaches (RCGNN, PG-RCGNN, RCGNN-CRL, RCGNN-CRL-IW). We
selected a learning rate of 0.001, a weight decay of 0.0005, a hidden
dimension of 256 for the predictive layers, and a dropout rate of 0.2.