The upper region of the South Platte river basin (Figure 1) serves as a water source for numerous communities. Colorado front range communities draw water from the South Platte and other nearby rivers to meet their water needs throughout the year. Data inputs were used from the major reservoirs (Stronita Springs, Cheesman, and Chatfield) along the South Platte, tunnels including ones that bring water across the Continental Divide (Roberts and Moffat tunnels) or move water within the drainage (homestake), and gauges representative of the many un-dammed streams (south fork of the South Platte, middle fork of Prince, and South Platte above Spinney), as well as gauges that are located below reservoirs (Eleven mile outflow, Antero outflow, below Brush creek, and Deckers), and the pipes (conduit 20 and 26) that send water from the reservoirs into Denver. Gauges were also included in the model that are outside of the South Platte drainage but are a part of the Denver Water system including South Boulder creek which supplies water to Denver and the Blue River below Dillon Reservoir which reflects releases of water into the Colorado that are therefore unavailable for transfer to the front range. Finally, gauges were included that do not directly supply water to the front range but that are connected indirectly because they are managed by Denver Water and governed by the Colorado River Compact of 1922 (MacDonnell, 2023) which sets the total amount of water that must be sent to downstream states. Transfers of water eastward across the Divide are often offset by releases into the Colorado River to comply with this compact.

Gauges were selected for the final model by examining which gauges in this region showed strong correlations to each other within this river basin (Figure 2). Some strong correlations are due to close geography. For example, the majority of the flow out of Strontia Spring reservoir flows directly into Chatfield reservoir and these gauges have a 0.969203 Pearson correlation coefficient. In contrast, a correlation of −0.126896 is seen between the Roberts tunnel which brings water from Dillon reservoir across the continental divide to be used in Denver and the inflow to Cheesman reservoir, which feeds Denver from other sources east of the divide. When water is available east of the divide, water is not transported from west of the divide.

The homestake pipeline stands out as being negatively correlated with many of the other gauges. This project takes water from the South Platte basin and delivers it further south, to the city of Colorado Springs. Similarly, two other gauges that were geographically within this region but delivered water to the town of Aurora likewise had weaker correlations.

Data from 1996 through 2021 was collected using the Colorado Division of Water Resources (https://dwr.state.co.us/Rest/) Application Programming Interface. Outlier datapoints 2 times above the 95th percentile or below the 10th percentile in flow or capacity were removed. Historical weather data was obtained from the National Centers for Environmental Information (https://www.ncdc.noaa.gov/cdo-web/datasets) for the city of Lakewood, Colorado which is close to the region of study, receives water from this river system, and has a long and continuous set of temperature readings through the period of study. Training was performed on all days available where 21 prior days are included in the model and 7 future days are predicted.

While any river flow data is likely to have seasonal variation, the data used here presents unique challenges. In a river system heavily influenced by anthropogenic activity, not all gauges follow a similar pattern of seasonality. Some gauges on unmodified streams have seasonality expected for natural rivers with increased flow during spring snowmelt. Other gauges are inactive for long periods of time, and some have absolutely no variability for the majority of the year. Some of the gauges that are largely inactive become active during spring runoff, but others are countercyclical, changing little during the spring, but becoming more active at other times of year. For example, several tunnels bring water across the continental divide and are often most active in the fall when water demands in the front range cannot be met with natural flow from streams into reservoirs.

Winter data was removed because flow is negligible at most of the gauges, however, this creates a discontinuity in the data. While careful selection of start and end dates can avoid using data that crosses this discontinuity, it is desirable to have flexibility in data handling and training of the neural network. Ideal flexibility would allow for selecting gauges, past and future window lengths, and setting weights prior to editing out data that contains such discontinuities. Under these circumstances, a final step must be performed to remove data that crosses the discontinuity. Since there are no column labels on data used for training and validation, another method is needed for identifying data that crosses this discontinuity. For this reason, a unique identifier is created by multiplying two arbitrarily selected column values. This product acts as a unique identifier of discontinuities and it is used to identify datasets containing discontinuities to be removed immediately before training.

After combining all seasons, the data was split into three sets. The first set was used for training the neural networks. A second validation set used to periodically check the network against data that had not yet been seen by the model to protect against overfitting; when validation data did not improve the fit, fitting was stopped regardless of whether the training data was still improving the quality of fit. Finally, a third set of data was set aside and not used for training or validation but rather for testing the model. All evaluation reported here is based on how well the model handled this testing data that was not used for training or for validation. All the data was Z-score normalized and the mean and standard deviation for the training set was used to normalize the training, validation, and testing data sets to avoid risk of information creeping into the validation and testing set. A final practical consideration centered on what years to include in each set. Given that later years had more gauges available, setting aside only later years for validation and testing was not practical. Instead, the training set was constructed from a mixture of older data (1996–2012) and newer data (2017–2021), validation was based on 2013–2014 and testing was based on 2015–2016. The total training data contained 4,558 elements and the validation data contained 398 elements where each element corresponds to three prior weeks of data connected to labels composed of the actual data (flow, temperature, and reservoir capacity) for the following 7 days.

Changes to water flow at some gauges are based on water managers and water users making predictions about their anticipated needs. For example, a farmer might place a water rights call when the forecast temperatures are high in upcoming days. For this reason, the forecast high temperatures for the next 4 days were included into the model. With the goal of forecasting based on current conditions, it is important to keep in mind that more accurate temperature forecasts are available now than in the past. And records of past forecasts are not available for as many regions with as many time points as are currently available. For these reasons, the past forecast data that was used was simulated. Specifically, the actual high temperatures recorded were used as a stand-in for what would have been an accurate prediction for those temperatures on previous days. With this method, any location that has accurate historical weather data can be used also as a source of simulated forecast data. This does not account for prior inaccurate forecasts that might have impacted changes to water use since these simulated forecasts are artificially perfect.

The base and per-gauge neural networks were built using a tensorflow ConvLSTM2D with data inputs arranged into a two-dimensional grid. The gauge data was spread across the x and y axes normally used as image vectors, creating an artificial image of the data. Two ConvLSTM2D layers were used with a padded 3 by 3 kernel. The first layer encoded input data and the second layer included future predictions from the input labels. These were followed by a time-distributed Dense layer.

Optimization of the ConvLSTM2D base model and each per-gauge model was done with Root Mean Squared Propagation, a learning rate of 0.003, and a clip value of 1.0. Kernel regularizers, recurrent regularizers, bias regularizers, and activity regularizer were all set to 1x10-6. A dropout and a recurrent dropout of 0.3 was used.

The classifier was based on a simple 3D CNN using a padded kernel of 3 with sigmoidal activation. The categorization model used a glorot uniform kernel initializer, an orthogonal recurrent initializer, and a bias initializer of zeros. Kernel regularizers, recurrent regularizers, bias regularizers, and activity regularizer were all set to 1x10-6. The recurrent activation was hard sigmoid, and the model was optimized using adam while minimizing losses calculated with sparse categorical cross entropy.

Training was performed in multiple steps (Figure 3), starting with the base ConvLSTM2D using all available data. Next, this base model was used to evaluate the training data a second time with the goal of determining in which cases the base model performed better than a simple persistence model. The ratio of mean average error (MAE) for the LSTM model and the persistence model was recorded and then used as a label for training the classification model. The classification models for each gauge were also built with CNNs but with a final sigmoidal activation producing a single value between 0 and 1 to describe each data input rather than predicting the future flow. Using min-max normalization of the ratio of MAE values for each individual gauge, a label between 0 and 1 was produced for each input that conveys whether that data was best described using a persistence model or the base model. A separate classification model was made for each gauge of interest (Table 1). The final gauge-specific models were trained on only data that was rated by the classifier as being modeled better by the base (convolutional LSTM) model than a persistence model. Each gauge had its own LSTM model with its weights optimized during training.

The slowest training step was creating the persistence labels which took 89 min. However, once created, these could be used for many different attempts at optimizing the per-gauge LSTM models.

Various methods are commonly used to evaluate river flow forecast accuracy. Both the Nash and Sutcliffe (1970) method (Eq 1) and the dimensionless version of the Willmott et al. (2011) technique (Eq 2) were utilized here. The shortcomings of the most common methods including these two have been documented (Legates and McCabe, 1999; Jackson et al., 2019). One particular concern that has been raised before for the Willmott method is that it is benchmarked against the variance in past data. For the type of data examined here, this is especially problematic as the past variance can be zero for gauges near reservoirs and pipes. Unlike natural streams, reservoir levels, tunnel flows, and gauges directly connected to reservoir outflows can report zero past fluctuation for substantial periods of time. The Nash method suffers from lacking a lower bound and being strongly influenced by extreme outliers. While neural networks have used Nash as a minimization function, both Nash and Willmott were unsatisfying as loss functions for multiple reasons. MAE was found to be a more satisfying minimization and evaluation criteria although both Nash:

and Willmott error analysis:

MAE was found to be a more satisfying both Nash (Eq 1) and Willmott (Eq 2) analysis of the model using the set aside test data set are provided here. While the MAE is very direct and easily understandable, one drawback is that it is not scaled. As a result, gauges that have higher average flows tend to have predictions with larger MAE values. The impact of this scaling issue on using MAE for optimization is minimal. Most of the models described here are specific to a single gauge, so minimizing them based on errors on that gauge is not problematic. In the base model, all gauges are weighted equally, and higher flow gauges may disproportionately impact the gradient descent, but this base model is not used directly for predictions. Another consequence of using MAE is that when comparing predictions on the test data it can appear that gauges with higher average flow have larger errors. One common method to solve this is to scale errors using the mean absolute scaled error (MASE) (Hyndman and Koehler, 2006). However, the MASE would be infinite or undefined when recent historical observations are all equal in value, which is common for this type of data.

Each final LSTM model used a different number of days of data (out of a possible 4558) for training. While the initial training was performed on all 4,558 data sets, the classification allowed for curating the training set to contain only that data expected to help the most. The classification method when applied to selecting data used in the final model for each gauge ranged from 2,739 to 4,524 (Figure 4) out of the possible 4,558 days available. The gauges nearer water projects used the persistence model more frequently than the natural streams.

When a 7-day prediction was run on 210 days of test data, the prediction 1 day into the future was close to that of the persistence model (Figure 5). In cases where there were rapid changes in flow, typically a rise, the convolutional LSTM model often was slower to show the increase. The convolutional LSTM models had mixed results when reservoirs or pipelines had long periods without any change. In some cases, such as below Dillon reservoir, this model was very similar to the persistence model. In other cases, such as homestake pipeline and cheesman reservoir, the convolutional LSTM model predicted flow at large sections of time when there was no flow recorded. Natural undammed streams such as below brush were reasonably accurate although the convolutional LSTM model was slower to respond to rapid changes. When run on the test data, the number of weeks labeled as being better modeled with persistence ranged from 98 to 133 out of the 210 days evaluated (yellow points in Figure 5).

In the same 7-day prediction, the prediction 7 days into the future was not unexpectedly, less accurate than 1 day into the future (Figure 6). Again, when changes in flow were rapid, the convolutional LSTM model often was slower and did not rise as high as the actual flow (blue) or the persistence prediction (red). The convolutional LSTM models again overestimated flow at times when reservoirs or pipelines had long periods without any change such as in homestake, underestimated at other times as in Roberts tunnel, and identified small peaks below cheesman reservoir when the actual data had no such spike. Although the reservoir was undergoing repairs during this time period that may have caused it to deviate from the model prediction.

The fully trained gauge specific models utilized approximately half of the available data provided from the classifier. These models trained on smaller datasets performed well at predicting flow from all the gauges. The median total error (sum of the absolute value of the difference between true value and predicted value for 7 future days combined) ranged from 0 cfs to as high as 413.7 cfs (Figure 7). The models with the higher flow gauges showed larger errors, as expected.

When analyzing the test data (Figure 3, right) using this combination of classification CNN and convolutional LSTM predictor, the median Willmott values for all gauges ranged from 0.801 to 1.0 (Figure 8). Here the lowest values were for pipeline and reservoir gauges rather than natural river and stream flows.

The median Nash values for each gauge ranged from 0.65 to −5.39 with gauges close to reservoirs showing worse performance. Here the Nash values are graphed as a boxplot after removing outliers one standard deviation beyond the IQR (Figure 9).

Predictions were also made using the base model to allow for an analysis of the benefit of including the active learning and transfer learning into the overall model. In almost all cases, the final model performed better than the base model (Table 2) despite both being based on convolutional LSTM architecture, demonstrating that this approach does add to the forecast accuracy. The total error (over 7 days) was 11 to 196 cubic feet per second (cfs) better for the final model, all but two gauges had a better Willmott score, and all but 5 had a better Nash-Sutcliffe score. The propensity of Nash-Sutcliffe to overweight outliers likely impacted this result.