A two-stage hierarchical runoff bias correction framework for observation-constrained river flow reconstruction

Mapping historical river flows across all stream orders is of vital importance for water resources management, flood risk assessment, energy production, and ecosystem preservation^1,2. Beyond historical analysis, these high-resolution estimates can be used to predict river flows with appropriate lead times, providing critical information hydrological extreme forecasting, such as flood awareness³ and drought monitoring⁴. However, most river reaches remain ungauged due to the sparse distribution of in situ gauging stations⁵. While remote sensing offers a promising data source, it still lacks the long-term temporal consistency required for historical flow reconstruction⁶. Therefore, the most common approach to estimate discharge at the reach scale is to combine a runoff generation model with a river routing model.

Existing modeling approaches for flow reconstruction generally fall into process-based and data-driven categories. The first pathway utilizes process-based models to simulate water fluxes, together with in situ, reanalysis, and satellite observations^7,8. A prominent example is the GloFAS-ERA5 reanalysis, a key component of the European Commission’s Copernicus Service⁹. By coupling the HTESSEL land surface model with the LISFLOOD routing model, GloFAS-ERA5 generates long-term, global river discharge at a 0.05° daily resolution, providing a critical baseline for flood awareness¹⁰. Other efforts include the Global Reach-Level A Prior Discharge Estimates (GRADES) dataset¹¹, which paired the VIC hydrological model with RAPID routing to estimate daily flow for nearly 3 million river reaches for 1979–2014, and the subsequent dataset Global Reach-Level Flood Reanalysis (GRFR) GRDR which pushed these estimates to a higher temporal frequency at 3-h¹². The second pathway relies on data-driven techniques, predominantly machine learning (e.g., Long Short-Term Memory. LSTM) to learn intricate runoff generation patterns directly from massive datasets without explicit hydrological processes^13,14. Previous researchers have often trained and applied LSTM to simulate river flows at the basin outlet by treating each catchment as a single hydrological response unit and utilizing basin-averaged meteorological forcing as input^15,16,17. To consider the internal spatial variability within the basin, recent studies have transitioned toward grid-based LSTM modeling. For instance, Yang et al.¹⁸ proposed a Grid LSTM-RAPID scheme that trains LSTM on small basins but applies it to rectangular grids, allowing for distributed runoff estimation and network-based flow routing over vast river systems.

Despite these advancements, simulated river flows often exhibit biases^19,20. To further mitigate these errors, post-processing or bias correction (BC) methods are frequently employed. Gauge-specific post-processing can improve agreement at individual stations, but it does not automatically propagate the correction to intervening reaches or define how corrections at nested gauges should be reconciled within one connected routing calculation (e.g.,¹⁹). Therefore, a more robust strategy is to correct the runoff before it enters the routing process. A notable example is the VIC-BC dataset¹¹, which applies grid-by-grid correction based on ML-derived global runoff characteristic maps²¹. However, these global BC products have several key drawbacks for regional applications. The coarse spatial resolution of 0.25 degree often fails to capture hydrological processes at a scale that meets locally societal needs. These products also suffer from low update frequencies and are not available for the most recent years. More importantly, these methods typically treat grid cells as independent units and overlook the hierarchical topology of the river network. Consequently, they cannot easily distinguish whether a downstream bias originates from the headwaters or the local drainage area. This lack of spatial connectivity and structural consistency ultimately limits the overall simulation accuracy across the full stream orders.

To provide an improved estimate of river flows across all stream orders, we propose the Two-Stage Hierarchical Bias Correction framework. This study utilizes ERA5 runoff (0.1° resolution, daily updates;⁹) as the runoff forcing and a high-resolution vector-based network MERIT Hydro at a 5 km scale²². By incorporating in situ observations from 18 stations in the case study Yangtze River, the proposed framework identifies the water balance discrepancy between simulated and observed discharge and linearly attributes these errors back to the grid-based runoff. Specifically, we first establish a runoff baseline in headwater basins and then correct the runoff incrementally between successive downstream stations. These corrected runoff fields are then routed through the mizuRoute routing model²³ to reconstruct high-accuracy daily discharge from 1980 to 2024 for approximately 40,000 river reaches. This approach respects the hierarchical topology of the river system and is scalable to any global basin for constructing locally relevant, high-resolution flow datasets.

Study area

The Yangtze River Basin (YRB, Fig. 1), encompassing a drainage area of approximately 1.8 million km², serves as a representative and complex environment for testing the proposed framework. As the largest river in China, it originates from the Tibetan Plateau at elevations exceeding 5000 m and traverses diverse topographical zones before discharging into the East China Sea. This 6300 km streamline creates a highly heterogeneous hydrological system characterized by a distinct subtropical monsoon climate. Annual precipitation varies significantly across the basin, ranging from less than 400 mm in the semi-arid headwaters of the Jinsha River to over 1600 mm in the humid downstream reaches. Managing water resources here is critical due to the presence of massive infrastructure and the frequent threat of catastrophic flooding. To capture this complexity, our modelling framework accounts for the entire basin’s heterogeneity, ensuring that both high-altitude headwaters and lowland mainstems are represented with physical rigor.

Data sources

ERA5 runoff forcing

The primary runoff forcing used in this research is the ERA5 reanalysis dataset from the European Centre for Medium-Range Weather Forecasts (ECMWF). We specifically extract total runoff, which is the sum of surface and subsurface runoff, at a spatial resolution of 0.1° (~ 11 km). ERA5 represents a significant advancement over previous reanalyses by incorporating the sophisticated HTESSEL land-surface scheme and improved data assimilation techniques. A critical advantage of ERA5 for our study is its long-term temporal consistency and near-real-time availability. Unlike legacy products that stop in 2014 or 2019, ERA5 provides daily updates to the present year. This allows our reconstruction to span from 1980 to 2024, effectively capturing recent extreme events in the Yangtze Basin. By utilizing the raw 0.1° runoff as our baseline, we provide a globally standardized forcing that ensures the scalability of the proposed framework to other transboundary basins.

VIC-BC runoff (benchmark dataset)

To evaluate the relative performance of our proposed framework, we incorporate the VIC-BC (Variable Infiltration Capacity—Bias Corrected) runoff dataset as a primary benchmark¹¹. VIC-BC was generated by applying a grid-based correction to the output of the VIC hydrological model using a machine learning approach. Specifically, it utilized Random Forest algorithms to develop global runoff characteristic maps (Qc maps) based on a large sample of gauged catchments²¹.

While VIC-BC has been widely used for global flow reconstruction, it originally operates at a coarser spatial resolution of 0.25° but downscaled to 0.05° and its temporal coverage is typically limited to the 1979–2019 period. Furthermore, as an independent grid-based correction, VIC-BC does not account for the hierarchical routing constraints between upstream and downstream reaches. In this study, we route the VIC-BC runoff through the same river routing setup to quantify how our framework improves upon this established global standard, particularly in capturing the spatial heterogeneities of a large-scale basin like the Yangtze.

MERIT hydro vector-based river network

The representation of the river channel system is based on MERIT Hydro²², which is currently the most accurate global hydrography dataset. MERIT Hydro significantly improves upon the SRTM (Shuttle Radar Topography Mission) by correcting for vegetation-induced height errors and speckle noise, which often lead to misidentified flow paths in densely forested or low-relief areas. From this product, we derived a vector-based river network for the YRB consisting of approximately 40,000 reaches. Each reach is defined with an average length of 5 km, providing a much finer representation than traditional gridded models. This vector-based approach treats each reach as a discrete longitudinal segment, preserving the physical topology and preventing pixel-stepping errors in flow direction. The 5 km scale ensures that the resulting dataset is locally relevant, providing enough detail for regional water management while remaining computationally efficient for a 45-year simulation.

In situ observations and metadata

The observational backbone of this study consists of daily discharge data from 18 in situ gauging stations strategically distributed to cover the major drainage divides of the Yangtze (Fig. 1). This includes control stations on the mainstem, such as Cuntan, Yichang, and Datong, and key stations at the outlets of major tributaries like the Min, Jialing, and Han Rivers. Metadata for each station includes precise geographical coordinates, the drainage area, and the associated MERIT Hydro reach ID. We categorize these 18 nodes into 10 headwater stations, which constrain the initial water yield, and 8 downstream stations, which constrain the incremental runoff. These observations act as the physical anchors for our Two-Stage Hierarchical Bias Correction, providing the necessary constraints to solve the mass balance equation across the basin’s hierarchical structure.

The mizuRoute routing model

To transform the spatially distributed runoff into time-variant river discharge, we employ the mizuRoute model²³. mizuRoute is a stand-alone, vector-based routing tool specifically designed to bridge the gap between grid-based land surface models and vector-based river networks. The modeling process is bifurcated into two distinct stages: hillslope routing and channel routing. In the hillslope component, mizuRoute performs a spatial aggregation by mapping the gridded runoff (from ERA5 or VIC-BC) onto the unique local drainage areas of each reach. This mapping uses area-weighted interpolation to ensure that water mass is strictly conserved during the transfer from the grid to the vector domain. This step is crucial for the proposed framework, as it allows us to attribute discharge errors back to the specific grid cells that contribute to a gauged reach. For the channel routing component, we utilize the Impulse Response Function (IRF), based on a unit hydrograph approach. The IRF simulates the movement of water through the network by considering both translation and attenuation of the flood wave. The same mizuRoute configuration was used for all runoff scenarios, and no scenario-specific calibration of routing parameters was performed. Channel-routing parameters, including those controlling flow translation and attenuation, followed the globally applicable default parameterization adopted in the default model configuration²³. Holding the routing setup fixed enables the relative differences among the experiments to be associated with the runoff correction within this configured modeling system, although it does not demonstrate that the routing parameters themselves are optimal for individual reaches or events.

Two-stage hierarchical bias correction (TSH-BC)

The TSH-BC framework (Fig. 2) assumes that persistent volumetric discrepancies between simulated and observed discharge can be represented, to first order, through a multiplicative adjustment of the grid-based runoff. This adjustment is an observation-constrained source-term correction rather than a physical diagnosis of the causes of model error. It does not imply that all discharge discrepancies originate from runoff generation. Instead, the corrected runoff acts as an effective source term that may collectively represent uncertainties associated with meteorological forcing, runoff generation, routing parameters, human regulation, and discharge observations. The principal purpose of the framework is to allocate these observation-constrained adjustments according to the hierarchical structure of the river network. Corrections are first applied to headwater drainage areas and are subsequently isolated within local incremental areas between successive control stations. Once introduced into mizuRoute, the adjusted source volumes are consistently accounted for and propagated through the connected downstream network. The resulting consistency therefore refers to the longitudinal coherence of the adjusted runoff volumes within the configured routing system, rather than to closure of the complete terrestrial water balance.　This hierarchical source-term treatment provides a practical advantage over independent gauge-level post-processing, which may improve discharge locally without introducing corresponding adjustments to runoff generated within the intervening drainage areas. Such local corrections can therefore be inconsistent with the surrounding routed flow field. The TSH-BC formulation instead produces a spatially continuous, observation-constrained reconstruction while retaining the upstream–downstream relationships represented by the river network. The multiplicative correction is primarily intended to reduce persistent volumetric biases at monthly to multi-year timescales. It is not designed to explicitly resolve event-scale flow dynamics, including flood-wave timing, channel attenuation, or low-flow recession behaviour, which also depend on routing parameters, catchment storage, reservoir operations, and other model structural processes. .

The first stage targets headwater basins (Fig. 2), which are defined as catchments with no upstream gauging stations. For a headwater station i, we identify all the contributing 0.1° ERA5 runoff grid cells that drain into that specific node. The correction factor k is derived from the volumetric ratio between the observed and simulated discharge. Crucially, this correction is applied uniformly to all upstream runoff grids within the headwater sub-basin. By correcting headwaters first, we establish a valid upstream boundary condition. This step ensures that any errors found at downstream stations can be isolated to the downstream landscape, rather than being confused with errors propagated from the source regions.

For basins situated between two or more successive gauging stations, the TSH-BC framework employs a cascading downstream correction logic. The downstream bias is strictly isolated to the local incremental drainage area, i.e., the specific grid cells that contribute flow between an upstream calibrated station and the next downstream station. To perform this isolation, we use the already corrected runoff from Stage 1 (and any preceding Stage 2 steps) to drive mizuRoute and calculate the simulated discharge at the downstream station (Q_{sim, downstream}). Any remaining discrepancy at downstream station is then attributed solely to the local incremental runoff. The correction factor k is calculated by comparing the observed incremental volume against the simulated incremental volume. By only applying this factor to the local grids, we avoid “over-writing” the valid runoff from the headwaters. This process is repeated reach-by-reach down the river hierarchy, ensuring that the accuracy cascades through the entire network while respecting the mass balance of each sub-segment.

To explore the sensitivity of the correction to temporal variability, the TSH-BC framework is implemented using two distinct strategies for calculating the correction factor k.

• Long-term Volume Correction (hereafter, ERA5 BC long-term): A single, static diagnostic bias ratio ((beta)) is calculated for the entire study period (1980–2024). This factor represents the total multi-year volumetric bias ({beta}_{i} = frac{sum {Q}_{sim,i}}{sum {Q}_{obs,i}})) for the station i. This approach effectively corrects systematic errors in the reanalysis while preserving the original daily temporal dynamics of the ERA5 forcing.

• Multi-year Monthly Correction (hereafter, ERA5 BC monthly): Recognizing that hydrological biases often exhibit strong seasonality due to monsoon patterns or vegetation cycles, we also calculate 12 distinct monthly correction factors (k₁, k₂, …, k₁₂) for each station. Each monthly factor is derived by aggregating all observations and simulations for that specific month across the 45-year period. This allows the model to adjust for seasonally varying errors. For instance, correcting a dry-season baseflow deficit differently than a wet-season peak-flow bias, thereby providing a more nuanced alignment with the basin’s natural regime.

It is essential to distinguish between the diagnostic bias ratio ((beta)) and the multiplicative correction factor ((k)) applied to the runoff grids. For a given station (i), the diagnostic bias ratio is evaluated as ({beta}_{i}=sum {Q}_{sim,i}/sum {Q}_{obs,i}). However, the applied correction factor (k) ensures volume matching and is the inverse: ({k}_{i}=sum {Q}_{obs,i}/sum {Q}_{sim,i}=1/{beta}_{i}). For Stage 2 (local incremental areas), the correction factor strictly isolates the lateral inflow. The local target volume (({V}_{loc,target,j})) is calculated by subtracting the already-corrected upstream routed flow (({V}_{up,j})) from the downstream observation (({V}_{obs,j})). The local incremental correction factor is then defined as: ({k}_{loc,j}=left({V}_{obs,j}-{V}_{up,j}right)/left({V}_{sim_all,j}-{V}_{up,j}right)), where the denominator represents the simulated local contribution.

Experimental design and evaluation metrics

To evaluate the effectiveness of TSH-BC, we designed four comparative simulation scenarios: ERA5 (uncorrected), ERA5 BC Monthly (seasonal correction), ERA5 BC Long-term, and the VIC-BC Benchmark. All simulations cover the period from 1980 to 2024. The performance is quantified using four widely recognized metrics. The Kling-Gupta Efficiency (KGE) balances correlation, variability, and bias. The Nash–Sutcliffe Efficiency (NSE) evaluates the ability to capture high-flow magnitudes. The Logarithmic NSE (logNSE) is calculated to emphasize performance during low-flow and drought periods, which is critical for water security analysis. Finally, Relative Bias (RB) measures the systematic percentage error in total water volume. These metrics provide a comprehensive evaluation of the TSH-BC framework’s ability to reconstruct the Yangtze’s hydrology across multiple scales and stream orders.

Overall performance and spatiotemporal patterns of correction factors

The evaluation of the TSH-BC framework across the 18 stations reveals a substantial improvement in hydrological simulation accuracy, as summarized in the statistical distributions shown in Fig. 3. In the baseline ERA5 simulation, the median KGE is only 0.45, with the interquartile range (IQR) spanning from 0.37 to 0.59. This poor performance is characterized by significant outliers in the upper Yangtze, such as Panzhihua, where the KGE is − 0.07, indicating that the raw reanalysis fails to provide even a basic representation of the observed flow. This systematic failure is primarily attributed to a massive wet bias in the high-altitude headwaters. Figure 3d illustrates this through the RB metric, where ERA5 shows a median overestimation of 40.1%, peaking at over 88% at Panzhihua. The bias ratio of 1.89 at this station suggests that the ERA5 significantly overestimates runoff generation on the Tibetan Plateau, likely due to the misrepresentation of complex orographic precipitation or the simplified treatment of snowmelt and permafrost dynamics.

By applying the TSH-BC framework, the accuracy of the simulation is substantially elevated. The ERA5 BC monthly scenario exhibits a median KGE of 0.74, which is similar to the state-of-the-art VIC-BC benchmark. However, a critical distinction lies in the volumetric consistency; while VIC-BC exhibits a negative median bias (− 4%), our monthly and long-term correction scenarios effectively bring the median RB to near-zero levels (1% and 2%, Fig. 3d). This correction is visually evident in the narrowing of the boxplots in Fig. 3, representing a more robust and spatially consistent simulation across the basin. The improvement in the NSE from a median of 0.20–0.64 (Fig. 3b) further confirms that the hierarchical bias correction not only captures the volume but also significantly improves the model’s ability to simulate the timing and magnitude of the Yangtze’s monsoon-driven flood peaks.

The spatial heterogeneity of the correction factors (k) provides deep insight into the regional biases of the ERA5 reanalysis. While we selected representative stations for detailed analysis (shown in Fig. 4), the overall trend across the 18 stations (not shown here) shows that the Yangtze mainstem requires more complex seasonal adjustments than its lower tributaries. At the Panzhihua headwater station, which monitors the high-altitude Jinsha River, the long-term bias ratio is 1.89, revealing that the raw reanalysis produces nearly double the observed annual water volume. This overestimation is not temporally uniform; during the dry season months of November and December, the bias ratios spike to 2.65 and 2.69, respectively. This suggests that the land surface model within ERA5 fails to accurately simulate the recession processes or overestimates groundwater contributions during the winter months. In the middle and lower reaches, the mainstem bias ratios at Cuntan (1.51) and Yichang (1.49) remain elevated, yet they reflect the cumulative effect of the corrected upstream flows. By implementing the monthly correction factor, we can observe the TSH-BC’s ability to normalize these seasonal deviations. At Cuntan, the monthly ratios are brought significantly closer to unity, ranging from 0.95 in June to 1.28 in November, thereby smoothing the artificial seasonality introduced by the raw forcing. The long-term correction, while effective at a multi-year scale, often fails to address these critical monthly discrepancies, leading to remaining seasonal biases that the monthly framework successfully mitigates.

Evaluation of reach-scale hydrographs and flow dynamics

The reconstruction of daily discharge time series provides a granular view of how the hierarchical correction improves flow dynamics across different stream orders. For the four representative stations (i.e., Panzhihua, Cuntan, Wulong, and Yichang shown in Fig. 4), the hydrograph analysis reveals that the simulation from ERA5 consistently overestimates both flood peaks and baseflows, particularly in the upper Yangtze. At Panzhihua, the KGE for ERA5 is actually negative (− 0.06), indicating weak overall agreement with the observed discharge. Nevertheless, this value remains above the mean-flow benchmark of approximately − 0.41 and therefore indicates an improvement over using the observed mean flow as a constant predictor²⁴. However, after the headwater correction, the monthly BC KGE soars to 0.9, representing a near-perfect alignment with the observed hydrograph’s volume and timing (Fig. 4a). This improvement is also evident in the Cuntan and Yichang mainstem stations, where KGE values rise from 0.45 to 0.86 and 0.84, respectively (Fig. 4b, d). The visual evidence in the hydrographs confirms that the monthly correction successfully restores the recession limb of the hydrograph, which is typically over-sustained in the raw reanalysis. This is particularly noticeable during the transition periods between the wet and dry seasons, where the TSH-BC framework captures the rapid drop in discharge following the monsoon’s cessation. By effectively reconciling the water balance at a monthly resolution, the framework mitigates the flashy response of the raw reanalysis to minor precipitation events while preserving the magnitude and translation speed of major flood waves. The alignment of the monthly BC hydrograph during extreme years, such as the 1998 floods, demonstrates the framework’s robustness in capturing high-magnitude events across large-scale drainage networks.

Furthermore, the logNSE results emphasize the framework’s superior performance during low-flow periods by penalizing volumetric errors more heavily at lower magnitudes. The median logNSE across the basin increases from 0.36 in ERA5 to 0.72 in the monthly BC scenario, a significant jump compared to the 0.62 achieved by the long-term correction. This improvement is critical because it suggests that the monthly factors better represent the nonlinear relationship between precipitation and baseflow during the dry season. At stations like Wulong on the tributary (Fig. 4c), the logNSE gain is particularly pronounced, as the monthly correction compensates for the systematic over-prediction of winter baseflow common in ERA5. Because the correction is applied to the source runoff before it enters the 40,000-reach network, this correction is inherited by even the smallest ungauged tributaries. This should ensure that the entire river hierarchy maintains a more realistic representation of flow regimes, assuming a reasonably representative model structure, while maintaining consistency between local water balances and downstream connectivity..

Impacts of hierarchical local area correction

The isolation of incremental drainage area errors through Stage 2 of the TSH-BC framework is an important factor in maintaining downstream simulation accuracy and physical consistency. This hierarchical logic ensures that errors identified at a downstream station are attributed exclusively to the local sub-segment of the river network, rather than simply adjusting the already-validated upstream flow. To evaluate the efficacy of this stage, we analyzed the performance gains across eight key downstream stations representing both the mainstem and major tributaries. As shown in Fig. 5, the incremental correction provides substantial performance leaps across all metrics, providing evidence that downstream biases are associated not only with headwater contributions but also with processes and lateral inflows within the middle reaches. Looking spatially from upstream to downstream, the framework first stabilizes the mainstem at Zhutuo, which receives flow from the Panzhihua headwater and the Gaochang tributary. Despite a vast local drainage area of 274,575 km², the local correction at Zhutuo yields a KGE Gain of 0.14 and an NSE Gain of 0.13 from our monthly BC. Moving further downstream, Cuntan (local area: 14,649 km²) and Yichang (local area: 55,581 km²) show continued improvements, with KGE Gains of 0.09 and 0.07 from our monthly BC, respectively. The smaller gains at these mid-reach mainstem points compared to Zhutuo suggest that while local runoff biases exist, the overall accuracy is increasingly dominated by the high-quality upstream boundary conditions established in Stage 1. This cascading accuracy validates the hierarchical structure: once the massive headwater volumes are corrected, the downstream correction can focus on refining the specific hydrological signatures of local lateral inflows.

The tributary systems, particularly the Han River (i.e., river reaches covering Baihe Xiangyang, and Huangzhuang stations), further demonstrate the necessity of local-scale refinement. The Han River represents a heavily managed and hydrologically complex sub-basin. At Xiangyang (local area: 46,890 km²), the logNSE gain reaches a remarkable 0.22 from our monthly BC, the highest among all downstream stations. This dramatic improvement reflects the framework’s ability to correct for complex local water balances that the ERA5 simulation fails to capture, particularly the low-flow dynamics. Similarly, at Huangzhuang (local area: 24,222 km²), we observe an NSE gain of 0.11. By isolating these reaches, the TSH-BC framework ensures that the unique biases of the Hanjiang do not contaminate the mainstem simulation until the physical confluence.

As we reach the middle and lower mainstem, the complexity increases as Yichang and the Han tributary converge at Hankou. With a massive local incremental area of 326,249 km², Hankou exhibits substantial KGE and NSE gains of 0.13 and 0.16 from our monthly BC, respectively. These gains are critical, as they account for the enormous lateral inflows from the middle reaches of the Yangtze, including the vast humid regions surrounding the Dongting Lake. Finally, at the Datong station, the ultimate basin outlet which aggregates flow from Jiujiang and three major southern tributaries (Waizhou, Meigang, and Lijiadu), the framework achieves an NSE gain of 0.19 and a KGE gain of 0.13. The fact that the gains using monthly correction at Datong are consistently higher than those using long-term correction (e.g., NSE gain of 0.19 vs. 0.17) underscores that the local correction must be seasonally adaptive to manage the monsoon-driven variance of the terminal drainage area. Ultimately, this segment-by-segment solution of the mass balance equation prevents the “over-correction” trap common in lumped models. In a non-hierarchical system, a model might achieve high accuracy at Datong by using physically unrealistic runoff values in the Tibetan headwaters to compensate for lowland errors. In contrast, the TSH-BC framework forces each reach to be correct. By the time the water reaches the sea, the framework has reconciled the water balance across a diverse hierarchy of stream orders and local climates. This improvement enables the TSH-BC framework to produce a hydrologically coherent dataset that performs better across the evaluated metrics and provides a useful foundation for reach-level water management across the entire 1.8 million km² basin.

The comprehensive reconstruction of the Yangtze River Basin’s hydrological history, as visualized in the reach-level discharge maps in Fig. 6, demonstrates the capacity of the TSH-BC framework to generate high-fidelity datasets that are both spatially continuous and locally relevant. By synthesizing the MERIT Hydro vector river network which provides a globally consistent 5-km routing topology with the ERA5 reanalysis runoff, we have established a high-resolution simulation framework that remains updated to the present day. Unlike traditional global modeling efforts that often struggle with regional biases, the results shown in Fig. 6 represent a functional bridge between 0.1° global atmospheric forcing and 5-km reach-level hydrology. This approach is designed for immediate local application; while global products often suffer from regional biases, our method allows local water authorities to utilize their own in situ records or publicly available datasets like GRDC and the CAMELS series^25,26 to anchor global models to local realities. The ability to bridge the gap between 0.1° global atmospheric forcing and 5-km reach-level hydrology ensures that the final product is not just a scientific exercise but a functional tool for regional water management.

A primary strength of this methodology lies in its unprecedented efficiency in data utilization and its respect for physical connectivity. Unlike traditional post-processing techniques that treat river stations as isolated points, the TSH-BC framework maintains the longitudinal integrity of the river system. By bias-correcting the runoff at the source grid-cell level before it enters the routing scheme, we ensure that the entire river hierarchy from the smallest headwater streams to the terminal outlet inherits a topology-aware flow regime. This is a significant advancement over simple “end-of-pipe” statistical corrections, which often achieve numerical accuracy at a gauge by sacrificing the mass balance of the upstream network. Our segment-by-segment hierarchical approach forces the simulation to be correct, accounting for the hydrodynamic relationships between upstream reaches and downstream lateral inflows.

Furthermore, the framework demonstrates that high-fidelity historical series can be constructed with minimal local data. While this study utilized 18 gauging stations to maximize accuracy across the Yangtze, the methodology is inherently scalable. In basins where gauging stations are abundant, the correction process can be further refined; in data-scarce regions, even a single downstream gauge can significantly constrain the simulation through our hierarchical logic. This simplicity and transferability mean the TSH-BC framework can be deployed to virtually any river basin globally, offering a robust alternative to complex deep-learning models that require decades of training data. By leveraging the near-real-time update cycle of ERA5, we provide a live historical record that supports immediate decision-making for contemporary hydrological challenges.

Looking forward, there are several avenues to further elevate the precision of this work. First, future iterations could utilize even finer hydrographic datasets, such as the Hydrography90m network²⁷, to model river dynamics at an even higher spatial resolution. However, shifting to such a granular scale will require advanced adjustments in routing algorithms to manage increased computational complexity and potential DEM-induced artifacts. Second, the current retrospective framework can be extended into a predictive mode. By applying the calculated monthly bias-correction factors to short-term and seasonal meteorological forecasts, we can provide high-skill, reach-level discharge predictions.

The resulting 45-year reach-level dataset (Fig. 6) is more than just a record of the past; it is a foundational resource for diverse applications. From assessing the environmental flow requirements of endangered species in small tributaries to diagnosing the propagation of flood waves during catastrophic monsoons, the spatial continuity of this data offers insights that point-based observations cannot. As we enter an era of increasingly volatile climate extremes, such high-resolution, physically consistent hydrological reconstructions will be indispensable for building resilient water infrastructure and managing the transboundary water resources of our planet’s great river systems.

While the TSH-BC framework substantially improves the volumetric agreement of the historical river flow reconstruction, several assumptions and limitations should be recognized. First, the correction factors and evaluation metrics were derived from the same available record for 1980–2024. This full-period, in-sample design was selected because the objective of the present study is to construct an observation-constrained historical baseline using all available information, rather than to develop a forecasting or extrapolation model. The use of the complete record is consistent with recent evidence that final hydrological model parameterization can benefit from maximizing the information content of the available observations²⁸. Nevertheless, the resulting performance metrics represent reconstruction agreement rather than independent predictive skill. In particular, the near-zero relative bias is an expected consequence of the volume-ratio correction and should not be interpreted as independent evidence of model skill. Temporal extrapolation, spatial transferability, and performance in ungauged basins remain to be independently evaluated. Second, the framework operationally represents the residual discharge discrepancy through an adjusted runoff source term. It therefore does not independently close the complete terrestrial water balance or identify the individual causes of model error. Uncertainties associated with meteorological forcing, runoff generation, routing parameters, model structure, and discharge observations may be jointly represented in the inferred correction factors. Streamflow observations are used as reference constraints without an explicit representation of rating-curve or measurement uncertainty. Consequently, improved discharge agreement may partly result from compensation among different sources of error. Third, the gauge observations contain the effects of human interventions, including reservoir regulation, lake storage changes, water transfers, abstractions, and return flows. The present routing configuration does not explicitly simulate reservoir-specific operations or consumptive water use. The resulting product should therefore be interpreted as an observation-constrained historical flow reconstruction rather than a reconstruction of purely naturalized discharge. The correction factors may absorb both natural hydrological errors and anthropogenic alterations into the effective runoff source term. Explicit separation of these influences would require time-varying reservoir storage and release observations, historical operating policies, water-use records, and additional process constraints. Because reservoir operations differ among facilities and evolve in response to multiple management objectives, the use of a generic operating scheme without adequate historical information could introduce additional uncertainty. Future work should therefore couple the hierarchical correction framework with reservoir-specific operating information, water-use modules, and process-oriented routing evaluations. Such extensions would support a clearer attribution of discharge variability to natural runoff generation, river routing, and human regulation.

This study presents the Two-Stage Hierarchical Bias Correction (TSH-BC) framework for observation-constrained, reach-level reconstruction of historical river discharge across large river networks. Rather than correcting discharge independently at individual gauging stations, TSH-BC represents persistent discharge discrepancies through adjustments to the runoff source term in headwater and local incremental drainage areas. These adjusted runoff volumes are subsequently routed through the connected river network, thereby maintaining longitudinal coherence between upstream contributions, local lateral inflows, and downstream discharge within the configured routing system.

Application to the Yangtze River Basin improved discharge agreement across the evaluated gauging stations. The median KGE increased from 0.45 for the uncorrected ERA5 simulation to 0.74 for the monthly correction, while the median relative bias decreased from 40.1% to near zero. The hierarchical treatment of local incremental drainage areas further improved the evaluated performance metrics at downstream stations by preventing downstream corrections from being indiscriminately applied to already adjusted upstream drainage areas. These results demonstrate the value of explicitly accounting for river-network hierarchy when constructing a spatially continuous historical discharge reconstruction.

The reported performance represents full-period reconstruction agreement based on the available station-specific observations, rather than independent predictive or spatial-transferability skill. In addition, the inferred runoff adjustments may collectively represent errors associated with meteorological forcing, runoff generation, routing parameters, human regulation, and discharge observations. The resulting dataset should therefore be interpreted as an observation-constrained historical flow reconstruction rather than a purely naturalized discharge product.

The modular structure of TSH-BC provides a practical basis for application to other large river systems where sufficient discharge observations and a consistent river-network representation are available. Future work should evaluate temporal extrapolation and spatial transferability, incorporate reservoir operations and water-use processes, and examine event-scale flow characteristics under alternative routing configurations. Within these stated limitations, the framework offers a scalable approach for integrating global runoff forcing with regional observations to produce spatially continuous historical discharge estimates across connected river networks.