A time-varying distributed unit hydrograph method considering soil moisture

The distributed unit hydrograph (DUH) method has been widely used for flow routing in a watershed because it adequately characterizes the underlying surface characteristics and varying rainfall intensity. Fundamental to the calculation of DUH is flow velocity. However, the currently used velocity formula assumes a global equilibrium of the watershed and ignores the impact of time-varying soil moisture content on flow velocity, which thus leads to a larger flow velocity. The objective of this study was to identify a soil moisture content factor, which, based on the tension water storage capacity curve, was derived to investigate the response of DUH to soil moisture content in unsaturated areas. Thus, an improved distributed unit hydrograph, based on time-varying soil moisture content, was obtained. The proposed DUH considered the impact of both time-varying rainfall intensity and soil moisture content on flow velocity, assuming the watershed to be not in equilibrium but varying with soil moisture. The Qin River basin and Longhu River basin were selected as two case studies, and the synthetic unit hydrograph (SUH), the time-varying distributed unit hydrograph (TDUH) and the current DUH methods were compared with the proposed method. Then, the influence of time-varying soil moisture content on flow velocity and flow routing was evaluated, and results showed that the proposed method performed the best among the four methods. The shape and duration of the unit hydrograph (UH) were mainly related to the soil moisture content at the initial stage of a rainstorm, and when the watershed was approximately saturated, the grid flow velocity was mainly dominated by excess rainfall. The proposed method can be used for the watersheds with sparse gauging stations and limited observed rainfall and runoff data.

1 Introduction

Flow routing is an essential component of a hydrological model, whose accuracy directly affects runoff prediction and forecasting. Different types of flow routing techniques are available, such as hydraulic and hydrologic methods (Akram et al., 2014). Since hydraulic methods are usually computationally intensive, hydrologic methods are widely used all over the world. The unit hydrograph (UH), proposed by Sherman (1932), is one of the methods most widely used in the development of flood prediction and warning systems for gauged basins with observed rainfall and runoff data (Singh et al., 2014). However, the UH method has inherent problems, such as areal lumping of catchment and rainfall characteristics, as well as the utilization of linear system theory (Singh, 1988; James and Johanson, 1999). Moreover, current routing methods usually require numerous rainfall and runoff data. For watersheds with sparse gauging stations, it is difficult to develop an adequate relationship between physical watershed characteristics and unit hydrograph shape. The unit hydrograph estimation in small and ungauged basins is still a challenge in hydrological studies (Petroselli and Grimaldi, 2018).

The UH, which is a surface runoff hydrograph resulting from one unit of rainfall excess uniformly distributed spatially and temporally over the watershed for the specified rainfall excess duration (Chow 1964), can be categorized into four major types (Singh, 1988), including traditional, probability-based, conceptual and geomorphologic methods (Bhuyan et al., 2015).

Synthetic UH methods establish the relationships between watershed characteristic for describing the UH (e.g., peak flow, time to peak and time base) and parameters used to describe the basin. Snyder (1938), Mockus (1957) and U.S. Soil Conservation Service (SCS) (1972) proposed some of these methods, which are still used. The disadvantages of these methods are that they do not yield adequately satisfactory results, and their application to practical engineering problems is tedious and cumbersome (Nigussie et al., 2016).

Since most UHs have rising limbs steeper than their receding sides, and their shape resembles typical probability distribution functions (PDFs), many PDFs have been used for the derivation of UHs. The difficulty of this method is that the PDFs are diverse, and their parameters depend on numerous hydrological data (Bhuyan et al., 2015).

Conceptual methods are another technique for deriving UHs. Nash (1957) proposed a conceptual model composed of n linear reservoirs connected in series (or a cascade) with the same storage coefficient K for the derivation of the instantaneous unit hydrograph (IUH). Dooge (1959) proposed a generalized IUH based on linear reservoirs, linear channels, and time–area concentration diagrams. Bhunya et al. (2005) and Singh et al. (2007) represented a hybrid method and an extended hybrid method based on a linear reservoir. Singh (2015) proposed a new simple two-parameter IUH with conceptual and physical justification. Khaleghi et al. (2018) suggested a new conceptual model, namely, the inter-connected linear reservoir model (ICLRM), which, however, neglects the impact of uneven basin surface on the UH.

Rodriguez-Iturbe (1979) proposed a geomorphologic instantaneous unit hydrograph (GIUH) method, which couples the hydrologic characteristics of a catchment with geomorphologic parameters (Singh, 1988; Kumar et al., 2007). In this method, the IUH corresponds to the probability density function of travel times from the locations of runoff production to the watershed outlet (Gupta et al., 1980; Singh, 1988). With the development of digital elevation models (DEMs) and geographic information system (GIS) technology, the width-function-based geomorphological IUH method has been formulated. However, in capacity it is unable to properly account (i.e., to respect the geometry) for the spatial distribution of rainfall (Rigon et al., 2016).

The UH method assumes the watershed response to be linear and time-invariant and rainfall to be spatially homogeneous. Contrary to the linearity assumption, basins have been shown to exhibit nonlinearity in the transformation of excess rainfall to storm flow (Bunster et al., 2019). For a small watershed, Minshall (1960) showed that significantly different UHs were produced by different rainfall intensities. To cope with this nonlinearity, Rodríguez-Iturbe et al. (1982) extended the GIUH to the geomorphoclimatic IUH (GcIUH) by incorporating excess rainfall intensity. Lee et al. (2008) proposed a variable kinematic-wave GIUH accounting for time-varying rainfall intensity, which may be applicable to ungauged catchments that are influenced by high-intensity rainfall. Du et al. (2009) proposed a GIS-based routing method to simulate storm runoff with the consideration of spatial and temporal variability of runoff generation and flow routing through hillslope and river network. A similar work was done by Muzik (1996), Gironás et al. (2009), and Bunster et al. (2019).

The traditional, probabilistic, conceptual, and geomorphologic methods have not been able to fully consider the geomorphic characteristics of the watershed and incorporate time-varying rainfall intensity.

The spatially distributed unit hydrograph (DUH) method conceptualizes that the unit hydrograph can be derived from the time–area curve of the watershed using the S-curve method (Muzik, 1996). It is a type of geomorphoclimatic unit hydrograph, since its derivation considers watershed geomorphology (Du et al., 2009) and spatially distributed flow celerity, and temporally varying excess rainfall intensities can be considered in DUH (Bunster et al., 2019). In this method, the travel time of each grid cell can be calculated by dividing the travel distance of a cell to the next cell by the velocity of flow generated in that cell (Paul et al., 2018). The travel time is then summed along the flow path to obtain the total travel time from each cell to the outlet. The DUH is thus derived using the distribution of travel time from all grid cells in a watershed (Bunster et al., 2019). Some DUH methods assumed a time-invariant travel time field and ignored the dependence of travel time on excess rainfall intensity (Melesse and Graham, 2004; Noto and La Loggia, 2007; Gibbs et al., 2010), while others suggested various UHs corresponding to different storm events, namely the time-varying distributed unit hydrograph (TDUH) (Martinez et al., 2002; Sarangi et al., 2007; Du et al., 2009). Compared to the fully distributed methods based on the momentum equation, the DUH is a more efficient method because it allows for the use of distributed terrain information and is an alternative to semi-distributed and fully distributed methods for rainfall-runoff modeling (Bunster et al., 2019).

Besides excess rainfall intensity, the upstream contributions to the travel time estimation have also been considered in the time-varying DUH method. For instance, Maidment et al. (1996) defined the velocity in the cell as a function of the contributing area to take into account the velocity increase observed downstream in river systems (Gironás et al., 2009). Gad (2014) applied a grid-based method using stream power to relate flow velocity to the hydrologic parameters of the upstream watershed area. Similar work was done by Saghafian and Julien (1995), Bhattacharya et al. (2012) and Chinh et al. (2013). A major drawback of this method is the assumption that the watershed is near global equilibrium. Bunster et al. (2019) developed a spatially time-varying DUH method that accounts for dynamic upstream contributions and characterized the temporal behavior of upstream contributions and their impact on travel times in the basin. However, this time-varying DUH also assumed that equilibrium in each individual grid cell was reached before the end of the rainfall excess pulse. When continuous excess rainfall accrues in a watershed, the soil moisture content and surface runoff increase, and the infiltration rate decreases, leading to an acceleration of flow routing velocity, until the entire basin is saturated, and the routing velocity reaches its maximum. This assumption of equilibrium globally or in grid cells yields faster travel flow velocities, smaller travel time and higher peak discharge. However, these approximations neglect the impact of dynamic changes of soil moisture exchange and water storage in unsaturated regions.

The objective of this study was therefore to propose a time-varying distributed unit hydrograph method for runoff routing that accounts for dynamic rainfall intensity and soil moisture content based on the Xinanjiang (XAJ) model, namely the time-varying distributed unit hydrograph considering soil moisture content (TDUH-MC). The main contributions of the present study are as follows. First, a soil moisture content proportional factor in the unsaturated area was identified and expressed based on the Pareto distribution function. Second, the travel time function based on the kinematic wave theory was modified by considering the soil moisture content proportional factor. Besides rainfall intensity, the influence of time-varying soil moisture storage on flow velocity in the watershed was considered, where runoff generation was dominated by the saturation-excess mechanism. Finally, the Qin River basin and Longhu River basin in Guangdong Province, China, were selected as two case studies. The flow forecast method mainly consisted of the calculation of excess rainfall and the derivation of DUH. A new routing method was developed to incorporate the dynamic changes of soil moisture content and rainfall intensity, and the XAJ model was adopted to calculate excess rainfall. The synthetic unit hydrograph (SUH), DUH and TDUH methods were compared with the TDUH-MC method, and sensitivity analysis of parameters was conducted.

2 Improvement of flow routing method

2.1 Calculation of flow velocity considering time-varying soil moisture content

The DUH relies on the computation of travel time in the basin. Grimaldi et al. (2010) found that the Soil Conservation Service (SCS) formula, given by Eq. (1), can be used to adequately define the basin flow time. This formula was also used by NRCS (1997) and Grimaldi et al. (2012), but this formula is time-invariant, and the time-varying rainfall intensity should be considered, as given by Eq. (2), which was used by Wong (1995), Muzik (1996), Gironas et al. (2009), Du et al. (2009) and Kong et al. (2019).

$$\begin{array}{}\text{(1)}& {\displaystyle }V=k\cdot S^{\frac{\mathrm{1}}{\mathrm{2}}}\text{(2)}& {\displaystyle }V=k\cdot S^{\frac{\mathrm{1}}{\mathrm{2}}}\cdot {\left({\displaystyle \frac{I_t}{I_c}}\right)}^{\frac{\mathrm{2}}{\mathrm{5}}},\end{array}$$

where V (m s⁻¹) is the flow velocity, k (m s⁻¹) is the land use or flow type coefficient, S (m m⁻¹) is the slope of the grid cell, I_t (mm h⁻¹) represents the excess rainfall intensity at time t and I_c (mm h⁻¹) represents the reference excess rainfall intensity of the basin.

These formulas assume that equilibrium in individual grid cell can be reached before the end of the rainfall excess (Bunster et al., 2019), which leads to larger flow velocity, shorter travel time and higher peak discharge. Actually, the hillslope flow velocity in each grid is related to soil moisture content. Fast subsurface velocities and quick runoff responses to precipitation have been observed on many hillslopes (Hutchinson and Moore, 2000; Peters et al., 1995; Tani, 1997). The exact mechanisms that cause water to move through the preferential flow path network are not well quantified, but it is often assumed that saturated soil provides the connection between preferential features (Sidle et al., 2001; Steenhuis et al., 1988). Studies have also shown that antecedent moisture condition, precipitation intensity, precipitation amount, topography and so on play a significant role in flow configuration (Sidle et al., 2000; Tsuboyama et al., 1994; Anderson et al., 2009).

To that end, a soil moisture factor θ_t was introduced to characterize the soil moisture content in unsaturated areas. Because the flow velocity will reach its maximum value when the entire basin is saturated, this new factor (θ_t) was added to the current time-varying flow velocity formula as

$$\begin{array}{}\text{(3)}& V=k\cdot S^{\frac{\mathrm{1}}{\mathrm{2}}}\cdot {\left({\displaystyle \frac{I_t}{I_c}}\right)}^{\frac{\mathrm{2}}{\mathrm{5}}}\cdot {\left({\mathit{\theta }}_t\right)}^{\mathit{\gamma }},\end{array}$$

where θ_t (unitless) represents the soil moisture content of unsaturated areas at time t, and γ (unitless) is an exponent smaller than unity, which represents the nonlinear relationship between soil moisture content and flow velocity.

Factor θ_t was defined as the ratio of w_t and w_max,t, which is expressed by

$$\begin{array}{}\text{(4)}& {\mathit{\theta }}_t={\displaystyle \frac{w_t}{w_{max,t}}},\end{array}$$

where w_t (mm) represents the mean tension water storage of the unsaturated region, and w_max,t (mm) represents the maximum tension water storage of the unsaturated region at time t.

Specifically, w_t and w_max,t were calculated based on the Pareto distribution function in this study. The Pareto distribution function has mostly been used to express the spatial variability of soil moisture capacity (Moore, 1985). As shown in Fig. 1, the area below the curve represents the mean tension water capacity of the entire basin.

Figure 1 Watershed storage capacity curve.

For the tension water storage capacity curve, the specific formula is given by

$$\begin{array}{}\text{(5)}& \mathit{\alpha }=\mathrm{1}-{\left(\mathrm{1}-{\displaystyle \frac{\mathrm{WM}}{\mathrm{WMM}}}\right)}^b,\end{array}$$

where α (unitless) represents the proportion of the basin area where the tension water capacity is less than or equal to the value of the ordinate WM (mm). The tension water capacity at a point, WM, varies from 0 to a maximum WMM (mm) according to Eq. (5).

Since the soil moisture content in a basin varies with time, the state of the catchment at any time t can be represented by a point x(α_t,WM_t) on the curved line of Fig. 1 (Zhao, 1992). The area to the right and below the point x is proportional to the areal mean tension water storage (not capacity). Thus, WM_t, the ordinate of the point x, represents the tension water storage capacity in the basin at time t, w_t (mm) can be assumed to represent the mean tension water storage of the unsaturated region and w_max,t (mm) represents the maximum tension water storage of the unsaturated region at time t. Their expressions are given by

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathit{\alpha }}_t=\mathrm{1}-{\left(\mathrm{1}-{\displaystyle \frac{{\mathrm{WM}}_t}{\mathrm{WMM}}}\right)}^b\text{(7)}& {\displaystyle }{w}_t=\left(\mathrm{1}-{\mathit{\alpha }}_t\right)\cdot {\mathrm{WM}}_t\text{(8)}& {\displaystyle }{w}_{max,t}=\underset{{\mathit{\alpha }}_t}{\overset{\mathrm{1}}{\int }}\mathrm{WMM}\left[\mathrm{1}-{\left(\mathrm{1}-\mathit{\alpha }\right)}^{\frac{\mathrm{1}}{b}}\right]\mathrm{d}\mathit{\alpha }.\end{array}$$

Combining Eqs. (4), (7) and (8), the soil moisture content can be written as

$$\begin{array}{}\text{(9)}& {\mathit{\theta }}_t={\displaystyle \frac{w_t}{w_{\max ,t}}}={\displaystyle \frac{\left(\mathrm{1}-{\mathit{\alpha }}_t\right)\cdot {\mathrm{WM}}_t}{\underset{{\mathit{\alpha }}_t}{\overset{\mathrm{1}}{\int }}\mathrm{WMM}\left[\mathrm{1}-{\left(\mathrm{1}-\mathit{\alpha }\right)}^{\frac{\mathrm{1}}{b}}\right]\mathrm{d}\mathit{\alpha }}}.\end{array}$$

Substituting Eq. (6) into (9),

$$\begin{array}{}\text{(10)}& \begin{array}{rl}{\mathit{\theta }}_t& ={\displaystyle \frac{\left(\mathrm{1}-{\mathit{\alpha }}_t\right)\cdot {\mathrm{WM}}_t}{\mathrm{WMM}\left[\mathrm{1}-{\mathit{\alpha }}_t-\frac{b}{b+\mathrm{1}}{\left(\mathrm{1}-{\mathit{\alpha }}_t\right)}^{\mathrm{1}+\frac{\mathrm{1}}{b}}\right]}}\\ & ={\displaystyle \frac{\left(b+\mathrm{1}\right){\mathrm{WM}}_t}{\mathrm{WMM}+b{\mathrm{WM}}_t}}.\end{array}\end{array}$$

It can be seen from Eq. (10) that as rainfall continues, the soil moisture content in the unsaturated area continues to increase, whereas the non-runoff area continues to decrease. The range of θ_t is (0, 1], and with the gradual increase of soil moisture, θ_t tends to 1.

2.2 Calculation of runoff routing based on DUH

The GIS-derived DUH method was employed for runoff routing calculations, which allowed the velocity to be calculated on a grid cell basis over the watershed. The DUH routing method is a semi-analytical form of the width-function-based IUH enumerated by Rigon et al. (2016). The DUH has been used for small ungauged basins. To remove the linearity assumption, fully distributed models use routing methods which are usually computationally intensive because they solve the St. Venant equations (Bunster et al., 2019), so they are usually limited to small basins. Therefore, the DUH method is an alternative method that allows for the use of distributed information in a much more efficient manner, and we applied it to different sizes of watersheds.

The core of the DUH method is to equate the probability density function of time at which the rainfall flows to the basin outlet to form the instantaneous unit hydrograph, in which the time–area relationship is derived using the velocity field with spatial distribution characteristics. The traditional DUH method can route the time-variant spatially distributed rainfall to the watershed outlet, but such a method is a lumped linear model of watershed response (Grimaldi et al., 2010). The schematic diagram of the DUH method is shown in Fig. 2.

Figure 2 Schematic diagram of the DUH method considering time-varying rainfall intensity and soil moisture content, in which Eqs. (1), (2) and (3) are the time-invariant flow velocity, time-varying flow velocity considering excess rainfall intensity and time-varying flow velocity considering both excess rainfall intensity and soil moisture content. The unit hydrograph derived from the three flow velocity equations corresponds to DUH, TDUH and the TDUH-MC method respectively.

The steps of the DUH method are summarized as follows:

1. The drainage network based on the advanced DEM preprocessing method is identified. More details can be found in Grimaldi et al. (2012).

2. Estimate the flow path, which is measured for each grid cell along the flow directions to the basin outlet.

3. Calculate the flow velocity based on watershed characteristics and the spatial–temporal distribution characteristics of rainfall. Several flow velocity formulas are commonly used for deriving the spatially distributed unit hydrograph, such as Manning's formula (Chow et al., 1988), the SCS formula (Haan et al., 1994), the Darcy–Weisbach formula (Katz et al., 1995), and the uniform flow equation by Maidment et al. (1996).

4. To compute the total travel time τ_i of flow from each cell i to the outlet, we added travel times along the R_i cells belonging to the flow path that starts at that cell, given by Eq. (11) (Muzik, 1996). The travel time for each grid cell can be calculated by Eq. (12):

   $$\begin{array}{}\text{(11)}& {\displaystyle }{\mathit{\tau }}_i=\sum _{i\in R_i}\mathrm{\Delta }{\mathit{\tau }}_i\text{(12)}& {\displaystyle }\mathrm{\Delta }{\mathit{\tau }}_i={\displaystyle \frac{L_i}{V}}\mathrm{or}\mathrm{\Delta }{\mathit{\tau }}_i={\displaystyle \frac{\sqrt{\mathrm{2}}{L}_i}{V}},\end{array}$$

   where Δτ_i is the retention time in grid cell i, τ_i is the total travel time along the flow path in grid cell i and L_i is the grid cell size. Travel length in a specific grid cell is the cell size L_i when the rasterized flow is flowing along the edges of the grid, whereas the travel length is $\sqrt{\mathrm{2}}{L}_i$ when it is flowing diagonally.

5. Develop a cumulative travel time map of the watershed based on cell-by-cell estimates for hillslope velocities. The cumulative travel time map is further divided into isochrones, which can be used to generate a time–area curve and the resulting unit hydrograph (Kilgore, 1997).

Figure 3 Schematic diagram of the XAJ model.

Figure 4 Locations of meteorological stations and gauging stations in the Qin River basin.

3 Calculation of runoff generation

The Xinanjiang (XAJ) model was used for the calculation of excess rainfall in this study. It is a conceptual hydrologic model proposed by Zhao et al. (1980) for flood forecasts in the Xinan River basin. The XAJ model has been widely used in humid and semi-humid watersheds all over the world (Zhao, 1992). It mainly consists of four modules, namely the evapotranspiration module, runoff generation module, runoff partition module and runoff routing module (Zhou et al., 2019). Usually, a large watershed is divided into several sub-basins to capture the spatial variability of underlying surface, precipitation, and evaporation. In each sub-basin, the inputs of the XAJ model are the average areal rainfall as well as evaporation, and the output is streamflow. The schematic diagram of the XAJ model is shown in Fig. 3.

First, for the evapotranspiration module, the soil profile of each sub-basin is divided into three layers, the upper, lower and deeper layers, and only when water in the layer above it has been exhausted does evaporation from the next layer occur. Second, for runoff generation in the XAJ model, a catchment is divided into two parts by the percentage of impervious and saturated areas, namely pervious and impervious areas, respectively. Since the soil moisture deficit is heterogeneous, runoff distribution is usually nonuniform across the basin. Thus, a storage capacity curve was adopted by the XAJ model to accommodate the nonuniformity of soil moisture deficit or the tension water capacity distribution. Third, the runoff partition in the XAJ model divides the total runoff into three components by a free reservoir, which consists of surface runoff (RS), interflow runoff (RI) and groundwater runoff (RG). More details can be found in Zhao et al. (1980).

Finally, the SUH was selected as the runoff routing approach in the XAJ model. Specifically, the Nash instantaneous unit hydrograph model (Nash, 1957) was used to derive the SUH in this study. For the Nash IUH model, a catchment was assumed to be made up of a series of n identical linear reservoirs, each with the same storage constant K. The magnitudes of n and K were estimated based on the observed excess rainfall hyetograph and corresponding direct runoff hydrograph using the method of moments. Details can be found in Singh (1988) and Chow et al. (1988).

The Muskingum method was employed to produce streamflow from each sub-basin to the outlet of the entire basin. For the SUH, the basin was taken as a whole. The parameters of the Muskingum method, including the Muskingum time constant KE and Muskingum weighting factor XE, were calibrated with those of the XAJ model. The Shuffled Complex Evolution Algorithm (SCE-UA) method was used to calibrate the parameters of XAJ model (Chu et al., 2009; Moghaddam et al., 2016). For the DUH, the basin was divided into several sub-basins. Since natural rivers are multiple inflow–single outflow runoff systems with different travel times from the sub-basins to the outlet, we adopted the physical and numerical principles established by Cunge to calculate the routing parameters of the Muskingum method, which is suitable for ungauged watersheds (Ponce et al., 1996). The Muskingum parameters for each sub-basin were determined based on flow and channel characteristics, such as the top width of the river, wave celerity, reach length and reach slope, as described in Chow (1959) and Wilson and Ruffin (1988).

4 Study area and data

The Qin River basin and Longhu River basin were selected as two case study watersheds. One is a large watershed, and the other is a small watershed. The applicability of the TDUH-MC method to different size watersheds was verified, and parameter sensitivity analysis was done to evaluate the performance of the TDUH-MC method (Chen et al., 2022a).

The Qin River is a tributary of the Mei River, which originates from Guangdong Province, China. The river is 91 km long, with a basin area of 1578 km². The mean slope of the basin is 1.1 ‰. There are 21 meteorological stations and 1 flow station (the Jianshan station) in the basin, as shown in Fig. 4. Using the DEM data of the Qin River basin, the whole basin was divided into nine sub-basins, namely sub-basins 1–9 from upstream to downstream, as shown in Fig. 5. Details of each sub-basin are given in Table 1.

Figure 5 Sub-basins of the Qin River basin. (Note that the satellite images for the study area are available at http://www.gscloud.cn, last access: 18 November 2020.)

Table 1 Detailed information of each sub-basin.

The Longhu River basin is a small watershed, which has a drainage area of 102.7 km², located in Guangdong Province, China. The length of the river is 17.4 km.

The rainfall and evaporation data from meteorological stations for the two basins were collected from 1959 to 2018, and the simultaneous hourly runoff data for the Jianshan Station and Longhu Station were collected as well. A total of 64 isolated storms with the observed runoff responses from 1959 to 2018 were selected to calibrate and verify the established model, of which 35 events were collected from the Qin River basin and 29 from the Longhu River basin. In total, 25 and 23 flow events were used for model calibration in the Qin and Longhu River basins respectively, and 10 and 6 flow events were used for model validation in the two basins.

The statistics of flow events used for model calibration and validation are shown in Fig. 6. The average peak flows of the two basins were 1311 and 118 m³ s⁻¹, and the average flood durations were about 50 and 13 h, respectively. The antecedent precipitation was calculated based on the daily recession coefficient of water storage in the basin.

Figure 6 Statistics of flow events used for model calibration and validation.

Figure 7 Schematic of the calibration procedure.

5 Results and discussion

5.1 Calibration of parameters

5.1.1 Model calibration

The runoff generation model (XAJ model) and the routing model were calibrated separately in this study. First, the SUH and several distributed unit hydrographs (DUH, TDUH and MC-TDUH) were derived. Second, the Shuffled Complex Evolution Algorithm (SCE-UA) method, developed by the University of Arizona (Duan et al., 1992), was used to optimize the XAJ model parameters (Vrugt et al., 2006; Beskow et al., 2011; Zhou et al., 2019). The SUH was selected as the runoff routing method. As the SUH was derived from observed rainfall and runoff, the flow routing model corrected some inconsistencies of the hydrological model. Therefore, the parameters of excess runoff were calibrated. Third, the performances of XAJ + SUH and XAJ + DUHs (DUH, TDUH and MC-TDUH) were compared. Since the XAJ model parameters were determined by combining with SUH routing method, this calibration method would be more inclined to optimize the performance of XAJ + SUH model. When combined with other confluence models, the accuracy of results may be affected to some extent. The schematic of the calibration procedure is given below.

The steps of parameter calibration can be summarized as follows:

1. The XAJ model was used to calculate the excess rainfall, in which the SUH derived from observed runoff was selected as the runoff routing method. The SCE-UA method was used to optimize the XAJ model parameters in this study. In total, 25 and 23 flow events in the Qin River basin and Longhu River basin were used for the calibration of the XAJ + SUH model.

2. The SUH was derived using 25 and 23 flow events in the Qin River basin and Longhu River basin, respectively. The DUH, TDUH and MC-TDUH were derived, based on physical characteristics and rainfall intensities of the watersheds. The parameters' determination method is given in Sect. 5.1.3.

3. Since the objective of this study was to propose a new flow routing method, the runoff production model and its parameters were not changed in order to discuss the performance of flow routing models. The XAJ model with calibrated parameters in Step (1) and DUH, TDUH as well as MC-TDUH determined in Step (2) were used for the validation period. A total of 10 and 6 flow events of the two basins were then used for the validation of the XAJ + (SUH, DUH, TDUH and MC-TDUH) model.

The Nash–Sutcliffe efficiency (E_NS) (Nash and Sutcliffe, 1970; Chen et al., 2015), the Kling–Gupta efficiency (E_KG) (Gupta et al., 2009) and the root-mean-squared error to standard deviation ratio (R_SR) were chosen as criteria. Moreover, the new aggregated objective function (Brunner et al., 2021) targeted at optimizing flow characteristics was composed of these three metrics, in which E_KG focuses on high flows (Mizukami et al., 2019), log(E_NS) emphasizes low flows and R_SR quantifies volume errors. A similar method has been used by Chen et al. (2022b, c). Three metrics and the aggregated objective function are expressed by

$$\begin{array}{}\text{(13)}& {\displaystyle }{E}_{\mathrm{NS}}=\mathrm{1}-{\displaystyle \frac{{\sum }_{t=\mathrm{1}}^T\left|Q_{\mathrm{s}}^t-Q_{\mathrm{o}}^t\right|}{{\sum }_{t=\mathrm{1}}^T\left|Q_{\mathrm{o}}^t-\stackrel{\mathrm{‾}}{Q_{\mathrm{o}}}\right|}}\text{(14)}& {\displaystyle }{E}_{\mathrm{KG}}=\mathrm{1}-\sqrt{{\left(r-\mathrm{1}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{s}}}{{\mathit{\sigma }}_{\mathrm{o}}}}-\mathrm{1}\right)}^{\mathrm{2}}+{\left({\displaystyle \frac{{\mathit{\mu }}_{\mathrm{s}}}{{\mathit{\mu }}_{\mathrm{o}}}}-\mathrm{1}\right)}^{\mathrm{2}}}\text{(15)}& {\displaystyle }{R}_{\mathrm{SR}}=\sqrt{{\displaystyle \frac{{\sum }_{t=\mathrm{1}}^T{\left(Q_{\mathrm{o}}^t-Q_{\mathrm{s}}\right)}^{\mathrm{2}}}{{\sum }_{t=\mathrm{1}}^T{\left(Q_{\mathrm{o}}^t-\stackrel{\mathrm{‾}}{Q_{\mathrm{o}}}\right)}^{\mathrm{2}}}}}\text{(16)}& {\displaystyle }\begin{array}{rl}M& =\mathrm{0.5}\times \left(\mathrm{1}-E_{\mathrm{NS}}\right)+\mathrm{0.25}\\ & \times \left(\mathrm{1}-E_{\mathrm{KG}}\right)+\mathrm{0.15}\times \left(\mathrm{1}-\log \left(E_{\mathrm{NS}}\right)\right)+\mathrm{0.1}\times R_{\mathrm{SR}},\end{array}\end{array}$$

where $Q_{\mathrm{o}}^t$ is the observed discharge at time t, $Q_{\mathrm{s}}^t$ is the simulated discharge at time t, $\stackrel{\mathrm{‾}}{Q_{\mathrm{o}}}$ is the mean of observed discharge, T is the duration of the flow event, r is the correlation coefficient between observed and simulated floods, σ_s and σ_o are the standard deviation values for the simulated and observed responses, respectively, and μ_s and μ_o are the corresponding mean values.

Table 2 Calibrated parameters of the XAJ model.

Table 3 Specific values of k for different vegetational types.

Figure 8 (a) Slope distribution, (b) land types and (c) rasterized flow direction of the Qin River basin.

5.1.2 Calibrated parameters of runoff generation using the XAJ model

Since the Qin River basin and Longhu River basin are in a humid area of southern China, the saturation-excess mechanism with three-source runoff separation of the XAJ model was adopted to calculate excess rainfall. The initial condition of the XAJ model was considered by calculating the antecedent precipitation index before each flow event (Linsley et al. 1949). The synthetic unit hydrograph, derived by historical rainfall-runoff data, was used for flow routing in the process of model calibration. The time interval was 1 h. Several studies have shown that UH which is derived by considering antecedent soil moisture is more consistent than UH which ignores that (Yue and Hashino, 2000; Nourani et al., 2009). Therefore, the antecedent precipitation was calculated and considered in this study. In order to obtain the SUH, we defined excess rainfall and separated direct runoff and baseflow hydrographs in advance. The final SUH used for calibration is the average value deduced by multiple historical flow events. The parameter n of the Qin River basin and Longhu River basin was 4 and 3, and the parameter K for the two basins was 3.4 and 2.1, respectively. Then, the flow peak, flow volume and the occurrence time of flow peak are three main basic elements for describing the flow hydrograph, and Eq. (16) was used as the aggregated objective function. The average Nash–Sutcliffe efficiency, relative flood peak error and peak occurrence time error obtained in the calibration period of the XAJ model were 0.84 %, 10.4 % and 4.96 h, respectively, for the Qin River basin. Accordingly, for the Longhu River basin, they were 0.86 %, 8.81 % and 2.75 h respectively, indicating a good performance of the XAJ model. Detailed information on the calibrated parameters of the XAJ model is shown in Table 2.

5.1.3 Calibrated parameters of the TDUH-MC flow routing method

As mentioned in Sect. 2.2, the core of the DUH is the calculation of the grid flow velocity. As shown in Eq. (3), the parameters that needed to be calibrated were k, S, I_c and γ, in which I_c was determined using hourly mean rainfall intensity and flow forecast of the target basin. For the Qin River basin, I_c was set at 20 mm h⁻¹ because the mean rainfall intensity of multiple flows was about 20 mm h⁻¹, and this parameter was 10 mm h⁻¹ for the Longhu River basin. Additionally, parameter γ reflected the influence of soil moisture content in unsaturated regions on flow velocity. The smaller the parameter γ was, the smaller the influence of soil moisture content on the flow velocity was. When the value of γ was equal to 1, the flow velocity of grid cell was proportional to the soil moisture content factor θ_t. The parameter γ of soil moisture content was determined to be 0.5 to reflect the influence of soil moisture content on the flow velocity for the two basins. Furthermore, sensitivity analysis for this parameter was conducted in Sect. 5.6. In order to get the grid cell slope S, the slope distribution of the study areas was obtained from the DEM data of the target basin. Figure 8a plots the slope distribution of the Qin River basin. Parameter k is the velocity coefficient, which was determined based on different underlying surface types or different flow states (Ajward and Muzik, 2000). Parameter k changed with different land types, and the k values used in this study are given in Table 3. The land types of the Qin River basin are shown in Fig. 8b. Then the k values of each grid cell were determined by combining Fig. 8b and Table 3.

The grid flow velocity was calculated by Eq. (3) with the above parameter values. Then, the flow travel time was determined by Eqs. (11) and (12). It is noteworthy that the raster size of the Qin River basin was divided into 1 km × 1 km, and the rasterized flow direction of each sub-basin is shown in Fig. 8c. For the Longhu River basin, the difference was that its cell size was divided into 30 m × 30 m to evaluate the performance of the TDUH-MC method in this small watershed.

5.2 Calculation of the TDUH-MC

After determining the parameters, flow routing was calculated based on the proposed DUH, considering the time-varying soil moisture content. In order to improve the effectiveness of the routing method, the rainfall intensity and soil moisture content parameters were discretized. Then, a simplified TDUH considering time-varying soil moisture content and the TDUH were obtained in a certain range of rainfall intensities or soil moisture contents; these ranges are presented in Tables 4 and 5. To evaluate the performance of the TDUH-MC method, the traditional SUH, DUH and TDUH methods were used for comparison.

Table 4 The ratio of I_t to I_c of each period corresponds to the discrete rain intensity I_s.

Table 5 The soil moisture content θ_t of each period corresponds to the discrete soil moisture content θ_s.

Table 6 Comparison of four routing methods for the Qin River basin. The values in bold mean the result performed the best among the four methods.

The DUH without considering rainfall intensity and soil moisture was obtained using Eq. (1). Results of the DUH for each sub-basin of the Qin River basin are shown in Fig. 9. There is only one DUH for a specific sub-basin due to the simplification of the underlying surface, such as slope and land covers. The differences among the DUHs were mainly reflected in flow peaks and their occurrence times. It can also be seen from Fig. 9 that the peak of DUHs in sub-basins 4 and 6 was significantly lower than in others. The reason may be that the smaller mean slope values of sub-basins 4 and 9 lead to lower flow velocity, resulting in a lower peak of the DUH.

Figure 9 DUH for the Qin River basin.

The TDUHs corresponding to different rainfall intensities of nine sub-basins are shown in Fig. 10. It can be seen from Fig. 10 that different rainfall intensities corresponded to different TDUHs. The increased rainfall intensity led to higher peak and earlier peak occurrence time of the UH. This is because a larger rainfall intensity caused a larger flow velocity according to Eq. (2). In the practical use of TDUH, the UHs need to be selected according to rainfall intensities.

Figure 10 The TDUH for the Qin River basin. (a) Sub-basin 1. (b) Sub-basin 2. (c) Sub-basin 3. (d) Sub-basin 4. (e) Sub-basin 5. (f) Sub-basin 6. (g) Sub-basin 7. (h) Sub-basin 8. (i) Sub-basin 9.

The TDUH of each sub-basin was further divided according to the soil moisture content. The TDUHs considering soil moisture contents of sub-basin 1 are shown in Fig. 11. Obviously, under the same rainfall intensity, the soil moisture content was of great importance to the shape, peak value and duration of the TDUH. Specifically, when the proportion of soil moisture content θ_t increased, the TDUH-MC method considering soil moisture content was accompanied by a steeper rising limb, a higher peak and shorter duration. After the whole basin was saturated, the TDUH considering the soil moisture content was the same as the TDUH.

Figure 11 The TDUH considering soil moisture content for sub-basin 1 of the Qin River basin. (a) I_s=0.5. (b) I_s=1. (c) I_s=1.5. (d) I_s=2.

Similarly, the TDUHs considering the soil moisture content for the Longhu River basin are shown in Fig. 12. The grey line in Fig. 12b is the DUH, where I_s is equal to 1, and w_s is 0.85. Four grey unit hydrographs in Figs. 12a to d make up the TDUH without considering the soil moisture content.

Figure 12 The TDUH considering soil moisture content for the Longhu River basin. (a) I_s=0.5. (b) I_s=1. (c) I_s=1.5. (d) I_s=2.

Figure 13 Comparison of flow hydrographs obtained by the four methods. (a) Flow event no. 20130720. (b) Flow event no. 20130817. (c) Flow event no. 20150709. (d) Flow event no. 20160128. (e) Flow event no. 20161021. (f) Flow event no. 20180916.

5.3 Comparisons of flood routing methods

The runoff generation module of the calibrated XAJ model was used to calculate the excess rainfall, and the SUH, DUH, TDUH and improved TDUH considering soil moisture content were employed for flow routing calculations, respectively. The Muskingum parameters for each sub-basin are given in Table 1. Dozens of flow events were applied for model validation. Simulated results of the four methods for the Qin River basin are shown in Table 6. Three criteria were used for model performance evaluation, which included the Nash–Sutcliffe efficiency (E_NS), the ratio between the simulated and observed peak discharges ($Q_{\mathrm{p}}^{\mathrm{s}}/Q_{\mathrm{p}}^{\mathrm{o}}$), and the error between simulated and observed times to peak $\left(\left|t_{\mathrm{p}}^{\mathrm{s}}-t_{\mathrm{p}}^{\mathrm{o}}\right|\right)$. The ratio between simulated and observed peak discharges of the TDUH-MC method ranged from 0.97 to 1.10. The average peak occurrence time error of the TDUH-MC method was 1.4 h, which was the smallest among the four methods, and the mean E_NS coefficients of the 10 flow events for validation were above 0.8. Figure 13 shows the flow hydrographs of the four routing methods for part of the flow events (event no. 20130720, 20130817, 20150709, 20160128, 20161021 and 20180916). It is demonstrated that the TDUH-MC method outperformed the remaining three routing methods.

In addition, the forecast results of six flow events in the Longhu River basin using the SUH, DUH, TDUH and the TDUH-MC method are presented in Table 7. Results of the TDUH-MC method generally showed the best performance, which also verified the TDUH-MC formula for the small watershed. In general, the TDUH-MC method was capable of better simulation in this watershed than in the Qin River basin.

Figure 14 Distributions of time-varying θ_t at different times in each sub-basin using the TDUH-MC method. (a) Flow event no. 20130817. (b) Flow event no. 20150709. (c) Flow event no.20160128. (d) Flow event no. 20180916. θ_t represents the ratio of current soil moisture storage to the corresponding maximum soil moisture capacity in the unsaturated region.

Figure 15 Time-varying velocity values of a grid cell in different storm events. (a) Time-varying velocity in storm event no. 20130817. (b) Time-varying velocity in storm event no. 20150709. The rainfall content is I_s, and the soil moisture content is θ_s.

Table 7 Comparison of four routing methods for the Longhu River basin. The values in bold mean the result performed the best among the four methods.

For flow event no. 20161021, the simulation result of the TDUH-MC method was basically consistent with that of the TDUH method. This was because the antecedent rainfall was close to saturation under this flow event. As a result, the TDUH-MC method performed the same as the TDUH method when the watershed was saturated. For flow event no. 20180916, the simulation accuracy of the TDUH-MC method was lower than that of the TDUH. The possible reason for the inaccurate flow simulation is that the antecedent rainfall was relatively small. Because the runoff generation was not dominated by the saturation excess, it was not appropriate to calculate runoff with the XAJ model.

5.4 Influence of time-varying soil moisture content on flow forecasts

In order to evaluate the influence of time-varying soil moisture content on flow forecasts, three typical flow forecast results of the TDUH-MC method were selected for comparison in the Qin River basin. Specifically, compared with the forecasting results using TDUH, results of flow event no. 20130817 using the TDUH-MC method were relatively similar, results of flow events no. 20150709 and 20160128 had a better performance, and results of flow event no. 20180916 were poor. Their corresponding temporal evolution of soil moisture content in unsaturated regions was obtained. The box-and-whisker plots of soil moisture contents of all sub-basins for flow events no. 20130817, 20150709, 20160128 and 20180916 are shown in Fig. 14. It can be seen from Fig. 14 that the soil moisture content of each sub-basin was initially low, and then the soil moisture content of the sub-basin gradually increased. Meanwhile, it was obvious that θ_t was hard to reach the maximum value. For all flow events, nine sub-basins eventually reached saturation only under the condition of flow event no. 20130817. The mean values of θ_t for flow events no. 20150709, 20160128 and 20180916 ranged from 0.5 to 0.8, and the soil moisture content did not reach the maximum during the flow events. As shown from the observed flow in Fig. 13, the peak discharge of the flow event no. 20130817 was larger than those of other flow events, reaching 3500 m³ s⁻¹, which meant that the watershed more probably reached saturation during the flow period.

As discussed in Sect. 5.3, results of flow event no. 20130817 using the TDUH-MC routing method showed the same behavior as did TDUH. This was because the simulation performance of the TDUH-MC method considering time-varying soil moisture content was the same as that of TDUH when the soil moisture content was closer to 1. Additionally, the forecast results of flow events no. 20150709 and 20160128 with the TDUH-MC routing method were obviously better than those of DUH and TDUH. The reason can be summarized as follows. The mean values of θ_t ranged from 0.5 to 0.6 for the two flow events, and the θ_t values were initially low as shown in Fig. 14. Thus, the soil moisture content had a significant impact on the shape of the hydrograph. For flow event no. 20180916, the sub-basins did not reach a global saturation, and the time-varying values of θ_t were generally high, which led to lower flow velocity than in the TDUH method. The peak occurrence times of unit hydrographs used for runoff routing calculations were in general later, leading to a lag time between maximum rainfall intensity and peak discharge for the forecast result of flow event no. 20180916.

5.5 Comparison of velocity calculated by three DUH methods

The routing method considering both time-varying rainfall intensity and soil moisture content was more accurate as discussed in Sect. 5.3. To evaluate the effect of time-varying soil moisture content on flow velocity, we selected a grid cell in sub-basin 3, in which slope and land type parameters were constant. Then, the flow velocity was calculated under different storm conditions. The storm events no. 20130817 and 20150709 were selected and compared because storm event no. 20130817 had a high intensity and long duration, and storm event no. 20150709 had a short period of heavy rainfall. Thus, soil moisture contents during the two storm events were significantly different. Figure 15 shows the time-varying velocity values of a grid cell for storm events no. 20130817 and 20150709. For the two storm events, the mean velocity of the DUH method was the largest among the three methods, followed by the TDUH method. The velocity calculated by the TDUH-MC method considering soil moisture content was the smallest. The velocity of DUH method was constant in the two storms, and that of the TDUH method varied with the change of excess rainfall. Meanwhile, the flow velocity of the TDUH-MC method was not only dominated by rainfall intensity, but was also related to soil water content.

For storm event no. 20130817, the initial soil moisture content was large, and it reached the maximum rapidly. The flow velocity of the TDUH-MC method was slightly smaller than that of the TDUH method at the initial stage of storm events. When the whole basin reached saturation, the flow velocities of the two methods became equal. Therefore, the differences between hydrographs were small when using the TDUH method and the TDUH-MC method for flow routing calculation, which led to similar forecast results.

For storm event no. 20150709, the initial soil moisture content was small, and the entire basin could not reach saturation after the rainstorm. Therefore, the grid velocity in the early stage of the storm was greatly affected by the soil moisture content. In the later stage of the rainstorm, θ_t of the watershed did not reach the maximum but was nearly close to 1. Thus, the impact of later soil moisture content on the flow velocity was small. From the above analyses, it can be concluded that the shape and duration of the unit hydrograph were mainly related to the soil moisture content at the initial stage of a storm, and when the watershed was approximately saturated, the grid flow velocity was mainly dominated by the excess rainfall.

5.6 Sensitivity analysis for the TDUH-MC method

A sensitivity analysis for the proposed formula was done in the Longhu River basin. The improved method only has two additional parameters, compared with the current model. The objective of this study was to explore the influence of the soil moisture content factor on the performance of the DUH model. Parameter γ in Eq. (3) significantly affected the significant degree of influence over how large soil moisture content will be. Thus, sensitivity analysis for parameter γ was necessary. A specific grid cell in the Longhu River basin was taken as an example, where the slope of the grid cell was set to 0.22 m m⁻¹. The coefficient of flow velocity k and the ratio of rainfall intensity to the reference rainfall intensity I_s were assumed to be 1.5 m s⁻¹ and 1, respectively. When parameter γ was 0.1, 0.5 and 1, respectively, the hillslope flow velocity values corresponding to different rainfall and soil moisture contents using the proposed formula are given in Fig. 16.

Figure 16 Time-varying flow velocity values corresponding to different parameters.

It can be seen from Fig. 16 that when θ_t was equal to 1, the proposed Eq. (3) turned to Eq. (2). The flow velocity values in the last column were the same and only changed with rainfall intensities. When I_t was equal to the reference rainfall I_c, Eq. (2) turned to Eq. (1), and the flow velocity was 0.704 m s⁻¹. After introducing a soil moisture content factor into the flow velocity formula, the flow velocity values ranged from 0.107 to 0.928 m s⁻¹ when γ was equal to 1. The flow velocity values were significantly different corresponding to different values of parameter γ. Thus, the parameter γ significantly affected the performance of the new routing method.

Moreover, the mean flow velocity of the Longhu River basin was calculated under different rainfall intensities (e.g., $\frac{I_t}{I_c}=\mathrm{0.5}$, 1, 1.5, 2, respectively). Figure 17 plots the theoretical curve of mean velocity and soil moisture content.

Figure 17 The theoretical curve of mean velocity and soil moisture content for the Longhu River basin. (a) $\frac{I_t}{I_c}=\mathrm{0.5}$. (b) $\frac{I_t}{I_c}=\mathrm{1}$. (c) $\frac{I_t}{I_c}=\mathrm{1.5}$. (d) $\frac{I_t}{I_c}=\mathrm{2}$.

Figure 17 reveals that the mean flow velocity ranged from 0.6 to 1 under different rainfall intensities without considering the influence of soil moisture content. After introducing this new factor into the current flow velocity formula, the mean flow velocity was significantly influenced by exponent γ. In addition, when the soil moisture content exceeded 0.7, the variation range of mean flow velocity decreased sharply. Results showed that the influence of parameter γ on the flow velocity decreased gradually with the increase of soil moisture content.

6 Conclusions

An improved distributed unit hydrograph method considering time-varying soil moisture content was proposed for flow routing. The TDUH-MC method comprehensively considered the changes of time-varying soil moisture content and rainfall intensity. The response of the underlying surface to the soil moisture content was considered an important factor. The Qin River basin and Longhu River basin were selected as two case studies. The SUH, DUH, TDUH and TDUH-MC routing methods were used for flow forecasting, and simulated results were compared. The sensitivity analysis was conducted for parameter γ. The main conclusions can be summarized as follows.

1. The TDUH-MC runoff routing method, considering both time-varying rainfall intensity and soil moisture content, was proposed, and the influence of the inhomogeneity of runoff generation on the routing process was considered. It was found that the soil moisture content was a significant factor affecting the accuracy of flow forecast, especially in the catchment dominated by saturation-excess runoff, and the flow velocity increased gradually with more surface runoff after considering the soil moisture content in unsaturated regions.

2. The time-varying characteristics of the DUH can be further considered by introducing both rainfall intensity and soil moisture content into the flow velocity formula, which can effectively improve the accuracy of flow forecasts. Simulation hydrographs and criteria of the two case studies showed that the accuracy of the TDUH-MC method was the highest, followed by the SUH and TDUH methods and finally the DUH method.

3. The shape and duration of the improved TDUH considering soil moisture were mainly affected by rainfall intensity. Meanwhile, soil moisture content in the initial stage of a storm also played a significant role in the characteristics of the improved TDUH. When the watershed was approximately saturated, the grid flow velocity was mainly dominated by excess rainfall.

4. Results of sensitivity analysis showed that the accuracy of the TDUH-MC method was mainly affected by soil moisture content. The influence of parameter γ on the flow velocity decreased gradually with the increase of soil moisture content.