A methodology to examine a depth-discharge constitutive equation for rainfall-runoff modelling

A number of rainfall-runoff models have been developed for hydraulic and hydrological engineering with an empha‐ sis on reproducing river discharge time series. Physicallybased rainfall-runoff models have recently reached a certain level of achievement following the advancement of com‐ puters and the development of various geographical and meteorological datasets. However, it has been pointed out that the current physically-based models do not properly reflect observed hillslope water dynamics. The present paper proposes a methodology to examine the capability of a depth-discharge constitutive equation for physicallybased rainfall-runoff modelling to simulate hillslope water dynamics. An application of the methodology suggested that 1) the targeted constitutive equation was capable of representing the depth-discharge relationship on hillslopes under the assumed conditions, 2) the runoff simulations with the constitutive equation described hillslope water flows, at least in the downward direction, and 3) there was a possibility that the parameters in the constitutive equation was determined from the internal structure of hillslope water dynamics.


INTRODUCTION
Rainfall-runoff processes are a major research topic of hydraulic and hydrological engineering and a number of rainfall-runoff models have been developed with an emphasis on reproducing river discharge time series for practical demands. Physically-based rainfall-runoff models, in particular, have been studied in the past few decades following the advancement of computers and the development of various geographical and meteorological datasets (Paniconi and Putti, 2015).
However, it has been pointed out that the current physically-based models do not properly reflect observed hillslope water dynamics. Tani (2016) raised a concern that current physically-based rainfall-runoff models such as kinematic wave models do not consider the effects of saturated-unsaturated flow and macropore flow originating from soil layer heterogeneity, although both flows are suggested to have significant influences on hillslope rainfallrunoff dynamics. Beven and Germann (2013) argued in their comprehensive review on macropores and water flow in soils that an adequate physical theory was not yet available for linking all types of flow in soils. Kirchner (2006) mentioned an implicit and unrealistic upscaling premise used in most physically-based models, and pointed out the need to develop physically-based governing equations working at the catchment or hillslope scale, which may look different from the equations for the small-scale physics.
The issues raised above give an important perspective for improving physically-based rainfall-runoff models. The validity of rainfall-runoff models is usually assessed based on the discharge hydrograph that is the final output of the models (e.g. An et al., 2010). There is, however, the possibility that a rainfall-runoff model will be essentially improved by verifying whether or not the model properly reflects hillslope water dynamics with a special focus on the model internal structure rather than the model outputs. The purposes of the present paper are to propose a methodology to examine the capability of a physically-based rainfallrunoff model to simulate hillslope water dynamics and to analyze the performance of a constitutive equation used for rainfall-runoff modelling with the proposed methodology.
An advantage of the methodology proposed is the avoidance of the model equifinality. A rainfall-runoff model can be calibrated so that its output shows a good agreement with observed discharge; however, many different parameter sets often give similar model outputs (e.g. Her et al., 2019). Given this undesirable situation, the model validity cannot be fairly judged by simple comparisons between model outputs and observations. The methodology proposed in the present study examines the model internal structure in a direct manner, not by diagnosing the model capability indirectly from the model end products. kinematic wave models generally calculate water flows using a continuity equation and a constitutive equation representing a relationship between water depth and discharge. A constitutive equation results from various factors, including hillslope topography and soil hydraulic properties, and is exactly the core of physically-based rainfall-runoff modelling. When a constitutive equation properly expresses the hillslope depth-discharge relationship, the model can be judged to be valid.
The validity of a constitutive equation should ideally be assessed with a depth-discharge relationship estimated from the observed soil moisture and flow velocity. However, substantial difficulties are expected in spatially detailed measurements that can clarify the depth-discharge relationship on hillslopes. We propose an alternative methodology to assess the validity of a constitutive equation, in which detailed simulations of hillslope water dynamics are conducted to derive probable depth-discharge relationships, and a constitutive equation is compared with the derived depth-discharge relationships to check its validity.

Hillslope water dynamics model
The proposed method utilizes simulations of hillslope water dynamics as an alternative to field observations, and accordingly, the simulations should reflect hillslope water dynamics in as much detail as possible. Water flow on a hillslope consists of subsurface flow and surface flow. Subsurface flow is further divided into saturated-unsaturated flow in soil micropores and relatively fast flow that passes in soil macropores (macropore flow). The water dynamics model to be used needs to simulate all these different flows and their interactions.
The present study constructs a model comprehensively simulating surface flow, saturated-unsaturated flow, and macropore flow at hillslopes by incorporating the effects of soil macropores into the model coupling Richards' equation for the saturated-unsaturated subsurface flow and the diffusion wave equation for the surface flow, developed by An and Yu (2014). The macropore effects are modeled with the method of Tani (2008), which gives large hydraulic conductivity to macropores. The model is hereafter called "the detailed model", and the model equations are given below.
Richards' equation (1) for the saturated-unsaturated flow and macropore flow is given as: where θ is the volumetric water content, ψ is the pressure head, z is the gravity head, K is the hydraulic conductivity, t is the time, q sub is the general source term for subsurface flow, and e sub is the exchange rate with surface flow. The diffusion wave equation (2) for the surface flow is given as: where h is the water depth, n is Manning's roughness coefficient, S is the average local bed slope (S = ∂z/ ∂x 2 + ∂z/ ∂ y 2 ), x and y are horizontal spatial coordinates, b is the ground surface elevation, q sur is the general source term for surface flow, and e sur is the exchange rate with subsurface flow. The exchange rate between surface flow and subsurface flow (e sub in Equation (1) and e sur in Equation (2)) is given by Darcy's equation with the assumption that the surface water depth is consistent with the pressure head at the ground surface for subsurface flow. Note that e sub and e sur have different dimensions ([T -1 ] and [LT -1 ]) but are essentially identical.
The present study uses Equations (3) and (4) proposed by Kosugi (1996) for θ -ψ and K -ψ relationships as in the below: where θ s is the saturated volumetric water content, θ r is the residual volumetric water content, ψ m is the pressure head corresponding to the median of the pore diameter distribution, σ is the standard deviation in the lognormal distribution of the pore diameter, and K s is the saturated hydraulic conductivity. The function Q is given by the following equation, When soil water content is relatively low, most soil water exists in micropores, and it migrates as unsaturated flow. When soil around the macropores becomes nearly saturated, water begins to enter the macropores, and then the pores start to function as a path through which soil water flows rapidly (Tani, 2016). In order to incorporate such a mechanism into the model, the hydraulic conductivity in macropores is given by the following equation proposed by Tani (2008), where K e is the hydraulic conductivity in macropores, ε is a nondimensional magnification coefficient (> 1) and ψ t is the value of pressure head where the effect of macropores begins to appear (< 0).

Derivation of depth-discharge relationship and verification of constitutive equation
The detailed model simulates hillslope water dynamics and outputs hydrological quantities such as soil water content, surface flow depth and cross-sectional average flow velocity of each computation cell. The water depth (the total amount of water existing in a cross section orthogonal to the slope flow direction) and the discharge passing through the cross section are obtained by aggregating the simulated hydrological quantities of the cross section. A depth-discharge relationship at the cross section is derived by putting together the pairs between the depth and discharge obtained at each computation time step.
A derived depth-discharge relationship should cover as wide a range as possible because it is compared for validation with the constitutive equation used in rainfall-runoff modelling. The depth and discharge for a downstream cross

EXAMINING A DEPTH-DISCHARGE CONSTITUTIVE EQUATION
section generally show a wider range of changes than that for an upstream cross section. The depth-discharge relationships should therefore be derived for cross sections located in the downstream of a hillslope while excluding the downstream end and its vicinity where the relationships may be directly affected by the given boundary conditions.
The suitability of a constitutive equation is verified by comparison with the depth-discharge relationship derived from simulated hydrological quantities. If a constitutive equation properly describes the derived depth-discharge relationship, the constitutive equation is deemed to be capable of simulating the hillslope water dynamics, at least, in the downward direction of the hillslope.

Constitutive equation to be validated
The method proposed in the present study is applied to assess the validity of the constitutive equation developed by Tachikawa et al. (2004) which integrates the saturatedunsaturated flow and surface flow (Equation (7)). The equation has been widely used for runoff simulations in many river basins and has shown high performance in reproducing river discharge time series (e.g. Sayama et al., 2012;Tanaka and Tachikawa, 2015;Tanaka et al., 2020).
where q is the discharge per unit width, h is the water depth (the total water amount of a cross section), d c is the soil micropore depth of a cross section, v c is the actual flow velocity in micropores (v c = k c i, where i is the slope gradient, and k c is the actual saturated hydraulic conductivity for micropores, defined as the micropore saturated hydraulic conductivity divided by the micropore porosity), β is a decay constant, a is the actual flow velocity in macropores (a = k a i, where k a is the actual saturated hydraulic conductivity for macropores, defined as the macropore saturated hydraulic conductivity divided by the macropore porosity), d is the total pore depth of a cross section, α = n/i where n is Manning's roughness coefficient, and m = 5/3. This equation has been developed to model runoff mechanism ranging from dry to wet soil conditions. Water flow in a hillslope soil layer is modelled as a combination of capillary water flow in micropores and non-capillary water flow in macropores. The parameter d c works as a threshold to determine the subsurface flow phase; the macropore flow appears when the water depth h becomes larger than d c . The surface flow happens in addition to the subsurface flow when h becomes larger than d. It should be noted that the continuous macropores such as a soil pipe is not explicitly supposed in the equation. The fast flow in such a configuration is expressed as the non-capillary water flow.

Supposed hillslopes and simulation conditions
In order to validate the constitutive equation (Equation (7)), water dynamics simulations were performed with the detailed model, assuming a vertically two-dimensional hillslope with the longitudinal section as shown in Figure 1 to derive the depth-discharge relationships. The slope length, the thickness of the soil layer, and the slope angle were set to be 8 m, 1 m, and 30°, respectively. The soil layer was divided into computation cells of 0.1 m length and height. Two simulation cases were conducted, i.e. without and with the soil pipe. The case with the soil pipe assumed that the pipe was continuously located directly above the bottom of the soil layer from the upper to the lower ends of the slope.
The parameter values of Equations (3) and (4) were given by referring Nakamura and Kikuzawa (2018) as follows: θ r = 0.2, θ s = 0.38, K s = 5.0 × 10 -5 m/s, ψ m = -0.2 m, and σ = 1.6. The parameters ε and ψ t in Equation (6) for the case with the soil pipe were set with reference to Tani (2016) to be 10 and -0.02 m, respectively. Manning's roughness coefficient, n, was set to be 0.185 s/m 1/3 . The initial pressure head was set to -0.1 m for all computation cells. The water depth corresponding to this pressure head value was calculated to be 0.32 m from Equation (3), which means that the soil layer was in a considerably wet condition. The bottom and the upper end of the soil layer were regarded as impermeable boundaries, and the seepage flow boundary condition was given to the hillslope lower end. A constant rainfall at 20 mm/h was given to the hillslope ground surface in the direction perpendicular to the slope from the beginning till the end of the simulation. Figure 2 shows volumetric soil water contents of the assumed hillslopes simulated by the detailed model with the above-mentioned conditions. The left and right panels present the results for the cases without and with the soil pipe, respectively. The saturated region, shown in dark red, largely evolved across the whole of the soil layer in the case without the soil pipe, while the extent of the saturated region was limited to the lower half of the soil layer in the case with the soil pipe, because the soil water was promptly discharged through the pipe. The depth-discharge relationships derived from the simu- lations are shown in Figure 3. The vertical axis represents the discharge per unit width, q [m 2 /s], and the horizontal axis represents the water depth, h [m]. The left and right panels show, similar to Figure 2, the cases without and with the soil pipe, respectively. The depth-discharge relationships were obtained at the cross section two meters upstream from the lower slope end, since it was located downstream while the water dynamics were almost unaffected by the lower end boundary condition. The dot in the panels represents a pair of depth and discharge at a certain computation time in 5-minute intervals for 10 hours from the beginning of the simulation. It was observed that the discharge increasing rate, dq/dh, became large when the depth reached around 0.36 m in both cases. This tendency was particularly remarkable in the case with the soil pipe. This is because the pressure head at the soil pipe exceeded the threshold ψ t (in Equation (6); -0.02 m) just around a depth of 0.36 m, resulting in generation of fast flow through the pipe. In the case without the soil pipe, the saturated region developed over a wide area of the soil layer, generating the large discharge, and the surface flow occurred when the depth reached around 0.38 m, accelerating the increase in the discharge rate further.

Derived depth-discharge relationship and validation of constitutive equation
Each solid line in Figure 3 represents the fitted constitutive equation so that it gives the best approximation to the depth-discharge relationship derived for each case. The fitted equation showed high agreement with the depth-discharge relationship of the case without the soil pipe. The first row of Equation (7) described the flow condition with the water depth from 0.32 m to 0.36 m, in which the micropore flow was predominant. The second row of the equation expressed the region of water from 0.36 m to 0.38 m depth. The third row of the equation captured the flow features where the depth exceeded 0.38 m and the surface flow happened as well as the subsurface flow. The differences between the depth-discharge relationship and the fitted constitutive equation became slightly large for the case with the soil pipe. Flow conditions with depth smaller than 0.36 m were almost the same as those without the soil pipe, and the first row of Equation (7) gave a fairly good approximation. The second row of the equation expressed the rapid increase in discharge caused by the fast flow through the soil pipe when the depth exceeded 0.36 m. Table I summarizes the parameter values obtained by fitting Equation (7) to the depth-discharge relationships. The actual saturated hydraulic conductivity at macropores, k a , predominantly controlling the non-capillary water flow, showed a large difference between the cases without and with the soil pore. The value of k a for the case with the soil pipe was more than three times larger than that without the soil pipe, reflecting the combined effects of relatively large pores and the soil pipe.  The results of the detailed model and the constitutive equation were qualitatively in good agreement, especially in terms of the observation that the discharge gradually increased during the first two hours and then rapidly increased to reach the maximum rate. In the case without the soil pipe (left panel), the surface flow pattern was reproduced well by the constitutive equation although there were some differences in its occurrence time and the maximum rate.

DISCUSSION
The results described above suggest the following points: 1) the constitutive equation (Equation (7)) can represent the depth-discharge relationship on hillslopes under the assumed conditions, 2) the runoff simulations with the constitutive equation can describe the water flows on hillslopes, at least in the downward direction, and 3) there is a possibility that the parameters in the constitutive equation can be determined not by the calibration based on the discharge but by being more directly associated with the hillslope water dynamics.
A trigger of the present study was the concern that the physically-based rainfall-runoff models may not reflect hillslope water dynamics properly; it motivated us to examine a constitutive equation itself rather than the model outputs. Given the first two points above suggested by the results, it is deduced that the targeted constitutive equation describes, if not completely, hillslope water flows and their transitions depending upon soil moisture conditions using a piecewise form.
It should be noted that the present results are valid under limited conditions; they do not necessarily present the universal applicability of the rainfall-runoff models based on the constitutive equations. For example, most constitutive equations assume single-value depth-discharge relationships, but in reality, multiple different discharges are possible for a single value of water depth, depending on hillslope topography, soil properties, soil water profiles, and rainfall conditions. Even in such cases, it may be possible to match the simulated hydrographs with the observed ones by calibrating (tuning) parameter values of rainfallrunoff models. However, this impairs the benefits of physically-based rainfall-runoff modelling. The methodology proposed in the present study can help unbiased assessment of how good or bad physically-based rainfall-runoff models are, and its essential goal is to contribute to model improvements by clarifying problems existing in the model structures.

CONCLUSIONS
We proposed a methodology to examine the capability of a physically-based rainfall-runoff model to simulate hillslope water dynamics, in which the validity of a constitutive equation for rainfall-runoff modelling was assessed using the depth-discharge relationship derived from the detailed simulations of hillslope water dynamics. The proposed methodology is fundamentally different from a conventional verification based on the simulated hydrographs because the methodology targets the internal structure of the rainfall-runoff model.
The proposed methodology was applied to assess the validity of the constitutive equation which integrated the saturated-unsaturated flow and surface flow. The results suggested that 1) the constitutive equation was capable of representing the depth-discharge relationship on hillslopes under the assumed conditions, 2) the runoff simulations with the constitutive equation described hillslope water flows, at least in the downward direction, and 3) there was a possibility that the parameters in the constitutive equation can be determined from the internal structure of hillslope water dynamics. The proposed methodology can clarify the conditions under which the physically-based rainfall-runoff models work properly and should contribute to improving the models by exploring the problems existing in the model structures.