**Research article**
11 Aug 2020

**Research article** | 11 Aug 2020

# A coupled atmospheric–hydrologic modeling system with variable grid sizes for rainfall–runoff simulation in semi-humid and semi-arid watersheds: how does the coupling scale affects the results?

Jiyang Tian Jia Liu Yang Wang Wei Wang Chuanzhe Li and Chunqi Hu

^{1},

^{1},

^{1,2},

^{1,3},

^{1},

^{4}

**Jiyang Tian et al.**Jiyang Tian Jia Liu Yang Wang Wei Wang Chuanzhe Li and Chunqi Hu

^{1},

^{1},

^{1,2},

^{1,3},

^{1},

^{4}

^{1}State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing, 100038, China^{2}Swanson School of Engineering, University of Pittsburgh, Pittsburgh, PA 15261, USA^{3}College of Hydrology and Water Resources, Hohai University, Nanjing, 210098, China^{4}Bureau of Water Resources Survey of Hebei, Shijiazhuang, 050031, China

^{1}State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing, 100038, China^{2}Swanson School of Engineering, University of Pittsburgh, Pittsburgh, PA 15261, USA^{3}College of Hydrology and Water Resources, Hohai University, Nanjing, 210098, China^{4}Bureau of Water Resources Survey of Hebei, Shijiazhuang, 050031, China

**Correspondence**: Jia Liu (hettyliu@126.com)

**Correspondence**: Jia Liu (hettyliu@126.com)

Received: 01 Nov 2019 – Discussion started: 20 Jan 2020 – Revised: 15 May 2020 – Accepted: 09 Jul 2020 – Published: 11 Aug 2020

The coupled atmospheric–hydrologic modeling system is an effective way to improve the accuracy of rainfall–runoff modeling and extend the lead time in real-time flood forecasting. The aim of this study is to explore the appropriate coupling scale of the coupled atmospheric–hydrologic modeling system, which is established by the Weather Research and Forecasting (WRF) model and the gridded Hebei model with three different sizes (1 km×1 km, 3 km×3 km and 9 km×9 km). The Hebei model is a conceptual rainfall–runoff model designed to describe a mixed runoff generation mechanism, including both storage excess and infiltration excess, in the semi-humid and semi-dry area of northern China. The soil moisture storage capacity and infiltration capacity of different grids in the gridded Hebei model are obtained and dispersed using the topographic index. The lumped Hebei model is also used to establish the lumped atmospheric–hydrologic coupled system as a reference system. Four 24 h storm events occurring at two small- and medium-scale sub-watersheds in northern China are selected as case studies. Contrastive analyses of the flood process simulations from the gridded and lumped systems are carried out. The results show that the flood simulation results may not always be improved with higher-dimension precision and more complicated system, and the grid size selection has a strong relationship with the rainfall evenness. For the storm events with uniform spatial distribution, the coupling scale has less impact on flood simulation results, and the lumped system also performs well. For the storm events with uneven spatiotemporal distribution, the corrected rainfall can improve the simulation results significantly, and higher resolution leads to better flood process simulation. The results can help to establish the appropriate coupled atmospheric–hydrologic modeling system to improve the flood forecasting accuracy.

Compared to the traditional flood forecast models that take gauge precipitation as inputs, the coupling of a mesoscale numerical weather prediction (NWP) model with a hydrological model has been proven to be an effective way to improve the accuracy of rainfall–runoff modeling and extend the lead time in real-time flood forecasting (Wu et al., 2014; Wan and Xu, 2011). The atmospheric–hydrologic coupled system has become the current world technology frontier for flood forecasting, and correlation studies have received extensive attention (Benoit et al., 2000; Tian et al., 2019). Given the multiscale nature of watershed hydrologic models and rainfall data output by numerical atmospheric models, one key to obtaining accurate simulations results is the establishment of a suitable atmospheric–hydrologic coupled system.

Most atmospheric–hydrologic coupling studies tend to use lumped or
distributed hydrological models (constructed in research watersheds) to run
atmospheric models based on hydrological model requirements by inputting
meteorological information such as rainfall (Wagner et al., 2016). Liu et al. (2015) simulated the floods by a lumped hydrological model integrated with the Weather Research And Forecasting (WRF) model at a 10 km scale at a small catchment (135.2 km^{2}). Rogelis et al. (2018) focused on comparison of streamflows obtained from three different lumped hydrological models driven by WRF precipitation forecast at the Tunjuelo River basin (380 km^{2}). The spatial resolution of 1.67 km was set for finest domain in the WRF model. Li et al. (2017) extended the flood forecasting lead time in the Liujiang River basin (58 270 km^{2}) by WRF precipitation forecast with a distributed hydrological model. The WRF precipitation forecast with the resolution being 20 km × 20 km was downscaled to 200 m × 200 m, which was suited for the resolution of the hydrological model. However, most studies have no consideration for the spatial-scale matching issue between atmospheric and hydrological models, which has significant impact on the flood forecast accuracy.

Model complexity is another impact factor that must be considered during model construction. When investigating the complexity of data-driven models, Myung et al. (2009) found that the relationship between model complexity and error can be described as in Fig. 1. It is outlined that model selection and construction should not solely be based on the goodness of fit but must also consider the model complexity (whether the model's description of observed data is achieved in the simplest possible manner). This is because a highly complex model can provide a good fit without necessarily bearing any interpretable relationship with the underlying process. It is shown that model selection based solely on the fit to observed data will result in the choice of an unnecessarily complex model that overfits the data and, thus, generalizes poorly (Myung, 2000).

This could also be true when dealing with the coupled atmospheric–hydrologic systems. Although the models in the coupled system are physically based with deterministic structures and parameters, the complexity is to a great extent decided by choosing an appropriate spatial scale (grid size) for coupling (Verbunt et al., 2006; Hostache et al., 2011; Rogelis and Werner, 2018). A smaller coupling scale between the atmospheric and the hydrologic models helps utilize higher-resolution rainfall information, but a higher grid resolution results in a greater computational load, leading in turn to greater cumulative errors in the simulated discharge at the watershed outlet. On the other hand, a larger coupling scale allows a smaller computational load, but the accuracy of the rainfall information decreases such that the final simulation results may be subject to larger errors. Therefore, finding the underlying law of how the performance of the coupled atmospheric–hydrologic system is impacted by the coupling scale is of great importance in enhancing the accuracy of rainfall–runoff simulation. Within this scope, important issues that should be addressed in depth include how to effectively use the high-resolution rainfall information provided by atmospheric models, how to select a suitable spatial scale for constructing a coupled atmospheric–hydrologic system, and how to make tradeoffs between the computational efficiency and accuracy. Moreover, to date, the spatial scale of atmospheric–hydrologic coupling has not yet been investigated with respect to different rainfall processes.

In this study, a coupled atmospheric–hydrologic system with variable grid scales is firstly developed based on the construction of a gridded hydrologic model in combination with a mesoscale NWP model to provide rainfall driving data at different resolutions; this allows the atmospheric–hydrologic coupling at different spatial scales. The most recently developed mesoscale NWP model, i.e., the Weather Research & Forecasting (WRF) model, is adopted to provide the downscaled rainfall information and the Hebei model (Tian et al., 2019), a conceptual model with mixed runoff generation mechanisms of both saturation-excess and infiltration-excess, is used to construct the gridded hydrologic model. Three different coupling scales are considered in this study, which is 1 km×1 km, 3 km×3 km, and 9 km×9 km. Next, the coupled atmospheric–hydrologic system with different spatial scales to reflect the degree of model complexity, is applied to four storm events in semi-humid and semi-arid watersheds of northern China. The storm events are characterized by different representatives of the rainfall distribution evenness in space and time. The performance of the coupled system in different degree of complexity (i.e., with different coupling scale or grid size) is fully investigated with regard to different storm events, and the relationship between the rainfall evenness and the selection of the most appropriate coupling scale is further discussed.

## 2.1 Study site

In this study, two mountainous sub-watersheds of the Daqing River basin (Fuping of the south branch and Zijingguan of the north branch) are chosen
as the study area (Figs. 2 and 3). The Fuping sub-watershed has a total area of 2219 km^{2} and is located in the upper reaches of the Shahe River, a south branch of the Daqing River. The Zijingguan sub-watershed has a total area of 1760 km^{2} and is located in the upper reaches of the Juma River, a north branch of the Daqing River. The two sub-watersheds are the most concentrated and typical cinnamon regions. The land use mainly includes grassland, farmland, and forestland. The soil erosion is severe. Due to the dry soil conditions and the groundwater overexploitation, the river has extensive seepage during the storm season. More information about the two sub-watersheds is shown in Table 1. The study area embodies the representative rainfall–runoff characteristics of the sub-humid and sub-arid area in northern China. Rainfall in northern China is characterized by
summer storms with short durations and large intensities, which are likely
to result in severe flood disasters, especially in mountainous areas like Fuping and Zijingguan.

## 2.2 Storm events

Based on the analyses of the historical floods in the study area, four representative 24 h duration storm events with relatively high flow peaks (Fig. 4) are selected and used to test the performance of the coupled atmospheric–hydrologic system for rainfall–runoff modeling. The start and end times together with the 24 h cumulative rainfall amounts and the peak discharges are listed in Table 2. Among the four events, three occur in Fuping sub-watershed and one occurs in Zijingguan sub-watershed. Before the establishment of the coupled atmospheric–hydrologic system, quality controls of the observational rainfall–runoff data are carried out. Rainfall observations from the rain gauges are verified by the weather radar, and interpolations are done to guarantee the continuity of the flow observations (Liu et al., 2018).

The coefficient of variation *C*_{v} is used as an indicator to verify the representativeness of the selected storm events and to characterize their
temporal and spatial distributions. With regard to the spatial distribution,
the coefficient reflects the variations in the 24 h cumulative rainfall among different rainfall stations across the watershed. As for the temporal
distribution, it shows the variations in the average areal rainfall at
different time steps during the whole storm event. The coefficient of
variation *C*_{v} can thus be calculated as (Yue et al., 2014; Jha et al., 2015; Tian et al., 2017a, b)

For spatial distribution calculations, *x*_{i} is the 24 h cumulative rainfall at a certain rainfall station *i*, $\stackrel{\mathrm{\u203e}}{x}$ represents the average cumulative rainfall of all stations, and *n* refers to the number of rainfall stations in the sub-watershed. For temporal distribution calculations, *x*_{i} is the areal rainfall at the *i*th time step, $\stackrel{\mathrm{\u203e}}{x}$ is the average of the areal rainfall of all time steps, and *n* refers to total time steps the rainfall lasts for, which is 24 h in this study. The calculated spatial and temporal values for the four selected events can be found in Table 3. The smaller the value of *C*_{v} is, the more even the rainfall distribution is in space or time. According to Table 3, the ranking of the distribution evenness of rainfall in space is event 2 > event 1 > event 4 > event 3, and that in time is event 1 > event 2 > event 4 > event 3. Event 1 has relatively even rainfall in both space and time, and the rainfall of event 2 is evenly distributed in space but unevenly distributed in time. Both event 3 and event 4 have unevenly distributed rainfall in space and time, while event 4, which can be noted with much larger *C*_{v} values and higher cumulative rainfall, is considered to be an extreme case. The storm happened
in a very short period of time and in a very concentrated area with high
rainfall intensity.

## 3.1 Atmospheric modeling

### 3.1.1 WRF model setup

The WRF model with a variety of physical and numerical options has widely
been used for numerical rainfall prediction in the catchment scale (Liu et
al., 2012; Cassola et al., 2015; Tian et al., 2017b). In this study, three
nested domains are adopted, centered over the Fuping sub-watershed (at
39^{∘}04^{′}15^{′′} N, 113^{∘}59^{′}26^{′′} E) or the Zijingguan sub-watershed (39^{∘}25^{′}59^{′′} N, 114^{∘}46^{′}01^{′′} E). The two sub-watersheds are totally covered by the innermost domain, with the interaction between the parent and child domain being allowed. The horizontal grid spacing is set to 1, 3, and 9 km from the innermost to the outermost domain, which is consistent with the resolution of the hydrological model and allows easy operation during coupling (Tian et al., 2017a). The NCEP/NCAR final operational global analysis (FNL) data with spatial resolution of 1^{∘} × 1^{∘} and temporal resolution of 1 h are used to provide the lateral and boundary conditions of the WRF model (C. C. Wang et al., 2013; H. Wang et al., 2013). In order to eliminate the discrepancy of the initial and boundary conditions with the driven data, another outermost domain is set beyond the WRF three nested domains to downscale the FNL data to a spatial resolution of 27 km. The integration time step of the model is set to 6 s, and the output data interval is set to 1 h (Hong and Lee, 2009). Forty layers are considered in the three nested domains, with a top-layer pressure of 50 hPa. The projection method is chosen to be Lambert, which has good applicability in the midlatitude regions. The other basic settings of the WRF model adopted in this study can be found in Table 4.

### 3.1.2 Physical parameterizations

The WRF rainfall simulation is sensitive to the selection and the combination of its physical parameterizations (Di et al., 2015). In this study, the well-performing and extensively used parameterizations in northern China were chosen (Miao et al., 2011; Di et al., 2015; Tian et al., 2017a), which include two microphysics parameterizations, i.e., Purdue-Lin (Lin) (Lin et al., 1983) and WRF Single-Moment 6 (WSM6) (Hong et al., 2006); two cumulus parameterizations, i.e., Kain–Fritsch (KF) (Kain, 2004) and Grell–Devenyi (GD) (Grell and Freitas, 2014); and two PBL (planetary boundary layer) parameterizations, i.e., Mellor–Yamada–Janjic (MYJ) (Hong et al., 2006) and Yonsei University (YSU) (Janjić, 1994). Besides, Rapid Radiative Transfer Model (RRTM) and Dudhia (Evans et al., 2012) usually cooperate well as the long- and shortwave radiation parameterizations, and Noah is chosen to be the land surface model (Chen et al., 2014). The most suitable physical parameterizations resulting the best rainfall simulations (Tian et al., 2017a) are used for each of the four storm events, as shown in Table 5. It should be clarified that the comparison of different coupling scales is carried out for each of storm event under the same model parameterizations (MPs) and thus using different MPs for different events will not cause difficulties in analyzing the final results.

## 3.2 Hydrologic modeling

### 3.2.1 Construction of a grid-based hydrologic model

A mesoscale NWP model can produce downscaled rainfall data with different horizontal resolutions through the multilayer nested grids. In order to use the spatial information of the WRF rainfall output data, a loosely structured hydrological model is constructed in each of the two sub-watersheds. The sub-watersheds are firstly divided into grids, based on which a rainfall–runoff model is established on each grid cell to calculate the runoff generation and concentration within the grid cell. After that, a spatially distributed river network model is used to calculate the confluence of the river flow from the grids to obtain the total discharge at the outlet of the watershed. As the scale of grid-based hydrologic models is dependent on the grid size, it is possible to conduct the grid division of the watershed according to the horizontal resolution of the WRF output rainfall data.

In this study, the 30 m digital elevation model (DEM) provided by Geospatial Data Cloud (http://www.gscloud.cn/, last access: 30 January 2019) assists in establishing the grid-based hydrologic model. Based on the DEMs of the two sub-watersheds and the GIS spatial analysis tools, the flow direction of each grid cell is determined, the river network is generated, and the computing orders of the grids are then obtained. Considering the horizontal resolution of the WRF model in this study is set to be 1, 3, and 9 km from the innermost to the outermost domain, the two sub-watersheds, Fuping and Zijingguan, are divided into grids with the same three resolutions. The divided 1 km×1 km, 3 km×3 km, and 9 km×9 km grids of the two sub-watersheds with a distinction between channel and non-channel cells are shown in Figs. 5 and 6.

### 3.2.2 Runoff generation and water exchange processes

In this study, a conceptual rainfall–runoff model, named Hebei model, is built on each grid cell of the two sub-watersheds. The Hebei model is specially developed for rainfall–runoff modeling in the semi-humid and
semi-dry area of northern China and has widely been applied in Hebei Province by considering both infiltration-excess and saturation-excess
mechanisms of the runoff generation (Tian et al., 2019). Due to the
perennial water shortage and groundwater overexploitation, both storage excess and infiltration excess are found in the study area with extensive
seepage along the river channel during the storm season. The obvious
advantage of the Hebei model is the consideration of both storage-excess and
infiltration-excess mechanisms for rainfall–runoff generation. The model is
easily applied and can be used in other semi-humid and semi-arid watersheds.
In the Hebei model, the description for the storage excess part is the same
as that in the Xin'anjiang model (Zhao, 1992). On the other hand, the
infiltration capacity across the watershed is described by a distribution
curve described below, and the Horton model is applied to calculate the
seepage along the river channel during the river routing (Horton, 1933). The
structure of the Hebei rainfall–runoff model is shown in Fig. 7 and
calibrated parameters are listed in Table 6. *γ* is the proportion of the area with the infiltration capability within a certain value, *χ* is
the proportion of no runoff generation area to entire basin area, *W*^{′} is the storage capacity of the point in the basin, and WMM is the storage capacity of the maximum point in the basin. The surface runoff is generated when the rainfall intensity exceeds the infiltration capacity, and when the infiltrated water volume matches the soil's water shortage capacity, the surplus water generates the groundwater runoff. The total flow at the watershed outlet is composed of both the surface and the underground runoff.

The infiltration curve of the model can be expressed as

where *f* is infiltration rate (mm h^{−1}), *i* is rainfall intensity (mm h^{−1}), *n* is an exponent, *f*_{c} is the stable infiltration rate (mm h^{−1}), *m* is the surface soil moisture (mm), *u* is an exponent indicating the decreasing speed of the infiltration rate with the increase in the soil moisture, and *f*_{m} is the infiltration capacity inside the grid (mm h^{−1}).

When the calculation period is *t* hours, the infiltration volume in that time period can be calculated as

where $\stackrel{\mathrm{\u203e}}{{f}_{\mathrm{i}}}$ is the accumulation infiltration volume in *i*th time period.

Thus, surface runoff *R*_{s} in a given time period can be expressed as

where *P*_{t} is the precipitation volume (mm), *E*_{t} is the evaporation
volume (mm), and *F*_{t} is the infiltration volume (mm) in the time period *t*. The calculation of the underground runoff (*R*_{g}) is the same with the Xin'anjiang model. Detailed descriptions are shown by Zhao (1992).

According to the principles of infiltration volume calculation in the time period and the spatial distribution of the water storage capacity in the grid, the groundwater runoff in the time period can be calculated as follows.

When ${P}_{\mathrm{a}}+{F}_{t}<{W}_{m}^{\prime}$,

When ${P}_{\mathrm{a}}^{\prime}+{F}_{t}\ge {W}_{\mathrm{m}}^{\prime}$,

where *R*_{g} is the groundwater runoff in the time period (mm), *E*_{t} is the evaporation volume in the time period, *P*_{a} is the influenced precipitation during the early phase of the grid (mm), *W*_{m} is the average storage capacity of the grid (mm), ${W}_{\mathrm{m}}^{\prime}$ is the storage capacity of the grid at the maximum point (in mm), and *b* is the exponent of the water storage capacity curve.

### (1) Exchange of water in a non-channel grid cell

The exchange of water between non-channel grid cells is illustrated in Fig. 8, where *P* and *E* are the rainfall and the evaporation in the grid cell, *Q*_{Si} and *Q*_{Gi} are the surface runoff inflow and the underground runoff inflow from the upstream grid cell, *Q*_{S} and *Q*_{G} are the surface runoff and the groundwater runoff, and *Q*_{So} and *Q*_{Go} are the surface and the groundwater runoff outflow to the downstream grid cell.

As mentioned, when the rainfall intensity in the grid cell is greater than the infiltration capacity, the surface runoff *Q*_{S} occurs in this cell. If the upstream cell also generates the surface runoff *Q*_{Si}, the surface runoff outflow *Q*_{So} can be calculated as

As part of the groundwater runoff from the upstream grid cell supplements
the soil water storage, the ratio of which is set as “sc”; the groundwater
runoff outflow *Q*_{Go} is calculated as

The soil water storage of the grid cell ${P}_{{\mathrm{a}}_{t}}$ at the end of the period *t* is

If the soil water storage ${P}_{{\mathrm{a}}_{t}}$ is greater than the storage capacity of the grid cell, then let ${P}_{{\mathrm{a}}_{t}}$ equal to the soil moisture capacity.

### (2) Exchange of water in a channel grid cell

For the exchange of water in a channel grid cell, the generated surface and groundwater runoff is no longer supplied directly to the downstream cell, but imported to the river channel and passes through the channel routing to the next cell. The processes are as illustrated in Fig. 9; the symbols of which have the same meanings as Fig. 8.

With the rainfall *P* and the evaporation *E*, the surface runoff *Q*_{S} and the groundwater runoff *Q*_{G} are generated. Based on the flow routing, the surface and the groundwater inflow to the river channel *Q*_{Sr} and *Q*_{Gr} can be obtained. The generated total inflow *Q*_{c} of the channel grid is then calculated as

where *ε* is the proportional coefficient of the surface and the
groundwater runoff into the river channel. By adding *Q*_{in} and *Q*_{c}, the outflow of the channel grid cell *Q*_{out} to the downstream cell can be obtained as

### 3.2.3 Confluence calculation and channel seepage

The storage–discharge relation of grid cell can be simplified to a single-line form. The storage–discharge equation can be obtained as

where *S* is water storage (m^{3}), *Q* is the outflow rate (m^{3} s^{−1}), *A* is the confluence parameter, and *ω* is the shape parameter. According to the water balance equation

the outflow rate of the grid cell *Q*_{t} at the end of the period *t*
(Goutal and Sainte-Marie, 2011) is

where *I* is the inflow rate (m^{3} s^{−1}), *I*_{t−τ} is the inflow rate (m^{3} s^{−1}) at time *t*−*τ* (*τ* is the travel time of flood wave), and *D* can be expressed as

where Δ*t* is the calculation time interval (*s*).

Due to the perennial water shortage and extensive channel seepage in study site, Horton infiltration model is applied to obtain the infiltration volume. When the calculation interval is 1 h, the infiltration volume can be calculated by the following equation (Horton, 1933):

where *f*_{0} is the initial infiltration rate (mm h^{−1}), and *k* is a parameter about soil. The channel seepage should be deducted before channel confluence is calculated with Eq. (14).

### 3.2.4 Discretization method for the soil moisture storage capacity and the infiltration capacity

In a grid-type hydrological model, it is necessary to determine the soil
moisture storage capacity and infiltration capacity in the grid cells in order to calculate the runoff generated within the grid and to simulate the
generation and convergence of runoff in the river basin. Therefore, the
determination of soil moisture storage capacity and infiltration capacity are key to model construction. Beven and Kirky (1979) proposed a topography-based semi-distributary hydrological model (TOPMODEL) that fully considered the effects of topography on the formation and change of runoff areas, using the spatial distribution of the terrain index ln (*α*∕tan *β*) to reflect the spatial distribution of saturated and deficit water volumes in the river basin. Based on this theory, it can be assumed that areas with similar topographic indices have the same hydrological response.

According to a statistical analysis of topographic indices for the 1 km×1 km, 3 km×3 km, and 9 km×9 km grids in the Fuping and Zijingguan basins, the cumulative distribution curves of the different grids' topographic indices in the same area have the same shape (parabolic). However, the cumulative distribution curve of topographic indices between different areas is also very similar. Experimentation showed that the soil moisture storage capacity and infiltration capacity of different grids can be obtained and dispersed using the topographic indices as follows:

where *W*_{i} is the moisture storage capacity of a certain grid cell (mm),
WMM is the maximal moisture storage capacity of a certain grid cell, TI_{i} is the topographic index of the grid, TI_{min} is the topographic index of the minimum point in the basin, *α* is the scale parameter of the grid, and *β* is the shape parameter of the grid.

where *f*_{i} is the infiltration capacity of a certain grid cell in the river basin, and *f*_{m} is the infiltration capacity of the maximum point in the basin.

## 3.3 Establishment of coupled atmospheric–hydrologic systems

Either one-way or two-way coupling can be achieved. For the latter, it is necessary to establish a communication mechanism between the atmospheric and hydrologic models to allow both to respond and combine dynamically. However, this process requires more complicated mechanisms and a higher amount of supportive data, making this approach difficult to use widely in actual forecasting, so this study focuses on one-way coupling. In order to research coupling scales, it is necessary to establish a gridded atmospheric–hydrologic coupled system. The WRF model is used to generate gridded rainfall data, and the gridded Hebei model is used as the land surface hydrologic model. The latter is set up with different grid divisions to allow the retrieval of rainfall data with a corresponding precision from the output data of the WRF model. Differences in the output processes of the hydrological models at different scales are then analyzed, and the relationships between the differences and rainfall patterns are considered. The lumped Hebei model is also used to establish the lumped atmospheric–hydrologic coupled system. The gridded rainfall data from the WRF model are averaged over each sub-watershed, which is regarded as the input of the lumped Hebei model. The coupled atmospheric–hydrologic system is illustrated in Fig. 10.

The gridded Hebei model is divided into the three different grid sizes noted above. The coordinates of the grid cell centers could be used to retrieve the corresponding output data from the WRF model for driving the hydrologic model. The SCE-UA (shuffled complex evolution) method (Duan et al., 1994) is used to calibrate the parameters, and the calibrated values are shown in Table 6. Due to the limited observational data, seven storm events in Fuping and six storm events in Zijingguan are selected and used to calibrate the Hebei model, and another two from each sub-watersheds are used for model validation. In order to guarantee reasonable values for the initial model conditions, the storm events are not independently used but rather with an antecedent period of data with the length of 15 days before the start of the event. The validation results show an average Nash–Sutcliffe efficiency coefficient (NSE) value of up to 0.686, indicating the calibrated models are reliable for further applications. It should be noted that the four storm events in Sect. 2.2 are different from those used for calibration and validation.

## 3.4 Evaluation statistics

The flood simulation results are evaluated using three statistics: the
relative flood peak error (*R*_{f}), the relative flood volume error (*R*_{l}), and the Nash–Sutcliffe efficiency coefficient (NSE).

where ${Q}_{\mathrm{f}}^{\prime}$ and *Q*_{f} are the simulated and observed flood peak flow, ${Q}_{\mathrm{l}}^{\prime}$ and *Q*_{l} are the simulated and the observed flood volume, ${Q}_{i}^{\prime}$ and *Q*_{i} are the simulated and observed flow discharge at the *i*th time step, *N* is the total time steps of the flood event, and $\stackrel{\mathrm{\u203e}}{Q}$ is the average value of *Q*_{i} at each time step.

## 4.1 Simulation results of the coupled atmospheric–hydrologic systems

The coupled atmospheric–hydrologic systems mentioned in Sect. 3.3 are used to simulate the four storm events. The flood processes at different grid divisions are retrieved by the gridded atmospheric–hydrologic coupled system, and the comparison is also made with the results from the lumped atmospheric–hydrologic coupled system. As shown in Fig. 11 and Table 7, the lumped system performs worse than the gridded system in most cases. The simulation results of the lumped system are close to the gridded system with different grid divisions for storm events 1 and 2, while the simulation results have a significant difference between lumped system and grid system for storm events 3 and 4. The gridded system, especially with the 1 km×1 km and 3 km×3 km grid, can improve accuracy and reduce error in the flood simulation for event 4.

The simulated results based on different grid sizes are similar for storm event 1, though the 9 km×9 km grid has the optimal simulation performance. All the coupled atmospheric–hydrologic systems at different grid divisions perform well for storm event 2, and the 3 km×3 km grid produces the optimal simulation result. For storm event 3 (a short-duration heavy rainfall), although the peak flow error and flood volume error are not significant, the simulation results fail to truly reflect the runoff process. The 1 km×1 km grid generates the optimal simulation result for storm event 4.

It can be easily found that different rainfall events have different
simulation results with the three grid sizes. The absolute value of
*R*_{f} is within 3.69 %–8.60 % for the system with 1 km resolution, 4.69 %–10.15 % for the system with 3 km resolution, and 6.09 %–9.08 % for the system with 9 km resolution. The absolute value of *R*_{l} ranges from 4.40 % to 11.81 % for the system with 1 km resolution, ranges from 9.12 % to 10.36 % for the system with 3 km resolution, and ranges from 8.62 % to 13.43 % for the system with 9 km resolution. The NSE is from 0.3669 to 0.8986 for a 1 km×1 km grid, from 0.3232 to 0.8573 for a 3 km×3 km grid, and from 0.3105 to 0.8302 for a 9 km×9 km grid. Considering the rainfall simulation errors from the WRF model for storm event 3, the 1 km×1 km and 3 km×3 km grids provide better simulation results overall, while the 9 km×9 km grid leads to unstable variations over a wide area. However, for a specific storm event, higher resolution may not lead to better flood process simulation. For example, the 9 km×9 km grid provides the optimal simulation results for storm event 1, while the best grid size is 3 km×3 km for rainfall event 2.

## 4.2 Relationship between the rainfall evenness and the performance of the gridded model with different grid size

Considering the characteristics of the spatial rainfall distributions, the
gridded rainfall simulation results from the WRF model are more suitable for the
application of the gridded Hebei model than lumped Hebei model. In order to
analyze the relationship between the rainfall evenness and the performance
of the gridded model with different grid size, the *C*_{v} values are
calculated for the rainfall simulations from the WRF model. Table 8 shows that
there is not much difference among the *C*_{v} values with different grid size for the same storm event, while the difference is obvious for different storm events. The ranking of the spatial evenness of the rainfall
simulations is event 1 > event 2 > event 3 > event 4.

The average flood simulation results (${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{l}}$, ${\stackrel{\mathrm{\u203e}}{R}}_{\mathrm{f}}$, and $\stackrel{\mathrm{\u203e}}{\mathrm{NSE}}$) for different grid size are calculated and used to compare with the flood simulation results (*R*_{l-lumped}, *R*_{f-lumped}, and NSE_{-lumped}) of the lumped system (shown in Table 9). Due to the large error in the rainfall simulation, storm event 3 is used as a reference event for further comparisons in this study. For the other three storm events, the difference between the gridded system and lumped system becomes more significant as the heterogeneity of the rainfall spatial distribution increases. In other words, when the spatial distribution of the rainfall is uniform, as events 1 and 2, both the lumped and gridded systems achieve good simulation results. The simulation results of the former are even superior to those of the latter in some cases, such as the NSE for event 2. However, when the spatial distribution of the rainfall is uneven, the lumped system in most cases could only simulate the overall runoff situation but fails to fully describe the runoff process, while the gridded system could obtain a better simulation, such as for event 4. Furthermore, the 3 km×3 km grid produces the most stable simulation results for events 1 and 2, which have the evenly distributed rainfall in space. It means that the simulation results may not always be improved with higher-dimension precision and a more complicated system. As the spatial distribution of rainfall becomes uneven, the simulated effect with 9 km×9 km grid declines quickly, and other grid sizes lead to better flood process simulation. The simulation results tend to improve with higher-dimension precision.

The errors in the rainfall simulation from the WRF model have a significant influence on the coupled atmospheric–hydrologic modeling system. Do the above conclusions change if the rainfall errors were eliminated? It is necessary to analyze the flood simulation results driven by the simulated rainfall, which is corrected by the measured values from the rain gauges. We assume that the WRF model is able to capture the spatial patterns of the simulated rainfall. The Thiessen polygon method is used to divide the control region of each rain gauge. The areal rainfall of the control region based on the rain gauge is regarded as the true rainfall value, which is allocated by the spatial distribution ratio from the simulated rainfall. The average value of the corrected gridded rainfall should be equal to the measured value of the rain gauge in a specific control region.

The corrected gridded rainfall is used to drive the coupled atmospheric–hydrologic modeling system, and the simulation results are shown
in Table 10. The system obtains similar simulation results of the three
different grid sizes for events 1 and 2. It is hard to say which grid division is the most outstanding, and higher resolution may not lead to better flood process simulation. The corrected rainfall leads to much better flood process simulations than the uncorrected rainfall from the WRF model. The simulations have significant improvement by using the corrected rainfall for event 3, and the system can reproduce the flood process. The
NSE_{-corrected} values are all above 0.88, and the *R*_{f-corrected}values are all lower than 3 % at different grid sizes. The system with low *R*_{f-corrected} and
*R*_{l-corrected} and high NSE_{-corrected} also performs well with different grid divisions for event 4. Higher resolution can lead to better flood process simulation for events 3 and 4. The system with 1 km×1 km grid size has the best simulation result. The conclusions are similar for both the corrected rainfall and the uncorrected rainfall. It is not always better for higher resolution, and the grid size selection has a strong relationship with the rainfall evenness. Considering the study only focusing on two semi-humid and semi-dry watersheds with limited storm events involved, the results in this study should be verified by more case studies before more general conclusions can be achieved.

Comparing the results from Tables 7 and 10, the system errors from the rainfall simulations (as shown in Table 11) can be easily obtained by the
subtraction of the corresponding values in Tables 7 and 10. For event 1, the average $\left|{R}_{\mathrm{l}}\right|-\left|{R}_{\mathrm{l}\text{-}\mathrm{corrected}}\right|$,
$\left|{R}_{\mathrm{f}}\right|-\left|{R}_{\mathrm{f}\text{-}\mathrm{corrected}}\right|$, and NSE_{-corrected}−NSE of the three different grid sizes caused by the rainfall simulations are 7.26 %, 7.00 %, and 0.1469. In the same way, the average $\left|{R}_{\mathrm{l}}\right|-\left|{R}_{\mathrm{l}\text{-}\mathrm{corrected}}\right|$, $|{R}_{\mathrm{f}}-|{R}_{\mathrm{f}\text{-}\mathrm{corrected}}|$, and NSE_{-corrected}−NSE of the three grid sizes are
7.47 %, 6.34 %, and 0.1116 for event 2. A notable case is event 3.
$\left|{R}_{\mathrm{l}}\right|-\left|{R}_{\mathrm{l}\text{-}\mathrm{corrected}}\right|$ of event 3 with
the grid size 3 km×3 km (7.96 %) is the highest among the three grid sizes, and the highest $\left|{R}_{\mathrm{f}}\right|-\left|{R}_{\mathrm{f}\text{-}\mathrm{corrected}}\right|$ (3.56 %) comes from the grid size 9 km×9 km. Due to the errors in the rainfall simulations, all the NSEs decline by more than 0.5 for the three grid sizes. For event 4, the average $\left|{R}_{\mathrm{l}}\right|-\left|{R}_{\mathrm{l}\text{-}\mathrm{corrected}}\right|$, $\left|{R}_{\mathrm{f}}\right|-\left|{R}_{\mathrm{f}\text{-}\mathrm{corrected}}\right|$, and NSE_{-corrected}−NSE of the three grid sizes
caused by the rainfall simulations are 7.32 %, 6.58 %, and 0.0991. It can
easily be found that the magnitudes of most errors in Table 11 are higher
than those of Table 10, which indicates that the accuracy of the simulated
rainfall is the main factor affecting the performance of the coupled system.
In order to improve the rainfall simulation in small- and medium-scale
catchments, radar data with high spatiotemporal resolution should be a good
choice, such as radar quantitative precipitation estimates (QPE) or quantitative precipitation forecasts (QPF) and radar data assimilation for
the NWP model (Xiao and Sun, 2007; Harader et al., 2012; Dai et al., 2019).

It should be mentioned that there is a necessity to incorporate parameter uncertainty analysis in this study. However, this will need a considerable set of the observational data (Hughes et al., 2010). Due to the lack of sufficient historical storm–flood processes, it is impossible to carry out such analyses. Nevertheless, parameter uncertainty estimations and ensemble simulations with perturbed parameters are suggested for future study when sufficient observational data are available.

This study establishes a coupled atmospheric–hydrologic modeling system with variable grid sizes in sub-humid and sub-arid area in northern China. The choice of coupling scales (1 km×1 km, 3 km×3 km, 9 km×9 km) is discussed in depth. The WRF model and the gridded Hebei model are used to establish the gridded atmospheric–hydrologic coupled system. The lumped Hebei model is also used to establish the lumped atmospheric–hydrologic coupled system, and the simulation result serves as a reference.

Contrastive analyses of the flood process simulations from the gridded atmospheric–hydrologic coupled system and the lumped atmospheric–hydrologic coupled system are carried out. Four main conclusions can be drawn: (1) the lumped system performs worse than the gridded system in most cases, while the simulation results are close to the gridded system for the storm events with uniform spatial distribution; (2) the coupled atmospheric–hydrologic systems at different grid divisions obtain similar simulation results and perform well for the storm events with uniform spatial distribution; (3) the simulation results may not always be improved with higher-dimension precision, and the grid size selection should considering the rainfall evenness; and (4) for the storm events with uneven spatiotemporal distribution, the corrected rainfall can improve the simulation results significantly, and higher resolution can lead to better flood process simulation. Flood forecasting is a stress and difficulty in the sub-humid and sub-arid area in northern China. The coupled atmospheric–hydrologic modeling system is influenced by the rainfall forecast accuracy and physics mechanisms of a hydrologic model. Further research considering the improvement of rainfall forecast and hydrologic models should be carried out to further verify the conclusions of this study.

The hydrological data used in this study are provided by Bureau of Water Resources Survey of Hebei. The NCEP/NCAR final operational global analysis (FNL) data are downloaded from https://rda.ucar.edu/datasets/ds083.2/ (last access: 30 January 2019) (NCAR, 2019). Access to the 30 m digital elevation model (DEM) can be requested through the website, http://www.gscloud.cn/sources/?cdataid=302&pdataid=10 (last access: 30 January 2019) (Geospatial Data Cloud, 2019).

All the authors contributed to the conception and the development of this paper. JT and JL contributed to establish the coupled atmospheric–hydrologic modeling system. YW and WW assisted in the calculations and the analyses. CL and CH helped with the figure production and the paper writing.

The authors declare that they have no conflict of interest.

We thank the HESS editor and two reviewers for constructive comments that help to improve the study.

This research has been supported by the National Natural Science Foundation of China (grant no. 51822906), the National Key Research and Development Project (grant no. 2017YFC1502405), the Major Science and Technology Program for Water Pollution Control and Treatment (grant no. 2018ZX07110001), and the IWHR Research & Development Support Program (grant no. WR0145B732017).

This paper was edited by Thomas Kjeldsen and reviewed by two anonymous referees.

Benoit, R., Pellerin, P., Kouwen, N., Ritchie, H., Donaldson, N., Joe, P., and Soulis, E. D.: Toward the use of coupled atmospheric and hydrologic models at regional scale, Mon. Weather Rev., 128, 1681–1706, https://doi.org/10.1175/1520-0493(2000)128<1681:TTUOCA>2.0.CO;2, 2000.

Beven, K. J. and Kirkby, M. J.: A physically based variable contributing area model of basin hydrology, Hydrol. Sci. Bull., 24, 43–69, https://doi.org/10.1080/02626667909491834, 1979.

Cassola, F., Ferrari, F., and Mazzino, A.: Numerical simulations of Mediterranean heavy precipitation events with the WRF model: a verification exercise using different approaches, Atmos. Res., 164–165, 210–225, https://doi.org/10.1016/j.atmosres.2015.05.010, 2015.

Chen, F., Liu, C., and Dudhia, J.: A sensitivity study of high-resolution regional climate simulations to three land surface models over the western United States, J. Geophys. Res.-Atmos., 119, 7271–7291, https://doi.org/10.1002/2014JD021827, 2014.

Dai, Q., Yang, Q., Han, D., Rico‐Ramirez, M. A., and Zhang, S.: Adjustment of radar-gauge rainfall discrepancy due to raindrop drift and evaporation using the Weather Research and Forecasting model and dual-polarization radar, Water Resour. Res., 55, 9211–9233, https://doi.org/10.1029/2019WR025517, 2019.

Di, Z., Duan, Q., Gong, W., Wang, C., Gan, Y., Quan, J., Li, J., Miao, C., Ye, A., and Tong, C.: Assessing WRF model parameter sensitivity: A case study with 5 day summer precipitation forecasting in the Greater Beijing Area, Geophys. Res. Lett., 42, 579–587, https://doi.org/10.1002/2014GL061623, 2015.

Duan, Q., Sorooshian, S., and Gupta, V. K.: Optimal use of the SCE-UA global optimization method for calibrating watershed models, J. Hydrol., 158, 265–284, https://doi.org/10.1016/0022-1694(94)90057-4, 1994.

Evans, J. P., Ekström, M., and Ji, F.: Evaluating the performance of a WRF physics ensemble over South-East Australia, Clim. Dynam., 39, 1241–1258, https://doi.org/10.1007/s00382-011-1244-5, 2012.

Geospatial Data Cloud: Global Digital Elevation Model DEM 30M, available at: http://www.gscloud.cn/sources/?cdataid=302&pdataid=10, last access: 30 January 2019.

Goutal, N. and Sainte-Marie, J.: A kinetic interpretation of the section-averaged Saint-Venant system for natural river hydraulics, Int. J. Numer. Meth. Fl., 67, 914–938, https://doi.org/10.1002/fld.2401, 2011.

Grell, G. A. and Freitas, S. R.: A scale and aerosol aware stochastic convective parameterization for weather and air quality modeling, Atmos. Chem. Phys., 14, 5233–5250, https://doi.org/10.5194/acp-14-5233-2014, 2014.

Harader, E., Borrell-Estupina, V., Ricci, S., Coustau, M., Thual, O., Piacentini, A., and Bouvier, C.: Correcting the radar rainfall forcing of a hydrological model with data assimilation: application to flood forecasting in the Lez catchment in Southern France, Hydrol. Earth Syst. Sci., 16, 4247–4264, https://doi.org/10.5194/hess-16-4247-2012, 2012.

Hong, S. Y. and Lee, J. W.: Assessment of the WRF model in reproducing a flash-flood heavy rainfall event over Korea, Atmos. Res., 93, 818–831, https://doi.org/10.1016/j.atmosres.2009.03.015, 2009.

Hong, S. Y., Noh, Y., and Dudhia, J.: A new vertical diffusion package with an explicit treatment of entrainment processes, Mon. Weather Rev., 134, 2318–2341, https://doi.org/10.1175/MWR3199.1, 2006.

Horton, R. E.: The role of infiltration in the hydrologic cycle, Eos. T. Am. Geophys. Un., 14, 446–460, https://doi.org/10.1029/TR014i001p00446, 1933.

Hostache, R., Matgen, P., Montanari, A., Montanari, M., Hoffmann, L., and Pfister, L.: Propagation of uncertainties in coupled hydro-meteorological forecasting systems: A stochastic approach for the assessment of the total predictive uncertainty, Atmos. Res., 100, 296–274, https://doi.org/10.1016/j.atmosres.2010.09.014, 2011.

Hughes, D. A., Kapangaziwiri, E., and Sawunyama, T.: Hydrological model uncertainty assessment in southern Africa, J. Hydrol., 387, 221–232, https://doi.org/10.1016/j.jhydrol.2010.04.010, 2010.

Janjić, Z. I.: The step-mountain eta coordinate model: further developments of the convection, viscous sublayer, and turbulence closure schemes, Mon. Weather Rev., 122, 927–945, https://doi.org/10.1175/1520-0493(1994)122<0927:TSMECM>2.0.CO;2, 1994.

Jha, S. K., Zhao, H., Woldemeskel, F. M., and Sivakumar, B.: Network theory and spatial rainfall connections: An interpretation, J. Hydrol., 527, 13–19, https://doi.org/10.1016/j.jhydrol.2015.04.035, 2015.

Kain, J. S.: The Kain-Fritsch convective parameterisation: an update, J. Appl. Meteorol., 43, 170–181, https://doi.org/10.1175/1520-0450(2004)043<0170:TKCPAU>2.0.CO;2, 2004.

Li, J., Chen, Y., Wang, H., Qin, J., Li, J., and Chiao, S.: Extending flood forecasting lead time in a large watershed by coupling WRF QPF with a distributed hydrological model, Hydrol. Earth Syst. Sci., 21, 1279–1294, https://doi.org/10.5194/hess-21-1279-2017, 2017.

Lin, Y. L., Farley, R. D., and Orville, H. D.: Bulk parameterization of the snow field in a cloud model, J. Clim. Appl. Meteorol., 22, 1065–1092, https://doi.org/10.1175/1520-0450(1983)022<1065:BPOTSF>2.0.CO;2, 1983.

Liu, J., Bray, M., and Han, D.: Sensitivity of the weather research and forecasting (WRF) model to downscaling ratios and storm types in rainfall simulation, Hydrol. Process., 26, 3012–3031, https://doi.org/10.1002/hyp.8247, 2012.

Liu, J., Wang, J., Pan, S., Tang, K., Li, C., and Han, D.: A real-time flood forecasting system with dual updating of the NWP rainfall and the river flow, Nat. Hazards, 77, 1161–1182, https://doi.org/10.1007/s11069-015-1643-8, 2015.

Liu, J., Tian, J., Yan, D., Li, C., Yu, F., and Shen, F.: Evaluation of Doppler radar and GTS data assimilation for NWP rainfall prediction of an extreme summer storm in northern China: from the hydrological perspective, Hydrol. Earth Syst. Sci., 22, 4329–4348, https://doi.org/10.5194/hess-22-4329-2018, 2018.

Miao, S., Chen, F., and Li, Q.: Impacts of urban processes and urbanization on summer precipitation: A case study of heavy rainfall in Beijing on 1 August 2006, J. Appl. Meteorol. Clim., 50, 806–825, https://doi.org/10.1175/2007MWR2106.1, 2011.

Myung, J. I.: The importance of complexity in model selection, J. Math. Psychol., 44, 190–204, https://doi.org/10.1006/jmps.1999.1283, 2000.

Myung, J. I., Tang, Y., and Pitt, M. A.: Evaluation and comparison of computation models, Meth. Enzymol., 454, 287–304, https://doi.org/10.1016/S0076-6879(08)03811-1, 2009.

NCAR: NCEP FNL Operational Model Global Tropospheric Analyses, continuing from July 1999, available at: https://rda.ucar.edu/datasets/ds083.2, last access: 30 January 2019.

Rogelis, M. C. and Werner, M.: Streamflow forecasts from WRF precipitation for flood early warning in mountain tropical areas, Hydrol. Earth Syst. Sci., 22, 853–870, https://doi.org/10.5194/hess-22-853-2018, 2018.

Tian, J., Liu, J., Wang, J., Li, C., Yu, F., and Chu, Z.: A spatio-temporal evaluation of the WRF physical parameterisations for numerical rainfall simulation in semi-humid and semi-arid catchments of Northern China, Atmos. Res., 191, 141–155, https://doi.org/10.1016/j.atmosres.2017.03.012, 2017a.

Tian, J., Liu, J., Yan, D., Li, C., and Yu, F.: Numerical rainfall simulation with different spatial and temporal evenness by using a WRF multiphysics ensemble, Nat. Hazards Earth Syst. Sci., 17, 563–579, https://doi.org/10.5194/nhess-17-563-2017, 2017b.

Tian, J., Liu, J., Yan, D., Li, C., and Yu, F.: Ensemble flood forecasting based on a coupled atmospheric-hydrological modeling system with data assimilation, Atmos. Res., 224, 127–137, https://doi.org/10.1016/j.atmosres.2019.03.029, 2019.

Verbunt, M., Zappa, M., Gurtz, J., and Kaufmann, P.: Verification of a coupled hydrometeorological modelling approach for alpine tributaries in the Rhine basin, J. Hydrol., 324, 224–238, https://doi.org/10.1016/j.jhydrol.2005.09.036, 2006.

Wagner, S., Fersch, B., Yuan, F., Yu, Z., and Kunstmann, H.: Fully coupled atmospheric-hydrological modeling at regional and long-term scales: Development, application, and analysis of WRF-HMS, Water Resour. Res., 52, 3187–3211, https://doi.org/10.1002/2015WR018185, 2016.

Wan, Q. and Xu, J.: A numerical study of the rainstorm characteristics of the June 2005 flash flood with WRF/GSI data assimilation system over south-east China, Hydrol. Process., 25, 1327–1341, https://doi.org/10.1002/hyp.7882, 2011.

Wang, C. C., Kuo, H. C., Yeh, T. C., Chung, C. H., Chen, Y. H., Huang, Y. S., Wang, Y., and Liu, C. H.: High-resolution quantitative precipitation forecasts and simulations by the Cloud-Resolving Storm Simulator (CReSS) for Typhoon Morakot (2009), J. Hydrol., 506, 26–41, https://doi.org/10.1016/j.jhydrol.2013.02.018, 2013.

Wang, H., Sun, J., Fan, S., and Huang, X.: Indirect assimilation of radar reflectivity with WRF 3D-Var and its impact on prediction of four summertime convective events, J. Appl. Meteorol. Clim., 52, 889–902, https://doi.org/10.1175/JAMC-D-12-0120.1, 2013.

Wu, J., Lu, G., and Wu, Z.: Flood forecasts based on multi-model ensemble precipitation forecasting using a coupled atmospheric-hydrological modeling system, Nat. Hazards, 74, 325–340, https://doi.org/10.1007/s11069-014-1204-6, 2014.

Xiao, Q. and Sun, J.: Multiple-Radar Data Assimilation and Short-Range Quantitative Precipitation Forecasting of a Squall Line Observed during IHOP_2002, Mon. Weather Rev., 135, 3381–3404, https://doi.org/10.1175/MWR3471.1, 2007.

Yucel, I., Onen, A., Yilmaz, K. K., and Gochis, D. J.: Calibration and evaluation of a flood forecasting system: Utility of numerical weather prediction model, data assimilation and satellite-based rainfall, J. Hydrol., 523, 49–66, https://doi.org/10.1016/j.jhydrol.2015.01.042, 2015.

Yue, B. J., Shi, Z. H., and Fang, N. F.: Evaluation of rainfall erosivity and its temporal variation in the Yanhe River catchment of the Chinese Loess Plateau, Nat. Hazards, 74, 585–602, https://doi.org/10.1007/s11069-014-1199-z, 2014.

Zhao, R.: The Xinanjiang model applied in China, J. Hydrol., 135, 371–381, https://doi.org/10.1016/0022-1694(92)90096-E, 1992.