Ratio limits of water storage and outflow in a rainfall–runoff process

Flash floods typically occur suddenly within hours of heavy rainfall. Accurate forecasting of flash floods in advance using the two-dimensional (2D) shallow water equations (SWEs) remains a challenge, due to the governing SWEs being difficult-to-solve partial differential equations (PDEs). Aiming at shortening the computational time and gaining more time for issuing early warnings of flash floods, constructing a new relationship between water storage and outflow in the rainfall–runoff process is attempted by assuming the catchment as a water storage system. Through numerical simulations of the diffusion wave (DW) approximation of SWEs, the water storage and discharge are found to be limited to envelope lines, and the discharge/water-depth process lines during water rising and falling showed a grid-shaped distribution. Furthermore, if a catchment is regarded as a semi-open water storage system, then there is a nonlinear relationship between the inside average water depth and the outlet water depth, namely, the water storage ratio curve, which resembles the shape of a plume. In the case of an open channel without considering spatial variability, the water storage ratio curve is limited to three values (i.e., the upper, the steady, and the lower limits), which are found to be independent of meteorological (rainfall intensity), vegetation (Manning's coefficient), and terrain (slope gradient) conditions. Meteorological, vegetation, and terrain conditions only affect the size of the plume without changing its shape. Rainfall, especially weak rain (i.e., when rainfall intensity is less than 5.0 mm h⁻¹), significantly affects the fluctuations of the water storage ratio, which can be divided into three modes: Mode I (inverse S-shape type) during the rainfall beginning stage, Mode II (wave type) during the rainfall duration stage, and Mode III (checkmark type) during rainfall end stage. Results indicate that the determination of the nonlinear relationship of the water storage ratio curve under different geographical scenarios will provide new ideas for simulation and early warning of flash floods.

1 Introduction

Flood disaster is a significant global health and economic threat. Disastrous floods have caused millions of fatalities in the 20th century and billions of dollars in direct economic losses each year (Merkuryeva et al., 2015; Merz et al., 2021; Ruidas et al., 2022). According to statistics (Lee et al., 2020), from 2001 to 2018, over 2900 floods caused over 93 000 deaths and over USD 490 billion in economic damages worldwide. Based on 250 m resolution daily satellite images of 913 major flood events during the same period, the total area inundated by floods is estimated to be 2.23×10⁶ km², and the directly affected population is estimated to be 255 to 290 million (Tellman et al., 2021). With the influence of climate change and extreme El Niño events (Ward et al., 2014; Cai et al., 2014), flood events caused by extreme precipitation are occurring frequently in many regions around the world (Kirezci et al., 2020; Najibi and Devineni, 2018; Almazroui, 2020). From 2020 to 2023, catastrophic floods caused by several extreme rainfall events were reported in Germany (Tradowsky et al., 2023), China (Hsu et al., 2021), Italy (Valente et al., 2023), Japan (Kobayashi et al., 2023), Pakistan (Nanditha et al., 2023), and other developed or developing countries and regions, even in some desert areas, e.g., in the Taklimakan Desert and the Atacama Desert, as reported by Li and Yao (2023) and by Cabré et al. (2023), respectively. Research shows that under a high emissions scenario, at latitudes above 40° N, compound flooding could become more than 2.5 times as frequent by 2100 compared to the present (Bevacqua et al., 2020). This means that in the future, the fraction of the global population at risk of floods will be increasing.

Flood simulation provides an effective means of flood forecasting to reduce the loss of property and life in flood-threatened areas around the world. In particular, weather-prediction-based distributed hydrological/hydraulic models are considered to be an effective strategy for flood simulation (Ming et al., 2020). Hence, a large number of scholars are committed to improving the simulation efficiency or simulation accuracy of distributed hydrological/hydraulic models. Accordingly, they have developed many forms of hydrological models and hydrodynamic models in the past decades. Among them, the hydrological models include Stanford Watershed Model IV (SWM) (Crawford and Linsley, 1966), the SHE/MIKESHE model (Abbott et al., 1986), the TANK model (Sugawara, 1995), Soil and Water Assessment Tool (SWAT) (Arnold and Williams, 1987), and TOPMODEL (Beven and Kirkby, 1979). The hydrodynamic models include the one-dimensional (1D) Saint-Venant equations (Köhne et al., 2011), the two-dimensional (2D) shallow water equations (SWEs) (Camassa et al., 1994), and the three-dimensional (3D) integrated equations of runoff and seepage (Mori et al., 2015). In addition, a variety of hydrological–hydrodynamic coupling models have also been proposed by Kim et al. (2012), Liu et al. (2019), Hoch et al. (2019), and other scholars. In particular, SWEs are the main governing equations for simulating floods. However, flood simulation based on SWEs is a time-consuming process due to its governing equations being a hyperbolic system of first-order nonlinear partial differential equations (PDEs) (Li and Fan, 2017). Therefore, many scholars attempted to improve the efficiency and accuracy of flood simulation through computer technology, e.g., applying GPU parallel computing (Crossley et al., 2010) or advanced numerical schemes (Sanders et al., 2010). For hydrological studies, the performance of hydrological modeling is usually challenged by model calibration and uncertainty analysis during modeling exercises (Wu et al., 2021).

Efficient and stable solutions for hydrodynamic models has long been an important issue in flood forecasting. Since the SWEs are nonlinear hyperbolic PDEs, the increase in the calculation domain and the increase in the degree of discreteness will greatly increase the difficulty of solving SWEs. In addition, when using high-resolution terrain to improve model calculation accuracy, non-physical phenomena such as false high-flow velocity in steep terrain will also occur, resulting in calculation distortion and a sharp increase in calculation time. Hence, we try to ignore the complex exchange/transfer process of mass and momentum (hydrodynamic models), and we also abandon the empirical relationships (hydrological models) between the input (precipitation), the transmission (flow rate), and the output (discharge) in the catchment area. A catchment is regarded as a semi-open water storage system, and the complex problem is simplified into three megascopic variables, i.e., inflow, water storage, and outflow. For one watershed, the complex internal flow processes could be ignored if the physical mechanism between inflow, water storage, and outflow can be found under different meteorological, geographical, and geological conditions. In other words, if we can give a physical-based relationship between the three megascopic variables, then flood forecasting will become much simpler.

2 Methods

An arbitrary catchment (Fig. 1b) could be assumed to be a conceptual water tank (Fig. 1a). In this water tank, according to the law of conservation of mass, the complex confluence process of surface runoff could be neglected, and it can be described only by the relationship between input, storage, and output, which can be expressed as Eq. (1),

$$\begin{array}{}\text{(1)}& \underset{\text{storage}}{\underbrace{A\cdot {\displaystyle \frac{\mathrm{d}H}{\mathrm{d}t}}}}=\underset{\mathrm{rainfall}}{\underbrace{R\cdot A}}-\underset{\mathrm{infiltration}}{\underbrace{I\cdot A}}+\underset{\mathrm{exfiltration}}{\underbrace{F\cdot A}}-\underset{\mathrm{evaporation}}{\underbrace{E\cdot A}}-\underset{\mathrm{discharge}}{\underbrace{{\displaystyle \frac{Q}{A}}\cdot A}},\end{array}$$

where A is catchment area (m²), t is time (s), H is internal average water depth (m), R is rainfall intensity (m s⁻¹), I is infiltration (m s⁻¹), F is exfiltration (m s⁻¹), E is evaporation (m s⁻¹), and Q is discharge (m³ s⁻¹).

Figure 1 Conceptual schematic of the distributed runoff model (DRM) and numerical model. (a) Conceptual water tank; (b) conceptual catchment; (c) impermeable conceptual slope model; (d) design rainfall.

In this section, attention is focused on the surface flow of runoff, so the runoff–atmosphere moisture exchange (evaporation) and runoff–soil moisture exchange (infiltration and/or exfiltration) are not considered. Zhu et al. (2020) validated the effectiveness of a diffusion wave (DW) approximation of shallow water equations by numerical simulations for simulating ground surface runoff,

$$\begin{array}{}\text{(2)}& {\displaystyle \frac{\partial h}{\partial t}}-\mathrm{\nabla }\left({\displaystyle \frac{h^{\mathrm{5}/\mathrm{3}}}{n_m\sqrt{\left|S\right|}}}\mathrm{\nabla }\left(h+z\right)\right)=R,\end{array}$$

where h is water depth (m), z is elevation (m), n_m is Manning's coefficient (s m${}^{-\mathrm{1}/\mathrm{3}}$), and S is the slope gradient.

To improve the computational efficiency of the hydrodynamic model, after strict mathematical derivation according to the basic hydrodynamic equation and the law of conservation of mass, Zhu et al. (2022) proposed a hydrological–hydrodynamic integrated model, i.e., distributed runoff model (DRM), as

$$\begin{array}{}\text{(3)}& \left\{\begin{array}{l}\frac{\mathrm{d}H}{\mathrm{d}t}=R-q,\\ H=\mathit{\eta }h=\mathit{\eta }{\left(\frac{n_m}{\sqrt{S}}\right)}^{\mathrm{0.6}}{q}^{\mathrm{0.6}}{\left(\frac{A}{B}\right)}^{\mathrm{0.6}},\end{array}\right.\end{array}$$

where $q=Q/A$ is conceptual outflow (m s⁻¹), η is the water storage ratio, and B is the outlet width (m).

3 Limits and plume shape of water storage ratio curve

The conceptual hydrological model takes the inside average water depth (H) in the catchment area as the independent variable (Eq. 1). However, the hydrodynamic equations take the water depth at any outlet (h) as an independent variable (Eq. 2). If a relationship between the inside average water depth (H) and outlet water depth (h) can be established, then this relationship will have both hydrodynamic and hydrological characteristics. Therefore, to find the H–h relationship, an impermeable conceptual slope model was built as shown in Fig. 1c, and numerical simulations were performed using the diffusion wave (DW) approximation (Eq. 2) of shallow water equations (SWEs). The water storage ratio is defined as the inside average water depth (H) divided by the outlet water depth (h). Firstly, the numerical simulations are performed under a designed rainfall condition; that is, rainfall intensity is 10 mm h⁻¹, and rainfall duration is 90 min with a total time of 180 min as shown in Fig. 1d. From the time-dependent water storage ratio ($H/h$) under different catchment area (Fig. 2a), it can be seen that continuous rainfall will cause the water storage ratio ($H/h$) to gradually decrease from the initial value of 1.0 (upper limit) to a stable value, which is approximately 0.625 (steady limit). When the rainfall ends, the value of the water storage ratio ($H/h$) decreases first and then increases, showing a U-shaped curve with a lower limit, which is approximately 0.4125. Afterward, the water storage ratio curves are under 10 kinds of catchment area (Fig. 2b), three kinds of Manning's coefficient (Fig. 2c), four kinds of slope gradient (Fig. 2d), and four kinds of rainfall intensity (Fig. 2e); conditions are obtained from parametric analyses and collected in Fig. 2f.

Figure 2 Water storage ratio curves. (a) Time-dependent water storage ratio under different catchment areas with 10 mm h⁻¹; (b) water storage ratio curves under 10 kinds of catchment area; (c) water storage ratio curves under three kinds of Manning's coefficient; (d) water storage ratio curves under four kinds of slope gradient; (e) water storage ratio curves under four kinds of rainfall intensity; (f) collection of the above 21 water storage ratio curves. Three limit lines envelop all water storage ratio curves, i.e., upper limit ($H/h=\mathrm{1.0}$), steady limit ($H/h=\mathrm{0.625}$), and lower limit ($H/h=\mathrm{0.4125}$).

Finally, it is found that water storage ratio curves resemble the shape of a plume. When the water outlet depth is the same, the water storage ratio ($H/h$) of the water-rising limb is higher than that of the water-falling limb. Furthermore, in the case of an open channel without considering spatial variability, there are three limits (the upper, the steady, and the lower limits) of the water storage ratio curves, which are found to be independent of meteorological (rainfall intensity), vegetation (Manning's coefficient), and terrain (slope gradient) conditions. Meteorological, vegetation, and terrain conditions only affect the size of the plume without changing its shape which is anchored by three limits. This means that the three limits and the water storage ratio curves provide a key to establishing a relationship between the hydrodynamic models and the hydrological models.

Figure 3 Discharge/water-depth process lines during water rising and falling. (a) Discharge process lines during water rising and falling; (b) gradient lines of discharge process line during water rising and falling; (c) schematic diagram of looped rating curve; (d) water depth process lines during water rising and falling; (e) gradient lines of water depth process lines during water rising and falling; (f) change of water depth envelope line under different rainfall intensity (R), Manning's coefficient (n_m), and slope gradient (S).

4 Grid-shaped cross-distribution of the discharge/water-depth process lines during water rising and falling

To obtain further insights into the causes of the formation of the water-rising limb and the water-falling limb of the water storage ratio curve, the ratio of discharge (i.e., the ratio of the total outflows (Q_out) to the total inflows (Q_in)) and the water depth (h) along the slope are discussed in Fig. 3a and d, respectively. Results indicate that there is an envelope line that controls the distribution of the discharge and water depth along the slope, respectively. The discharge envelope line is a straight line with a slope of 1 % (Fig. 3a), while the water depth envelope line is a nonlinear curve controlled by a power function of the general form h=kx^a (Fig. 3d). This means that if the duration of rainfall with a constant intensity is long enough, the catchment system will eventually reach an equilibrium state between inflow and outflow.

On the other hand, the process lines of discharge and water depth during water rising and falling present a grid-shaped cross-distribution (Fig. 3a and d). Similarly, from the view of the gradient of the discharge and water depth process lines during water rising and falling, the discharge gradient curves (Fig. 3b) and the water depth gradient curves (Fig. 3e) also present a grid-shaped cross-distribution during water rising and falling, which might be the cause of the looped rating curve (Fig. 3c), i.e., higher discharges for the rising limb (Q_u) than for the recession limb (Q_f) at the same stage (Petersen-Øverleir, 2006). After fitting the value of parameter k and a under different rainfall intensity (R), Manning's coefficient (n_m), and slope gradient (S) conditions (Fig. 3f), it is found that the parameter a is a constant, while the change of parameter k is positively correlated with the change of rainfall intensity (R) and Manning's coefficient (n_m), but it is negatively correlated with the change of slope gradient (S).

Based on the water storage ratio curve, a hydrological–hydrodynamic integrated model, namely, the distributed runoff model (DRM), is established with the governing equations in Eq. (3). To check the effectiveness and applicability of the DRM, a comparative analysis of the numerical results obtained from the DRM and the DW model is implemented. We found that the DRM quickly reproduces the calculation results of the time-consuming DW model under different rainfall intensities (Fig. 4a and b), different Manning's coefficients (Fig. 4c), and different slope gradients (Fig. 4d). meaning that the water storage ratio curve will provide new ideas for simulation and early warning of floods. In addition, due to the governing equations of the DRM being ordinary differential equations (ODEs), the computational efficiency of the DRM is much higher than the DW model, which is governed by nonlinear partial differential equations (PDEs). More attention should be paid to the determination of the nonlinear relationship of the water storage ratio curve under different geographical scenarios, which will be beneficial to the proposal of more efficient flood-forecasting methods or early-warning systems.

Figure 4 Comparative analyses of discharge calculated by DW and the DRM under designed rainfall. (a) controlled group; (b) compared with (a), only the rainfall intensity is changed; (c) compared with (a), rainfall intensity and Manning coefficient are changed; (d) compared with (a), rainfall intensity and slope gradient are changed.

5 Validation of the DRM by considering infiltration calculated by the Horton infiltration method

In the above section, the simulations of DW and the DRM are based on an impermeable conceptual slope model as shown in Fig. 1c. After considering infiltration in DW and the DRM, Eqs. (2) and (3) become

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial h}{\partial t}}-\mathrm{\nabla }\left({\displaystyle \frac{h^{\frac{\mathrm{5}}{\mathrm{3}}}}{n_m\sqrt{\left|S\right|}}}\mathrm{\nabla }\left(h+z\right)\right)=R-I;\text{(5)}& {\displaystyle }\left\{\begin{array}{l}\frac{\mathrm{d}H}{\mathrm{d}t}=R-q-I,\\ H=\mathit{\eta }h=\mathit{\eta }{\left(\frac{n_m}{\sqrt{S}}\right)}^{\mathrm{0.6}}{q}^{\mathrm{0.6}}{\left(\frac{A}{B}\right)}^{\mathrm{0.6}},\end{array}\right.\end{array}$$

Infiltration (I) is calculated by Horton's infiltration model (Horton, 1933), which suggests an exponential equation for modeling the soil infiltration capacity f_p (m s⁻¹):

$$\begin{array}{}\text{(6)}& f_p\left(t\right)=f_{\mathrm{c}}+(f_{\mathrm{0}}-f_{\mathrm{c}})e^{-kt},\end{array}$$

where f₀ is the initial infiltration capacity (m s⁻¹), f_c is the final infiltration capacity (m s⁻¹), and k represents the rate of decrease in the capacity (s⁻¹). The infiltration parameter sets are listed in Table 1.

Table 1 Infiltration parameter sets.

A rainfall event begins with a weak precipitation intensity. When the rainfall intensity is less than the infiltration capacity, all the rainwater will infiltrate into the soil. When the rainfall intensity exceeds the soil infiltration capacity, the surface water is generated, and Horton law (Eq. 6) applies:

$$\begin{array}{}\text{(7)}& I=\left\{\begin{array}{l}R\left(t\right)\\ f_p\left(t\right)\end{array}\right.\begin{array}{l}\text{if}R\left(t\right)\le f_p\left(t\right),\\ \text{if}R\left(t\right)>f_p\left(t\right).\end{array}\end{array}$$

Results of outlet discharge (Q) and runoff volume (ROV) calculated by DW and DRM are compared with the reference results adopted from Fernández-Pato et al. (2016) as shown in Fig. 5. Figure 5a shows the comparison of results under a uniform design rainfall. In this case, the rain volume is 75 000 m³ with a duration of 250 min. Figure 5b shows the comparison of results under a non-uniform rainfall. Rain volume is 75 000 m³ with a duration of 250 min. From Fig. 5, it can be recognized that after considering infiltration, except that the calculation results of the DRM are a little small at the end stage of rainfall, the calculation results of the DRM are still highly consistent with the calculation results of the DW model and reference results adopted from Fernández-Pato et al. (2016).

Figure 5 Outlet discharge (Q) and runoff volume (ROV) calculated by DW and the DRM vs. reference results adopted from Fernández-Pato et al. (2016).

6 Fluctuation of water storage ratio under natural rainfall conditions

After implementing a real rainfall event in the impermeable conceptual slope model (Fig. 1c), the change in the water storage ratio is calculated as shown in Fig. 6. Rainfall data were recorded from 9 August 2022 00:00 JST (Japan standard time, UTC+9) to 10 August 2022 00:00 JST in Aomori Prefecture, Japan, and from 29 August 2016 01:00 JST to 31 August 2016 09:00 JST in Nissho Pass, Japan (https://www.data.jma.go.jp, last access: 10 August 2023). The total simulation times are 30 and 56 h, respectively. Results show that in addition to the fluctuations of water storage ratio in the beginning and end stages of rainfall, there are mainly 10 fluctuation periods of water storage ratio during the rainfall duration stage, identified as 1^#, 2^#, 3^#, 4^#, and 5^# in Fig. 6a and 6^#, 7^#, 8^#, 9^#, and 10^# in Fig. 6b. The fluctuations are found to be mainly caused by weak rainfall (i.e., rainfall intensity near 5.0 mm h⁻¹) as pointed out by the red arrows in Fig. 6a and b. The magnitude of the fluctuations appears to be positively correlated with the difference between rainfall intensity and 5.0 mm h⁻¹. When the rainfall intensity continues to be greater than 5.0 mm h⁻¹, the fluctuation of the water storage ratio is not obvious. The water storage ratio is stable near the steady limit, even if there is heavy rainfall during this period.

Figure 6 The fluctuation of water storage ratio and the effectiveness of the DRM in natural rainfall events. (a) Aomori Prefecture; (b) Nissho Pass.

Besides, the fluctuations of the water storage ratio can be divided into three modes: Mode I identified as the inverse S-shape type during the rainfall beginning stage (Fig. 7a), Mode II identified as a wave type during the weak rainfall duration stage (Fig. 7b), and Mode III identified as a checkmark type during the rainfall end stage (Fig. 7c). Among them, Mode I describes how water storage ratio drops from the upper limit to the steady limit in an inverse S-shape. Mode II represents the water storage fluctuations around the steady limit. Mode III happens when the water storage ratio first drops from the steady limit to the lower limit and then rises to the upper limit. This means that the certainty of the fluctuation modes will provide the possibility for quantitative analysis of the fluctuation of the water storage ratio induced by the change in the rainfall intensity.

Figure 7 Three kinds of water storage ratio fluctuation modes in natural rainfall events. (a) Mode I during the rainfall beginning stage; (b) Mode II during the weak rainfall duration stage; (c) Mode III during the rainfall end stage.

Figure 8a and b show the simulation results of discharge calculated by the DRM and DW model using the rainfall data recorded in Aomori Prefecture and Nissho Pass, Japan, respectively. Results suggest that after the determination of the water storage ratio fluctuations, the calculation results of the DRM are in good agreement with those of the DW model, meaning that the DRM provides a new and more effective theoretical scheme for flood prediction.

Figure 8 Time-dependent discharge calculated by the DRM and DW model. (a) Aomori Prefecture; (b) Nissho Pass.

7 Discussions and conclusions

Based on a conceptual slope model, numerical simulations of the rainfall–runoff process are performed by using the diffusion wave (DW) approximation of SWEs. A plume-shaped nonlinear relationship between water storage and outflow, defined as the water storage ratio, is found between the inside average water depth and the outlet water depth in a catchment. The water storage ratio is controlled by three limits, namely, upper limit, steady limit, and lower limit with values of approximately 1.0, 0.625, and 0.4125, respectively. Under the control of the three limits, meteorological, vegetation, and terrain conditions only affect the size of the plume without changing its shape. The regular curve shape of the water storage ratio provides the possibility to construct a correlation between the water storage in the catchment area and the outlet discharge.

Based on the water storage ratio, a hydrological–hydrodynamic integrated model (DRM) is established, which shows high calculation accuracy and computational efficiency. This is because the governing equations of the DRM are ordinary differential equations (ODEs), which are much easier to solve than nonlinear partial differential equations (PDEs). However, the calculations of the DRM and DW only involve the confluence part of surface water and infiltration, while the interbasin groundwater flow as inputs to the watershed (exfiltration) and evaporation are not considered. This is inconsistent with the real rainfall–runoff process in the watershed and may lead to deviations in the calculation results. Therefore, the flow exchange between surface water and groundwater during the existence and extinction of runoff also needs to be further realized by establishing a dynamic coupling model of surface water and groundwater.

In addition, the water storage and discharge are limited to envelope lines, and the discharge/water-depth process lines during water rising and falling showed a grid-shaped distribution, which might be the cause of the looped rating curve, i.e., higher discharges for the rising limb than for the recession limb at the same stage. Rainfall, especially weak rainfall (i.e., rainfall intensity less than 5.0 mm h⁻¹), significantly affects the fluctuations of the water storage ratio. The fluctuations of water storage ratio during a real rainfall event can be divided into three modes: Mode I identified as inverse S-shape type during the rainfall beginning stage, Mode II identified as Wave type during weak rainfall duration stage, and Mode II identified as checkmark type during rainfall end stage. It is worth noting that a qualitative determination of the three fluctuation modes of water storage ratio during rainfall events is obtained, but the quantitative analysis still needs to be further carried out in the future.

The findings in this study provide a key to establishing a simpler prediction model for flash floods. The water storage ratio has been proven to be effective in improving the effectiveness and efficiency of flood forecasting. Therefore, the determination of the nonlinear relationship of the water storage ratio curve under different geographical scenarios will provide new ideas for simulation and early warning of flash floods.