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.

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.

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

Nimmo (2010) developed an approach (Equation (1)) to model the flux of diffuse and preferential flow by combining the Darcy–Buckingham law and the continuity equation:
1
where 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; qD(z,t) is the diffuse-flow portion of flux density; qS(z,t) is the source-responsive flux density; K(θD) is the hydraulic conductivity; θD is the diffuse-flow domain moisture content; ψ 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 qSmax(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. qSmax(z) is a constitutive property of the medium and varies with z as the size or number of preferential flow paths varies with z.
From the geometry of flow domains the sources-responsive flow, the sources-responsive water content is θ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, Lmax(z), therefore, the source-responsive formula for flux density in films is (Nimmo 2007, 2010). Equation (1) can be transformed into:
2
where g is gravitational acceleration (9.8 m s−2); υ is the kinematic viscosity of water (υ = 1.0 × 10−6 m2 s−1 at 20 °C); Lmax is the uniform and maximum supportable free-surface film thickness (m); Lmax 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:

  1. 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.

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

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

There are two forms of water infiltrating to the soil. One is the diffuse-flow domain and the other is source-responsive domain. The total cumulative depth of infiltration is F = F1 + F2, where F1 and F2 are the cumulative depth of infiltration diffuse-flow domain and source-responsive domain, respectively. F1 can be determined by F1 = zf(θsθi), where zf is the thickness of the soil layer. F2 is unknown and must be derived. The infiltration rate f is the first derivative of F with respect to time t.
3
4
5
According to the Green–Ampt model for the diffuse-flow domain, F1 can be calculated from
6
For determining F2, 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
7
which equals unity (indicative of maximal preferential flow) at the starting time t0 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). F2 can be calculated by
8
9
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

The surface soil layer is assumed as a spatially lumped system for the following equations:
10
where 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 Wmax is the maximum water capacity of the soil layer. Wmax can be estimated with:
11
The drainage rate gg(t) is determined by the following relation (Famiglietti & Wood 1994)
12
where λ is a pore size distribution index linked to the structure of the soil layer.
The evapotranspiration component e(t) is represented as a linear relationship, between the potential evapotranspiration ETp(t) and the soil saturation (Brocca et al. 2011):
13
ETp(t) is computed with the empirical relation of Blaney and Criddle as modified by Doorenbos & Pruitt (1977):
14
where Ta(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)/Wmax, 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.

If a given catchment is divided into Nb elements that represent the subcatchment with outlets along the main channel and assumed homogeneous to constitute a lumped system, the rainfall excess ɛj(t) for element j (j = 1,…, Nb) is estimated by SCS-CN formula (Chow et al. 1988):
15
where Rj(t) is the rainfall depth from the beginning of the storm, Sj is the soil potential maximum retention at the beginning of the storm, and λ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:
16
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):
17
where η is a parameter to be calibrated.
The direct runoff hydrograph, Q(t), at the catchment outlet is estimated by a diffusive linear approach (Troutman & Karlinger 1985):
18
with p(t) the diffusive routing function given by
19
where Lc is the distance between the element and catchment outlet, c is the celerity, and D is the diffusivity parameter.
We tested the new infiltration and preferential flow model by applying it to the Elder Creek basin, located in the South Fork Eel River basin in northern California, USA (Figure 1). The watershed is within the Angelo Coast Range Reserve, which since 1959 has been protected by the University of California for research purposes. This 16.84 km2 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 m3 s−1 in September to 1.9 m3 s−1 in January, with 1,290 mm as the average annual runoff (Nourani & Zanardo 2014).
Figure 1

The Elder Creek River basin.

Figure 1

The Elder Creek River basin.

Close modal

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 from the MISDc model whose infiltration is estimated by the Green–Ampt equation

To analyze infiltration, the Green–Ampt equation from the MISDc hydrological model was used. The observed soil moisture data were used to calibrate the parameters of the Green–Ampt equation (Table 1). The parameter Wmax (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.
Table 1

Parameter values of the SWB model based on different formulations for infiltration

The Green–Ampt equation model
The model combing the source-responsive model and the Green–Ampt model
ψfKs (mm h−1)λbψfKs (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
ψfKs (mm h−1)λbψfKs (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.

Figure 2

The temporal evolution of rainfall in Elder Creek River basin.

Figure 2

The temporal evolution of rainfall in Elder Creek River basin.

Close modal
Figure 3

Comparison between observed and simulated soil moisture without accounting for diffuse and preferential flow.

Figure 3

Comparison between observed and simulated soil moisture without accounting for diffuse and preferential flow.

Close modal
Figure 4

Comparison between observed and simulated runoff without accounting for diffuse and preferential flow, less base flow (6.27 m3/s).

Figure 4

Comparison between observed and simulated runoff without accounting for diffuse and preferential flow, less base flow (6.27 m3/s).

Close modal

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

With the diffuse and preferential flow taken into account in the soil moisture simulation, observed and simulated soil content and discharge/runoff for the selected flood events are as shown in Figures 5 and 6, respectively. For both formulations and their simulation results in SWB models with (Figure 2) or without (Figure 5) accounting for the preferential flow, their parameters were estimated by minimizing the mean square error and are presented in Table 1. The pore size distribution index λ 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.
Figure 5

Comparison between observed and simulated soil moisture with accounting for diffuse and preferential flow.

Figure 5

Comparison between observed and simulated soil moisture with accounting for diffuse and preferential flow.

Close modal
Figure 6

Comparison between observed and simulated runoff with accounting for diffuse and preferential flow, less base flow (6.27 m3/s).

Figure 6

Comparison between observed and simulated runoff with accounting for diffuse and preferential flow, less base flow (6.27 m3/s).

Close modal

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.

To demonstrate that the goodness of fit is improved more than would be expected by simply adding more model parameters in Equations (3) or (11), an 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:
20
where RMSEi is the root mean square error of each model i; N is the number of the data set; k1 is the number of the MISDc model parameters; and k2 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.

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.

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).

Ahuja
L. R.
Hebson
C.
1992
Root Zone Water Quality Model
.
GPSR Technical Report No. 2
,
USDA, ARS, Fort Collins
,
CO
,
USA
.
Anctil
F.
Michel
C.
Perrin
C.
Andreassian
V.
2004
A soil moisture index as an auxiliary ANN input for stream flow forecasting
.
J. Hydrol.
286
,
155
167
.
Beven
K. J.
Germann
P.
1982
Macropores and water flow in soils
.
Water Resour. Res.
18
,
1311
1325
.
Beven
K.
Germann
P.
2013
Macropores and water flow in soils revisited
.
Water Resour. Res.
49
(
6
),
3071
3092
.
Chow
V. T.
Maidment
D. R.
Mays
L. W.
1988
Applied Hydrology
.
McGraw-Hill
,
New York
,
USA
.
Cislerova
M.
Vogel
T.
Votrubova
J.
Robovska
A.
2002
Searching below thresholds: tracing the origins of preferential flow within undisturbed soil samples
. In:
Environmental Mechanics: Water, Mass and Energy Transfer in the Biosphere
(
Raats
P. A. C.
Smiles
D.
Warrick
A. W.
, eds).
Geophysical Monograph Series 129. American Geophysical Union
,
Washington, DC
,
USA
, pp.
265
274
.
Corradini
C.
Melone
F.
Ubertini
L.
1995
A semi-distributed model for direct runoff estimate
. In:
Applied Simulation and Modelling
(
Hamza
M. H.
, ed.).
IASTED ACTA Press
,
Anaheim, CA
,
USA
, pp.
541
545
.
Corradini
C.
Morbidelli
R.
Saltalippi
C.
Melone
F.
2002
An adaptive model for flood forecasting on medium size basins
. In:
Applied Simulation and Modelling
(
Ubertini
L.
, ed.).
IASTED Acta Press
,
Anaheim, CA
,
USA
, pp.
555
559
.
Doorenbos
J.
Pruitt
W. O.
1977
Background and development of methods to predict reference crop evapotranspiration (ETo). In: FAO-ID-24, Appendix II
, pp.
108
119
.
Famiglietti
J. S.
Wood
E. F.
1994
Multiscale modeling of spatially variable water and energy balance processes
.
Water Resour. Res.
11
,
3061
3078
.
Gerke
H. H.
2006
Preferential flow descriptions for structured soils
.
J. Plant Nutr. Soil Sci.
169
,
382
400
.
Gerke
H. H.
Germann
P.
Nieber
J.
2010
Preferential and unstable flow: from the pore to the catchment scale
.
Vadose Zone J.
9
,
207
212
.
Gupta
V. K.
Waymire
E.
Wang
C. T.
1980
A representation of an instantaneous unit hydrograph from geomorphology
.
Water Resour. Res.
16
,
855
862
.
Hendrickx
J. M.
Flury
M.
2001
Uniform and Preferential Flow Mechanisms in the Vadose Zone
.
National Academy Press
,
Washington, DC, USA, pp
.
149
188
.
Jarvis
N. J.
1994
The MACRO Model (Version 3.1). Technical Description and Sample Simulations
.
Reports and Dissertations 19
.
Department of Soil Science, Swedish University of Agricultural Science
,
Uppsala
,
Sweden
, p.
51
.
Krzeminska
D. M.
Bogaard
T. A.
van Asch
W. J.
van Beek
P. H.
2012
A conceptual model of the hydrological influence of fissures on landslide activity
.
Hydrol. Earth Syst. Sci.
16
,
1561
1576
.
Mast
M. A.
Clow
W. D.
2000
Environmental Characteristics and Water Quality of Hydrologic Benchmark Network stations in the Western United States. US Geological Survey Circular 1173-D
, p.
115
.
NIST/SEMATECH e-Handbook of Statistical Methods
2003
,
updated 2010. http://www.itl.nist.gov/div898/handbook/eda/section3/eda3673.htm (accessed 10 September 2010)
.
Ogawa
S.
Baveye
P.
Parlange
J. Y.
Steenhuis
T.
2002
Preferential flow in the field soils
.
Forma
17
,
31
53
.
Shao
W.
Bogaard
T. A.
Bakker
M.
Greco
R.
2015
Quantification of the influence of preferential flow on slope stability using a numerical modelling approach
.
Hydrol. Earth Syst. Sci.
19
,
2197
2212
. doi:10.5194/hess-19-2197-2015.
Snehota
M.
Sobotkova
M.
Cislerova
M.
2008
Impact of the entrapped air on water flow and solute transport in heterogeneous soil: Experimental set-up
.
J. Hydrol. Hydromech.
56
,
247
256
.
Swamee
P. K.
Rathie
P. N.
Ozelim
L. C. S. M.
Cavalcante
A. L. B.
2014
Recent advances on solving the three-parameter infiltration equation
.
J. Hydrol.
509
,
188
192
.
Vogel
T.
Brezina
J.
Dohnal
M.
Jaromir
D.
2010
Physical and numerical coupling in dual-continuum modeling of preferential flow
.
Vadose Zone J.
9
,
260
267
.