The WRF-Driven Grid-Xin’anjiang Model and Its Application in Small and Medium Catchments of China

Abstract

The distributed Grid-Xin’anjiang (Grid-XAJ) model is very sensitive to the spatial and temporal distribution of data when used in humid and semi-humid small and medium catchments. We used the successive correction method to merge the gauged rainfall with rainfall forecasted by the Weather Research and Forecasting (WRF) model to enhance the spatiotemporal accuracy of rainfall distribution. And we used the Penman–Monteith equation to calculate the potential evapotranspiration (PE_PM). Then, we designed two forcing scenarios (WRF-driven rainfall (Wr) + PE_PM, WRF-merged rainfall (Wm) + PE_PM) to drive the Grid-XAJ model for flood forecasting. We found the WRF-driven Grid-XAJ model held significant potential in flood forecasting. The Grid-XAJ model provided only an approximation of flood hygrographs when driven by scenario Wr + PE_PM. The results in scenario Wm + PE_PM showed a high degree-of-fit with observed floods with mean Nash–Sutcliffe efficiency coefficient (NSE) values of 0.94 and 0.68 in two catchments. Additionally, scenario Wm + PE_PM performed better flood hygrographs than scenario Wr + PE_PM. The flood volumes and flow peaks in scenario Wm + PE_PM had an obvious improvement compare to scenario Wr + PE_PM. Finally, we observed that the model exhibited superior performance in forecasting flood hydrographs, flow peaks, and flood volumes in humid catchments compared with semi-humid catchments.

1. Introduction

Floods are one of the most serious natural disasters, and in recent years, they have received increasing attention [1,2]. The primary economic and demographic focal points of China are primarily concentrated in humid and semi-humid regions, thus emphasizing the paramount significance of studying flood forecasting in these areas [3]. In recent years, the risk of extreme floods in the large basins of China has been substantially reduced due to improvements in flood forecasting [4]. However, the floods that occur in small and medium catchments (SMCs) are characterized by rapidity and suddenness, and usually have a short flood duration and cause a rapid rise in the flow peak, which makes flood prediction difficult [5]. Moreover, flood forecasting in both humid SMCs and semi-humid SMCs faces additional challenges due to the fact that flood forecasting in SMCs is more sensitive to the spatial and temporal distribution of data.

Hydrological models serve as important tools in flood forecasting and have been widely used in various regions [5,6]. Hydrological models can be classified into lumped, semi-distributed, and distributed models on the basis of the hydrological processes and watershed attributes [7,8,9]. In this study, we used the Grid-XAJ model [10] for flood forecasting. It is a distributed version of the well-known Xin’anjiang model. This model has been widely used in numerous catchments throughout China [11,12]. It operates at a high resolution (≤1 km) and conducts flood calculations for each grid cell, thereby incorporating the spatial variability of rainfall and the underlying surface conditions more effectively. Traditionally, the Grid-XAJ model is driven by gauged rainfall and pan evaporations [10], which overlook the comprehensive influence of meteorology on hydrological processes. Additionally, distributed hydrological models require grid data driving [13]. The grid data for the Grid-XAJ model are usually generated using spatial interpolation methods. The accuracy of station observations is relatively high, and the uncertainty in expressing the characteristics of a meteorological element at the basin scale using station observations is largely due to interpolation or data assimilation processes.

The precise calibration of parameters in the hydrological model is crucial for accurately simulating floods [14,15]. Some parameters of the Grid-XAJ model can be a priori estimated using underlying surface data. For instance, the channel bed slope (S_oc) can be estimated utilizing the digital elevation model (DEM), and the leaf area index (LAI) can be estimated based on land use and land cover (LULC) data. On the other hand, other parameters necessitate calibration. Over time, the Grid-XAJ model has undergone continuous development and refinement.

Acquiring reliable rainfall spatial information on high temporal resolution is essential to improve the ability of flood forecasting for the Grid-XAJ model. Typically, researchers interpolate gauged rainfall onto computational units to drive the Grid-XAJ model [16]. However, rainfall spatial interpolation is a non-trivial task due to many challenging factors, such as the rainfall’s irregular distribution in space, and the accuracy of the interpolation results is sometimes unsatisfactory [17,18]. Several methods have been employed to improve the interpolation results, such as graph neural networks (GNNs) [19,20] and graph for spatial interpolation (GSI) [21]. Additionally, various numerical weather prediction (NWM) models are used to achieve high spatial resolution in rainfall data [22]. For example, the Weather Research and Forecasting (WRF) model has a proven potential in depicting the spatial distribution of rainfall [23]. It has been used by Yao [24] to generate rainfall data for the Grid-XAJ model. However, WRF-driven rainfall often underestimates both the rain area and its intensity [25,26]. To enhance the reliability of rainfall predictions, a range of methods have been employed to merge WRF rainfall with multiple sources [27]. In this study, we used a successive correction method [28] to merge the WRF-driven rainfall with gauged rainfall to improve the spatio-temporal accuracy of rainfall distribution.

Pan evaporation is another type of input gauged data for the Grid-XAJ model, which is converted into potential evapotranspiration (PE) using the pan coefficient (KC). Evapotranspiration (ET) plays a crucial role in facilitating water–energy exchange between the land surface and atmosphere [29,30,31]. The accurate assessment of PE is significant for studying hydrological cycle processes and achieving full coupling between hydrological models and meteorological models. The hourly PE in the Grid-XAJ model is calculated for each grid by evenly distributing the daily pan evaporation across each hour (PE_pan). However, usually only one or two pans are available within a catchment, and these pans exclusively capture daily evaporation. So the PE_pan may not resonate with real-world conditions [32,33]. Various methods can improve the spatiotemporal accuracy of PE distribution. With the advancement of remote sensing technology, a growing number of datasets now directly offer high-spatial-resolution PE data with reliable accuracy at no cost. Notable examples include the Moderate Resolution Imaging Spectroradiometer (MODIS) [34] and the Global Land Evaporation Amsterdam Model (GLEAM) [35]. But their temporal resolution is sometimes not satisfactory. Acquiring high-spatiotemporal-resolution meteorological variables such as temperature and wind speed is considerably easy. Hence, machine learning methods with the inputs of meteorological data have been extensively employed for the best PE estimates in recent years [36,37]. In addition, using energy-based methods to estimate PE with high-spatiotemporal-resolution meteorological data is reliable. For example, the Penman–Monteith equation, which is based on energy balance, has been used by Meng [38] to effectively estimate the temporal and spatial distribution of PE in the Grid-XAJ model. We employed the Penman–Monteith equation to evaluate the potential evapotranspiration (PE_PM) using meteorological data generated by the WRF model in this paper.

The objective of this study is to enhance the accuracy of the spatial and temporal distribution of rainfall and PE in the Grid-XAJ model via the utilization of the WRF model. This enhancement aims to improve the flood simulations of the Grid-XAJ model. To accomplish this objective, we initially compared three rainfall products based on WRF-driven rainfall and gauged rainfall. Then, we analyzed the characteristics of PE_PM. Finally, to comprehensively assess the WRF-driven Grid-XAJ model’s performance, we conducted a comparative analysis of flood simulations in both the humid Tunxi catchment and the semi-humid Chenhe catchment.

2. Methodology

2.1. The WRF Model

The governing equations of the WRF model encompass principles such as energy conservation and momentum conservation. The model parameterization schemes include parameterizations of microphysics and Cumulus parameterizations, among others [39]. The WRF modeling system mainly consists of the WRF Preprocessing System (WPS), the Advanced Research WRF (ARW) solver, post-processing, visualization tools, and the WRF Data Assimilation (WRF-DA) system. The functions of the WPS involve defining the simulation domains, interpolating terrestrial data to the simulation domain, and processing meteorological data to this simulation domain. The ARW solver is the key component of the modeling system, and comprises several initialization programs for idealized and real-data simulations, and the numerical integration program. The post-processing and visualization tools are responsible for analyzing the output data. The WRF-DA can be utilized to update the initial conditions in cycling mode when the WRF model is running. It can use observational data to estimate the initial fields as accurately as possible to enhance the forecast quality [40].

The WRF model (version 4.0) was utilized to generate the rainfall and PE for the Grid-XAJ model across the two study catchments, with each catchment covered by three-nested domains of 25 km, 5 km, and 1 km resolution. In the Tunxi catchment, the center of the three-nested domain was located at (29°58′ N, 117°58′ E), consisting of 111 × 91 grids for the inner domain, 121 × 111 grids for the middle domain, and 143 × 139 grids for the outer domain. For the Chenhe catchment, the center of the three-nested domain was located at (33°51′ N, 108°2′ E), with 91 × 61 grids for the inner domain, 141 × 121 grids for the middle domain, and 99 × 95 grids for the outer domain. The 6 h, 1° × 1° Final Operational Global Analysis data (FNL data, see Section 3.3), which were firstly processed by WPS to the simulation domains and then used in program “real.exe” of the WRF model, were used to generate boundary layer conditions (25 km × 25 km) for outer domain and initial conditions (25 km × 25 km, 5 km × 5 km, and 1 km × 1 km) for all three domains in this study.

Table 1. Physical parameterizations of the WRF model.

Category	Parameterization Selected	Reference
Microphysics option	Thompson	Thompson [41]
Cumulus option	Kain–Fritsch	Kain and Fritsch [42]
Planetary boundary layer	YSU	Hong and Pan [43]
Radiation physics option	RRTMG	Shin [44]

presents the main physical parameterizations of the WRF model.

2.2. The Grid-XAJ Model

The Grid-XAJ model [10] assumes that each DEM grid within the catchment is a computational unit, assuming the uniform spatial distribution of rainfall and the underlying conditions within each grid. It only considers the variability of various elements between different grids. Firstly, the model calculates the precipitation intercepted by the canopy (Icum). Secondly, the channel precipitation (Ich) that directly falls into the channel is treated as direct runoff. Thirdly, an evapotranspiration (ET) calculation is performed using a three-layer soil moisture module. Fourthly, the runoff (R) generation in each grid is simulated based on the saturation-excess runoff. Then, R is divided into the surface runoff (RS), interflow (RI), and groundwater runoff (RG) using the free-water sluice reservoir method. Finally, flow routing is accomplished using a grid-based Muskingum method. The static data required for the Grid-XAJ model include the DEM and slope direction, among others. Figure 1 shows the structure of the Grid-XAJ model.

We employed the trial-and-error method to optimize the model parameters. In detail, different parameter combinations were set in the Grid-XAJ model to simulate floods. The parameters are global optimal when evaluation metrics reach the global optimal values.

In this study, gauged rainfall was interpolated to each grid using the inverse distance weighting (IDW) method (Ir). The rainfall driving the Grid-XAJ model contains WRF-driven rainfall (Wr) and WRF-merged rainfall (Wm), which were obtained by merging Wr with gauged rainfall. PE for the Grid-XAJ model was calculated using the Penman–Monteith equation (PE_PM). The input data for the Penman–Monteith equation, including radiation, temperature, pressure, humidity, and wind speed, were acquired from the WRF simulation.

2.3. Successive Correction Method

The successive correction method [45] is employed as a technique to refine forecast fields by incorporating gauged data. This approach can be traced back to the objective analysis method, initially proposed by Cressman in 1959. In this research, we utilized the successive correction method to merge the WRF-driven rainfall with gauged rainfall. The calculation steps are as follows.

The initial step of the successive correction method was to determine the initial field, for which we used the WRF-driven rainfall in our research. Subsequently, selecting an appropriate weight function becomes crucial. Notably, a variety of weight functions are offered to the successive correction method. We adopted Cressman’s circular scheme [45] as the weight function in this paper.

$$W_{ij,k}=\left\{\begin{array}{cc}\frac{R^2-{d_{ij,k}}^2}{R^2+{d_{ij,k}}^2}\hfill & d_{ij,k}<R\hfill \\ 0\hfill & d_{ij,k}\ge R\hfill \end{array}\right.,$$

(1)

where W_ij,k is the weight factor of grid (i, j) for rain station k; R is the influence radius, and is defined as the minimum value between the number of rows and columns within the WRF domain for each catchment; and d_ij,k is the distance from grid (i, j) to rainfall station k.

The subsequent step was acquiring the correction field, which was computed based on the gauged data and the forecasted data at the grid:

$$C_{ij}=\frac{1}{n}{\sum }_{k=1}^n{W}_{ij,k}{E}_k,$$

(2)

$$E_k=O_k-G_k,$$

(3)

where C is the correction field, n is the total number of stations, E_k is the correction value of station k, O_k are the gauged data in station k, and G_k is the initial value in station k.

Finally, we obtained the analysis field based on the initial field and the correction field. Sometimes multiple iterations are necessary. In such instances, the current analysis field serves as the initial field for subsequent iteration. This study involved a total of five iterative corrections.

$$A_{ij}^m=A_{ij}^{m-1}+C_{ij},$$

(4)

where A is the analysis field, and m is the number of iterative corrections, m = 5.

2.4. The Penman–Monteith Equation

Penman first established an evaporation model based on an energy balance equation [46]. Then, Monteith enhanced the model by using aero-dynamic conductance, which depends not only on wind speed, but also on surface ruggedness, atmospheric stability, and vegetation height [47]. This modified model is known as the Penman–Monteith model. In 1998, the experts and researchers of the Food and Agriculture Organization (FAO) updated and simplified the Penman–Monteith model by using some assumed constant parameters, such as a grass reference crop of height 0.12 m, a surface resistance of 70 s/m, and an albedo value of 0.23 [48].

The new equation for hourly potential evapotranspiration [49] is as follows:

$${PE}_{PM}=\frac{0.408\Delta \left(R_n-G\right)+\gamma \frac{37.0}{T_{hr}+273}{u}_2(e_{T_{hr}}^0-e_a)}{\Delta +\gamma (1+0.34u_2)},$$

(5)

where PE_PM is the hourly potential evapotranspiration (mm/h); Rn is the net radiation at the crop surface (MJ/(m h)); G is the soil heat flux density (MJ/(m h)); T_hr is the mean hourly air temperature at a 2 m height (°C); u₂ is the wind speed at a 2 m height (m/s); $e_{T_{hr}}^0$ is the saturation vapor pressure (kPa); e_a is the actual vapor pressure (kPa); Δ is the slope of the vapor pressure curve (kPa/°C); and γ is the psychrometric constant, whose value is 0.066 kPa/°C. The parameters of the above equation are described below [49]:

$$R_n=\left(1-\alpha \right)RS+\epsilon (RL-\sigma T_s^4),$$

(6)

where α is the surface albedo; ε is the surface emissivity; RS is the incoming short-wave radiation, W/m²; RL is the incoming long-wave radiation, W/m²; σ is the Stefan–Boltzmann constant, whose value is 5.67 × 10⁻⁸ W/(m² K⁴); and Ts is the land surface temperature, K.

G is estimated as follows:

$$G=\left\{\begin{array}{cc}0.1R_n& \mathrm{daytime}(8:00~20:00)\\ 0.5R_n\hfill & \mathrm{nighttime}(20:00~8:00)\hfill \end{array}\right.,$$

(7)

The WRF provides the wind speed (u₁₀) at height h = 10 m, which can be converted to a 2 m wind speed:

$$u_2=u_{10}\frac{4.87}{\ln (67.8h-5.42)},$$

(8)

$e_{T_{hr}}^0$ is estimated as follows:

$$e_{T_{hr} }^0=0.6108\exp \frac{17.27T_{hr}}{T_{hr}+237.3},$$

(9)

e_a is estimated as follows [50]:

$$q=0.622\frac{e_a}{p-0.378e_a},$$

(10)

where q is the humidity, g/g, and p is the air pressure at 2 m, kPa.

Δ is estimated as follows:

$$\Delta =\frac{4098(0.6108\exp (\frac{17.27T_{hr}}{T_{hr}+273}))}{{(T_{hr}+273)}^2},$$

(11)

2.5. Evaluation Methods and Metrics

The accuracy of rainfall, evapotranspiration, and streamflow is characterized using five assessment metrics: the correlation coefficient (RR) [51], the Nash–Sutcliffe efficiency coefficient (NSE) [52], the percent bias (PB) [53], Shannon entropy (SE) [54], and the time error of flow peak (TEP).

$$\mathrm{RR}=\frac{{\sum }_{i=1}^n\left(S_i- \overline{S}\left)\right(O_i- \overline{O}\right)}{\sqrt{{\sum }_{i=1}^n{(S_i-\overline{S})}^2\sqrt{{\sum }_{i=1}^n{(O_i -\overline{O})}^2}}}\overline{S}=\frac{1}{n}{\sum }_{i=1}^n{S}_i,$$

(12)

$$\mathrm{NSE}=1 - \frac{{\sum }_{i=1}^n{\left(S_i - \overline{S}\right)}^2}{{\sum }_{i=1}^n{\left(O_i -\overline{O}\right)}^2}\overline{o}=\frac{1}{n}{\sum }_{i=1}^n{O}_i,$$

(13)

$$\mathrm{PB}=\frac{{\sum }_{i=1}^n(S_i -O_i)}{{\sum }_{i=1}^n{O}_i}\times 100\%,$$

(14)

$$\mathrm{SE}\left(x\right)= -{\sum }_{k=1}^{n}p\left(x_k\right){\log }_2\left(p\left(x_k\right)\right)x={(x_1,x_2,\dots ,x_N)}^T,$$

(15)

where S_i are the simulated results for each time step i, O_i is the observed value, n is the total number of time series, x is a discrete random variable recording all different values in the grids, p(x_k) is the frequency of the value of x_k appearing in the grids, and N is the length of the x. The RR and NSE quantify the degree-of-fit between simulation and observation. PB measures the errors of the simulated results. In the case of single-value (such as rainfall peak) simulation output, PB represents the relative error between this value and the gauged value. The SE is regarded as a measure for spatial variability, and the higher its values are, the more complex the related spatial information is. The TEP is used to evaluate the time accuracy of the simulated flow peak.

The rainfall data were characterized in terms of three aspects: rainfall series, cumulative rainfall, and rainfall spatial distribution. Evaluation metrics including RR, NSE, and the PB of the rainfall peak (PB_pr) were selected for analyzing the rainfall series [25]. For assessing the cumulative rainfall, the PB of the cumulative rainfall (PB_cr) was selected as an evaluation metric. Due to the distinct spatial distributions of the three types of rainfall, SE was selected to evaluate their spatial characteristics [25]. For PE_PM, we mainly focused on its spatial and temporal variations and possible influencing factors. The PB, NSE, PB of the flow peak (PB_pf), and TEP were used to evaluate the flood simulations.

According to the Standard for Hydrological Information and Hydrological Forecasting (GB/T22482-2008) [55], if the PB ≤ 20%, the flood volume simulation is considered acceptable. Similarly, if the PB_pf ≤ 20%, the flow peak simulation is considered acceptable. Additionally, the TEP within a three-hour window is deemed acceptable for peak occurrence time forecasting. The qualification rate (QR) of flood simulations can be determined as follows:

$$\mathrm{QR}=\frac{n}{m}\times 100\%$$

(16)

where n are the qualified events, and m are the total simulated events.

3. Study Area and Data

3.1. Study Catchments

The Tunxi catchment (29°27′ N~30°06′ N, 117°38′ E~118°30′ E) is located in the southeastern part of Anhui Province, with an area of 2670 km². It has an elevation ranging from 92 m to 1622 m. The catchment is situated within a subtropical monsoon climate zone and experiences an average annual rainfall of approximately 1600 mm, making it a representative humid catchment. Rainfall is highly unevenly distributed throughout the year, with the months from April to June being rainy, which often leads to flood disasters, while the period from July to September experiences frequent drought. Farmland and woodland collectively occupy over 60% of the Tunxi catchment, with clay loam being the predominant soil type.

The Chenhe catchment (33°40′ N~34°01′ N, 107°42′ E~108°20′ E) is located in the central-southern part of Shaanxi Province, with an area of approximately 1380 km² and elevations ranging from 600 m to 3800 m. This catchment represents a typical semi-humid catchment, with an annual runoff depth ranging from 100 to 500 mm and a multi-year average rainfall of 780 mm. The land use in the Chenhe catchment primarily consists of grassland, with a minor proportion of mixed forest and deciduous broadleaf forest. The predominant soil type is loam. The distribution of rivers and stations in the two catchments is shown in Figure 2.

3.2. Gauge Data

The hydrometeorological data for the Tunxi catchment, including daily evaporation, hourly discharge, and hourly precipitation, were provided by the Hydrology Bureau of Anhui Province. The catchment is equipped with eight rainfall stations and the Tunxi flow station. The Tunxi flow station, located at the catchment outlet, provided both rainfall and discharge data, while the rainfall stations only provided rainfall data. The evaporation data were obtained from the evaporation pan at the Tunxi flow station. We used twelve flood events for parameter calibration and six flood events for validation between 2012 and 2017. Then, four flood events were chosen to show the flood simulation performance of the WRF-driven Grid-XAJ model.

The Hydrology Bureau of Shaanxi Province provided the daily evaporation, hourly discharge, and 6 h precipitation data for the Chenhe catchment. Within this catchment, there are also eight rainfall stations and the Chenhe flow station. Some rainfall data and discharge data were sourced from the Tunxi flow station, and the rainfall stations only provided rainfall data. Evaporation was obtained via the evaporation pan at the Chenhe flow station. To determine the optimal model parameters in the Chenhe catchment, we utilized eight flood events for parameter calibration and four flood events for validation between 2003 and 2012. Similar to the Tunxi case study, we selected four flood events to show the flood simulation performance of the WRF-driven Grid-XAJ model. The initial state of the Grid-XAJ model was acquired after a warm-up period of a minimum of four months of daily flood simulation.

Table 2. The set of flood events.

Catchment	Year	Events	Simulation Period	Peak Discharge (m3/s)
Tunxi	2012	120423	23 April, 14:00–27 Aapril, 06:00	3170
	2013	130606	6 June, 14:00–11 June, 20:00	3610
	2015	150607	7 June, 08:00–10 June, 17:00	3010
	2017	170623	23 June, 08:00–28 June, 00:00	4210
Chenhe	2003	030903	3 September, 02:00–8 September, 20:00	740
	2005	050928	28 September, 08:00–3 October, 20:00	1740
	2011	110916	16 September, 14:00–21 September, 08:00	1200
	2012	120830	30 August, 20:00–3 September, 18:00	1710

lists the selected events in the Tunxi and Chenhe catchments for driving the WRF-driven Grid-XAJ model.

3.3. FNL Data

The FNL(Final) Operational Global Analysis data provided by the National Centers for Environmental Prediction (NCEP) boast a high resolution and incorporate various observational and satellite-derived data. The FNL data are collected four times daily, the observations spanning the 6 h prior at least, and undergo global data analysis. They are widely utilized in numerical simulations and a range of weather analysis research [56,57]. For this study, FNL data from 2003 to 2017 were selected to provide the initial conditions for running the WRF model, with the FNL data being processed by the WPS.

4. Results and Discussion

4.1. Parameters Calibration

As explained before, a priori estimates of most model parameters were obtained from the underlying surface data. Then, the trial-and-error method was used to update the values of the parameters. The parameters that needed to be calibrated included K, C, Ci, Cg, Cs, and Lag, whose detailed description can be seen in reference [10]. For the Tunxi catchment, when all the flood evaluation metrics reached the global optimal values, the mean NSE for the eighteen flood events was 0.93 and the mean PB was 0.04. The QR of the flow peak was 94% and the QR of the peak occurrence time was 100%. For the Chenhe catchment, the mean NSE for the twelve flood events was 0.77 and the mean PB was 0.11. The QR of the flow peak was 94% and the QR of the peak occurrence time was 100%. The optimal values of K, C, Ci, Cg, Cs, and Lag are shown in

Table 3. The optimal values of the parameters in two catchments.

Parameter	Description	Optimal Estimate of the Tunxi Catchment	Optimal Estimate of the Chenhe Catchment
K	Ratio of potential evapotranspiration to pan evaporation	0.98	0.7
C	Evapotranspiration coefficient of deeper layer	0.18	0.08
Ci	Recession constant of interflow storage	0.3	0.55
Cg	Recession constant of groundwater storage	0.98	0.87
Cs	Recession constant in the lag and route technique	0.93	0.89
Lag	Lag time	1.0	2.0

. In addition, the simulation performance of the Grid-XAJ model in the Tunxi catchment was obviously better than that in the Chenhe catchment.

4.2. Evaluation of Two Rainfall Products

Evaluating the Wr and Wm is crucial when analyzing the WRF-driven Grid-XAJ model. In this study, we regarded the Ir as the observed areal precipitation. The analysis focused on the rainfall series, cumulative rainfall, and rainfall spatial distribution.

Table 4. The values of the four evaluation metrics of Wr and Wm.

Catchment	Events	Wr				Wm
		RR	NSE	PBpr%	PBcr%	RR	NSE	PBpr%	PBcr%
Tunxi	120423	0.66	0.32	−40.1	−60.1	0.99	0.98	−2.6	−0.4
	130606	0.90	0.80	−10.3	−14.3	0.99	0.99	3.5	−2.0
	150607	0.73	0.18	33.9	−4.6	0.99	0.98	−4.8	3.3
	170623	0.69	0.36	−35.8	13.1	0.99	0.99	−8.1	6.7
Chenhe	030903	0.39	−1.29	142.6	51.2	0.95	0.75	88.1	36.8
	050928	0.18	−1.62	83.7	15.9	0.98	0.91	10.3	15.8
	110916	0.45	0.18	−58.0	−16.0	0.99	0.96	17.3	13.8
	120830	0.70	0.42	−47.8	−38.3	0.99	0.92	23.0	21.1

displays the four evaluation metrics of the Wr and Wm. Figure 3 illustrates the hourly rainfall and cumulative rainfall. From Figure 3 and Table 4, Wr adequately simulates the real rainfall series but exhibits limitations in terms of the rainfall peak, cumulative rainfall, and degree-of-fit. Conversely, Wm significantly reduces the simulation errors, achieving a level of accuracy comparable to Ir.

For the Tunxi catchment, the WRF model provided a reasonably accurate simulation of the rainfall series, with Wr showing a strong correlation with Ir (RR > 0.65). However, the degree-of-fit between Wr and Ir was unsatisfactory, as only event 130606 achieved a NSE value of 0.80, while the other three events achieved NSE values below 0.40. Notably, the WRF model successfully simulated all rainfall peaks, with all PB_pr values being below 50%. Event 130606 exhibited a relatively large PB_cr value of 60.1%, while the PB_cr values for the other three events remained within 20%. Merging Wr with the gauged rainfall significantly improved rainfall accuracy, as indicated by Wm achieving higher RR values (>0.99). Furthermore, Wm demonstrated notable improvements in the degree-of-fit with Ir, exceeding all NSE values above 0.98. Additionally, the PB_pr of Wm was substantially reduced to below 10%. In terms of cumulative rainfall analysis, there was a significant decrease in PB_cr values across all four events, with values falling below 10%.

For the Chenhe catchment, the accuracy of Wr was relatively poor. Although Wr captures the rainfall series, its correlation with Ir was low, with only event 120830 achieving a RR value of 0.70. However, the RR values for the other three events were all below 0.50. Additionally, Wr exhibited a poor degree-of-fit with Ir, as indicated by NSE values below 0.50 for all four events. Moreover, the WRF model failed to accurately simulate rainfall peaks, for it generated non-existent peaks for events 030903 and 050928, and only partially simulated peaks for events 110916 and 120830. The PB_pr of Wr for all four events exceeded 40%. Furthermore, the PB_cr values were above 30% for events 030903 and 120830. After merging Wr and gauged rainfall, Wm achieved RR values of at least 0.95, with the NSE values all being above 0.75. The PB_pr of Wm was noticeably reduced compared to that of Wr, and the PB_cr values were also significantly decreased (<37%). In addition, the cumulative rainfall of Wm was consistently higher than that of Ir.

We further compared the rainfall of the Tunxi catchment and the Chenhe catchment. As shown in Table 4, the rainfall series, rainfall peaks, and cumulative rainfall in the Tunxi catchment outperformed the results in the Chenhe catchment. In the Tunxi catchment, the values of mean RR, mean NSE, mean PB_pr, and mean PB_cr were 0.75, 0.42, 30.0%, and 23.0% for Wr; and 0.99, 0.99, 4.8%, and 3.1% for Wm. In the Chenhe catchment, they were 0.43, −0.58, 83.0%, and 30.4% for Wr; and 0.98, 0.89, 34.7%, and 21.9% for Wm.

Figure 4 and Figure 5 present spatial distribution maps of the maximum cumulative rainfall in the two catchments. They demonstrate that Wm combines the advantages of both Ir and Wr. (1) The cumulative rainfall of Wm falls within most of the ranges of Ir and Wr. (2) The rainfall center of Ir is situated at the rain stations due to limitations imposed by the density and quality of gauged data. Wr can capture the rainfall center and spatial distribution details that rain stations fail to capture, but its spatial distribution and rainfall volume are not very accurate. Wm can incorporate rainfall centers from the other two rainfall products. (3) In both catchments, the spatial distribution of Ir is relatively smooth, and is closely related to the distribution of rainfall stations, resulting in simple spatial information. Wr exhibits a complex multi-banded spatial distribution. Wm also displays a multi-banded spatial distribution, but provides more detailed information compared to Ir, leading to a higher SE than that of Ir (0.3%~49.2%). When SE(Ir) < SE(Wr), merging Wr with Ir will simplify the spatial complexity of Wr. Conversely, when SE(Ir) > SE(Wr), merging Wr with Ir will enhance the spatial information.

4.3. Characteristics of PE_PM

Due to the lack of observed areal PE, we made a concise analysis comparing PE_PM and PE_pan. Then, we mainly focused on some influence factors for PE_PM and its spatiotemporal distribution characteristics.

Figure 6 presents a comparison between the PE_PM and PE_pan series. As shown in the figure, PE_pan remains stable throughout the day, with sudden changes occurring only at the transition to the next day. Contrarily, PE_PM exhibits hourly variations and changes gradually over time. In addition, the range of PE_PM (0.0 mm~0.5 mm for the Tunxi catchment, 0.0 mm~0.4 mm for the Chenhe catchment) is greater than that of PE_pan (0.0 mm~0.3 mm for the Tunxi catchment, 0.0 mm~0.25 mm for the Chenhe catchment).

As shown in Figure 6, rainfall significantly influences PE_PM. PE_PM notably decreases during rainy periods in both catchments. Moreover, the influence of different climatic zones is evident in the PE_PM values. The PE_PM values in the Tunxi catchment are considerably higher compared to those in the Chenhe catchment. Specifically, in the Tunxi catchment, the PE_PM values range from 0.0 to 0.5 mm. In the Chenhe catchment, the PE_PM values range from 0.0 to 0.2 mm for events 050928 and 110916, and from 0.0 to 0.4 mm for the other two events.

PE_PM exhibits distinct temporal distribution characteristics. As illustrated in Figure 6, evapotranspiration mainly occurs between 6:00 a.m. and 6:00 p.m., with the peak PE_PM occurring around noon at 12:00 p.m. During the nighttime, the PE_PM values approach zero.

PE_PM exhibits distinct spatial distribution characteristics. Figure 7 depicts the spatial distribution of the mean PE_PM, where the mean PE_PM refers to the average spatial distribution of PE_PM over all time periods. The SE values exceed 2.0 for the Tunxi catchment and surpass 3.0 for the Chenhe catchment. Figure 8 illustrates the relationship between the mean PE_PM and elevation. As shown in Figure 8, there is a relationship between PE_PM and elevation: regions with higher elevations display relatively lower mean PE_PM values, while areas with moderate elevations exhibit the highest mean PE_PM values, and slightly lower values are observed in lower elevation areas. Here are the specifics: (1) In the lower elevation areas of the catchment, the mean PE_PM is moderately high. Within the 0~400 m elevation range in the Tunxi catchment, except for event 170623, there is a slight increase in the mean PE_PM with higher elevations. Similarly, within the 600~800 m elevation range in the Chenhe catchment, except for event 120830, there is also a slight increase in the mean PE_PM with higher elevations. Although the mean PE_PM values of event 170623 in the Tunxi catchment and event 120830 in the Chenhe catchment decreased with increasing elevation, it is worth noting that some lower elevation areas, such as rivers near the catchment outlets, exhibit a relatively low mean PE_PM. (2) In the middle elevation areas of the catchment, the mean PE_PM is the highest. Within the 400~1200 m elevation range of the Tunxi catchment, there is a gradual decrease in the mean PE_PM as the elevation increases. Similarly, within the 800~2400 m elevation range of the Chenhe catchment, although there is a fluctuating trend in the mean PE_PM with increasing elevation for event 110916, for the other three events, there is a gradual decrease in the mean PE_PM with higher elevations. (3) In the high elevation areas of the catchment, the mean PE_PM is the lowest. In the regions with an elevation above 2400 m in the Chenhe catchment and above 1200 m elevation in the Tunxi catchment, the mean PE_PM exhibits a fluctuating pattern. For all events, as the mean PE_PM decreases to a certain extent with increasing elevation, it then shows an upward trend. The possible reason for this phenomenon is that PE is mainly affected by wind speed in the middle and low elevations of the catchment, while PE is mainly affected by temperature in the high elevations of the catchment.

4.4. Evaluation of Discharge

Different from the traditional Grid-AXJ model, Wr and Wm were used as input rainfall to drive the Grid-XAJ model for flood forecasting in this paper. Furthermore, we considered substituting PE_pan with PE_PM in the model. The results are presented in Figure 9.

Table 5. The evaluation metrics of the simulated floods.

Catchment	Event	Wr + PEPM				Wm + PEPM
		NSE	PB%	PBpf%	TEP (Hour)	NSE%	PB%	PBpf%	TEP (Hour)
Tunxi	120423	0.12	−67.26	−77.55	−3	0.94	3.85	−7.60	0
	130606	0.93	−5.42	−25.38	−2	0.91	4.18	−16.24	1
	150607	0.82	−1.51	36.41	−2	0.94	1.71	14.51	−2
	170623	0.68	11.96	31.36	0	0.96	4.14	−5.59	0
Chenhe	030903	0.46	31.97	−1.20	−9	0.66	23.71	21.15	4
	050928	−1.65	24.48	−1.49	−42	0.54	26.08	−9.06	−3
	110916	0.54	−16.43	−29.75	20	0.87	19.20	16.27	0
	120830	0.54	−42.58	−48.40	−5	0.66	41.69	7.10	−2

shows the values of the metrics used for the simulated floods. It is evident from Figure 9 and Table 5 that directly using Wr and PE_PM to simulate floods yields unsatisfactory results. However, incorporating Wm as input rainfall significantly enhances the forecasting performance.

To evaluate the flood forecasting capability of the WRF-driven Grid-XAJ model, Wr + PE_PM and Wm + PE_PM forcing scenarios were designed for the Grid-XAJ model. In the Tunxi catchment, for scenario Wr + PE_PM, only one flood event exhibits an NSE value exceeding 0.90. The PB values are relatively small for events 130606 and 150607. However, for the other two events, particularly event 120423, the PB values are significantly larger. All simulated flow peaks fail to meet the qualification criteria, with all PB_pf values exceeding 20%, despite having a small TEP. For scenario Wm + PE_PM, all NSE values surpass 0.90. The flood volume simulation is robust, with all PB values below 5%. For flow peaks, the PB_pf values for all four events are below 20%, with events 130606 and 170623 achieving PB_pf values below 10%. The TEP for the four floods remain within a 3 h margin, demonstrating significant reductions compared to scenario Wr + PE_PM.

In the Chenhe catchment, for scenario Wr + PE_PM, the NSE values for all events remain below 0.55, indicating a suboptimal degree-of-fit between the simulation and observation. Except for event 110916, all PB values exceed 20%. Additionally, notable deficiencies are observed in peak flows. Although the PB_pf values of events 030903 and 050928 meet the criteria, significant temporal discrepancies exist for the TEP (−9 h, −42 h). For events 110916 and 120830, neither the flow peaks nor the peak occurrence times align with the gauged discharges. For scenario Wr + PE_PM, all four events achieve NSE values surpassing 0.5. Although event 030903 does not meet the criteria for PB_pf and TEP, these deviations are relatively small. The other three events yield satisfactory results in terms of PB_pf and TEP.

To further illustrate the impact of merging the method of rainfall on flood forecasting, we compared the streamflow simulations of the Grid-XAJ model under two forcing scenarios. It is clear that scenario Wm + PE_PM significantly improves the degree-of-fit, flood volume, and flow peaks over scenario Wr + PE_PM. In the Tunxi catchment, compared to scenario Wr + PE_PM, the mean NSE value of scenario Wm + PE_PM increased by 0.3, the mean PB value decreased from 21.5% to 3.5%, the QR of PB_pf increased from 0 to 100%, and the TEP decreased one hour. In the Chenhe catchment, compared to scenario Wr + PE_PM, the NSE value of scenario Wm + PE_PM increased from −1.65~0.54 to 0.54~0.87. The PB values exhibit minor fluctuations, with reductions of 8.26% and 0.89% for events 030903 and 120830, and a slight increase of 1.6% and 2.77% for events 050928 and 110916. The QR of PB_pf increased from 50% to 100% and the TEP decreased by seventeen hours.

We further compared the flood simulations of the WRF-driven Grid-XAJ model in different climatic catchments. For scenario Wr + PE_PM, the mean NSE value in the Tunxi catchment is 66.4% higher than that in the Chenhe catchment, and the mean PB is 7.3% lower than that in the Chenhe catchment. Additionally, the QR of the TEP is 100% in the Tunxi catchment but 0% in the Chenhe catchment. For scenario Wr + PE_PM, all NSE values exceed 0.90 in the Tunxi catchment, whereas in the Chenhe catchment, all NSE values fall below 0.90. In the Tunxi catchment, the QRs for PB, PB_pf, and TEP are all 100%, while in the Chenhe catchment, the QRs for the PB, PB_pf, and TEP are 50%, 75%, and 75%, respectively. These outcomes indicate that the results of the model in the semi-humid catchment can only serve as a reference.

4.5. Discussion

We built a WRF-driven Grid-XAJ model, fully considering the influence of meteorology to enhance the accuracy of flood simulations. The traditional forcing scenario of spatial interpolation rainfall and pan evaporation renders the Grid-XAJ model inapplicable in regions with limited rainfall data and neglects other meteorological variables. Additionally, their spatiotemporal accuracy may be affected by many factors. In this study, we fully considered the influence of meteorology on the Grid-XAJ model compared with Yao [24] and Meng [38], though the forcing scenarios in their studies do improve the accuracy of flood forecasting. Using the successive correction method actually enhanced the spatiotemporal accuracy of Wr distribution (Table 4), which was consistent with other studies that used this method to improve rainfall accuracy [58]. PE_PM, which exhibited a reasonable spatiotemporal distribution (Figure 7 and Figure 8), was also consistent with other research that used it with NWM output for PE calculation [59].

Table 4 shows that Wr and Wm performed better in the Tunxi catchment than in the Chenhe catchment. This can be attributed to the fact that the humid catchment typically experiences more intense rainfall with shorter durations, and the WRF model performs better in short-term heavy rainfall forecasting rather than long-term light rainfall. Additionally, errors in the WRF-driven rainfall can partly reflect the performance of the WRF-driven Grid-XAJ model (Figure 3 and Figure 9). The NSE values of Wm are higher than those of Wr, suggesting that Wm may exhibit better performance in degree-of-fit for simulated floods. Errors in the peak rainfall can lead to inaccurate predictions of the peak flood. The application of Wm enhances the PB_pr of Wr, which positively influences the simulation of flow peaks. An accurate estimate for the cumulative rainfall will improve the precision of simulated flood volume. The error in the cumulative rainfall of Wm is smaller than that of Wr, which indicates that Wm is more suitable to simulate flood volume than Wr. In addition, the error of each metric for Wr and Wm in the Tunxi catchment is lower than that in the Chenhe catchment, suggesting that the Grid-XAJ model may yield more accurate flood simulations in the Tunxi catchment. Figure 6 shows that PE_PM exhibits a correlation with rainfall. This association can be attributed to the concurrent decrease in air temperature and increase in moisture content during rainy periods. Furthermore, PE_PM in the Tunxi catchment was higher compared to that in the Chenhe catchment (Figure 6). This may be attributed to the elevated soil moisture content and a greater availability of water vapor within the humid catchment. Research has shown a strong correlation between PE_PM and elevation [60], and that was proved in our study (Figure 7).

As mentioned in Section 2.2, the Grid-XAJ model adopts the saturation-excess runoff module to calculate runoff generation, which represents the main runoff generation mode in humid catchments [61]. The runoff generation mode of a semi-humid catchment is intricate and potentially encompasses various other modes of runoff generation. Thus, the Grid-XAJ model appears to be better suited for a humid watershed.

5. Conclusions

Based on gauged rainfall, WRF-driven rainfall (Wr), and meteorological data output by WRF, WRF-merged rainfall (Wm) and potential evapotranspiration calculated using the Penman–Monteith equation (PE_PM) were utilized in this paper. Then, two forcing scenarios (Wr + PE_PM and Wm + PE_PM) were designed to drive the Grid-XAJ model in both the Tunxi and the Chenhe catchments, aiming to enhance the accuracy of flood simulations. The main conclusions are as follows.

The accuracy of Wm surpasses that of Wr in both humid and semi-humid catchments, indicating a more rational spatial distribution than Ir. Wm also incorporates rainfall centers from both Wr and Ir. PE_PM demonstrates significant spatiotemporal heterogeneity, exhibiting hourly variations compared to PE_pan. Additionally, PE_PM can be influenced by rainfall, and its values are higher in the humid catchment than in the semi-humid catchment. The spatial distribution of PE_PM exhibits a gradual decrease from moderately elevated regions towards areas characterized by both higher and lower elevations.

The WRF-driven Grid-XAJ model holds significant potential in flood forecasting. The Grid-XAJ model, when driven by scenario Wr + PE_PM, shows a broad capacity to simulate floods. In contrast, under scenario Wm + PE_PM, the model demonstrates a superior flood forecasting performance. Additionally, the WRF-driven Grid-XAJ model demonstrates superior flood forecasting performance in the humid catchment compared to the semi-humid catchment.

6. Study Limitations and Prospects

In this study, the rainfall inputs into the Grid-XAJ model included WRF-driven rainfall and WRF-merged rainfall, which are based on FNL data and gauged data, thereby not extending the lead time of flood forecasting. In addition, we examined the influence of meteorological factors on hydrological processes, but the feedback of hydrological factors on meteorological processes was not considered. In future studies, we will try to extend the lead time by using forecasted data (i.e., GFS data) to drive the WRF model and improve the accuracy of the WRF-forecast rainfall before it is used to drive the Grid-XAJ model. Furthermore, we will explore methods to realize the two-way coupling between the WRF model and the Grid-XAJ model.