Quantify the Pore Water Velocity Distribution by a Celerity Function

Key Laboratory of Meteorological Disaster, Ministry of Education / Joint International Research Laboratory of Climate and Environment Change / Collaborative Innovation Centre on Forecast and Evaluation of Meteorological Disasters / School of Hydrology and Water Resources, Nanjing University of Information Science and Technology, Nanjing, Jiangsu 210044, China Water Resources Section, Faculty of Civil Engineering and Geosciences, Delft University of Technology, Stevinweg 1, PO Box 5048, 2600 GA Delft, Netherlands Department of Physical Geography and Geoecology, Faculty of Science, Charles University, Albertov 6, 128 43 Prague 2, Czech Republic Key Laboratory of Mountain Hazards and Surface Process, Institute of Mountain Hazards and Environment, Chinese Academy of Sciences, Chengdu 610041, China


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][2][3].In hydrological studies, it is impractical to deterministically represent the pore-scale fluid dynamics and transport processes in soil, for the variability of pore water velocity is difficult to observe or predict [4][5][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 well-known Darcy-Richards 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 advection-diffusion equation to simulate solute/thermal transport [10][11][12].Applications of the Darcy-Richards 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][14][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 perturbation-induced 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 nearly-instant 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 cause-effect 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 first-order 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][2][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][33][34], then the experimentally measured water retention curve is used to infer pore size distribution and pore water velocity distribution [35][36][37][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., Young-Laplace equation) and a pipe/cylinder liquid flow equation (i.e., the Hagen-Poiseuille law).The hydraulic conductivity functions integrated from pore velocity distribution can indirectly reflect the pore-scale hydraulic properties [39,40].As the first-order 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][42][43], which have been captured either by an experimentally obtained breakthrough curve [5,[44][45][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 high-intensity 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][52][53][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 Brooks-Corey model or the van Genuchten model for a single-permeability 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 unithydraulic-gradient 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 Brooks-Corey model and both unmodified and modified Mualem-van 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 dual-permeability 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 Mualem-van Genuchten model and the Brooks-Corey model.Additionally, the equifinal parameter sets of a bimodal soil hydraulic function were distinguished by the celerity function.The follow-up discussion focused on the mechanism of pressure propagation and tracer transport in a natural subsurface flow system.

Theory
2.1.Definitions.In subsurface, the vertical component q (LT −1 ) of a specific discharge vector (i.e., the volume flux of water per unit cross-sectional area) is formulated using Darcy's law as where K (LT −1 ) is the hydraulic conductivity, h (L) is the pressure head, and z (L) is the vertical coordinate (positive downward).The average pore water velocity v (LT −1 ) is where θ (−) is the volumetric water content and θ r (−) is the residual water content.
The continuity equation for one-dimensional vertical flow is The derivative of specific discharge q with respect to z is written as Substitution of (4) for ∂q/∂z in (3) results in an advection equation where c (LT −1 ) is the celerity [22] c = ∂q ∂θ 6 Theoretically, the advection equation represents the advection of the soil-water content θ with speed c.When c is a constant value, it would mean that an arbitrarily shaped pulse of θ moves with constant speed c 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 α K (−), 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][33][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 (i = 1) to the largest (i = N).Under unsaturated conditions, Mof the N tube groups are filled with water.The pore water velocity v i (LT −1 ) in one specific equivalent (saturated) cylinder tube i can be determined by its specific discharge where k i (LT −1 ) is a coefficient relating the head gradient and the average velocity in pore i.
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 N 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 water-filled 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 v i follows a distribution ordering from smallest pore group (with velocity v 1 ) to largest pore group (with velocity v M ).At a certain water content value, v θ (LT −1 ) may denote the maximum pore water flow velocity taking place in the water-filled pore with the largest equivalent radius-then the specific discharge q (LT −1 ) and the hydraulic conductivity K (LT −1 ) can be obtained through a summation: 3 Geofluids where M is the largest tube group filled with water.Substitution of (8) for v i in (9) gives where the hydraulic conductivity K is defined as The unit hydraulic gradient condition considers a gravity-driven vertical infiltrating flow (i.e., when ∇h = 0 and K = q), and the hydraulic conductivity function K θ therefore can be obtained by summing up the pore water velocity v i 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 v θ 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][33][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.

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 ∂ h − z /∂z = −1 (when h is constant), the specific discharge is equal to the hydraulic conductivity (q = K) [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 q, the increments of Δq to the upper boundary correspondingly induce an increment of soil water content Δθ that generates a piston-shaped advancing wetting front with velocity v WF (LT −1 ) expressed as [59] The velocity of wetting front v WF 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 Δq.
On the other hand, the celerity in unsaturated flow can be illustrated by a perturbation analysis; namely, a small increment of δq at the surface boundary will activate water flow in a larger-size 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 k θ , and the flow velocity v θ in the newly 4 Geofluids activated water flow path will be larger than other water-filled tubes.Assuming that the fluid dynamics are independent among tubes, the celerity under variably saturated condition express the one-to-one 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 first-order 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 first-order derivative of the hydraulic conductivity function.The two most widely used soil hydraulic conductivity functions of the Brooks-Corey model and the Mualem-van 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 Mualem-van Genuchten model in near-saturated soil is approaching an infinite value, which is caused by a mathematical artefact [38].To overcome this problem, the modified Mualem-van Genuchten model was used to prevent the unrealistic variations in unsaturated hydraulic conductivity for the near-saturation range [38,[61][62][63].
The mathematical formulations of water retention curve, hydraulic conductivity, average pore velocity, celerity, and kinematic ratio based on the Brooks-Corey and Modified Mualem-van Genuchten parameterization are shown in Table 1 [35,38].The modified Mualem-van 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 Mualem-van Genuchten model sets the air entry pressure h s to zero, so that the parameter ε equals one.The formulation of the kinematic ratio derived from the Brooks-Corey model is a constant value of 2/n BC + 3 that equals to the power index in the hydraulic conductivity function.
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][14][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 open-channel 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.

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 one-dimensional 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 θ w , M of the N tube groups are filled with water, so that the water content θ w may be written as , n (−), and m (−) are the fitting parameters for the Brooks-Corey model (detonated with subscription of "BC") and the Mualem-van Genuchten model (detonated with subscription of "VG"); Θ is the effective saturation (−) that is defined as h s is the minimum capillary height, which becomes zero when θ m = θ s ; and l VG is the pore connectivity parameter and is usually assumed to be 0.5.

Geofluids
Considering a soil column with length L (L) with a uniform distributed water content θ w , 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 t i needed for tracer breakthrough from a pore group with a velocity v i can be expressed as Tracer breakthrough sequentially occurs from the largest water-filled pores to the smallest one.Especially, the first arrival time of the tracer transport is associated with the maximum pore water velocity L/v θ w and thus can be calculated with L/c θ w .
When the tracer concentration in the surface flow is changing from C 0 to C 0 + ΔC, the tracer breakthrough will be achieved at the time of L/c θ w , and then the concentration of the drainage flow will be increased gradually.At time t i (after the first arrival time), water in all tube groups j > i 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

21
where where v c is pore water flow velocity in a certain flow path at the arrival time of t c for tracer breakthrough θ c is the water content corresponding to a flow velocity v c 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 C t based on the Mualem-van 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 Mualem-van Genuchten models.We hereafter demonstrate an example of deriving the breakthrough curve C t c of a nonreactive tracer from the celerity function based on the Brooks-Corey model, to be more explicitly linking the tracer transport with the pore water velocity distribution [58].
Substitution of the celerity equation for the Brooks-Corey model in Table 1 into (25) and rearrangement of terms giving the pore water velocity v θ c can be calculated as If the initial tracer concentration is zero, the breakthrough curve for the Brooks-Corey model can be written explicitly as where L/c θ w 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 above-given 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 α K that simplified (28) by replacingn BC : with t * as a dimensionless time defined as When t * 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 6 Geofluids average pore water velocity.Such a uniform pore water velocity distribution leads to a piston-shaped 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 Brooks-Corey 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 Brooks-Corey model (see Table 1), which means that the breakthrough curve as a function of dimensionless time t * in variably saturated soil is independent of the specific discharge or effective saturation for the Brooks-Corey 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 dual-continuum 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][70][71][72][73][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 dual-permeability 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 dual-continuum 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 f and m denote the preferential flow and matrix flow domains, respectively.The total water content and total specific discharge of a soil with a dual-continuum 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 dual-continuum pore system is Under the steady flow condition, the celerity of each a domain can be defined as follows: Finally, the celerity of the dual-continuum 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 near-saturation 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.

Soil Hydraulic Functions under a Single-Domain
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 7 Geofluids values in the modified Mualem-van Genuchten model were adopted from those of the Brooks-Corey 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.
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 Mualem-van Genuchten model reaches an infinite value; instead, the slopes from the Brooks-Corey model and the modified Mualem-van 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, q/K s ranges from 0 to 1 as expected.v/K s is between 0 and 3 (because θ s − θ r is around 0.35; see Table 2).c/K s ranges from 0 to 70. α K also ranges from 0 to 30, except for the unmodified Mualem-van Genuchten model under near-saturated 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, c/K s 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,c/K s is much smaller.The value of K/K s is below 0.5 for θ = 0 85 for all soil types.The value of K/K s 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 K/K s curves are near-exponentially 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 near-saturated condition, in fine-textured soil (i.e., sandy clay and clay), the celerity responds to the effective saturation more significantly than that in coarse soil textures.c/K s in coarse-textured 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., c/K s of silty soil ranges from 0 to 70).Under the low saturation state (e.g., Θ < 0 7), 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, K/K s is below 0.5 in all soil types (sand K/K s = 0 5 and clay K/K s = 0 2), 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 α K curve is used as an indicator, shown in Figure 3.The kinematic ratios derived from the unmodified Mualem-van Genuchten model approach the infinite values, which were induced by a mathematical artefact.However, the kinematic ratios under the modified Mualem-van Genuchten model and the Brooks-Corey model consistently have (nearly) constant values for all the saturation ranges (specific values are given in Table 3).α K indicates that the celerity is proportionally increasing with the average pore water velocity (see Figure 3(a)).

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 Brooks-Corey model and a numerical approach for the unmodified and modified van-Genuchten models.We adopted a dimensionless time (see (29)) to exclude the impact of K s , θ s , and θ r 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 t * 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   9 Geofluids unsaturated conditions.The piston-shaped 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 coarse-textured soil than in fine-textured soil; therefore, the dimensionless first arrival time of fine-textured soil is earlier.The first arrival time described by the Brooks-Corey model and the modified Mualem-van Genuchten model is significantly larger than that calculated by the unmodified Mualem-van Genuchten model for saturated condition.The kinematic ratio under the unmodified Mualem-van 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 t * = 1, 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  10 Geofluids bundles.When t * < 1, the relative concentration C t * /ΔC 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 t * > 1, the relative concentration C t * /ΔC 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 Brooks-Corey model, the kinematic ratio is constant, and the breakthrough curves are independent of the effective saturation.In the modified Mualem-van 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 (θ = 0 8) conditions is shown in Figure 5.Under saturated conditions (Figure 5(a)), the dimensionless first arrival time is zero for the unmodified Mualem-van Genuchten model (h s = 0), 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 Mualem-van Genuchten model in Figure 6.Except for the saturated case (Θ = 1 0), the other three breakthrough curves are very similar.

Analysis of the Dual-Permeability
Model.The bimodal soil hydraulic function describes the hydraulic properties of a conceptualized dual-continuum 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 single-continuum soil hydraulic function of the water retention curve and soil hydraulic conductivity might result in the "equifinal" parameter sets.For example, based on the Mualem-van 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 Mualem-van 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 dual-continuum 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 (w = 0 1 and w = 0 025) for a dualcontinuum porous medium (shown as "Total" in Figure 7) agrees well with that of the single-continuum 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 single-continuum 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 w for conducting the equivalent amount of total specific discharge.
The proposed constitutive relation between celerity and pressure head c h (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 (w = 0 025) in c h is approximately twice as large as with the first parameter set (w = 0 1) (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., c h ) is able to distinguish the implicitly formulated pore water velocity distribution of equifinal parameter sets.

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 REV-scale variable with the pore-scale 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 13 Geofluids 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 near-saturated 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 v 0 with an effective precipitation rate i 0 (a spatially and temporarily averaged precipitation applied to the surface boundary) and suggested that the ratio of v 0 /i 0 is an essentially constant value which equaled to 18 based on an empirical estimation from their geometric mean.The obtained  14 Geofluids ratio with an order-of-magnitude accuracy facilitates a fast prediction of worst-case contaminate travel times.Nimmo illustrated that the low variability of the ratio v 0 /i 0 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 v 0 /i 0 in Nimmo's paper has the same physical meaning as the ratio of c/q that can be defined in this study.The value of c/q is related to both kinematic ratio and soil water content as c/q = c/ vθ = α K /θ.The kinematic ratio of unsaturated flow based on the Brooks-Corey model is constant and ranges from 3 to 16 for the soils of Table 2, and q is equal to K under unit gradient condition.
The ratio c/K is plotted versus the effective saturation Θ in Figures 9(a The value of c/K 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, c/K is fairly constant and only weakly depends on the specific discharge (q/K s ) as shown in Figures 9(d)-9(f).Assuming that the soil-water 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 c/K which approximately ranges between 6 and 80.The geometric mean of c/K 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.

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 water-filled 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 first-order 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) Mualem-van Genuchten and Brooks-Corey models with standard parameter sets.The results of all soil textures present nonlinear relations of effective saturation with hydraulic conductivity and celerity.3.

Geofluids
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 pore-scale hydraulic characteristics and the complex subsurface transport processes.The celerity-effective 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 θ h , hydraulic conductivity function K θ , and celerity function c θ ) 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 c/q is fairly constant for specific soil texture, in accordance with published experiment findings.

Figure 1 :
Figure 1: A numerical integration of pore water velocities from pore groups.

Figure 2 :
Figure2: Constitutive relations of the effective saturation to logarithmic pressure head (a-c), dimensionless specific discharge (d-f), dimensionless average pore water velocity (g-i), dimensionless celerity (j-l), and kinematic ratio (m-o) in 5 typical soils under models of Brook-Corey, unmodified Mualem-van Genuchten, and modified Mualem-van Genuchten (illustrated in each column).

Figure 4 :
Figure 4: Breakthrough curves and the piston-shaped breakthrough curves under the unit hydraulic gradient condition for 5 typical saturated soils in (a) the Brook-Corey model, (b) the unmodified Mualem-van Genuchten model, and (c) the modified Mualem-van Genuchten model.

Figure 5 :1 1 CFigure 6 :
Figure 5: Effect of air entry pressure and effective saturation on breakthrough curves for sandy loam soil in the (a) unmodified Mualem-van Genuchten model and (b) the modified Mualem-van Genuchten model.

c
(cm/day) h (cm) w f = 0.1 matrix w f = 0.025 matrix w f = 0.1 PF w f = 0.025 PF w f = 0.1 total w f = 0.025 total (d) Celerity

Figure 7 :
Figure 7: Under parameter settings of a single-continuum soil hydraulic model in Table 4, bimodal soil hydraulic functions of (a) water retention θ h , (b) hydraulic conductivity K h , (c) average pore velocity v h , and (d) celerity c h for matrix flow (matrix), preferential flow (PF), single-continuum model, and dual-continuum model (dual) are given, respectively.
)-9(c).The value of c/K decreases with Θ for the Brooks-Corey and the modified van Genuchten model.

Figure 8 :
Figure 8: Relative difference of (a) unsaturated hydraulic conductivity and (b) celerity of total soils under the two parameter sets of Table3.

Table 1 :
Constitutive relationships under the unit hydraulic gradient condition formulated by the Brooks-Corey model and the modified Mualem-van Genuchten model.

Table 2 :
Parameter sets of 5 typical soils in the Brooks-Corey model and the van-Genuchten model.
Note. h s is the air entry pressure value of the Brooks-Corey model and the modified Mualem-van Genuchten model.

Table 4 :
Equifinal parameter sets of the bimodal soil hydraulic function for a hypothetical single-continuum soil.