*2.1. Study Area*

We studied two mountainous catchments of the Daqing River (i.e., Fuping (2210 km2) and Zijingguan (1760 km2)). Two catchments lie along the upper reach of the Shazhi River on the south branch and the upper reach of the Juma River on the north branch, respectively (Figure 1). Low vegetation coverage, bare hills, mountains, and thin soil layers cause uneven infiltration capacity distribution across the study areas, and the surface often features extra-infiltration flow, which means surface flow occurs when the intensity of infiltration exceeds the intensity of precipitation. Rainfall occurs primarily during the flood season from late May to early September. The combined effects of the western Pacific subtropical high, the westerly cold vortex, and the western Taihang Mountain uplift influence the heavy convective rainfall that prevails in the region. Additionally, due to steep terrain, it is easy to form high intensity, short duration floods with a large peak flow. These floods cause damage in the region. According to statistics of torrential rain and flood data from 1958 to 2015 in both catchments, floods that occurred more than once a decade occurred five times in Fuping and six times in Zijingguan, with increasing frequency in extreme events in recent years.

**Figure 1.** Locations of the Fuping and Zijingguan catchments covered by the weather research forecasting (WRF) nested domains and the weather radar.

#### *2.2. Storm Events*

Four storm events in the Fuping and Zijingguan catchments were screened. The events were divided according to their spatial and temporal distribution. In comparison to the southern part of China, it is difficult to find absolute homogeneous rainfall events in northern China, either in spatial or temporal dimensions. Herein, we calculated the storm events' coefficient of variance (*Cv*) values to indicate the relative response to the spatiotemporal homogeneousness of the flood events [34]. *Cv* values were calculated as follows:

$$Cv = \sqrt{\frac{\sum\_{i=1}^{N} \left(\frac{x\_i}{N} - 1\right)^2}{N}} \tag{1}$$

where *x* is the average of *xi*. When calculating *Cv* in the temporal dimension, *xi* is the catchment areal rainfall at the *i*th hour and *N* is the total number of hours of the storm event. In this case, time *Cv* represents the temporal deviation of the catchment areal rainfall at each time step. When calculating *Cv* in the spatial dimension, *xi* is the 24 h rainfall accumulation at the *i*th gauge and *N* is the number of rain gauges. In this case, the *Cv* value reflects the spatial deviation of the rainfall accumulation at each rain gauge.

The results are shown in Table 1, where larger *Cv* values represent a more uneven spatio-temporal distribution. Based on this criterion, the spatio-temporal consistency of the rainfall-runoff events listed in Table 1 were classified. Events 1 and 2 had relatively uniform spatial and temporal distributions, whereas Events 3 and 4 featured rainfall with uneven spatial and temporal distributions. As shown in Figure 2, the rainfall-runoff events had a flood recession time with different lengths, so we used a 72-h time window for Event 1 and Event 2, a 60-h time window for Event 3, while for Event 4, a 36-h time window was adopted. Different time windows were used to calculate the statistics when evaluating the forecasting results. Event 4 occurred on 21 July 2012, bringing the largest flood disaster in 60 years to Beijing and the Daqing River. The rainfall at the largest monitoring point corresponded to an event that occurred once in more than 500 years. In the Zijingguan catchment, its 24 h maximum single station cumulative rainfall was 355 mm, and the maximum flow at the outlet reached 2580.0 m3/s. Figure 2 and Table 1 show other rainfall and flood characteristics.

**Table 1.** Storm events and spatio-temporal rainfall evenness characterized by *Cv.*


**Figure 2.** Observed rainfall-runoff processes for the four selected storm events.

#### **3. Three Atmospheric–Hydrologic Coupled Systems**

*3.1. The Numerical Weather Prediction (NWP) Model*

The mesoscale NWP model has limited skill in convective precipitation forecast, even in the WRF model. One of the considerable reasons for the poor performance is due to its nonlinear and chaotic nature, since solving quasistationary meso-β and meso-γ processes

is complicated [35,36], and the power spectrum of the turbulence in convective motion has a resolution of a few kilometers [37,38]. 3DVar data assimilation could merge the fine rainfall information into the modeling system to obtain accurate high-precision rainfall data. The following are detailed descriptions of the WRF model configurations and the 3DVar data assimilation techniques used in this study.

#### 3.1.1. Weather Research Forecasting (WRF) Model Configurations

As a next-generation mesoscale forecasting model and data assimilation system, the simulated scale used in the WRF spans tens of meters to thousands of kilometers, and is mainly used to enhance understanding and forecasting of mesoscale weather, assist operational forecasting, and improve atmospheric research. Detailed descriptions of specific WRF model settings are shown in Table 2 and Figure 1. The inner research area includes the Taihang Mountains, Bohai Sea, and major cities in Beijing, Tianjin, and Hebei. A warm-up time period is necessary for both WRF and WRF-Hydro. A seven-day length warm-up time period was set to generate the restart file, which serves mainly as the starting condition for WRF, and was also used to provide the initial soil moisture condition for WRF-Hydro to generate more realistic runoff. The WRF has an output time step of 1 h and an integration step of 6 s. Global Forecast System (GFS) data in each 6 h window was to provide the initial side boundary and atmospheric background conditions.



The parameterization scheme chosen can affect the mode of operation of the WRF; different parameterization schemes have different emphases on physical processes. The diversity of parameterization schemes and the corresponding differences in rainfall events in specific regions result in difficulties inherent to the accurate simulation of spatial and temporal rainfall distributions and cumulative rainfall. For the four rainfall processes in this study, we chose the optimal physical parameterization based on the relevant experimental research on the selection of sensitive parameterization schemes for ensemble rainfall forecasting (more details can be found in Tian et al. [39,40]). The microphysical parameterization schemes chosen for our forecasts included Purdue-Lin (Lin) [41] and WRF Single-Moment 6 (WSM6) [42]; the cumulus convection schemes included the Kain-Fritsch (KF) [43] and Grell-Devenyi (GD) [44] ensembles; and the planetary boundary layer schemes included the Yonsei University (YSU) [42] and Mellor-Yamada-Janjic (MYJ) [45] schemes. The specific parameterization details of the four studied storms are shown in Table 3. It should be noted that the cumulus parameterization was not used in the innermost domain, where the convective rainfall generation was assumed to be explicitly resolved.

**Table 3.** Optimal physical parameterization for the studied storm events.


#### 3.1.2. Data Assimilation with WRF-3DVar

3D variational assimilation has numerous advantages, for example, that the objective function contains physical processes and utilizes the model itself as a dynamic constraint (i.e., it can efficiently represent complex nonlinear constraint relationships). The objective function of WRF-3DVar assimilation can be expressed as follows:

$$f(X) = \frac{1}{2}(X - X\_b)^T B^{-1} (X - X\_b) + \frac{1}{2}(H(X) - \Upsilon\_0)^T R^{-1} (H(X) - \Upsilon\_0) \tag{2}$$

where *X* is a parameter reflecting the atmosphere and surface conditions; *Xb* is the background field at the time of change; *B* is the corresponding background field error covariance matrix; and *H* is the observation function containing *X* variables. *H* can change the variables in the atmospheric model from model space projected into the observation space. *Y*0 represents the assimilated observation vector and *R* is the error covariance matrix from observation.

Global Telecommunications System (GTS) data and Doppler radar data were used as the assimilation data in this study. In previous studies, we have demonstrated that the assimilation of radar reflectivity can improve local rainfall processes, especially for data <500 m, and can thus improve the reliability of assimilation information for forecasting systems [46,47]. Taking into account the spatial resolution of data from different sources and maintaining the stability of the atmospheric motion, we chose to assimilate <500 m radar reflectivity data in the inner domain (Dom2) and GTS data in the outer domain (Dom1). This assimilation method was supported by multi-source data assimilation experiments carried out in the same study area [47].

GTS data are a collection of traditional ground and upper air meteorological data including barometric pressure, wind direction and velocity, temperature, and humidity. GTS data have wide coverage in both the vertical and horizontal directions. These data are updated at 6 h intervals and are therefore widely used in large-scale atmospheric studies. GTS data provided by the National Center for Atmospheric Research (NCAR) were assimilated for the outer domain, sourced from sounders, ground-based weather stations, pilot balloons, and aircraft observations. An observation preprocessor with defined observation error covariance was employed to ensure quality control of the GTS data prior to assimilation, and a default U.S. Air Force (AFWA) OBS error file was also used for processing the GTS data and identifying instrumental and sensor errors.

The forecasting radar site of the Shijiazhuang s-band Doppler weather radar and specific coverage is shown in Figure 1. Its reflectance spatial resolution was 1 km × 1◦ with a scan radius of 250 km (i.e., it covered both the Fuping and Zijingguan catchments). The radar completed a cycle every 6 min at nine different scan elevation angles (0.5◦, 1.5◦, 2.4◦, 3.4◦, 4.3◦, 6.0◦, 9.9◦, 14.6◦, and 19.5◦). Raw data for the radar were obtained from the National Meteorological Administration of China and, following quality control, the radar reflectivity was programmed to conform to the WRF-3DVar data format.

The absorption of radar reflectivity data by the WRF-3DVar module presupposes that the total water mixing ratio is used as a control variable and that a warm rainfall process is simultaneously introduced into the assimilation module. Assuming a Marshall-Palmer raindrop size distribution and no ice relative reflectivity effect, the following equation can be derived for radar reflectivity *Z* by introducing a rainwater mixing ratio *qr*, in which *ρ* is the density of air in kg m<sup>−</sup>3:

$$Z = 43.1 + 17.5 \log(\rho q\_r) \tag{3}$$

For the four storm events, we selected a uniform 24 h time window that completely described the entire process of rainfall. The duration of the WRF is required to cover the start and end of the rainfall time window; therefore, the assimilation data started to be merged 6 h before the storm events, then assimilated every 6 h until five assimilation cycles were completed, for a total running duration of 36 h (6 × 6 h). Figure 3 presents a schematic diagram of the assimilation using Event 4 for further explanation. As illustrated in Figure 3, there was a seven day warm-up time period (indicated with a dashed line) and a 6 h spin-up time period (indicated with dotted lines) for WRF, and run1 presents the initial WRF run with no data assimilation. Data assimilation started at 00:00 on 21 July 2012 and, subsequently, run2, run3, run4, run5, and run6 were executed at 6 h intervals (i.e., at 00:00, 06:00, 12:00, and 18:00 on 21 July 2012 and 00:00 on 22 July 2012, respectively). The first output file generated in the previous run was used to provide the initial and lateral boundary conditions for subsequent runs and was treated as a benchmark reflecting improvements in rainfall simulation after data assimilation.

**Figure 3.** The cycling data assimilation process for Event 4.

#### *3.2. Hydrological Models with Differing Complexity*

#### 3.2.1. The Lumped Hebei Model

ÿ

¨

¨

The lumped Hebei model was developed based on rainfall-runoff mechanisms identified in semi-humid and semi-arid catchments in northern China and has been applied to real-time flood forecasting in Hebei Province. The model comprehensively considers the two inflow mechanisms of infiltration excess and saturation excess (shown in Figure 4). *W'* and *WMM* represent the point storage capacity and the storage capacity of the maximum point in catchments. *γ* and *λ* are the proportion of the infiltration capability beyond a certain value and the proportion of runoff generation in the catchment. When rainfall intensity was greater than infiltration intensity during the studied period, surface runoff occurs, whereupon continuous infiltration supplements the soil moisture deficiency until the system reaches the point soil water capacity, generating underground runoff. The model divides soil into two layers, in which soil depth should be determined according to the conditions of the catchment. The volume of infiltration (*ft*), surface (*rs*), and groundwater (*rg*) can be calculated using Equations (2)–(4) as follows:

$$f\_t = \sum\_{i=1}^t \overline{f}\_i = \left( i - \frac{i^{(1+n)}}{(1+n)f\_m^n} \right) e^{-um} + f\_c \tag{4}$$

$$r\_s = p\_t - f\_t \tag{5}$$

where *pt* (mm/h) is the rainfall rate during the *i*th hour and *f i* (mm/h) is the *i*th hour empirical infiltration rate computed by a variant of Horton equation; *m* (mm) is the surface soil moisture; and *u* is an index accounting for the decreasing infiltration rate with

increasing soil moisture. *fc* (mm/h) and *fm* (mm/h) are the stable infiltration rate and the infiltration capacity, respectively. It follows that:

$$r\_{\mathcal{S}} = \begin{cases} \begin{array}{l} f\_t + p\_a - e\_t - w\_m \\ f\_t + p\_a - e\_t - w\_m + w\_m \left(1 - \frac{f\_t + p\_a}{w\_m}\right)^{1 + b} \end{array} & p\_a + f\_t < w\_m \end{cases} \tag{6}$$

where *et* is the evaporation volume; *pa* (mm) represents the rainfall influenced; *wm* (mm) is the average storage capacity; and *b* is a coefficient of the water storage capacity curve.

**Figure 4.** Rainfall-runoff mechanisms for the grid-based Hebei model.

#### 3.2.2. The Grid-Based Hebei Model

The grid-based Hebei model retains the runoff generation concepts behind the lumped Hebei model. Furthermore, it develops terrain index *T* [48,49] to obtain a grid-type representation of soil moisture storage and infiltration capacity. *T* = ln(*α*/tan*β*), where *α* and *β* are the grid scale and shape parameters, respectively. The value of *T* in the Fuping and Zijingguan catchments were statistically analyzed. Experimentation demonstrated that the cumulative distribution curves of *T* values between different regions are always similar [50]. The *T* values of different grids in the same region can be expressed by a parabolic empirical formula. Hence, the soil water storage capacity and infiltration capacity of different grids can be determined as followed:

$$\frac{\mathcal{W}\_{i}}{\mathcal{W}MM} = \exp\left\{-\left[\frac{\ln(T\_{i} - T\_{\text{min}} + 1)}{\alpha}\right]^{\beta}\right\} \tag{7}$$

$$\frac{f\_i}{f\_m} = \left\{ 1 - \left[ 1 - \exp\left[ -\frac{1}{a} [\ln(T\_i - T\_{\rm min} + 1)]^\beta \right]^b \right]^{1/n} \right\}^{1/n} \tag{8}$$

where *i* is a certain grid cell, and *Ti* and *Wi* (mm) are the corresponding moisture storage capacity and terrain index of the *i*th grid. *T*min is the minimum value of the grid terrain index. In a non-channel grid, surface runoff outflow *qso* = *qsi* + *qs*, where *qs* is the grid surface runoff. *qs* occurs when the rainfall intensity exceeds *fm*. *qsi* denotes the runoff values generated via the surrounding upstream grid. Similarly, the generated groundwater runoff from upstream grid (i.e., *qg*, is transmitted to the groundwater runoff outflow according to *qgo*, *qgo* = *qgi*·(<sup>1</sup> − *μ*) + *qg*, where *μ* is the ratio coefficient). The grid-based Hebei model adopts a simplified form of the Saint Venant equations for confluence calculation. Due to perennial channel water shortage and substantial channel seepage in the study area, an additional Horton infiltration equation was considered in the grid-based Hebei model [50].

#### 3.2.3. The WRF-Hydro Modeling System

The open-code WRF-Hydro model is a multi-core parallel, multi-physics, multiscale, fully distributed hydraulic model. It can be operated only with meteorological data or coupled with the WRF model to constitute the atmospheric-hydrological model system. Its multiscale 3D land surface hydrological simulation mode can improve the one-dimensional vertical generalization of water transport using the original LSM Noah and Noah-MP. WRF-Hydro can use additional modules to achieve lateral flow exchange between the surface and subsurface, thereby improving upon the "column-only" one-dimensional vertical structure. It uses a disaggregation-aggregation solution module between the land model and terrain routing grid, which enables convergence processes occurring at a smaller resolution. WRF-Hydro can set the variables of the scale factor including soil moisture and excess infiltration into the fine grid values.

Noah or Noah-MP LSM are connected to WRF-Hydro to calculate water and heat flux exchange processes. In Noah-MP, a column cylinder structure is used to substitute soil layers into thicknesses, from top to bottom, of 0–10 cm, 10–40 cm, 40–100 cm, and 100–200 cm. The inputs required by WRF-Hydro include meteorological forcing data such as rainfall, air temperature, relative air humidity, surface pressure, wind speed, downward longwave radiation, and downward solar radiation. Gochis and Chen [51] described a sub-grid, spatially weighted disaggregation-aggregation solution to coordinate the land model and terrain routing grid in WRF-Hydro. Its runoff generation method uses a simple water balance (SWB) [52]. Similar to the Hebei model, when the rainfall capacity exceeds the infiltration capacity, a surface infiltration excess occurs in the top soil layer and the corresponding change in surface water depth *h* (m) is:

$$\frac{\partial h}{\partial t} = \frac{\partial p\_c}{\partial t} \left\{ 1 - \frac{\left[\sum\_{k=1}^4 Z\_i(\delta\_s - \delta\_k)\right] \left[1 - \exp\left(-S \frac{R\_{df}}{R\_{fd}} \frac{\Delta t}{86400}\right)\right]}{p\_c + \left[\sum\_{k=1}^4 Z\_i(\delta\_s - \delta\_k)\right] \left[1 - \exp\left(-S \frac{R\_{df}}{R\_{fd}} \frac{\Delta t}{86400}\right)\right]} \right\} \tag{9}$$

where *h* (m) refers to the change in surface water depth; Δ*t* (s) is the model time step; *k* is the integer number of the soil layer (i.e., 1–4), *Zk* (m) and *δk* (m3m−3) are the depth and grid soil moisture of the *k*th soil layer; *δs* (m3m−3) is the maximum soil moisture content; *S* is a coefficient given by Richards' equation to regulate runoff infiltration; and *Rdt* and *Rfd* represent the tunable surface infiltration coefficient and saturated hydraulic conductivity, respectively.

Subsurface routing followed the approach of Wigmosta and Lettenmaier [53] by using a quasi-3D flow taking into account topography, saturated soil depth, and saturated hydraulic conductivity with soil depth. Overland routing is a fully unsteady finite-difference and diffuse wave approach, implemented as described in Downer et al. [54]. River routing is similar to the overland case, using an explicit, one-dimensional, variable time-step diffusion wave equation. Details of specific rainfall-runoff processes can be obtained from the official user guide [55].

#### *3.3. Establishment of Three Coupled Atmospheric-Hydrological Systems*

Figure 5 shows a flood prediction diagram for coupled systems comprising three types of hydrological models, WRF models, and WRF-3DVar data assimilation modules. Rainfall in the WRF model was enhanced through optimal parameter scenarios and assimilation. GTS data and radar reflectivity were assimilated every 6 h to generate 3 km grid data in the inner layer Dom2. To eliminate the deviation of rainfall, for each storm, the best performance in the assimilated rainfall ensemble was selected as the forecast forcing data. Subsequently, predicted rainfall was converted into discharges using the aforementioned three models, with parameters calibrated using 17 historical rainfall-runoff events from the two studied catchments since the 1980s. For the calibration of WRF-Hydro, previous parameter sensitivity analysis in the study area noted several key parameters that can

cause large fluctuations in flooding prediction including the runoff infiltration parameter (*REFKDT*), the channel Manning roughness parameter (*MannN*), the surface retention depth scaling parameter (*RETDEPRTFAC*), and the overland flow roughness scaling parameter (*OVROUGHRTFAC*). To evaluate the impact of data assimilation on the obtained runoff, WRF outputs under the selected optimal parameterization scheme were also used to drive the model for runoff forecasting; the results of the three systems before and after data assimilation were compared as described below.

**Figure 5.** Framework for the coupled atmospheric-hydrologic modeling systems used for rainfallrunoff prediction.

The three coupling systems tested showed three types of rainfall input at different resolutions. The lumped model gives selected grid-averaged rainfall data across the sub-basins, whereas the grid-based Hebei model and WRF-Hydro directly pass the grid coordinates to locate forcing data. The former has a rainfall resolution of 3 km, and the latter has a grid division factor of 10 to allow further downscaling from the WRF assimilation data to 300 m of routing data.

The Nash efficiency coefficient (*NSE*), relative flood peak (*Rp*), and relative flood volume (*Rv*) were used to analyze the flow discharge forecasts as follows:

$$NSE = 1 - \frac{\sum\_{i=1}^{\overline{n}} \left( q\iota\_i - q\_i \right)^2}{\sum\_{i=1}^{\overline{n}} \left( q\_i - \overline{q} \right)^2} \tag{10}$$

$$R\_p = \left(q\iota\_p - q\_p\right) / q\_p \tag{11}$$

$$R\_{\upsilon} = (r t\_{\upsilon} - r\_{\upsilon}) / r\_{\upsilon} \tag{12}$$

where *i* is the time step; *n* is the total *i* of flood processes; *qi*, *qi*, and *q* are the simulated, observed, and averaged flood discharge, respectively; *qp* and *qp* are the simulated and

observed discharge; respectively; and *<sup>r</sup>v* and *rv* are the simulated and observed flood volumes, respectively.
