Coupling WRF with HEC-HMS and WRF-Hydro for flood forecasting in typical mountainous catchments of northern China

The atmospheric–hydrological coupling systems are essential to flood forecasting because they allow for more improved and comprehensive prediction of flood events with an extended forecast lead time. Achieving this goal requires a reliable hydrological model system that enhances both rainfall predictions and hydrological forecasts. This study evaluates the potential of coupling the mesoscale numerical weather prediction model, i.e., the weather research and forecasting (WRF) model, with different hydrological modeling systems to improve the accuracy of flood simulation. The fully distributed WRF-Hydro modeling system and the semi-distributed Hydrological Engineering Center Hydrological Modeling System (HEC-HMS) were coupled with the WRF model, and the lumped HEC-HMS model was also adopted using the observed gauge precipitation as a benchmark to test the model uncertainty. Four distinct storm events from two mountainous catchments in northern China characterized by varying spatial and temporal rainfall patterns were selected as case studies. Comparative analyses of the simulated flooding processes were carried out to evaluate and compare the performance of the coupled systems with different complexities. The coupled WRF–HEC-HMS system performed better for long-duration storm events and obtained optimal performance for storm events uniformly distributed both temporally and spatially, as it adapted to more rapid recession processes of floods. However, the coupled WRF–HEC-HMS system did not adequately capture the magnitude of the storm events as it had a larger flow peak error. On the other hand, the fully distributed WRF–WRF-Hydro system performed better for shorter-duration floods with higher flow peaks as it can adapt to the simulation of flash floods. However, the performance of the system became poor as uniformity decreased. The performance of the lumped HEC-HMS indicates some source of uncertainty in the hydrological model when compared with the coupled WRF–HEC-HMS system, but a larger magnitude error was found in the WRF output rainfall. The results of this study can help establish an adaptive atmospheric–hydrologic coupling system to improve flood forecasting for different watersheds and climatic characteristics.

1 Introduction

Floods are frequent and widespread natural hazards that result in substantial annual losses to human lives and properties worldwide (Jonkman, 2005). Due to climate change, the future is expected to bring more intense precipitation, which might potentially lead to an increase in extreme-rainfall-induced flood events and elevated flood risk (Mirza, 2003). Flood forecasting is essential to mitigating the impact of floods by providing timely warnings and enabling proactive measures that help to safeguard lives, property, and infrastructure in vulnerable areas (Merz et al., 2020). Improving the ability to predict flood risks ahead of time is essential to the premise of promoting forecast accuracy. To improve the simulation accuracy and extend the forecast lead time, there is a growing trend in favor of substituting the conventional throughfall with the mesoscale numerical weather prediction (NWP; Ozkaya, 2023; Trinh et al., 2023; Kaufmann et al., 2003). An effective strategy to do this involves coupling the hydrological model with a high-resolution regional NWP model. This approach has demonstrated capabilities to not only improve the accuracy of flood forecasting but also extend real-time forecast lead time, compared to conventional flood forecasting that relies only on gauge observations as inputs (Seid et al., 2021).

In recent years, coupling hydrological models with high-resolution NWP models covering different scales has emerged as a promising approach to improve flood simulation (Jasper et al., 2002; Bartholmes and Todini, 2005; Nam et al., 2014; Wu et al., 2014, 2020; Cattoën et al., 2016; Li et al., 2017; Chen et al., 2020; Ming et al., 2020; Dasgupta et al., 2023; Patel and Yadav, 2023; Bacelar et al., 2023). Jasper et al. (2002) coupled the WaSim-ETH model with surface observations, forecast data from five high-resolution NWP models, and weather radar data for seven extreme flood events. They concluded that future simulation improvement hinged on the NWP model development. Bartholmes and Todini (2005) analyzed the effects of coupling meteorological mesoscale quantitative precipitation forecasts at various scales with the TOPKAPI model to extend the flood simulation horizon. The results highlighted the limited reliability of quantitative simulation precipitation generated by meteorological models. Cattoën et al. (2016) coupled the NZLAM regional NWP model at high and low resolutions with the TopNet hydrological model for flood simulation, and their findings indicated the advantage of utilizing a high-resolution convection-permitting time-lagged ensemble simulation over a lower-resolution large-scale model. Li et al. (2017) coupled the Weather Research and Forecasting quantitative precipitation forecast (WRF QPF) with the Liuxihe model to extend the flood-forecasting lead in southern China and found that as the lead time increased, both the accuracy of the WRF QPF and the flood simulation capability decreased. Chen et al. (2020) coupled the GRAPE_MESO model, a two-dimensional hydrodynamical flood model, and a rainstorm prediction reconstruction method for urban flood simulation. Their results showed that the coupled modeling system achieved accurate predictions with high resolution and an extended lead time. Ming et al. (2020) produced high-resolution catchment-scale rainfall–runoff and flood forecasting by coupling the WRF model with a GPU-accelerated hydrodynamic model. The system provided a 34 h lead time based on weather forecasts available 36 h in advance. Patel and Yadav (2023) researched the coupling of hybrid ensemble linear regression, the HEC-HMS model, and the Bayesian numerical weather model to simulate hourly reservoir inflows. Their results demonstrated the effectiveness of these coupled systems in accurately predicting reservoir inflows in the Sabarmati River basin in India.

Although recent studies have been conducted to improve flood forecasting by coupling NWP models with hydrological models, few have addressed the implications of choosing between fully distributed and semi-distributed models of varying complexities in constructing these coupling systems. The complexity of the hydrological model plays an important role in determining the generation of the streamflow, which should not be neglected when establishing and evaluating the atmospheric–hydrological coupling system (Ahmed et al., 2023). A fully distributed model divides a watershed into smaller spatial units, allowing for a detailed representation of the entire area. A semi-distributed model groups similar sub-basins, providing a balance between detail and computational efficiency (Valiya Veettil et al., 2021). By integrating meteorological data from the NWP models into these hydrological models, it is possible to create a holistic understanding of how meteorological inputs impact the generation of the streamflow. Understanding the source of uncertainties involved in this process and how to eliminate them is also of paramount importance in improving the performance of the atmospheric–hydrological systems for flood forecasting. Extreme weather events, particularly intense storms, pose significant challenges for hydrological modeling due to their complex interactions with surface and subsurface processes.

The main objective of this study is to evaluate the potential of coupling the mesoscale numerical weather prediction model, i.e., the weather research and forecasting (WRF) model, with different hydrological modeling systems to improve the accuracy of flood forecasting. The fully distributed WRF-Hydro modeling system and the semi-distributed Hydrological Engineering Center Hydrological Modeling System (HEC-HMS) were coupled with the WRF model. Additionally, the lumped HEC-HMS model was adopted using observed gauge precipitation as a benchmark to test the rainfall input uncertainty. This approach examines the effectiveness of the HEC-HMS model in replicating observed conditions at the gauge locations as shown in Fig. 1. Four distinct storm events, characterized by varying spatial and temporal rainfall patterns, were selected as case studies. These events occurred in two mountainous catchments along the Daqing River, where precise flood prediction is urgently required to mitigate the risks associated with construction in northern China's downstream area. The WRF model stands out as the predominant mesoscale NWP model for simulating and forecasting rainfall in hydrology and water resource-related disciplines (Done et al., 2004; Lo et al., 2008; Liu et al., 2012; Haghroosta et al., 2014; Chawla et al., 2018; Yáñez-Morroni et al., 2018; Huang et al., 2023). The advancement of the WRF-Hydro modeling system, built upon the research on the WRF model, has enhanced the efficiency of utilizing WRF for hydrological simulation (Gochis et al., 2013). This innovation addresses the issue of misalignment between the resolution of the atmospheric model and the resolution of the hydrological model. In recent years, the development of the WRF-Hydro modeling system by the National Center for Atmospheric Research (NCAR) and its collaboration partners has been coupled with the WRF model in various hydrological research projects. This coupling has enhanced flood-forecasting accuracy (Senatore et al., 2015) and improved the representation of streamflow dynamics (Ryu et al., 2017). It has demonstrated effectiveness in simulating extreme weather events and their hydrological impacts (Wang et al., 2020; Sun et al., 2020) and in capturing spatial variability in hydrological responses (Quenum et al., 2022; Wang et al., 2022). Additionally, it has proven robust across various climates and geographies, ensuring reliable hydrological predictions (Liu et al., 2023; Naabil et al., 2023). The HEC-HMS model is a widely used hydrological model with a flexible structure to be built in either a lumped mode or a semi-distributed mode. In recent years, there have also been efforts to couple HEC-HMS with the WRF model (Herath et al., 2016; Givati et al., 2016; Niyogi et al., 2022; Tien Thanh et al., 2023; Ting et al., 2023). The study of Herath et al. (2016) demonstrates the coupled WRF and HEC-HMS to be a helpful tool for the flood warning system in Polgolla Barrage in Sri Lanka. Givati et al. (2016) coupled the HEC-HMS model using a 3 km hourly precipitation generated by the WRF model for flood forecasting in the Mediterranean region. Niyogi et al. (2022) also coupled the WRF model of 1.5 km resolution with the HEC-HMS model and HEC-RAS 2D hydraulic model of 10 m resolution.

Figure 1 Locations of the two study sub-catchments in the Daqinghe catchment.

This study establishes two atmospheric–hydrological systems by coupling the WRF model with WRF-Hydro and the HEC-HMS model to forecast four typical storm events in the study basin. This study utilizes a one-way coupling approach between the two hydrological model structures and the WRF model. This means that WRF output drove the hydrological models without reciprocally influencing the atmospheric modeling processes. The 1 km×1 km output rainfall from the WRF model is used to drive the WRF-Hydro model and the HEC-HMS model to produce flood forecasting. Also, the lumped HEC-HMS driven by the observed garage precipitation is used to produce flood forecasting. Comparative analyses were carried out to evaluate the performance of the forecast processes of the four storm events by the coupled atmospheric–hydrological and lumped systems and to identify the source of uncertainties further.

The analysis in Sect. 4 is structured as follows:

• in Sect. 4.1, we evaluate and compare the performance of the coupled systems, i.e., WRF–WRF-Hydro and WRF–HEC-HMS, with different complexities;

• in Sect. 4.2, we evaluate and compare the performance of the coupled WRF–HEC-HMS and lumped HEC-HMS model driven by observed rainfall to analyze the model uncertainty;

• in Sect. 4.3, we evaluate and analyze the error of the WRF model output rainfall and its resulting uncertainty in the atmospheric–hydrological coupling systems.

2 Study area and events

2.1 Study area

This study adopts the Fuping and Zijingguan sub-catchments within the Daqing River basin, located in northern China, as study areas. Fuping (2219 km²) and Zijingguan (1760 km²) are typical mountainous sub-catchments located in the upper Shahe River of the southern branch and the upper Juma River of the northern branch of the Daqing River, respectively, as shown in Fig. 1. The Fuping sub-catchment has a longitudinal river slope of 5.7 % and a residential area of 0.63 %. The Zijingguan sub-catchment has a longitudinal river slope of 5.5 % and a residential area of 0.52 %. The predominant land use in the Daqing River basin is farmland, forestland, and grassland, with a granitic gneiss type of geology. The Fuping catchment has a total of 8 gauged stations, and the Zijingguan catchment has a total of 11 gauged stations as shown in Fig. 1. Hydrological stations measured the flow at the outlets of the two catchments. The Daqing River basin experiences severe soil erosion attributed to dry soil conditions and excessive groundwater exploitation. Additionally, during the storm season, typically from June to September, the river undergoes substantial seepage. The mean annual rainfall is approximately 490 and 650 mm for Fuping and Zijingguan, respectively, with most rain occurring between late May and early September. Summer storms with high intensities and short durations are typical of the rainfall found in China's mountainous regions, such as Fuping and Zijingguan. Therefore, they mostly result in severe flood disasters in the downstream Daqing River basin.

2.2 Storm events

Four storm events with 24 h durations and relatively high flow peaks, as shown in Fig. 2, are selected to test the performance of the coupled hydrological rainfall–runoff modeling systems constructed in this study. Three events occurred in the Fuping sub-catchment, and one event occurred in the Zijingguan sub-catchment. Storm events in the Fuping and Zijingguan catchments in the Daqinghe basin are driven by the East Asian monsoon bringing moist air and intense rainfall. These storms often form into mesoscale convective systems (MCSs), which are large clusters of thunderstorms with sustained heavy rain, leading to rapid river rises and potential flooding. The gauged rainfall and flow data were provided by the Ministry of Water Resources of the People's Republic of China. Table 1 shows the duration, cumulative rainfall, and peak discharges of the storm events. The four 24 h storm events are categorized according to their spatial and temporal distributions. The coefficient of variance (C_v) of the storm event is calculated to designate the different homogenous characteristics (Hosking and Wallis, 1997) as

$$\begin{array}{}\text{(1)}& C_{\mathrm{v}}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N{\left(\frac{x_i}{\stackrel{\mathrm{‾}}{x}}-\mathrm{1}\right)}^{\mathrm{2}}}{N}}}.\end{array}$$

When calculating C_v for the spatial distribution, the 24 h cumulative rainfall at any ith rain gauge is x_i, the cumulative average rainfall of all stations is $\stackrel{\mathrm{‾}}{x}$, and N is the total number of all the stations. For C_v, in temporal distribution calculation, the catchment areal rainfall at any hour i time step is x_i, the average areal rainfall of all time steps is $\stackrel{\mathrm{‾}}{x}$, and N is the total duration of the storm events (24 h).

Figure 2 The rainfall–runoff observations of the four 24 h storm events.

Table 1 The selected four 24 h storm events in Fuping and Zijingguan sub-catchments.

The C_v values of the selected storm events shown in Table 2 reflect the spatiotemporal derivation of the catchment accumulative gauge rainfall of each station and the average rainfall at each time step, respectively. A smaller C_v value indicates that the rainfall distribution is more uniform or even in space or time. As seen in Table 1, storm event 2, with the smallest spatial distribution C_v value, is the most uniform in space, followed by event 1, event 4, and event 3. The most uniform in time is storm event 1, with the smallest temporal distribution C_v value, followed by event 2, event 4, and event 3. Generally, from the categorization, event 1 has rainfall that is uniform in spatial and temporal distributions, and event 2 has rainfall that has a uniform spatial distribution but non-uniform temporal distribution. Events 3 and 4 have non-uniform rainfall in spatial and temporal distributions.

Table 2 The coefficient of variance (C_v) of the four 24 h storm events.

The selected rainfall–runoff storm events have different lengths and flood recession times. We employ a 71 h duration for events 1 and 3, a 67 h duration for event 2, and a 36 h duration for event 4, as shown in Fig. 2. Event 4 can be noted as an extreme situation; it has larger C_v values and occurred in a shorter duration with the highest flow peak and rainfall intensity. Event 4 corresponds to the events at the most significant monitoring point that occur once every 500 years and is regarded as one of the largest flood disasters (Du et al., 2016). These different rainfall–runoff characteristics are used in the evaluation of the simulation results.

3 Hydrological models and experimental design

3.1 WRF model setup

The Weather Research and Forecasting (WRF) non-hydrostatic model is a widely used numerical weather prediction system that simulates and forecasts atmospheric processes (Powers et al., 2017). The WRF model version 3.7 is used in this study. Its structure includes initialization, dynamics, physics, grid options, output tools, boundary conditions, nesting capabilities, and parallel-computing support, making it versatile for various applications and domains. The versatile and flexible abilities of the WRF model predict weather patterns at different spatial and temporal scales, from global to regional and even local levels (Cassola et al., 2015). Additional details about the WRF model can be explored in Skamarock and Klemp (2008). The performance of the WRF model largely depends on the parameterization schemes, which can be effective for some storm events but not for others (Liu et al., 2013). Due to the difficulty in determining the best schemes for future storms, these are often preset in operational uses (Liu et al., 2015). In this study, we utilize the most widely used physical parameterizations for northern China. Details of the parameterizations that significantly influence precipitation generation are provided in Table 3 and further elaborated by Tian et al. (2017a).

Table 3 Main WRF model physical schemes used in this study.

This study employs the most widely used parameter setup option from previous studies in northern China because determining the optimal parameters for future storm events poses a significant challenge (Tian et al., 2017b). More information on the parameter setup can be explored in Tian et al. (2017a). Table 4 shows the main parameterization configurations of the WRF model for the two sub-catchments that have more influence on precipitation generation. The initial boundary conditions for simulation are derived from the 1°×1° FNL driving data at 6 h intervals, with the integration time step set at 6 s (Zhu et al., 2022). The FNL data are the Final Operational Global Analysis meteorological data, which are provided by the National Centers for Environmental Prediction (NCEP; available at http://rda.ucar.edu/datasets/ds083.2/, last access: 25 January 2023). The WRF output data interval is set at 1 h with a spin-up time of 6 h. Three nested domains are set up over the Fuping and Zijingguan sub-catchments. The innermost domain of the WRF is set up at a 1 km horizontal resolution, and the nesting ratio of the three layers is configured at 1:3. The grid center of the Fuping sub-catchment is at $\mathrm{39}\mathit{°}{\mathrm{04}}^{\prime }{\mathrm{15}}^{\prime \prime }$ N and $\mathrm{113}\mathit{°}{\mathrm{59}}^{\prime }{\mathrm{26}}^{\prime \prime }$ E, and the nesting grid divisions from domain 1 to domain 3 are 252 km×234 km, 144 km×126 km, and 96 km×84 km. The grid center of the Zijingguan sub-catchment is at $\mathrm{39}\mathit{°}{\mathrm{25}}^{\prime }{\mathrm{59}}^{\prime \prime }$ N and $\mathrm{114}\mathit{°}{\mathrm{46}}^{\prime }{\mathrm{01}}^{\prime \prime }$ E, and the nesting grid divisions from domain 1 to domain 3 are 216 km×198 km, 108 km×90 km, and 72 km×42 km. With a Lambert projection, a 40 vertical discretization up to a 50 hPa top-layer pressure is set up for the three nested domains (Tian et al., 2020). The downscaled output precipitation from the WRF model serves as input to drive the HEC-HMS and the WRF-Hydro models.

Table 4 Main WRF model parameterization configurations for the two sub-catchments.

3.2 WRF-Hydro model

The Weather Research and Forecasting Hydrological (WRF-Hydro) model is a widely used, fully distributed hydrological modeling system that integrates atmospheric, land surface, and hydrological processes to simulate and predict surface and subsurface water fluxes. Its structure is designed to represent the complex interactions between the atmosphere and the land surface, including precipitation, runoff, streamflow, and soil moisture. The WRF-Hydro model can only be run by coupling with the WRF model or utilizing meteorological data to establish an atmospheric–hydrological model system. This study implements a one-way run utilizing the WRF-Hydro modeling system version 3.0 with the WRF model (Gochis et al., 2015).

The hydrological components of WRF-Hydro model are described below.

• Routing. This simulates water movement through river networks and channels, accounting for flow routing and storage dynamics. WRF-Hydro uses a simplified Muskingum–Cunge routing equation for river routing,

  $$\begin{array}{}\text{(2)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial Q}{\partial t}}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{\Delta }x}}(Q^n-Q^{n-\mathrm{1}})\\ {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{S^n+S^{n-\mathrm{1}}}{\mathrm{2}}}\right)\left({\displaystyle \frac{\mathrm{\Delta }{Q}^n}{\mathrm{\Delta }t}}\right),\end{array}\end{array}$$

  where Q is the river discharge, t is time, Δx is the river reach length, S is the channel storage, n and n−1 represent the current and previous time steps, and ΔQ is the change in discharge over time $\left(\frac{\partial Q}{\partial t}\right)$.

• Runoff generation. The WRF-Hydro model generates runoff using a simple water balance (SWB) method. The topsoil layer experiences a surface infiltration excess when the precipitation capacity surpasses the infiltration capacity, resulting in a corresponding alteration in the surface water depth h (m),

  $$\begin{array}{}\text{(3)}& {\displaystyle }\begin{array}{rl}{\displaystyle }& {\displaystyle \frac{\partial h}{\partial t}}=\\ {\displaystyle }& {\displaystyle \frac{\partial p_e}{\partial t}}\left\{\mathrm{1}-{\displaystyle \frac{\left\{\begin{array}{c}\left[{\sum }_{K=\mathrm{1}}^{\mathrm{4}}{Z}_i({\mathit{\delta }}_s-{\mathit{\delta }}_k)\right]\\ \times \left[\mathrm{1}-\exp \left(-s\frac{R_{\text{dt}}}{R_{\text{fd}}}\frac{\mathrm{\Delta }t}{\mathrm{86}\mathrm{400}}\right)\right]\end{array}\right\}}{p_e+\left\{\begin{array}{c}\left[{\sum }_{k=\mathrm{1}}^{\mathrm{4}}{Z}_i({\mathit{\delta }}_s-{\mathit{\delta }}_k)\right]\\ \times \left[\mathrm{1}-\exp \left(-S\frac{R_{\text{dt}}}{R_{\text{fd}}}\frac{\mathrm{\Delta }t}{\mathrm{86}\mathrm{400}}\right)\right]\end{array}\right\}}}\right\},\end{array}\end{array}$$

  where h (m) represents the change in surface water depth, k is an integer representing the soil layer (ranging from 1 to 4), δ_k (m³ m⁻³) and Zk (m) are the soil moisture grid and depth of the kth soil layer, δ_s (m³ m⁻³) represents the maximum soil moisture content, Δt (s) represents the model time step, S denotes the coefficient from the regulating runoff infiltration Richards' equation, and R_fd and R_dt denote the saturated hydraulic conductivity and the tunable coefficients for surface infiltration, respectively.

• Groundwater flow. The baseflow module typically describes groundwater flow. The bucket model uses a conceptual storage equation. This equation depicts groundwater storage changes over time, influenced by input recharge and outflow proportional to the difference between current storage and baseflow threshold,

  $$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{d}S}{\mathrm{d}t}}=R-K_{\mathrm{b}}\cdot (S-S_{\mathrm{0}}),\end{array}$$

  where $\frac{\mathrm{d}s}{\mathrm{d}t}$ is the rate of change of storage with respect to time, R is the recharge to the groundwater storage (from excess soil moisture), K_b is the baseflow recession coefficient, S is the storage of groundwater, and S₀ is the baseflow threshold storage.

The horizontal resolution of WRF-Hydro is specified by segmenting the inner domain of WRF into a grid spacing of 100 m in this study. The horizontal routing grids for catchments in the WRF-Hydro model are computed using the Muskingum–Cunge method, which handles channel routing with time-varying parameter estimates and neglects the backwater effect (Wang et al., 2022). In this study, the vertical routing process integrates the Noah-MP land surface model (LSM), which includes four soil layers (10, 30, 60, and 100 cm) spanning a 2 m soil column from top to bottom. Enabling the fully coupled option initiates the involvement of the hydrological module, disaggregation–aggregation module, and LSM components of WRF-Hydro when running WRF. It should be noted that the default Noah configurations in WRF-Hydro were employed rather than using site-specific settings. Also, the baseflow bucket model is switched off for simulation periods; the WRF-Hydro model primarily accumulates subsurface runoff and redistributes it to the channel, effectively increasing river flow (Xue et al., 2000). It should be noted that the findings presented in this study should be considered a benchmark for the WRF-Hydro fundamental model performance. The intention is to offer valuable insights for future users of the model operating in particular basins within northern China and comparable regions, as well as to provide guidance for prospective model enhancements in future years.

3.3 HEC-HMS model

The Hydrologic Engineering Center's Hydrologic Modeling System (HEC-HMS) is a comprehensive software tool developed by the United States Army Corps of Engineers for hydrological modeling and the simulation of watershed runoff. It is used for a wide range of applications, including flood forecasting, reservoir management, and water resource planning. The HEC–HMS model version 4.10 is employed in this study (Bartles et al., 2006).

The key HEC-HMS model structure components are described below.

• Data processing. This component includes tools for GIS connection, data import, watershed delineation, and data transformation.

• Meteorologic data. This component allows input of historical or synthetic rainfall data in various formats.

• Hydrologic models.

  • The loss model estimates losses due to interception, depression storage, and infiltration. Direct runoff (Q) is calculated as

    $$\begin{array}{}\text{(5)}& \left(Q\right)={\displaystyle \frac{(P-I_{\mathrm{a}}{)}^{\mathrm{2}}}{(P-I_{\mathrm{a}}+S)}},\end{array}$$

    where Q is direct runoff, P is precipitation, I_a is initial abstraction, and S is potential maximum retention after runoff begins.

  • Routing model. The model routes flow through river channels using methods like the Muskingum–Cunge method:

    $$\begin{array}{}\text{(6)}& Q(k+\mathrm{1})=(\mathrm{1}-x)\times Q\left(k\right)+x(P-P_{\text{loss}}),\end{array}$$

    where Q(k+1) is the outflow at the next time step, Q(k) is the current outflow, x is the routing parameter, P is inflow, and P_loss represents losses (Niazkar and Zakwan, 2022).

• Unit hydrograph.

  • The ModClark unit hydrograph method uses a linear convolution equation:

    $$\begin{array}{}\text{(7)}& Q\left(t\right)={\displaystyle \frac{P\left(t\right)}{A}}\times \text{UH}\left(t\right),\end{array}$$

    where Q(t) is the discharge (runoff) at time t, P(t) is the precipitation at time t, A is the area of the watershed, and UH(t) is the ModClark unit hydrograph at time t (Che et al., 2014).

  • The SCS (Soil Conservation Service) unit hydrograph method uses a time–area approach:

    $$\begin{array}{}\text{(8)}& Q\left(t\right)={\displaystyle \frac{P\left(t\right)}{\mathrm{12}}}\times \text{UH}\left(t\right),\end{array}$$

    where Q(t) is the discharge (runoff) at time t, P(t) is the precipitation at time t, and UH(t) is the SCS unit hydrograph at time t (Shatnawi and Ibrahim, 2022).

3.4 Experimental design

In atmospheric–hydrological coupling, two-way coupling is essential for studying complex climate interactions, while one-way coupling is often employed for practical meteorological or climate prediction models to simplify computational demands and focus on specific phenomena (Wu et al., 2016). Therefore, this study adopts the one-way coupling when constructing the atmospheric–hydrological coupling systems. The experimental structure for the gridded coupled atmospheric–hydrological systems and the lumped system, which comprise three types of hydrological models, i.e., WRF, WRF-Hydro, and HEC-HMS models, is shown in Fig. 4. The WRF model is used to downscale the 1°×1° FNL data, which will generate gridded rainfall to drive the HEC-HMS and the WRF-Hydro model for flood simulation. The lumped HEC-HMS model is also used for flood simulation using the observed rainfall. The purpose is to set a benchmark for the coupled WRF–HEC-HMS system in order to analyze the source of model uncertainties in the coupled atmospheric–hydrologic system.

Figure 3 HEC-HMS schematic DEM map of Fuping and Zijingguan catchments.

Figure 4 Experimental structure for the coupled atmospheric–hydrological and lumped modeling system.

When calibrating the WRF-Hydro model, we employed manual calibration by carefully considering several crucial parameters that have the potential to impact flood-forecasting accuracy significantly. These parameters were identified through prior research on parameter sensitivity analysis conducted in the study region (Liu et al., 2021b). Because of the limited data in our study area, a systematic adjustment of model parameters to minimize discrepancies between observed and simulated streamflow was done for each storm event to enhance model accuracy. These parameters include the scaling parameter for overland flow roughness (OVROUGHRTFAC), the scaling parameter for surface retention depth (RETDEPRTFAC), the channel Manning roughness parameter (MannN), and the parameter for runoff infiltration (REFKDT).

For the calibration of the HEC-HMS model, the studied watersheds were delineated into 8 and 11 sub-basins for Fuping and Zijingguan, respectively, as shown in Fig. 3. The 1 km×1 km output gridded rainfall from the WRF was imported and interpolated using the bilinear resampling method. The ModClark and Soil Conservation Service (SCS) unit hydrograph methods were employed t simulate excess precipitation into direct surface runoff for the gridded and lumped models, respectively. The initial and constant method was employed to model infiltration loss, and the exponential recession model was used to model baseflow for the gridded and lumped HEC-HMS. Typically, in standard situations, the model calibration process involves making subjective adjustments to its parameters through a trial-and-error approach. While it is possible to calibrate the model manually, HEC-HMS additionally provides an inherent automatic optimization procedure designed to assess the suitability and practicality of parameter values and their respective ranges for the intended use of the model (Feldman and US Army Corps of Engineers, 2000).

Comparative analyses of the forecast flood processes were carried out to evaluate and compare the performance of the coupled and lumped systems. The flood-forecasting results are evaluated using the following statistical criteria: the Nash–Sutcliffe efficiency coefficient (NSE), root mean square error (RMSE), relative flood peak (R_f), and the relative flood volume (R_v).

$$\begin{array}{}\text{(9)}& {\displaystyle }\text{NSE}=\mathrm{1}-{\displaystyle \frac{{\sum }_{i=\mathrm{1}}^N(y_i^{\prime }-y_i{)}^{\mathrm{2}}}{{\sum }_{i=\mathrm{1}}^N(y_i-\stackrel{\mathrm{‾}}{y}{)}^{\mathrm{2}}}}\text{(10)}& {\displaystyle }\text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N(y_i^{\prime }-y_i{)}^{\mathrm{2}}}\text{(11)}& {\displaystyle }{R}_{\mathrm{f}}={\displaystyle \frac{(y_{\mathrm{f}}^{\prime }-y_{\mathrm{f}})}{y_{\mathrm{f}}}}\text{(12)}& {\displaystyle }{R}_{\mathrm{v}}={\displaystyle \frac{(y_{\mathrm{v}}^{\prime }-y_{\mathrm{v}})}{y_{\mathrm{v}}}}\end{array}$$

Here, $y_i^{\prime }$ and y_i represents the simulated and observed discharge flow at a specific time step denoted by i; N denotes the total number of flood event time steps; $\stackrel{\mathrm{‾}}{y}$ is the calculated average discharge; $y_{\mathrm{f}}^{\prime }$ and y_f represent the simulated and observed flood peaks, respectively; and $y_{\mathrm{v}}^{\prime }$ and y_v represents the simulated and observed flood volumes, respectively.

4 Results

4.1 Results from the coupled WRF–HEC-HMS and WRF–WRF-Hydro systems

The simulation results of the coupled WRF–HEC-HMS and WRF–WRF-Hydro modeling systems are shown in Fig. 5 and Table 5. As demonstrated in Sect. 3.4, the 1 km×1 km output rainfall from the WRF model is used to simulate the four storm events using the coupled atmospheric–hydrological systems. Comparing the simulation results, it can be seen that the coupled WRF–HEC-HMS model performs better for storm events 1 and 3 with NSE of 0.84 and 0.79, respectively, as these events have longer durations, demonstrating that the model adapted well in modeling prolonged floods. Also, the coupled WRF–HEC-HMS system's best performance is event 1 (NSE=0.84), characterized by long duration and relatively small flood magnitude, demonstrating the ability to model floods subjected to rapid recession. The coupled WRF–WRF-Hydro system performs better for events 2 and 4 with NSEs of 0.63 and 0.62, respectively, which are floods with very high magnitudes and shorter durations, demonstrating abilities in modeling flash floods.

Figure 5 Simulated flood hydrographs of the coupled WRF–HEC-HMS and WRF–WRF-Hydro systems for the four storm events.

Table 5 Simulation results of the coupled WRF–HEC-HMS and WRF–WRF-Hydro system for the four storm events.

By comparing the observed and simulated hydrographs of the coupled WRF–HEC-HMS and WRF–WRF-Hydro modeling systems (Fig. 5), it can be seen that the flow peaks were underestimated for both models. The WRF–HEC-HMS system has the largest flow peak error, ranging from −10.78 % to −58.48 %, but a better RMSE, ranging from 4.26 to 404.55. The worst flow peak error (−58.48 %) was for event 4 (highest flow peak and short duration), and it only has a better flow peak error for storm event 3 (−10.73 %) when compared to the coupled WRF–WRF-Hydro system. The WRF–WRF-Hydro system has a better flow peak error ranging from −5.75 % to −54.90 % but higher RMSE ranging from 7.11 to 513.59. Its worst flow peak error was also for event 4 (−54.9 %), for which it had a better performance. However, the coupled WRF–HEC-HMS model with an NSE value ranging from 0.51 to 0.84 and the coupled WRF–WRF-Hydro model with an NSE value ranging from 0.34 to 0.63 indicate that the HEC-HMS model results in a better average performance when coupled with WRF.

4.2 Results from the lumped HEC-HMS model driven by the observed rainfall

Considering the unsatisfactory performance of the coupling systems, the four storm events were also simulated with the lumped HEC-HMS using the observed rainfall (demonstrated in Sect. 3.4). Lumped models do not possess the necessary spatial data to depict hydrological processes accurately because they only have temporal inputs. The simulation results of the lumped HEC-HMS are shown in Table 6 and Fig. 6.

Table 6 Simulation results of the lumped HEC-HMS model driven by observed rainfall for the four storm events.

Figure 6 Flood-forecasting hydrographs of the coupled WRF–HEC-HMS and lumped HEC-HMS systems for the four storm events.

As shown in Table 6, the best performance of the lumped HEC-HMS is event 1 with NSE=0.90, which has the best uniform temporal distribution with a long duration. The worst performance is found in event 3 with NSE=0.48, which has the least uniform temporal distribution. Generally, the lumped HEC-HMS obtained a late flood peak for some storm events, i.e., −9 and −1 h for events 1 and 3, respectively. Also, it underestimated the flow peak for events 2 and 4, indicating some uncertainty in the model. The lumped HEC-HMS model shows an average performance that is not too great compared to the coupled WRF–HEC-HMS system, i.e., 0.70 (average lumped HEC-HMS NSE value) minus 0.68 (average coupled WRF–HEC-HMS value) equals 0.02.

The simulation results (R_f, RMSE, NSE) of the coupled WRF–HEC-HMS and lumped HEC-HMS model were compared (R_p–R_p lumped, RMSE–RMSE lumped, NSE–NSE lumped) as shown in Fig. 6 and Table 7. This comparison analyzes the level of uncertainty in the hydrological model. The largest error in performance was found in event 3 (NSE=0.31), which has the largest non-uniform temporal distribution, showing a better performance for the coupled WRF–HEC-HMS system. For the other three events, the difference between the coupled WRF–HEC-HMS model and the lumped HEC-HMS model increased as the rainfall temporal distribution heterogeneity increased. This shows that when the rainfall temporal distribution is uniform, as in event 1, event 2, and event 4, the performance error decreases to −0.07, −0.08, and −0.22, respectively. Also, there is a large difference in flow peak error and RMSE, indicating the level of underestimation of the flow peak from the coupled WRF–HEC-HMS system, more notably in storm event 4. These comparisons indicate some uncertainty from the coupled WRF–HEC-HMS model.

Table 7 Comparison between the coupled WRF–HEC and lumped HEC-HMS model for flood forecasting of the four storm events.

4.3 Error in the simulated WRF rainfall

The simulated WRF rainfall is analyzed to determine the error level, which influences the forecasting results of the coupled systems. The relative errors (R_v) of the 24 h rainfall accumulations of the observed and simulated four storm events are shown in Table 8. Figure 7 displays the temporal variations of the observed and simulated rainfall in the accumulative curves and time series bars for the four storm events.

Table 8 The four 24 h storm events observed and simulated WRF rainfall.

Figure 7 The 24 h accumulative curves of the observed and the simulated rainfall of the four storm events.

It can be seen in Table 8 and Fig. 7 that event 1, which has the most uniform temporal distribution with C_v=0.60 (as shown in Table 2), has a better simulation result compared to the other storm events, having the lowest R_v value of 9.53 %. The largest relative errors were found in storm events 3 and 4, with the most non-uniform temporal rainfall distributions. Generally, both coupled systems performed poorly, with the largest flow peak error for storm event 4 having the largest relative error of −66.73 %. Therefore, it can be seen that the simulation uncertainty of the storm events is directly proportional to the temporal rainfall non-uniformity, especially for events with very high peaks.

The spatial variation of the 24 h accumulations of rainfall of the four storm events is further analyzed in the study catchments. Figure 8 shows the spatial patterns of the accumulation rainfall distribution from the observed rain gauges, simulated WRF output, the coupled WRF–WRF-Hydro system, and the spatial differences of the accumulation rainfall distribution of simulated WRF output and the coupled WRF–WRF-Hydro model (i.e., the coupled WRF–WRF-Hydro system minus WRF), respectively. As seen in the figures, the largest spatial variation was found in event 3, with the most non-uniform spatial distribution C_v=0.74 (Table 2) compared to the other events. The simulations of the coupled WRF–WRF-Hydro system exhibited more noticeable variations for event 3 compared to the WRF model. This spatial difference is evident in subfigure (c) for event 3 in Fig. 8, illustrating a better performance of the coupled WRF–WRF-Hydro model compared to the WRF model at the crucial storm center. This improvement is attributed to the model's efficient rainfall spatial redistribution.

Figure 8 The 24 h accumulative spatial distributions of the four storm events in the study catchments: gauge observation, (a) WRF model, (b) coupled WRF–WRF-Hydro system, and (c) coupled WRF–WRF-Hydro system minus WRF.

Storm event 4, characterized by intense convective rainfall, experienced a rapid and substantial increase in rainfall intensity within a shorter duration, with maximum gauge cumulative values reaching 355 mm. Both the WRF model and the coupled WRF–WRF-Hydro system underestimated this storm, with the WRF model exhibiting even poorer performance. This suggests a failure to capture such intense storms in the Zijingguan catchment accurately. For storm events 1 and 2, the WRF and the coupled WRF–WRF-Hydro systems both exhibited nearly identical spatial patterns in cumulative rainfall distributions. The simulations for these events aligned more closely with cumulative rain gauge observations compared to events 3 and 4. The WRF and coupled WRF–WRF-Hydro models effectively captured the storm centers of events 1 and 2. However, for storm event 3, some areas with high rainfall accumulations within the catchment were not accurately represented by the WRF and coupled WRF–WRF-Hydro models compared to observations.

5 Discussion

The WRF atmospheric model, when coupled with hydrological models, significantly enhances the representation of precipitation and its impact on hydrological processes. The model's ability to account for localized variations improves accuracy in capturing spatial heterogeneity, thereby contributing to its effectiveness in various scenarios (Liu et al., 2021b). The results of the coupled WRF–HEC-HMS and WRF–WRF-Hydro systems highlight intriguing findings in flood simulation within mountainous catchments of northern China. These coupled models exhibit distinct characteristics in effectively capturing the complex interplay between precipitation and runoff. The coupled WRF–WRF-Hydro system, a fully distributed system, performs better for extreme storm events, as in the case of storm events 2 and 4. However, this coupled system performed poorly for storm events that exhibited faster surface runoff recession, exemplified in storm events 1 and 3. This may be because of the interaction between the land surface and rainfall–runoff generation, occurring within the WRF-Hydro model at a shorter integration time step. This increases the volume of infiltrated precipitation, leading to higher soil moisture and a decrease in runoff production (Branch et al., 2016). On the other hand, the coupled WRF–HEC-HMS model, a simpler semi-distributed model,performs better for prolonged floods, as in the case of storm events 1 and 3, as it can adapt to the rapid recession process. Nonetheless, the coupled WRF–HEC-HMS model may have oversimplified the hydrological processes in some storm events, resulting in less accurate predictions, particularly in extreme conditions such as events 2 and 4. The lumped parameter approach in HEC-HMS, i.e., assuming uniformity in parameters such as soil properties and land use, might not be able to represent the rapid hydrological changes across different sub-basins.

The performance error of the lumped system reflects the uncertainty inherent in the hydrological model. The lumped HEC-HMS is driven by observed data, representing a system with the true inputs of the studied catchments. Adopting the lumped HEC-HMS model allows us to assess the influence of rainfall errors on the coupled systems. The comparison between the lumped HEC-HMS model driven by observed rainfall and the coupled WRF–HEC-HMS model reveals that the error is primarily attributed to simulated WRF rainfall. The model structure also has some influence (lumped > gridded), although this influence is negligible compared to the simulated rainfall error. The results show minimal differences in average model performance between the lumped HEC-HMS and the coupled WRF–HEC-HMS model. The simulation results of the coupled WRF–HEC-HMS and WRF–WRF-Hydro modeling systems are influenced by the uncertainties in the hydrological models and the driven data. Previous studies have shown that the error in the WRF output rainfall, which depends on the quality of the driven data (i.e., FNL data), necessarily causes parallel uncertainty in hydrological forecasts (Merino et al., 2022). We analyzed the temporal and spatial errors inherent in the simulated WRF output rainfall, influencing the simulation results of the coupled systems. The temporal variations scrutinize how well the WRF model captures the timing and intensity of rainfall events, providing a comprehensive understanding of its predictive capabilities in the study catchments. The most significant temporal variation between the 24 h accumulations of observed and simulated rainfall was found in storm event 4, characterized by the highest rainfall intensity. Additionally, we assess the 24 h accumulations of rainfall spatial variation, identifying regions where the models may struggle to predict precipitation patterns accurately in the study catchments. The most significant spatial variation was found in event 3, where both the WRF model and coupled WRF–WRF-Hydro systems failed to capture areas with high rainfall accumulations within the catchment accurately compared to the observations.

Even though there are some uncertainties in the coupled atmospheric–hydrological systems, it is evident that a larger magnitude of error is attributed to the WRF output rainfall. This indicates that the primary factor influencing the overall accuracy of the coupled systems is the accuracy of simulated WRF rainfall. To enhance simulation rainfall in small- and medium-scale mountainous catchments, such as those in northern China, here are some recommendations: the observed rainfall can be used to correct the simulated rainfall, implement radar data assimilation in numerical weather prediction (NWP) models, and integrate the simulated rainfall with radar quantitative precipitation forecasts (QPFs) or quantitative precipitation estimates (QPEs; Vendrasco et al., 2016; Tong et al., 2016; Liu et al., 2021a). In comparison with previous studies, our results show significant improvements in flood-forecasting accuracy by coupling the WRF model with different hydrological systems. This aligns with Jasper et al. (2002) on NWP model development, Bartholmes and Todini (2005) on high-resolution NWP data benefits, and Cattoën et al. (2016) on forecast reliability. Our findings also support Li et al. (2017) on decreased accuracy with longer lead times, mitigated here by high-resolution models. By integrating fully distributed and semi-distributed models, we extend the work of Ming et al. (2020) and Chen et al. (2020), who demonstrated the advantages of coupled NWP and hydrodynamic models. Similar to the research of Giannaros et al. (2021) and Varlas et al. (2024), our approach shows that high-resolution NWP data integrated with hydrological models can enhance real-time flood prediction. Implementing these coupled systems in operational frameworks can provide more accurate and timely flood warnings.

It should be noted that the unique topographic and climatic features of our study area compound the inherent uncertainty in simulations. Coupling the WRF model with hydrological models introduces additional challenges, such as model parameterization, spatial static data, downscaling resolution, and integration time step, all of which can significantly influence simulation outcomes. Furthermore, uncertainties in meteorological inputs, such as precipitation forecasts, and the choice of model coupling approaches, such as the one-way coupling used in this study, may affect the reliability of flood predictions. However, the general conclusions of this study aim to provide valuable insights into the performance and potential enhancements of this modeling approach in the face of complex topographical and meteorological conditions.

6 Conclusions

This study coupled the semi-distributed HEC-HMS and fully distributed WRF-Hydro models with a 1 km×1 km rainfall output from the WRF model, alongside the lumped HEC-HMS using observed gauge precipitation. We conduct a comparative analysis of the forecast processes from these coupled hydrological systems for storm events with varying characteristics. The lumped HEC-HMS model is adopted as a benchmark to compare and analyze the level of uncertainty in the coupled WRF–HEC-HMS model system. Additionally, we analyze the error in simulated WRF output rainfall. From the results, we conclude the following.

• The coupled WRF–HEC-HMS system exhibits better performance in predicting prolonged storm events with optimal accuracy in uniformly distributed spatial and temporal patterns (e.g., event 1). It effectively adapts to rapid recession processes, despite challenges in accurately capturing flood magnitudes, leading to larger flow peak errors.

• In contrast, the coupled WRF–WRF-Hydro system performs better for shorter-duration floods characterized by higher flow peaks, demonstrating strong capability in flash flood forecasting. However, its performance is reduced as uniformity in storm events decreases, which might be due to incorrect representation of the spatial rainfall.

• The performance of the lumped model driven by the gauge-observed rainfall indicates the uncertainty in the hydrological models when compared with the coupled WRF–HEC-HMS system, but a larger magnitude error was found in the WRF output simulated rainfall.

The conclusions from this study verify the notion that the coupled atmospheric–hydrological modeling systems are influenced by the rainfall simulation accuracy and complexity of the hydrological model. We hope this study will encourage further research on improving the simulated rainfall and hydrological models to verify the conclusions of this study.