Enhancing hydrological hazard early warning: a 60 d streamflow forecasting framework integrating deep learning and process-based modeling

Reliable medium- and long-term streamflow forecasting is a cornerstone of hydrological hazard early warning and water resources management, yet achieving accurate predictions with sufficient lead time remains a formidable challenge. This study proposes a 60 d streamflow forecasting framework to strengthen early warning capabilities by systematically integrating a convolutional neural network (CNN) for bias correction of precipitation forecasts from the UK Met Office (UKMO) numerical weather prediction model, the Geomorphology-Based Eco-Hydrological Model (GBEHM) for streamflow simulation, and an autoregressive with exogenous input (ARX) model for statistical post-processing. Applying the proposed framework to the Upper Yangtze River Basin, results indicate that the CNN model reduces the areal-averaged precipitation root mean square error (RMSE) by around 35 % and elevates the temporal correlation coefficient (TCC) from 0.62 to 0.74 against raw UKMO forecasts across the 60 d horizon, with performance gains amplifying at longer lead times. Subsequently, when driving the GBEHM with corrected precipitation and applying ARX post-processing, the streamflow forecasts exhibit substantial enhancements with a reduction in RMSE of 36 %, a decrease in relative error (RE) from 48.2 % to 17.4 %, and an increase in Nash–Sutcliffe efficiency (NSE) from 0.33 to 0.72 compared to those driven by raw forecasts in terms of 60 d mean performance. Error decomposition identifies precipitation forecast errors which intensify with lead time as the dominant source of uncertainty for medium- and long-term streamflow forecasting, while confirming that hydrological model uncertainty remains a significant component, highlighting that the selection of a robust hydrological model is crucial for enhancing the reliability and predictive skill of the streamflow forecasts. By systematically leveraging the CNN to mitigate drifting meteorological biases, the GBEHM to capture physical catchment dynamics, and the ARX to minimize residual errors, the proposed framework extends the effective early warning horizon to 60 d with high volumetric accuracy and temporal consistency, providing vital decision support for flood and drought risk management and regional water security.

1 Introduction

The escalating frequency and intensity of extreme hydrological events, such as devastating floods and prolonged droughts, present significant challenges to global water security and socioeconomic stability within the context of climate change (Swain et al., 2025; Sutanto et al., 2025; Kreibich et al., 2022; Tabari, 2020). Reliable medium- and long-term streamflow forecasting, typically referring to predictions with lead times ranging from 3 d to 1 year (Shao et al., 2024), is recognized as a critical measure for hydrological hazard early warning by providing more substantial lead times for proactive decision-making in flood and drought mitigation, hydropower optimization and water supply security (Kondal et al., 2024; Lee et al., 2022; Jackson-Blake et al., 2022; Xu et al., 2020). However, constrained by the inherent complexity of hydrological processes and atmospheric uncertainties, the accuracy and lead time of current forecasts remain insufficient to fully meet the rigorous demands of disaster early warning (Slater et al., 2023; Neri et al., 2020). Consequently, enhancing the predictive skill and extending the reliable horizon of streamflow forecasts are imperative for strengthening regional disaster resilience and safeguarding water security (Koh and Galelli, 2024; Sutanto et al., 2020; Pendergrass et al., 2020).

Medium- and long-term streamflow forecasting methods are typically divided into three main categories: time series analysis, data-driven methods and dynamical methods. Time series analysis forecasts streamflow by extrapolating statistical patterns solely from historical observations (Wang et al., 2020; Huang et al., 2020; Luo et al., 2019). For instance, Guo et al. (2023) applied a hybrid time-series framework based on signal decomposition to forecast medium- and long-term streamflow in the lower Yellow River. However, inherently limited to coarser temporal scales and devoid of atmospheric drivers, this univariate approach struggles to capture both daily streamflow variability and dynamic hydrological responses to climate change (Jamei et al., 2024; Kratzert et al., 2018). Data-driven methods forecast streamflow by establishing direct quantitative mappings from meteorological predictors (Wang et al., 2023; Hunt et al., 2022; Adnan et al., 2020; Yuan et al., 2024; Zhou et al., 2024). For example, Cheng et al. (2020) utilized long short-term memory (LSTM) networks to map precipitation and antecedent flow to future streamflow for long lead-time forecasting. However, lacking explicit representation of physical hydrological processes, these data-driven methods are fundamentally limited in capturing non-stationary environmental impacts thereby hindering the prediction of unprecedented hydrological hazards (Reichstein et al., 2019; Shen, 2018). Dynamical methods generate streamflow forecasts by driving process-based hydrological models with meteorological predictions (Lee et al., 2024; Greuell and Hutjes, 2023; Quedi and Fan, 2020). For instance, Tian et al. (2018) employed the THREW model driven by bias-corrected European Centre for Medium-Range Forecasts (ECMWF) precipitation forecasts for seasonal prediction in the Upper Hanjiang River Basin. Among these three categories, dynamical methods leverage their explicit physical mechanisms and interpretability to offer superior transferability and the capability to capture climate change impacts, thereby making them indispensable for robust hydrological hazard early warning (Falck et al., 2025; Zhang et al., 2024; Andrade et al., 2024).

As the primary meteorological forcing, precipitation forecasts play a decisive role in determining the accuracy of medium- and long-term streamflow forecasts (Ghimire et al., 2021). Currently, medium- and long-term precipitation forecasts rely predominantly on numerical weather prediction models (Adams III and Dymond, 2019; Bauer et al., 2015). However, these raw outputs typically exhibit large discrepancies and systematic biases (Siqueira et al., 2020; Lin et al., 2019). As the forecast lead time extends, their reliability diminishes significantly, making them unsuitable for direct application in hydrological hazard assessment (Du et al., 2025; Bogner et al., 2022). Consequently, bias correction is essential to render these outputs usable for streamflow forecasting (Vernon et al., 2025; Monhart et al., 2019; Anghileri et al., 2019). Traditional correction methods, such as quantile mapping, have been shown to effectively correct systematic biases in precipitation forecasts (Cannon et al., 2015; Teutschbein and Seibert, 2012). Recently, deep learning approaches have demonstrated superior performance by leveraging their ability to model complex nonlinear relationships within multidimensional datasets (Yin et al., 2023; Lyu et al., 2023). Nie and Sun (2024) proposed a method combining deep learning and dynamical-statistical projection model to correct ECMWF sub-seasonal precipitation forecasts, demonstrating significant skill enhancements in Southwest China for lead times of up to 30 d. Lyu et al. (2024) utilized a hybrid CSG-UNET model to refine ECMWF ensemble precipitation forecasts over the China mainland during summer, extending the effective forecast skill to approximately 4 weeks. In general, current research predominantly targets ECMWF forecasts due to their superior predictive skill (Falck et al., 2025; Dong et al., 2025; Andrade et al., 2024; Quedi and Fan, 2020). However, most studies are confined to the first 30 d, as results beyond this point are either omitted or unreliable for practical use (Nie and Sun, 2024; Lyu et al., 2024; Yin et al., 2023). In contrast, alternative models such as that from the UK Met Office (UKMO), which provide 60 d forecasts, remain under-explored due to the lower raw accuracy of their precipitation forecasts. To address this gap, this study aims to apply bias correction to 60 d precipitation forecasts, thereby enhancing reliability and extending the effective forecast horizon.

In addition to errors in precipitation forecasts, inherent uncertainties stemming from the hydrological model structural generalization and parameter estimation can also propagate into streamflow predictions (Donegan et al., 2021; Dion et al., 2021), rendering the correction of hydrological simulation errors imperative. Common correction strategies range from autoregressive methods and data assimilation to deep learning approaches (Tanguy et al., 2025; Sabzipour et al., 2023; Siqueira et al., 2021). Among these, the autoregressive with exogenous input (ARX) model is particularly suitable for addressing the systematic bias and serial correlation inherent in hydrological residuals, providing a computationally efficient and effective solution for post-processing (McInerney et al., 2021; Sharma et al., 2019).

Therefore, this study proposes a 60 d streamflow forecasting framework that integrates a convolutional neural network (CNN) to correct precipitation forecasts from the UKMO numerical weather prediction model, the geomorphology-based eco-hydrological model (GBEHM) to simulate hydrological processes, and an ARX model to minimize residual errors. Applying the proposed framework in the Upper Yangtze River Basin (UYRB), we expect this research to offer a practical and reliable tool for decision support with extended lead times for hydrological hazard early warning.

2 Study area and data

2.1 Study area

The Upper Yangtze River Basin (UYRB) is highly vulnerable to extreme hydrological events driven by uneven streamflow distribution, frequently triggering devastating floods and droughts that threaten downstream infrastructure and regional water security (Liang et al., 2023; Wang et al., 2022). In addition, the basin serves as a strategic hub for China's hydropower generation, playing a vital role in regional energy security (Zhong et al., 2020). Consequently, improving forecast accuracy at medium- and long-range horizons in this region is of paramount importance for integrated disaster mitigation and reservoir operation (Su et al., 2017).

This study focuses on the UYRB above the Shigu hydrological station, which spans 90–101° E and 26–36° N with a drainage area of approximately 2.2×10⁵ km² (Fig. 1). The terrain is complex with elevation gradually decreasing from over 6500 m in the headwaters to about 1800 m at the outlet. It is a typical alpine cold mountainous zone dominated by grassland, shrubland, and barren areas. The region has a mean annual precipitation of around 670 mm and a mean annual streamflow of around 1300 m³ s⁻¹ observed at the Shigu station.

Figure 1 The upper Yangtze River Basin controlled by the Shigu hydrological station.

2.2 Data sources

2.2.1 Observed meteorological and hydrological data

To support model training and validation, the China Gauge-Based Daily Precipitation Analysis (CGDPA), provided by the National Meteorological Information Center, is used as the reference observational dataset. This 0.25° daily precipitation gridded product is derived from over 2400 gauge stations across China and covers the period from 1960 to 2019. Its reliability in capturing precipitation characteristics over complex terrain has been extensively validated, and detailed station distribution density maps along with gauge information can be found in Shen and Xiong (2016).

To drive the hydrological model, additional meteorological inputs, including air temperature, wind speed, relative humidity, and sunshine duration, are obtained from meteorological stations operated by the China Meteorological Administration (CMA, http://data.cma.cn, last access: 21 May 2026). To align with the hydrological model configuration, these gauge observations are spatially interpolated onto the 8 km ×8 km simulation grid using the elevation-adjusted angular distance weighting technique (Yang et al., 2004). Furthermore, to meet the temporal requirements, the daily data are disaggregated into hourly forcing series following the algorithms proposed by Gao et al. (2015).

To calibrate the hydrological model and evaluate forecast performance, daily streamflow observations from the Shigu hydrological station for 1960 to 2019 are collected from the hydrological yearbook.

2.2.2 Underlying surface data

The distributed model GBEHM requires topographic data and underlying surface properties. Topographically, a 90 m resolution digital elevation model (DEM) from the Shuttle Radar Topography Mission (SRTM) database (Jarvis et al., 2008) is utilized to delineate the river network and basin boundaries. Regarding soil properties, soil texture distributions are derived from Shangguan et al. (2014), while key soil hydraulic parameters, including saturated hydraulic conductivity, saturated and residual water contents, and van Genuchten parameters, are sourced from the China Soil Hydraulic Parameters Dataset (Dai et al., 2013). For land surface characterization, this study uses a 100 m resolution land use map from the Resource and Environment Data Cloud Platform (http://www.resdc.cn/, last access: 21 May 2026). Finally, vegetation dynamics, specifically leaf area index (LAI) and fraction of photosynthetically active radiation (FPAR), are parameterized using the GIMMS NDVI3g-based datasets developed by Zhu et al. (2013).

2.2.3 Meteorological forecast data

Meteorological forecast data are sourced from the Sub-seasonal to Seasonal Prediction Project (Vitart et al., 2017) which provides medium- and long-range forecasts generated by numerical weather prediction models from various international operational centers. Among the models providing a 60 d lead time (including those from CMA and UKMO), the UKMO model is selected as the source of raw precipitation forecasts due to its superior performance in the study region (Li et al., 2019). The dataset, provided at a 0.25° spatial resolution, covers the hindcast period of 1997–2016 with 44 initialization dates per year. The forecast variables include the target precipitation forecasts together with a suite of predictors for correction: convective precipitation, 2 m dew point temperature, 2 m air temperature, total cloud cover, and multi-level variables (geopotential height, specific humidity, temperature, and wind speed at 200, 500, and 850 hPa). The selection of these specific predictors is based on previous studies which have demonstrated their efficacy in improving forecast accuracy (Zhang et al., 2024, 2023; Li et al., 2023).

Figure 2 Flowchart of the proposed forecasting framework.

3 Methods

3.1 Overview

In this study, we propose a framework to improve medium- and long-term forecasts of daily precipitation and streamflow with lead times of up to 60 d to enhance hydrological hazard early warning capabilities. The flowchart of the proposed framework is illustrated in Fig. 2. First, a LeNet-based CNN model is developed to perform bias correction on the precipitation forecasts by modeling the relationship between local grid precipitation and predictor variables from both the target and surrounding grids. Subsequently, the gridded precipitation forecast serves as input to drive the distributed hydrologic model GBEHM, generating daily streamflow forecasts for the 1–60 d horizon. Finally, the simulated streamflow undergoes post-processing via an ARX model to yield the final forecasts. The performance of these forecasts is quantitatively evaluated using metrics including the root mean squared error (RMSE), temporal correlation coefficient (TCC), relative error (RE) and the Nash–Sutcliffe efficiency (NSE).

3.2 Convolutional neural networks

A modified LeNet-based convolutional neural network (CNN) is designed to model the complex non-linear relationship between local grid precipitation and large-scale atmospheric conditions. The model input consists of a multivariate 3D tensor with dimensions of $\mathrm{20}\times \mathrm{9}\times \mathrm{9}$ representing a 9×9 spatial neighborhood centered on the target 0.25° grid cell. The input comprises 20 predictor channels, including the DEM, precipitation forecasts and other meteorological variables detailed in Sect. 2.2.3. Structurally, the network integrates three primary components: embedding layers, a convolutional backbone, and fully connected layers. To explicitly capture spatial location dependencies, the latitude and longitude indices of each grid are encoded via two independent embedding layers. The convolutional backbone consists of four cascaded layers (configured with 64, 32, 16, and 8 filters, respectively), each followed by batch normalization to enhance training stability. The extracted spatial feature maps are flattened and concatenated with the spatial embeddings before passing through two fully connected hidden layers.

The final layer of the network produces three real numbers (z₁, z₂, z₃) that govern the parameters of a censored shifted gamma (CSG) distribution. This probabilistic formulation is specifically adopted to address the mixed discrete-continuous nature of precipitation, effectively characterizing zero-inflation (dry days), skewness, and the transition from light to heavy rainfall (Towler et al., 2025; Zhang et al., 2017). To ensure mathematical validity, the raw network outputs are transformed into the physical distribution parameters of shift (γ), mean (μ), and standard deviation (σ) via specific activation functions: $\mathit{\gamma }=-\sqrt{z_{\mathrm{1}}^{\mathrm{2}}+\mathit{ϵ}}$, μ=exp(z₂), and $\mathit{\sigma }=\sqrt{\exp \left(z_{\mathrm{3}}\right)+\mathit{ϵ}}$, where ϵ is a small constant for numerical stability. Mathematically, the CSG distribution is defined by a latent variable $X=Z+\mathit{\gamma }$, where $Z\sim \mathrm{Gamma}(\mathit{\alpha },\mathit{\beta })$ with shape $\mathit{\alpha }=(\mathit{\mu }/\mathit{\sigma }{)}^{\mathrm{2}}$ and scale $\mathit{\beta }={\mathit{\sigma }}^{\mathrm{2}}/\mathit{\mu }$. The precipitation Y is then obtained by censoring X at zero, such that $Y=max(\mathrm{0},X)$. The probability of zero precipitation is given by $P(Y=\mathrm{0})=F_{\mathrm{\Gamma }}(-\mathit{\gamma };\mathit{\alpha },\mathit{\beta })$, where F_Γ is the cumulative distribution function of the Gamma distribution. Finally, the model generates a deterministic forecast by constructing a large-scale pseudo-ensemble from the predicted CSG distribution at equal quantiles and calculating the ensemble mean (see Sect. S1 in the Supplement for the detailed sampling procedure). This approach enables a comprehensive quantification of uncertainty and ensures that low-probability, high-impact extreme events in the tail of distribution are adequately captured, thereby mitigating the smoothing effect inherent in deterministic forecasts.

The model is trained using the continuous ranked probability score (CRPS) as the loss function to optimize probabilistic performance:

$$\begin{array}{}\text{(1)}& \mathrm{CRPS}\left(F,y\right)={\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}{\left(F\left(x\right)-H_{\left\{x\mathit{>=}y\right\}}\right)}^{\mathrm{2}}\mathrm{d}x\end{array}$$

where F is the predicted cumulative distribution function and H is the Heaviside step function that equals 1 if x≥y and 0 otherwise. To ensure a robust and unbiased evaluation, the 20-year UKMO hindcast dataset (1997–2016) is partitioned into training (1997–2008, 12 years), validation (2013–2016, 4 years), and testing (2009–2012, 4 years) subsets. This specific testing period is strategically aligned with the final streamflow evaluation to ensure the CNN model is assessed on entirely unseen data. The optimization is performed using the Adam algorithm, with hyperparameters tuned to minimize overfitting.

3.3 Geomorphology-based eco-hydrological model

The GBEHM is a physically-based, distributed model designed to simulate hydrological processes in topographically complex catchments. To address spatial heterogeneity, the model employs a hierarchical sub-grid parameterization scheme (Yang et al., 2015), where the hillslope-valley system serves as the fundamental computational unit. Specifically, the study region is first delineated into sub-basins based on the DEM to construct the river network topology for flow routing. Within each sub-basin, the landscape is first discretized into the grid system to integrate meteorological forcing data and capture the spatial heterogeneity of land surface properties. Finally, these grids are further subdivided into hillslope-valley units. This allows the model to solves vertical water and energy balances at the hillslope scale and dynamically aggregates the runoff to the sub-basin scale for lateral river routing. Runoff generation is explicitly resolved through three primary pathways: surface overland flow, lateral subsurface flow, and groundwater discharge. Specifically, vertical soil water movement in the unsaturated zone is governed by the one-dimensional Richards' equation, while lateral flow in the saturated zone and groundwater-river exchange are quantified using Darcy's law and mass balance principles (Cong et al., 2009). To close the water and energy budget, the hydrological module is coupled with the Simple Biosphere Model 2 (SiB2), which estimates evapotranspiration losses, including canopy interception and soil evaporation, based on energy transfer within the soil–plant–atmosphere continuum (Sellers et al., 1996). To accurately simulate these complex hydrological processes, several key model parameters require careful calibration. These primarily include the evaporation parameters (C₁, C₂, C₃), soil saturated hydraulic conductivity (K_s), groundwater transmissivity (K_g) and storage coefficient, the snowmelt factor (M_f), and the hillslope shape parameter (f_ss).

A distinctive feature of GBEHM is its enhanced representation of cryosphere hydrology, making it particularly robust for cold regions. The model integrates a rigorous coupled heat and water balance equation (Flerchinger and Saxton, 1989) to simulate soil freeze-thaw cycles, which critically alter soil hydraulic conductivity and infiltration capacity. The vertical soil profile is discretized into a multi-layer structure extending to a depth of 50 m. Crucially, the active soil layer (top 1–3 m) features a refined mesh resolution to accurately capture the dynamics of the active layer thickness and the maximum frozen depth (Guo and Wang, 2013).

The selection of GBEHM for this study is driven by its proven efficacy in simulating complex hydrological processes within the Tibetan Plateau and other high-altitude regions (Shi et al., 2020; Gao et al., 2018). The study area is characterized by rugged terrain and a cold climate, where glacier and snow melt as well as soil freeze-thaw processes exert a dominant control on the hydrological regime. Since GBEHM explicitly accounts for phase changes in the soil and the delayed runoff response caused by the cryosphere, it is well-suited for accurate runoff forecasting in this catchment. In this application, the UYRB is discretized into an 8 km ×8 km grid system and further delineated into 479 sub-basins based on the DEM.

3.4 Autoregressive with exogenous input

While the GBEHM captures the fundamental physical processes of runoff generation, hydrological simulations inevitably contain systematic biases and persistent errors due to uncertainties in model structure and parameters. To mitigate these discrepancies and improve forecast accuracy, this study implements a statistical post-processing technique. Instead of modeling the streamflow directly, an autoregressive with exogenous input (ARX) model is constructed to simulate and correct the hydrological residuals.

Let Q_obs(t) and Q_sim(t) denote the observed and GBEHM-simulated streamflow at time t, respectively. To stabilize the variance and eliminate seasonal scale effects, both series are first standardized using the mean (μ) and standard deviation (σ) derived strictly from the calibration period. The standardized hydrological error, E_std(t), is defined as the deviation of the simulated flow from the observations:

$$\begin{array}{}\text{(2)}& E_{\mathrm{std}}\left(t\right)=Q_{{\mathrm{obs}}_{\mathrm{std}}}\left(t\right)-Q_{{\mathrm{sim}}_{\mathrm{std}}}\left(t\right)\end{array}$$

Hydrological errors typically exhibit strong temporal autocorrelation and dependence on flow magnitude. To capture these dynamics, the error at the current time step is modeled as a function of the current simulated state, antecedent simulated states, and antecedent errors. The ARX model formulation is expressed as:

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{E}_{\mathrm{std}}^{\ast }& \left(t\right)={\mathit{\beta }}_{\mathrm{0}}+{\sum }_{i=\mathrm{1}}^p{\mathit{\varphi }}_i{E}_{\mathrm{std}}\left(t-i\right)\\ & +{\sum }_{j=\mathrm{0}}^k{\mathit{\gamma }}_j{Q}_{{\mathrm{sim}}_{\mathrm{std}}}\left(t-j\right)+\mathit{\epsilon }\left(t\right)\end{array}\end{array}$$

where p and k are the lag orders for the autoregressive error term and the exogenous simulated streamflow term respectively, which are optimally determined by identifying the combination that yields the minimum values of the Akaike Information Criterion (AIC) and the Bayesian Information Criterion (BIC) (Hipel and McLeod, 1994); ϕ_i represents the autoregressive coefficients describing the persistence of model errors; γ_j are coefficients for the exogenous input, accounting for magnitude-dependent bias; β₀ is the intercept standing for the independent bias, and ε(t) is the residual white noise. The model parameters are calibrated using the least squares method. The corrected error $E_{\mathrm{std}}^{\ast }\left(t\right)$ is generated recursively and added to the standardized forecast simulation. Finally, the post-processed streamflow is obtained by applying inverse standardization to the result.

3.5 Evaluation metrics

The accuracy of the corrected precipitation forecasts is evaluated using the RMSE to quantify error magnitude and the TCC to measure phase consistency. These metrics are defined as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }\mathrm{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{n}}{\sum }_{i=\mathrm{1}}^n(S_i-O_i{)}^{\mathrm{2}}}\text{(5)}& {\displaystyle }\mathrm{TCC}={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(S_i-\stackrel{\mathrm{‾}}{S})(O_i-\stackrel{\mathrm{‾}}{O})}{\sqrt{{\sum }_{i=\mathrm{1}}^n(S_i-\stackrel{\mathrm{‾}}{S}{)}^{\mathrm{2}}}\sqrt{{\sum }_{i=\mathrm{1}}^n(O_i-\stackrel{\mathrm{‾}}{O}{)}^{\mathrm{2}}}}}\end{array}$$

where S_i and O_i represent the simulated (or forecasted) and observed values at time step i, respectively; $\stackrel{\mathrm{‾}}{S}$ and $\stackrel{\mathrm{‾}}{O}$ denote their corresponding means; and n is the total number of samples.

For streamflow simulations, model performance is appraised using three standard hydrological efficiency criteria: the RMSE (same as above), RE, and NSE. The RE and NSE are calculated as:

$$\begin{array}{}\text{(6)}& {\displaystyle }\mathrm{RE}={\displaystyle \frac{\stackrel{\mathrm{‾}}{S}-\stackrel{\mathrm{‾}}{O}}{\stackrel{\mathrm{‾}}{O}}}\times \mathrm{100}\mathit{\%}\text{(7)}& {\displaystyle }\mathrm{NSE}=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n(S_i-O_i{)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^n(O_i-\stackrel{\mathrm{‾}}{O}{)}^{\mathrm{2}}}}\end{array}$$

4 Results

4.1 Calibration and validation of GBEHM and ARX model

The GBEHM is calibrated for the period of 1990–2005 and validated for the period of 2006–2015, leveraging extended historical gauge records to robustly capture complex long-term hydrological dynamics. Figure 3 illustrates the comparison of the simulated daily streamflow and the observations at the Shigu station. The RMSE, RE and NSE values for the calibration period are 330 m³ s⁻¹, 20.2 % and 0.93, respectively, and for the validation period are 347 m³ s⁻¹, 20.2 % and 0.90, which indicating the high reliability of the GBEHM for daily streamflow simulation in the study area.

The ARX model is calibrated during 2000–2008 and validated during 2009–2012, striking a balance between avoiding early-stage GBEHM warm-up uncertainties and maintaining sufficient data length for robust parameterization. After testing different model orders, p=3 and k=3 is selected based on validation results (see Sect. S2 for details). For the validation period, the raw GBEHM simulation shows the RMSE of 351 m³ s⁻¹, RE of 18.9 %, and NSE of 0.91. Post-correction with the ARX model, these metrics improve to 267 m³ s⁻¹, 12.5 %, and 0.95, respectively. This confirms that the ARX model effectively mitigates errors in streamflow simulation.

Figure 3 Comparison of observed and simulated daily streamflow at the Shigu station during the calibration (1990–2005) and validation (2006–2015) periods.

Figure 4 (a) RMSE and (b) TCC of the UK and UK-CNN precipitation forecasts, calculated for aggregated 10 d lead times.

4.2 Evaluation of precipitation forecasts

Figure 4 illustrates the performance of the areal-averaged raw and corrected (denoted as UK and UK-CNN) precipitation forecasts in terms of RMSE and TCC across different lead time ranges for the 2009–2012 period.

In terms of overall performance averaged over all lead times, the areal-averaged UK raw forecasts exhibit an RMSE of 2.6 mm d⁻¹ and a TCC of 0.62. The UK-CNN effectively improve these metrics, reducing the average RMSE to 1.7 mm d⁻¹ (a decrease of 35 %) and increasing the TCC to 0.74 (an improvement of 18 %). Regarding the temporal evolution, the skill of UK raw forecasts naturally deteriorates with increasing lead time. The UK-CNN follow a similar declining trend initially but stabilize at a relatively constant level after approximately 20 d. Consequently, the relative improvement gap widens as the lead time extends. For instance, the RMSE reduction increases from 33 % in the first 20 d to 35 % and 38 % in the middle and last 20 d periods, respectively. Similarly, the TCC improvement rises significantly from 5 % to 17 % and 30 % over the same intervals.

Figure 5 (a) RMSE, (b) RE and (c) NSE of streamflow forecasts driven by raw (UK) and corrected (UK-CNN) precipitation, and the streamflow subsequently post-processed by the ARX model (UK-ARX), calculated for aggregated 10 d lead times.

4.3 Evaluation of streamflow forecasts

Figure 5 shows the RMSE, RE, and NSE for streamflow forecasts driven by raw UK and corrected UK-CNN precipitation, as well as the UK-CNN driven forecasts further post-processed by the ARX model (denoted as UK-ARX). In terms of overall performance averaged over all lead times, the streamflow forecasts driven by raw UK precipitation yield an RMSE of 930 m³ s⁻¹, an RE of 48.2 %, and an NSE of 0.33. Using UK-CNN precipitation significantly improves these metrics to 706 m³ s⁻¹, 21.0 %, and 0.61, representing relative improvements of 24 %, 56 %, and 87 %, respectively. The UK-ARX achieves the best performance, further reducing the RMSE and RE to 596 m³ s⁻¹ and 17.4 %, and increasing the NSE to 0.72. This corresponds to additional improvements of 16 %, 17 %, and 18 % over the UK-CNN benchmark.

As the forecast lead time extends, the skill of streamflow forecasts naturally declines with increasing lead time; however, the correction methods notably mitigate the rate of this deterioration. For instance, the RE of the UK-driven forecasts rises sharply from 25.2 % to 63.2 % by lead day 51–60, while the NSE plummets from 0.78 to 0.09. Conversely, the RE of the UK-CNN forecasts increases only marginally from 17.6 % to 24.5 %, with the NSE showing a more gradual decline from 0.87 to a moderate level of 0.43. The UK-ARX further alleviates this degradation trend, maintaining the RE below 20.9 % and sustaining the NSE above 0.59 throughout the entire 60 d period.

Moreover, the extent of improvement exhibits distinct trends. The performance gap between UK and UK-CNN becomes more pronounced at longer lead times. For instance, the reduction in RE achieved by UK-CNN (compared to UK) increases from around 44 % in the first 20 days to 58 % and 61 % in the middle and last 20 d periods, respectively. Conversely, the incremental improvement from UK-ARX (relative to UK-CNN) remains stable, providing RE reductions of around 21 %, 16 %, and 16 % across these lead time ranges.

5 Discussion

5.1 Efficacy of the CNN-based precipitation bias correction

This study extends the horizon of reliable precipitation forecasting to 60 d by employing a CNN-based deep learning model for bias correction. Specifically, the model achieves a TCC of approximately 0.70 in the extended range (days 21–30), representing a substantial advantage over the results reported by Lyu et al. (2023) who achieved a TCC of about 0.35 for the same period when correcting summer precipitation over Southern China using ECMWF data with lead times up to 30 d. More importantly, the model maintains a TCC above 0.66 even at long lead times of 51–60 d, with the RMSE stabilizing after the initial 20 d. These findings validate the model's robustness in mitigating error accumulation and sustaining forecast reliability up to 60 d.

This substantial enhancement in performance is primarily attributed to the deep learning approach and the probabilistic output scheme. First, unlike traditional point-to-point statistical correction methods (e.g., quantile mapping) that process each grid cell independently, the CNN architecture effectively extracts spatial dependencies from surrounding atmospheric conditions using a 9×9 grid neighborhood, and models the complex non-linear interactions between the multiple meteorological predictors and local precipitation patterns (Baño-Medina et al., 2020; Pan et al., 2019). Second, instead of generating deterministic point predictions, this model explicitly corrects the CSG distribution optimized via CRPS. This approach not only precisely captures the mixed zero-inflated and heavy-tailed characteristics of precipitation (Ghazvinian et al., 2022) but also mitigates the smoothing effect common in deep learning by deriving forecasts from a pseudo-ensemble (Ravuri et al., 2021; Rasp and Lerch, 2018), thereby preserving signal variability and ensuring high reliability even at extended lead times.

5.2 Added value of ARX-based hydrological post-processing

Our framework extends the valid streamflow forecast horizon to 60 d with the ARX post-processing model providing significant added value in mitigating hydrological residuals. Specifically, while the skill of UK-CNN forecasts naturally degrades at extended horizons (with NSE dropping to 0.43 by day 60), the ARX model effectively alleviates this decline, sustaining a reliable NSE above 0.57 and keeping RE below 22 % throughout the entire period. Furthermore, the ARX model achieves a consistent performance boost (about 16 %–21 % reduction in RE) across all lead times, demonstrating its unique capability to mitigate intrinsic structural errors within the hydrological model and ensuring the validity of forecasts up to 60 d.

The efficacy of the ARX model in achieving these results stems from its capacity to mitigate intrinsic hydrological modeling uncertainties while remaining a computationally efficient solution. While precipitation bias correction improves inputs, hydrological models inevitably introduce intrinsic biases due to simplified parameterizations, structural deficiencies, or uncertain initial conditions. The ARX model addresses this by exploiting error autocorrelation to effectively remove systematic and temporal discrepancies, ensuring stable correction capabilities across the entire forecast horizon. Furthermore, compared to data assimilation approaches which aim to reduce these uncertainties by continuously updating internal model states via computationally intensive ensemble simulations (Nearing et al., 2022; Liu et al., 2012), the ARX model provides a more efficient alternative by bypassing these complex internal adjustments entirely, thereby serving as a straightforward post-processing solution characterized by interpretability and minimal data requirements.

5.3 Attribution of the streamflow forecast error

An error decomposition framework is adopted to evaluate the relative contributions of error sources. The total streamflow forecast error is decomposed into hydrological model error (MSE_m), defined as the MSE between observed precipitation driven simulations and observed streamflow, and precipitation forecast error (MSE_p), defined as the MSE between observed and forecast precipitation driven simulations. It should be noted that MSE_m encompasses not only intrinsic hydrological model uncertainties but also the errors introduced by interpolating 0.25° precipitation data onto the 8 km grid. Additionally, the error reduction contribution from the ARX post-processing (MSE_ARX) is quantified as the MSE between the streamflow forecasts before and after ARX correction. Notably, the total forecast error is generally smaller than the arithmetic sum of MSE_m and MSE_p. This non-additivity is attributed to the interaction between error sources, where a compensation effect between precipitation biases and hydrological model deficiencies helps mitigate the overall error.

Figure 6 Decomposition of streamflow forecast errors into hydrological model (MSE_m) and precipitation forecast (MSE_p) components before and after ARX post-processing across lead times.

Figure 6 depicts the absolute values and relative proportions of MSE_m and MSE_p before and after ARX post-processing. Generally, the magnitudes of MSE_m and the ARX-induced error reduction MSE_ARX show little variation with lead time, while MSE_p exhibits a steady increase. Prior to ARX post-processing, MSE_m is the dominant error source within the first week of lead time, accounting for over 50 % of the total. This proportion drops rapidly and stabilizes at around 0.15 after 30 d, while the share of MSE_p inversely rises to stabilize at around 0.85. With the application of ARX, the dominance of the hydrological model error is significantly curtailed; it ceases to be the primary error source after a lead time of only 5 d and stabilizes at a proportion of roughly 0.10 after 20 d. Notably, the MSE displays a pronounced 8 d periodicity with peaks occurring at the third phase of each cycle. This pattern stems from the systematic alignment of validity dates caused by the 8 d forecast interval, where errors from sparse extreme events are systematically projected onto specific lead times, thus driving the peak formation (details provided in Sect. S3).

The error decomposition analysis reveals a clear shift in dominance: short-term skill is constrained by hydrological modeling, while medium-to-long-term skill is limited by precipitation forecasts. In our study, hydrological model error (MSE_m) stabilizes at a notably low contribution of approximately 15 % at long lead times, outperforming the about 30 % reported by Dong et al. (2025) who used a hybrid deep learning-conceptual model. This superior performance is attributed to the GBEHM's capacity to resolve the spatial heterogeneity of the underlying surface and explicitly model the cryosphere hydrological processes, thereby representing the complex runoff generation mechanisms in this cold mountainous region more accurately. The ARX post-processing further compresses this error to about 10 % and shortens its dominance period to just 5 d, confirming the high reliability of process-based models augmented by post-processing. Nevertheless, the upper limit of forecast skill remains constrained by meteorological forcing, which dominates the long-term error contribution (>85 %). This highlights the indispensability of the proposed framework: effective medium-to-long-range forecasting demands both a high-precision hydrological model to minimize internal uncertainty and advanced precipitation bias correction to mitigate the overwhelming external forcing errors.

5.4 Limitations of this study

One limitation of this study lies in the use of an identical set of predictors across the entire 1–60 d forecast horizon. Although separate models are trained for different lead times, allowing the weights of these predictors to adjust dynamically, previous studies indicate that the dominant atmospheric drivers often shift as the lead time extends (Lyu et al., 2023). Therefore, incorporating predictor importance ranking to select lead-time-specific predictor subsets could likely further enhance the model's fitting capability.

Beyond this specific constraint, future research could expand in several directions. First, it is essential to evaluate the potential of multi-model ensembles (e.g., incorporating the CMA 60 d forecast) and verify the generalizability of our framework across more diverse catchments. Second, although this study uses the mathematical expectation of probabilistic forecasts to balance computational efficiency and extreme event retention, future research should transition to fully probabilistic streamflow forecasting, such as employing surrogate models to ultimately maximize the informational value of early warnings. Finally, while the current framework is optimized for overall streamflow simulation accuracy, it holds significant potential for extreme hydrological hazard forecasting. To fully realize its utility for flood early warning, future efforts should focus on specializing the CNN and ARX models for extreme-value optimization (e.g., by integrating extreme-value-guided loss functions in the CNN and employing threshold-based autoregression in the ARX), and adjusting GBEHM calibration objectives to prioritize high-flow precision.

6 Summary and conclusions

This study constructs a robust 60 d streamflow forecasting framework by coupling a CNN for correcting UKMO precipitation forecasts, the GBEHM for process-based hydrological simulation, and an ARX model for mitigating hydrological modelling residuals. This framework is applied to the Upper Yangtze River Basin (UYRB).

The results demonstrate a significant extension of the valid streamflow forecast horizon to 60 d. Compared to raw forecasts, the proposed method reduced the RE from a range of 25.2 %–63.2 % to 13.7 %–20.9 %, while elevating the NSE from 0.78–0.09 to a reliable 0.91–0.59. In terms of contributions, this performance boost is driven by two key components. First, the CNN-based model significantly improved meteorological inputs by reducing precipitation RMSE by 35 % and elevating TCC from 0.62 to 0.74, particularly at longer lead times; this enhancement in precipitation accounts for approximately two-thirds of the total improvement in streamflow forecasts. Meanwhile, the ARX post-processing contributed the remaining 33 % to the total streamflow error reduction by effectively mitigating intrinsic hydrological residuals.

Our framework highlights the distinct advantages of integrating a physically robust hydrological model with a dual-stage error correction strategy. First, the distributed GBEHM demonstrates exceptional precision in characterizing complex catchment dynamics, maintaining notably low residual errors (15 %) even at extended lead times. Second, by strategically coupling deep learning (CNN) for precipitation input correction with statistical post-processing (ARX) for hydrological output refinement, the framework systematically mitigates both external meteorological biases and internal simulation uncertainties. This synergy yields forecasts that are both volumetrically accurate and temporally consistent. From an operational perspective, the pre-trained nature of the framework ensures high computational efficiency for real-time deployment, and the requisite real-time UKMO forecasts can be readily secured through professional institutional data acquisition. Therefore, the proposed framework provides a highly reliable solution with an extended 60 d horizon for hydrological hazard early warning and proactive flood and drought risk mitigation.