Contribution of PoreScale Approach to Macroscale Geofluids Modelling in Porous Media
View this Special IssueResearch Article  Open Access
Quantify the Pore Water Velocity Distribution by a Celerity Function
Abstract
Celerity and velocity of subsurface flow in soil porous medium are intimately linked with tracer transport, while only few current studies focused on their relations. This study conducted theoretical analyses based on a pore bundle model, under which celerity in unsaturated flow is equivalent to maximum velocity. Furthermore, under 3 models BrooksCorey model and unmodified and modified Mualemvan Genuchten models, we used a celerity function to derive breakthrough curves aiming to quantify the advective tracer transport for 5 typical soil textures. The results showed that the celerity under nearsaturated conditions can be 5 up to 100 times larger than the saturated hydraulic conductivity, and a small volumetric fraction (<15%) of pores contributed more than half of the specific discharge. First arrival time and extensive tailing of the breakthrough curves were controlled by maximum velocity and velocity distribution, accordingly. As the kinematic ratio in the BrooksCorey model remained constantly for each specific soil, we used it to quantify the ratio of maximum tracer velocity over average tracer velocity. We also found that a bimodal soil hydraulic function (for a dualpermeability model) may result in similar soil hydraulic conductivity functions for different parameter sets, but their celerities are different. The results showed that the proposed celerity function can assist in investigating subsurface flow and tracer transport, and the kinematic ratio could be used to predict the first arrival time of a conservative tracer.
1. Introduction
Subsurface flow is conducted in the soil porous medium, in which the pore structure (i.e., pore connectivity, tortuosity, and pore size distribution) controls the permeabilities and pore water velocities in numerous flow paths [1–3]. In hydrological studies, it is impractical to deterministically represent the porescale fluid dynamics and transport processes in soil, for the variability of pore water velocity is difficult to observe or predict [4–6]. For describing the subsurface flow processes, the continuum approach was often used to conceptualize the discrete materials and water particles as one single continuous mass. The soil hydraulic properties (e.g., the lumped hydraulic conductivity and moisture capacity) and state variables (e.g., volumetric water content and capillary pressure) are often quantified in a representative elementary volume (REV) [7, 8]. The continuum subsurface flow equation (i.e., the wellknown DarcyRichards equation) calculates specific discharge and average pore water velocity for transport phenomenon. Specifically, the specific discharge quantifies the volumetric flow [9, 10], and the average pore water velocity can be coupled with an advectiondiffusion equation to simulate solute/thermal transport [10–12]. Applications of the DarcyRichards equation relay on the indicators of velocities and specific discharge.
For describing the water fluid dynamics, the terminologies of velocity (i.e., pore water velocity) and celerity have been commonly used in various hydrological systems, such as river channels [13–15], estuaries [16, 17], and soil porous medium [8, 18, 19]. Theoretically, velocity expresses mass flux, while celerity represents the speed of pressure propagation through a flow domain [20]. In surface water, for instance, the pressure wave is a flood wave propagating through river channels. Similarly, the celerity in subsurface flow describes a perturbationinduced pressure wave that is affected by precipitation, evaporation, or fluid injection and extraction [21, 22]. In saturate soil, the difference between velocity and celerity can be illustrated by a virtual experimental [20], and the authors pointed out: “in a cylinder full of sand and saturated with water, changing the flow rate or head at the input boundary will immediately cause a change in flow at the output boundary. While the water flow velocity through the sand is slow, the celerity in this case is (theoretically) instant, hence the immediate response. At larger scales, this case is analogous to a confined aquifer with incompressible water and rock. Allowing for the compressibility will slow the celerity a little, but the velocities of flow will still be much less than the celerities.”
In confined saturated flow, the nearlyinstant celerity is caused by the low compressibility of water and soil porous medium. While in unsaturated soil, the pressure propagation was achieved by the variation of capillary pressure that is driven by different mechanisms. For example, a response of the pore water pressure in unsaturated soil can be caused by either preferential flow [23] or pressure wave through the entrapped air [24]. However, how to distinguish the initial causeeffect mechanism of pressure response is still a challenge [22, 25]. Many studies used the kinematic wave equation to simulate unsaturated flow [8, 18, 26, 27], in which the celerity was defined as a firstorder derivative of hydraulic conductivity (i.e., equal to specific discharge) with respect to the water content [22, 28, 29]. These studies illustrated the celerity as “mobility of water in soil” [29, 30], and the mechanisms of pressure propagation in unsaturated zone was, however, not discussed.
On the other hand, the pore water velocity in numerous flow paths is dictated by the pore structure such as pore connectivity, tortuosity, and pore size distribution [1–3]. Especially, the soil hydraulic functions can be derived as an integration of pore water velocities in all flow paths if considers a unit gradient condition [10, 12, 31]. Based on a pore bundle model, the hydraulic properties of each flow path can be equivalent to that of a cylinder tube with a small circular cross section and a sufficient length [32–34], then the experimentally measured water retention curve is used to infer pore size distribution and pore water velocity distribution [35–38]. Under such condition, each tube conducts the viscous flow; therefore, the relations of the equivalent tube radius to capillarity pressure and pore hydraulic conductivity can be, respectively, described by a capillary rise equation (i.e., YoungLaplace equation) and a pipe/cylinder liquid flow equation (i.e., the HagenPoiseuille law). The hydraulic conductivity functions integrated from pore velocity distribution can indirectly reflect the porescale hydraulic properties [39, 40]. As the firstorder derivative of hydraulic conductivity can be defined as a celerity function for unsaturated soil [22, 24], the celerity can be related to the pore water velocity distribution.
The complex solute transport behaviours are caused by the pore water velocity distribution in a subsurface flow system [28, 41–43], which have been captured either by an experimentally obtained breakthrough curve [5, 44–46] or by natural and artificial tracers in a hydrological system [28, 47, 48]. The breakthrough curve expresses the time series of the tracer concentration of drainage flow, which was often characterized by a coexisting early initial breakthrough and extensive tailing [44]. The tracer experiment in a hillslope or a catchment assists in inferring the travel time distribution [49, 50]. However, many tracer studies found a paradox behaviour; that is, during a highintensity rainstorm a large amount of “old” water resident in soil or groundwater is expelled into the stream, while the labelled “new” water of the rainfall on the hillslope immediately starts to appear in the stream after a fast infiltration through a fast flow path [51–54]. Those complex tracer transport behaviours address the importance of quantifying the pore water velocity distribution for achieving the more detailed physical interpretations of unsaturated flow and transport [28, 55].
The celerity and velocity are intimately linked with tracer transport, while most of the current studies focused on one or two aspects. For example, mathematical derivations for the celerity and kinematic ratio have been provided by Rasmussen et al. [22], while they did not find the relation between tracer transport and pressure propagation. Wang et al. [24] and Mohammadi et al. [56] found the derivation of the tracer breakthrough curve through the soil hydraulic conductivity functions, and their method have been validated through their experiments. However, the studies from Wang et al. [24] and Mohammadi et al. [56] were based on the soil hydraulic functions of either the BrooksCorey model or the van Genuchten model for a singlepermeability model, and more soil hydraulic functions still need further studies. There is no report existing that explored the relations among celerity (pressure propagation), velocity (water transport), and tracer transport (pore velocity distribution) yet.
The objective of this paper was to (i) illustrate the concept of celerity in unsaturated flow, (ii) propose a new soil hydraulic function (i.e., celerity function) to express the pore water velocity distribution, and (iii) quantify the advective breakthrough curves using the proposed celerity function. Specifically, in Section 2.1, the definitions of velocity and celerity were given to illustrate the mechanisms of pressure propagation in unsaturated soils. Based on a pore bundle model, Section 2.2 mathematically proved the equivalence between celerity in unsaturated soil and maximum pore water velocity among all flow paths that conduct water. Furthermore, a theoretical analysis was conducted to illustrate the difference between celerity and velocity of wetting front under a unithydraulicgradient condition, showing that the celerity function derived from a perturbation analysis can reveal the pore water velocity distribution (see Section 2.3). The celerity functions of the BrooksCorey model and both unmodified and modified Mualemvan Genuchten models are provided in Section 2.4. The impact of velocity distribution on tracer transport can be illustrated by mathematically deriving a breakthrough curve through a celerity function (Section 2.5). Functions for the celerity in dualpermeability models are derived in Section 2.6. Section 3 analysed the hydraulic characteristic and breakthrough curves of typical soil textures by using the proposed celerity function. Hereto, the pore water velocity distribution and breakthrough curve of typical soils were analysed using the existing parameter sets of the Mualemvan Genuchten model and the BrooksCorey model. Additionally, the equifinal parameter sets of a bimodal soil hydraulic function were distinguished by the celerity function. The followup discussion focused on the mechanism of pressure propagation and tracer transport in a natural subsurface flow system.
2. Theory
2.1. Definitions
In subsurface, the vertical component (LT^{−1}) of a specific discharge vector (i.e., the volume flux of water per unit crosssectional area) is formulated using Darcy’s law as where (LT^{−1}) is the hydraulic conductivity, (L) is the pressure head, and (L) is the vertical coordinate (positive downward).
The average pore water velocity (LT^{−1}) is where (−) is the volumetric water content and (−) is the residual water content.
The continuity equation for onedimensional vertical flow is
The derivative of specific discharge with respect to is written as
Substitution of (4) for in (3) results in an advection equation where (LT^{−1}) is the celerity [22]
Theoretically, the advection equation represents the advection of the soilwater content with speed . When is a constant value, it would mean that an arbitrarily shaped pulse of moves with constant speed without changing its shape. The celerity is, however, not a constant but a function of water content. Hence, the celerity is an approximate advection speed of a small pulse of .
A ratio between celerity and average velocity is called the kinematic ratio (−), as defined by Rasmussen et al. [22]:
2.2. Celerity and Maximum Pore Velocity
The pore sizes and their intrinsic permeabilities vary throughout the unsaturated zone. Based on a pore bundle model, many previous studies considered soil as a bundle of nonintersecting, parallel, and cylindrical tubes with varying radii to quantify the water and solute transport [32–34, 57]. When water enters soil flowing through a pore of certain size, each pore is either entirely filled with water or entirely empty. Considering variably saturated conditions, pores are filled from the smallest tube group () to the largest (). Under unsaturated conditions, of the tube groups are filled with water. The pore water velocity (LT^{−1}) in one specific equivalent (saturated) cylinder tube can be determined by its specific discharge where (LT^{−1}) is a coefficient relating the head gradient and the average velocity in pore .
The pores with the same size and the same pore water velocity can be classified as one group, based on which the water content (i.e., all the volumetric fractions of the pore space filled with water) is discretised in fragments of equal water content and ordered from small pores to larger pores. Assuming a unique, nonhysteretic relationship between water content and capillary pressure, water can flow through the waterfilled saturated smaller pores with higher capillary pressure, while the larger pores will remain empty and inactivated for conducting unsaturated flow [58].
By discretizing the flow domain into numerous flow tubes that are vertically homogeneous, the hydraulic characteristics of each tube can be formulated as functions of the volumetric water content . Considering the pore size distribution, the pore water flow velocity follows a distribution ordering from smallest pore group (with velocity ) to largest pore group (with velocity ). At a certain water content value, (LT^{−1}) may denote the maximum pore water flow velocity taking place in the waterfilled pore with the largest equivalent radius—then the specific discharge (LT^{−1}) and the hydraulic conductivity (LT^{−1}) can be obtained through a summation: where is the largest tube group filled with water. Substitution of (8) for in (9) gives where the hydraulic conductivity is defined as
The unit hydraulic gradient condition considers a gravitydriven vertical infiltrating flow (i.e., when and ), and the hydraulic conductivity function therefore can be obtained by summing up the pore water velocity from a limited number of pore groups as shown in Figure 1.
The integral equation of (9) is and the integral equation of (11) is
The celerity is now obtained as where is the maximum velocity corresponding to the water content . Hence, the celerity is equal to the maximum velocity corresponding to a certain water content .
Many previous studies using the pore bundle model to derive the hydraulic conductivity function by upscaling pore water velocities to REV scale were similar to (12)–(13) [32–34], while an inverse derivation in (14) demonstrates that the celerity in unsaturated flow is equal to the maximum pore water flow velocity. Moreover, Figure 1 implies that the piecewise linear approximation of the soil hydraulic conductivity function can infer the flow paths (or tubes groups) in which the water flow has distinct velocities.
2.3. Expressing Pore Water Velocity Distribution by a Celerity Function
The effects of pore water velocity distribution on soil hydrology can be visualized by a perturbation analysis with a hypothetical infiltration experiment. Considering soil under unit hydraulic gradient condition as (when is constant), the specific discharge is equal to the hydraulic conductivity () [58]. Under this condition, an increment of water flow on the surface boundary will cause downward pressure propagation. The pressure propagation may be linked to the celerity as maximum pore water velocity, while it can also be described by a similar terminology of the velocity of wetting front. A hypothetical analysis is provided hereafter to exemplify their different roles in quantifying the pressure wave propagation under unit hydraulic gradient condition.
Considering a steady flow condition driven by an infiltration flow , the increments of to the upper boundary correspondingly induce an increment of soil water content that generates a pistonshaped advancing wetting front with velocity (LT^{−1}) expressed as [59]
The velocity of wetting front expresses the average pore water velocity of the newly activated pore system with the volumetric fraction of , which is induced by a change of the specific discharge .
On the other hand, the celerity in unsaturated flow can be illustrated by a perturbation analysis; namely, a small increment of at the surface boundary will activate water flow in a largersize pore with water content of :
Furthermore, under unit hydraulic gradient condition, the celerity can be mathematically linked to the maximum pore hydraulic conductivity
A perturbation (i.e., an small increment) of water flow at the surface boundary will activate water flow in the “new” tube group with larger size and higher pore hydraulic conductivity , and the flow velocity in the newly activated water flow path will be larger than other waterfilled tubes. Assuming that the fluid dynamics are independent among tubes, the celerity under variably saturated condition express the onetoone correspondence between the maximum pore water velocities with the water content [58]. When exerting the perturbation analysis in saturated/unsaturated soil, the celerity as a firstorder derivative of the hydraulic conductivity function can quantify the pore hydraulic conductivity distribution (equivalent to the pore water velocity distribution) under unit hydraulic gradient condition.
2.4. Different Soil Hydraulic Models
The pore water velocity distribution can be obtained by a celerity function that was derived as the firstorder derivative of the hydraulic conductivity function. The two most widely used soil hydraulic conductivity functions of the BrooksCorey model and the Mualemvan Genuchten model [35, 37, 60] have been chosen by Rasmussen et al. [22] to formulate celerity functions. However, the kinematic ratio and celerity under the Mualemvan Genuchten model in nearsaturated soil is approaching an infinite value, which is caused by a mathematical artefact [38]. To overcome this problem, the modified Mualemvan Genuchten model was used to prevent the unrealistic variations in unsaturated hydraulic conductivity for the nearsaturation range [38, 61–63].
The mathematical formulations of water retention curve, hydraulic conductivity, average pore velocity, celerity, and kinematic ratio based on the BrooksCorey and Modified Mualemvan Genuchten parameterization are shown in Table 1 [35, 38]. The modified Mualemvan Genuchten model here includes a nonzero minimum capillary height of air entry pressure to overcome the mathematical artefact, and the results will be given in Section 3.1. The unmodified Mualemvan Genuchten model sets the air entry pressure to zero, so that the parameter equals one. The formulation of the kinematic ratio derived from the BrooksCorey model is a constant value of that equals to the power index in the hydraulic conductivity function.
 
Notation of table: (L^{−1}), (−), and (−) are the fitting parameters for the BrooksCorey model (detonated with subscription of “BC”) and the Mualemvan Genuchten model (detonated with subscription of “VG”); is the effective saturation (−) that is defined as ; is the correction factor to modify the vanGenuchten model , with ; is the minimum capillary height, which becomes zero when ; and is the pore connectivity parameter and is usually assumed to be 0.5. 
The kinematic ratio is a constant value if deriving it from a power function that describes unsaturated hydraulic conductivity [64, 65], and its value usually falls in a range between 2.5 and 24.5 [65]. Similarly, the celerity in surface water has been extensively investigated for describing fluid dynamics [13–15, 17]. The kinematic ratio for surface flow is approximately equal to 1.67 derived from the Manning resistance coefficient in a kinematic wave equation for openchannel flow [15]. Since the celerity in unsaturated subsurface flow is analogous to the wave celerity in surface flow [8], the kinematic ratio of unsaturated flow can be 1.5–14.7 times larger than that of surface flow, which implies that the velocity distribution in unsaturated subsurface flow is highly nonlinear.
2.5. Derivation of a Breakthrough Curve through the Celerity Function
Transport of a conservative tracer in a porous medium is governed by physical mechanisms of advection, dispersion, and molecular diffusion, among which advection and dispersion are a function of the velocity distribution. The tracer advection is governed by the average pore water velocity, and the tracer dispersion is caused by the velocity distribution [66, 67].
Considering a onedimensional and uniform flow under the unit hydraulic gradient condition, the transport of conservative tracer in a subsurface flow system can be described with the pore bundle model, and the breakthrough curve can be derived from velocity distribution [58]. If the water content is constant and equal to , of the tube groups are filled with water, so that the water content may be written as
Considering a soil column with length (L) with a uniform distributed water content , the specific discharge is equal to the specified constant water flow at the upper boundary under unit hydraulic gradient condition. Considering the varying pore water velocities in numerous vertical flow paths in a soil column, the time needed for tracer breakthrough from a pore group with a velocity can be expressed as
Tracer breakthrough sequentially occurs from the largest waterfilled pores to the smallest one. Especially, the first arrival time of the tracer transport is associated with the maximum pore water velocity and thus can be calculated with .
When the tracer concentration in the surface flow is changing from to , the tracer breakthrough will be achieved at the time of , and then the concentration of the drainage flow will be increased gradually. At time (after the first arrival time), water in all tube groups has travelled from the top of the column to the bottom of the column; the concentration of drainage flow can be formulated through a summation [57]:
The integral equivalent of (21) is where where is pore water flow velocity in a certain flow path at the arrival time of for tracer breakthrough
is the water content corresponding to a flow velocity that is equal to the celerity and thus
Under unit hydraulic gradient condition, the integrals in (21) are equal to the specific discharge (see (12)), which are equal to the hydraulic conductivity, so that it becomes
The mathematical expression for the breakthrough curve based on the Mualemvan Genuchten parameterization requires an analytical approach [56]. In this study, the numerical approach was used to derive the breakthrough curve of the unmodified and modified Mualemvan Genuchten models. We hereafter demonstrate an example of deriving the breakthrough curve of a nonreactive tracer from the celerity function based on the BrooksCorey model, to be more explicitly linking the tracer transport with the pore water velocity distribution [58].
Substitution of the celerity equation for the BrooksCorey model in Table 1 into (25) and rearrangement of terms giving the pore water velocity can be calculated as
If the initial tracer concentration is zero, the breakthrough curve for the BrooksCorey model can be written explicitly as where approximates the first arrival time of the tracer. Eq. (27) indicates that the first arrival time of a tracer is determined by the maximum pore water velocity, and the tailing is determined by the pore water velocity distribution. The abovegiven mathematical derivations are similar with the previous work that also formulated the breakthrough curve by soil hydraulic conductivity functions under the saturated condition [57]. In this study, we extended the formulation to unsaturated soils and further formulated the tracer transport as a function of kinematic ratio that simplified (28) by replacing: with as a dimensionless time defined as
When is equal to 1, the flow through a soil column is equal to one pore volume.
The breakthrough curve expressed with (28) shows that the kinematic ratio can directly determine the first arrival time and tailing of a conservative tracer in drainage flow. For an extreme case, when the kinematic ratio is equal to 1, the pore water velocity in all the flow paths is equal to the average pore water velocity. Such a uniform pore water velocity distribution leads to a pistonshaped breakthrough curve:
Eq. (30) expresses only the tracer advection that is controlled by the average pore water velocity, while excluding the dispersion effect caused by the pore water velocity distribution.
For the BrooksCorey model, substituting the hydraulic conductivity function (from Table 1) into (29), the dimensionless time can be formulated as a function of either soil water content or specific discharge:
The kinematic ratio is a constant in the BrooksCorey model (see Table 1), which means that the breakthrough curve as a function of dimensionless time in variably saturated soil is independent of the specific discharge or effective saturation for the BrooksCorey model.
2.6. Celerity of the Bimodal Soil Hydraulic Functions
This section derives a celerity function through the bimodal soil hydraulic conductivity functions for quantifying the pore water velocity distribution in structured soil. Conceptualizing a soil porous medium as two overlapping continuums of the matrix domain and the preferential flow domain as a dualcontinuum pore system [41, 68], the composite bimodal soil hydraulic functions use two water retention curves and hydraulic conductivity functions to, respectively, describe the soil hydraulic properties [69–74]. The preferential flow domain consists of the pores with relatively large size (often taken as larger than 0.3 mm in equivalent tube diameter) and low tortuosity, such as worm burrows, root channels, tension cracks, and interaggregate pores [41, 75, 76]. The remaining micropore volume is classified as the matrix domain. The bimodal soil hydraulic function has been widely used to parameterize the dualpermeability models that simulate two groups of state variables (i.e., different pressure heads, water contents, and tracer concentrations between two domains) for describing the nonequilibrium phenomena [77]. The hydraulic characteristics of the preferential flow domain and the matrix domain can be distinguished by their parameterization in bimodal soil hydraulic functions [70, 78]. Below, we derive the soil hydraulic functions of average pore water velocity and celerity for a dualcontinuum system under unit hydraulic gradient condition.
First, the volumetric ratio of the preferential flow domain and the matrix flow domain sums up to 1: where the subscripts and denote the preferential flow and matrix flow domains, respectively. The total water content and total specific discharge of a soil with a dualcontinuum pore system are calculated as the weighted average of the two domains water content and specific discharge:
The average pore velocity of each a domain is defined as
Then, the average pore velocity of the dualcontinuum pore system is
Under the steady flow condition, the celerity of each a domain can be defined as follows:
Finally, the celerity of the dualcontinuum pore system can be expressed as
When the matrix domain is unsaturated, the matrix conducts water, and the preferential flow domain remains inactive. Therefore, pressure propagation will take place in the matrix domain as well. When the matrix domain approaches nearsaturation condition, the preferential flow domain becomes activated. At this stage, the water flow in the preferential domain quickly takes over the majority of water flow and consequently the pressure propagation of celerity.
3. Analysis
3.1. Soil Hydraulic Functions under a SingleDomain Conceptualization
This section combines the theoretical formulation (given in Section 2) with the existing parameterization of 5 typical soils to analyse the soil hydraulic behaviours. Table 2 shows the standard parameter sets of 5 typical soils obtained from the UNSODA database [79, 80]. The detailed values of each soil type have slightly different values of saturated hydraulic conductivity and residual/saturated water content, which might be related to the different fitting algorithms and sampling data. Considering that the values of these parameters are drastically varied for different soil textures, we then decided to use dimensionless indicators, for example, the effective saturation. Furthermore, the values of specific discharge, average velocity, and celerity were all divided by the saturated hydraulic conductivity to convert those variables to dimensionless values. Consequently, all the values of each individual indicator for the different soil types fall in a similar range (Figure 2). The air entry pressure values in the modified Mualemvan Genuchten model were adopted from those of the BrooksCorey model for all 5 typical soils and listed in the last column of Table 2. Solely using either specific discharge or average pore water velocity to quantify the transport processes is theoretically insufficient [20, 22, 81]; we included the celerity and kinematic ratio in Table 2 to provide a more complete illustration of the dependence of pore water flow velocity and pressure propagation on soil effective saturation.
 
Note. is the air entry pressure value of the BrooksCorey model and the modified Mualemvan Genuchten model. 
The water retention curves (first row in Figure 2) for 5 typical soils have distinct curvatures attributed to the different pore size distributions [38, 76]. Specifically, in coarse texture soil (e.g., sandy soil), a flatter curve might manifest a relatively uniform pore size distribution and homogeneous soil hydraulic characteristic of the coarse texture. When approaching the fully saturated state, the slopes of the water retention curve from the Mualemvan Genuchten model reaches an infinite value; instead, the slopes from the BrooksCorey model and the modified Mualemvan Genuchten model (Figures 2(a) and 2(c)) are close to 0, which are attributed to the inclusion of an air entry pressure head stemmed from the observation of various types of soils [80].
Celerity, average pore water velocity, and specific discharge all increase as the effective saturation increases, but the ranges of their values differ. Specifically, ranges from 0 to 1 as expected. is between 0 and 3 (because is around 0.35; see Table 2). ranges from 0 to 70. also ranges from 0 to 30, except for the unmodified Mualemvan Genuchten model under nearsaturated condition. All four dimensionless variables reach their maximum value when soils are saturated. The celerity curves can manifest the maximum pore water velocity at a certain saturation state, and Figure 2 shows that the celerity is over 20 times larger than the average pore water velocity and 10 up to 100 times larger than the specific discharge.
The celerity curves show different patterns for different soil types. Near saturation, can reach above 50 for fine textured soils, while it reaches around 20 for coarsetextured soils. On the contrary, when the effective saturation drops below 0.8–0.9, is much smaller. The value of is below 0.5 for for all soil types. The value of is highest for sand (= 0.5) and lowest for clay (= 0.2), which means that more than 50% of the specific discharge flows through only 15% of the pore space when flow is saturated.
The curves are nearexponentially increased with the effective saturation increasing till the soils reach the fully saturated state. However, the celerity curves show different patterns under the parameters set for different soil types. Under nearsaturated condition, in finetextured soil (i.e., sandy clay and clay), the celerity responds to the effective saturation more significantly than that in coarse soil textures. in coarsetextured soil (e.g., sandy soil) only changes from 0 to around 15 (Figures 2(i) and 2(j)), while the celerity curve of fine texture soil has a much larger range (e.g., of silty soil ranges from 0 to 70). Under the low saturation state (e.g., ), specific discharge and average pore water velocity behave similarly with the relative low values between different soil textures. Till the effective saturation reaches 0.9, is below 0.5 in all soil types (sand and clay ), which implicitly manifests that more than 50% of specific discharge is conducted in pores with a volumetric fraction less than 10%.
The slope of the curve is used as an indicator, shown in Figure 3. The kinematic ratios derived from the unmodified Mualemvan Genuchten model approach the infinite values, which were induced by a mathematical artefact. However, the kinematic ratios under the modified Mualemvan Genuchten model and the BrooksCorey model consistently have (nearly) constant values for all the saturation ranges (specific values are given in Table 3). indicates that the celerity is proportionally increasing with the average pore water velocity (see Figure 3(a)).
(a)
(b)
(c)

3.2. Breakthrough Curves
The pore water velocity distribution can further dictate the tracer transport, which can be formulated and illustrated with a breakthrough curve (Figure 4). The functions of the breakthrough curve in typical soils used the same parameter set given in Table 2. Breakthrough curves were plotted for three soil hydraulic models, using an analytic approach for the BrooksCorey model and a numerical approach for the unmodified and modified vanGenuchten models. We adopted a dimensionless time (see (29)) to exclude the impact of , , and on breakthrough curve. It is recalled that only advective transport is considered based on the pore bundle model, where water particles remain in the same pore group (and hence travel with the same velocity) for the entire length of the column. Here, we used dimensionless time to eliminate the influence of different specific discharges in various soil textures, and the kinematic ratio is directly correlated with the shape of breakthrough curves (following (28)).
The generated breakthrough curves are presented in Figure 4 for saturated conditions and Figures 5 and 6 for unsaturated conditions. The pistonshaped breakthrough curves (black lines in Figure 4, based on (31)) were computed for tracer advection driven by flow with uniform velocity distribution. The dimensionless first arrival time can be determined as a reciprocal of the kinematic ratio (see (28)). Specifically, the kinematic ratio is lower in coarsetextured soil than in finetextured soil; therefore, the dimensionless first arrival time of finetextured soil is earlier. The first arrival time described by the BrooksCorey model and the modified Mualemvan Genuchten model is significantly larger than that calculated by the unmodified Mualemvan Genuchten model for saturated condition. The kinematic ratio under the unmodified Mualemvan Genuchten model is expressed as an infinite value under saturated condition, which is caused by the overestimated celerity leading to a nearly instant breakthrough of the tracer.
After the first arrival time, the increment of relative concentration over the dimensionless time in each soil is caused by pore water velocity distribution. When , the relative concentration is above 0.7 for all soil types and all soil hydraulic models. All breakthrough curves show very long tailing, caused by the low velocity in the smaller pore bundles. When , the relative concentration approximately ranges between 0 and 0.7; the water particle transport process is dominated by the pore water velocity that is larger than the average pore water velocity. When , the relative concentration generally approaches 1.0, and the slopes of the breakthrough curves turn to be flatter. The long tailing is determined by the pore water velocity distribution, because a significant fraction of pores with the low pore water velocity needs much longer time for tracer breakthrough.
In the BrooksCorey model, the kinematic ratio is constant, and the breakthrough curves are independent of the effective saturation. In the modified Mualemvan Genuchten model, the shapes of the breakthrough curves for saturated soil are affected by the value of the air entry pressure, while the shapes of the breakthrough curves for unsaturated soil are not affected by the effective saturation.
The effect of the air entry pressure on the breakthrough curves for sandy loam soil under saturated and unsaturated () conditions is shown in Figure 5. Under saturated conditions (Figure 5(a)), the dimensionless first arrival time is zero for the unmodified Mualemvan Genuchten model (), and it approaches 0.2 when the air entry pressure is increased from 3 cm to 30 cm. In contrast, breakthrough curves do not depend significantly on the air entry pressure in unsaturated sandy loam (Figure 5(b)).
Breakthrough curves are shown for different effective saturation values for sandy loam under the unmodified Mualemvan Genuchten model in Figure 6. Except for the saturated case (), the other three breakthrough curves are very similar.
3.3. Analysis of the DualPermeability Model
The bimodal soil hydraulic function describes the hydraulic properties of a conceptualized dualcontinuum pore system as a composite function. However, parameterization of a bimodal soil hydraulic function is difficult as two conceptualized domains cannot be experimentally separated. Fitting the parameter set of the bimodal soil hydraulic function according to the measurable singlecontinuum soil hydraulic function of the water retention curve and soil hydraulic conductivity might result in the “equifinal” parameter sets. For example, based on the Mualemvan Genuchten model, Köhne et al. [70] obtained 5 sets of the equifinal parameters of bimodal soil hydraulic functions, all of which can perfectly fit a singlecontinuum water retention curve and hydraulic conductivity function. Different soil hydraulic behaviours described by those 5 parameter sets are difficult to distinguish. Therefore, this section demonstrates a way to distinguish the soil hydraulic characteristics of the equifinal parameter sets in the bimodal soil hydraulic function, by quantifying the celerity as a function of pressure head.
Table 4 shows the hydraulic properties of each pore domain specified by the Mualemvan Genuchten model under two “equifinal” parameter settings [70]. Under unit hydraulic gradient condition, it is often assumed that the pressure heads in the matrix flow domain and in the preferential flow domain are the same when facilitating the illustration of the hydraulic behaviour of structured soil [70]. For a dualcontinuum conceptualization, the hydraulic characteristics in the matrix domain, preferential flow domain, and total domain were obtained according to (32)–(38).

Figure 7(a) shows the water retention curve for the bulk soil; the fitted bimodal water retention curve using two groups of parameters ( and ) for a dualcontinuum porous medium (shown as “Total” in Figure 7) agrees well with that of the singlecontinuum soil hydraulic model (shown as “Single” in Figure 7). Figures 7(b) and 7(c) further show that the specific discharge and average pore water velocity in a bimodal hydraulic conductivity function under two groups of “equifinality” parameter sets are different from that in singlecontinuum soil hydraulic function. However, the parameter sets of a bimodal soil hydraulic function generate an equifinal phenomenon in the formulated total hydraulic conductivity and total average pore water velocity. Thus, under the equifinal parameter sets of the bimodal soil hydraulic model, the higher hydraulic conductivity of the preferential flow domain can be compensated by a smaller volumetric fraction of for conducting the equivalent amount of total specific discharge.
(a) Water content
(b) Hydraulic conductivity
(c) Average pore velocity
(d) Celerity
The proposed constitutive relation between celerity and pressure head (in Figure 7(d)) is able to distinguish the implicitly formulated pore water velocity distribution of equifinal parameter sets. The celerity in the matrix domain is nearly the same for the two parameter sets. When the pressure head is smaller than −35 cm, the celerity is controlled by the matrix flow, and the preferential flow domain can be regarded as an insignificant component contributing to the pressure propagation, while when the pressure head is larger than −35 cm, the celerity is in turn controlled by the preferential flow with a much higher pore water velocity than the matrix flow, and the pressure propagation is induced by the preferential flow with much higher pore water velocity than that of the matrix flow.
The relative difference of the hydraulic conductivity and the celerity for the total domain is shown in Figure 8. For the two parameter sets, the relative difference of the hydraulic conductivity is less than 5% (Figure 8(a)). However, the celerity for the second parameter set () in is approximately twice as large as with the first parameter set () (Figure 8(b)). This is an important result, as it can be decided which of the two parameter sets is the better one, if the celerity can be measured. The celerity can be measured with a tracer experiment by using a conservative tracer and measuring the first arrival time of the tracer. In summary, the proposed constitutive relation between celerity and pressure head (i.e., ) is able to distinguish the implicitly formulated pore water velocity distribution of equifinal parameter sets.
(a) Hydraulic conductivity
(b) Celerity
4. Discussion
4.1. Relations of Soil Hydraulic Properties to the Breakthrough Curve
The pore water velocity distribution is feasible by measuring the variation of both specific discharge and water content under a perturbation analysis (i.e., by generally increasing the water fluxes on the soil surface) [33]. The pore bundle model is linking the measurable REVscale variable with the porescale hydraulic characteristics. On such a basis, the proposed celerity function enables the derivation of the breakthrough curve to express the influence of pore water velocity distribution on tracer transport. Consequently, the mathematical derivation (in Section 2) implies that the same parameter set for the water retention curve or hydraulic conductivity function can also be used to parameterize pore water velocity distribution and the breakthrough curve. Inversely, the measured breakthrough curve might also be used to infer the parameter for soil hydraulic function. Recent studies show that the breakthrough curve derived from a tracer experiment could assist in obtaining a reasonably accurate parameter estimation for describing the water retention curve and the hydraulic conductivity function [56, 57]. Thus, mutually fitting the breakthrough curve, water retention curve, and hydraulic conductivity function by an inverse parameter estimation through experiments can provide a more accurately mathematical formulation to express the soil hydraulic properties.
The obtained pore water velocity distribution in typical soil textures is highly nonlinear [31, 38], which manifests that the nonuniform pore size distribution influences both water and tracer transport. In the breakthrough curve, the first arrival time is dictated by the maximum pore water velocity, and the long tailing is dictated by the relatively low pore water velocity. This phenomenon implies that the transport processes are significantly different in the preferential flow path within the matrix. A large volumetric fraction of water stored in micropores is nearly stagnant or with extremely low velocity [50], being classified as the matrix domain, which causes an extensive tailing of tracer concentration in the breakthrough curve. The matrix flow has a relatively insignificant contribution in delivering specific discharge and propagating pressure wave under nearsaturated condition, while flow occurs in a small volumetric fraction of pores (e.g., macropores) that can be categorized as the preferential flow path dominating the transport processes and pressure propagation.
4.2. Implementing the Kinematic Ratio to Predict the Maximum Tracer Velocity
The fast tracer transport in heterogeneous soil will be dominated by the preferential flow with the high flow velocity; thus, the effect of absorption, reaction, or diffusion might be negligible. Nimmo [23] analysed 64 field experiments and determined that the maximum tracer velocities varied within a small range, which could be predicted with by a simple model. These tracer experiments were conducted in various types of soil or bedrock with a transport distance ranging from 0.3 to 1300 m, and the maximum tracer velocities (with a 90% probability) spread between 0.8 and 200 m/day. Nimmo [23] proposed a predictive model using a ratio to link the fastest tracer velocity with an effective precipitation rate (a spatially and temporarily averaged precipitation applied to the surface boundary) and suggested that the ratio of is an essentially constant value which equaled to 18 based on an empirical estimation from their geometric mean. The obtained ratio with an orderofmagnitude accuracy facilitates a fast prediction of worstcase contaminate travel times. Nimmo illustrated that the low variability of the ratio could be caused by a natural speed limit of the preferential flow in terms of frictional force and water exchange (e.g., via absorption) from the macropore to the matrix. However, the predicting ability of Nimmo’s model might be hampered due to lack of the parameter constrains.
Based on our analyses of celerity, a theoretical proof and parameter constrain for Nimmo’s model can be achieved by using a kinematic ratio to predict the maximum tracer velocity. Considering a hillslope, the process of water particle and tracer transport can be summarized as follows. For vertical infiltration flow in an unsaturated zone, celerity represents the maximum pore water velocity as well as the maximum tracer velocity, and the kinematic ratio approximately represents the ratio between maximum tracer velocity and average tracer velocity. For lateral flow in a saturated zone, celerity represents a nearly instant pressure propagation, yet the actual flow velocity in each pore is related with the intrinsic permeability depending on the size and tortuosity of the pore [31, 38]. Therefore, the celerity function derived for unsaturated flow can express the pore water velocity distribution even in saturated soil. In Section 3.1, we found that the kinematic ratio of unsaturated flow is likely to be a constant value for all ranges of saturation. Therefore, it can be used as an indicator to quantify the ratio between maximum tracer velocity and average tracer velocity for both the saturated soil and unsaturated soil.
The presented analysis of celerity in unsaturated soils can be compared to Nimmo’s model. The ratio in Nimmo’s paper has the same physical meaning as the ratio of that can be defined in this study. The value of is related to both kinematic ratio and soil water content as . The kinematic ratio of unsaturated flow based on the BrooksCorey model is constant and ranges from 3 to 16 for the soils of Table 2, and is equal to under unit gradient condition.
The ratio is plotted versus the effective saturation in Figures 9(a)–9(c). The value of decreases with for the BrooksCorey and the modified van Genuchten model. The value of approaches infinity when the effective saturation decreases from 0.3 to 0, which is not shown here because the corresponding specific discharge is very small. Except for the unmodified van Genuchten model, is fairly constant and only weakly depends on the specific discharge () as shown in Figures 9(d)–9(f). Assuming that the soilwater content in the natural system is within a range of 0.2–0.5 (effective saturation is within a range of 0.4–1.0), we can also derive which approximately ranges between 6 and 80. The geometric mean of is 23 based on this rough estimation, which is close to the value of 18 from Nimmo’s experimental finding. Moreover, the value of the kinematic ratio can be approximately inferred by the soil textures, water retention curve, soil hydraulic function, or maximum tracer velocity, and using the estimated kinematic ratio can provide an additional parameter constrain.
5. Conclusion
We conducted a theoretical study attempting to investigate the mechanism of pressure propagation and tracer transport in a subsurface flow system by analyzing the celerity in unsaturated flow. The celerity was firstly defined according to the continuum equations. Furthermore, based on the conceptualization of a pore bundle model, a mathematical derivation showed that the celerity in unsaturated flow can be illustrated as the maximum pore water velocity among all the waterfilled flow paths, and the kinematic ratio can be illustrated as the ratio of the maximum pore water velocity to the average pore water velocity. Under unit hydraulic gradient condition, the celerity function is formulated as a firstorder derivative of the soil hydraulic conductivity function to manifest the pore water velocity distribution, which can further be used to derive a breakthrough curve for tracking the tracer and water particle transport under steady flow condition.
The soil hydraulic characteristic of typical soil textures is reanalysed by using the (un/modified) Mualemvan Genuchten and BrooksCorey models with standard parameter sets. The results of all soil textures present nonlinear relations of effective saturation with hydraulic conductivity and celerity. The results show that water in a small volumetric fraction (around 15%) of pores has a much higher velocity than the remaining pore volume. The derived breakthrough curve can manifest the influence of pore water velocity distribution on the advective tracer transport. The first arrival time of the tracer is determined by the maximum pore water velocity, and the long tailing is caused by the nearly stagnate matrix flow in micropores with a large volumetric fraction. Furthermore, the pore water velocity distribution formulated in a composite bimodal soil hydraulic function can also be described by the proposed celerity function. Different parameter sets may result in similar water retention curves and soil hydraulic functions, but their celerities differ significantly which shows that the equifinal parameter sets of a bimodal soil hydraulic function can be distinguished by using the celerity function.
Analyzing the celerity in unsaturated flow facilitates the understanding and quantification of the porescale hydraulic characteristics and the complex subsurface transport processes. The celerityeffective saturation curve may be applicable to the bimodal soil hydraulic functions to demonstrate the distributions of pore water velocities for natural soils. Furthermore, the following possible approaches for implementation are suggested and discussed in Section 4. Under unit hydraulic gradient, the three constitutive relations (i.e., water retention curve , hydraulic conductivity function , and celerity function ) and the derived breakthrough curve share the same parameter set, which could facilitate the formulation and parameterization of soil hydraulic models. Finally, we found that the ratio of celerity and specific discharge is fairly constant for specific soil texture, in accordance with published experiment findings.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was financially supported by the National Natural Science Foundation of China (Grant nos. 41771021 and 41471012), the China Postdoctoral Science Foundation (nos. 2017M621783 and 2018T110527), the International Postdoctoral Exchange Fellowship Program by the China Postdoctoral Council (year 2017), and the Startup Foundation for Introducing Talent of Nanjing University of Information Science and Technology (no. 2017r045). The authors would also like to thank Professor Huub Savenije, Professor Mark Bakker, and Dr. Thom Bogaard for their comments and suggestions for this manuscript.
References
 J. Bear, Dynamics of Fluids in Porous Media, Dover, 1988.
 K. Beven, “Changing ideas in hydrology — the case of physicallybased models,” Journal of Hydrology, vol. 105, no. 12, pp. 157–172, 1989. View at: Publisher Site  Google Scholar
 D. R. Nielsen, M. Th. Van Genuchten, and J. W. Biggar, “Water flow and solute transport processes in the unsaturated zone,” Water Resources Research, vol. 22, no. 9S, pp. 89S–108S, 1986. View at: Publisher Site  Google Scholar
 X. Bai, B. He, X. Li et al., “First assessment of sentinel1A data for surface soil moisture estimations using a coupled water cloud model and advanced integral equation model over the Tibetan plateau,” Remote Sensing, vol. 9, no. 7, p. 714, 2017. View at: Publisher Site  Google Scholar
 W. Durner and H. Flühler, “Multidomain model for poresize dependent transport of solutes in soils,” Geoderma, vol. 70, no. 2–4, pp. 281–297, 1996. View at: Publisher Site  Google Scholar
 L. Wei, J. Dong, M. Gao, and X. Chen, “Factors controlling temporal stability of surface soil moisture: a watershedscale modeling study,” Vadose Zone Journal, vol. 16, no. 10, 2017. View at: Publisher Site  Google Scholar
 H. Flühler, W. Durner, and M. Flury, “Lateral solute mixing processes — a key for understanding fieldscale transport of water and solutes,” Geoderma, vol. 70, no. 2–4, pp. 165–183, 1996. View at: Publisher Site  Google Scholar
 V. P. Singh, Kinematic Wave Modeling in Water Resources, Environmental Hydrology, John Wiley & Sons, 1997.
 P. D'Odorico, S. Fagherazzi, and R. Rigon, “Potential for landsliding: dependence on hyetograph characteristics,” Journal of Geophysical Research: Earth Surface, vol. 110, 2005. View at: Publisher Site  Google Scholar
 S. K. Kampf and S. J. Burges, “A framework for classifying and comparing distributed hillslope and catchment hydrologic models,” Water Resources Research, vol. 43, no. 5, 2007. View at: Publisher Site  Google Scholar
 O. Cainelli, A. Bellin, and M. Putti, “On the accuracy of classic numerical schemes for modeling flow in saturated heterogeneous formations,” Advances in Water Resources, vol. 47, pp. 43–55, 2012. View at: Publisher Site  Google Scholar
 M. van Genuchten, C. NaveiraCotta, T. Skaggs, A. Raoof, and E. Pontedeiro, “The Use of Numerical Flow and Transport Models in Environmental Analyses,” in Application of Soil Physics in Environmental Analyses, W. G. Teixeira, M. B. Ceddia, M. V. Ottoni, and G. K. Donnagema, Eds., pp. 349–376, Springer International Publishing, 2014. View at: Publisher Site  Google Scholar
 M. J. Lighthill and G. B. Whitham, “On kinematic waves. I. Flood movement in long rivers,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 229, no. 1178, pp. 281–316, 1955. View at: Publisher Site  Google Scholar
 V. P. Singh, Kinematic Wave Modeling in Water Resources, SurfaceWater Hydrology, John Wiley & Sons, 1996.
 T. S. W. Wong, “Influence of upstream inflow on wave celerity and time to equilibrium on an overland plane,” Hydrological Sciences Journal, vol. 41, no. 2, pp. 195–205, 1996. View at: Publisher Site  Google Scholar
 H. H. G. Savenije, “3  Tidal Dynamics,” in Salinity and Tides in Alluvial Estuaries, H. H. G. Savenije, Ed., pp. 69–107, Elsevier Science Ltd, Amsterdam, 2005. View at: Google Scholar
 H. H. G. Savenije and E. J. M. Veling, “Relation between tidal damping and wave celerity in estuaries,” Journal of Geophysical Research: Oceans, vol. 110, no. C4, 2005. View at: Publisher Site  Google Scholar
 P. F. Germann and K. Beven, “Kinematic wave approximation to infiltration into soils with sorbing macropores,” Water Resources Research, vol. 21, no. 7, pp. 990–996, 1985. View at: Publisher Site  Google Scholar
 V. P. Singh and E. S. Joseph, “Kinematicwave model for soilmoisture movement with plantroot extraction,” Irrigation Science, vol. 14, no. 4, pp. 189–198, 1994. View at: Publisher Site  Google Scholar
 J. J. McDonnell and K. Beven, “Debates—the future of hydrological sciences: a (common) path forward? A call to action aimed at understanding velocities, celerities and residence time distributions of the headwater hydrograph,” Water Resources Research, vol. 50, no. 6, pp. 5342–5350, 2014. View at: Publisher Site  Google Scholar
 J. Davies and K. Beven, “Comparison of a multiple interacting pathways model with a classical kinematic wave subsurface flow solution,” Hydrological Sciences Journal, vol. 57, no. 2, pp. 203–216, 2012. View at: Publisher Site  Google Scholar
 T. C. Rasmussen, R. H. Baldwin, J. F. Dowd, and A. G. Williams, “Tracer vs. pressure wave velocities through unsaturated saprolite,” Soil Science Society of America Journal, vol. 64, no. 1, pp. 75–85, 2000. View at: Publisher Site  Google Scholar
 J. R. Nimmo, “Simple predictions of maximum transport rate in unsaturated soil and rock,” Water Resources Research, vol. 43, no. 5, 2007. View at: Publisher Site  Google Scholar
 Z. Wang, J. Feyen, D. R. Nielsen, and M. T. van Genuchten, “Twophase flow infiltration equations accounting for air entrapment effects,” Water Resources Research, vol. 33, no. 12, pp. 2759–2767, 1997. View at: Publisher Site  Google Scholar
 Y. Su, J. Langhammer, and J. Jarsjö, “Geochemical responses of forested catchments to bark beetle infestation: evidence from high frequency instream electrical conductivity monitoring,” Journal of Hydrology, vol. 550, pp. 635–649, 2017. View at: Publisher Site  Google Scholar
 A. Alaoui, P. Germann, N. Jarvis, and M. Acutis, “Dualporosity and kinematic wave approaches to assess the degree of preferential flow in an unsaturated soil,” Hydrological Sciences Journal, vol. 48, no. 3, pp. 455–472, 2003. View at: Publisher Site  Google Scholar
 R. E. Smith, “Approximate soil water movement by kinematic characteristics,” Soil Science Society of America Journal, vol. 47, no. 1, pp. 3–8, 1983. View at: Publisher Site  Google Scholar
 K. Beven and P. Germann, “Macropores and water flow in soils revisited,” Water Resources Research, vol. 49, no. 6, pp. 3071–3092, 2013. View at: Publisher Site  Google Scholar
 R. J. Charbeneau, “Kinematic models for soil moisture and solute transport,” Water Resources Research, vol. 20, no. 6, pp. 699–706, 1984. View at: Publisher Site  Google Scholar
 V. P. Singh, “Is hydrology kinematic?” Hydrological Processes, vol. 16, no. 3, pp. 667–716, 2002. View at: Publisher Site  Google Scholar
 A. Peters, W. Durner, and G. Wessolek, “Consistent parameter constraints for soil hydraulic functions,” Advances in Water Resources, vol. 34, no. 10, pp. 1352–1365, 2011. View at: Publisher Site  Google Scholar
 N. T. Burdine, “Relative permeability calculations from pore size distribution data,” Journal of Petroleum Technology, vol. 5, no. 3, pp. 71–78, 1953. View at: Publisher Site  Google Scholar
 E. C. Childs and N. CollisGeorge, “The permeability of porous materials,” Proceedings of the Royal Society A Mathematical, Physical and Engineering Sciences, vol. 201, no. 1066, pp. 392–405, 1950. View at: Publisher Site  Google Scholar
 Y. Mualem, “A new model for predicting the hydraulic conductivity of unsaturated porous media,” Water Resources Research, vol. 12, no. 3, pp. 513–522, 1976. View at: Publisher Site  Google Scholar
 R. H. Brooks and A. T. Corey, Hydraulic Properties of Porous Media, Colorado State University, 1964.
 M. Tuller and D. Or, “Hydraulic conductivity of variably saturated porous media: film and corner flow in angular pore space,” Water Resources Research, vol. 37, no. 5, pp. 1257–1276, 2001. View at: Publisher Site  Google Scholar
 M. T. van Genuchten, “A closedform equation for predicting the hydraulic conductivity of unsaturated soils,” Soil Science Society of America Journal, vol. 44, no. 5, pp. 892–898, 1980. View at: Publisher Site  Google Scholar
 H. J. Vogel, “A numerical experiment on pore size, pore connectivity, water retention, permeability, and solute transport using network models,” European Journal of Soil Science, vol. 51, no. 1, pp. 99–105, 2000. View at: Publisher Site  Google Scholar
 M. A. Celia, P. C. Reeves, and L. A. Ferrand, “Recent advances in pore scale models for multiphase flow in porous media,” Reviews of Geophysics, vol. 33, no. S2, pp. 1049–1057, 1995. View at: Publisher Site  Google Scholar
 H. Dahle, M. Celia, and S. M. Hassanizadeh, “BundleofTubes Model for Calculating Dynamic Effects in the CapillaryPressureSaturation Relationship,” in Upscaling Multiphase Flow in Porous Media, D. B. Das and S. M. Hassanizadeh, Eds., pp. 5–22, Springer Netherlands, 2005. View at: Publisher Site  Google Scholar
 K. Beven and P. Germann, “Macropores and water flow in soils,” Water Resources Research, vol. 18, no. 5, pp. 1311–1325, 1982. View at: Publisher Site  Google Scholar
 J. Freer, J. J. McDonnell, K. J. Beven et al., “The role of bedrock topography on subsurface storm flow,” Water Resources Research, vol. 38, no. 12, pp. 51–516, 2002. View at: Publisher Site  Google Scholar
 T. Uchida, Y. Asano, T. Mizuyama, and J. J. McDonnell, “Role of upslope soil pore pressure on lateral subsurface storm flow dynamics,” Water Resources Research, vol. 40, no. 12, 2004. View at: Publisher Site  Google Scholar
 M. L. Brusseau and P. S. C. Rao, “Modeling solute transport in structured soils: a review,” Geoderma, vol. 46, no. 1–3, pp. 169–192, 1990. View at: Publisher Site  Google Scholar
 F. A. Dullien, Porous Media: Fluid Transport and Pore Structure, Academic press, 1991.
 A. G. Williams, J. F. Dowd, D. Scholefield, N. M. Holden, and L. K. Deeks, “Preferential flow variability in a wellstructured soil,” Soil Science Society of America Journal, vol. 67, no. 4, pp. 1272–1281, 2003. View at: Publisher Site  Google Scholar
 C. Kendall and J. J. McDonnell, “Preface,” in Isotope Tracers in Catchment Hydrology, C. K. J. McDonnell, Ed., pp. vii–vix, Elsevier, Amsterdam, 1998. View at: Publisher Site  Google Scholar
 J. M. Köhne, S. Köhne, and J. Šimůnek, “A review of model applications for structured soils: (a) water flow and tracer transport,” Journal of Contaminant Hydrology, vol. 104, no. 1–4, pp. 4–35, 2009. View at: Publisher Site  Google Scholar
 J. Beckers and Y. Alila, “A model of rapid preferential hillslope runoff contributions to peak flow generation in a temperate rain forest watershed,” Water Resources Research, vol. 40, no. 3, 2004. View at: Publisher Site  Google Scholar
 J. Renée Brooks, H. R. Barnard, R. Coulombe, and J. J. McDonnell, “Ecohydrologic separation of water between trees and streams in a Mediterranean climate,” Nature Geoscience, vol. 3, no. 2, pp. 100–104, 2010. View at: Publisher Site  Google Scholar
 M. Hrachowitz, H. H. G. Savenije, G. Blöschl et al., “A decade of predictions in ungauged basins (PUB)—a review,” Hydrological Sciences Journal, vol. 58, no. 6, pp. 1198–1255, 2013. View at: Publisher Site  Google Scholar
 J. W. Kirchner, “A double paradox in catchment hydrology and geochemistry,” Hydrological Processes, vol. 17, no. 4, pp. 871–874, 2003. View at: Publisher Site  Google Scholar
 F. M. Phillips, “Hydrology: soilwater bypass,” Nature Geoscience, vol. 3, no. 2, pp. 7778, 2010. View at: Publisher Site  Google Scholar
 M. Weiler, B. L. McGlynn, K. J. McGuire, and J. J. McDonnell, “How does rainfall become runoff? A combined tracer and runoff transfer function approach,” Water Resources Research, vol. 39, no. 11, 2003. View at: Publisher Site  Google Scholar
 J. Seibert, T. Grabs, S. Köhler, H. Laudon, M. Winterdahl, and K. Bishop, “Linking soil and streamwater chemistry based on a riparian flowconcentration integration model,” Hydrology and Earth System Sciences, vol. 13, no. 12, pp. 2287–2297, 2009. View at: Publisher Site  Google Scholar
 M. H. Mohammadi, M. R. Neishabouri, and H. Rafahi, “Predicting the solute breakthrough curve from soil hydraulic properties,” Soil Science, vol. 174, no. 3, pp. 165–173, 2009. View at: Publisher Site  Google Scholar
 Q. Wang, R. Horton, and J. Lee, “A simple model relating soil water characteristic curve and soil solute breakthrough curve 1,” Soil Science, vol. 167, no. 7, pp. 436–443, 2002. View at: Publisher Site  Google Scholar
 W. Shao, Numerical Modeling of the Effect of Preferential Flow on Hillslope Hydrology and Slope Stability, Doctoral thesis, TU Delft Water Resources, 2017.
 J. L. Nieber, R. Z. Dautov, A. G. Egorov, and A. Y. Sheshukov, “Dynamic capillary pressure mechanism for instability in gravitydriven flows; review and extension to very dry conditions,” Transport in Porous Media, vol. 58, no. 12, pp. 147–172, 2005. View at: Publisher Site  Google Scholar
 S. C. Iden and W. Durner, “Freeform estimation of the unsaturated soil hydraulic properties by inverse modeling using global optimization,” Water Resources Research, vol. 43, no. 7, 2007. View at: Publisher Site  Google Scholar
 P. Nasta, J. A. Vrugt, and N. Romano, “Prediction of the saturated hydraulic conductivity from Brooks and Corey’s water retention parameters,” Water Resources Research, vol. 49, no. 5, pp. 2918–2925, 2013. View at: Publisher Site  Google Scholar
 A. Peters and W. Durner, “A simple model for describing hydraulic conductivity in unsaturated porous media accounting for film and capillary flow,” Water Resources Research, vol. 44, no. 11, 2008. View at: Publisher Site  Google Scholar
 M. G. Schaap and M. T. van Genuchten, “A modified Mualem–van Genuchten formulation for improved description of the hydraulic conductivity near saturation,” Vadose Zone Journal, vol. 5, no. 1, pp. 27–34, 2006. View at: Publisher Site  Google Scholar
 G. S. Campbell, “A simple method for determining unsaturated conductivity from moisture retention data,” Soil Science, vol. 117, no. 6, pp. 311–314, 1974. View at: Publisher Site  Google Scholar
 Y. Mualem and G. Dagan, “Hydraulic conductivity of soils: unified approach to the statistical models,” Soil Science Society of America Journal, vol. 42, no. 3, pp. 392–395, 1978. View at: Publisher Site  Google Scholar
 D. Jacques, J. Šimůnek, A. Timmerman, and J. Feyen, “Calibration of Richards' and convection–dispersion equations to fieldscale water flow and solute transport under rainfall conditions,” Journal of Hydrology, vol. 259, no. 1–4, pp. 15–31, 2002. View at: Publisher Site  Google Scholar
 J. Vanderborght, M. Vanclooster, A. Timmerman et al., “Overview of inert tracer experiments in key belgian soil types: relation between transport and soil morphological and hydraulic properties,” Water Resources Research, vol. 37, no. 12, pp. 2873–2888, 2001. View at: Publisher Site  Google Scholar
 W. Shao, T. Bogaard, M. Bakker, and M. Berti, “The influence of preferential flow on pressure propagation and landslide triggering of the Rocca Pitigliana landslide,” Journal of Hydrology, vol. 543, pp. 360–372, 2016. View at: Publisher Site  Google Scholar
 W. Durner, “Hydraulic conductivity estimation for soils with heterogeneous pore structure,” Water Resources Research, vol. 30, no. 2, pp. 211–223, 1994. View at: Publisher Site  Google Scholar
 J. M. Köhne, S. Köhne, and H. H. Gerke, “Estimating the hydraulic functions of dualpermeability models from bulk soil data,” Water Resources Research, vol. 38, no. 7, pp. 261–2611, 2002. View at: Publisher Site  Google Scholar
 D. Mallants, P.H. Tseng, N. Toride, A. Tinunerman, and J. Feyen, “Evaluation of multimodal hydraulic functions in characterizing a heterogeneous field soil,” Journal of Hydrology, vol. 195, no. 1–4, pp. 172–199, 1997. View at: Publisher Site  Google Scholar
 B. P. Mohanty, R. S. Bowman, J. M. H. Hendrickx, and M. T. van Genuchten, “New piecewisecontinuous hydraulic functions for modeling preferential flow in an intermittentfloodirrigated field,” Water Resources Research, vol. 33, no. 9, pp. 2049–2063, 1997. View at: Publisher Site  Google Scholar
 N. Romano, P. Nasta, G. Severino, and J. W. Hopmans, “Using bimodal lognormal functions to describe soil hydraulic properties,” Soil Science Society of America Journal, vol. 75, no. 2, pp. 468–480, 2011. View at: Publisher Site  Google Scholar
 T. Zurmühl and W. Durner, “Determination of parameters for bimodal hydraulic functions by inverse modeling,” Soil Science Society of America Journal, vol. 62, no. 4, pp. 874–880, 1998. View at: Publisher Site  Google Scholar
 J. M. Hendrickx and M. Flury, “Uniform and Preferential Flow Mechanisms in the Vadose Zone,” in Conceptual Models of Flow and Transport in the Fractured Vadose Zone, pp. 149–188, National Academies Press, Washington, DC, USA, 2001. View at: Google Scholar
 N. J. Jarvis, “A review of nonequilibrium water flow and solute transport in soil macropores: principles, controlling factors and consequences for water quality,” European Journal of Soil Science, vol. 58, no. 3, pp. 523–546, 2007. View at: Publisher Site  Google Scholar
 W. Shao, J. Ni, A. K. Leung, Y. Su, and C. W. W. Ng, “Analysis of plant root–induced preferential flow and porewater pressure variation by a dualpermeability model,” Canadian Geotechnical Journal, vol. 54, no. 11, pp. 1537–1552, 2017. View at: Publisher Site  Google Scholar
 W. Shao, T. A. Bogaard, M. Bakker, and R. Greco, “Quantification of the influence of preferential flow on slope stability using a numerical modelling approach,” Hydrology and Earth System Sciences, vol. 19, no. 5, pp. 2197–2212, 2015. View at: Publisher Site  Google Scholar
 F. J. Leij, The UNSODA Unsaturated Soil Hydraulic Database: User’s Manual, National Risk Management Research Laboratory, Office of Research and Development, US Environmental Protection Agency, 96, 1996.
 M. T. Van Genuchten, F. Leij, and S. Yates, The RETC Code for Quantifying the Hydraulic Functions of Unsaturated Soils, Robert S. Kerr Environmental Research Laboratory, 1991.
 G. W. Waswa, A. D. Clulow, C. Freese, P. A. L. Le Roux, and S. A. Lorentz, “Transient pressure waves in the vadose zone and the rapid water table response,” Vadose Zone Journal, vol. 12, no. 1, 2013. View at: Publisher Site  Google Scholar
Copyright
Copyright © 2018 Wei Shao et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.