Articles | Volume 27, issue 4
https://doi.org/10.5194/hess-27-895-2023
https://doi.org/10.5194/hess-27-895-2023
Research article
 | 
24 Feb 2023
Research article |  | 24 Feb 2023

Mixed formulation for an easy and robust numerical computation of sorptivity

Laurent Lassabatere, Pierre-Emmanuel Peyneau, Deniz Yilmaz, Joseph Pollacco, Jesús Fernández-Gálvez, Borja Latorre, David Moret-Fernández, Simone Di Prima, Mehdi Rahmati, Ryan D. Stewart, Majdi Abou Najm, Claude Hammecker, and Rafael Angulo-Jaramillo
Abstract

Sorptivity is one of the most important parameters for the quantification of water infiltration into soils. Parlange (1975) proposed a specific formulation to derive sorptivity as a function of the soil water retention and hydraulic conductivity functions, as well as initial and final soil water contents. However, this formulation requires the integration of a function involving hydraulic diffusivity, which may be undefined or present numerical difficulties that cause numerical misestimations. In this study, we propose a mixed formulation that scales sorptivity and splits the integrals into two parts: the first term involves the scaled degree of saturation, while the second involves the scaled water pressure head. The new mixed formulation is shown to be robust and well-suited to any type of hydraulic function – even with infinite hydraulic diffusivity or positive air-entry water pressure heads – and any boundary condition, including infinite initial water pressure head, h-. Lastly, we show the benefits of using the proposed formulation for modeling water into soil with analytical models that use sorptivity.

1 Introduction

Soil sorptivity represents the capacity of soil to absorb water by capillarity (Cook and Minasny2011). The accurate estimation of soil sorptivity is crucial for the modeling of water infiltration into soils and the hydraulic characterization of soils (Angulo-Jaramillo et al.2016; Stewart and Abou Najm2018). Several models and methods make use of this variable, such as the Beerkan Estimation of Soil Transfer parameters (BEST) methods (Lassabatere et al.2006; Yilmaz et al.2010; Bagarello et al.2014a; Lassabatere et al.2019; Angulo-Jaramillo et al.2019) and related simplified Beerkan approaches (Bagarello et al.2014b; Di Prima et al.2020; Yilmaz2021). Sorptivity is also required for the computation of several hydraulic parameters, like the macroscopic capillary length (Bouwer1964; White and Sully1987).

The squared sorptivity is related to the flux-concentration function, F(θ), as follows (Philip and Knight1974):

(1) S 2 θ 0 , θ 1 = 2 θ 0 θ 1 θ - θ 0 F ( θ ) D ( θ ) d θ ,

where D(θ)=K(θ)dh/dθ is the hydraulic diffusivity function, and θ0 and θ1 stand for the initial and final water contents. In the context of water infiltration into soils, the initial water content, θ0, refers to the water content along the soil profile before water infiltration, and the final water content, θ1, corresponds to that imposed at the soil surface (at the water source). Several studies have investigated the definition of the flux-concentration functions, depending on the type of soils (Angulo-Jaramillo et al.2016, Table 1, p. 33). Ross et al. (1996) suggested the use of the approximation proposed by Parlange (1975) for most soils leading to two main forms of sorptivity, a diffusivity form, SD, and a conductivity form, SK, both summarized below:

(2)F(θ)=2θ-θ0θ1+θ-2θ0(3)SD2θ0,θ1=θ0θ1θ1+θ-2θ0D(θ)dθ(4)SDθ0,θ1=θ0θ1θ1+θ-2θ0D(θ)dθ(5)SKh0,h1=h0h1θ1+θ(h)-2θ0K(h)dh,

where the initial and the final values of the water pressure heads, h0 and h1, correspond to the water contents θ0=θ(h0) and θ1=θ(h1).

These two forms, the diffusivity form, SD, and the conductivity form, SK, each have their own shortcomings. For certain hydraulic models, D(θ) tends towards infinity when θ0θs, making it difficult to compute the right-hand side of Eq. (4). Moreover, when the surface water pressure head exceeds the air-entry water pressure head, ha, SD misses the saturated part of sorptivity, hah1(θ1+θ(h)-2θ0)K(h)dh (Ross et al.1996). The conductivity form SK must be used when it is necessary to account for the two parts of sorptivity, i.e., the unsaturated and saturated parts, as indicated by the following relationship (Lassabatere et al.2021):

(6) S K h 0 , h 1 h a = S D 2 θ 0 , θ s + 2 θ s - θ 0 h 1 - h a K s = S D 2 θ h 0 , θ h 1 + 2 θ h 1 - θ h 0 h 1 - h a K s .

We can thus conclude that the conductivity form, SK, is the most general equation. However, SK can also be difficult to handle when the initial conditions are very dry. In particular, for very dry initial conditions, the initial water pressure head corresponding to θr corresponds to h0-. Then, the calculation of SK requires the evaluation of an integral that involves an infinite lower bound: -h1θh1+θ(h)-2θrK(h)dh.

In this study, we propose a new mixed formulation that overcomes these problems. We compare it to the approaches commonly used to compute sorptivity, i.e., Eqs. (4) and (5). The proposed mixed formulation automatically accounts for the saturated and unsaturated parts of sorptivity. It also allows for easy computation under any initial condition, including the extreme case of an initial water content equal to the residual water content, θ0=θr (corresponding to a negative initial water pressure head that tends towards infinity, h0-), and a final water pressure head higher than the air-entry water pressure head, h1ha, even including positive values. In addition to proposing a new robust formulation for the computation of sorptivity, we aim to demonstrate the following points: (i) the computation of sorptivity with classic approaches may be challenging depending on the hydraulic models chosen for describing the water retention (WR) and hydraulic conductivity (HC) functions; (ii) the proposed mixed formulation is an ideal estimator for sorptivity and performs well at all times; (iii) the usual methods, based on the use of SD (Eq. 4) or SK (Eq. 5), do not necessarily provide accurate estimations of the nominal sorptivity; and (iv) those misestimations of sorptivity may have substantial impacts on the prediction of water infiltration into soils when inserting sorptivity into analytical models.

The paper is organized as follows. The Theory section presents the proposed mixed formulation. Next, the paper analyzes the precision of the mixed formulation by comparing it with the exact analytical formulation for the case of the maximum sorptivity, S(-,0)=-0θs+θ(h)-2θrK(h)dh. The maximum sorptivity, S(-,0), encompasses the two types of problems, i.e., infinite negative initial water pressure head and infinite diffusivity function close to water saturation (θθs), and also the omission of the saturated part of sorptivity by regular approaches when h1>ha. We considered three commonly used hydraulic models, for which Lassabatere et al. (2021) proposed analytical formulations for S(-,0): Brooks and Corey (BC), van Genuchten–Burdine (vGB), and van Genuchte–Mualem (vGM). The second part of the paper compares the accuracy of the mixed formulation with the current strategies for the same three hydraulic models plus the Kosugi (KG) model and demonstrates the risk of serious misestimations with prior approaches. By presenting a new formulation that is applicable to any type of condition, this paper completes the study of Lassabatere et al. (2021), who proposed a scaling procedure for the approximation of SKh0,h1=0 with the condition of null water pressure head at the surface, i.e., h1=0. Lastly, we show how using the proposed approach can improve the accuracy of sorptivity estimation and, consequently, the modeling of water infiltration into soils with analytical models that make use of sorptivity.

2 Theory

2.1 Proposed new mixed formulation for computing sorptivity

To build the mixed formulation, SM, we start with the conductivity form of sorptivity, SK, since it includes both unsaturated and saturated parts. Then, we define an intermediate water pressure head between the initial and final water pressure heads, hch0,h1, smaller than the air-entry pressure, hc<ha0, and we split the integral into two separate parts as follows:

(7) S M h 0 , h 1 = S K h 0 , h 1 = h 0 h 1 θ h 1 + θ ( h ) - 2 θ h 0 K ( h ) d h = h 0 h c θ h 1 + θ ( h ) - 2 θ h 0 K ( h ) d h + h c h 1 θ h 1 + θ ( h ) - 2 θ h 0 K ( h ) d h = θ h 0 θ h c θ h 1 + θ - 2 θ h 0 D ( θ ) d θ + h c h 1 θ h 1 + θ ( h ) - 2 θ h 0 K ( h ) d h .

In Eq. (7), the integral h0hcθh1+θ(h)-2θh0K(h)dh is transformed into θh0θhcθh1+θ-2θh0D(θ)dθ thanks to the change of variable hθ. This operation requires that the function θ(h) is bijective over the whole interval [h0,hc], which is valid so long as hc<ha. Alternatively, the mixed formulation, SM, may be written alternatively as follows:

(8) S M h 0 , h 1 = θ 0 θ c θ 1 + θ - 2 θ 0 D ( θ ) d θ A + h c h 1 θ 1 + θ ( h ) - 2 θ 0 K ( h ) d h B ,

where θc=θ(hc), θ1=θ(h1), and θ0=θ(h0). The constraint hc<ha ensures that the computation of A in Eq. (8) avoids the challenging integration of infinite diffusivity close to saturation, Dθs=+, since θc<θs. In addition, hc is bounded to a finite value to avoid integration over infinite intervals for part B. Then, the two integrals involved in Eq. (8), A and B, only involve bounded functions over finite intervals, ensuring an easy numerical computation. An illustration of the procedure is depicted in Fig. 1.

https://hess.copernicus.org/articles/27/895/2023/hess-27-895-2023-f01

Figure 1Concept of the mixed formulation, SM(h0,h1): the integration of h0h1(θ1+θ(h)-2θ0)K(h)dh (=SM2(h0,h1)) is converted into the sum of the integration of two bounded functions over bounded intervals, θ0θc(θ1+θ-2θ0)D(θ)dθ and hch1(θ1+θ(h)-2θ0)K(h)dh. Note that the data are depicted on a log scale for clarity, but the integration is performed directly on the integrands instead of their log-scaled counterparts. Consequently, the integrals do not correspond to the areas below the curves, conversely to the case of linear scales. An illustrative case of the computation of sorptivity, SM(h0,h1), for a synthetic loamy soil with an initial water pressure head of −103 mm and a final water pressure head of -10-1 mm.

Download

Next, we scale sorptivity to separate the respective contributions of scale and shape parameters, as suggested by Lassabatere et al. (2021). We consider the following scaling relationships for hydraulic variables and sorptivity:

(9) S e = θ - θ r θ s - θ r h * = h | h g | K r = K K s ,

where Se is the saturation degree, h* is the scaled water pressure head, Kr is the relative hydraulic conductivity, θr and θs are the residual and the saturated soil water contents, hg is the scale parameter for the water pressure head, and Ks is the saturated hydraulic conductivity. The application of scaling relationships of Eq. (9) to the dimensional sorptivity expressions leads to the following equation (Ross et al.1996):

(10) S = | h g | K s θ s - θ r S * ,

where S and S* are respectively the dimensional and the scaled sorptivities. The application of the scaling equations Eqs. (9) and (10) to the mixed formulation SM, defined by Eq. (7), leads to the final expression proposed in our study:

(11) S M h 0 , h 1 = | h g | K s θ s - θ r S M * h 0 * , h 1 * S M * h 0 * , h 1 * = S e h 0 * S e h c * S e h 1 * + S e - 2 S e h 0 * D * S e d S e + h c * h 1 * S e h 1 * + S e h * - 2 S e h 0 * K r h * d h * 1 / 2 ,

where SM*h0*,h1* is the scaled version of the proposed mixed formulation, SM(h0,h1), with h0*=h0/|hg| and h1*=h1/|hg|. Equation (11) can be demonstrated by changing the integration variable θSe in the first and hh* in the second integral of Eq. (7). Since Eq. (11) cannot be evaluated with coding software that does not allow infinite values, h0*=-, we replace the input h0 (that may be infinite) with Se,0=Seh0*, which always remains bounded in the final expression of the mixed formulation SM*:

(12) S M * S e , 0 , h 1 * = S e , 0 S e h c * S e h 1 * + S e - 2 S e , 0 D * S e d S e + h c * h 1 * S e h 1 * + S e h * - 2 S e , 0 K r h * d h * .

Several options exist for the choice of the intermediate water pressure head hc* and intermediate saturation degree Se,c=Sehc*. In this study, our preferred option is to set the intermediate saturation degree as the average between the initial and the final saturation degrees, Se,c=Se,0+Se,1/2. However, under certain circumstances (e.g., for soils with gradual water retention functions; see Results section), the value of h*(Se,c) may reach very large values, leading to numerical instabilities. Therefore, we use the following criteria to ensure that h*(Se,c) remains finite:

(13) h c * = - min h * S e , 0 + S e , 1 2 , 10 z z Z S e , c = S e h c * .

When necessary, the value of z is varied until the two integrals in SM* (Eq. 12) converge. In most cases, z{-2,-1,0,1,2} ensures convergence regardless of soil type and situation. Note that for hydraulic models with non-null water entry pressure head, ha<0, z should be fixed with z≥0 so as to ensure -10z-1 and thus hc*ha*=-1. This condition is necessary to ensure the bijectivity of the function Se(h*) over the interval h0*,hc*, which is required for the use of Eq. (12).

In the following section, the mixed formulation SM* (Eqs. 12 and 13) is compared to several strategies previously proposed in the literature to cope with situations of numerical indeterminacy, e.g., at saturation θ1=θs (or, Se,1=1) for a null water pressure head at the surface, h1=0, and for very dry initial conditions θ0θr (or, h0*-).

2.2 Usual methods for computing sorptivity based on SD and SK

2.2.1 Computing sorptivity with SK for very dry initial conditions, h0-

Regarding the computation of sorptivity for very dry initial conditions with SK, one of the strategies found in the literature applies the regular definitions of sorptivity (Parlange1975), Eq. (5), to the case of very low values of h0. Such an approach was used by Di Prima et al. (2020) to estimate SK(h0,h1) for very dry soils. In this case, the maximum sorptivity, S-,h1, is approached as follows:

(14) S - , h 1 = - h 1 θ 1 + θ ( h ) - 2 θ r K ( h ) d h = lim h 0 - h 0 ¯ h 1 θ 1 + θ ( h ) - 2 θ h 0 ¯ K ( h ) d h = lim h 0 - S K h 0 ¯ , h 1 .
https://hess.copernicus.org/articles/27/895/2023/hess-27-895-2023-f02

Figure 2Illustration of the regular strategies for the estimation of the limits for the case of very dry conditions with the estimation of S(-,h1) using either SK(h0,h1), Eq. (5), (a) or SK-V2(h0,h1), Eq. (16), (b). The integration proceeds from the given value of h0 to the right, hh0, corresponding to the solid zones, and the value of h0 is lowered step by step to reach the limit S(-,h1) (see the arrow and the extension of the solid zones). In the equations, the variables that are adjusted to reach the targeted limits are underlined, and figures are zoomed in on in the vicinity of the limits. The differences in equations between the target and the integrated integrands are in red. An illustrative case of the computation of the limits of SK(h0,h1) or SK-V2(h0,h1) for a synthetic loamy soil with a final water pressure head of h1=-1 mm; successive values are considered for the initial water pressure heads, h0{-101,-102,-103} mm.

Download

Note that, for the sake of clarity, the underline in Eq. (14) shows the variables that are adjusted until the solution converges. Note also that the preceding equations are valid since limf(x)=limf(x), with xx being a continuous function. In practical applications of this method, SK(h0,h1) is computed for decreasing values of h0 until stabilization is reached (Fig. 2a). The last value obtained in this way is considered to represent the sorptivity at extremely dry conditions, i.e., S-,h1. This option is quite easy to implement since it requires the user to only code the regular function SK before applying it to very negative values of h0.

We propose an alternative procedure that employs a specific integrand to compute the same limit (Fig. 2b). In this case, the water content θ0 is set equal to θr in the integrand so as to correspond to the targeted initial conditions θh0=-=θr:

(15) S - , h 1 = - h 1 θ 1 + θ ( h ) - 2 θ r K ( h ) d h = lim h 0 - h 0 ¯ h 1 θ 1 + θ ( h ) - 2 θ r K ( h ) d h = lim h 0 - S K - V 2 h 0 ¯ , h 1

with the specific function SK-V2 defined as follows:

(16) S K - V 2 h 0 , h 1 = h 0 h 1 θ 1 + θ ( h ) - 2 θ r K ( h ) d h .

In comparison to SK, the water content θ0 is replaced with θr in the integrand for SK-V2. We expect this modification to improve numerical convergence towards the lower integration limit since SK-V2 directly integrates the right integrand (Fig. 2b). Conversely, SK integrates a distinct integrand, i.e., θ1+θ(h)-2θh0K(h)θ1+θ(h)-2θrK(h), thus involving an additional source of error (Fig. 2a). Briefly, SK combines the error due to the integral bound h0>- and the difference between its integrand and the targeted integrand. Note that the function SK-V2 should be restricted to the evaluation of S-,h1 and never used for the computation of other cases, i.e., Sh0-,h1, since the water content in the integrand is fixed at θr, which corresponds exclusively to the case of h0=-.

The application of the scaling, Eqs. (9) and (10), to the preceding definitions, Eqs. (14) and (15), leads to scaled versions of those equations, which will be used in the computations below:

(17) S * - , h 1 * = lim h 0 * - h 0 * ¯ h 1 * S e , 1 + S e h * - 2 S e h 0 * ¯ K r h * d h * = lim h 0 * - S K * h 0 * ¯ , h 1 *

(18) S * - , h 1 * = lim h 0 * - h 0 * ¯ h 1 * S e , 1 + S e h * K r h * d h * = lim h 0 * - S K - V 2 * h 0 * ¯ , h 1 *

with SK* and SK-V2*, the scaled versions of SK and SK-V2, defined as follows:

(19) S K * h 0 * , h 1 * = h 0 * h 1 * S e , 1 + S e h * - 2 S e h 0 * K r h * d h * S K - V 2 * h 0 * , h 1 * = h 0 * h 1 * S e , 1 + S e h * K r h * d h * .

The derivation of these equations involved scaling Eqs. (9) and (10) along with the change of variable θSe.

2.2.2 Computing sorptivity with SD for null water pressure head at surface, h1=0

A similar approach is often used with the SD formulation to avoid numerical indeterminacy close to saturation. The first option considers SD(θ0,θ1) with θ1θs as suggested, for instance, by Fernández-Gálvez et al. (2019):

(20) S u θ 0 , θ s = θ 0 θ s θ s + θ - 2 θ 0 D ( θ ) d θ = lim θ 1 θ s θ 0 θ 1 ¯ θ 1 ¯ + θ - 2 θ 0 D ( θ ) d θ = lim θ 1 θ s S D θ 0 , θ 1 ¯ .

Note that with this method, we can only account for the unsaturated part of sorptivity, Suθ0,θs=θ0θsθs+θ-2θ0D(θ)dθ, and we systematically miss the saturated portion of sorptivity 2θs-θ0Ksh1-ha, as mentioned in Sect. 1. The total sorptivity corresponds to the sum of its two components (see Eq. 6): Sh0,h1>ha=Su2θ0,θs+2θs-θ0Ksh1-ha. The subscript “u” in Su(θ0,θs) stands for “unsaturated” and serves as a reminder of that limitation (Ross et al.1996). This point will be further illustrated and discussed in Sect. 3.

The integrand specified by SD corresponds to θ1+θ(h)-2θ0D(θ). Consequently, SD combines the error due to the discrepancy between the integrated and the targeted integrands with the error resulting from the restriction of the integration to [θ0,θ1] instead of [θ0,θs] (Fig. 3a, SD). To correct this problem, we define a different estimator, SD-V2, to integrate directly the targeted integrand, θs+θ(h)-2θ0 D(θ) (Fig. 3b, SD-V2), and the following developments emerge:

(21) S u θ 0 , θ s = θ 0 θ s θ s + θ - 2 θ 0 D ( θ ) d θ = lim θ 1 θ s θ 0 θ 1 ¯ θ s + θ - 2 θ 0 D ( θ ) d θ = lim θ 1 θ s S D - V 2 θ 0 , θ 1 ¯

with the function SD-V2 defined as follows:

(22) S D - V 2 θ 0 , θ 1 = θ 0 θ 1 θ s + θ - 2 θ 0 D ( θ ) d θ .
https://hess.copernicus.org/articles/27/895/2023/hess-27-895-2023-f03

Figure 3Illustration of the regular strategies for estimating the limits for the case of saturation θ1θs, with the estimation of Su(θ0,θs) using either SD(θ0,θ1), Eq. (4), (a) or SD-V2(θ0,θ1), Eq. (22), (b). The integration proceeds from the given value of θ1 to the left, θθ1, defining the solid zones. The values of θ1 are increased step by step to reach Su(θ0,θs) (see the arrow and the extension of the solid zones). In the equations, the variables that are varied to reach the targeted limits are underlined, and figures are zoomed in on in the vicinity of the limits. The differences in equations between the target and the integrated integrands are in red. An illustrative case of the computation of the limits of SD(θ0,θ1) or SD-V2(θ0,θ1) for a synthetic loamy soil with an initial water content of θ0=0.25; successive values are considered for the final water contents, θ1{0.39,0.40,0.41}.

Download

As mentioned above, for SK-V2, SD-V2 should be only used for the determination of Su(θ0,θs) and not for the computation of sorptivity corresponding to other values of final water contents, since SD-V2 integrates the right integrand if and only if θ1=θs. The scaled version of these equations can be easily found by applying the scaling equations, Eqs. (9) and (10), to the previous equations, Eqs. (20) and (21), leading to their scaled versions:

(23)Su*Se,0,1=limSe,11Se,0Se,1¯Se,1¯+Se-2Se,0D*SedSe=limSe,11SD*Se,0,Se,1¯(24)Su*Se,0,1=limSe,11Se,0Se,1¯1+Se-2Se,0D*SedSe=limSe,11SD-V2*Se,0,Se,1¯.

The scaled versions of the formulation SD and SD-V2, SD* and SD-V2*, can be written as

(25) S D * S e , 0 , S e , 1 = S e , 0 S e , 1 S e , 1 + S e - 2 S e , 0 D * S e d S e S D - V 2 * S e , 0 , S e , 1 = S e , 0 S e , 1 1 + S e - 2 S e , 0 D * S e d S e .

In the following sections, we compare these previously used strategies, based on the use of SK*, and SD*, and the improved versions designed for the purpose of this study, SK-V2* and SD-V2* with the proposed mixed formulation SM*, by quantifying the accuracy and efficiency of each.

2.3 Validation of estimates against the nominal sorptivity for the selected hydraulic models

2.3.1 Hydraulic models and nominal sorptivity

The validation of the computation of sorptivity with the proposed mixed formulation SM (Eq. 11) and the usual strategies (see Sect. 2.2) was performed for hydraulic models that are commonly used for the hydraulic characterization of soils and also present different challenging features:

  • The Brooks and Corey (BC) model (Brooks and Corey1964) was among the first hydraulic models of soil physics (Hillel1998). It uses power law relationships to define the water retention (WR) and hydraulic conductivity (HC) functions and is often considered for integrating sorptivity and finding analytical solutions for water infiltration into soils (e.g., Varado et al.2006). The BC model reads as follows:

    (26) θ BC ( h ) = θ s h h BC θ r + θ s - θ r h BC h λ BC h < h BC K BC ( θ ) = K s θ - θ r θ s - θ r η BC .
  • The van Genuchten–Burdine (vGB) model combines the van Genuchten (1980) model with the Burdine condition m=1-2n for the WR function and the Brooks and Corey (1964) model for the HC function. It was the basis of the development of BEST methods and was often considered for the hydraulic characterization of soils (Lassabatere et al.2006; Yilmaz et al.2010; Bagarello et al.2014a). These formulations are considered to be one of the most consistent to use for modeling water infiltration into soils (Fuentes et al.1992). The vGB model reads as follows:

    (27) θ vGB ( h ) = θ r + θ s - θ r 1 + h h vGB n vGB - m vGB m vGB = 1 - 2 n vGB K vGB ( θ ) = K s θ - θ r θ s - θ r η vGB .
  • The van Genuchten–Mualem (vGM) model combines the van Genuchten (1980) model with Mualem's condition m=1-1n for the WR function and the Mualem (1976) capillary model for the HC function. The vGM model is among the most widely used models, in particular for the numerical modeling of flow in the vadose zone (Šimůnek et al.2003). The vGM model reads as follows:

    (28) θ vGM ( h ) = θ r + θ s - θ r 1 + h h vGM n vGM - m vGM m vGM = 1 - 1 n vGM K vGM ( θ ) = K s θ - θ r θ s - θ r l vGM 1 - 1 - θ - θ r θ s - θ r 1 m vGM m vGM 2 .
  • The Kosugi (KG) model relates the WR function to the soil pore size distribution assuming lognormal distributions (Kosugi1996). It is quite popular as the consequence of its physical meaning and soundness (Pollacco et al.2013; Nasta et al.2013) and was also recently implemented into BEST methods for the hydraulic characterization of soils (Fernández-Gálvez et al.2019). The KG model reads as follows:

    (29) θ KG ( h ) = θ r + θ s - θ r 2 erfc ln h h KG 2 σ KG K KG ( θ ) = K s θ - θ r θ s - θ r l KG 1 2 erfc erfc - 1 2 θ - θ r θ s - θ r + σ KG 2 2 ,

    where erfc stands for the complementary error function.

These models involve the following common scale hydraulic parameters: residual water content, θr; saturated water content, θs; scale parameter for the water pressure head, hg (hBC, hvGB, hvGM, or hKG); and saturated hydraulic conductivity, Ks. The BC models involve a non-null air-entry water pressure head, hBC, meaning that air needs a given suction to enter into the soil and to desaturate the soil. For the sake of simplicity, the scale parameter for the water pressure head is often fixed at the air-entry pressure head, so that hg=hBC. In addition, these hydraulic models involve one or two shape parameters for each set of WR and HC functions, specifically λBC and ηBC for the BC model; mvGB, nvGB, and ηvGB for the vGB model; mvGM, nvGM, and lvGM for the vGM model; and, lastly, σKG and lKG for the KG model. In order to simplify these equations and to reduce the risk of equifinality and non-unique optimization when inverting (Pollacco et al.2013), the following capillary model has been proposed to link these shape parameters (Haverkamp et al.2005):

(30) η = 2 λ + 2 + p ,

where λ=λBC for the BC model and λ=mn for the vGB models. In addition, the shape parameters lvGM and lKG are usually fixed at 12. In this study, the computations are performed considering the relationship given by Eq. (30).

The application of the scaling procedure Eq. (9) to these hydraulic models, i.e., Eqs. (26)–(29), leads to the following scaled hydraulic models (Lassabatere et al.2021).

  • BC model:

    (31) S e , BC h * = 1 - H 1 + h * | h * | - λ BC + H 1 + h * K r , BC S e = S e η BC .
  • vGB model:

    (32) S e , vGB h * = 1 + | h * | n vGB - m vGB with m vGB = 1 - 2 n vGB K r , vGB S e = S e η vGB .
  • vGM model:

    (33) S e , vGM h * = 1 + | h * | n vGM - m vGM with m vGM = 1 - 1 n vGM K r , vGM S e = S e l vGM 1 - 1 - S e 1 m vGM m vGM 2 .
  • KG model:

    (34) S e , KG h * = 1 2 erfc ln | h * | 2 σ KG K r , KG S e = S e l KG 1 2 erfc erfc - 1 2 S e + σ KG 2 2 .

These hydraulic models have the following hydraulic diffusivity functions, D*Se=KSedh*dSe (Lassabatere et al.2021):

(35) D BC * S e = 1 λ BC S e η BC - 1 λ BC + 1 D vGB * S e = 1 - m vGB 2 m vGB S e η vGB - 1 + m vGB 2 m vGB 1 - S e 1 m vGB - 1 + m vGB 2 D vGM * S e = 1 - m vGM m vGM S e l vGM - 1 m vGM 1 - S e 1 m vGM - m vGM + 1 - S e 1 m vGM m vGM - 2 D KG * S e = 1 2 π 2 σ KG S e l KG erfc erfc - 1 2 S e + σ KG 2 2 e erfc - 1 2 S e 2 + 2 σ KG erfc - 1 2 S e .

These equations are needed for the computation of the dimensionless sorptivity with the proposed mixed formulation SM*2 and the regular formulations SK*2, SK-V2*2, SD*2, and SD-V2*2.

The studied hydraulic models exhibit contrasting and challenging features for the computation of sorptivity, including non-null water pressure heads ha*<0 and infinite hydraulic diffusivity close to saturation limSe1D*Se=+. The complexity may also increase with the values of related shape parameters. Therefore, Lassabatere et al. (2021) defined a shape index, x, to characterize the spread of the WR functions around Seh*=12. Regardless of the chosen hydraulic model, the values of x close to zero correspond to a large spread of WR functions with a smooth descent of the saturation degree, Se, with the increase of |h*| (see Fig. 2, in Lassabatere et al.2021, and also Sect. 3). Conversely, when x gets close to unity, WR functions approach the stepwise function with an abrupt decrease of Se with the increase of |h*|. Lassabatere et al. (2021) defined the WR shape index x as follows:

(36) x BC = λ BC 2 + λ BC x vGB = m vGB x vGM = m vGM x KG = 1 1 + σ KG .

Lassabatere et al. (2021) also analytically determined the maximum squared scaled sorptivity S*2(-,0), also referred to as the parameter cp, as a function of the WR shape index x, for the BC, vGB, and vGM models:

(37) c p , BC ( x ) = 2 + 1 - x 5 x + 1 + 1 - x 7 x + 1 c p , vGB ( x ) = Γ 3 - x 2 Γ 1 + 5 x 2 Γ 1 + 2 x + Γ 1 + 7 x 2 Γ 1 + 3 x c p , vGM ( x ) = Γ ( 2 - x ) Γ 3 2 x 3 2 x - 1 Γ 1 2 x + Γ 5 2 x 5 2 x - 1 Γ 3 2 x + 1 - x Γ 3 2 x Γ ( 1 + x ) 3 2 x - 1 Γ 5 2 x + Γ 5 2 x Γ ( 1 + x ) 5 2 x - 1 Γ 7 2 x - 2 1 3 2 x - 1 + 1 5 2 x - 1 .

Note that no analytical formulation can be found for the case of Kosugi's hydraulic model, so cp,KG(x) must be computed numerically (Lassabatere et al.2021). Note also that in Eq. (37), nominal sorptivities are defined with the use of the capillary model Eq. (30) for the BC and vGB models, and lvGM=lKG=12.

2.3.2 Paper methodology and computations

In this study, we stress the following points: (i) the studied hydraulic models for WR and HC functions exhibit challenging conditions for the computation of sorptivity; (ii) the proposed mixed formulation is an ideal estimator for sorptivity; (iii) the usual methods, based on the use of SK and SD (Eqs. 5 and 4), or their improved versions, SK-V2 and SD-V2 (Eqs. 16 and 22), do not necessarily provide accurate estimations of the targeted nominal sorptivity; and (iv) errors in sorptivity estimation may drastically impact its use for modeling water infiltration into soils. To demonstrate these points, we consider the following conditions. Firstly, we only investigate the case of scaled sorptivity. For any estimator any estimator S^ of the nominal dimensional sorptivity S, and the related scaled variables, S^* and S*, the following relations emerge between the relative errors, Er, related to dimensional and scaled sorptivity:

(38) E r ( S ) = S ^ - S S = S ^ * | h g | K s θ s - θ r - S * | h g | K s θ s - θ r S * | h g | K s θ s - θ r = S ^ * - S * S * = E r S * ,

proving that the accuracy of the scaled estimator corresponds exactly to the accuracy of the dimensional estimator. Next, we consider the maximum scaled sorptivity, S*(-,0), since it involves at the same time the two types of challenges, i.e., very dry initial conditions with infinite water pressure head, h0=-, and saturated final conditions with null water pressure head, h1=0. Then, the modeling of water infiltration into soils is illustrated for a synthetic loamy soil and relies on the use of dimensional sorptivity (Sect. 3.4).

The first step, point (i), involves the study of the selected models with regard to the shapes of WR and HC functions. We computed the WR and HC functions for the four selected models, considering the following values of the WR shape index: x{0.01, 0.02, …, 0.99}. For the second goal of the study, point (ii), we compared the values provided by the proposed mixed formulation SM* with the nominal (error-free) values of sorptivity, S*=S*(-,0)=cp, provided by the exact analytical formulations (Eq. 37) for BC, vGB, and vGM models. The computations were performed for all the values of the WR shape index x, and the accuracy of SM* is discussed as a function of x.

For the third goal, point (iii), the estimations provided by the usual strategies were compared to the estimates provided by the proposed mixed formulation, SM*, to evaluate the efficiency of those previously used strategies. We considered several scenarios for the use of SK*, SK-V2*, SD*, and SD-V2*, with several values of the lower water pressure head, h0*, used in Eqs. (17) and (18) and several values of the final saturation degree, Se,1, used in Eqs. (23) and (24). We also investigated the estimation error as a function of the shape parameters for the four studied hydraulic models (BC, vGB, vGM, and KG).

Finally, regarding point (iv), we computed the dimensional sorptivity for a synthetic loamy soil, S(h0,h1), using both the mixed formulation and the regular methods, to investigate its dependency on both initial and final conditions and hydraulic models (BC, vGM, and KG). For the sake of clarity, we omitted the case of the vGB model and considered the other three models since their behaviors provided the greatest contrasts. Then, we investigated its use for the modeling of water infiltration into loamy soil. For that final purpose, we applied the different calculated sorptivity values to analytical models for the modeling of cumulative infiltrations corresponding to a set of water pressure heads imposed at the soil surface (hf=-150 or hf=30 mm) along with an initial water pressure head of −10 m (hi=-104 mm). The analytical models used to model the cumulative infiltrations are described in Sect. 3. The BC, vGM, and KG models were fitted to the hydraulic functions of the loamy soil, as defined by Carsel and Parrish (1988), before computation of sorptivity to investigate the impact of the choice of the hydraulic models on sorptivity estimation and related impact on cumulative infiltrations. The different synthetic cumulative infiltrations were discussed with regard to the sorptivity estimation method, the choice of hydraulic models for the WR and HC functions, and the choice of the analytical model. All the computations were performed using Scilab software (Campbell et al.2010) and are available on Zenodo (Lassabatere2022).

3 Results

3.1 Analysis of the selected hydraulic models and related challenging features

This section presents the four sets of models for describing the water retention and hydraulic conductivity functions, their features, and their dependency on related shape parameters. The features of the WR-HC functions are determinant with regard to the estimation of sorptivity. The BC model has a non-null air-entry water pressure head (Fig. 4a, shown by the plateau for |h*|1, i.e., h*-1), meaning that the sorptivity has a non-null saturated part that must be accounted for. This condition is one of the challenging features of the diffusivity form of sorptivity, Eq. (4). Conversely, the three other hydraulic models do not have any air-entry water pressure head values (Fig. 4e, i, and m do not have plateaus) but rather have infinite values of the hydraulic diffusivity close to saturation, thus posing potential problems of convergence (Fig. 4h, l and p). The use of these models allows us to characterize the improvements offered by SM* compared to the usual use of SD* for the cases of problematic computation close to saturation (i.e., h1→0 and θ1θs). Similarly, the accuracy of the commonly used SK* version can be challenged when integrating over infinite intervals (−∞, 0], particularly for hydraulic models that are characterized by a slight decrease in the saturation degree, Se, for quasi-infinite water pressure heads, h*. The chosen hydraulic models are thus expected to be challenging, in particular the BC, vGB, and vGM models, which keep high values of saturation degrees even for quasi-infinite water pressure heads (Fig. 4a, e, and i). The KG model, which is more symmetrical and characterized by a larger decrease in Se when h* decreases (Fig. 4m), is expected to be less challenging. Regarding those challenges, the shapes of the WR, HC, and hydraulic diffusivity functions are also of importance. We expect the small values of x to be more problematic, in particular with the use of SK*, due to the smooth decrease in Se with h* (Fig. 4, first column). Conversely, we expect more problems with the use of SD* for large values of x, with quasi-infinite values for the hydraulic diffusivity close to saturation, i.e., Se→1 (Fig. 4h, l, and p).

https://hess.copernicus.org/articles/27/895/2023/hess-27-895-2023-f04

Figure 4Water retention (WR) and hydraulic conductivity (HC) curves for different values of the WR shape index x. Panels (a, e, i, m) show WR as Se(h*), panels (b, f, j, n) show HC as Kr(Se), panels (c, g, k, o) show HC as Kr(h*), and panels (d, h, l, p) show diffusivity as D*(Se); the four tested models include Brooks and Corey (BC) (a–d), van Genuchten–Burdine (vGB) (e–h), van Genuchten–Mualem (vGM) (i–l), and Kosugi (KG) models (m–p). The arrows indicate the trends with increasing WR shape index x. The hydraulic parameters λBC, mvGM, mvGB, and σKG were computed as a function of x using Eq. (36) with lvGM=lKG=12. Adapted from Lassabatere et al. (2021).

3.2 Validation of the proposed mixed formulation, SM*, against the nominal sorptivity, S*(-,0)

The computations using SM* (Eqs. 12 and 13) with hc*=h*Se,0+Se,12 were efficient in most cases, regardless of the value of the WR shape index x. The use of the threshold 10z was necessary for the first value of the WR shape index x for the BC model (xBC=0.01), the first two values for the vGB model (xvGB{0.01,0.02}), and the first 17 values for the vGM model (xvGM{0.01,0.02, …,,0.17}). The value of z=0 was enough to allow the computation in all cases, apart from the case of the vGM model for which the value of z=-1 had to be considered for xvGM{0.02, …, 0.07} and z=-2 for xvGM=0.01. Note that, as long as convergence is obtained, the values of SM*(-,0) do not depend on z. With this strategy involving a threshold, the proposed mixed formulation, SM*, provides estimates for all cases, i.e., for all hydraulic models and all values of the WR shape index, x. A sensitivity analysis was also performed for the KG model and led to the same success with z=0 regardless of the value of the WR shape index xKG (data not shown).

The relative error, Er, between the estimates provided by the proposed mixed formulation, SM*(-,0), and the targeted sorptivity, S*(-,0)=cp, was analyzed in terms of means, standard deviations, and minimum and maximum values (Table 1):

(39) E r = S M * ( - , 0 ) - c p c p .

Table 1Absolute values of relative errors, |Er|, between the proposed mixed formulation, SM*, and the targeted scaled sorptivity, S*(-,0)=cp, with the mean |Er|, the standard deviation (σ|Er|), the minimum and the maximum values for the three hydraulic models whose sorptivity is analytically tractable: Brooks and Corey (BC), van Genuchten–Burdine (vGB), van Genuchten–Mualem (vGM). Note that 10−16 corresponds to the relative precision of numbers in Scilab.

Download Print Version | Download XLSX

The accuracy of the proposed mixed formulation, SM*, Eqs. (12) and (13), is excellent for all the models and all values of the WR shape index x (Table 1). The average relative errors are on the order of 10−13 for the BC and vGB models and on the order of 10−9 for the vGM model (Table 1Er). The minimum errors were <10-15 for all the models (Table 1, Min). The maximum errors are 10-12 for the BC and the vGB models and 10-7 for the vGM model (Table 1, Max). In other words, the mixed formulation, SM*, provides extremely accurate estimations of the targeted scaled sorptivity S*(-,0). The proposed mixed formulation, SM*, can therefore be considered to be an excellent estimator of sorptivity in all cases, regardless of the choice of the hydraulic model and related values of shape parameters.

3.3 Study of the usual strategies for estimating the nominal sorptivity, S*(-,0)

In this section, we compare the estimates provided by the strategies commonly considered (i.e., SD*, SD-V2*, SK*, and SK-V2*) with the reference values of the targeted sorptivity S*=S*(-,0). For these comparisons, we consider the analytical formulations of Eq. (37) for the BC, vGB, and vGM models and use the proposed mixed formulation SM* for the KG model. Indeed, no analytical expressions are available for this last model, whereas the estimates provided by S*=SM* are very accurate for the BC, vGB, and vGM models (see Sect. 3.2), and thus SM* is assumed to be similarly accurate for the KG model.

3.3.1 Illustrative example

In this example, we investigate the accuracy of the limits of the functions SK* and SK-V2* towards S*. These functions, in particular, SK*, are often used without special attention regarding their accuracy. SK* and SK-V2* converge to the limit S* when |h0*| becomes large enough (Fig. 5, left column). For instance, for the BC model with λ=0.56, the use of SK* with h0*=-10 and h0*=-100 has respective relative errors of Er=-4.1 % and Er=-15 % (Fig. 5a, dashed red line). The second estimator, SK-V2*, converges much faster than SK*. With the same values of h0*, the relative errors drop below 0.01 % in absolute values (Fig. 5a, continuous blue line). In this case, the convergence towards the target can be achieved by integrating from h0*=-10, with high accuracy. Such an improvement results from setting the initial saturation degree, Se,0, at its target value, Se,0=0, in the integrand (see Eq. 18 versus Eq. 17). The same conclusions hold regardless of the hydraulic model (Fig. 5c, e, and g).

https://hess.copernicus.org/articles/27/895/2023/hess-27-895-2023-f05

Figure 5Relative errors of the regular functions SD*(0,Se,1) and SD-V2*(0,Se,1) (b, d, f, h) and SK*(h0*,0) and SK-V2*(h0*,0) (a, c, e, g) towards S*(-,0) for the four hydraulic models, for the following specific cases: BC model with λBC=56 (xBC=0.22), vGB model with nvGB=3 (xvGB=1/3), vGM model with nvGM=2 (xvGM=1/2) with lvGM=0.5, and KG model with σKG=1.5 (xKG=0.4) and lKG=0.5.

Download

The convergence of the two functions, SD* and SD-V2*, is depicted for the case of BC functions in Fig. 5b. For both estimators, the estimates are far from the target values, S*, with |Er| close to 50 %. This large error results from the omission of the saturated part of sorptivity, as explained above (Sect. 2.2). By design, SD* and SD-V2* converge towards the unsaturated part of sorptivity, Su(0,1)=011+SeD*SedSe, and miss the additional saturated part of the scaled sorptivity, 2Se,1-Se,0KrSe,1h1*-ha*=2, given that h1*=0 (Se,1=1), h0*=- (Se,0=0), and ha*=1. When the saturated portion is added to the estimators, the convergence becomes excellent (Fig. 5b, green line). These results show that the computation of sorptivity using the integration with regards to saturation degree, like Eqs. (4) and (25), leads to erroneous estimations. Accounting for the term 2(θs-θ0)Ks|ha| in Eq. (6) is thus essential. Note that Eq. (4) is often considered and used in most studies, even though it can lead to a large underestimation of sorptivity for the case of non-null air-entry water pressure heads. That factor is of prime importance regarding the accurate estimation of sorptivity.

For the other hydraulic models, the estimators SD* and SD-V2* are significantly better (Fig. 5d, f, and h). Indeed, the air-entry water pressure head is null in these models, which removes the large underestimation due to the omission of the saturated part of sorptivity, as observed for the BC model. However, the underestimation remains substantial, with relative errors on the order of 15 %–20 % for Se,1=0.99 (Fig. 5d, f and h). Even when Se,1=0.999, the relative errors are larger than 10 %. From these results, we conclude that determining the sorptivity, S*, by setting the final saturation degree to near-unity does not provide good estimates, regardless of the estimators considered. These poor estimates result from the fact that the diffusivity and, thus, the integrand are infinite. The convergence of the integrals defined by Eq. (25) and thus Eq. (4) is slow and prevents accurate estimations. Conversely, in the case of SK* and SK-V2*, changing the integrand by fixing it to the targeted integrand does not significantly improve the convergence.

3.3.2 Relative errors as a function of the WR shape index x

The previous results were presented for specific values of the shape parameters in Fig. 5. One may wonder if the accuracy of estimators SK* and SK-V2*, on the one hand, and SD* and SD-V2*, on the other, varies with the shape parameters and indices. Figure 6 depicts the relative errors of the estimators SK* and SK-V2* as a function of the shape parameters for different values of h0* (left column of Fig. 6) and of the estimators SD* and SD-V2* for different values of Se,1 (right column of Fig. 6). For all hydraulic models except the BC model, the best estimates for SD* and SD-V2* are obtained for intermediate values of WR shape indices (Fig. 6, right column). Discrepancies increase for both small or large values of the WR shape index. As detailed above (Sect. 3.3.1), accurate predictions require the practitioner to fix the upper integration boundary, Se,1, to a minimum of 0.999, to get errors less than 10 %. However, estimates remain poor, even with Se,1=0.9999 or Se,1=0.99999, when the WR shape indices tend towards unity (Fig. 6d, f, and h).

https://hess.copernicus.org/articles/27/895/2023/hess-27-895-2023-f06

Figure 6Relative errors of estimators SD* and SD-V2* (b, d, f, h) and SK* and SK-V2* (a, c, e, g) for the four selected hydraulic models, i.e., BC, vGB, vGM, and KG models, as a function of the WR shape indices.

Download

For the BC model, the estimators SD* and SD-V2* provide poor estimates under all circumstances, since they miss the saturated part of sorptivity, as explained above (Fig. 6b, SD*SD-V2*). Adding the saturated part of sorptivity substantially improves the computation, with very accurate estimations (Fig. 6b, 2+SD*2). Regarding the estimators SK* and SK-V2* (Fig. 6, left column), the estimates improve when the WR shape index converges towards 1. It can be noted that SK-V2* always provides quite accurate predictions with very small relative errors (Fig. 6a, c, e, and f).

The results obtained in this study also revealed particular behaviors for the different formulations used to compute sorptivity, specifically SK*, SK-V2*, SD*, and SD-V2*. The approaches that compute sorptivity by integration with regards to h* (i.e., SK* and SK-V2*) have errors that come from the choice of the lower integration boundary h0* and from the value of Se,0 in the integrand. Solutions that compute sorptivity by integration with regards to Se (i.e., SD* and SD-V2*) have more error associated with the omission of the saturated part, as well as slow convergence of the integration procedure when dealing with infinite diffusivity values. Among the usual procedures, the SK* estimator is better than SD*, and the SK-V2* estimate is superior to both. Indeed, SK-V2* does not miss the saturated part of sorptivity, like SK*, but provides much lower discrepancies than the other estimators regardless of the selected hydraulic model (Fig. 6). At the same time, the adaptation of the integrand, by fixing the values of the initial saturation degree, Se,0, to its final values, substantially improves the quality of the estimator. This formulation thus always provides estimates with |Er|<10 %, regardless of the selected hydraulic model and the value of the shape index. However, these errors are still much larger than those obtained with the mixed formulation, SM*. Consequently, the use of the mixed formulation, SM*, should be favored as far as possible.

3.4 Computation of sorptivity and accuracy of water infiltration models based on sorptivity

In the following, we illustrate the need to use accurate sorptivity estimates when modeling water infiltration into soils. We consider the case of one-dimensional (1D) water infiltration into a loamy soil as a function of initial conditions, i.e., hi, and final conditions, i.e., hf. In addition to the illustration of the practical use of the proposed mixed formulation, SM, we demonstrate that (i) SM provides better estimates for sorptivity than the usual estimator, SD, because it always includes the saturated part of sorptivity (when applicable); (ii) this improvement avoids significant misestimations of sorptivity and related errors for the modeling of the cumulative infiltration into soils; and (iii) the proper analytical model is needed to account for positive water pressure head at surface. Note that, for the sake of clarity, we compare SM to SD, only, with no comparison with the other methods since SD-V2 gives similar results to SD and since SK and SK-V2 do not miss the saturated part of sorptivity and involve much lower errors. Lastly, the impact of the selection of the hydraulic models for the description of the soil hydraulic functions on the computation of sorptivity and on related modeling of the cumulative infiltration is also discussed.

3.4.1 Illustrative example of computation of sorptivity and dependency upon initial and final conditions

The studied synthetic loamy soil was defined by Carsel and Parrish (1988) considering the vGM hydraulic model (Eq. 28) with the hydraulic parameters tabulated in Table 2 (column “vGM”). In addition, we used the BC and KG models (Eqs. 26 and 29), with the parameters tabulated in Table 2 (column “BC” and “KG”). Related hydraulic parameters were optimized by fitting BC and KG models to the vGM model and are available in the HYDRUS software suite (Radcliffe and Simunek2010). Figure 7 depicts the water retention curve (Fig. 7a) and unsaturated hydraulic conductivity (Fig. 7b) and demonstrates the proper alignment of BC and KG models on the vGM model, despite the slightly greater discrepancy for the BC model close to saturation for the hydraulic conductivity (Fig. 7b).

https://hess.copernicus.org/articles/27/895/2023/hess-27-895-2023-f07

Figure 7Water retention (a) and hydraulic conductivity (b) curves described by the three sets of hydraulic models (BC, vGM, and KG) and sorptivity as a function of initial (c) and final (d) water pressure heads with the sorptivity estimated by the proposed mixed formulation SM and the regular model SD.

Download

Table 2Hydraulic parameters, i.e., values of the residual and saturated water contents, θr and θs; the scale parameter for water pressure head hg; the air-entry water pressure head ha; the saturated hydraulic conductivity Ks; and the shape parameters involved in the hydraulic models used for the description of the water retention and the hydraulic conductivity curves, i.e., Eqs. (26), (28), and (29).

Download Print Version | Download XLSX

At first, we illustrate the use of the mixed formulation, SM, for the computation of sorptivity for a final water pressure head at the surface of hf=-150 mm and an initial water pressure head of hi=-104 mm =−10 m. This initial condition corresponds to the isostatic water pressure head in equilibrium with a water table positioned at 10 m below the ground. The value of −150 mm, considered at the surface, is often used for water infiltration with tension disk experiments to deactivate the macropores and infiltrate only into the matrix (e.g., Timlin et al.1994; Malone et al.2004; Lassabatere et al.2014). Practically, the computation of dimensional sorptivity with the proposed procedure SM can be performed using the codes developed for this paper and downloadable from Zenodo (Lassabatere2022). Alternatively, computations may be performed as follows:

  1. Compute the initial and final saturation degrees, Se,i=Se(hi) and Se,f=Se(hf).

  2. Compute the intermediate water pressure head, hc*=-minh*Se,i+Se,f2,10z with z{-2,-1,0,1,2}.

  3. Compute the intermediate saturation degree from the intermediate water pressure head, Se,c=Se(hc*).

  4. Integrate the lower part of the squared sorptivity

    A=Se,iSe,c(Se,f+Se-2Se,i)D*(Se)dSe.
  5. Integrate the upper part of the squared sorptivity

    B=hc*hf*(Se,f+Se(h*)-2Se,i)Kr(h*)dh*.
  6. Combine the two parts to compute the mixed formulation SM*=A+B.

  7. Upscale by multiplying the scaled sorptivity,

    SM=|hg|Ks(θs-θr)SM*.

For the studied case, with hi=-10 m and hf=-150 mm, it leads to the following results:

  1. The related scaled water pressure heads are hi*=-36 and hf*=-0.54, with related saturation degrees of Se,i=0.134 and Se,f=0.89.

  2. The intermediate water pressure head takes the value of hc*=-1, which corresponds to the maximum between h*Se,i+Se,f2=-2.96 and -100=-1.

  3. The corresponding saturation degree takes the value of Se,c*=0.780.

  4. The computation of the lower part of the squared sorptivity gives A=0.0344.

  5. The computation of the upper part of the squared sorptivity gives B=0.0515.

  6. The combination of the two parts gives the scaled sorptivity: SM*=0.2931.

  7. The upscaling finalizes the computation of sorptivity: SM=0.156 mm s-1/2.

Then, we computed sorptivity with the mixed formulation, SM, and the regular estimator, SD, as a function of the initial water pressure head, hi, for the water pressure head at the surface of hf=-150 mm (Fig. 7c) and as a function of final water pressure head, hf, for an initial water pressure head of hi=-104 mm =−10 m (Fig. 7d). We considered three different models (vGM, BC, and KG). We first analyzed the values of sorptivity as estimated by the mixed formulation SM. The values of sorptivity decrease with the initial water pressure head, hi (Fig. 7c), and increase with the final water pressure head, hf (Fig. 7d), regardless of the chosen hydraulic model, vGM, BC, or KG. This trend was expected since the integrand is always positive, and the integration of positive functions decreases with regard to its lower boundary and increases with regard to its upper boundary. Those trends have already been documented by many authors (e.g., Stewart et al.2013). The values of sorptivity were also computed with the regular estimator SD (Fig. 7c and d, dashed lines). In most cases, the two estimators, SM and SD, provided the same estimations, except when the final pressure head was higher than the air-entry water pressure head, hfha. In this case, as explained above, SD missed the saturated part of the sorptivity, which explains why the values predicted with SD no longer increase with hf but remain equal to the unsaturated part of sorptivity: SDθi,hfha=SDθi,θs=Suθi,θs (Fig. 7d, plateaus with dashed lines). Note for vGM and KG models, the air-entry water pressure head is zero, then the plateaus begin at hf=0, whereas the air-entry water pressure head is lower than zero for the BC model, ha=-111.5 mm, inducing an earlier beginning of the plateau (Fig. 7d, BC versus vGM and KG).

Another interesting point is the difference between the three models. The BC model exhibits much larger sorptivity values, followed by the vGM model and then the KG model. However, the three models are meant to represent the same water retention and hydraulic conductivity functions, i.e., the same soil. That result shows that the selection of the hydraulic model for the water retention and hydraulic conductivity functions substantially impacts the value of sorptivity, even when the models are fitted to the same data and represent the same soil. This point was already raised by Lassabatere et al. (2021). This feature is counterintuitive since the use of close hydraulic models is expected to provide close values of sorptivity. Such a discrepancy may result from the different mathematical behaviors of the three different hydraulic models close to saturation. A closer look at the hydraulic functions shows that the unsaturated hydraulic conductivity differs close to saturation with the following ranking – BC  vGM  KG (Fig. 7b) – which may explain the observed discrepancy between values of sorptivity. We suspect that the ranking of the values of sorptivity between models reflects the ordering between unsaturated hydraulic conductivity close to saturation. Stewart and Abou Najm (2018) also attributed higher sorptivity values with the use of the BC model to the extra area under the K(h) curve.

3.4.2 Impact of misestimations of sorptivity on cumulative infiltrations

Given the dependence of sorptivity on the choice of the estimator and the hydraulic model, we expected some implications for water infiltration. To demonstrate this, we modeled water infiltration into the studied loamy soil, with a 30 mm water pressure head at the surface (hf=30 mm) and an initial water pressure head of −10 m (hi=-104 mm). To model water infiltration, we used the analytical model developed for hfha by Haverkamp et al. (1994) with its extension to hf>ha developed by Haverkamp et al. (1990). The first model, referred to as the quasi-exact implicit (QEI) model by many authors (e.g. Fernández-Gálvez et al.2019), reads as follows:

(40) 2 Δ K 2 S 2 t = 1 1 - β 2 Δ K S 2 I 1 D - K i t - ln e 2 β Δ K S 2 I 1 D - K i t + β - 1 β ,

where ΔK stands for the difference between final and initial values of hydraulic conductivity ΔK=Kf-Ki, and β stands for an infiltration constant, usually fixed at 0.6. Varado et al. (2006) followed by Lassabatere et al. (2009) suggested introducing the following scaling procedure for time t and cumulative infiltrations I:

(41) γ I = S 2 2 Δ K I = γ I I * + K i t γ t = S 2 2 Δ K 2 t = γ t t * .

These scaling equations lead to the following relationship between the scaled and dimensional versions of the QEI model (Lassabatere et al.2009):

(42) I QEI ( t ) = γ I I QEI * t γ t + K i t

with the following equations for the scaled QEI model, IQEI*t*:

(43) t * = 1 1 - β I * - ln e β I * + β - 1 β .

Note that in both equations Eqs. (40) and (43), the cumulative infiltration is defined implicitly, since time is defined as a function of the cumulative infiltration, instead of defining cumulative infiltration as a function of time. On the same basis, Ross et al. (1996) scaled the model developed by Haverkamp et al. (1990) for the case of hf>ha, leading to

(44)I*=σq*-1+1-σβln1+βq*-1(45)t*=1-σβ(1-β)ln1+βq*-1+σq*-1-1-σβ1-βln1+1q*-1,

where q*=qγq is the scaled infiltration rate and γqK. In this case, there is no direct relationship between the time t* and the cumulative infiltration I* but instead two relationships, with the first relating time and the infiltration rate, t*=t*q*, defined by Eq. (45), and the second relating the cumulative infiltration and the infiltration rate, I*=I*q*, defined by Eq. (44). To compute the cumulative infiltration as a function of time, the infiltration rate q* must be retrieved as the root of the Eq. (45) before being inserted into Eq. (44) to define the scaled infiltration as a function of the scaled time, IQEI_ext*t*. In the following, this relationship is referred to as the extended version of the QEI model and denoted as the QEI_ext model. This set of scaled models, Eqs. (44) and (45), must be upscaled considering the same scaling factors, as defined in Eq. (41) with the additional parameter, σ, which quantifies the relative contribution of the saturated part of sorptivity to the total sorptivity:

(46) σ = 2 θ s - θ i K s h f - h a S 2 .

The QEI_ext model addresses all cases. When hfha, the saturated part of sorptivity is null, leading to σ=0, and the QEI_ext model reduces to the regular QEI model. When hf>ha, the regular QEI model is no longer valid, and the extended QEI_ext model should be used instead. In other words, the QEI_ext model is always valid and should be preferred to the QEI model in any case. In the following example, we consider that two different errors may arise: (i) the wrong choice of the estimator for the sorptivity and (ii) the wrong choice of the model.

https://hess.copernicus.org/articles/27/895/2023/hess-27-895-2023-f08

Figure 8Cumulative infiltrations as a function of initial (a) and final (b) water pressure heads, computed with the correct estimate for sorptivity (SM) and the right model and illustration of the dependency of cumulative infiltration on the model selection (QEI model versus extended QEI_ext model) and on the selected estimator for the sorptivity (SM versus SD) for (c) the case of unsaturated conditions imposed at surface hfha (BC model, hi=-10 m and hf=-150 mm) and (d) the case of saturated conditions, hfha (BC model, hi=-10 m and hf=30 mm).

Download

Before investigating the impact of errors, we evaluated the dependence of the cumulative infiltrations as a function of initial and final conditions. For this analysis, we use the mixed formulation for sorptivity with the extended QEI_ext model. The results are depicted in Fig. 8a and b. Cumulative infiltration decreases with the initial water pressure head hi (Fig. 8a, arrows downwards) and increases with the final water pressure head hf (Fig. 8b, arrows upwards), as a consequence of the variations of sorptivity as a function of hi and hf (Fig. 7c and d). It can be noted that the impact of the water pressure head at surface hf is much more important, with much wider curve beams (Fig. 8b versus a; see shaded areas). Again, these trends are in line with previous studies (e.g., Lassabatere et al.2009, 2014; Angulo-Jaramillo et al.2016). The driver of water infiltration into the soil is the hydraulic gradient and, thus, the difference between the final and the initial water pressure heads. The higher the final water pressure head and the lower the initial water pressure head, the higher the hydraulic gradient and, thus, the higher the cumulative infiltration. As for the sorptivity, the selection of a hydraulic model for the water retention and hydraulic conductivity functions drastically impacts the results, with the same ordering as for sorptivity, i.e., largest cumulative infiltration for BC, then vGM, and, lastly, KG models. As a consequence, the fit of specific hydraulic models on the same water retention and hydraulic conductivity functions may lead to contrasting results in terms of water infiltration. To the authors' knowledge, this is the first time that the influence of the choice of the hydraulic model on sorptivity estimation and water infiltration has been so clearly demonstrated.

As stated above, either a wrong estimation of sorptivity or a wrong choice of the analytical model (QEI instead of QEI_ext model) may lead to erroneous estimates of the cumulative infiltration. The four combinations of the two models and the two estimates of sorptivity were used to model water infiltration for the case of an initial water pressure head of hi=-104 mm and a final water pressure head of hf=-150 mm (Fig. 8c) and for the case of the same initial water pressure head, hi=-104 mm, but a positive surface water pressure head hf=30 mm (Fig. 8d). For the first case, there was no difference at all (Fig. 8c). Both SD and SM accurately quantified the sorptivity, and the two models, QEI and QEI_ext were similar. In contrast, for the second case, the difference was quite significant. The use of the more appropriate estimator, i.e., SM instead of SD, increased the amount of modeled cumulative infiltration, regardless of the selected model. In addition, the use of the correct model, i.e., the QEI_ext instead of the QEI model, substantially increased the cumulative infiltration, with a slightly lower increase in comparison to that brought by the choice of SM for the computation of sorptivity. These results demonstrate the importance of using the proper estimator for sorptivity in combination with the best possible infiltration model.

4 Conclusions

The proper calculation of sorptivity is crucial to model accurately water infiltration into soils. However, in some cases, e.g., when the initial state is very dry, or the final state corresponds to saturated conditions, the numerical computation of sorptivity using Eq. (4) or Eq. (5) from the hydraulic functions may be a source of numerical errors and difficulties. Indeed, the integration procedure typically used involves either an infinite boundary or unbounded integrands. Many previous studies had attempted to alleviate these problems by fixing arbitrarily finite limits for the integration interval. In this study, we investigated the accuracy of these approaches and demonstrated the potentially massive misestimation of sorptivity that is possible when these arbitrary corrections are used. To alleviate those problems, we proposed a mixed formulation that was validated against analytical expressions of sorptivity for specific hydraulic models. The proposed mixed formulation proved highly accurate for all hydraulic models and shape parameters tested, with negligible relative errors (<10-7). Conversely, the use of regular estimates for sorptivity leads to large underestimates in many instances, in particular when the final water pressure head exceeds the air-entry water pressure head.

This study demonstrates that, through the use of the new mixed formulation, it is possible to compute sorptivity easily and very accurately. The proposed formulation presents a very practical tool that may be applied to any type of hydraulic model and any value of initial and final water pressure heads and water contents. The proposed approach allows sorptivity to be computed in all cases, thus improving the modeling of water infiltration into soils and the estimation of soil hydraulic properties. In addition, we used the proposed formulation to investigate the sensitivity of sorptivity to initial and final water pressure heads and to the choice of the hydraulic model chosen to quantify the water retention and unsaturated hydraulic curves. This analysis clearly demonstrated that sorptivity increases with the final water pressure head and decreases with the initial water pressure head. We also showed that a proper estimate of sorptivity is crucial with regard to the modeling of water infiltration into soils and that the selection of the model for the hydraulic functions drastically impacts the computation of sorptivity and, consequently, the final amounts of cumulative infiltrations.

Code availability

All computations were done using Scilab free software. The scripts for the computation of all the results and figures presented in this paper, and in particular the script for the proposed mixed formulation, can be downloaded online: https://doi.org/10.5281/zenodo.7044940 (Lassabatere2022).

Data availability

No data sets were used in this article.

Author contributions

LL established the question, performed the analytical developments, computed the numerical results, and provided the first draft of the manuscript. PEP confirmed the analytical and numerical developments and wrote the first draft with LL. DY, BL, DMF, SDP, and MR verified parts of the numerical computations. JP and JFG helped with the use of the Kosugi model. RDS and MAN helped with the editing, the layout of the manuscript, and the presentation of the results. All the authors contributed to the editing of the manuscript.

Competing interests

The contact author has declared that none of the authors has any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

This work was performed within the INFILTRON project supported by the French National Research Agency (ANR-17-CE04-010).

Financial support

This research has been supported by the Agence Nationale de la Recherche (grant no. ANR-17-CE04-010).

Review statement

This paper was edited by Nunzio Romano and reviewed by two anonymous referees.

References

Angulo-Jaramillo, R., Bagarello, V., Iovino, M., and Lassabatere, L.: Infiltration measurements for soil hydraulic characterization, Infiltration Measurements for Soil Hydraulic Characterization, Springer, https://doi.org/10.1007/978-3-319-31788-5, 2016. a, b, c

Angulo-Jaramillo, R., Bagarello, V., Di Prima, S., Gosset, A., Iovino, M., and Lassabatere, L.: Beerkan Estimation of Soil Transfer parameters (BEST) across soils and scales, J. Hydrology, 576, 239–261, https://doi.org/10.1016/j.jhydrol.2019.06.007, 2019. a

Bagarello, V., Di Prima, S., and Iovino, M.: Comparing alternative algorithms to analyze the beerkan infiltration experiment, Soil Sci. Soc. Am. J., 78, 724–736, https://doi.org/10.2136/sssaj2013.06.0231, 2014a. a, b

Bagarello, V., Di Prima, S., Iovino, M., and Provenzano, G.: Estimating field-saturated soil hydraulic conductivity by a simplified Beerkan infiltration experiment, Hydrol. Process., 28, 1095–1103, https://doi.org/10.1002/hyp.9649, 2014b. a

Bouwer, H.: Unsaturated flow in ground-water hydraulics, J. Hydraul. Div., 90, 121–144, https://doi.org/10.1061/JYCEAJ.0001098, 1964. a

Brooks, R. and Corey, A.: Hydraulic Properties of Porous Media, Hydrology Papers, Colorado State University, 1964.  a, b

Campbell, S. L., Chancelier, J. P., and Nikoukhah, R.: Modeling and Simulation in SCILAB, Modeling and Simulation in Scilab/Scicos with ScicosLab 4.4, Springer, https://doi.org/10.1007/978-1-4419-5527-2, 2010. a

Carsel, R. F. and Parrish, R. S.: Developing joint probability distributions of soil water retention characteristics, Water Resour. Res., 24, 755–769, https://doi.org/10.1029/WR024i005p00755, 1988. a, b

Cook, F. and Minasny, B.: Sorptivity of Soils, in: Encyclopedia of Earth Sciences Series, Springer, 824–826, https://doi.org/10.1007/978-90-481-3585-1_161, 2011. a

Di Prima, S., Stewart, R., Castellini, M., Bagarello, V., Abou Najm, M., Pirastru, M., Giadrossich, F., Iovino, M., Angulo-Jaramillo, R., and Lassabatere, L.: Estimating the macroscopic capillary length from Beerkan infiltration experiments and its impact on saturated soil hydraulic conductivity predictions, J. Hydrol., 589, 125159, https://doi.org/10.1016/j.jhydrol.2020.125159, 2020. a, b

Fernández-Gálvez, J., Pollacco, J., Lassabatere, L., Angulo-Jaramillo, R., and Carrick, S.: A general Beerkan Estimation of Soil Transfer parameters method predicting hydraulic parameters of any unimodal water retention and hydraulic conductivity curves: Application to the Kosugi soil hydraulic model without using particle size distribution data, Adv. Water Resour., 129, 118–130, 2019. a, b, c

Fuentes, C., Haverkamp, R., and Parlange, J.-Y.: Parameter constraints on closed-form soil water relationships, J. Hydrol., 134, 117–142, 1992. a

Haverkamp, R., Parlange, J.-Y., Starr, J., Schmitz, G., and Fuentes, C.: Infiltration under ponded conditions: 3. A predictive equation based on physical parameters, Soil Sci., 149, 292–300, 1990. a, b

Haverkamp, R., Ross, P. J., Smettem, K. R. J., and Parlange, J. Y.: 3-Dimensional analysis of infiltration from the disc infiltrometer. 2. Physically-based infiltration equation, Water Resour. Res., 30, 2931–2935, 1994. a

Haverkamp, R., Leij, F. J., Fuentes, C., Sciortino, A., and Ross, P.: Soil water retention, Soil Sci. Soc. Am. J., 69, 1881–1890, 2005. a

Hillel, D.: Environmental soil physics: Fundamentals, applications, and environmental considerations, Academic Press, ISBN 0-12-348525-8, 1998. a

Kosugi, K.: Lognormal Distribution Model for Unsaturated Soil Hydraulic Properties, Water Resour. Res., 32, 2697–2703, https://doi.org/10.1029/96WR01776, 1996. a

Lassabatere, L.: Scilab script for mixed formulation for an easy and robust computation of sorptivity (version 2), Zenodo [code], https://doi.org/10.5281/zenodo.7044940, 2022. a, b, c

Lassabatere, L., Angulo-Jaramillo, R., Soria Ugalde, J. M., Cuenca, R., Braud, I., and Haverkamp, R.: Beerkan estimation of soil transfer parameters through infiltration experiments – BEST, Soil Sci. Soc. Am. J., 70, 521–532, https://doi.org/10.2136/sssaj2005.0026, 2006. a, b

Lassabatere, L., Angulo-Jaramillo, R., Soria-Ugalde, J., ǎimůnek, J., and Haverkamp, R.: Numerical evaluation of a set of analytical infiltration equations, Water Resour. Res., 45, W12415, https://doi.org/10.1029/2009WR007941, 2009. a, b, c

Lassabatere, L., Yilmaz, D., Peyrard, X., Peyneau, P., Lenoir, T., Šimůnek, J., and Angulo-Jaramillo, R.: New analytical model for cumulative infiltration into dual-permeability soils, Vadose Zone J., 13, 1–15, https://doi.org/10.2136/vzj2013.10.0181, 2014. a, b

Lassabatere, L., Di Prima, S., Bouarafa, S., Iovino, M., Bagarello, V., and Angulo-Jaramillo, R.: BEST-2K Method for Characterizing Dual-Permeability Unsaturated Soils with Ponded and Tension Infiltrometers, Vadose Zone J., 18, 180124, https://doi.org/10.2136/vzj2018.06.0124, 2019. a

Lassabatere, L., Peyneau, P.-E., Yilmaz, D., Pollacco, J., Fernández-Gálvez, J., Latorre, B., Moret-Fernández, D., Di Prima, S., Rahmati, M., Stewart, R. D., Abou Najm, M., Hammecker, C., and Angulo-Jaramillo, R.: A scaling procedure for straightforward computation of sorptivity, Hydrol. Earth Syst. Sci., 25, 5083–5104, https://doi.org/10.5194/hess-25-5083-2021, 2021. a, b, c, d, e, f, g, h, i, j, k, l, m

Malone, R. W., Shipitalo, M. J., and Meek, D. W.: Relationship between herbicide concentration in percolate, percolate breakthrough time, and number of active macropores, T. ASAE, 47, 1453–1456, 2004. a

Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res., 12, 513–522, 1976. a

Nasta, P., Assouline, S., Gates, J. B., Hopmans, J. W., and Romano, N.: Prediction of unsaturated relative hydraulic conductivity from Kosugi's water retention function, Proced. Environ. Sci., 19, 609–617, 2013. a

Parlange, J.-Y.: On Solving the Flow Equation in Unsaturated Soils by Optimization: Horizontal Infiltration, Soil Sci. Soc. Am. J., 39, 415–418, 1975. a, b, c

Philip, J. and Knight, J.: On solving the unsaturated flow equation: 3. New quasi-analytical technique, Soil Sci., 117, 1–13, 1974. a

Pollacco, J., Nasta, P., Soria-Ugalde, J., Angulo-Jaramillo, R., Lassabatere, L., Mohanty, B., and Romano, N.: Reduction of feasible parameter space of the inverted soil hydraulic parameter sets for Kosugi model, Soil Sci., 178, 267–280, https://doi.org/10.1097/SS.0b013e3182a2da21, 27, 2013. a, b

Radcliffe, D. E. and Simunek, J.: Soil physics with HYDRUS: Modeling and applications, CRC Press, ISBN 978-1-4200-7380-5, 2010. a

Ross, P., Haverkamp, R., and Parlange, J.-Y.: Calculating parameters for infiltration equations from soil hydraulic functions, Transp. Porous Media, 24, 315–339, 1996. a, b, c, d, e

Šimůnek, J., Jarvis, N. J., van Genuchten, M. T., and Gärdenäs, A.: Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone, J. Hydrol., 272, 14–35, 2003. a

Stewart, R. D. and Abou Najm, M. R.: A comprehensive model for single ring infiltration II: Estimating field-saturated hydraulic conductivity, Soil Sci. Soc. Am. J., 82, 558–567, https://doi.org/10.2136/sssaj2017.09.0313, 2018. a, b

Stewart, R. D., Rupp, D. E., Abou Najm, M. R., and Selker, J. S.: Modeling effect of initial soil moisture on sorptivity and infiltration, Water Resour. Res., 49, 7037–7047, https://doi.org/10.1002/wrcr.20508, 2013. a

Timlin, D. J., Ahuja, L. R., and Ankeny, M. D.: Comparison of 3 field methods to characterize apparent macropore conductivity, Soil Sci. Soc. Am. J., 58, 278–284, 1994. a

van Genuchten, M.: A closed-form equation for predicting the hydraulic conductivity of unsaturated soils, Soil Sci. Soc. Am. J., 44, 892–898, 1980. a, b

Varado, N., Braud, I., Ross, P. J., and Haverkamp, R.: Assessment of an efficient numerical solution of the 1D Richards' equation on bare soil, J. Hydrol., 323, 244–257, 2006. a, b

White, I. and Sully, M. J.: Macroscopic and microscopic capillary length and time scales from field infiltration, Water Resour. Res., 23, 1514–1522, 1987. a

Yilmaz, D.: Alternative α* parameter estimation for simplified Beerkan infiltration method to assess soil saturated hydraulic conductivity, Eurasian Soil Sci., 54, 1049–1058, https://doi.org/10.1134/S1064229321070140, 2021.  a

Yilmaz, D., Lassabatere, L., Angulo-Jaramillo, R., Deneele, D., and Legret, M.: Hydrodynamic characterization of basic oxygen furnace slag through an adapted best method, Vadose Zone J., 9, 107–116, https://doi.org/10.2136/vzj2009.0039, 2010. a, b

Short summary
Sorptivity is one of the most important parameters for quantifying water infiltration into soils. In this study, we propose a mixed formulation that avoids numerical issues and allows for the computation of sorptivity for all types of models chosen for describing the soil hydraulic functions and all initial and final conditions. We show the benefits of using the mixed formulation with regard to modeling water infiltration into soils.