Using Data-Driven Prediction of Downstream 1D River Flow to Overcome the Challenges of Hydrologic River Modeling

Abstract

Methods for downstream river flow prediction can be categorized into physics-based and empirical approaches. Although based on well-studied physical relationships, physics-based models rely on numerous hydrologic variables characteristic of the specific river system that can be costly to acquire. Moreover, simulation is often computationally intensive. Conversely, empirical models require less information about the system being modeled and can capture a system’s interactions based on a smaller set of observed data. This article introduces two empirical methods to predict downstream hydraulic variables based on observed stream data: a linear programming (LP) model, and a convolutional neural network (CNN). We apply both empirical models within the Colorado River system to a site located on the Green River, downstream of the Yampa River confluence and Flaming Gorge Dam, and compare it to the physics-based model Streamflow Synthesis and Reservoir Regulation (SSARR) currently used by federal agencies. Results show that both proposed models significantly outperform the SSARR model. Moreover, the CNN model outperforms the LP model for hourly predictions whereas both perform similarly for daily predictions. Although less accurate than the CNN model at finer temporal resolution, the LP model is ideal for linear water scheduling tools.

1. Introduction

Rivers play a keystone role in the functioning of human societies, where they provide wildlife habitat, freshwater, food supply, transportation, and energy [1]. The security of critical water resources like rivers is increasingly threatened by climate change, population growth, and pollution, among other factors [2]. As a result, water resources managers must consider a complex and diverse slate of stakeholders including farmers and other irrigators, power suppliers, fishermen, floodplain dwellers, ecologists, and hobbyists. Effective water management strategies that understand the effects of human influences on river flows can aim to minimize threats like flood, drought, species endangerment, or hydropower failure. The development of models for river flow prediction, for example river stream discharge and stage prediction, is therefore a crucial component for maintaining water security, safeguarding environmental quality, and promoting sustainable development [3]. Hydrologic models can simulate river flow in 1D (e.g., along the river path), 2D (e.g., in planar representation), or 3D (e.g., accounting for vertical and lateral dynamics). 2D and 3D models are generally more accurate, but computational costs increase with the number of dimensions being modeled [4]. This study addresses the challenge of 1D river flow modeling using data-driven methods where observation data within a watershed enables the prediction of downstream hydraulic variables.

Characterizing the driving mechanisms of water storage and movement of water within watersheds, a process known as the hydrologic cycle, is a core focus of hydrology. Accurate prediction of downstream river flow is complicated by a variety of influencing factors including precipitation and its heterogenous spatial distribution, other weather factors, fluxes via ground water exchange and irrigation, and ungauged stream inflows. Over the years, numerous methods have been developed for water flow prediction to address these challenges with differing success. In this brief literature review we categorize some of the most common methods for river flow prediction into empirical approaches and physics-based approaches.

Empirical approaches leverage statistical modeling, machine learning or other data-driven methods to predict hydraulic variables. These methods rely on their ability to capture a system’s interactions based on observed data rather than through formulating relationships between underlying physical principles. As such, they require less information about the system being modeled and may be more flexible in adapting to new data or changes in the system. Statistical techniques based on stochastic process modeling are the most traditional methods for empirical water flow prediction [5,6,7,8]. Generally, these include the autoregressive model (AR) [9,10], moving average model (MA) [11,12], autoregressive and moving average model (ARMA) [13], autoregressive integrated moving average model (ARIMA) [14], and seasonal decomposition of time series. The optimization of statistical parameters is performed using historical data, producing time series models capable of reasonable short-term forecasting of stationary hydraulic variables using linear combinations of past values.

Newer data-driven methods inhabit a rapidly growing body of research including linear programming, machine learning, and artificial intelligence. Machine learning models such as Support Vector Machines (SVMs) [15,16] and Artificial Neural Networks (ANNs) [14,17] can learn complex patterns in data and capture non-linear relationships between variables. Empirical approaches also can impose inductive biases that guide the learning of parameters within a set of structural assumptions leading to diverse subtypes of ANNs, for example time series models that use convolutions and recurrent connections [18,19]. For this reason, such models are increasingly used for processes within the hydrologic cycle and specifically for river flow prediction.

Physics-based models attempt to simulate the hydrologic cycle based on known characteristics of the watershed. Some considerations include soil moisture, infiltration, land use, and surface friction [20], each of which plays key roles in water partitioning, routing, and flow speed. Depending on the objective, these factors may be represented using known physical relationships. For example, the Gauckler–Manning formula estimates velocity within an open channel based on cross-sectional geometry, topography, and friction [21]. Physics-based models widely used by the community include Soil and Water Assessment Tool (SWAT) [22], the Hydrologic Engineering Center’s Hydrologic Modeling System (HEC-HMS) [23], and the Variable Infiltration Capacity (VIC) model [24], each with varying applications and features.

One disadvantage of physics-based models is that they can be more complex and computationally expensive to develop and use. These models require a thorough understanding of the underlying physical principles and the relevant governing equations, as well as a detailed knowledge of the system’s geometry, boundary conditions, and material properties. Another disadvantage of physics-based models is that they may not always accurately capture all the complex interactions and feedback mechanisms that can occur within a system. In contrast, empirical models are often simpler and easier to develop because they are based on observed data or relationships. They generally require less information about the system being modeled and may be more flexible in adapting to new data or changes in the system. Empirical models are also able to capture a system’s complex interactions and feedback mechanisms more accurately than a physics-based model because they are based on observed data or relationships, rather than on assumptions about the underlying physical principles.

However, the main disadvantage of empirical models is the need for large training datasets. Because training datasets are key to empirical models, they are not suitable for cases where historical data are limited or absent. For example, physics-based models are more suitable for river systems where measuring gages have only recently been active, where the watershed structure was radically and suddenly modified (e.g., due to extreme events such as flooding or earthquake), or for hypothetical/future river systems that do not exist yet. Another disadvantage of empirical models is their lack of interpretability. Like most machine learning models, the models described in this manuscript cannot provide a deep insight into the links between the learned features and the physical conditions of the system, whereas physics-based models aim to explicitly describe the physical mechanisms at play.

In this paper, two empirical methods inspired by traditional statistical techniques are introduced and compared for their ability to predict river flow downstream of the Flaming Gorge Dam. The first model employs Linear Programming (LP) optimization to learn discrete unit hydrographs at multiple flow magnitudes. The next model is based on a 1D convolutional neural network (CNN) and encodes upstream hydrograph signals to use in downstream prediction. Both models use the same National Water Information System (NWIS) data sources within the Middle Green River basin to provide defensible comparisons between the two approaches. Alongside this comparison, simulations from the physics-based Streamflow Synthesis and Reservoir Regulation (SSARR) model [25] are used to supplement the comparison.

2. Materials and Methods

2.1. The SSARR Model

The SSARR model is a physics-based model that is designed for operational use in hydrologic engineering studies and daily streamflow forecasting [25]. The model has been regularly updated in the past decades and is still in use today in several water management models. For example, the SSARR model is used by the Bureau of Reclamation (BoR) in their RiverWare implementation of the Colorado River Simulation System (CRSS) [26]. It is also used by the Western Area Power Administration (WAPA) in the hydropower scheduling model GTMaxSL [27] to predict the water flow and level at the Jensen gage (USGS 09261000) [28], downstream the Flaming Gorge dam. However, because SSARR is a physics-based model, river flow and stage prediction are calculated based on a large set of physical equations (e.g., generalized snowmelt, channel routing) and watershed data (e.g., drainage area, temperature index, evapotranspiration index, soil moisture index) [25] and several equation parameters are “determined by trial-and-error” [25].

The model introduced in this study and used by WAPA for hydropower scheduling uses SSARR’s basic 1D routing method [25] to route water through a river system. The law of continuity in storage provides the numerical scheme for flow through a channel:

$$I=O+\frac{dS}{dt}$$

where $I$, $O$, and $S$ are the inflow, outflow, and storage of the channel over a given duration $t$.

The magnitude of flow through a channel dictates the amount of time that flow takes to traverse the channel, referred to as the duration of storage. Duration of storage ($T_s$) is given by:

$$T_s=\frac{K_{T_s}}{Q^n}$$

where $Q$ is discharge, $K_{T_s}$ is a constant determined from physical measurements of flow and routing times, and $n$ is a coefficient. The river reach between the Flaming Gorge Dam and the Jensen, UT station is represented as five reaches in SSARR, each with multiple sub reaches, i.e., channel components (25, 35, 40, 25, and 23 sub reaches for each respective reach). Routing for every sub reach within a reach is calculated using the same $K_{T_s}$ and $n$ values.

Contrary to SSARR, the linear programming (LP) and convolutional neural network (CNN) methods presented in the following sections are both empirical models. These new models are designed to learn how to predict the downstream water flow and/or water level of a river system solely based on observational data. Compared to the SSARR model or other physics-based models, these models do not need to characterize the dynamics of the physical systems and their predictive parameters can be automatically updated and learned from the most recent historical water flow and water level data.

2.2. The Discrete Convolution Approach

The two empirical models introduced here are both based on variants of a discrete convolution function, where it is assumed that the downstream river flow can be modeled as the convolution of upstream river flows with their corresponding unit hydrographs [29]. It should be noted that the approach used here is different from the autoregressive (AR) approach traditionally used in water flow prediction (which can also be interpreted as a discrete convolution). Specifically, the AR approach predicts flow at a river point based on flow at the same point in preceding hours. The convolution approach used here instead considers upstream flow when predicting the downstream flow rate.

Consider sequences $f$ and $g$, where $f$ is an endogenous time series and a filter is given by sequence $g=(g_1,\cdots ,g_M)$. The $k$-th element of the discrete convolution of $f$ and $g$ is written as:

$${\left(f\ast g\right)}_k=\sum _{m=1}^M{g_{m}f}_{k-m},$$

(1)

This discrete convolution formula is used in more advanced empirical implementations presented below.

2.3. The Linear Programming Model

Linear programming (LP), or linear optimization, is a mathematical method used to identify the minimum or maximum value of an objective function with respect to requirements that are modeled by a set of linear equality and inequality constraints [30]. In the LP model introduced here, we aim to minimize the error in predicting the downstream river discharge. The LP model assumes that the water discharge profile at the downstream point of a river network is equal to the sum of the individual discharge impacts from its upstream sources. Each individual discharge impact, in turn, can be estimated in the medium-term range (i.e., a few days to a month) as the convolution between a unit hydrograph and the water discharge at the upstream source. Based on these assumptions, the LP model identifies the optimal set of unit hydrographs [29] that minimizes the error in predicting the discharge level at a downstream river point based on the discharge level at the upstream sources. The LP model also simultaneously maximizes the smoothness of the identified unit hydrographs based on a predefined smoothness coefficient. An important contribution of this LP model is that it allows the unit hydrograph of each reach to slowly evolve over time according to a predefined set of evolving linear coefficients. These evolving linear coefficients can be guided by long-term trends such as seasons, average temperatures, or seasonal inflow levels.

The LP model described here is inspired by the one introduced in [31]. The main contributions of the new LP model proposed here are the following:

• Instead of predicting the downstream discharge by identifying the optimal unit hydrograph of a single upstream source, our model simultaneously identifies the optimal unit hydrograph of multiple contributing upstream sources;
• Instead of assuming the unit hydrograph of each upstream source to be static, the model allows the identification of dynamic unit hydrographs;
• Apart from minimizing the error in predicting the downstream flow, the model also maximizes the smoothness of the identified unit hydrographs.

The LP model can be described by the compact mathematical formulation (2)–(8) below. A more complete formulation, together with a description of each set and variable, is provided in Appendix A.

$$Min\left({\sum }_{t=0}^{T-1}{\delta }_t+c{\sum }_{k=1}^K{s}_{k,l}\right),$$

(2)

$$\mathrm{s}.\mathrm{t}. \left|Q_t^d-\widehat{Q_t^d}\right|\le {\delta }_t$$

(3)

$$\widehat{Q_t^d}={\sum }_{k=1}^K\widehat{Q_{k,t}^d},$$

(4)

$$\widehat{Q_{k,t}^d}={\sum }_{l=1}^L{a}_{k,l,t}{q}_{k,l,t}^d$$

(5)

$$q_{k,l,t}^d=h_{k,l,t}\ast Q_{k,t}^u={\sum }_{u=0}^{T^H-1}{h}_{k,l,u}{Q}_{k,t-u}^u$$

(6)

$$d_{k,l,t}=h_{k,l,t+1}-2h_{k,l,t}+h_{k,l,t-1}$$

(7)

$$\left|d_{k,l,t}\right|\le s_{k,l}$$

(8)

Equation (2) is the objective function and simultaneously minimizes the estimation error of the predicted downstream flow and the variability of the unit hydrographs. More specifically, Equation (2) concurrently minimizes (a) the sum of absolute differences between the downstream flow estimator and its historical value, and (b) the sum of maximum absolute value of the second order derivative of each unit hydrograph. Equation (3) define ${\delta }_t$ as an upper bound of the error in approximating the downstream flow in time t. Equation (4) defines the downstream flow estimator as the sum of the individual upstream components. Equation (5) allows each upstream component $\widehat{Q_{k,t}^d}$ of the downstream flow estimator to be defined as an evolving linear combination of subcomponents $q_{k,l,t}^d$. Each subcomponent $q_{k,l,t}^d$ is defined as the convolution between the historical upstream flow $Q_{k,t}^u$ and an elementary unit hydrograph $h_{k,l,t}$, as described in Equation (6). The coefficients $a_{k,l,t}$ of the evolving linear combination are predefined by the user and can be based on the time of the year, weather conditions, or hydrology conditions. Equation (7) defines $d_{k,l,t}$ as the discrete second order derivative of the unit hydrographs. Equation (8), together with the minimization imposed by Equation (2), ensure that the maximum absolute value of $\left|d_{k,l,t}\right|$ is as low as possible. Note that minimizing the maximum absolute value of the second order derivative of each unit hydrograph is equivalent to maximizing their smoothness.

As seen in Equations (5) and (6), the presented formulation does not assume the unit hydrograph of a given reach $k$ to be static. Instead, the formulation assumes that the water flow component $\widehat{Q_{k,t}^d}$ of reach $k$ is an evolving linear combination of $L$ subcomponents $q_{k,l,t}^d$ that are each associated to a static unit hydrograph $h_{k,l,u}$. Mathematically, we can write:

$$\begin{array}{l}\widehat{Q_{k,t}^d}=\sum _{l=1}^L{a}_{k,l,t}{q}_{k,l,t}^d=\sum _{l=1}^L{a}_{k,l,t}\left(\sum _{u=0}^{T^H-1}{h}_{k,l,u}{Q}_{k,t-u}^u\right)\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} =\sum _{u=0}^{T^H-1}\left(\sum _{l=1}^L{a}_{k,l,t}{h}_{k,l,u}\right)Q_{k,t-u}^u=\sum _{u=0}^{T^H-1}{H}_{k,t,u}{Q}_{k,t-u}^u\end{array}$$

(9)

with $H_{k,t,u}={\sum }_{l=1}^L{a}_{k,l,t}{h}_{k,l,u}$.

In other words, the estimator component $\widehat{Q_{k,t}^d}$ can be interpreted as the convolution between the historical upstream flow $Q_{k,t}^u$ and an evolving unit hydrograph $H_{k,t,u}$, and this evolving unit hydrograph is defined as an evolving linear combination of elementary unit hydrographs $h_{k,l,u}$.

However, the convolution ${\sum }_{u=0}^{T^H-1}{H}_{k,t,u}{Q}_{k,t-u}^u$ only makes sense if the value of $H_{k,t,u}$ is almost constant when $t$ varies over a period of length $T^H$. In other words, the convolution makes sense if we assume that the parameter $a_{k,l,t}$ varies slowly over a time period of length $T^H$, i.e., $\underset{t\le u\le t+T^H-1}{\max }\left|a_{k,l,t}-a_{k,l,u}\right|\ll 1$.

Therefore, the input coefficients $a_{k,l,t}$ of this model must describe a slow long-term trend that has negligible variations over the reach’s average water travel time.

Note that the proposed LP model is designed to predict a river downstream discharge based on upstream discharges. However, it might be more relevant at times to predict the downstream river stage instead of the downstream river discharge. Contrary to the downstream river discharge, the relationship between the upstream river flow and the downstream river stage is highly non-linear, which prevents the river stage from being directly modeled into an LP formulation. However, the river stage of a specific point of a river is known to be a function of the river discharge [32]. Moreover, the relationship between the river stage and the river discharge can be modeled using historical data.

Consequently, the LP model described above can be used to identify/learn the optimal unit hydrographs that minimize the error in predicting the downstream river discharge based on upstream flows. These unit hydrographs can later be applied to new upstream flow data to predict future downstream river discharge. Finally, an empirical discharge-to-stage function can be applied to the predicted downstream discharge to predict the downstream stage.

The LP model is solved using the PuLP Python package with Gurobi (9.5.1) solver [33].

2.4. The Convolutional Neural Network Encoder

The CNN is a type of artificial neural network (ANN) that uses convolutional filters to learn space-dependent patterns. The building block of all ANNs is the multilayer perceptron (MLP), a series of stacked units formulated as [34]:

$$\begin{array}{l}{\mathbf{x}}_{}^{l+1}=\sigma \left({{\mathbf{w}}_{}^{l+1}}_{}^T{\mathbf{x}}_{}^l+{\mathbf{b}}_{}^{l+1}\right)\\ {\mathbf{x}}_{}^L={{\mathbf{w}}_{}^L}_{}^T{\mathbf{x}}_{}^{L-1}+{\mathbf{b}}_{}^L\end{array}$$

(10)

where ${\mathbf{w}}_{}^L$ and ${\mathbf{b}}_{}^L$ are learnable parameter vectors for layer $l$, ${\mathbf{x}}_{}^l$ is the output vector of layer $l$ (${\mathbf{x}}_{}^0$ is the input vector, ${\mathbf{x}}_{}^L$ is the MLP output), and $\sigma$ represents a non-linear activation, e.g., the sigmoid ($\sigma \left(x\right)=max(0,x)$) or tanh ($\sigma =tanh$) function. Depending on the prediction domain, deep learning often achieves improved performance by including restrictive assumptions (known as inductive bias) within model architectures. In image and time-series problems where data adheres to common spatial or temporal patterns, this is achieved through convolutions with learnable filters. The neural network used here, known as a fully convolutional encoder network, is made entirely from convolutional layers and formulated as [34]:

$$\begin{array}{l}{\mathbf{x}}_k^{l+1}=\sigma \left(b_k^{l+1}+\sum _{i=0}^{C^{l-1}}{{\mathbf{w}}_{ik}^{l+1}}_{}^T\ast {\mathbf{x}}_i^l\right)\\ {\mathbf{x}}_{}^L=b_{}^L+\sum _{i=0}^{C^{L-1}}{{\mathbf{w}}_i^L}_{}^T\ast {\mathbf{x}}_i^{L-1}\end{array}$$

(11)

where ${\mathbf{x}}_k^l$ is the output of neuron $k$ in layer $l$, $b_k^l$ is the bias term for neuron $k$ in layer $l$, and ${\mathbf{w}}_{ik}^l$ is the filter (kernel) of neuron $k$ in layer $l$ applied to the activations of $i$-th neuron (channel) of layer $l-1$. The number of layers $L$, number of filters (and output channels) per layer $C^l$, and the lengths of weight kernels ${\mathbf{w}}_{ik}^l$ are hyperparameters that may require optimization. Note that the equation presented here for 1D convolutions uses the neuron-wise formulation (in contrast to Equation (10) for MLP where vector-wise formulation is used) to emphasize flexibility in the number and size of kernels (${\mathbf{w}}_{ik}^l$).

For training and prediction of Jensen flow, 15 min discharge series from the Greendale and Deerlodge USGS stream gages are produced by applying a discharge-variant shift described in the following section. The two shifted series are then used as a 2-channel input into the CNN. The prediction of Jensen stage is given by the center-cropped output of the last convolution layer. Input series are 3-day windows and the model output is a 1-day window.

The TensorFlow Python package is used for deep learning [35]. For tuning, Bayesian optimization via the GPyOpt Python package is performed on the hyperparameter input space presented in

Table A1. Search space for Bayesian optimization and optimal values of CNN hyperparameters. Number of convolutional layers refers to the number of stacked convolution operations in the CNN, represented as L in Equation (11); number of convolution filters refers to the number of independent convolution designs in each convolution layer, represented as C in Equation (11); kernel size is the number of contiguous timesteps used in a convolution, represented as the length of vector w in Equation (11); kernel non-negative constraint refers to a numerical constraint added to the convolution filters; kernel regularizer L2 factor is L2 regularization factor applied to the kernel weights during training; learning rate is used to adjust the relative effect of model loss on parameter updates after each training epoch; and batch size is the number of data points used in one training step. Type indicates whether the variable has discrete types (e.g., “discrete” and “Boolean”) or can occupy a range of values (e.g., “continuous”). Domain defines the numerical space to search.

Parameter	Type	Domain	Optimized Value	Note
Number of convolution layers	Discrete	(1, 10)	5
Number of convolution filters	Discrete	{4, 8, 16, 32, 64}	16
Kernel size	Discrete	(4, 193)	24	Subject to number of convolutions
Kernel non-negative constraint	Boolean	{True, False}	False
Kernel regularizer L2 factor	Continuous	(0, 0.1)	0.03
Learning rate	Continuous	(0.00001, 0.01)	0.002
Batch size	Discrete	{64, 128, 256}	64

[36]. The optimal CNN presented here, identified using Bayesian optimization, employs 5 convolution layers, each containing 16 filters with a kernel size of 24. The model was trained using the Adam optimizer [37] with a learning rate of 0.002, mean squared error loss, and batch size of 64. Regularization with an L2 factor of 0.03 was applied to convolution kernel weights during training.

2.5. Discharge-Variant Water Travel Time

The travel time of a wave between two points along a river reach is assumed to be exponentially proportional to the volume of water transported by that wave [38]. This informs the water travel time (WTT) estimation algorithm presented here in which travel time is relative to upstream discharge. To develop this relationship, rolling windows are applied to the first-difference upstream and downstream discharge timeseries. The average discharge of the $i$-th upstream window is assumed to drive the travel time of impulses to the $i$-th downstream window. Thus, selection of window length $L$ is prone to calibration but must include sufficient time for an impulse to appear in both up and downstream windows. Travel time for the $i$-th window, $L_i$, is estimated as the cross-correlation time lag which maximizes similarity between the window pairs, bound by previous window’s lag:

$$L_i=\underset{n\in [L_{i-1}-\delta ,L_{i-1}+\delta ]}{max}\{\sum _{m=-\infty }^{\infty }\overline{f^{\prime }\left[m\right]}{g}^{\prime }\left[m+n\right]\}$$

(12)

where $\overline{f\left[t\right]}$ denotes the complex conjugate, $f{}^{\prime }$ is the first difference (first-order backshift) of $f$, $n$ is displacement order, $f\left[t\right]$ and $g\left[t\right]$ are the $i$-th upstream and downstream hydrograph windows. The algorithm is applied across every rolling window. An exponential regression is used to calculate the power law relationship between the upstream discharge and WTT.

Upstream hydrographs are augmented by reindexing according to the lag estimated by the WTT estimator (i.e., a discharge-variant shift is applied). Any discontinuities created by this operation are filled via linear interpolation, while discontinuities that pre-existed in the original hydrographs are maintained.

Based on observation, the WTT between both upstream gages (Deerlodge, UT and Greendale, UT) rarely exceeds three days ($T_{max}$), so the window length $L={2T}_{max}+1=7d$ is chosen.

2.6. Model Validation

A period of four water years (October 2018–September 2022, Figure 1) is selected to constrain modelling by a recent timeframe in which hydraulic relationships are expected to be relatively steady. Empirical models in this study use the first three water years for training. The last year, water year 2022, is reserved for model evaluation and is known as the test set. Due to the propensity of deep learning models to overfit on training data, the training set is divided into five cross-validation folds for parameter tuning of the CNN. When evaluating the CNN’s performance, the CNN’s prediction is given by the average of five models each independently trained on a training fold. Because the CNN is trained on 15 min flow, predicted hourly flow is generated by calculating 1 h average stage.

2.7. Site Description

The site used to assess the performance of the proposed methods is the river system located on the Green River and Yampa River and delimited by the Flaming Gorge Dam (first upstream source), the Deerlodge Park gage (second upstream source), and the Jensen gage (downstream point) (see Figure 2).

The Green River is a major tributary of the Colorado River, flowing through Utah, Colorado, Wyoming and Arizona. The Flaming Gorge Dam impounds the Upper Green River, creating the 91-mile (~146-km) Flaming Gorge Reservoir, 400 miles (~643 km) upstream of the major Colorado and Green River confluence. Constructed from 1958 to 1964, the reservoir and dam play crucial roles in flood control, hydropower generation, irrigation, and recreation. The Yampa River is a major tributary of the Green River, flowing through northwestern Colorado. The Yampa River meets the Green River at the confluence near the Dinosaur National Monument in Colorado and Utah. The Jensen gage is a river gage located on the Green River near Jensen, Utah, 28.6 miles downstream of the Yampa confluence. The Jensen gage is used to measure the flow of the Green River and monitor its water levels.

Per the most recent Flaming Gorge Environmental Impact Statement (EIS) [39], during the base flow period (from mid-July to the end of February), hourly release patterns from the Flaming Gorge Dam must be not produce more than a 0.1 m stage change at the Jensen gage within a 24 h period, except during emergency operations. This restriction is part of a broader set of environmental regulations imposed on the Flaming Gorge Dam operations to achieve the flow and temperature regimes recommended by the Upper Colorado River Endangered Fish Recovery Program (Recovery Program) [40]. One of the goals of these environmental regulations is to protect critical nursery habitats for endangered fish in the Green River downstream from Jensen.

In order to comply with the 0.1 m stage change restriction at the Jensen gage, both BoR and WAPA have been relying on SSARR in their river simulation model and hydro-scheduling tool [41]. As seen above, SSARR includes a physics-based model designed to predict downstream river flow based on upstream discharges. In particular, SSARR is used by BoR and WAPA to predict Jensen flow rate based on Flaming Gorge Dam releases and Yampa flow rate. The Jensen stage level is then deduced from Jensen flow rate based on a flow-to-stage curve regularly updated by USGS [32]. However, the geometry and geology of the Green River and Yampa River have noticeably evolved since the first time SSARR was used to predict Jensen flow, in the 1990s. As a physics-based model, SSARR parameters require to be regularly updated by conducting detailed studies about the new river shape and soil conditions, which is not always feasible for economic reasons. As a result, the predictive ability of SSARR at the Jensen gage has been declining in the past decade.

By opposition, the LP and CNN models presented in this article are both empirical models and are designed to “learn”, i.e., automatically update, their predictive parameters based on the most recent hydrology time series data, which are publicly available and constantly monitored by USGS gages. The Green River system is used as a case study to assess the performance of both models and compare it to the performance of SSARR.

2.8. Datasets

Stream gage information is accessed via the National Water Information System using the instantaneous value data retrieval REST API [28]. The USGS uses a station number to categorize sites from which hydrological data are measured, included here as “USGS 0#######”. Flaming Gorge Dam release is approximated by the Green River near Greendale, UT stream gage (USGS 09234500); Yampa River by the gage at Deerlodge Park, UT (USGS 09260050); and Jensen by the Green River near Jensen, UT stream gage (USGS 09261000).

3. Results

Models are evaluated for two measures which serve specific considerations for 1D flow prediction for the Jensen, UT stream reach: prediction of (1) hourly flow, and (2) daily maximum stage change. Hourly flow prediction is the 1 h discharge series and performance demonstrates abilities of each model to forecast high-resolution 1D flow. Often, operational considerations may value flow summary metrics over hourly flow. The second measure used to compare model performance is important in the Jensen reach where daily maximum stage change is used for evaluating dam compliance with the 0.1 m stage requirement. Each measure is evaluated with four performance metrics: mean absolute error (MAE), mean squared error (MSE), coefficient of determination (R²), and the maximum error residual. Performance metrics are calculated on the test dataset.

Flow predictions from SSARR, LP, and CNN were generated for water year 2022 (see example predictions in Figure 3) to compare the performance of each method.

Table 1. Time to train and performance metrics for hourly flow and daily minimum-to-maximum flow predicted by SSARR, LP, and CNN. The best performing metric is indicated in bold. Performance metrics for hourly flow and daily minimum-to-maximum flow include: mean absolute error (MAE) in m, mean squared error (MSE) in m², coefficient of determination (R²) on a scale of 0 to 1, and the maximum error residual (m). The LP model is trained on laptop Intel Core i7-11800H with 32 GB of RAM; The CNN model is trained on a server NVIDIA A100 40 GB. Training time is the number of seconds it takes to train one complete CNN or LP model.

	Training Time (Seconds)	Hourly Prediction				Daily Minimum-to-Maximum Prediction
Model		MAE (m)	MSE (m2)	R2	Max Error Residual (m)	MAE (m)	MSE (m2)	R2	Max Error Residual (m)
SSARR	-	0.0411	0.0030	0.987	0.295	0.027	0.0014	0.669	0.180
LP	17	0.0296	0.0020	0.991	0.210	0.016	0.00056	0.856	0.130
CNN	170	0.0171	0.00056	0.998	0.145	0.015	0.00046	0.877	0.142

shows training time and performance for both evaluation metrics across the two data-driven approaches (LP and CNN) compared to the SSARR hydrologic model. The data-driven LP and CNN models outperform SSARR under both hourly and daily min-to-max performance metrics. CNN has a mean absolute error approximately 42% smaller than LP for hourly flow, but mean squared error over 70% smaller, indicating a smaller number of larger errors in the hourly forecast. Despite performing better for hourly prediction, CNN and LP have relatively similar performance for the daily min-to-max measure. Notably, LP has a smaller maximum residual than CNN under the daily min-to-max measure. It is likely this discrepancy is due to chance rather than a structural difference that allows LP to out-perform CNN for daily min-to-max stage. The LP model can be computed about ten times faster than the CNN, but both models can be trained in under five minutes.

We are also interested in the performance of each stage model over the course of the model year and at different discharge magnitudes. Monthly error distributions from each model are shown in Figure 4. Due to large differences in flow magnitudes over the water year, the distribution of relative error (absolute error divided by monthly average stage) is also shown. While absolute error is higher when the stage near Jensen is larger, there is no evident pattern in relative error over the water year. Notably, during July through November, LP and CNN significantly outperform the SSARR model with smaller errors.

Both models are developed by applying convolution to upstream hydrographs to predict Jensen flow. LP is applied by considering Jensen flow as an evolving linear combination of five unit hydrographs for each of the two upstream sources. As described in the linear programming model section, the LP model is solved to learn the optimal unit hydrographs from the training set, and these optimal unit hydrographs are used to predict the downstream river discharge in the test set. The downstream river stage is then calculated using the Jensen flow-to-stage empirical function [32]. The input parameters $a_{k,l,t}$ describing the slowly evolving change in the unit hydrograph is based on the average monthly discharge at Yampa. Figure 5 below illustrates the five optimal unit hydrographs identified by the LP model for the Greendale upstream source (USGS 09234500). Note that the average time delay from Greendale to Jensen is decreasing as the average Yampa flow level increases. As Yampa flow speed increases, it increases flow speed on the Green River which reduces water travel time from Greendale.

The CNN considers Jensen flow as a convolution of the two series produced when applying a discharge-variant lag to Greendale and Deerlodge hydrographs. The discharge-variant lag applied to upstream sources, modeled as an exponential function, would be analogous to the LP’s varying unit hydrographs at different flow magnitudes. To fit the exponential function, the maximum cross-correlation algorithm described in the Discharge-Variant Water Travel Time section is used to aggregate a sample dataset. Manual data cleansing is required to fit a clean exponential function. Figure 6 demonstrates the exponential function fit for one training fold. As expected, the water travel time increases as discharge decreases for both upstream sources.

An important consideration for water managers and modelers using empirical techniques for 1D flow prediction is the training dataset size, which impacts the accuracy of the model. To examine the effect of training dataset size on accuracy, the LP and CNN models were trained by varying the training dataset length. Since water year 2022 is used as the evaluation dataset, training datasets with 1, 2, and 3 years of data prior to water year 2022 were used to train both models (the models trained on 3 years of data are the same models presented previously). The MAE and MSE from these experiments are shown for the LP and CNN models in Figure 7. Both CNN and LP models are impacted by training dataset length, though CNN impressively has very similar performance between models trained on 2 prior years and 3 prior years.

4. Discussion

The data-driven approaches used by the LP and CNN models demonstrate considerable prediction capabilities for addressing the challenge of 1D flow forecasting. Though modelers must consider a variety of stakeholder needs when considering methods for 1D flow, a major advantage of empirical modeling is the ability to quickly relearn hydrologic relationships at a low cost when model performance inevitably degrades. Indeed, contrary to physics-based models that require an accurate and thorough description of numerous watershed characteristics, including site geometry, and evolving soil moisture, infiltration, land use, and surface friction, the data-driven approaches presented here only require time series data from constantly monitoring and relatively low-cost gages. Because of this, the adoption of such data-driven approaches for river discharge and stage prediction has significant cost savings potential for future hydrology/hydropower-related studies.

Selecting a data-driven approach is another serious consideration. This manuscript introduces two viable solutions: a linear mathematical optimization (LP) model and one that uses artificial intelligence with non-linearities (CNN). Both models can learn optimal prediction parameters based on historical flow data and they both use convolution methods when predicting output results. We demonstrate that both models have similar performances when predicting downstream river stage at the Jensen gage. The CNN model demonstrates higher prediction accuracy than the LP model at an hourly resolution. However, the difference in prediction accuracy becomes negligeable when comparing daily min-to-max results. This implies that both models exhibit similar performances when predicting whether the water level at the Jensen gage complies with the environmental 0.1 m stage change restriction. The strength of the CNN model is based on its ability to explicitly model non-linearities between the discharge and stage data of the system, owing to its neural network structure. Contrary to the CNN model, the LP model is only able to model linear relationships, which limits its predictive ability to downstream discharge based on upstream discharges. To predict downstream stage, an additional step is required with the LP model, in which predicted discharge is converted to predicted stage owing to an empirical flow-to-stage function. This implies that its ability to predict downstream stage is limited by the accuracy of the flow-to-stage function. However, the strength of the LP model is based on its ability to identify linear predictive functions (i.e., the unit hydrographs) that are mathematically guaranteed to ensure the minimum possible prediction error. Moreover, powerful state-of-the-art LP solving methods such as the interior-point method [42] enable the optimal predictive functions to be identified in just a few seconds, even when considering time series of thousands of data points. The linear predictive functions identified by the LP model also make them more suitable to be integrated in linear water management tools or hydropower scheduling models such as GTMaxSL [27].

As demonstrated, the performance of empirical models also relies primarily on the quality of the observation data used to train them. While the LP and CNN models both perform with reasonable accuracy when decreasing the training dataset size, the performance degradation is significant when only using 1 year of data. It should also be expected that when river conditions (e.g., magnitude of flow, riverbed characteristics, or river geometry) deviate from the conditions represented by the training dataset, prediction accuracy will be hampered.

In summary, the structure of the CNN model and its explicit modeling of non-linearities makes it ideal to accurately predict downstream measures at a high time resolution (hours, minutes) that may include water discharge, water level, but also other environmental measures (e.g., temperature, gas concentration). On the other hand, the formulation of LP model and its ability to identify optimal linear predictive parameters makes it ideal to be integrated in linear water management tools or hydropower scheduling models.

5. Conclusions

This paper introduces two novel empirical models to predict downstream river flow and stage based on upstream flow data. Both models are based on convolution methods and are designed to learn predictive parameters using a training dataset (measured hydrologic data). The LP model solves a linear mathematical optimization problem to identify slowly evolving unit hydrographs that describe the impact of the upstream water flows on the downstream flow. The predicted downstream flow can later be translated into a river stage via an empirical flow-to-stage function. Conversely, the CNN model learns to directly predict the downstream stage based on upstream flows owing to a non-linear neural network. When tested on a real-world case study, both models show promising results and outperform the physics-based model currently used by federal agencies on the selected river system. Further work is needed to identify the ideal size of the training period, as there is a tradeoff between having a large enough training dataset and recent historical data. Alternatively, the prediction error of the hydrologic variables could be weighted based on how recent they are.