Mixed formulation for an easy and robust numerical computation of sorptivity

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\to -\mathrm{\infty }$. 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 Minasny, 2011). 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 Najm, 2018). 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; Yilmaz, 2021). Sorptivity is also required for the computation of several hydraulic parameters, like the macroscopic capillary length (Bouwer, 1964; White and Sully, 1987).

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

$$\begin{array}{}\text{(1)}& S^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{1}}\right)=\mathrm{2}\underset{{\mathit{\theta }}_{\mathrm{0}}}{\overset{{\mathit{\theta }}_{\mathrm{1}}}{\int }}{\displaystyle \frac{\left(\mathit{\theta }-{\mathit{\theta }}_{\mathrm{0}}\right)}{F\left(\mathit{\theta }\right)}}D\left(\mathit{\theta }\right)d\mathit{\theta },\end{array}$$

where $D\left(\mathit{\theta }\right)=K\left(\mathit{\theta }\right)dh/d\mathit{\theta }$ is the hydraulic diffusivity function, and θ₀ and θ₁ stand for the initial and final water contents. In the context of water infiltration into soils, the initial water content, θ₀, refers to the water content along the soil profile before water infiltration, and the final water content, θ₁, 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, S_D, and a conductivity form, S_K, both summarized below:

$$\begin{array}{}\text{(2)}& {\displaystyle }& {\displaystyle }F\left(\mathit{\theta }\right)={\displaystyle \frac{\mathrm{2}\left(\mathit{\theta }-{\mathit{\theta }}_{\mathrm{0}}\right)}{\left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)}}\text{(3)}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{D}}^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{1}}\right)=\underset{{\mathit{\theta }}_{\mathrm{0}}}{\overset{{\mathit{\theta }}_{\mathrm{1}}}{\int }}\left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)D\left(\mathit{\theta }\right)d\mathit{\theta }\text{(4)}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{D}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{1}}\right)=\sqrt{\underset{{\mathit{\theta }}_{\mathrm{0}}}{\overset{{\mathit{\theta }}_{\mathrm{1}}}{\int }}\left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)D\left(\mathit{\theta }\right)d\mathit{\theta }}\text{(5)}& {\displaystyle }& {\displaystyle }{S}_{\mathrm{K}}\left(h_{\mathrm{0}},h_{\mathrm{1}}\right)=\sqrt{\underset{h_{\mathrm{0}}}{\overset{h_{\mathrm{1}}}{\int }}\left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)K\left(h\right)dh},\end{array}$$

where the initial and the final values of the water pressure heads, h₀ and h₁, correspond to the water contents θ₀=θ(h₀) and θ₁=θ(h₁).

These two forms, the diffusivity form, S_D, and the conductivity form, S_K, each have their own shortcomings. For certain hydraulic models, D(θ) tends towards infinity when θ₀→θ_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, h_a, S_D misses the saturated part of sorptivity, $\underset{h_{\mathrm{a}}}{\overset{h_{\mathrm{1}}}{\int }}({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }(h)-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}})K\left(h\right)dh$ (Ross et al., 1996). The conductivity form S_K 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):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}{S}_{\mathrm{K}}& \left(h_{\mathrm{0}},h_{\mathrm{1}}\ge h_{\mathrm{a}}\right)=\sqrt{S_{\mathrm{D}}^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{s}}\right)+\mathrm{2}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}\right)\left(h_{\mathrm{1}}-h_{\mathrm{a}}\right)K_{\mathrm{s}}}\\ & =\sqrt{S_{\mathrm{D}}^{\mathrm{2}}\left(\mathit{\theta }\left(h_{\mathrm{0}}\right),\mathit{\theta }\left(h_{\mathrm{1}}\right)\right)+\mathrm{2}\left(\mathit{\theta }\left(h_{\mathrm{1}}\right)-\mathit{\theta }\left(h_{\mathrm{0}}\right)\right)\left(h_{\mathrm{1}}-h_{\mathrm{a}}\right)K_{\mathrm{s}}}.\end{array}\end{array}$$

We can thus conclude that the conductivity form, S_K, is the most general equation. However, S_K 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 $h_{\mathrm{0}}\to -\mathrm{\infty }$. Then, the calculation of S_K requires the evaluation of an integral that involves an infinite lower bound: $\underset{-\mathrm{\infty }}{\overset{h_{\mathrm{1}}}{\int }}\left(\mathit{\theta }\left(h_{\mathrm{1}}\right)+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{r}}\right)K\left(h\right)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, θ₀=θ_r (corresponding to a negative initial water pressure head that tends towards infinity, $h_{\mathrm{0}}\to -\mathrm{\infty }$), and a final water pressure head higher than the air-entry water pressure head, h₁≥h_a, 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 S_D (Eq. 4) or S_K (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(-\mathrm{\infty },\mathrm{0})=\sqrt{\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}\left({\mathit{\theta }}_{\mathrm{s}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{r}}\right)K\left(h\right)dh}$. The maximum sorptivity, $S(-\mathrm{\infty },\mathrm{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 h₁>h_a. We considered three commonly used hydraulic models, for which Lassabatere et al. (2021) proposed analytical formulations for $S(-\mathrm{\infty },\mathrm{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 $S_{\mathrm{K}}\left(h_{\mathrm{0}},h_{\mathrm{1}}=\mathrm{0}\right)$ with the condition of null water pressure head at the surface, i.e., h₁=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, S_M, we start with the conductivity form of sorptivity, S_K, since it includes both unsaturated and saturated parts. Then, we define an intermediate water pressure head between the initial and final water pressure heads, $h_{\mathrm{c}}\in \left[h_{\mathrm{0}},h_{\mathrm{1}}\right]$, smaller than the air-entry pressure, $h_{\mathrm{c}}<h_{\mathrm{a}}\le \mathrm{0}$, and we split the integral into two separate parts as follows:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& S_{\mathrm{M}}\left(h_{\mathrm{0}},h_{\mathrm{1}}\right)=S_{\mathrm{K}}\left(h_{\mathrm{0}},h_{\mathrm{1}}\right)\\ & =\sqrt{\underset{h_{\mathrm{0}}}{\overset{h_{\mathrm{1}}}{\int }}\left(\mathit{\theta }\left(h_{\mathrm{1}}\right)+\mathit{\theta }\left(h\right)-\mathrm{2}\mathit{\theta }\left(h_{\mathrm{0}}\right)\right)K\left(h\right)dh}\\ & =\sqrt{\underset{h_{\mathrm{0}}}{\overset{h_{\mathrm{c}}}{\int }}\left(\mathit{\theta }\left(h_{\mathrm{1}}\right)+\mathit{\theta }\left(h\right)-\mathrm{2}\mathit{\theta }\left(h_{\mathrm{0}}\right)\right)K\left(h\right)dh+\underset{h_{\mathrm{c}}}{\overset{h_{\mathrm{1}}}{\int }}\left(\mathit{\theta }\left(h_{\mathrm{1}}\right)+\mathit{\theta }\left(h\right)-\mathrm{2}\mathit{\theta }\left(h_{\mathrm{0}}\right)\right)K\left(h\right)dh}\\ & =\sqrt{\underset{\mathit{\theta }\left(h_{\mathrm{0}}\right)}{\overset{\mathit{\theta }\left(h_{\mathrm{c}}\right)}{\int }}\left(\mathit{\theta }\left(h_{\mathrm{1}}\right)+\mathit{\theta }-\mathrm{2}\mathit{\theta }\left(h_{\mathrm{0}}\right)\right)D\left(\mathit{\theta }\right)d\mathit{\theta }+\underset{h_{\mathrm{c}}}{\overset{h_{\mathrm{1}}}{\int }}\left(\mathit{\theta }\left(h_{\mathrm{1}}\right)+\mathit{\theta }\left(h\right)-\mathrm{2}\mathit{\theta }\left(h_{\mathrm{0}}\right)\right)K\left(h\right)dh}.\end{array}\end{array}$$

In Eq. (7), the integral $\underset{h_{\mathrm{0}}}{\overset{h_{\mathrm{c}}}{\int }}\left(\mathit{\theta }\left(h_{\mathrm{1}}\right)+\mathit{\theta }\left(h\right)-\mathrm{2}\mathit{\theta }\left(h_{\mathrm{0}}\right)\right)K\left(h\right)dh$ is transformed into $\underset{\mathit{\theta }\left(h_{\mathrm{0}}\right)}{\overset{\mathit{\theta }\left(h_{\mathrm{c}}\right)}{\int }}\left(\mathit{\theta }\left(h_{\mathrm{1}}\right)+\mathit{\theta }-\mathrm{2}\mathit{\theta }\left(h_{\mathrm{0}}\right)\right)D\left(\mathit{\theta }\right)d\mathit{\theta }$ thanks to the change of variable h→θ. This operation requires that the function θ(h) is bijective over the whole interval [h₀,h_c], which is valid so long as h_c<h_a. Alternatively, the mixed formulation, S_M, may be written alternatively as follows:

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& S_{\mathrm{M}}\left(h_{\mathrm{0}},h_{\mathrm{1}}\right)=\\ & \sqrt{\underset{A}{\underbrace{\underset{{\mathit{\theta }}_{\mathrm{0}}}{\overset{{\mathit{\theta }}_{\mathrm{c}}}{\int }}\left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)D\left(\mathit{\theta }\right)d\mathit{\theta }}}+\underset{B}{\underbrace{\underset{h_{\mathrm{c}}}{\overset{h_{\mathrm{1}}}{\int }}\left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)K\left(h\right)dh}}},\end{array}\end{array}$$

where θ_c=θ(h_c), θ₁=θ(h₁), and θ₀=θ(h₀). The constraint h_c<h_a ensures that the computation of A in Eq. (8) avoids the challenging integration of infinite diffusivity close to saturation, $D\left({\mathit{\theta }}_{\mathrm{s}}\right)=+\mathrm{\infty }$, since θ_c<θ_s. In addition, h_c 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.

Figure 1 Concept of the mixed formulation, S_M(h₀,h₁): the integration of $\underset{h_{\mathrm{0}}}{\overset{h_{\mathrm{1}}}{\int }}({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }(h)-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}})K\left(h\right)dh$ ($=S_{\mathrm{M}}^{\mathrm{2}}(h_{\mathrm{0}},h_{\mathrm{1}})$) is converted into the sum of the integration of two bounded functions over bounded intervals, $\underset{{\mathit{\theta }}_{\mathrm{0}}}{\overset{{\mathit{\theta }}_{\mathrm{c}}}{\int }}({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}})D\left(\mathit{\theta }\right)d\mathit{\theta }$ and $\underset{h_{\mathrm{c}}}{\overset{h_{\mathrm{1}}}{\int }}({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }(h)-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}})K\left(h\right)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, S_M(h₀,h₁), for a synthetic loamy soil with an initial water pressure head of −10³ mm and a final water pressure head of $-{\mathrm{10}}^{-\mathrm{1}}$ mm.

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:

$$\begin{array}{}\text{(9)}& \left\{\begin{array}{l}{S}_{\mathrm{e}}=\frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}\\ h^{*}=\frac{h}{\left|h_{\mathrm{g}}\right|}\\ K_{\mathrm{r}}=\frac{K}{K_{\mathrm{s}}}\end{array}\right.,\end{array}$$

where S_e is the saturation degree, h^* is the scaled water pressure head, K_r is the relative hydraulic conductivity, θ_r and θ_s are the residual and the saturated soil water contents, h_g is the scale parameter for the water pressure head, and K_s 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):

$$\begin{array}{}\text{(10)}& S=\sqrt{\left|h_{\mathrm{g}}\right|K_{\mathrm{s}}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)}{S}^{*},\end{array}$$

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 S_M, defined by Eq. (7), leads to the final expression proposed in our study:

$$\begin{array}{}\text{(11)}& \left\{\begin{array}{ll}{S}_{\mathrm{M}}\left(h_{\mathrm{0}},h_{\mathrm{1}}\right)=& \sqrt{\left|h_{\mathrm{g}}\right|K_{\mathrm{s}}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)}{S}_{\mathrm{M}}^{*}\left(h_{\mathrm{0}}^{*},h_{\mathrm{1}}^{*}\right)\\ S_{\mathrm{M}}^{*}\left(h_{\mathrm{0}}^{*},h_{\mathrm{1}}^{*}\right)=& \left(\underset{S_{\mathrm{e}}\left(h_{\mathrm{0}}^{*}\right)}{\overset{S_{\mathrm{e}}\left(h_{\mathrm{c}}^{*}\right)}{\int }}\left(S_{\mathrm{e}}\left(h_{\mathrm{1}}^{*}\right)+S_{\mathrm{e}}-\mathrm{2}{S}_{\mathrm{e}}\left(h_{\mathrm{0}}^{*}\right)\right)D^{*}\left(S_{\mathrm{e}}\right)d{S}_{\mathrm{e}}\right.\\ & {\left.\mathrm{\dots }+\underset{h_{\mathrm{c}}^{*}}{\overset{h_{\mathrm{1}}^{*}}{\int }}\left(S_{\mathrm{e}}\left(h_{\mathrm{1}}^{*}\right)+S_{\mathrm{e}}\left(h^{*}\right)-\mathrm{2}{S}_{\mathrm{e}}\left(h_{\mathrm{0}}^{*}\right)\right)K_{\mathrm{r}}\left(h^{*}\right)d{h}^{*}\right)}^{\mathrm{1}/\mathrm{2}}\end{array}\right.,\end{array}$$

where $S_{\mathrm{M}}^{*}\left(h_{\mathrm{0}}^{*},h_{\mathrm{1}}^{*}\right)$ is the scaled version of the proposed mixed formulation, S_M(h₀,h₁), with $h_{\mathrm{0}}^{*}=h_{\mathrm{0}}/\left|h_{\mathrm{g}}\right|$ and $h_{\mathrm{1}}^{*}=h_{\mathrm{1}}/\left|h_{\mathrm{g}}\right|$. Equation (11) can be demonstrated by changing the integration variable θ→S_e in the first and $h\to h^{*}$ in the second integral of Eq. (7). Since Eq. (11) cannot be evaluated with coding software that does not allow infinite values, $h_{\mathrm{0}}^{*}=-\mathrm{\infty }$, we replace the input h₀ (that may be infinite) with $S_{\mathrm{e},\mathrm{0}}=S_{\mathrm{e}}\left(h_{\mathrm{0}}^{*}\right)$, which always remains bounded in the final expression of the mixed formulation $S_{\mathrm{M}}^{*}$:

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& S_{\mathrm{M}}^{*}\left(S_{\mathrm{e},\mathrm{0}},h_{\mathrm{1}}^{*}\right)=\\ & \sqrt{\underset{S_{\mathrm{e},\mathrm{0}}}{\overset{S_{\mathrm{e}}\left(h_{\mathrm{c}}^{*}\right)}{\int }}\left(S_{\mathrm{e}}\left(h_{\mathrm{1}}^{*}\right)+S_{\mathrm{e}}-\mathrm{2}{S}_{\mathrm{e},\mathrm{0}}\right)D^{*}\left(S_{\mathrm{e}}\right)d{S}_{\mathrm{e}}+\underset{h_{\mathrm{c}}^{*}}{\overset{h_{\mathrm{1}}^{*}}{\int }}\left(S_{\mathrm{e}}\left(h_{\mathrm{1}}^{*}\right)+S_{\mathrm{e}}\left(h^{*}\right)-\mathrm{2}{S}_{\mathrm{e},\mathrm{0}}\right)K_{\mathrm{r}}\left(h^{*}\right)d{h}^{*}}.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(13)}& \left\{\begin{array}{lc}{h}_{\mathrm{c}}^{*}=-\min \left(\left|h^{*}\left(\frac{S_{\mathrm{e},\mathrm{0}}+S_{\mathrm{e},\mathrm{1}}}{\mathrm{2}}\right)\right|,{\mathrm{10}}^z\right)& z\in \mathbb{Z}\\ S_{\mathrm{e},\mathrm{c}}=S_{\mathrm{e}}\left(h_c^{*}\right)& \end{array}\right..\end{array}$$

When necessary, the value of z is varied until the two integrals in $S_{\mathrm{M}}^{*}$ (Eq. 12) converge. In most cases, $z\in \mathit{\{}-\mathrm{2},-\mathrm{1},\mathrm{0},\mathrm{1},\mathrm{2}\mathit{\}}$ ensures convergence regardless of soil type and situation. Note that for hydraulic models with non-null water entry pressure head, h_a<0, z should be fixed with z≥0 so as to ensure $-{\mathrm{10}}^z\le -\mathrm{1}$ and thus $h_{\mathrm{c}}^{*}\le h_{\mathrm{a}}^{*}=-\mathrm{1}$. This condition is necessary to ensure the bijectivity of the function S_e(h^*) over the interval $\left[h_{\mathrm{0}}^{*},h_{\mathrm{c}}^{*}\right]$, which is required for the use of Eq. (12).

In the following section, the mixed formulation $S_{\mathrm{M}}^{*}$ (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 θ₁=θ_s (or, $S_{\mathrm{e},\mathrm{1}}=\mathrm{1}$) for a null water pressure head at the surface, h₁=0, and for very dry initial conditions θ₀→θ_r (or, $h_{\mathrm{0}}^{*}\to -\mathrm{\infty }$).

2.2 Usual methods for computing sorptivity based on S_D and S_K

2.2.1 Computing sorptivity with S_K for very dry initial conditions, $h_{\mathrm{0}}\to -\mathrm{\infty }$

Regarding the computation of sorptivity for very dry initial conditions with S_K, one of the strategies found in the literature applies the regular definitions of sorptivity (Parlange, 1975), Eq. (5), to the case of very low values of h₀. Such an approach was used by Di Prima et al. (2020) to estimate S_K(h₀,h₁) for very dry soils. In this case, the maximum sorptivity, $S\left(-\mathrm{\infty },h_{\mathrm{1}}\right)$, is approached as follows:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}S\left(-\mathrm{\infty },h_{\mathrm{1}}\right)& =\sqrt{\underset{-\mathrm{\infty }}{\overset{h_{\mathrm{1}}}{\int }}\left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{r}}\right)K\left(h\right)dh}\\ & =\underset{h_{\mathrm{0}}\to -\mathrm{\infty }}{\lim }\sqrt{\underset{\underset{\mathrm{¯}}{h_{\mathrm{0}}}}{\overset{h_{\mathrm{1}}}{\int }}\left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }\left(h\right)-\mathrm{2}\mathit{\theta }\left(\underset{\mathrm{¯}}{h_{\mathrm{0}}}\right)\right)K\left(h\right)dh}\\ & =\underset{h_{\mathrm{0}}\to -\mathrm{\infty }}{\lim }{S}_{\mathrm{K}}\left(\underset{\mathrm{¯}}{h_{\mathrm{0}}},h_{\mathrm{1}}\right).\end{array}\end{array}$$

Figure 2 Illustration of the regular strategies for the estimation of the limits for the case of very dry conditions with the estimation of $S(-\mathrm{\infty },h_{\mathrm{1}})$ using either S_K(h₀,h₁), Eq. (5), (a) or S_K-V2(h₀,h₁), Eq. (16), (b). The integration proceeds from the given value of h₀ to the right, h≥h₀, corresponding to the solid zones, and the value of h₀ is lowered step by step to reach the limit $S(-\mathrm{\infty },h_{\mathrm{1}})$ (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 S_K(h₀,h₁) or S_K-V2(h₀,h₁) for a synthetic loamy soil with a final water pressure head of $h_{\mathrm{1}}=-\mathrm{1}$ mm; successive values are considered for the initial water pressure heads, $h_{\mathrm{0}}\in \mathit{\{}-{\mathrm{10}}^{\mathrm{1}},-{\mathrm{10}}^{\mathrm{2}},-{\mathrm{10}}^{\mathrm{3}}\mathit{\}}$ mm.

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 $\sqrt{\lim f\left(x\right)}=\lim \sqrt{f\left(x\right)}$, with $x\to \sqrt{x}$ being a continuous function. In practical applications of this method, S_K(h₀,h₁) is computed for decreasing values of h₀ 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\left(-\mathrm{\infty },h_{\mathrm{1}}\right)$. This option is quite easy to implement since it requires the user to only code the regular function S_K before applying it to very negative values of h₀.

We propose an alternative procedure that employs a specific integrand to compute the same limit (Fig. 2b). In this case, the water content θ₀ is set equal to θ_r in the integrand so as to correspond to the targeted initial conditions $\mathit{\theta }\left(h_{\mathrm{0}}=-\mathrm{\infty }\right)={\mathit{\theta }}_{\mathrm{r}}$:

$$\begin{array}{}\text{(15)}& \begin{array}{rl}S\left(-\mathrm{\infty },h_{\mathrm{1}}\right)& =\sqrt{\underset{-\mathrm{\infty }}{\overset{h_{\mathrm{1}}}{\int }}\left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{r}}\right)K\left(h\right)dh}\\ & =\underset{h_{\mathrm{0}}\to -\mathrm{\infty }}{\lim }\sqrt{\underset{\underset{\mathrm{¯}}{h_{\mathrm{0}}}}{\overset{h_{\mathrm{1}}}{\int }}\left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{r}}\right)K\left(h\right)dh}\\ & =\underset{h_{\mathrm{0}}\to -\mathrm{\infty }}{\lim }{S}_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}\left(\underset{\mathrm{¯}}{h_{\mathrm{0}}},h_{\mathrm{1}}\right)\end{array}\end{array}$$

with the specific function S_K-V2 defined as follows:

$$\begin{array}{}\text{(16)}& S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}\left(h_{\mathrm{0}},h_{\mathrm{1}}\right)=\sqrt{\underset{h_{\mathrm{0}}}{\overset{h_{\mathrm{1}}}{\int }}\left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{r}}\right)K\left(h\right)dh}.\end{array}$$

In comparison to S_K, the water content θ₀ is replaced with θ_r in the integrand for S_K-V2. We expect this modification to improve numerical convergence towards the lower integration limit since S_K-V2 directly integrates the right integrand (Fig. 2b). Conversely, S_K integrates a distinct integrand, i.e., $\left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }\left(h\right)-\mathrm{2}\mathit{\theta }\left(h_{\mathrm{0}}\right)\right)K\left(h\right)\ne \left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{r}}\right)K\left(h\right)$, thus involving an additional source of error (Fig. 2a). Briefly, S_K combines the error due to the integral bound $h_{\mathrm{0}}>-\mathrm{\infty }$ and the difference between its integrand and the targeted integrand. Note that the function S_K-V2 should be restricted to the evaluation of $S\left(-\mathrm{\infty },h_{\mathrm{1}}\right)$ and never used for the computation of other cases, i.e., $S\left(h_{\mathrm{0}}\ne -\mathrm{\infty },h_{\mathrm{1}}\right)$, since the water content in the integrand is fixed at θ_r, which corresponds exclusively to the case of $h_{\mathrm{0}}=-\mathrm{\infty }$.

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:

$$\begin{array}{}\text{(17)}& \begin{array}{rl}{S}^{*}& \left(-\mathrm{\infty },h_{\mathrm{1}}^{*}\right)\\ & =\underset{h_{\mathrm{0}}^{*}\to -\mathrm{\infty }}{\lim }\sqrt{\underset{\underset{\mathrm{¯}}{h_{\mathrm{0}}^{*}}}{\overset{h_{\mathrm{1}}^{*}}{\int }}\left(S_{\mathrm{e},\mathrm{1}}+S_{\mathrm{e}}\left(h^{*}\right)-\mathrm{2}{S}_{\mathrm{e}}\left(\underset{\mathrm{¯}}{h_{\mathrm{0}}^{*}}\right)\right)K_{\mathrm{r}}\left(h^{*}\right)d{h}^{*}}\\ & =\underset{h_{\mathrm{0}}^{*}\to -\mathrm{\infty }}{\lim }{S}_{\mathrm{K}}^{*}\left(\underset{\mathrm{¯}}{h_{\mathrm{0}}^{*}},h_{\mathrm{1}}^{*}\right)\end{array}\end{array}$$

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{S}^{*}\left(-\mathrm{\infty },h_{\mathrm{1}}^{*}\right)& =\underset{h_{\mathrm{0}}^{*}\to -\mathrm{\infty }}{\lim }\sqrt{\underset{\underset{\mathrm{¯}}{h_{\mathrm{0}}^{*}}}{\overset{h_{\mathrm{1}}^{*}}{\int }}\left(S_{\mathrm{e},\mathrm{1}}+S_{\mathrm{e}}\left(h^{*}\right)\right)K_{\mathrm{r}}\left(h^{*}\right)d{h}^{*}}\\ & =\underset{h_{\mathrm{0}}^{*}\to -\mathrm{\infty }}{\lim }{S}_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}\left(\underset{\mathrm{¯}}{h_{\mathrm{0}}^{*}},h_{\mathrm{1}}^{*}\right)\end{array}\end{array}$$

with $S_{\mathrm{K}}^{*}$ and $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$, the scaled versions of S_K and S_K-V2, defined as follows:

$$\begin{array}{}\text{(19)}& \left\{\begin{array}{l}{S}_{\mathrm{K}}^{*}\left(h_{\mathrm{0}}^{*},h_{\mathrm{1}}^{*}\right)=\sqrt{\underset{h_{\mathrm{0}}^{*}}{\overset{h_{\mathrm{1}}^{*}}{\int }}\left(S_{\mathrm{e},\mathrm{1}}+S_{\mathrm{e}}\left(h^{*}\right)-\mathrm{2}{S}_{\mathrm{e}}\left(h_{\mathrm{0}}^{*}\right)\right)K_{\mathrm{r}}\left(h^{*}\right)d{h}^{*}}\\ S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}\left(h_{\mathrm{0}}^{*},h_{\mathrm{1}}^{*}\right)=\sqrt{\underset{h_{\mathrm{0}}^{*}}{\overset{h_{\mathrm{1}}^{*}}{\int }}\left(S_{\mathrm{e},\mathrm{1}}+S_{\mathrm{e}}\left(h^{*}\right)\right)K_{\mathrm{r}}\left(h^{*}\right)d{h}^{*}}\end{array}\right..\end{array}$$

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

2.2.2 Computing sorptivity with S_D for null water pressure head at surface, h₁=0

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

$$\begin{array}{}\text{(20)}& \begin{array}{rl}{S}_{\mathrm{u}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{s}}\right)& =\sqrt{\underset{{\mathit{\theta }}_{\mathrm{0}}}{\overset{{\mathit{\theta }}_{\mathrm{s}}}{\int }}\left({\mathit{\theta }}_{\mathrm{s}}+\mathit{\theta }-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)D\left(\mathit{\theta }\right)d\mathit{\theta }}\\ & =\underset{{\mathit{\theta }}_{\mathrm{1}}\to {\mathit{\theta }}_{\mathrm{s}}}{\lim }\sqrt{\underset{{\mathit{\theta }}_{\mathrm{0}}}{\overset{\underset{\mathrm{¯}}{{\mathit{\theta }}_{\mathrm{1}}}}{\int }}\left(\underset{\mathrm{¯}}{{\mathit{\theta }}_{\mathrm{1}}}+\mathit{\theta }-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)D\left(\mathit{\theta }\right)d\mathit{\theta }}\\ & =\underset{{\mathit{\theta }}_{\mathrm{1}}\to {\mathit{\theta }}_{\mathrm{s}}}{\lim }{S}_{\mathrm{D}}\left({\mathit{\theta }}_{\mathrm{0}},\underset{\mathrm{¯}}{{\mathit{\theta }}_{\mathrm{1}}}\right).\end{array}\end{array}$$

Note that with this method, we can only account for the unsaturated part of sorptivity, $S_{\mathrm{u}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{s}}\right)=\sqrt{\underset{{\mathit{\theta }}_{\mathrm{0}}}{\overset{{\mathit{\theta }}_{\mathrm{s}}}{\int }}\left({\mathit{\theta }}_{\mathrm{s}}+\mathit{\theta }-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)D\left(\mathit{\theta }\right)d\mathit{\theta }}$, and we systematically miss the saturated portion of sorptivity $\mathrm{2}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}\right)K_{\mathrm{s}}\left(h_{\mathrm{1}}-h_{\mathrm{a}}\right)$, as mentioned in Sect. 1. The total sorptivity corresponds to the sum of its two components (see Eq. 6): $S\left(h_{\mathrm{0}},h_{\mathrm{1}}>h_{\mathrm{a}}\right)=\sqrt{S_{\mathrm{u}}^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{s}}\right)+\mathrm{2}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}\right)K_{\mathrm{s}}\left(h_{\mathrm{1}}-h_{\mathrm{a}}\right)}$. The subscript “u” in S_u(θ₀,θ_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 S_D corresponds to $\left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)D\left(\mathit{\theta }\right)$. Consequently, S_D 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 [θ₀,θ₁] instead of [θ₀,θ_s] (Fig. 3a, S_D). To correct this problem, we define a different estimator, S_D-V2, to integrate directly the targeted integrand, $\left({\mathit{\theta }}_{\mathrm{s}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)$ D(θ) (Fig. 3b, S_D-V2), and the following developments emerge:

$$\begin{array}{}\text{(21)}& \begin{array}{rl}{S}_{\mathrm{u}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{s}}\right)& =\sqrt{\underset{{\mathit{\theta }}_{\mathrm{0}}}{\overset{{\mathit{\theta }}_{\mathrm{s}}}{\int }}\left({\mathit{\theta }}_{\mathrm{s}}+\mathit{\theta }-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)D\left(\mathit{\theta }\right)d\mathit{\theta }}\\ & =\underset{{\mathit{\theta }}_{\mathrm{1}}\to {\mathit{\theta }}_{\mathrm{s}}}{\lim }\sqrt{\underset{{\mathit{\theta }}_{\mathrm{0}}}{\overset{\underset{\mathrm{¯}}{{\mathit{\theta }}_{\mathrm{1}}}}{\int }}\left({\mathit{\theta }}_{\mathrm{s}}+\mathit{\theta }-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)D\left(\mathit{\theta }\right)d\mathit{\theta }}\\ & =\underset{{\mathit{\theta }}_{\mathrm{1}}\to {\mathit{\theta }}_{\mathrm{s}}}{\lim }{S}_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{0}},\underset{\mathrm{¯}}{{\mathit{\theta }}_{\mathrm{1}}}\right)\end{array}\end{array}$$

with the function S_D-V2 defined as follows:

$$\begin{array}{}\text{(22)}& S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{1}}\right)=\sqrt{\underset{{\mathit{\theta }}_{\mathrm{0}}}{\overset{{\mathit{\theta }}_{\mathrm{1}}}{\int }}\left({\mathit{\theta }}_{\mathrm{s}}+\mathit{\theta }-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)D\left(\mathit{\theta }\right)d\mathit{\theta }}.\end{array}$$

Figure 3 Illustration of the regular strategies for estimating the limits for the case of saturation θ₁→θ_s, with the estimation of S_u(θ₀,θ_s) using either S_D(θ₀,θ₁), Eq. (4), (a) or S_D-V2(θ₀,θ₁), Eq. (22), (b). The integration proceeds from the given value of θ₁ to the left, θ≥θ₁, defining the solid zones. The values of θ₁ are increased step by step to reach S_u(θ₀,θ_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 S_D(θ₀,θ₁) or S_D-V2(θ₀,θ₁) for a synthetic loamy soil with an initial water content of θ₀=0.25; successive values are considered for the final water contents, ${\mathit{\theta }}_{\mathrm{1}}\in \mathit{\{}\mathrm{0.39},\mathrm{0.40},\mathrm{0.41}\mathit{\}}$.

As mentioned above, for S_K-V2, S_D-V2 should be only used for the determination of S_u(θ₀,θ_s) and not for the computation of sorptivity corresponding to other values of final water contents, since S_D-V2 integrates the right integrand if and only if θ₁=θ_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:

$$\begin{array}{}\text{(23)}& {\displaystyle }\begin{array}{rl}{S}_{\mathrm{u}}^{*}\left(S_{\mathrm{e},\mathrm{0}},\mathrm{1}\right)& =\underset{S_{\mathrm{e},\mathrm{1}}\to \mathrm{1}}{\lim }\sqrt{\underset{S_{\mathrm{e},\mathrm{0}}}{\overset{\underset{\mathrm{¯}}{S_{\mathrm{e},\mathrm{1}}}}{\int }}\left(\underset{\mathrm{¯}}{S_{\mathrm{e},\mathrm{1}}}+S_{\mathrm{e}}-\mathrm{2}{S}_{\mathrm{e},\mathrm{0}}\right)D^{*}\left(S_{\mathrm{e}}\right)d{S}_{\mathrm{e}}}\\ & =\underset{S_{\mathrm{e},\mathrm{1}}\to \mathrm{1}}{\lim }{S}_{\mathrm{D}}^{*}\left(S_{\mathrm{e},\mathrm{0}},\underset{\mathrm{¯}}{S_{\mathrm{e},\mathrm{1}}}\right)\end{array}\text{(24)}& {\displaystyle }\begin{array}{rl}{S}_{\mathrm{u}}^{*}\left(S_{\mathrm{e},\mathrm{0}},\mathrm{1}\right)& =\underset{S_{\mathrm{e},\mathrm{1}}\to \mathrm{1}}{\lim }\sqrt{\underset{S_{\mathrm{e},\mathrm{0}}}{\overset{\underset{\mathrm{¯}}{S_{\mathrm{e},\mathrm{1}}}}{\int }}\left(\mathrm{1}+S_{\mathrm{e}}-\mathrm{2}{S}_{\mathrm{e},\mathrm{0}}\right)D^{*}\left(S_{\mathrm{e}}\right)d{S}_{\mathrm{e}}}\\ & =\underset{S_{\mathrm{e},\mathrm{1}}\to \mathrm{1}}{\lim }{S}_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}\left(S_{\mathrm{e},\mathrm{0}},\underset{\mathrm{¯}}{S_{\mathrm{e},\mathrm{1}}}\right).\end{array}\end{array}$$

The scaled versions of the formulation S_D and S_D-V2, $S_{\mathrm{D}}^{*}$ and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$, can be written as

$$\begin{array}{}\text{(25)}& \left\{\begin{array}{l}{S}_{\mathrm{D}}^{*}\left(S_{\mathrm{e},\mathrm{0}},S_{\mathrm{e},\mathrm{1}}\right)=\sqrt{\underset{S_{\mathrm{e},\mathrm{0}}}{\overset{S_{\mathrm{e},\mathrm{1}}}{\int }}\left(S_{\mathrm{e},\mathrm{1}}+S_{\mathrm{e}}-\mathrm{2}{S}_{\mathrm{e},\mathrm{0}}\right)D^{*}\left(S_{\mathrm{e}}\right)d{S}_{\mathrm{e}}}\\ S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}\left(S_{\mathrm{e},\mathrm{0}},S_{\mathrm{e},\mathrm{1}}\right)=\sqrt{\underset{S_{\mathrm{e},\mathrm{0}}}{\overset{S_{\mathrm{e},\mathrm{1}}}{\int }}\left(\mathrm{1}+S_{\mathrm{e}}-\mathrm{2}{S}_{\mathrm{e},\mathrm{0}}\right)D^{*}\left(S_{\mathrm{e}}\right)d{S}_{\mathrm{e}}}\end{array}\right..\end{array}$$

In the following sections, we compare these previously used strategies, based on the use of $S_{\mathrm{K}}^{*}$, and $S_{\mathrm{D}}^{*}$, and the improved versions designed for the purpose of this study, $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$ and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$ with the proposed mixed formulation $S_{\mathrm{M}}^{*}$, 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 S_M (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 Corey, 1964) was among the first hydraulic models of soil physics (Hillel, 1998). 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:

  $$\begin{array}{}\text{(26)}& \left\{\begin{array}{l}{\mathit{\theta }}_{\mathrm{BC}}\left(h\right)=\left\{\begin{array}{ll}{\mathit{\theta }}_{\mathrm{s}}& h\ge h_{\mathrm{BC}}\\ {\mathit{\theta }}_{\mathrm{r}}+\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right){\left(\frac{h_{\mathrm{BC}}}{h}\right)}^{{\mathit{\lambda }}_{\mathrm{BC}}}& h<h_{\mathrm{BC}}\end{array}\right.\\ K_{\mathrm{BC}}\left(\mathit{\theta }\right)=K_{\mathrm{s}}{\left(\frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}\right)}^{{\mathit{\eta }}_{\mathrm{BC}}}\end{array}\right..\end{array}$$

• The van Genuchten–Burdine (vGB) model combines the van Genuchten (1980) model with the Burdine condition $\left(m=\mathrm{1}-\frac{\mathrm{2}}{n}\right)$ 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:

  $$\begin{array}{}\text{(27)}& \left\{\begin{array}{l}{\mathit{\theta }}_{\mathrm{vGB}}\left(h\right)={\mathit{\theta }}_{\mathrm{r}}+\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right){\left(\mathrm{1}+{\left(\frac{h}{h_{\mathrm{vGB}}}\right)}^{n_{\mathrm{vGB}}}\right)}^{-m_{\mathrm{vGB}}}\\ m_{\mathrm{vGB}}=\mathrm{1}-\frac{\mathrm{2}}{n_{\mathrm{vGB}}}\\ K_{\mathrm{vGB}}\left(\mathit{\theta }\right)=K_{\mathrm{s}}{\left(\frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}\right)}^{{\mathit{\eta }}_{\mathrm{vGB}}}\end{array}\right..\end{array}$$

• The van Genuchten–Mualem (vGM) model combines the van Genuchten (1980) model with Mualem's condition $\left(m=\mathrm{1}-\frac{\mathrm{1}}{n}\right)$ 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:

  $$\begin{array}{}\text{(28)}& \left\{\begin{array}{l}{\mathit{\theta }}_{\mathrm{vGM}}\left(h\right)={\mathit{\theta }}_{\mathrm{r}}+\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right){\left(\mathrm{1}+{\left(\frac{h}{h_{\mathrm{vGM}}}\right)}^{n_{\mathrm{vGM}}}\right)}^{-m_{\mathrm{vGM}}}\\ m_{\mathrm{vGM}}=\mathrm{1}-\frac{\mathrm{1}}{n_{\mathrm{vGM}}}\\ K_{\mathrm{vGM}}\left(\mathit{\theta }\right)=K_{\mathrm{s}}{\left(\frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}\right)}^{l_{\mathrm{vGM}}}{\left(\mathrm{1}-{\left(\mathrm{1}-{\left(\frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}\right)}^{\frac{\mathrm{1}}{m_{\mathrm{vGM}}}}\right)}^{m_{\mathrm{vGM}}}\right)}^{\mathrm{2}}\end{array}\right..\end{array}$$

• The Kosugi (KG) model relates the WR function to the soil pore size distribution assuming lognormal distributions (Kosugi, 1996). 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:

  $$\begin{array}{}\text{(29)}& \left\{\begin{array}{l}{\mathit{\theta }}_{\mathrm{KG}}\left(h\right)={\mathit{\theta }}_{\mathrm{r}}+\frac{\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)}{\mathrm{2}}\mathrm{erfc}\left(\frac{\ln \left(\frac{h}{h_{\mathrm{KG}}}\right)}{\sqrt{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{KG}}}\right)\\ K_{\mathrm{KG}}\left(\mathit{\theta }\right)=K_{\mathrm{s}}{\left(\frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}\right)}^{l_{\mathrm{KG}}}{\left(\frac{\mathrm{1}}{\mathrm{2}}\mathrm{erfc}\left({\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}\frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}\right)+\frac{{\mathit{\sigma }}_{\mathrm{KG}}}{\sqrt{\mathrm{2}}}\right)\right)}^{\mathrm{2}},\end{array}\right.\end{array}$$

  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, h_g (h_BC, h_vGB, h_vGM, or h_KG); and saturated hydraulic conductivity, K_s. The BC models involve a non-null air-entry water pressure head, h_BC, 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 h_g=h_BC. 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; m_vGB, n_vGB, and η_vGB for the vGB model; m_vGM, n_vGM, and l_vGM for the vGM model; and, lastly, σ_KG and l_KG 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):

$$\begin{array}{}\text{(30)}& \mathit{\eta }={\displaystyle \frac{\mathrm{2}}{\mathit{\lambda }}}+\mathrm{2}+p,\end{array}$$

where λ=λ_BC for the BC model and λ=mn for the vGB models. In addition, the shape parameters l_vGM and l_KG are usually fixed at $\frac{\mathrm{1}}{\mathrm{2}}$. 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:

  $$\begin{array}{}\text{(31)}& \left\{\begin{array}{l}{S}_{\mathrm{e},\mathrm{BC}}\left(h^{*}\right)=\left(\mathrm{1}-H\left(\mathrm{1}+h^{*}\right)\right)|h^{*}{|}^{-{\mathit{\lambda }}_{\mathrm{BC}}}+H\left(\mathrm{1}+h^{*}\right)\\ K_{\mathrm{r},\mathrm{BC}}\left(S_{\mathrm{e}}\right)=S_{\mathrm{e}}^{{\mathit{\eta }}_{\mathrm{BC}}}.\end{array}\right.\end{array}$$

• vGB model:

  $$\begin{array}{}\text{(32)}& \left\{\begin{array}{l}{S}_{\mathrm{e},\mathrm{vGB}}\left(h^{*}\right)={\left(\mathrm{1}+|h^{*}{|}^{n_{\mathrm{vGB}}}\right)}^{-m_{\mathrm{vGB}}}\\ \mathrm{with}{m}_{\mathrm{vGB}}=\mathrm{1}-\frac{\mathrm{2}}{n_{\mathrm{vGB}}}\\ K_{\mathrm{r},\mathrm{vGB}}\left(S_{\mathrm{e}}\right)=S_{\mathrm{e}}^{{\mathit{\eta }}_{\mathrm{vGB}}}.\end{array}\right.\end{array}$$

• vGM model:

  $$\begin{array}{}\text{(33)}& \left\{\begin{array}{l}{S}_{\mathrm{e},\mathrm{vGM}}\left(h^{*}\right)={\left(\mathrm{1}+|h^{*}{|}^{n_{\mathrm{vGM}}}\right)}^{-m_{\mathrm{vGM}}}\\ \mathrm{with}{m}_{\mathrm{vGM}}=\mathrm{1}-\frac{\mathrm{1}}{n_{\mathrm{vGM}}}\\ K_{\mathrm{r},\mathrm{vGM}}\left(S_{\mathrm{e}}\right)=S_{\mathrm{e}}^{l_{\mathrm{vGM}}}{\left(\mathrm{1}-{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m_{\mathrm{vGM}}}}\right)}^{m_{\mathrm{vGM}}}\right)}^{\mathrm{2}}.\end{array}\right.\end{array}$$

• KG model:

  $$\begin{array}{}\text{(34)}& \left\{\begin{array}{l}{S}_{\mathrm{e},\mathrm{KG}}\left(h^{*}\right)=\frac{\mathrm{1}}{\mathrm{2}}\mathrm{erfc}\left(\frac{\ln \left(\left|h^{*}\right|\right)}{\sqrt{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{KG}}}\right)\\ K_{\mathrm{r},\mathrm{KG}}\left(S_{\mathrm{e}}\right)=S_{\mathrm{e}}^{l_{\mathrm{KG}}}{\left(\frac{\mathrm{1}}{\mathrm{2}}\mathrm{erfc}\left({\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}{S}_{\mathrm{e}}\right)+\frac{{\mathit{\sigma }}_{\mathrm{KG}}}{\sqrt{\mathrm{2}}}\right)\right)}^{\mathrm{2}}.\end{array}\right.\end{array}$$

These hydraulic models have the following hydraulic diffusivity functions, $D^{*}\left(S_{\mathrm{e}}\right)=K\left(S_{\mathrm{e}}\right)\frac{\mathrm{d}{h}^{*}}{\mathrm{d}{S}_{\mathrm{e}}}$ (Lassabatere et al., 2021):

$$\begin{array}{}\text{(35)}& \left\{\begin{array}{l}{D}_{\mathrm{BC}}^{*}\left(S_{\mathrm{e}}\right)=\frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{BC}}}{S}_{\mathrm{e}}^{{\mathit{\eta }}_{\mathrm{BC}}-\left(\frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{BC}}}+\mathrm{1}\right)}\\ D_{\mathrm{vGB}}^{*}\left(S_{\mathrm{e}}\right)=\frac{\mathrm{1}-m_{\mathrm{vGB}}}{\mathrm{2}{m}_{\mathrm{vGB}}}{S}_{\mathrm{e}}^{{\mathit{\eta }}_{\mathrm{vGB}}-\frac{\mathrm{1}+m_{\mathrm{vGB}}}{\mathrm{2}{m}_{\mathrm{vGB}}}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m_{\mathrm{vGB}}}}\right)}^{-\frac{\mathrm{1}+m_{\mathrm{vGB}}}{\mathrm{2}}}\\ D_{\mathrm{vGM}}^{*}\left(S_{\mathrm{e}}\right)=\frac{\mathrm{1}-m_{\mathrm{vGM}}}{m_{\mathrm{vGM}}}{S}_{\mathrm{e}}^{l_{\mathrm{vGM}}-\frac{\mathrm{1}}{m_{\mathrm{vGM}}}}\left({\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m_{\mathrm{vGM}}}}\right)}^{-m_{\mathrm{vGM}}}\right.\\ \left.+{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m_{\mathrm{vGM}}}}\right)}^{m_{\mathrm{vGM}}}-\mathrm{2}\right)\\ D_{\mathrm{KG}}^{*}\left(S_{\mathrm{e}}\right)=\frac{\mathrm{1}}{\mathrm{2}}\sqrt{\frac{\mathit{\pi }}{\mathrm{2}}}{\mathit{\sigma }}_{\mathrm{KG}}{S}_{\mathrm{e}}^{l_{\mathrm{KG}}}{\left(\mathrm{erfc}\left({\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}{S}_{\mathrm{e}}\right)+\frac{{\mathit{\sigma }}_{\mathrm{KG}}}{\sqrt{\mathrm{2}}}\right)\right)}^{\mathrm{2}}\\ e^{{\left({\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}{S}_{\mathrm{e}}\right)\right)}^{\mathrm{2}}+\sqrt{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{KG}}{\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}{S}_{\mathrm{e}}\right)}.\end{array}\right.\end{array}$$

These equations are needed for the computation of the dimensionless sorptivity with the proposed mixed formulation $S_{\mathrm{M}}^{*\mathrm{2}}$ and the regular formulations $S_{\mathrm{K}}^{*\mathrm{2}}$, $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*\mathrm{2}}$, $S_{\mathrm{D}}^{*\mathrm{2}}$, and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*\mathrm{2}}$.

The studied hydraulic models exhibit contrasting and challenging features for the computation of sorptivity, including non-null water pressure heads $h_{\mathrm{a}}^{*}<\mathrm{0}$ and infinite hydraulic diffusivity close to saturation $\underset{S_{\mathrm{e}}\to \mathrm{1}}{\lim }{D}^{*}\left(S_{\mathrm{e}}\right)=+\mathrm{\infty }$. 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 $S_{\mathrm{e}}\left(h^{*}\right)=\frac{\mathrm{1}}{\mathrm{2}}$. 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, S_e, with the increase of $\left|h^{*}\right|$ (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 S_e with the increase of $\left|h^{*}\right|$. Lassabatere et al. (2021) defined the WR shape index x as follows:

$$\begin{array}{}\text{(36)}& \left\{\begin{array}{l}{x}_{\mathrm{BC}}=\frac{{\mathit{\lambda }}_{\mathrm{BC}}}{\mathrm{2}+{\mathit{\lambda }}_{\mathrm{BC}}}\\ x_{\mathrm{vGB}}=m_{\mathrm{vGB}}\\ x_{\mathrm{vGM}}=m_{\mathrm{vGM}}\\ x_{\mathrm{KG}}=\frac{\mathrm{1}}{\mathrm{1}+{\mathit{\sigma }}_{\mathrm{KG}}}\end{array}\right..\end{array}$$

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

$$\begin{array}{}\text{(37)}& \left\{\begin{array}{l}{c}_{\mathrm{p},\mathrm{BC}}\left(x\right)=\mathrm{2}+\frac{\mathrm{1}-x}{\mathrm{5}x+\mathrm{1}}+\frac{\mathrm{1}-x}{\mathrm{7}x+\mathrm{1}}\\ c_{\mathrm{p},\mathrm{vGB}}\left(x\right)=\mathrm{\Gamma }\left(\frac{\mathrm{3}-x}{\mathrm{2}}\right)\left[\frac{\mathrm{\Gamma }\left(\frac{\mathrm{1}+\mathrm{5}x}{\mathrm{2}}\right)}{\mathrm{\Gamma }\left(\mathrm{1}+\mathrm{2}x\right)}+\frac{\mathrm{\Gamma }\left(\frac{\mathrm{1}+\mathrm{7}x}{\mathrm{2}}\right)}{\mathrm{\Gamma }\left(\mathrm{1}+\mathrm{3}x\right)}\right]\\ c_{\mathrm{p},\mathrm{vGM}}\left(x\right)=\mathrm{\Gamma }(\mathrm{2}-x)\left[\frac{\mathrm{\Gamma }\left(\frac{\mathrm{3}}{\mathrm{2}}x\right)}{\left(\frac{\mathrm{3}}{\mathrm{2}}x-\mathrm{1}\right)\mathrm{\Gamma }\left(\frac{\mathrm{1}}{\mathrm{2}}x\right)}+\frac{\mathrm{\Gamma }\left(\frac{\mathrm{5}}{\mathrm{2}}x\right)}{\left(\frac{\mathrm{5}}{\mathrm{2}}x-\mathrm{1}\right)\mathrm{\Gamma }\left(\frac{\mathrm{3}}{\mathrm{2}}x\right)}\right]\\ +\left(\mathrm{1}-x\right)\left[\frac{\mathrm{\Gamma }\left(\frac{\mathrm{3}}{\mathrm{2}}x\right)\mathrm{\Gamma }(\mathrm{1}+x)}{\left(\frac{\mathrm{3}}{\mathrm{2}}x-\mathrm{1}\right)\mathrm{\Gamma }\left(\frac{\mathrm{5}}{\mathrm{2}}x\right)}+\frac{\mathrm{\Gamma }\left(\frac{\mathrm{5}}{\mathrm{2}}x\right)\mathrm{\Gamma }(\mathrm{1}+x)}{\left(\frac{\mathrm{5}}{\mathrm{2}}x-\mathrm{1}\right)\mathrm{\Gamma }\left(\frac{\mathrm{7}}{\mathrm{2}}x\right)}\right.\\ \left.-\mathrm{2}\left(\frac{\mathrm{1}}{\frac{\mathrm{3}}{\mathrm{2}}x-\mathrm{1}}+\frac{\mathrm{1}}{\frac{\mathrm{5}}{\mathrm{2}}x-\mathrm{1}}\right)\right]\end{array}\right..\end{array}$$

Note that no analytical formulation can be found for the case of Kosugi's hydraulic model, so c_p,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 $l_{\mathrm{vGM}}=l_{\mathrm{KG}}=\frac{\mathrm{1}}{\mathrm{2}}$.

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 S_K and S_D (Eqs. 5 and 4), or their improved versions, S_K-V2 and S_D-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 $\widehat{S}$ of the nominal dimensional sorptivity S, and the related scaled variables, ${\widehat{S}}^{*}$ and S^*, the following relations emerge between the relative errors, E_r, related to dimensional and scaled sorptivity:

$$\begin{array}{}\text{(38)}& \begin{array}{rl}{E}_{\mathrm{r}}\left(S\right)& ={\displaystyle \frac{\widehat{S}-S}{S}}\\ & ={\displaystyle \frac{{\widehat{S}}^{*}\cdot \sqrt{\left|h_{\mathrm{g}}\right|K_{\mathrm{s}}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)}-S^{*}\cdot \sqrt{\left|h_{\mathrm{g}}\right|K_{\mathrm{s}}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)}}{S^{*}\cdot \sqrt{\left|h_{\mathrm{g}}\right|K_{\mathrm{s}}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)}}}\\ & ={\displaystyle \frac{{\widehat{S}}^{*}-S^{*}}{S^{*}}}\\ & =E_{\mathrm{r}}\left(S^{*}\right),\end{array}\end{array}$$

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^{*}(-\mathrm{\infty },\mathrm{0})$, since it involves at the same time the two types of challenges, i.e., very dry initial conditions with infinite water pressure head, $h_{\mathrm{0}}=-\mathrm{\infty }$, and saturated final conditions with null water pressure head, h₁=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 $S_{\mathrm{M}}^{*}$ with the nominal (error-free) values of sorptivity, $S^{*}=S^{*}(-\mathrm{\infty },\mathrm{0})=\sqrt{c_{\mathrm{p}}}$, 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 $S_{\mathrm{M}}^{*}$ 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, $S_{\mathrm{M}}^{*}$, to evaluate the efficiency of those previously used strategies. We considered several scenarios for the use of $S_{\mathrm{K}}^{*}$, $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$, $S_{\mathrm{D}}^{*}$, and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$, with several values of the lower water pressure head, $h_{\mathrm{0}}^{*}$, used in Eqs. (17) and (18) and several values of the final saturation degree, S_e,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(h₀,h₁), 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 ($h_{\mathrm{f}}=-\mathrm{150}$ or h_f=30 mm) along with an initial water pressure head of −10 m ($h_{\mathrm{i}}=-{\mathrm{10}}^{\mathrm{4}}$ 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 (Lassabatere, 2022).

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 $\left|h^{*}\right|\le \mathrm{1}$, i.e., $h^{*}\ge -\mathrm{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 $S_{\mathrm{M}}^{*}$ compared to the usual use of $S_{\mathrm{D}}^{*}$ for the cases of problematic computation close to saturation (i.e., h₁→0 and θ₁→θ_s). Similarly, the accuracy of the commonly used $S_{\mathrm{K}}^{*}$ 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, S_e, 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 S_e 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 $S_{\mathrm{K}}^{*}$, due to the smooth decrease in S_e with h^* (Fig. 4, first column). Conversely, we expect more problems with the use of $S_{\mathrm{D}}^{*}$ for large values of x, with quasi-infinite values for the hydraulic diffusivity close to saturation, i.e., S_e→1 (Fig. 4h, l, and p).

Figure 4 Water retention (WR) and hydraulic conductivity (HC) curves for different values of the WR shape index x. Panels (a, e, i, m) show WR as S_e(h^*), panels (b, f, j, n) show HC as K_r(S_e), panels (c, g, k, o) show HC as K_r(h^*), and panels (d, h, l, p) show diffusivity as D^*(S_e); 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, m_vGM, m_vGB, and σ_KG were computed as a function of x using Eq. (36) with $l_{\mathrm{vGM}}=l_{\mathrm{KG}}=\frac{\mathrm{1}}{\mathrm{2}}$. Adapted from Lassabatere et al. (2021).

3.2 Validation of the proposed mixed formulation, $S_{\mathrm{M}}^{*}$, against the nominal sorptivity, $S^{*}(-\mathrm{\infty },\mathrm{0})$

The computations using $S_{\mathrm{M}}^{*}$ (Eqs. 12 and 13) with $h_{\mathrm{c}}^{*}=h^{*}\left(\frac{S_{\mathrm{e},\mathrm{0}}+S_{\mathrm{e},\mathrm{1}}}{\mathrm{2}}\right)$ were efficient in most cases, regardless of the value of the WR shape index x. The use of the threshold 10^z was necessary for the first value of the WR shape index x for the BC model (x_BC=0.01), the first two values for the vGB model ($x_{\mathrm{vGB}}\in \mathit{\{}\mathrm{0.01},\mathrm{0.02}\mathit{\}}$), and the first 17 values for the vGM model ($x_{\mathrm{vGM}}\in \mathit{\{}\mathrm{0.01},\mathrm{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=-\mathrm{1}$ had to be considered for x_vGM∈{0.02, …, 0.07} and $z=-\mathrm{2}$ for x_vGM=0.01. Note that, as long as convergence is obtained, the values of $S_{\mathrm{M}}^{*}(-\mathrm{\infty },\mathrm{0})$ do not depend on z. With this strategy involving a threshold, the proposed mixed formulation, $S_{\mathrm{M}}^{*}$, 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 x_KG (data not shown).

The relative error, E_r, between the estimates provided by the proposed mixed formulation, $S_{\mathrm{M}}^{*}(-\mathrm{\infty },\mathrm{0})$, and the targeted sorptivity, $S^{*}(-\mathrm{\infty },\mathrm{0})=\sqrt{c_{\mathrm{p}}}$, was analyzed in terms of means, standard deviations, and minimum and maximum values (Table 1):

$$\begin{array}{}\text{(39)}& E_{\mathrm{r}}={\displaystyle \frac{S_{\mathrm{M}}^{*}(-\mathrm{\infty },\mathrm{0})-\sqrt{c_{\mathrm{p}}}}{\sqrt{c_{\mathrm{p}}}}}.\end{array}$$

Table 1 Absolute values of relative errors, $\left|E_{\mathrm{r}}\right|$, between the proposed mixed formulation, $S_{\mathrm{M}}^{*}$, and the targeted scaled sorptivity, $S^{*}(-\mathrm{\infty },\mathrm{0})=\sqrt{c_{\mathrm{p}}}$, with the mean $\stackrel{\mathrm{‾}}{\left|E_{\mathrm{r}}\right|}$, the standard deviation (${\mathit{\sigma }}_{\left|{\mathrm{E}}_{\mathrm{r}}\right|}$), 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⁻¹⁶ corresponds to the relative precision of numbers in Scilab.

The accuracy of the proposed mixed formulation, $S_{\mathrm{M}}^{*}$, 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⁻¹³ for the BC and vGB models and on the order of 10⁻⁹ for the vGM model (Table 1, $\stackrel{\mathrm{‾}}{E_{\mathrm{r}}}$). The minimum errors were $<{\mathrm{10}}^{-\mathrm{15}}$ for all the models (Table 1, Min). The maximum errors are $\approx {\mathrm{10}}^{-\mathrm{12}}$ for the BC and the vGB models and $\approx {\mathrm{10}}^{-\mathrm{7}}$ for the vGM model (Table 1, Max). In other words, the mixed formulation, $S_{\mathrm{M}}^{*}$, provides extremely accurate estimations of the targeted scaled sorptivity $S^{*}(-\mathrm{\infty },\mathrm{0})$. The proposed mixed formulation, $S_{\mathrm{M}}^{*}$, 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^{*}(-\mathrm{\infty },\mathrm{0})$

In this section, we compare the estimates provided by the strategies commonly considered (i.e., $S_{\mathrm{D}}^{*}$, $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$, $S_{\mathrm{K}}^{*}$, and $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$) with the reference values of the targeted sorptivity $S^{*}=S^{*}(-\mathrm{\infty },\mathrm{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 $S_{\mathrm{M}}^{*}$ for the KG model. Indeed, no analytical expressions are available for this last model, whereas the estimates provided by $S^{*}=S_{\mathrm{M}}^{*}$ are very accurate for the BC, vGB, and vGM models (see Sect. 3.2), and thus $S_{\mathrm{M}}^{*}$ 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 $S_{\mathrm{K}}^{*}$ and $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$ towards S^*. These functions, in particular, $S_{\mathrm{K}}^{*}$, are often used without special attention regarding their accuracy. $S_{\mathrm{K}}^{*}$ and $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$ converge to the limit S^* when $\left|h_{\mathrm{0}}^{*}\right|$ becomes large enough (Fig. 5, left column). For instance, for the BC model with λ=0.56, the use of $S_{\mathrm{K}}^{*}$ with $h_{\mathrm{0}}^{*}=-\mathrm{10}$ and $h_{\mathrm{0}}^{*}=-\mathrm{100}$ has respective relative errors of $E_{\mathrm{r}}=-\mathrm{4.1}$ % and $E_{\mathrm{r}}=-\mathrm{15}$ % (Fig. 5a, dashed red line). The second estimator, $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$, converges much faster than $S_{\mathrm{K}}^{*}$. With the same values of $h_{\mathrm{0}}^{*}$, 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 $h_{\mathrm{0}}^{*}=-\mathrm{10}$, with high accuracy. Such an improvement results from setting the initial saturation degree, S_e,0, at its target value, $S_{\mathrm{e},\mathrm{0}}=\mathrm{0}$, in the integrand (see Eq. 18 versus Eq. 17). The same conclusions hold regardless of the hydraulic model (Fig. 5c, e, and g).

Figure 5 Relative errors of the regular functions $S_{\mathrm{D}}^{*}(\mathrm{0},S_{\mathrm{e},\mathrm{1}})$ and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}(\mathrm{0},S_{\mathrm{e},\mathrm{1}})$ (b, d, f, h) and $S_{\mathrm{K}}^{*}(h_{\mathrm{0}}^{*},\mathrm{0})$ and $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}(h_{\mathrm{0}}^{*},\mathrm{0})$ (a, c, e, g) towards $S^{*}(-\mathrm{\infty },\mathrm{0})$ for the four hydraulic models, for the following specific cases: BC model with λ_BC=56 (x_BC=0.22), vGB model with n_vGB=3 ($x_{\mathrm{vGB}}=\mathrm{1}/\mathrm{3}$), vGM model with n_vGM=2 ($x_{\mathrm{vGM}}=\mathrm{1}/\mathrm{2}$) with l_vGM=0.5, and KG model with σ_KG=1.5 (x_KG=0.4) and l_KG=0.5.

The convergence of the two functions, $S_{\mathrm{D}}^{*}$ and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$, is depicted for the case of BC functions in Fig. 5b. For both estimators, the estimates are far from the target values, S^*, with $\left|E_{\mathrm{r}}\right|$ close to 50 %. This large error results from the omission of the saturated part of sorptivity, as explained above (Sect. 2.2). By design, $S_{\mathrm{D}}^{*}$ and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$ converge towards the unsaturated part of sorptivity, $S_{\mathrm{u}}(\mathrm{0},\mathrm{1})=\sqrt{\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(\mathrm{1}+S_{\mathrm{e}}\right)D^{*}\left(S_{\mathrm{e}}\right)d{S}_{\mathrm{e}}}$, and miss the additional saturated part of the scaled sorptivity, $\mathrm{2}\left(S_{\mathrm{e},\mathrm{1}}-S_{\mathrm{e},\mathrm{0}}\right)K_{\mathrm{r}}\left(S_{\mathrm{e},\mathrm{1}}\right)\left(h_{\mathrm{1}}^{*}-h_{\mathrm{a}}^{*}\right)=\mathrm{2}$, given that $h_{\mathrm{1}}^{*}=\mathrm{0}$ ($S_{\mathrm{e},\mathrm{1}}=\mathrm{1}$), $h_{\mathrm{0}}^{*}=-\mathrm{\infty }$ ($S_{\mathrm{e},\mathrm{0}}=\mathrm{0}$), and $h_{\mathrm{a}}^{*}=\mathrm{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 $\mathrm{2}({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}})K_{\mathrm{s}}\left|h_{\mathrm{a}}\right|$ 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 $S_{\mathrm{D}}^{*}$ and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$ 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 $S_{\mathrm{e},\mathrm{1}}=\mathrm{0.99}$ (Fig. 5d, f and h). Even when $S_{\mathrm{e},\mathrm{1}}=\mathrm{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 $S_{\mathrm{K}}^{*}$ and $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$, 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 $S_{\mathrm{K}}^{*}$ and $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$, on the one hand, and $S_{\mathrm{D}}^{*}$ and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$, on the other, varies with the shape parameters and indices. Figure 6 depicts the relative errors of the estimators $S_{\mathrm{K}}^{*}$ and $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$ as a function of the shape parameters for different values of $h_{\mathrm{0}}^{*}$ (left column of Fig. 6) and of the estimators $S_{\mathrm{D}}^{*}$ and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$ for different values of S_e,1 (right column of Fig. 6). For all hydraulic models except the BC model, the best estimates for $S_{\mathrm{D}}^{*}$ and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$ 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, S_e,1, to a minimum of 0.999, to get errors less than 10 %. However, estimates remain poor, even with $S_{\mathrm{e},\mathrm{1}}=\mathrm{0.9999}$ or $S_{\mathrm{e},\mathrm{1}}=\mathrm{0.99999}$, when the WR shape indices tend towards unity (Fig. 6d, f, and h).

Figure 6 Relative errors of estimators $S_{\mathrm{D}}^{*}$ and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$ (b, d, f, h) and $S_{\mathrm{K}}^{*}$ and $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$ (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.

For the BC model, the estimators $S_{\mathrm{D}}^{*}$ and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$ provide poor estimates under all circumstances, since they miss the saturated part of sorptivity, as explained above (Fig. 6b, $S_{\mathrm{D}}^{*}\approx S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$). Adding the saturated part of sorptivity substantially improves the computation, with very accurate estimations (Fig. 6b, $\sqrt{\mathrm{2}+S_{\mathrm{D}}^{*\mathrm{2}}}$). Regarding the estimators $S_{\mathrm{K}}^{*}$ and $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$ (Fig. 6, left column), the estimates improve when the WR shape index converges towards 1. It can be noted that $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$ 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 $S_{\mathrm{K}}^{*}$, $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$, $S_{\mathrm{D}}^{*}$, and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$. The approaches that compute sorptivity by integration with regards to h^* (i.e., $S_{\mathrm{K}}^{*}$ and $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$) have errors that come from the choice of the lower integration boundary $h_{\mathrm{0}}^{*}$ and from the value of S_e,0 in the integrand. Solutions that compute sorptivity by integration with regards to S_e (i.e., $S_{\mathrm{D}}^{*}$ and $S_{\mathrm{D}\text{-}\mathrm{V}\mathrm{2}}^{*}$) 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 $S_{\mathrm{K}}^{*}$ estimator is better than $S_{\mathrm{D}}^{*}$, and the $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$ estimate is superior to both. Indeed, $S_{\mathrm{K}\text{-}\mathrm{V}\mathrm{2}}^{*}$ does not miss the saturated part of sorptivity, like $S_{\mathrm{K}}^{*}$, 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, S_e,0, to its final values, substantially improves the quality of the estimator. This formulation thus always provides estimates with $\left|E_{\mathrm{r}}\right|<\mathrm{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, $S_{\mathrm{M}}^{*}$. Consequently, the use of the mixed formulation, $S_{\mathrm{M}}^{*}$, 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., h_i, and final conditions, i.e., h_f. In addition to the illustration of the practical use of the proposed mixed formulation, S_M, we demonstrate that (i) S_M provides better estimates for sorptivity than the usual estimator, S_D, 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 S_M to S_D, only, with no comparison with the other methods since S_D-V2 gives similar results to S_D and since S_K and S_K-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 Simunek, 2010). 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).

Figure 7 Water 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 S_M and the regular model S_D.

Table 2 Hydraulic parameters, i.e., values of the residual and saturated water contents, θ_r and θ_s; the scale parameter for water pressure head h_g; the air-entry water pressure head h_a; the saturated hydraulic conductivity K_s; 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).

At first, we illustrate the use of the mixed formulation, S_M, for the computation of sorptivity for a final water pressure head at the surface of $h_{\mathrm{f}}=-\mathrm{150}$ mm and an initial water pressure head of $h_{\mathrm{i}}=-{\mathrm{10}}^{\mathrm{4}}$ 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 S_M can be performed using the codes developed for this paper and downloadable from Zenodo (Lassabatere, 2022). Alternatively, computations may be performed as follows:

1. Compute the initial and final saturation degrees, $S_{\mathrm{e},\mathrm{i}}=S_{\mathrm{e}}\left(h_{\mathrm{i}}\right)$ and $S_{\mathrm{e},\mathrm{f}}=S_{\mathrm{e}}\left(h_{\mathrm{f}}\right)$.

2. Compute the intermediate water pressure head, $h_{\mathrm{c}}^{*}=-\min \left(\left|h^{*}\left(\frac{S_{\mathrm{e},\mathrm{i}}+S_{\mathrm{e},\mathrm{f}}}{\mathrm{2}}\right)\right|,{\mathrm{10}}^{\mathrm{z}}\right)$ with $z\in \mathit{\{}-\mathrm{2},-\mathrm{1},\mathrm{0},\mathrm{1},\mathrm{2}\mathit{\}}$.

3. Compute the intermediate saturation degree from the intermediate water pressure head, $S_{\mathrm{e},\mathrm{c}}=S_{\mathrm{e}}\left(h_{\mathrm{c}}^{*}\right)$.

4. Integrate the lower part of the squared sorptivity

   $$A=\underset{S_{\mathrm{e},\mathrm{i}}}{\overset{S_{\mathrm{e},\mathrm{c}}}{\int }}(S_{\mathrm{e},\mathrm{f}}+S_{\mathrm{e}}-\mathrm{2}{S}_{\mathrm{e},\mathrm{i}})D^{*}\left(S_{\mathrm{e}}\right)d{S}_{\mathrm{e}}.$$

5. Integrate the upper part of the squared sorptivity

   $$B=\underset{h_{\mathrm{c}}^{*}}{\overset{h_{\mathrm{f}}^{*}}{\int }}(S_{\mathrm{e},\mathrm{f}}+S_{\mathrm{e}}(h^{*})-\mathrm{2}{S}_{\mathrm{e},\mathrm{i}})K_{\mathrm{r}}\left(h^{*}\right)d{h}^{*}.$$

6. Combine the two parts to compute the mixed formulation $S_{\mathrm{M}}^{*}=\sqrt{A+B}$.

7. Upscale by multiplying the scaled sorptivity,

   $$S_{\mathrm{M}}=\sqrt{\left|h_{\mathrm{g}}\right|K_{\mathrm{s}}({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}})}{S}_{\mathrm{M}}^{*}.$$

For the studied case, with $h_{\mathrm{i}}=-\mathrm{10}$ m and $h_{\mathrm{f}}=-\mathrm{150}$ mm, it leads to the following results:

1. The related scaled water pressure heads are $h_{\mathrm{i}}^{*}=-\mathrm{36}$ and $h_{\mathrm{f}}^{*}=-\mathrm{0.54}$, with related saturation degrees of $S_{\mathrm{e},\mathrm{i}}=\mathrm{0.134}$ and $S_{\mathrm{e},\mathrm{f}}=\mathrm{0.89}$.

2. The intermediate water pressure head takes the value of $h_{\mathrm{c}}^{*}=-\mathrm{1}$, which corresponds to the maximum between $h^{*}\left(\frac{S_{\mathrm{e},\mathrm{i}}+S_{\mathrm{e},\mathrm{f}}}{\mathrm{2}}\right)=-\mathrm{2.96}$ and $-{\mathrm{10}}^{\mathrm{0}}=-\mathrm{1}$.

3. The corresponding saturation degree takes the value of $S_{\mathrm{e},\mathrm{c}}^{*}=\mathrm{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: $S_{\mathrm{M}}^{*}=\mathrm{0.2931}$.

7. The upscaling finalizes the computation of sorptivity: S_M=0.156 mm s${}^{-\mathrm{1}/\mathrm{2}}$.

Then, we computed sorptivity with the mixed formulation, S_M, and the regular estimator, S_D, as a function of the initial water pressure head, h_i, for the water pressure head at the surface of $h_{\mathrm{f}}=-\mathrm{150}$ mm (Fig. 7c) and as a function of final water pressure head, h_f, for an initial water pressure head of $h_{\mathrm{i}}=-{\mathrm{10}}^{\mathrm{4}}$ 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 S_M. The values of sorptivity decrease with the initial water pressure head, h_i (Fig. 7c), and increase with the final water pressure head, h_f (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 S_D (Fig. 7c and d, dashed lines). In most cases, the two estimators, S_M and S_D, provided the same estimations, except when the final pressure head was higher than the air-entry water pressure head, h_f≥h_a. In this case, as explained above, S_D missed the saturated part of the sorptivity, which explains why the values predicted with S_D no longer increase with h_f but remain equal to the unsaturated part of sorptivity: $S_{\mathrm{D}}\left({\mathit{\theta }}_{\mathrm{i}},\left(h_{\mathrm{f}}\ge h_{\mathrm{a}}\right)\right)=S_{\mathrm{D}}\left({\mathit{\theta }}_{\mathrm{i}},{\mathit{\theta }}_{\mathrm{s}}\right)=S_{\mathrm{u}}\left({\mathit{\theta }}_{\mathrm{i}},{\mathit{\theta }}_{\mathrm{s}}\right)$ (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 h_f=0, whereas the air-entry water pressure head is lower than zero for the BC model, $h_{\mathrm{a}}=-\mathrm{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 (h_f=30 mm) and an initial water pressure head of −10 m ($h_{\mathrm{i}}=-{\mathrm{10}}^{\mathrm{4}}$ mm). To model water infiltration, we used the analytical model developed for h_f≤h_a by Haverkamp et al. (1994) with its extension to h_f>h_a 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:

$$\begin{array}{}\text{(40)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{2}\mathrm{\Delta }{K}^{\mathrm{2}}}{S^{\mathrm{2}}}}t& ={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\mathit{\beta }}}\left({\displaystyle \frac{\mathrm{2}\mathrm{\Delta }K}{S^{\mathrm{2}}}}\left(I_{\mathrm{1}\mathrm{D}}-K_{\mathrm{i}}t\right)\right.\\ & \left.-\ln \left({\displaystyle \frac{e^{\frac{\mathrm{2}\mathit{\beta }\mathrm{\Delta }K}{S^{\mathrm{2}}}\left(I_{\mathrm{1}\mathrm{D}}-K_{\mathrm{i}}t\right)}+\mathit{\beta }-\mathrm{1}}{\mathit{\beta }}}\right)\right),\end{array}\end{array}$$

where ΔK stands for the difference between final and initial values of hydraulic conductivity $\mathrm{\Delta }K=K_{\mathrm{f}}-K_{\mathrm{i}}$, 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:

$$\begin{array}{}\text{(41)}& \left\{\begin{array}{ll}{\mathit{\gamma }}_{\mathrm{I}}=\frac{S^{\mathrm{2}}}{\mathrm{2}\mathrm{\Delta }K}& I={\mathit{\gamma }}_{\mathrm{I}}{I}^{*}+K_{\mathrm{i}}t\\ {\mathit{\gamma }}_t=\frac{S^{\mathrm{2}}}{\mathrm{2}\mathrm{\Delta }{K}^{\mathrm{2}}}& t={\mathit{\gamma }}_t{t}^{*}\end{array}\right..\end{array}$$

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

$$\begin{array}{}\text{(42)}& I_{\mathrm{QEI}}\left(t\right)={\mathit{\gamma }}_{\mathrm{I}}{I}_{\mathrm{QEI}}^{*}\left({\displaystyle \frac{t}{{\mathit{\gamma }}_t}}\right)+K_{\mathrm{i}}t\end{array}$$

with the following equations for the scaled QEI model, $I_{\mathrm{QEI}}^{*}\left(t^{*}\right)$:

$$\begin{array}{}\text{(43)}& t^{*}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\mathit{\beta }}}\left(I^{*}-\ln \left({\displaystyle \frac{e^{\mathit{\beta }{I}^{*}}+\mathit{\beta }-\mathrm{1}}{\mathit{\beta }}}\right)\right).\end{array}$$

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 h_f>h_a, leading to

$$\begin{array}{}\text{(44)}& {\displaystyle }{I}^{*}={\displaystyle \frac{\mathit{\sigma }}{q^{*}-\mathrm{1}}}+{\displaystyle \frac{\mathrm{1}-\mathit{\sigma }}{\mathit{\beta }}}\ln \left(\mathrm{1}+{\displaystyle \frac{\mathit{\beta }}{q^{*}-\mathrm{1}}}\right)\text{(45)}& {\displaystyle }\begin{array}{rl}{t}^{*}& ={\displaystyle \frac{\mathrm{1}-\mathit{\sigma }}{\mathit{\beta }(\mathrm{1}-\mathit{\beta })}}\ln \left(\mathrm{1}+{\displaystyle \frac{\mathit{\beta }}{q^{*}-\mathrm{1}}}\right)+{\displaystyle \frac{\mathit{\sigma }}{q^{*}-\mathrm{1}}}\\ & -{\displaystyle \frac{\mathrm{1}-\mathit{\sigma }\mathit{\beta }}{\mathrm{1}-\mathit{\beta }}}\ln \left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{q^{*}-\mathrm{1}}}\right),\end{array}\end{array}$$

where $q^{*}=\frac{q}{{\mathit{\gamma }}_q}$ is the scaled infiltration rate and γ_q=ΔK. 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^{*}\left(q^{*}\right)$, defined by Eq. (45), and the second relating the cumulative infiltration and the infiltration rate, $I^{*}=I^{*}\left(q^{*}\right)$, 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, $I_{\mathrm{QEI}\mathrm{\_}\mathrm{ext}}^{*}\left(t^{*}\right)$. 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:

$$\begin{array}{}\text{(46)}& \mathit{\sigma }={\displaystyle \frac{\mathrm{2}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{i}}\right)K_{\mathrm{s}}\left(h_{\mathrm{f}}-h_{\mathrm{a}}\right)}{S^{\mathrm{2}}}}.\end{array}$$

The QEI_ext model addresses all cases. When h_f≤h_a, the saturated part of sorptivity is null, leading to σ=0, and the QEI_ext model reduces to the regular QEI model. When h_f>h_a, 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.

Figure 8 Cumulative infiltrations as a function of initial (a) and final (b) water pressure heads, computed with the correct estimate for sorptivity (S_M) 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 (S_M versus S_D) for (c) the case of unsaturated conditions imposed at surface h_f≤h_a (BC model, $h_{\mathrm{i}}=-\mathrm{10}$ m and $h_{\mathrm{f}}=-\mathrm{150}$ mm) and (d) the case of saturated conditions, h_f≥h_a (BC model, $h_{\mathrm{i}}=-\mathrm{10}$ m and h_f=30 mm).

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 h_i (Fig. 8a, arrows downwards) and increases with the final water pressure head h_f (Fig. 8b, arrows upwards), as a consequence of the variations of sorptivity as a function of h_i and h_f (Fig. 7c and d). It can be noted that the impact of the water pressure head at surface h_f 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 $h_{\mathrm{i}}=-{\mathrm{10}}^{\mathrm{4}}$ mm and a final water pressure head of $h_{\mathrm{f}}=-\mathrm{150}$ mm (Fig. 8c) and for the case of the same initial water pressure head, $h_{\mathrm{i}}=-{\mathrm{10}}^{\mathrm{4}}$ mm, but a positive surface water pressure head h_f=30 mm (Fig. 8d). For the first case, there was no difference at all (Fig. 8c). Both S_D and S_M 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., S_M instead of S_D, 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 S_M 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 ($<{\mathrm{10}}^{-\mathrm{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.