A DeepAR-Based Modeling Framework for Probabilistic Mid-Long Term Streamflow Prediction

1. Introduction

The combined impacts of climate change and human activities are profoundly altering the processes of runoff generation and confluence, resulting in increasing uncertainty in the evolution of water resources [1,2,3,4]. At the same time, socio-economic drivers such as population growth and urbanization are intensifying the demand for water supply [5,6,7,8]. Against this backdrop, mid-long term streamflow prediction (MLSP) has become increasingly vital for effective water resources management and integrated utilization, as it provides valuable insights into future runoff patterns [9,10,11,12,13,14]. Consequently, MLSP is attracting growing attention in both research and practical applications [11,12,13,14].

Many models have been developed and applied in MLSP to improve the predictive performance, and provide information for the comprehensive utilization of water resources [11,12,13,14,15,16,17]. These models can be broadly divided into physical-based models, which simulate the streamflow based on the runoff generation and confluence equations, and data-driven models, which directly simulate the relationship between streamflow and predictors including precipitation, temperature and other factors [13,14,18,19,20,21]. Along with the development of artificial intelligence (AI) methods, the AI-based data-driven models, including support vector regression (SVR), artificial neural network (ANN), gated recurrent unit neural network (GRU), long short-term memory network (LSTM) and so on, can obtain better predictive performance than traditional models and have become predominant in the MLSP [10,19,22,23,24,25,26]. For instance, many studies applied SVR models in MLSP and the results demonstrate that the SVR models can generate more accurate forecasts than linear models and ANN models [27,28,29,30,31]. Xie et al. (2024) compared five AI-based models and the results demonstrate that the LSTM model outperformed other models in forecasting monthly streamflow in 37 basins [14].

The proposed AI-based models in MLSP demonstrate significant improvements in predictive performance, but their limited capacity to characterize future water resource uncertainties constrains practical applications [32,33,34,35]. To overcome this issue, many post-processing methods are adopted to produce ensemble forecasting results capable of characterizing predictive uncertainties [34,35,36]. For example, Liang et al. (2018) proposes the hydrological uncertainty processor to post-process the deterministic outputs from the SVR model to quantify prediction uncertainties [33]. Mo et al. (2023) applies the generalized autoregressive conditional heteroskedasticity model to identify time-varying forecasting errors to improve the predictive performance [37].

Although the post-processing approach in MLSP demonstrates predictive capability, it exhibits two fundamental limitations: (1) inability to directly generate probability distributions, and (2) failure to preserve the inherent statistical characteristics of streamflow. To overcome these constraints in MLSP, this study adopts a DeepAR-based modeling framework. The DeepAR architecture exhibits two key capabilities: (1) direct prediction of the predictand's distribution parameters, and (2) flexible selection of probability distributions for the target variable [38,39]. Empirical validation across five time series datasets demonstrates its superior performance over other existing state-of-the-art methods. However, the DeepAR's robustness under heavy-tailed streamflow distributions requires further verification, and the selection criteria for appropriate probability distributions lack systematic guidance. Therefore, the objectives of this study are (1) to develop a DeepAR-based probabilistic forecasting framework for MLSP, (2) to validate modeling framework's applicability through implementation in two case studies, and (3) to systematically evaluate impacts of base model architecture and distribution type selection on predictive skill.

The remaining sections of this paper are organized as follows. The data, case studies and methods are introduced in section 2. The results will be demonstrated in section 3 and discussed in section 4. Finally, the main conclusions will be summarized in section 5.

2. Materials and Methods

2.1. DeepAR Model

DeepAR is a probabilistic forecasting framework developed by Amazon Research, which combines recurrent neural networks (RNNs) with parametric probability distributions to generate time-series predictions [38]. Unlike traditional point-forecasting models, DeepAR directly outputs the parameters of user-specified distributions (e.g., Gaussian for real-valued data, Negative Binomial for positive count data, Beta for data in the unit interval), enabling native uncertainty quantification.

Let $z_{\mathrm{t}}$ denotes the streamflow value at time t, and ${\mathit{x}}_{\mathit{t}}$ represents the vector of predictor variables at time t. The DeepAR model estimates the conditional distribution of future streamflow values:

$$P\left({\mathbf{z}}_{t_0+1:t_0+L}\right|{\mathbf{z}}_{t_0-H+1:t_0},{\mathit{x}}_{t_0-H+1:t_0+L})$$

(1)

where $t_0$ is the forecast initialization time, $H$ is the context length and $L$ is the prediction length.

For a trained DeepAR model, the distribution parameters ${\mathit{\theta }}_t$ at time t are computed as a function of the hidden state ${\mathit{h}}_{\mathit{t}}$ and model parameters $\mathit{\Theta }$:

$${\mathit{\theta }}_{\mathit{t}}=\mathrm{f}({\mathit{h}}_{\mathit{t}},\mathit{\Theta })$$

(2)

where $\mathrm{f}(·)$ is a function used to map the hidden state to distribution parameters, and the hidden state ${\mathit{h}}_{\mathit{t}}$ evolves recursively via:

$${\mathit{h}}_{\mathit{t}}=\mathrm{g}({\mathit{h}}_{\mathit{t}-1},z_{\mathrm{t}-1},{\mathit{x}}_{\mathit{t}},\mathit{\Theta })$$

(3)

where $\mathrm{g}(·)$ is a nonlinear transition function implemented as a multi-layer RNN (LSTM or GRU) parametrized by $\mathit{\Theta }$, and ${\mathit{h}}_{\mathit{t}-1}$ is the hidden state from the previous time step.

Then the simulated or forecasted streamflow value at time t can be sampled by:

$$\widetilde{z_t}~p(·|\mathrm{f}({\mathit{h}}_{\mathit{t}},\mathit{\Theta }))$$

(4)

2.1.1. Training

Given a streamflow time series ${\left\{{\mathbf{z}}_{\mathbf{t}}\right\}}_{\mathrm{t}=\mathrm{1,2},\dots ,\mathrm{T}}$ and associated predictor variables ${\mathbf{x}}_{\mathbf{t}}$, the DeepAR model parameters $\mathit{\Theta }$, including the parameters of both the $\mathrm{f}\left(·\right)$ and $\mathrm{g}(·)$, can be learned by maximizing the log-likelihood as below:

$$\mathcal{L}=\sum _{t_0}\sum _{t=t_0+1}^{t_0+L}\log p\left(z_t\right|\mathrm{f}({\mathit{h}}_{\mathit{t}},\mathit{\Theta }))$$

(5)

2.1.2. Prediction

Given an observed streamflow sequence ${\mathbf{z}}_{t_0-H+1:t_0}$ and corresponding predictor variables ${\mathit{x}}_{t_0-H+1:t_0+L}$, the trained DeepAR model generates probabilistic forecasts for future streamflow values ${\mathit{z}}_{t_0+1:t_0+L}$ through the following procedure:

1) The hidden state ${\mathit{h}}_{{\mathit{t}}_0}$ is obtained by recursively processing the historical streamflow ${\mathit{z}}_{t_0-H+1:t_0}$ and predictors ${\mathit{x}}_{t_0-H+1:t_0}$ through the RNN transition function in Equation (3);

2) Initial conditions are set as $\stackrel{~}{{\mathit{h}}_{{\mathit{t}}_0}}={\mathit{h}}_{{\mathit{t}}_0}$ and $\stackrel{~}{z_{{\mathrm{t}}_0}}=z_{{\mathrm{t}}_0}$;

3) For each subsequent time step $t=t_0+1$ to $t_0+L$, the hidden state $\widetilde{{\mathit{h}}_{\mathit{t}}}$ is updated using $\mathrm{g}(\widetilde{{\mathit{h}}_{\mathit{t}-1}},\widetilde{z_{\mathrm{t}-1}},{\mathit{x}}_{\mathit{t}},\mathit{\Theta })$, the forecast $\widetilde{z_t}$ is sampled from $p(·|\mathrm{f}(\widetilde{{\mathit{h}}_{\mathit{t}}},\mathit{\Theta }))$, and

4) Step 3) is repeated N times to produce an ensemble forecasts ${\left\{{\mathbf{z}}_{\mathbf{i},{\mathit{t}}_0+1:{\mathit{t}}_0+\mathit{L}}\right\}}_{\mathrm{i}=\mathrm{1,2},\dots ,\mathrm{N}}$, providing a Monte Carlo approximation of the predictive distribution.

2.1.3. Likelihood Model

The likelihood $p(·|\mathrm{f}({\mathit{h}}_{\mathit{t}},\mathit{\Theta }))$, which defines the target distribution, must be carefully selected to match the statistical characteristics of predictand. While Salinas et al. (2020) recommended Gaussian distributions for real-valued data and Negative Binomial distributions for positive count data, these choices prove suboptimal for streamflow forecasting due to: 1) heavy-tailed characteristics of streamflow, and 2) the strictly positive nature of streamflow values. Therefore, to better characterize the statistical properties of streamflow, this study employs Student's t-distribution which is suitable for extreme events, Gamma distribution which is suitable for skewed streamflow, and conventional Gaussian distribution.

For Student's t-distribution, the likelihood and function $\mathrm{f}(·)$ are as below:

$$\begin{gathered} p_{ST}\left(z|\mu ,\sigma ,\nu \right)={\displaystyle \frac{\mathsf{\Gamma }\left(\frac{\nu +1}{2}\right)}{\mathsf{\Gamma }\left(\frac{\nu }{2}\right)\sqrt{\nu \pi }\sigma }}{\left(1+{\displaystyle \frac{{\left(z-\mu \right)}^2}{\nu {\sigma }^2}}\right)}^{-{\displaystyle \frac{\nu +1}{2}}} \\ \mu \left({\mathit{h}}_{\mathit{t}}\right)={{\mathit{w}}_{\mathit{\mu }}^{\mathit{T}}\mathit{h}}_{\mathit{t}}+b_{\mu } \\ \sigma \left({\mathit{h}}_{\mathit{t}}\right)=\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{p}\mathrm{l}\mathrm{u}\mathrm{s}\left({{\mathit{w}}_{\mathit{\sigma }}^{\mathit{T}}\mathit{h}}_{\mathit{t}}+b_{\sigma }\right) \\ \nu \left({\mathit{h}}_{\mathit{t}}\right)=2+\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{p}\mathrm{l}\mathrm{u}\mathrm{s}\left({{\mathit{w}}_{\mathit{\nu }}^{\mathit{T}}\mathit{h}}_{\mathit{t}}+b_{\nu }\right) \end{gathered}$$

(6)

where the softplus activation function $\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{p}\mathrm{l}\mathrm{u}\mathrm{s}\left(x\right)=\log (1+e^x)$ is applied to enforce positivity constraints on the distribution parameters.

For Gamma distribution, the likelihood and function $\mathrm{f}(·)$ are as below:

$$\begin{gathered} p_{Gamma}\left(z|\alpha ,\beta \right)={\displaystyle \frac{{\mathsf{\beta }}^{\mathsf{\alpha }}}{\mathsf{\Gamma }\left(\alpha \right)}}{z}^{\alpha -1}{e}^{-\beta z} \\ \alpha \left({\mathit{h}}_{\mathit{t}}\right)=\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{p}\mathrm{l}\mathrm{u}\mathrm{s}\left({{\mathit{w}}_{\mathit{\alpha }}^{\mathit{T}}\mathit{h}}_{\mathit{t}}+b_{\alpha }\right) \\ \beta \left({\mathit{h}}_{\mathit{t}}\right)=\mathrm{s}\mathrm{o}\mathrm{f}\mathrm{t}\mathrm{p}\mathrm{l}\mathrm{u}\mathrm{s}\left({{\mathit{w}}_{\mathit{\beta }}^{\mathit{T}}\mathit{h}}_{\mathit{t}}+b_{\beta }\right) \end{gathered}$$

(7)

2.2. DeepAR-Based Modeling Framework

This study develops a DeepAR-based probabilistic forecasting framework for MLSP (Figure 1), which consists of four main steps: (1) data preparing, involving data collection and optimal distribution selection, (2) data splitting, which spites data into training, validation and test datasets, (3) model calibration, involving model training and selection, and (4) model evaluation, involving model performance assessment [40,41].

2.2.1. Data Preparing

Data preparing and analysis are proposed to collect streamflow and corresponding predictors (precipitation), and select optimal distribution for the streamflow. To select the most appropriate distribution from these candidates, both the Akaike Information Criterion (AIC) and Bayesian Information Criterion (BIC) are employed to evaluate the goodness-of-fit between the observed data and theoretical distributions. Then the optimal distribution is selected by minimizing both AIC and BIC values, which are computed as:

$$\begin{gathered} \mathrm{A}\mathrm{I}\mathrm{C}=-2\mathrm{l}\mathrm{n}\left(Li\right)+2k \\ \mathrm{A}\mathrm{I}\mathrm{C}=-2\mathrm{l}\mathrm{n}\left(Li\right)+k\mathrm{l}\mathrm{n}\left(n\right) \end{gathered}$$

(8)

where Li is the maximized likelihood value, k denoted parameters number, and n is length of streamflow time series.

2.2.2. Data Splitting

Data splitting is an important step in data-driven modeling process, through which the available data is divided into training, validation and test datasets [40,41,42]. In this study, the data after a specific time point is first separated as the test set to ensure no overlap between the test set data and other data. Subsequently, the remaining data is randomly shuffled and split into training and validation datasets in an specific ratio.

2.2.3. Model Calibration

Model calibration is used to optimize model architecture and parameters in order that the model can represent the underlying relationships between predictors and predictand. First, the optimal streamflow distribution obtained by minimizing the AIC/BIC metrics is adopted to define the model architecture (i.e. model’s output distribution type). Then, multiple model variants are generated by varying input conditions and parameters are optimized using training dataset by the ADAM optimizer [43]. Finally, the model variants are used in validation dataset and compared to select the best model in terms of their predictive performance.

2.2.4. Model Evaluation

Model evaluation is used to assess the predictive performance of the selected model over an independent dataset (i.e. test dataset) [44]. The deterministic and probabilistic predictions are generated and two metrics are proposed to evaluate the predictive performance: 1) root mean square error (RMSE), quantifying the accuracy of ensemble mean predictions against observations; and 2) continuous ranked probability score (CRPS), which measures the overall probabilistic predictive performance [45,46]. These two metrics can be calculated according to the following equations:

$$\begin{gathered} \mathrm{R}\mathrm{M}\mathrm{S}\mathrm{E}=\sqrt{{\displaystyle \frac{1}{n}}{\sum }_{i=1}^n{\left(\widehat{Q_i}-Q_i\right)}^2} \\ CRPS={\displaystyle \frac{1}{N}}{\sum }_{i=1}^N\int {\left(F_i\left(x\right)-H\left(x-Q_i\right)\right)}^{2}dx \\ H\left(x-Y_i\right)=\left\{\begin{array}{c}\begin{array}{cc}0& x<Y_i\end{array}\\ \begin{array}{cc}1& x\ge Y_i\end{array}\end{array}\right. \end{gathered}$$

(9)

where n is the sample size in test dataset, i is the sample index, $Q_i$ is the observed streamflow, $\widehat{Q_i}$ is the ensemble mean prediction, $F_i\left(\right)$ is the cumulative distribution function (CDF) of the probabilistic forecast.

2.3. Case Study and Data

The performance of the DeepAR-based modeling framework is examined in two basins upstream of the Wudongde (WDDR) and Sanxia (SXR) reservoirs, as illustrated in Figure 2. The available data and statistical characteristics for the two basins is presented in Table 1, including ten-day naturalized streamflow and areal mean precipitation records spanning January 1980 to September 2022. To ensure data compatibility with neural network requirements, all variables are first normalized before model processing, and then inversely transformed to their original scales.

2.4. Experiment Setup

Following the DeepAR-based modeling framework, three candidate distribution types—Gaussian, Gamma, Student’s t-distribution—are provided to account for streamflow characteristics during the data preparing stage. Then, the data after January 2017 is separated as the test dataset, while the other data is randomly split into training and validation datasets at an 8:2 ratio.

During the model calibration phase, different models for predicting streamflow in the next 18 ten-day periods, are trained based on varying input conditions that incorporate combinations of three precipitation input scenarios (temporal lags: [0], [0,1] or [0,1,2]) and two streamflow input scenarios (temporal lags: [1], [1,2]). The optimal input combination is then selected based on comparative evaluation of all model variants' predictive performance (i.e. RMSE) on the validation dataset.

After data preparing and model calibration, the final model combining the optimal input configuration and the probability distribution output is established and applied to produce deterministic predictions and probabilistic predictions with 100 members in this study. In order to evaluate the impact of output distribution and RNN structure, alternative models with different probability distribution outputs and different RNN structures (LSTM and GRU), named GRU-N, GRU-S, GRU-G, LSTM-N, LSTM-S, LSTM-G, are established and compared in terms of their predictive performance on the testing dataset. In the “GRU-N” naming structure, the former denotes the RNN structure (including GRU and LSTM), while the latter denotes the model output distribution (including normal distribution N, Student's t-distribution S, and Gamma distribution G).

3. Results

Following the modeling framework, the results are presented in three sections: optimal probability distribution selection (Section 3.1), input configuration optimization (Section 3.2), and testing performance evaluation (Section 3.3).

3.1. Optimal Probability Distribution Selection

The distribution fitting results for both study areas are presented in Figure 3, while the corresponding AIC and BIC values are listed in Table 2. The results reveals that the Gamma distribution provides superior statistical performance compared to Normal and Student's t-distributions. For the Upper WDDR area, the Gamma distribution achieves the lowest AIC (27321.09) and BIC (27337.03) values, substantially outperforming the Normal distribution (AIC: 28574.36, BIC: 28584.99) and Student's t-distribution (AIC: 28455.62, BIC: 28471.57). Similarly, in the Upper SXR area, the Gamma distribution exhibits superior performance with AIC (30848.71) and BIC (30864.66) values considerably lower than the alternative distributions. The probability density plots also illustrate that the Gamma distribution more accurately captures the right-skewed characteristics and tail behavior of the streamflow data in both study areas.

3.2. Input Configuration Optimization

The comparative analysis of different input configurations reveals consistent performance patterns across both study areas. For the Upper WDDR area, the precipitation input scenario with temporal lags [0,1,2] achieves the lowest RMSE (1199.21 m³/s), indicating superior predictive accuracy when incorporating precipitation data from the current time step and two previous time steps. In contrast, the precipitation-only scenario with no temporal lag [0] shows the poorest performance with the highest RMSE (1277.07 m³/s). Similarly, in the Upper SXR area, the [0,1,2] precipitation configuration demonstrates optimal performance with an RMSE value of 3481.18 m³/s, while the [0] configuration shows the highest RMSE (3634.35 m³/s). The performance gradient follows a consistent pattern across both areas, where increased temporal lag information progressively improves model accuracy.

Regarding streamflow input configurations, the comparison between temporal lags [1] and [1,2] shows that incorporating additional historical streamflow information ([1,2]) yields only marginal improvements, with RMSE decreasing modestly from 1245.33 to 1230.19 m³/s in the Upper WDDR area and from 3577.72 to 3577.31 m³/s in the Upper SXR area. This suggests that the contribution of additional streamflow lag information is relatively limited compared to the substantial performance gains observed with precipitation temporal lags.

3.3. Testing Performance Evaluation

The predictive performance metrics of the six models (i.e. GRU-N, GRU-S, GRU-G, LSTM-N, LSTM-S, LSTM-G) are presented in Table 4. Models with Gamma distribution output achieve the best performance in both study areas and evaluation metrics. In the Upper WDDR area, LSTM-G demonstrates the lowest RMSE (1016.54 m³/s) and CRPS (473.26 m³/s), followed by GRU-G with RMSE of 1098.98 m³/s and CRPS of 517.54 m³/s. For the Upper SXR area, LSTM-G also shows superior performance with RMSE of 4047.15 m³/s and CRPS of 1717.93 m³/s. The differences between GRU and LSTM architectures are relatively small and show no consistent patterns. In Upper WDDR, LSTM-G slightly outperforms GRU-G, while in Upper SXR, the differences are marginal. This suggests that RNN architecture choice has minimal impact compared to output distribution selection.

4. Discussion

The predictive performance of the six models (i.e. GRU-N, GRU-S, GRU-G, LSTM-N, LSTM-S, LSTM-G) is discussed in this section. First, Section 4.1 discusses the deterministic prediction performance (i.e., RMSE) of models with different distribution outputs and different RNN structures across various forecast horizons. Then, Section 4.2 provides a comparative analysis of their probabilistic prediction performance (i.e., CRPS). Finally, Section 4.3 demonstrates the overall predictive performance of the models across different forecast horizons.

4.1. Deterministic Prediction Performance of Different Models

The deterministic predictive performance (i.e., RMSE) of six models (i.e., GRU-N, GRU-S, GRU-G, LSTM-N, LSTM-S, LSTM-G) across 18 forecast horizons (10-day periods) for two study areas is presented in Table 5 and Figure 4. It is evident that the forecasting accuracy is generally observed to decline as the forecast horizon increases, regardless of model architecture or output distribution. Among the three distributions, the Gamma distribution consistently results in lower RMSE values, particularly at longer lead times. This suggests a better capacity of the Gamma-based models to capture the positively skewed or heteroscedastic characteristics often found in hydrological data. In the Upper WDDR area, the LSTM-G model achieves the lowest RMSE across most forecast horizons, with values ranging from 859.8 m³/s at the 1st forecast period to 1096.7 m³/s at the 18th. A similar pattern is observed in the Upper SXR area, where LSTM-G yields the minimum RMSE of 3687.7 m³/s at the 1st period and 3998.1 m³/s at the 18th.

Differences between LSTM and GRU structures are also observed but are found to be less consistent and less influential than those resulting from output distribution selection. In the Upper WDDR area, GRU-N and GRU-S outperform their LSTM counterparts, while the advantage shifts to LSTM only when the Gamma distribution is employed. In the Upper SXR area, the superiority of LSTM-G is observed only at shorter forecast horizons, with GRU-G performing better as the forecast horizon increases. Across all configurations, no consistent advantage is associated with either architecture, suggesting that model structure plays a secondary role relative to the output distribution in determining predictive performance.

Notably, the magnitude of performance gains brought by the Gamma distribution varies across regions. In the Upper WDDR area, the use of Gamma distribution yields substantial improvements compared to Normal and Student’s t assumptions, particularly when paired with LSTM. In contrast, in the Upper SXR area, although Gamma-based models still outperform others, the improvement is relatively marginal. For instance, the reduction in RMSE from LSTM-N to LSTM-G at the 18th horizon is only modest (from 4203.0 m³/s to 4070.5 m³/s), whereas the gain is more pronounced in the Upper WDDR region (from 1559.5 m³/s to 1096.7 m³/s). These findings demonstrate that the benefit of applying flexible, non-Gaussian output distributions such as Gamma is region-dependent and may vary with underlying hydrological complexity or data characteristics.

In summary, the results highlight the critical role of output distribution in deterministic streamflow forecasting, with the Gamma distribution consistently offering performance benefits, though to varying extents across regions. While model structure has some impact, it exerts less influence than distributional assumptions.

4.2. Probabilistic Prediction Performance of Different Models

The probabilistic predictive performance (i.e., CRPS) of six models (i.e., GRU-N, GRU-S, GRU-G, LSTM-N, LSTM-S, LSTM-G) across 18 forecast horizons (10-day periods) for two study areas is presented in Table 6 and Figure 5. Similar to the deterministic results, probabilistic forecast accuracy generally declines with increasing lead time, as reflected by rising CRPS values across all model configurations. Among the output distributions, the Gamma distribution again provides the most substantial improvements in forecast skill. This consistency between deterministic and probabilistic evaluations reinforces the robustness of the Gamma assumption in capturing the inherent characteristics of streamflow data. In the Upper WDDR area, LSTM-G yields the lowest CRPS across nearly all horizons, ranging from 340.4 m³/s at the 1^st period to 521.4 m³/s at the 18^th. In the Upper SXR area, LSTM-G performs best during the first three horizons, while GRU-G shows superior accuracy at longer horizons. This pattern echoes deterministic results, where GRU-G also performed better in long-term forecasting under complex conditions.

Compared to output distribution selection, differences due to model structure are less consistent and generally less influential. While LSTM tends to perform better in the Upper WDDR area, GRU offers competitive or better results in the Upper SXR region when combined with the Gamma distribution.

Regional differences are again observed, with CRPS values in the Upper SXR region consistently higher than those in the Upper WDDR area, reflecting greater predictive uncertainty. However, the relative benefit of using the Gamma distribution persists in both regions, though with different magnitudes. For instance, at the 18^th horizon, the CRPS reduction from GRU-N to GRU-G is more substantial in the Upper WDDR area (from 705.2 m³/s to 572.2 m³/s) than in the Upper SXR area (from 1843.3 m³/s to 1699.6 m³/s), consistent with patterns seen in deterministic forecasts.

In summary, probabilistic forecasting results confirm key findings from deterministic evaluations—particularly the consistent advantage of using the Gamma distribution across regions and horizons. However, nuanced differences are observed in the relative contributions of model structure, especially under probabilistic metrics. These findings suggest that while distributional assumptions remain the most critical component for improving streamflow forecast quality, model structure and regional hydrological characteristics jointly shape both deterministic and probabilistic forecasting performance.

4.3. Overall Predictive Performance

The predicted and observed streamflow at the 1^st forecast horizon are shown in Figure 6. It is evident that all models produce narrow prediction intervals that closely follow the observations during non-flood seasons, when streamflow is low and uncertainty is limited. In contrast, flood seasons are characterized by increased variability, leading to wider forecast intervals and larger deviations, particularly around peak flows. The choice of output distribution also significantly affects forecast reliability: models using the Normal distribution often generate overly wide intervals due to their poor fit to the skewed nature of streamflow, whereas Gamma-based models produce tighter and more stable intervals. While Gamma-based models generally provide better predictions in most cases, their advantage may not hold in every situation—for example, the GRU-G model perform worst in the Upper WDDR region. Nonetheless, their ability to constrain uncertainty remains evident, and such model-specific variability also underscores the importance of model fusion strategies for improving forecast robustness under diverse application scenarios [20,47,48,49,50].

The predicted and observed streamflow across 18 forecast periods for models using Gamma distribution in both the Upper WDDR and Upper SXR regions are presented in Figure 7. As the forecast horizon increases, streamflow uncertainty becomes more pronounced, and the width of the predictive intervals correspondingly expands. This widening of intervals effectively reflects the growing uncertainty associated with longer lead times, particularly during high-flow periods. The use of the Gamma distribution enables the models to better capture the skewed nature of streamflow, resulting in predictive intervals that not only reflect the asymmetry of the data but also reliably encompass the observed hydrographs across seasons and forecast periods. This demonstrates the capability of using suitable output distribution to improve both the reliability and calibration in probabilistic streamflow prediction.

5. Conclusions

In this study, a DeepAR-based modeling framework is developed to generate probabilistic streamflow forecasts by integrating distribution selection with the DeepAR model. The framework is applied to two case studies to evaluate both deterministic and probabilistic forecasting performance. The influence of model structure and output distribution choice on prediction accuracy is also examined. The main conclusions are as follows.

The proposed DeepAR-based modeling framework effectively identifies the most appropriate output distribution and input configuration, leading to superior deterministic and probabilistic forecasting performance in both case studies. By integrating distribution selection and input optimization, the framework ensures that the final models are well-calibrated to the characteristics of local streamflow dynamics.

Multiple factors influence predictive performance, among which the selection of output distribution has the most significant impact. Meanwhile, differences between RNN model structures (i.e., GRU vs. LSTM) are relatively minor and less consistent. Additionally, regional differences in forecasting performance suggest that the effectiveness of a given model configuration is influenced by underlying hydrological conditions.

As forecast lead time increases and during flood seasons, streamflow uncertainty grows substantially. The developed probabilistic forecasting framework is capable of capturing this variation in uncertainty, and models using the Gamma distribution demonstrate superior performance by better representing the skewed nature of streamflow. This contributes to more reliable and better-calibrated forecasts across both typical and high-variability conditions.

Given that different models exhibit varying performance under different scenarios, future work could explore model fusion approaches to leverage complementary model strengths and enhance forecast robustness. Moreover, the framework’s generalizability could be further validated by testing a wider range of output distributions—tailored to characteristics of daily or annual streamflow and other hydro-meteorological variables—thereby expanding its applicability across diverse forecasting contexts.