A Coupled Hydrologic–Hydraulic Model (XAJ–HiPIMS) for Flood Simulation

: To protect ecologies and the environment by preventing ﬂoods, analysis of the impact of climate change on water requires a tool capable of considering the rainfall-runo ﬀ processes on a small scale, for example, 10 m. As has been shown previously, hydrologic models are good at simulating rainfall-runo ﬀ processes on a large scale, e.g., over several hundred km 2 , while hydraulic models are more advantageous for applications on smaller scales. In order to take advantages of these two types of models, this paper coupled a hydrologic model, the Xinanjing model (XAJ), with a hydraulic model, the Graphics Processing Unit (GPU)-accelerated high-performance integrated hydraulic modelling system (HiPIMS). The study was completed in the Misai basin (797 km 2 ), located in Zhejiang Province, China. The coupled XAJ–HiPIMS model was validated against observed ﬂood events. The simulated results agree well with the data observed at the basin outlet. The study proves that a coupled hydrologic and hydraulic model is capable of providing ﬂood information on a small scale for a large basin and shows the potential of the research.


Introduction
Flood, as the most common natural hazard, causes huge causalities, heavy economic losses, and ecological and environmental risks worldwide [1,2]. A warmer climate creates more extreme storm events, combined with rapid land use/cover change and a growing population, the risk of floods to urban and farming areas is increasing and posing a threat to ecologies and the environment [2][3][4][5][6]. Extreme rainfall is considered as the main cause of flooding [7], however, inadequate flood risk management also contributes to the increase of flood risk [8]. There is an urgent requirement for flood simulation/forecasting and integrated flood risk management in our changing environment.
Hydrologic models are commonly implemented for rainfall-runoff simulations to predict the basin-scale flow of water passing through the basin outlet, normally on a daily/hourly temporal scale, e.g., Soil and Water Assessment Tool (SWAT) [9,10], Xinanjing (XAJ) [11,12], Stanford Watershed Model (SWM) [13], Sacramento Model (SAC) [14]. However, high spatial heterogeneity of a realistic basin raises the requirement for fine flood management characterized by spatial information inside of the complex basin that affects its response to extreme storms. Numerous distributed grid-based hydrologic models have been developed to represent spatial variability [15][16][17]. River channel flow is computed by the Muskingum method, e.g., [15], or the Muskingum method is coupled with a one-dimensional (1D) hydraulic model to generate water depth along a river channel, e.g., [18]. The overland flow is barely considered as a fully dynamic flow, and is normally solved for by hydrologic concepts such as linear reservoir [19] or is assumed to be a two-dimensional (2D) kinematic/diffusive wave equation [16].
Hence, the flood characteristics (i.e., flood extent, water depth, and velocity) cannot be represented inside the two-dimensional domain.
Currently, the two-dimensional hydraulic model provides a good solution for flood routing processes for fine flood risk assessment and management [20][21][22]. The fully 2D hydraulic model is a robust numerical methodology to provide accurate simulation of runoff flow on a minutely temporal scale, based on high-resolution data sources and GPU-accelerated techniques [23][24][25]. However, runoff generation is still a big concern for the application of the fully 2D hydraulic model in the large-scale of a basin, especially for small and medium floods. The typical solution of infiltration is to adopt the Kostiakov equation (e.g., [26]), Green-Ampt model (e.g., [27]), Horton equation (e.g., [28,29]), or Richard model (e.g., [30]).
In recent years, how to jointly use the hydrologic and hydraulic models for flood risk assessment and reservoir operation has been discussed [31][32][33]. As has been shown previously, hydrologic models are good at simulating rainfall-runoff processes for a basin on a large scale, e.g., several hundred km 2 , while hydraulic models are more advantageous for applications on smaller scales. In order to take advantage of these two types of models, this study attempts to couple a hydrologic model with a hydraulic model.
In this study, the Misai river basin (797 km 2 ), located in Zhejiang Province, China, is taken as the study area. For the hydrologic model, XAJ was selected as it has been widely used in China for a long time [34][35][36]. For the hydraulic model, the high-performance integrated hydraulic modelling system (HiPIMS) was chosen because of its high-speed calculations and wide application [23][24][25]. This paper introduces how the XAJ-HiPIMS is coupled from XAJ and HiPIMS and validates the new model against typical floods. Section 2 outlines the materials and methods, Section 3 presents results and discussion, and Section 4 gives conclusions.

XAJ
XAJ is a lumped rainfall-runoff model developed by the authors of [11,12]. Based on the watershed saturation-excess runoff theory, evapotranspiration is calculated in three soil layers. The soil's moisture storage capacity distribution curve is utilized to provide a non-uniform distribution of storage capacity in the hydrologic unit. The total runoff is composed of surface runoff, interflow, and groundwater runoff. The distribution of these three components can be decided by using a free water capacity distribution curve. The surface runoff is routed by the unit hydrograph. The interflow and groundwater flow are routed by the linear reservoir method. Hence, four modules compose the XAJ: evapotranspiration, runoff generation, runoff sources partition, and runoff routing ( [11,12], Figure 1). The parameters of XAJ are described in Table 3. The output of XAJ is outflow discharge at the outlet section of the basin.

HiPIMS
The governing equations of HiPIMS are fully 2D shallow water equations based on the Cartesian uniform grid [38], as follows,

HiPIMS
The governing equations of HiPIMS are fully 2D shallow water equations based on the Cartesian uniform grid [38], as follows, In Equation 1, q stands for the vector of conserved flow variables; f and g are the flux vectors in the x-and y-directions, respectively; and R, Sb, and Sf denote the source terms of rainfall, bed slope, and friction, respectively. In Equation 2, h is the water depth (h = η − zb, where η is the water surface elevation and zb is the bed elevation); u and v are the depth-averaged velocity components in the xand y-directions, respectively; g is the acceleration of gravity; ρ is the water density; r is the generated surface runoff (which is the , in Equation 2); ∂b / ∂x and ∂b / ∂y denote the bed slope in the x-and y-directions; τbx and τby are the bed friction stresses, calculated by the following equations: where Cf ( = gn 2 / h 1/3 ) stands for the bed roughness coefficient; and n is the Manning's coefficient.
In this study, HiPIMS employs a first-order finite volume Godunov-type numerical scheme to solve the above governing equations. Equation 1 is discretized by Equation 4 and the flow variables can be updated as follows. In Equation (1), q stands for the vector of conserved flow variables; f and g are the flux vectors in the x-and y-directions, respectively; and R, S b , and S f denote the source terms of rainfall, bed slope, and friction, respectively. In Equation (2), h is the water depth (h = η − z b , where η is the water surface elevation and z b is the bed elevation); u and v are the depth-averaged velocity components in the xand y-directions, respectively; g is the acceleration of gravity; ρ is the water density; r is the generated surface runoff (which is the r n i, j in Equation (2)); ∂b / ∂x and ∂b / ∂y denote the bed slope in the x-and y-directions; τ bx and τ by are the bed friction stresses, calculated by the following equations: where C f ( = gn 2 / h 1/3 ) stands for the bed roughness coefficient; and n is the Manning's coefficient. In this study, HiPIMS employs a first-order finite volume Godunov-type numerical scheme to solve the above governing equations. Equation (1) is discretized by Equation (4) and the flow variables can be updated as follows.
where the superscript m denotes the present time level; ∆t is the present time step (which is the ∆t n in Equation (5)); Ω i stands for the area of cell i; F k (q) presents the fluxes that are normal to the cell edges; k is the index of the cell edges (k = 1-4); l k is the corresponding cell edge length; and n = (n x ,n y ) denotes the unit vector of the outward normal direction. In HiPIMS, the flux terms and bed slope term are calculated by an explicit scheme. The interface fluxes are calculated by the HLLC (Harten-Lax-van Leer-Contact) approximate Riemann solver. The local Riemann problems at the cell interfaces are solved by implementing the surface reconstruction method [39]. The friction term is solved by an implicit scheme [40], which is capable of maintaining the numerical stability when dealing with very small water depth. The adaptive time step method is implemented and controlled by the Courant-Friedrichs-Lewy (CFL) criterion. The high-performance of this CPU/GPU-integrating model is achieved by adopting the OpenCL programming framework [38,41].

Coupling Framework
In the coupling framework, XAJ generates the surface runoff and produces interflow and underground flow while the hydraulic model is run for 2D confluence calculations of the surface runoff generated by XAJ. The generated surface runoff was presented in the governing equations of HiPIMS as the source term in the mass conservation equation. The temporal scale of XAJ is hourly while the computing time-step of HiPIMS is in seconds, and the spatial scale of XAJ is the whole basin, while computing grid-cell of HiPIMS is in meters. In order to deal with this temporal and spatial downscaling issue, functions (Equations (5) and (6)) were created as below.
In Equation (5), R XAJ denotes the surface runoff calculated by the XAJ; ∆t n is the time-step of HiPIMS; n means the presently considered time-step; and r * i,j stands for the downscaled surface runoff. In Equation (6), r n i,j means the surface runoff in terms of the source term in the mass conservation equation of HiPIMS and m is the amount of the computing grid cell.
Then, the surface runoff generated by XAJ inputs to the hydraulic model to produce two-dimensional flow over the basin. Thus, the procedure is proposed as follows (shown in Figure 2

Statistical Method
In this work, the statistical methods absolute relative error of flood-peak discharge (ARED) (Equation 7), difference of peak arrival time (DPAT) (Equation 8), and the Nash-Sutcliffe efficiency coefficient (NSE) (Equations 9 and 10) are employed to estimate the discrepancy of the discharge. Three types of data are considered herein, i.e., observation, simulation by XAJ, and simulation by XAJ-HiPIMS.
where MAX(Qs) stands for the simulated peak discharge and MAX(Qo) denotes the observed peak

Statistical Method
In this work, the statistical methods absolute relative error of flood-peak discharge (ARED) (Equation (7)), difference of peak arrival time (DPAT) (Equation (8)), and the Nash-Sutcliffe efficiency coefficient (NSE) (Equations (9) and (10)) are employed to estimate the discrepancy of the discharge. Three types of data are considered herein, i.e., observation, simulation by XAJ, and simulation by XAJ-HiPIMS.
Water 2020, 12, 1288 where MAX(Q s ) stands for the simulated peak discharge and MAX(Q o ) denotes the observed peak discharge. ARED∈ (−∞,+∞), where T MAX(Qs) is the arrival time of simulated peak discharge by XAJ and XAJ-HiPIMS and T MAX(Qo) is the arrival time of the observed peak discharge.
where NSE 1 is the Nash-Sutcliffe efficiency coefficient computed by the XAJ simulation and observation; NSE 2 is the Nash-Sutcliffe efficiency coefficient computed by the XAJ-HiPIMS simulation and observation; Q o is the observed discharge; Q 1 is the simulated discharge of XAJ; Q 2 is the simulated discharge of XAJ-HiPIMS; and Q o is the averaged observed discharge.

Study Area
The Misai basin is located in the Zhejiang province of China with an area of 797 km 2 . The bed elevation differs from 105 to 1253 m. In this study, the 30-m ASTER GDEM V2 elevation data were adopted for computing, as shown in Figure 3. Six precipitation gauging stations are located in this basin. One is located at the outlet of the basin, and it also served as stream gauging station. The 30-m Landsat TM image data were used to provide the land use information. Considering the land characteristics, the land use types were classified into seven categories in this work, i.e., forest, heavy brush, cultivated land, grass land, pond/river, bare land, and urban land (distribution shown in Figure 4). The proportional area and Manning's coefficients of land use types are summarized in Table 1.

Study Area
The Misai basin is located in the Zhejiang province of China with an area of 797 km². The bed elevation differs from 105 to 1253 m. In this study, the 30-m ASTER GDEM V2 elevation data were adopted for computing, as shown in Figure 3. Six precipitation gauging stations are located in this basin. One is located at the outlet of the basin, and it also served as stream gauging station. The 30-m Landsat TM image data were used to provide the land use information. Considering the land characteristics, the land use types were classified into seven categories in this work, i.e., forest, heavy brush, cultivated land, grass land, pond/river, bare land, and urban land (distribution shown in Figure 4). The proportional area and Manning's coefficients of land use types are summarized in Table 1.

Modelling Set
XAJ: The rainfall input was given at 1-h intervals from the historic records from the six gauging stations in the basin. The 13 coefficients were calibrated and are shown in Table 3.
HiPIMS: The considered basin was divided into 909,824 grid cells upon a regular uniform grid of resolution 30 m. The Manning's n are given in Table 1. The CFL number was set to 0.35 in the present work.

Results and Discussion
In this study, the simulations of XAJ-HiPIMS and XAJ were assessed by comparing their results with the observed discharge at the outlet gauging station, as presented in Figure 5. Good agreement of the simulation with the observation was demonstrated for the three typical floods. Better performance was found in the bigger flood category. Further assessment has been done by implementing the statistical methods of ARED, DPAT, and NSE, shown in Table 4. We found that 1) in terms of the peak discharge, the simulation of XAJ-HiPIMS presented good agreement with the measurement; 2) in terms of ARED, the XAJ-HiPIMS achieved a better flood peak simulation, especially for the bigger flood; 3) in terms of DPAT, the Further assessment has been done by implementing the statistical methods of ARED, DPAT, and NSE, shown in Table 4. We found that (1) in terms of the peak discharge, the simulation of XAJ-HiPIMS presented good agreement with the measurement; (2) in terms of ARED, the XAJ-HiPIMS achieved a better flood peak simulation, especially for the bigger flood; (3) in terms of DPAT, the differences between simulations and measurements were no more than 1 h and smaller differences were found in the medium flood; (4) in terms of NSE, XAJ-HiPIMS obtained good flood simulations, which were greater than 0.85. The above results demonstrate XAJ-HiPIMS could provide satisfactory simulation at different flood scales. According to the 2D distribution of flood information, inundated extent and maximum water depth were investigated. The threshold of inundated depth was set to 0.01 m in this study. The snapshots of inundated extent were produced by XAJ-HiPIMS at the peak moment (the moment of peak discharge) shown in Figure 6 (left column). The value of water depth in each grid cell was recorded during the calculation. Maximum water depth is defined as the maximum value of water depth in each grid cell in the entire flood event record. The maps of maximum water depth were analyzed and are presented in Figure 6 (right column).
We found that: (1) Inundated area was mostly located along the river network and in the low-lying urban areas and cultivated land. Inundated extent and maximum water depth were obviously increased due to the bigger rainfall height. (2) In three FPs, maximum water depth in low-lying areas of urban land exceed 1 meter. In FP(B), cultivated land was also inundated with the maximum water depth of more than 1 meter. The results suggest that spatial distribution of flood information could be adopted in flood risk assessment and also help with the planning of flooding mitigation and adaptation strategies.
These findings imply that the XAJ-HiPIMS (1) presents good agreement with the observations; (2) can provide 2D presentation of flood processes and more information of flood characteristics at the basin scale; (3) provides a new perspective and methodology for ecological and environmental protection, flood prevention and management, and urban water landscape design.

Conclusions
This study presented a coupled hydrologic-hydraulic model, XAJ-HiPIMS. The model was validated in the Misai basin, Zhejiang Province, China by comparing its performance with measurements and statistically assessing the absolute relative error of flood-peak discharge, the difference of peak arrival time, and the Nash-Sutcliffe efficiency coefficient.
The present study: (1) demonstrates that the XAJ-HiPIMS can provide a good performance in basin-scale flood simulation and deliver a 2D viewpoint of flood characteristics; (2) reveals the potential application of XAJ-HiPIMS in ecological and environmental protection, urban flood management, landscape design, and climate change analysis; and (3) brings forward the necessity and recommendations for further research.
We found that: 1) Inundated area was mostly located along the river network and in the lowlying urban areas and cultivated land. Inundated extent and maximum water depth were obviously increased due to the bigger rainfall height. 2) In three FPs, maximum water depth in low-lying areas of urban land exceed 1 meter. In FP(B), cultivated land was also inundated with the maximum water depth of more than 1 meter. The results suggest that spatial distribution of flood information could be adopted in flood risk assessment and also help with the planning of flooding mitigation and adaptation strategies. These findings imply that the XAJ-HiPIMS 1) presents good agreement with the observations; 2) can provide 2D presentation of flood processes and more information of flood characteristics at the basin scale; 3) provides a new perspective and methodology for ecological and environmental protection, flood prevention and management, and urban water landscape design.