Assessment of the Relationship between Land Use and Flood Risk Based on a Coupled Hydrological–Hydraulic Model: A Case Study of Zhaojue River Basin in Southwestern China

Abstract

As an ecological consequence of intensified anthropogenic activities, more frequent extreme rainfalls have resulted in significant increases in water levels and discharge in southwestern China. This phenomenon presents a significant challenge in flood risk and ecological management. Land use is one of the major factors significantly affecting the flooding process, and it is inextricably tied to the ecological risk of floods. Hence, flood risk estimates based on land use are essential for flood control and land use planning. In this study, a coupled hydrologic–hydraulic model was developed to analyze the relationship between flood ecological risk and land use in order to provide new insights into current flood risk management practices. Ten real flood events (of different magnitudes) in the Zhaojue river basin (650 km²) were chosen to evaluate the credibility and performance of the coupled model’s application. Promising results were obtained, with sufficient reliability for flood risk assessment purposes. The results of our flood risk analysis also indicated that the model effectively reproduced overland flow and competently accounted for flood evolution. This work is significant in the understanding of the mechanism of the flood process and its relationship with land use, and it can be used in decision support for the prevention and mitigation of flood disasters and for land use planning.

1. Introduction

Floods are considered the most common and damaging natural hazards around the world [1,2,3,4], in which they cause over 50% of all casualties and more than 30% of economic losses [5]. The frequency of extreme rainfall events has grown with increases in population and the acceleration of urbanization, and the damage caused by floods has gradually become one of the biggest risks to social development worldwide [6,7,8]. At present, the frequency, intensity, loss, and ranges of impacts of flood disasters in the world have shown significant upward trends [9,10], demonstrating that flood control and disaster mitigation still face challenges. China has some of the most severe flood disasters in the world [11,12,13,14]. Floods have become a comprehensive problem that cannot be disregarded in sustainable development, which has increased pressure on local water resource management departments in flood management and prevention.

Flood risk reduction requires flood risk assessment [15,16,17], which is an essential component of flood control and disaster mitigation and plays an important role in flood management as a non-engineering measure [18]. With the rapid development of 3S technology, numerical calculation techniques, and flood simulation technology, flood risk assessment has received more and more attention in flood control and emergency work [18]. Studies have demonstrated that reasonable flood risk assessment can provide crucial information (such as inundation extent, water depth, flow velocity, peak flow, and arrival time) for strategy development and planning adaptation [19]. Attems et al. emphasized that flood risk assessment should be applied to flood management as a non-structural measure alongside structural measures [20]. Costache et al. suggested that flood risk assessment is an essential and mandatory procedure that should be included in flood management strategies [21]. Wu et al. proposed a hydrodynamic model for flood risk assessment in different environments that could be used to provide support for urban flood control and flood risk management [22]. Lowe et al. combined multiple models to systematically test various flood risks [23]. Mehr et al. used a forecasting system for flood risk assessment and achieved good results [24]. Furthermore, numerous large-scale flood risk map projects have been established by national institutions, thus attaching serious importance to flood risk assessment [25].

Previous studies have shown that land use associated with anthropogenic activities is the main factor affecting runoff and flooding [26,27,28,29,30], causing more flood risks [31,32,33]. Related research has suggested that land use change is the strongest contributor to changes in the flood process that may be directly responsible for more than 70% of floods worldwide [34,35]. Various land use types, such as urbanization [36], cultivated land [37], and vegetation changes [38,39], have disparate hydrological consequences. Land use structures and patterns also play important roles in the flood process [40,41]. Overall, land use planning has a potentially significant effect on flood risk. Ran et al. pointed out that flood risk management should be integrated and investigated in tandem with land use planning [42]. According to related research [43], reasonable land use planning based on flood risk analysis could decrease future flood risks by 39–50%, so increased flood risks caused by unreasonable land use also require attention [44]. Hence, flood risk assessment guided by the land–flood dual ecosystem theory needs to be studied [45].

Flood risk assessment in southwestern China has become more urgent due to its complex and unique characteristics [46,47,48,49]. For instance, 10% of China’s area comprises hills and mountains, with slopes of greater than 15° accounting for 60.0% of the terrain, demonstrating that this region is at serious flood risk [50]. Hence, compared to other regions, the flood control and disaster mitigation of China may face greater challenges.

As the core contents of flood risk assessment, coupled hydrological–hydraulic models have been increasingly appreciated in recent decades [51,52,53,54]. Coupled models can consider the effects of land use on the flooding process, so these models have been widely applied as important tools for flood management and land use planning [55,56,57,58,59,60]. On the one hand, hydrological models have major advantages over current alternatives in terms of modeling complex situations, computational efficiency, and imputation accuracy. On the other hand, hydraulic models have the ability to characterize hydrological behavior information in detail, laying the foundation for refined flood risk analysis [61]. Hence, hydrological–hydraulic coupled models can be used when a single model cannot meet simulation requirements. By linking hydrology to hydrodynamics, it is possible to increase the advantages and compensate for the weaknesses of both [57]. Flood behavior can be described more accurately from the different perspectives with their respective advantages. Lai et al. presented a coupled numerical model and successfully reproduced the flow regime of the Yangtze river basin [62]. Paiva et al. coupled an MGB–IPH hydrological model and a hydrodynamic model to simulate the flood process in the Amazon basin, and their results showed that the coupled model could reproduce the main hydrological characteristics [63]. Yin et al. applied coupled models to predict the flow and inundation of rivers and obtained information on the spatial distribution of floods [64]. Lian et al. employed a hydrology–hydrodynamic coupling model to simulate the flow of the Illinois river basin at different time scales, and the results showed that the accuracy of the coupled model was higher than that of a single model [65]. The same conclusion was drawn in a related paper [66].

Though there have been numerous studies based on coupled models for flood events, combining the advantages of both on the basin scale still faces many problems [67,68] regarding, e.g., the space–time scale [54,55,69], numerical instabilities [57,70], realistic flood events [71], validation issues [68,72], regional parameterization [54,73], infiltration issues [74], stream flow allocation issues [75], and connectivity issues [76,77]. To date, fully coupled modeling studies in realistic basins with both rainfall and infiltration processes have been rare, and coupling lumped rainfall–runoff (RR) models and fully dynamic hydraulic models have not been frequently reported in the literature [77].

Therefore, to fill the current research and practical gaps, our work aimed to develop a flood risk assessment approach using an integrated lumped RR and two-dimensional hydraulic model (2D) in realistic flood events incorporating both rainfall and infiltration processes. Land use information was adopted to solve the regional parameterization of the basin. The Zhaojue river basin in southwestern China was selected as a case study to assess the model’s performance. We also investigated the relationship between land use and flood risk in the basin. The proposed analytical method and knowledge gained from this study could help land use decision-makers and planners to identify flooded risk areas prior to future land use development. The remainder of this paper is organized as follows: Section 2 describes the methods and study area, Section 3 summarizes the research results, discussions are provided in Section 4, and conclusions are drawn in Section 5.

2. Materials and Methods

2.1. Study Area and Data Description

The study was conducted in a typical southwestern mountainous basin: the Zhaojue river basin with different hillslopes, soil properties, and vegetative covers. The Zhaojue river basin is located in Liangshan Yi Autonomous Prefecture, in the southwest of Sichuan province, China, with an area of 650 km² (28.00°–28.51° N, 102.69°–102.91° E) (Figure 1). It is a second-order tributary of the Jinsha River. There are three rainfall stations and a hydrologic station in this basin: ShenGuoZhuang station, SuoLuo station, Bier station, and Zhaojue hydrologic station. Elevation in the basin ranges from 2047 to 3815 m, and its surface relief and slope are significantly steep. These topographies are more likely to trigger extreme rainfall and flooding events, resulting in considerable infrastructure damage and losses of life.

Extreme rainfall events frequently occur in this region due to climate change and the influence of the southwest monsoon [78,79]. Additionally, due to its complex terrain and geological, ecological, and social backgrounds, rainfall events can easily cause flooding and threaten the safety of people’s lives and property [80,81] (Figure 2). Furthermore, related research suggests that there will be considerably more flooding in this area in the future [82]. Hence, the flood risk analysis of this area is of practical importance.

The dataset used in this study included rainfall, evaporation, discharge, DEM, and land use information. Details of the various data types used in this study are given in

Table 1. Details of various types of data used for developing the coupled model.

Data	Description	Source	Purpose
Rainfall (mm) Discharge (m3/s)	Hourly: (Flood data (2005–2010)) Monthly: ERA5 (2005–2010)	Hourly: hydrological administrations Zhaojue station (28.00° N,102.85° E) Monthly: Climate Forecasting Application Network from the European Centre for Medium-Range Weather Forecast (ECMWF)	Hourly: Model calibration/ validation Monthly: Hydrology analysis (Section 3.1)
Evaporation (mm)	Hourly: (Flood data (2005–2010))	Hydrological administration Zhaojue station (28.00° N, 102.85° E)	Model calibration/ validation
Temperature (°C)	Daily: 2005–2010 0.25˚ × 0.25˚ ERA5 precipitation and temperature data with a temporal resolution of 1 h	Climate Forecasting Application Network from the European Centre for Medium-Range Weather Forecast (ECMWF)	Hydrology analysis (Section 3.1)
Regional map (River network stations)	Vector format GCS_WGS_1984 1:50,000	Bureau of hydrology	Correction of delineated catchment boundaries and river network
Digital Elevation Model (FABDEM)	30 m × 30 m Publication date: 17 Dec 2021 GCS_WGS_1984 31.52°–32.72° N, 113.25°–114.77° E	University of Bristol https://data.bris.ac.uk/data/dataset/ (accessed on 8 April 2022 )	FABDEM is new a global elevation map that removes building and tree height Delineation of catchment boundaries Extraction of river network and floodplain cross-sections Input data for the hydraulic module of the coupled model
Landsat 4-5 TM	30 m × 30 m 2005-2010 (June) GCS_WGS_1984	United States Geological Survey (USGS) (https://eartheplorer.usgs.gov accessed on 10 April 2022)	Land use classification Assigning the surface flow to basin grids Determining the Manning coefficient as the input data for the hydraulic module of the coupled model
	30 m × 30 m
Google Earth	Raster format 2006	Google LLC	The base map of flood hazard map

.

2.2. Description of the Coupled Model

The coupled model was constructed based on a lumped hydrological model (XAJ) and a hydraulic model (2D model). In the coupling framework, the XAJ model was adopted to generate runoff, and the routing scheme was replaced by the 2D model (discretized surface runoff as input) to better simulate flash floods at the basin scale. The whole coupling process is illustrated in Figure 3.

With hourly rainfall, evaporation, and discharge data, the XAJ model served as a rainfall–runoff generator that generates surface runoff, interflow and groundwater flow after evapotranspiration, runoff generation, and runoff source partition (including a priori calibration) data (Figure 3, left). After these processes, the routing processes of interflow and groundwater flow still used the XAJ model (linear reservoir method). However, the routing technique of surface runoff was replaced by the 2D model to analyze flood evolution. The discretized surface flow was moved from upstream to downstream over the entire basin through a topography-based cell-to-cell connectivity sequence. Additionally, the Manning coefficient (derived from land use) was used to estimate the flow resistance to adjust the effect of underlying surface on water flow velocity (Figure 3, upper right). The model was verified by comparing the sum of the discharge at the basin outlet from the 2D model and the underground discharge from the XAJ model with the measured discharge. Additionally, more flash flood information (flood depth, velocity, inundation area, and arrival time) could be obtained by the coupled model. The analysis of this information could contribute to reducing flood hazards and associated consequences (Figure 3, lower right). The remainder of this section details the components of the coupled model.

2.2.1. Hydrologic Component (XAJ)

The XAJ model was developed by Zhao et al. [83]; it is one of the most famous lumped hydrological models based on saturation excess runoff and is widely used in the humid and semi-humid basins of China [84]. The XAJ model uses four main components: evapotranspiration, runoff generation, runoff source partition, and runoff routing (Figure 4). Here, the three-layer evapotranspiration sub-model was utilized to calculate the evaporation, and a tension water storage capacity curve was adopted to generate the total runoff.

Runoff generation is one of the most important modules of the model, and a parabolic curve was used to calculate the runoff yield, which could be approximately represented by the following equation:

$$\frac{\mathrm{f}}{F}=1 - {\left(1 - \frac{W_m^{\prime }}{W_{mm}^{\prime }}\right)}^{\mathrm{B}},$$

(1)

$$\{\begin{array}{cc}R=P-E-WM+W_0+WM{\left[1-\left(P-E+A\right)/W_{mm}^{\prime }\right]}^{1+B}& PE+A<W_{mm}^{\prime }\\ R=P-E-WM+W_0& PE+A\ge W_{mm}^{\prime }\end{array},$$

(2)

where $W_m^{\prime }$ is the tension water storage capacity of a point in the basin (mm), f is the fraction of the basin area for which the storage capacity is less than $W_m^{\prime }$, F is the whole basin area, $W_{mm}^{\prime }$ is the maximum value of $W_m^{\prime }$, $\mathrm{B}$ is the shape parameter of the storage capacity distribution, $R$ is the total runoff, $P$ is the rainfall, E is the evapotranspiration, WM is the average tension water storage capacity, $W_0$ is the initial soil water storage, and $A$ is the initial state of the basin.

The total runoff is composed of surface runoff, interflow, and groundwater runoff by a free water capacity curve:

$$\frac{\mathrm{FS}}{FR}=1 - {\left(1 - \frac{SM{F}^{\prime }}{\mathrm{SMMF}}\right)}^{\mathrm{EX}},$$

(3)

where $SM{F}^{\prime }$ is the free water storage capacity of a point in the basin, $\mathrm{SMMF}$ is the maximum value of $SM{F}^{\prime }$, $\mathrm{FS}$ is the fraction of the basin area for which the free water storage is less than $SM{F}^{\prime }$, FR is the runoff-producing area proportion, and $\mathrm{EX}$ is the exponent of the free water storage capacity curve.

The surface runoff was set to be routed with the unit hydrograph or the lag-and-route method. The interflow and groundwater flow were set to be routed with the linear reservoir method [85]. The outflow of each sub-basin was set to be routed to the basin outlet with the Muskingum method [86].

2.2.2. Hydraulic Component (2D Model)

The two-dimensional hydraulic model was used to solve the 2D shallow water equations (SWEs) based on the principles of the conservation of mass and the conservation of momentum under steady-state flow conditions [61,87,88]. The method was developed for simulating two-dimensional overland flow with spatially variable infiltration and microtopography using hydrodynamic flow equations [89].

In order to ensure well-balanced properties when the grids near the dry and wet interface in the calculation domain are processed separately, water depth is treated as the water level in this model [88]. The two-dimensional hydraulic model can be represented as follows:

$$\frac{\partial q}{\partial t}+\frac{\partial f}{\partial x}+\frac{\partial g}{\partial y}=s,$$

(4)

where $q$ is the vector containing the conserved flow variables, $x$ and $y$ are the Cartesian directions, $f$ and $g$ are the flux vector terms in the two Cartesian directions, $t$ is time, and $s$ is a source term.

Furthermore,

$$q=\left[\begin{array}{c}\eta \\ hu\\ hv\end{array}\right], f=\left[\begin{array}{c}hu\\ u^{2}h+\frac{1}{2}\\ uvh\end{array}g\left({\eta }^2-2\eta z_b\right)\right], g=\left[\begin{array}{c}hv\\ hvu\\ h{v}^2+\frac{1}{2}g\left({\eta }^2-2\eta z_b\right)\end{array}\right],$$

(5)

$$s=\left[\begin{array}{c}\begin{array}{c}0\\ -\frac{{\tau }_{bx}}{\rho }-gh\frac{\partial z_b}{\partial x}\end{array}\\ -\frac{{\tau }_{by}}{\rho }-gh\frac{\partial z_b}{\partial y}\end{array}\right],$$

(6)

where $\eta$ is the water surface elevation, $h$ is the total water depth, $u$ and $v$ are the two Cartesian velocity components,$g$ is the acceleration due to gravity,$z_b$ is the bed elevation above datum,$\rho$ is the water density, and ${\tau }_{bx}$ and ${\tau }_{by}$ are the friction stresses that can be represented as follows:

$${\tau }_{bx}=\rho C_{f}u\sqrt{u^2+v^2},{\tau }_{by}=\rho C_{f}v\sqrt{u^2+v^2},$$

(7)

$$C_f=\frac{g{n}^2}{h^{\frac{1}{3}}},$$

(8)

where $C_f$ is the roughness coefficient and $n$ is the Manning coefficient.

2.2.3. Coupled Component (Coupled Model)

The 2D hydraulic model described in Section 2.2.2 only deals with movement in the horizontal direction, even though water flows vertically in addition to horizontally between soils in the entire process of the water cycle. The model additionally lacks the ability to convert rainfall into surface runoff (it cannot handle the infiltration process). Hence, we had to integrate the surface runoff generated from Equation (3) into Equation (6) to build the coupled model. Hence, in the coupled model, the s variable in Equation (6) needed to be adjusted to the following:

$$s=\left[\begin{array}{c}\begin{array}{c}{R}_{sur}\\ -\frac{{\tau }_{bx}}{\rho }-gh\frac{\partial z_b}{\partial x}\end{array}\\ -\frac{{\tau }_{by}}{\rho }-gh\frac{\partial z_b}{\partial y}\end{array}\right],$$

(9)

where $R_{sur}$ is the surface runoff. Ultimately, by integrating Equation (9) into Equation (4), the coupled model could be constructed. In the coupled model, the infiltration component is represented by the XAJ model, and the routing scheme of the XAJ model is replaced by the two-dimensional (2D) hydraulic model.

The $R_{sur}$ obtained by the XAJ model is only the mean value of the basin, so a certain theory was needed to assign $R_{sur}$ to each grid in the basin to obtain the discretized surface runoff. The spatially discretized surface runoff based on land use can be represented as follows:

$$R_{sur, }{}_i=f\left(p-e,s,l,r,SM,EX\right),$$

(10)

where $R_{sur,}{}_i$ is the surface runoff of the i-th grid;$p$ is the rainfall;$e$ is the evaporation;$s$ is the soil water content of free water;$l$ is the land use;$r$ is the contributing area of runoff;$SM$ is the free water capacity;$EX$ is the parameter of the free water capacity curve;$p$ and $e$ are observed values;$s$,$r, SM, \mathrm{and} EX$ can be obtained by the XAJ model; and $l$ can be obtained from the land use information. Overall, by integrating Equations (9) and (10) into Equation (4), the coupled model could be constructed.

2.3. Data Preprocessing and Model Calibration

According to the requirements of Equations (1) to (10), the main data preprocessing steps were established with the aid of GIS tools and remote sensing data. We used a digital elevation model (DEM) as the main input for topographic parameter extraction in both hydrological and hydraulic modeling for the preprocessing of data regarding, e.g., sub-basins, drainage networks, slope, directions of water flow, and cross-sectional profiles. FABDEM was adopted to identify and extract these features. The average areal rainfall in each sub-basin was calculated with the Thiessen polygonal method [90]. Manning’s coefficient was determined according to land use to reflect the surface roughness in this study (

Table 2. Manning’s coefficient.

Land Use	Manning’s Coefficient (s/m1/3)
Urban land	0.016
Water	0.035
Grassland	0.03
Cultivated land	0.035
Forest land	0.075

), and such values were supported by a number of previous works [61,91].

The model calibration aimed to align the model’s results with observed data by adjusting the model parameters, which was crucial for model application. The PEST (Parameter Estimation) automatic calibration algorithm is a flexible and generic calibration tool that can be used in any of these models and has been widely applied [92]. Hence, the PEST algorithm was linked to the coupled model in the MATLAB environment for model calibration and validation in our study. The calibration procedure was used to optimize 13 parameters (

Table 3. The parameters of the XAJ model.

	Parameters	Physical Meaning	Range and Units [93,94,95]	Optimized Value
Evapotranspiration	K	Ratio of potential evapotranspiration to pan evaporation	0.5–1.1 (−)	0.91
	C	Coefficient of the deep layer	0.1–0.3 (−)	0.2
	WUM	Averaged soil moisture storage capacity of the upper layer	5–100 (mm)	5
	WLM	Averaged soil moisture storage capacity of the lower layer	50–300 (mm)	86
	WDM	Averaged soil moisture storage capacity of the deeper layer	5–100 (mm)	35
Runoff generation	B	Exponent of the distribution to tension water capacity	0.1–2 (−)	0.34
	IMP	Percentage of impervious and saturated areas in the basin	0.01–0.1 (%)	0.01
Runoff sources partition	SM	Areal mean free water capacity of the surface soil layer	5–100 (mm)	85
	EX	Exponent of the free water capacity curve influencing the development of the saturated area	1–1.5 (−)	1.5
Runoff routing	KSS	Outflow coefficients of the free water storage to interflow relationships	0.01–0.7	0.23
	KG	Outflow coefficients of the free water storage to groundwater relationships	0.01–0.7 (−)	0.47
	KKI	Recession constants of the interflow storage	0.05–0.95 (−)	0.74
	KKG	Recession constants of the groundwater storage	0.9–0.999 (−)	0.998

) related to evapotranspiration (K, C, WUM, WLM, and WDM), runoff generation (B and IMP), runoff source partition (SM and EX), and runoff routing (KI, KG, KKI, and KKG) following the work of Zhao et al. [83]. Three main factors were considered during parameter optimization: the upper and lower bounds of each parameter, the objective function, and the iteration termination criteria. The upper and lower bounds were determined by the model’s recommended parameter range and their physical meaning to avoid the selection of localized parameter values (Table 3) [83,93,94]. The objective function was calculated as follows:

$$\mathrm{F}={{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Q}}_{\mathrm{s}}-{\mathrm{Q}}_{\mathrm{o}}\right)}^2,$$

(11)

where $\mathrm{F}$ is an objective function, ${\mathrm{Q}}_{\mathrm{s}}$ is the simulated discharge, and ${\mathrm{Q}}_{\mathrm{o}}$ is the observed discharge.

Ten significant flood events in the period of 2005–2010 were available for the study. We used six flood events to calibrate the model, and the others were used for validation purposes (the results are introduced in Section 3.2).

A uniform grid spacing of 30 m was used, so the entire basin had more than 700,000 grid points. At each grid point, hydraulic conductivity based on the remote sensing data was defined. According to the definition of a drainage watershed, all surface water located within the watershed ultimately converged to an outlet. Therefore, the coupled model only needed two types of “known” boundary conditions, i.e., non-transmissible (wall) for all boundaries except for an outlet and transmissible (open) boundary conditions for outlets. A time interval of 1 s was used for the simulation. The simulations were carried out in a continuous fashion and included the representation of subsurface and surface dynamics, and the output time interval was 1 h. The time step used for hydraulic routing was 0.5 s, which could be adaptively increased according to the local CFL conditions.

2.4. Statistical Method

The performance of the coupled model could be computed using statistical indices and graphical comparisons. The Nash–Sutcliffe efficiency (NSE), relative error (RE), peak discharge error (PE), and time difference of arrival peak discharge ($\Delta \mathrm{T}$) were employed to estimate the model’s performance [96].

$$\mathrm{NSE}=1 - \frac{{{\displaystyle \sum }}^{ }{{(\mathrm{Q}}_{\mathrm{o}}{-\mathrm{Q}}_{\mathrm{s}})}^2}{{{\displaystyle \sum }}^{ }{{(\mathrm{Q}}_{\mathrm{o}}-\overline{{\mathrm{Q}}_{\mathrm{o}}})}^2},$$

(12)

$$\mathrm{RE}=\frac{{{\displaystyle \sum }}^{ }{\mathrm{Q}}_{\mathrm{s}}-{{\displaystyle \sum }}^{ }{\mathrm{Q}}_{\mathrm{o}}}{{{\displaystyle \sum }}^{ }{\mathrm{Q}}_{\mathrm{o}}} \times 100\%,$$

(13)

$$\mathrm{PE}=\frac{\mathrm{MAX}\left({\mathrm{Q}}_{\mathrm{s}}\right){-\mathrm{MAX}(\mathrm{Q}}_{\mathrm{o}})}{{\mathrm{MAX}(\mathrm{Q}}_{\mathrm{o}})} \times 100\%,$$

(14)

$$\Delta \mathrm{T}=\left|T_{\mathrm{MAX}\left({\mathrm{Q}}_{\mathrm{s}}\right)}-T_{{\mathrm{MAX}(\mathrm{Q}}_{\mathrm{o}})}\right|,$$

(15)

where ${\mathrm{Q}}_{\mathrm{o}}$ is the observed discharge, ${\mathrm{Q}}_{\mathrm{s}}$ is the simulated discharge, $\overline{{\mathrm{Q}}_{\mathrm{o}}}$ is the mean of the observed discharge,$\mathrm{MAX}\left({\mathrm{Q}}_{\mathrm{s}}\right)$ is the simulated peak discharge, $\mathrm{MAX}\left({\mathrm{Q}}_{\mathrm{o}}\right)$ is the observed peak discharge, $T_{\mathrm{MAX}\left({\mathrm{Q}}_{\mathrm{s}}\right)}$ is the arrival time of the observed peak discharge, and $T_{{\mathrm{MAX}(\mathrm{Q}}_{\mathrm{o}})}$ is the arrival time of the simulated peak discharge.

3. Results

3.1. Hydrological Characteristics and Land Use Change in the Study Area

The consistency and stability of the basin’s hydrology impacted the simulation results of the coupled model. Climate and land use change are two important factors in a coupled model’s consistency and stability [46,97], so studies of land use and climate in the study area are significant for subsequent analyses of the flood simulation.

Climate change was found to have a significant impact on the relationship between rainfall and runoff in the basin. The temperature distribution in the Zhaojue river basin during the study period is shown in Figure 5a. We found that the temperature of the basin significantly varied from −10.0 to 20 °C, with a hot summer and cold winter. The annual average temperature was around 10 °C for many years. The highest temperature occurred around July, generally 19.0 °C, and the lowest temperature occurred around January, generally −5.0 °C. The temperature change correlated with the hydrological characteristics of the study area. The annual hydrological characteristics of the basin during the study period are shown in Figure 5b. The mean annual rainfall was about 1400 mm, the mean runoff depth was 600 mm, and the runoff coefficient was relatively stable at about 0.45. These results further suggest that the study area had significant hydrological consistency during the study period.

Land use change has an important impact on the flood process at the basin scale. The area percentages of land use in the basin during the study period are shown in

Table 4. Area percentages of land use.

Land Use Type	Cultivated Land	Forest	Grassland	Water	Unused	Urban Land
Area percentage (%)	15.03	41.63	43.04	0.1	0.01	0.19

, and the maps of land use are shown in Figure 6. There are six types of land use in the study area, the main ones being cultivated land, grassland, and forest, which account for over 99% in total, and it is worth noting that the degree of land use change is not obvious. Consequently, land use changes will have little impact on the consistency of the hydrology, which is beneficial for our research.

3.2. Performance Assessment of Coupled Model for Real Flood Events

The performance of the coupled model was evaluated by analyzing the discharge process at the outlet gauging station, as indicated in Section 2.4. Considering the availability and consistency of hydrological data, six flood events were selected for model calibration and four flood events were used for model validation. The results of the coupled model are presented in

Table 5. Simulation results of the coupled model.

	Flood	NSE	RE (%)	PE (%)	ΔT (h)
Calibration	20050707	0.77	−25.76	−42.90	1
	20050916	0.95	−4.02	−16.83	2
	20060706	0.87	6.74	−18.82	1
	20060917	0.67	28.50	−0.15	1
	20070605	0.80	33.95	12.64	1
	20070916	0.76	−10.06	3.10	1
Validation	20080703	0.77	22.37	−8.41	1
	20090726	0.94	1.22	−15.32	1
	20090920	0.90	16.75	−12.00	1
	20100821	0.48	81.92	3.88	2

.

According to these statistical indicators, NSE values were generally within an acceptable range between 0.48 and 0.95, with an average value of 0.79 in all flood events. Event no. 20050916 presented the highest values, whereas event no. 20100821 had the lowest values due to the complex topography and rainfall distribution [46]. The absolute RE in the calibration period ranged between 4.02% and 33.95%, with an average value of 18.17%, and between 1.22% and 81.92% in the validation period. The PE indices were rather satisfactory, with low error rates (within ±20%) in most events (except for event no. 20050707 due to antecedent precipitation). The investigation of the arrival time of the flood peak revealed that the model accurately simulated the timing of the confluence. A coupled model can simulate the peak discharge and routing process in enough detail to comprehensively consider the impacts of land use on flooding [98], so the arrival time of the flood peak from the coupled model was well matched with the observations.

In addition to the statistical indicators, simulated hydrographs from the coupled model were plotted with the observed discharge hydrographs to compare the observed and simulated discharge processes (Figure 7). Figure 7 shows that the simulated stage hydrograph agreed well with that observed at the outlet of the basin. The agreement indicates that the developed equation and method could adequately describe the observed surface runoff process. However, the peak discharges were underestimated for most of the flood events, which may have been caused by the infiltration excess runoff due to topography [81]. It is also noteworthy that the discharges were overestimated for the low flow; this is due to agricultural storage (15% of all land, Table 4) and the depression storage or the antecedent precipitation error [99]. During storms, the coupled model could accurately simulate the peak flow and the arrival time in most flood events, demonstrating the strengths of the hydraulic model for confluence.

In summary, although some simulated events (event nos. 20050707 and 20100821) had certain errors, the overall accuracy of the results based on the evaluation criteria in Section 2.4 was within an acceptable range. As such, the coupled model can be applied in flood risk assessment.

3.3. Flood Risk Response of Different Land Uses

A flood risk management method is essential for protecting downstream residents. This method can be used to prepare flood risk maps in order to analyze the impact of land use on the flood process [100,101,102]. These maps could provide more decision support for land use planners, authorities, and governments to develop appropriate policy responses to mitigate or prevent potential damage due to flooding [33].

For this paper, we selected the flood event 20060706 as an example for flood risk analysis because it showed a typical flow distribution. We set the threshold of inundated depth to 0.01 m. First, we investigated how the inundation dynamics evolved over time, and water depth maps at various time points are presented in Figure 8. The overland runoff flowed from the hillside to the channel and eventually along the channel. With the increase in rainfall intensity, the scope of surface flow also expanded, and the coverage of water flow widened (from 5 to 18 h). When the rainfall approached the peak moment (9 h), the inundation extent obviously increased. These results show that the spatial distribution of flood behavior calculated with the coupled model could be used for flooding mitigation and land use planning.

Based on the flood evolution process, the changes in the flood inundation areas and average water depth were calculated and are presented in Figure 9. Both of them gradually increased over time (Figure 9); at 20 h (near maximum cumulative rainfall), flood inundation areas reached their peak, accounting for 1.3% and reaching 8.45 km². As can be seen in the inundation extent information shown in Figure 8, the roads close to the river (as well as cultivated and urban land) contributed the most to the inundation areas. The results indicate that these risk loci were prone to flooding during significantly heavy rainfall. After 20 h, the spatial distribution of flooded areas was basically stable, indicating that the duration was long and the damage area was relatively large.

The maximum water depth is an important indicator for flood evaluation and hazard estimation. In this study, the maximum water depth was analyzed, and the maps of it are presented in Figure 10a (their zoomed-in views are also illustrated in Figure 10b,c). The maximum water depth was less than 1 m in most areas, and the maximum value was 1.22 m (Figure 10a). Taken as a whole, extensive inundation could be observed in main roads, intersections, and dense residential quarters (Figure 10b,c). Heavy rainstorms within these areas will cause more severe, more intense, and longer floods. Even worse, Zhaojue county has a low elevation and is located at the outlet of the basin, so flooding will pose significant threats to people’s lives and property when rainstorms occur. Therefore, we recommend implementing river regulation projects and structural reinforcement in the upper reaches of Zhaojue river.

Flood risk was divided into four categories based on water depth: 0.05–0.1 m, 0.1–0.2 m, 0.2–0.4 m, and >0.4 m. The inundation areas and percentages under different levels were calculated and are shown in

Table 6. Flood inundation in different grades.

Depth (m)	Inundation Area (km2)	Area Ratio (%)
[0.05, 0.1)	5.17	1.00
[0.1, 0.2)	2.80	0.54
[0.2, 0.4)	1.27	0.25
[0.4, +∞)	0.47	0.09

. The results show that a total area of 9.17 km² (about 1.88% of the total area) was at risk of flooding. Although the overall proportion was not large, it was mostly in the landslide-prone areas, in which the harm caused was significant. In addition, residential areas and cultivated land accounted for a large proportion in the inundation extent (Figure 10a), followed by roads, thus allowing for serious disasters.

The modeled velocity over highly irregular topography was of significant interest as it could clearly reflect the impact of the topography on the flood flow. The flow velocity in the basin was found to vary between 0.01 and 1.0 m/s—mostly below 0.1 m/s—and sometimes between 0.1 and 0.3 m/s (Figure 11a). The velocities in the narrow reach of the channel were higher than those in the wider reach. Furthermore, no non-physical (unrealistically high) velocity was discerned, which reflected the satisfactory performance of the present model from the numerical perspective. “Derivative” metrics, such as shear stress and vorticity, could be derived from velocity information, which would be helpful for understanding the effects of land use and basin hydrology on soil erosion and sediment transport. A map of the arrival time of the maximum water depth (relative to initial time of the event) is presented in Figure 11b, which clearly shows the flow process from the hillside to the river.

Furthermore, some high-risk sites are shown in Figure 12, which shows that the risk points were mostly located in the downstream river channel (Figure 12b) and surrounding villages (Figure 12c). The arrival time of the maximum water depth is an important variable, providing a reference for disaster relief time and relief route—for example, transferring people before a certain inundation depth. Hence, it is recommended to prioritize the evacuation of large numbers of people ahead of the arrival of some extreme flooding using information gained from the arrival time map (Figure 11b). Overall, the diverse spatial patterns of the changes in flood occurrence reflect various aspects of rainfall intensity, topography, and land use over the basin, demonstrating that the results were reasonable.

4. Discussions

4.1. The Coupled Hydrological–Hydraulic Model and Flood Process

Although the authors of numerous studies have investigated the flood process by integrating hydrological and 2D models (e.g., [103,104,105]), it is relatively rare for lumped rainfall–runoff (RR) models to be coupled with hydraulic models at the basin scale, even though some scholars have called for more research into this topic [77]. It is still a significant challenge to allocate overland flow in lumped RR to the boundary of a hydraulic model (usually uniform distribution [61,91,106]). In our coupled approach, the application effect of the XAJ model in runoff production has widely been proven possible [83,94], and 2D models have also been successfully adopted to assess anthropogenic impacts on the flood process at the basin scale [107,108,109]. Therefore, we are confident in the coupling method of this paper. Furthermore, compared to previous studies, we have expanded the scope of application (not only in experimental conditions [110]) and also addressed the computational domain problem with our coupled model [76,77].

Accordingly, the proposed coupled model has demonstrated its value in our study area. The results (Table 5) indicate the model’s strong performance during the calibration and validation periods, which is in line with related research with similar methods [54,105]. In terms of the NSE value (average NSE = 0.79), the results showed an acceptable accuracy that was consistent with other studies (e.g., average NSE = 0.80 in [58]). Compared to a single hydrological model (e.g., average NSE = 0.73 in Longchangqiao station, 0.65 in Huangguoshu station, and 0.68 in Huangmaocun station in [47]), the coupled model reported in this paper showed considerable improvements. Regarding the arrival time (average $\Delta \mathrm{T}$ = 1.2), the performance of our model was similar to or better than that of other studies (e.g., [71,98]). Wang et al. adopted the same coupled model to simulate floods at different magnitudes, but the results presented in this paper are better in terms of small and medium flood events due to the more precise surface flow assignment in our work (the surface flow was homogenized in that study) [91]. Similarly, the results of this study are better than those of our previous research with the same model [61] due to the fact that the study area of the previous research contained many dams that the coupled model failed to consider. However, we still observed some errors in our study (20100821 in Table 5) due to errors in the meteorological data or runoff production caused by the complex topography [81,111].

4.2. The Impact of Land Use and Topography on Flood Process

Munar et al. reported that the coupled model is a promising approach to investigate the impact of land use and topography on the flood process [57]. Related scholars have also proposed the necessity of applying the coupled model in flood research with complex terrains [88,112]. Therefore, we explored the relationship between topography and the flood process. Our flood evolution map (Figure 8) showed that the flow movement was basically consistent with the terrain. Additionally, due to the lower terrain, downstream areas were found to be more vulnerable to flooding than upstream areas (Figure 10), which puts more pressure on the flood control of downstream cities. These findings are consistent with those of Possa [113].

Land use was found to have significant impacts on runoff production and flood routing through different forms. The velocity map (Figure 11a) showed that different land uses had different effects on the dynamics of floods: the flow velocity sequentially increased in the vegetation, cultivated, and building zones. This phenomenon indicated that the interception rate of vegetation to flooding was high and that the urban land reduced the convergence time. The maximum water depth map (Figure 10) showed the impact of different land uses on floods: the vegetation zones were hardly submerged (less than 0.01 m), indicating the weakening effect of vegetation on flooding. However, there were more flood risk points in the urban area (Figure 10 and Figure 12), indicating that urbanization will lead to the aggravation of rainstorms and flood intensity.

4.3. Strengths and Weaknesses

The coupled model proposed in this paper combines the advantages of both hydrological and hydraulic models and could be used to simulate an entire basin. The backwater problem (i.e., the effects of building and tree heights) mentioned in previous research [114,115] has been addressed as efficiently as possible by FABDEM [116]. Due to the unevenness of underlying surfaces (microtopography), overland flow depths and speeds are highly variable and discontinuous in space and time [89]. FABDEM allowed us to significantly improve the consistency of the model and will allow for better representation of elevation, avoiding the “nonflow phenomenon” that occurs in some grids at low runoff production rates. Furthermore, in terms of accuracy, FABDEM is superior to existing elevation products, and its properties are well suited to hydraulic simulations [116].

There are three types of coupling: full coupling, internal coupling, and external coupling [117]. To facilitate such model selection and further promote the coupling of hydrological and hydraulic models, the proposed coupled model was designed with external coupling, which can incorporate more models to replace the XAJ model in an expandable way. Additionally, external coupling is the simplest and most common approach that can be highly extensible and combined without modifying the code of each model. We believe that we were able to better understand the flood process by joining the locally best performing hydrological and hydraulic models.

Despite these advantages, the coupled model has several limitations.

Begnudelli et al. pointed out that roughness has a vital role in routing time [118], and previous research has suggested that the Manning coefficient is a more sensitive parameter [98]. Hence, the impact of land use on roughness in our work requires further study. Additionally, Manning’s coefficient is related to terrain and water depth [89], so the construction of a comprehensive solution for the Manning coefficient would be helpful to improve the model’s performance. Furthermore, although the agreement between the coupled model and the observed discharge at the outlet was satisfactory, there is still a need to evaluate the velocities and water depths. This will greatly benefit from the linking of more data sources, such as the Surface Water and Ocean Topography (SWOT) mission or more detailed data sources (e.g., radar or more gauge stations) [119].

In addition, the scale problem of flood simulation needs to be seriously considered [120]. Coarser resolution datasets may impact model applications [57]. Although the resolution of the dataset used in this paper was consistent with that of previous studies [19,91,121] and can be used for flood simulation, its impact on the model’s performance was not considered; we need to further study the model’s performance under different resolutions. Furthermore, the uncertainty of input data impacts model results, and rainfall is considered the most important source of uncertainty [122]. Accordingly, to decrease the impact caused by rainfall uncertainty, the average rainfall of three rainfall stations was used to represent the rainfall distribution in this study. Topography is another source of uncertainty; although FABDEM was adopted to characterize the topographic changes in the basin, which is helpful for flood simulation [116], its application effect needs further evaluation.

5. Conclusions

The impact of land use on the hydrological process is particularly pronounced. In this paper, we have presented a systematic approach based on the effective coupling of hydrological and hydraulic models to assess the effect of land use on flood risk. In the coupled model, the XAJ model was adopted to generate overland flow and the floodplain evolution was implemented by a 2D hydraulic model. A strong agreement between the simulated and measured discharge results was achieved. This integrated approach could be a useful strategy for flood and land use management and could complement existing flood planning tools in southwestern China. The key findings are summarized as follows.

(1)

According to the analysis of hydrological characteristics (Section 3.1), it was concluded that the hydraulic and hydrologic features did not significantly change during the study period. The relationship between rainfall and runoff had a significant consistency. This was beneficial for coupled model application.

(2)

The results (Section 3.2) showed that the coupled model had a strong applicability to the basin in southwestern China. The model was able to reproduce observed hydrographs, and discharge was well simulated (Table 5 and Figure 7). We found that 80% of the flood events had an NSE value exceeding 0.70. In particular, the time of concentration and peak flows were well simulated: the simulated peak time matched well with the observed values with an error of less than ±3 h, and the PE values of most flood events were within 20%.

(3)

A unique advantage of the coupled model is that in addition to the flow hydrograph, it offers inundated area, water depth, and velocity information that is fundamental for reliable flood risk analysis. Our flood risk analysis results (Section 3.3) suggested that the coupled model is suitable for simulating inundation and could provide an important tool for flood management to reduce damage in terms of lives and property in the Zhaojue river basin.

(4)

The flood risk maps (Section 3.3) indicated that topography and land use played the most major roles in flood wave attenuation and delay. The flow velocity map (Figure 12a) showed that different land use types had different impacts on flood process. The order of flow velocity was as follows: forest < grassland < cultivated land < urban land. This order also reflected the degree of interception of floods by different land use types. These results showed that the strengthening of land management has played an important role in reducing flood risk.

Although our preliminary test demonstrates the relationship between the flood process and land use, there are still many questions to be addressed. Ongoing research can focus on the spatial verification of flood inundation extents with more measured data or radar data. Additionally, the Manning coefficient adopted in this paper was calculated from relevant research results; thus, the calculation of more accurate Manning coefficients to improve the coupled model’s performance is worth exploring. Furthermore, the effect of the resolution of remote sensing data on the coupled model’s results is also worth investigating.