Preferential flow is significant for its contribution to rapid response to hydrologic inputs at the soil surface and unsaturated zone flow, which is critical for flow generation in rainfall–runoff (RR) models. In combination with the diffuse and source-responsive flow equations, a new model for water infiltration that incorporates preferential flow is proposed in this paper. Its performance in estimating soil moisture at the catchment scale was tested with observed water content data from the Elder sub-basin of the South Fork Eel River, located in northern California, USA. The case study shows that the new model can improve the accuracy of soil water content simulation even at the catchment scale. The impacts of preferential flow on RR simulation were tested by the Modello Idrologico Semi-Distributio in continuo lumped hydrological model for the Elder River basin. Eleven significant floods events, which were defined as having flood peak magnitudes greater than ten times average discharge during the study period, were employed to assess runoff simulation improvement. The accuracy of the runoff simulation incorporating the preferential flow at the catchment scale improved significantly according to the likelihood ratio test.

## INTRODUCTION

Water infiltration plays a major role in rainfall–runoff (RR) generation. Thus, reliable infiltration equations are essential for RR modeling accuracy in both theoretical and practical hydrological applications (Swamee *et al.* 2014). Preferential flow, which refers to the phenomenon of faster than average water movement through only a fraction of the pore space in a soil column, occurs when water infiltrates through the vadose zone of soil. Preferential flow can influence the soil moisture distribution which has effects on runoff hydrographs, especially subsurface flow in highly heterogeneous slopes (Hendrickx & Flury 2001; Jarvis 2007). It is also an important component of flow generation, especially in many humid regions where overland flow is rarely observed.

Preferential flow occurs in many soils, and is related to worm holes, root channels, and inter-aggregate fissures (Bouma 1981; Beven & Germann 1982), and associated with textural differences (Cislerova *et al.* 2002; Snehota *et al.* 2008). Preferential flow has also been increasingly recognized as a process of great practical significance for the transport of water at different scales (Ogawa *et al.* 2002; Krzeminska *et al.* 2012). There has been a wide variety of approaches to characterizing preferential flow, including the analysis of the flow in macropores formed by shrinkage and biological processes, within or over structured textural heterogeneous porous formations, unstable flows due to gravity instabilities or viscous instabilities, flow in natural soil pipes or man-made drain tubes, and topographically controlled local soil and groundwater recharge (Gerke *et al.* 2010). Poiseuille's equation (Ahuja & Hebson 1992), the Green and Ampt or Philip infiltration models (Ahuja & Hebson 1992; Chen & Wagenet 1992), the kinematic wave equation (Germann & Beven 1985; Jarvis 1994), and the Richards equation (Gerke & van Genuchten 1993) are used to calculate water flow in macropores or inter-aggregate pores. These equations assume that the porous medium consists of two interacting regions, and are widely applied in hydrologic simulation (Šimůnek *et al.* 2003; Gerke 2006; Jarvis 2007; Köhne *et al.* 2009; Vogel *et al.* 2010).

These dual-permeability models neglect or simplify the processes and geometrics of preferential flow, which are significantly different from those of non-preferential flow (Nimmo 2010). For example, capillarity is a dominant influence in diffuse unsaturated flow. ‘Source-responsive’ model refers to the thesis that preferential flow may be triggered and modulated in response to the source of water to the unsaturated medium, with the model's dependence on local potential gradients much less than in traditional unsaturated flow models (Nimmo 2010). However, upscaling the preferential flow process in the pores at the micrometer-range scale to the profile scale or even hillslope scale is still limited. There is, as yet, no real consensus as to how to quantify the effects of preferential flow pathways at catchment scales (Beven 2010; Gerke *et al.* 2010). The widely used infiltration equations in hydrological models cannot effectively simulate preferential flow resulting in rapid infiltration (Nieber & Sidle 2010; Beven & Germann 2013; Shao *et al.* 2015).

The objectives of this study were to (1) upscale the preferential flow process to the catchment scale and (2) analyze the effect of preferential flow on simulating the runoff process. To achieve this, the source-responsive model was combined with the Green and Ampt model, leading to expressions for soil moisture simulation at a catchment scale. We also applied a lumped hydrological model (Modello Idrologico Semi-Distribuito in continuo (MISDc)) to include the preferential flow process into flow generation, with the expectation of improving predictive models compared to no accounting of preferential flow.

## METHODOLOGY

### Model for the flux of diffuse and preferential flow using the Green–Ampt model

*q*(

*z*,

*t*) is the total flux density (or the volumetric flow rate in the

*z*direction per unit cross-sectional area of medium, m s

^{−1});

*z*is the depth;

*t*is time;

*q*(

_{D}*z*,

*t*) is the diffuse-flow portion of flux density;

*q*(

_{S}*z*,

*t*) is the source-responsive flux density;

*K*(

*θ*) is the hydraulic conductivity;

_{D}*θ*is the diffuse-flow domain moisture content;

_{D}*ψ*is the suction (capillary) at wetting front relevant to diffuse flow;

*r*(

*z,t*) is a dimensionless factor referred to as the active area fraction in reference to the concept of film flow along some of the internal surface area of pores at given

*t*and

*z*; and

*q*

_{Smax}(

*z*) is the maximum volumetric flux density sustainable within the source-responsive domain at

*z*.

*r*(

*s,t*) ranges from 0 to 1 and depends on some yet-unspecified conditions of the flow system.

*q*

_{Smax}(

*z*) is a constitutive property of the medium and varies with

*z*as the size or number of preferential flow paths varies with

*z*.

*θ*=

_{D}*rML*. In the most general case, the volume of the source-responsive domain filled with actively flowing free-surface film and was assumed to be the relations between flow speed and flow thickness. The average vertical flow speed is , and the maximum source-responsive flow corresponds to all macropore facial areas having a maximum supportable free surface film thickness,

*L*

_{max}(

*z*), therefore, the source-responsive formula for flux density in films is (Nimmo 2007, 2010). Equation (1) can be transformed into: where

*g*is gravitational acceleration (9.8 m s

^{−2}); υ is the kinematic viscosity of water (υ = 1.0 × 10

^{−6}m

^{2}s

^{−1}at 20 °C);

*L*

_{max}is the uniform and maximum supportable free-surface film thickness (m);

*L*

_{max}corresponds to the situation of maximum source-responsive flow and its value is a characteristic of the general nature of the flow system; and

*M*(

*z*) is the facial area density (m

^{−1}).

Green & Ampt (1911) produced an elegant infiltration formula for predicting cumulative infiltration in dry or uniformly wet soil based on the Darcy–Buckingham law for fluxes. In the Green–Ampt model, if the following assumptions for the diffuse-flow domain are satisfied, the approximate infiltration model for fluxes of diffuse and preferential flow can be derived:

As rain continues to fall and water infiltrates, the wetting front advances at the same rate as depth, which produces a well-defined wetting front.

The volumetric water content remains constant above and below the wetting front as it advances.

The soil-water suction immediately below the wetting front remains constant with both time and location as the wetting front advances.

*F*=

*F*

_{1}+

*F*

_{2}, where

*F*

_{1}and

*F*

_{2}are the cumulative depth of infiltration diffuse-flow domain and source-responsive domain, respectively.

*F*

_{1}can be determined by

*F*

_{1}=

*z*(

_{f}*θ*–

_{s}*θ*), where

_{i}*z*is the thickness of the soil layer.

_{f}*F*

_{2}is unknown and must be derived. The infiltration rate

*f*is the first derivative of

*F*with respect to time

*t*. According to the Green–Ampt model for the diffuse-flow domain,

*F*

_{1}can be calculated from For determining

*F*

_{2}, the active-area fraction

*ra*(

*z,t*) in Equation (3) can be assumed to be a declining function of time, nearly uniform with depth

*z*. Therefore, a reasonable representation of

*ra*(

*z,t*) is which equals unity (indicative of maximal preferential flow) at the starting time

*t*

_{0}of the ponded interval, and then declines exponentially with a time constant

*τ*

_{1}, which will probably have a value of a few hours (Nimmo 2010).

*F*

_{2}can be calculated by Therefore, the total (cumulative) infiltration can be determined from Equations (6) and (9).

### A modified MISDc lump hydrological model which accounts for the diffuse and preferential flow in soil water balance

The MISDc is a semi-distributed RR model developed by Brocca *et al.* (2011) at the Research Institute for Geo-Hydrological Protection of Perugia, Italy. The model, which has a simple structure with few parameters and low computational effort, couples a soil water balance (SWB) model with an event-based RR model to obtain a simple and robust structure (Anctil *et al.* 2004; Manfreda *et al.* 2005; Sheikh *et al.* 2009). MISDc consists of two main components as follows.

#### The SWB model to simulate the soil moisture temporal pattern

*f*(

*t*) is estimated by Equation (3),

*e*(

*t*) is the evapotranspiration rate,

*gg*(

*t*) is the drainage rate due to the interflow and/or deep percolation, and

*W*

_{max}is the maximum water capacity of the soil layer.

*W*

_{max}can be estimated with:

*gg*(

*t*) is determined by the following relation (Famiglietti & Wood 1994) where

*λ*is a pore size distribution index linked to the structure of the soil layer.

*e*(

*t*) is represented as a linear relationship, between the potential evapotranspiration

*ET*(

_{p}*t*) and the soil saturation (Brocca

*et al.*2011):

*ET*(

_{p}*t*) is computed with the empirical relation of Blaney and Criddle as modified by Doorenbos & Pruitt (1977): where

*T*(

_{a}*t*) is the air temperature in °C,

*ξ*is the percentage of total daytime hours during the period considered out of total daytime hours of the year, and

*b*is a parameter to be calibrated.

The output of the SWB model is the degree of saturation, *W*(*t*)/*W*_{max}, that is used to determine the initial condition in the MISD model.

#### A semi-distributed event-based RR model for flood simulation

A semi-distributed event-based RR model (MISD) developed by Corradini *et al.* (1995) employed the Soil Conservation Service-Curve Number method for abstraction (SCS-CN) for estimation of losses; the geomorphological instantaneous unit hydrograph (IUH) for routing of rainfall in excess of subcatchments (Gupta *et al.* 1980); and the linear reservoir IUH, a diffusive linear approach, for routing of areas draining directly into the main channel.

*N*elements that represent the subcatchment with outlets along the main channel and assumed homogeneous to constitute a lumped system, the rainfall excess

_{b}*ɛ*(

_{j}*t*) for element

*j*(

*j*= 1,…,

*N*) is estimated by SCS-CN formula (Chow

_{b}*et al.*1988): where

*R*(

_{j}*t*) is the rainfall depth from the beginning of the storm,

*S*is the soil potential maximum retention at the beginning of the storm, and λ

_{j}_{1}is a parameter linked to the initial abstraction, assumed constant for all the elements. The direct runoff hydrograph,

*Y*(

*t*), at the element outlet is given by the convolution of

*ɛ*(

*t*) and the IUH,

*h*(

*t*), as: where

*τ*

_{2}is an auxiliary variable for time and

*A*is the element area. The dynamic parameter lag time,

*L*, is given by the empirical lag–area relationship proposed by Corradini

*et al.*(2002): where

*η*is a parameter to be calibrated.

*Q*(

*t*), at the catchment outlet is estimated by a diffusive linear approach (Troutman & Karlinger 1985): with

*p*(

*t*) the diffusive routing function given by where

*L*is the distance between the element and catchment outlet,

_{c}*c*is the celerity, and

*D*is the diffusivity parameter.

## STUDY AREA AND DATA

^{2}basin ranges in elevation from 410 to 1,284 m above sea level and drains a landscape of steep hill slopes and narrow canyons. The basin has a Mediterranean climate with annual average precipitation of approximately 2,150 mm y

^{−1}, the vast majority of which falls as rain, mostly concentrated during winter (October–April). The average annual air temperature was historically 10 °C, and the seasonal minimum and maximum temperatures were −8 °C and 32 °C, respectively (Kim

*et al.*2014). The basin is located in the Sierran Steppe-Mixed Forest-Coniferous Forest-Alpine Meadow eco-region, and vegetation grows in a mosaic of mixed forest, oak woodland, mixed chaparral, riparian, and grassland communities (http://angelo.berkeley.edu/). Soils in the basin are classified as Inceptisols and Ultisols and primarily belong to the Hugo and Josephine soil series (Mast & Clow 2000). The length of its main stem is about 8 km, with an average stream gradient of 80 m/km. Stream flows are highly variable during heavy storms due to excess runoff from steep hill slopes, but are well sustained by springs during low rainfall periods. Mean monthly discharge ranges from 0.03 m

^{3}s

^{−1}in September to 1.9 m

^{3}s

^{−1}in January, with 1,290 mm as the average annual runoff (Nourani & Zanardo 2014).

The stream hosts a US Geological Survey hydrologic benchmark station near its mouth (latitude 39 ° 43′ 47″, longitude 123 ° 38′ 34″) with 15-min continuous discharge records from 1967 to 2013. Rainfall data (from six stations), air temperature data (from 12 gauge stations), and soil moisture data (from two stations) are from Berkeley Sensor Database: Angelo Reserve Data. We used data from October 28, 2012 to August 31, 2014, converted to half hour (30 min) intervals. The areal data were calculated from every gauge station by using Thiessen polygon method.

For assessing the impacts of infiltration and preferential flow on RR, flood events characterized by a continuous rainfall pattern with greater than ten times average discharge () were selected; there were 11 such flood events in the study period. A time-domain reflectometry device was used to measure the soil moisture.

## RESULTS AND DISCUSSION

### Results from the MISDc model whose infiltration is estimated by the Green–Ampt equation

*W*

_{max}(the maximum water capacity of the soil layer) is the most sensitive parameter of the model and was fixed in accordance with the spatial distribution of the curve values computed from soil/land use characteristics (Chow

*et al.*1988; Brocca

*et al.*2011). The temporal evolution of rainfall in the study period is shown in Figure 2. Figure 3 illustrates the temporal variation of the observed and simulated soil moisture without accounting for the diffuse and preferential flow, while Figure 4 compares the observed and simulated discharge (along with the starting time of the 11 selected flood events) at the outlet gauge station. The variation of soil moisture is consistent with the change of the rainfall. However, most of the soil moisture simulation from the model during the heavy rainfall periods is more than the observed data while the simulation results during the light rainfall periods seem to be smaller than the observed value.

The Green–Ampt equation model | The model combing the source-responsive model and the Green–Ampt model | ||||||||
---|---|---|---|---|---|---|---|---|---|

ψ_{f} | Ks (mm h^{−1}) | λ | b | ψ_{f} | Ks (mm h^{−1}) | λ | b | τ | (m^{−1}) |

−0.800 | 12.505 | 0.640 | 0.761 | −0.800 | 12.505 | 0.520 | 0.827 | 550 | 40,500 |

The Green–Ampt equation model | The model combing the source-responsive model and the Green–Ampt model | ||||||||
---|---|---|---|---|---|---|---|---|---|

ψ_{f} | Ks (mm h^{−1}) | λ | b | ψ_{f} | Ks (mm h^{−1}) | λ | b | τ | (m^{−1}) |

−0.800 | 12.505 | 0.640 | 0.761 | −0.800 | 12.505 | 0.520 | 0.827 | 550 | 40,500 |

is the average facial area density.

The accuracy of the soil moisture and runoff simulated models is demonstrated with standard performance indicators. Nash–Sutcliffe efficiency (NSE) coefficients were found equal to 0.61, while the root mean squared error (RMSE) was less than 0.04 when the Green–Ampt equation was used to simulate the soil moisture. As can be seen in Figures 2 and 3, the variation of the simulated soil moisture reflects the response of the rainfall to soil. The relatively large differences between the observed and the simulated soil moisture in Figure 3 are probably related to the poor soil moisture data representation for the study catchment. Specifically, the water content from the Green–Ampt equation for high rainfall intensity is overestimating the infiltration. The poor soil moisture representation is because there are only two stations. Additionally, the soil characteristics might be quite variable at different depths and at the whole catchment scale. NSE is 0.71 for the selected flood events. This implies that the MISDc hydrological model can mediate the errors from the soil moisture model, especially for the peak flow simulation.

### Results from the MISDc model whose infiltration simulation has accounted for the diffuse and preferential flow

*λ*linked to the structure of the soil layer and the correction coefficients for potential evapotranspiration

*b*are different while the parameters related to the Green–Ampt infiltration equation are equal. As expected from Equations (9), (11), and (14), both the simulations of infiltration are based on the Green–Ampt infiltration equation while the simulation of the preferential flow is independent from the non-preferential process.

As shown in Figure 5, the performance of the soil moisture balance model incorporating preferential flow is better, with RMSE decreasing from 0.04 to 0.037 and NSE increasing from 0.61 to 0.66. Figure 6 compares the observed and simulated runoff with preferential flow in the 11 selected flood events, and shows a consistent temporal variation of the soil moisture balance simulation. Including preferential flow increased the NSE value to 0.73. We note that the preferential flow rate will decrease with the time *t* according to Equation (8). Therefore, the effect of preferential flow on soil moisture simulation should be reduced during a rainfall event.

*F*-test (NIST/SEMATECH e-Handbook of Statistical Methods, 2003), a likelihood ratio test, was employed in this study. The

*F*-test is calculated with: where RMSE

*is the root mean square error of each model*

_{i}*i*;

*N*is the number of the data set;

*k*

_{1}is the number of the MISDc model parameters; and

*k*

_{2}is the number of the MISDc model parameters incorporating preferential flow. Application of the

*F*-test shows that use of more parameters (that is, incorporating preferential flow) in the SWB model and the runoff simulation can statistically improve accuracy. This implies that the use of preferential flow at the catchment scale can improve simulation of soil moisture and discharge.

## CONCLUSIONS

This study demonstrated that soil moisture simulation accuracy in SWB models can be improved by incorporating preferential flow in the infiltration process.

The sensitivity of soil moisture to preferential flow is reduced over time during a rainfall event. However, the runoff simulation during flood events are highly correlated with the soil moisture evaluated by the SWB model, with low sensitivity to preferential flow. Incorporating preferential flow also significantly improves the accuracy of runoff simulation at the catchment scale through the likelihood ratio test (according to Occam's Razor rule, the more parameters, the higher accuracy of the simulation).

Moreover, the spatial heterogeneity of the soil characteristics and the real areal observed soil moisture data at the catchment scale can influence the accuracy of the model simulation and should be taken into account in future studies.

## ACKNOWLEDGEMENTS

The authors thank William E. Dietrich and his research team members for contributing both data and valuable insights to this study. This work was supported by the National Science Foundation CZP EAR-1331940 for the Eel River Critical Zone Observatory and National Natural Science Foundation of China (Grant No. 51379148 and 51579183). Dedi Liu was also supported by the China Scholarship Council (Grant No. 201308420310).