# Mixed formulation for an easy and robust numerical computation of sorptivity

Laurent Lassabatere, Pierre-Emmanuel Peyneau, Deniz Yilmaz, Joseph Pollacco, Jesús Fernández-Gálvez, Borja Latorre, David Moret-Fernández, Simone Di Prima, Mehdi Rahmati, Ryan Stewart, et al.

► To cite this version:

Laurent Lassabatere, Pierre-Emmanuel Peyneau, Deniz Yilmaz, Joseph Pollacco, Jesús Fernández-Gálvez, et al.. Mixed formulation for an easy and robust numerical computation of sorptivity. Hydrology and Earth System Sciences, 2023, 27 (4), pp.895-915. 10.5194/hess-27-895-2023. hal-04047636

HAL Id: hal-04047636

https://hal.inrae.fr/hal-04047636

Submitted on 27 Mar 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Mixed formulation for an easy and robust numerical computation of sorptivity

Laurent Lassabatere<fnref n="1" />, Pierre-Emmanuel Peyneau<fnref n="2" />, Deniz Yilmaz<fnref n="3" />, Joseph Pollacco<fnref n="4" />, Jesús Fernández-Gálvez<fnref n="5" />,
Borja Latorre<fnref n="6" />, David Moret-Fernández<fnref n="6" />, Simone Di Prima<fnref n="7" />, Mehdi Rahmati<fnref n="8" />,<fnref n="9" />, Ryan D. Stewart<fnref n="10" />,
Majdi Abou Najm<fnref n="11" />, Claude Hammecker<fnref n="12" />, and Rafael Angulo-Jaramillo<fnref n="1" />

**Correspondence:** Laurent Lassabatere (laurent.lassabatere@entpe.fr)

Received: 17 December 2021 – Discussion started: 23 December 2021

Revised: 3 September 2022 – Accepted: 15 January 2023 – Published: 24 February 2023

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

<footnote n="1">¹Univ Lyon, Université Claude Bernard Lyon 1, CNRS, ENTPE, UMR 5023 LEHNA, 69518 Vaulx-en-Velin, France</footnote>

<footnote n="2">²Univ Gustave Eiffel, GERS-LEE, 44344 Bouguenais, France</footnote>

<footnote n="3">³Civil Engineering Department, Engineering Faculty, Munzur University, Tunceli, Türkiye</footnote>

<footnote n="4">⁴Manaaki Whenua - Landcare Research, 7640 Lincoln, New Zealand</footnote>

<footnote n="5">⁵Department of Regional Geographic Analysis and Physical Geography, University of Granada, 18071 Granada, Spain</footnote>

<footnote n="6">⁶Departamento de Suelo y Agua, Estación Experimental de Aula Dei, Consejo Superior de Investigaciones Científicas (CSIC), P.O. Box 13034, 50080 Zaragoza, Spain</footnote>

<footnote n="7">⁷School of Agricultural, Forestry, Food and Environmental Sciences (SAFE), University of Basilicata, 85100 Potenza, Italy</footnote>

<footnote n="8">⁸Department of Soil Science and Engineering, Faculty of Agriculture, University of Maragheh, Maragheh, Iran</footnote>

<footnote n="9">⁹Forschungszentrum Jülich GmbH, Institute of Bio- and Geosciences: Agrosphere (IBG-3), Jülich, Germany</footnote>

<footnote n="10">¹⁰School of Plant and Environmental Sciences, Virginia Polytechnic Institute and State University, Blacksburg, VA, United States</footnote>

<footnote n="11">¹¹Department of Land, Air and Water Resources, University of California, Davis, CA 95616, United States</footnote>

<footnote n="12">¹²University of Montpellier, UMR LISAH, IRD, Montpellier, France</footnote>
tion 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(\theta)$, as follows (Philip and Knight, 1974):

$$S^2(\theta_0, \theta_1) = 2 \int_{\theta_0}^{\theta_1} \frac{(\theta - \theta_0)}{F(\theta)} D(\theta) d\theta, \quad (1)$$

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

$$F(\theta) = \frac{2(\theta - \theta_0)}{(\theta_1 + \theta - 2\theta_0)} \quad (2)$$

$$S_D^2(\theta_0, \theta_1) = \int_{\theta_0}^{\theta_1} (\theta_1 + \theta - 2\theta_0) D(\theta) d\theta \quad (3)$$

$$S_D(\theta_0, \theta_1) = \sqrt{\int_{\theta_0}^{\theta_1} (\theta_1 + \theta - 2\theta_0) D(\theta) d\theta} \quad (4)$$

$$S_K(h_0, h_1) = \sqrt{\int_{h_0}^{h_1} (\theta_1 + \theta(h) - 2\theta_0) K(h) dh}, \quad (5)$$

where the initial and the final values of the water pressure heads, $h_0$ and $h_1$, correspond to the water contents $\theta_0 = \theta(h_0)$ and $\theta_1 = \theta(h_1)$.

These two forms, the diffusivity form, $S_D$, and the conductivity form, $S_K$, each have their own shortcomings. For certain hydraulic models, $D(\theta)$ tends towards infinity when $\theta_0 \to \theta_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, $\int_{h_a}^{h_1} (\theta_1 + \theta(h) - 2\theta_0) K(h) 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):

$$S_K(h_0, h_1 \ge h_a) = \sqrt{S_D^2(\theta_0, \theta_s) + 2(\theta_s - \theta_0)(h_1 - h_a)K_s} \\ = \sqrt{S_D^2(\theta_0, \theta_1) + 2(\theta_1 - \theta_0)(h_1 - h_a)K_s}. \quad (6)$$

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 $\theta_r$ corresponds to $h_0 \to -\infty$. Then, the calculation of $S_K$ requires the evaluation of an integral that involves an infinite lower bound:

$$\int_{-\infty}^{h_1} (\theta(h_1) + \theta(h) - 2\theta_r) K(h) dh.$$

In this study, we propose a new mixed formulation that overcomes these problems. We compare it to the approaches commonly used to compute sorptivity, i.e., Eqs. (4) and (5). The proposed mixed formulation automatically accounts for the saturated and unsaturated parts of sorptivity. It also allows for easy computation under any initial condition, including the extreme case of an initial water content equal to the residual water content, $\theta_0 = \theta_r$ (corresponding to a negative initial water pressure head that tends towards infinity, $h_0 \to -\infty$), and a final water pressure head higher than the air-entry water pressure head, $h_1 \ge 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(-\infty, 0) = \sqrt{\int_{-\infty}^{0} (\theta_s + \theta(h) - 2\theta_r) K(h) dh}$. The maximum sorptivity, $S(-\infty, 0)$, encompasses the two types of problems, i.e., infinite negative initial water pressure head and infinite diffusivity function close to water saturation ($\theta \to \theta_s$), and also the omission of the saturated part of sorptivity by regular approaches when $h_1 > h_a$. We considered three commonly used hydraulic models, for which Lassabatere et al. (2021) proposed analytical formulations for $S(-\infty, 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_K(h_0, h_1 = 0)$ with the condition of null water pressure head at the surface, i.e., $h_1 = 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_c \in [h_0, h_1]$, smaller than the air-entry pressure, $h_c < h_a \le 0$, and we split the integral into two separate parts as follows:

$$ 
\begin{aligned}
S_M(h_0, h_1) &= S_K(h_0, h_1) \\
&= \sqrt{\int_{h_0}^{h_1} (\theta(h_1) + \theta(h) - 2\theta(h_0)) K(h) dh} \\
&= \sqrt{\int_{h_0}^{h_c} (\theta(h_1) + \theta(h) - 2\theta(h_0)) K(h) dh + \int_{h_c}^{h_1} (\theta(h_1) + \theta(h) - 2\theta(h_0)) K(h) dh} \\
&= \sqrt{\int_{\theta(h_0)}^{\theta(h_c)} (\theta(h_1) + \theta(h) - 2\theta(h_0)) D(\theta) d\theta + \int_{h_c}^{h_1} (\theta(h_1) + \theta(h) - 2\theta(h_0)) K(h) dh}. \quad (7)
\end{aligned}
$$

In Eq. (7), the integral $\int_{h_0}^{h_c} (\theta(h_1) + \theta(h) - 2\theta(h_0)) K(h) dh$ is transformed into $\int_{\theta(h_0)} (\theta(h_1) + \theta(h) - 2\theta(h_0)) D(\theta) d\theta$ thanks to the change of variable $h \to \theta$. This operation requires that the function $\theta(h)$ is bijective over the whole interval $[h_0, h_c]$, which is valid so long as $h_c < h_a$. Alternatively, the mixed formulation, $S_M$, may be written alternatively as follows:

$$ S_M(h_0, h_1) = \sqrt{\underbrace{\int_{\theta_0}^{\theta_c} (\theta_1 + \theta - 2\theta_0) D(\theta) d\theta + \int_{h_c}^{h_1} (\theta_1 + \theta(h) - 2\theta_0) K(h) dh}_{A}} \underbrace{\int_{h_c}^{h_1} (\theta_1 + \theta(h) - 2\theta_0) K(h) dh}_{B}, \quad (8) $$

where $\theta_c = \theta(h_c)$, $\theta_1 = \theta(h_1)$, and $\theta_0 = \theta(h_0)$. 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(\theta_s) = +\infty$, since $\theta_c < \theta_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.

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:

$$ 
\left\{ 
\begin{aligned}
S_e &= \frac{\theta - \theta_r}{\theta_c - \theta_r} \\
h^* &= \frac{h}{|h_g|} \\
K_r &= \frac{K}{K_s}
\end{aligned}
\right. , \quad (9) $$

where $S_e$ is the saturation degree, $h^*$ is the scaled water pressure head, $K_r$ is the relative hydraulic conductivity, $\theta_r$ and $\theta_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):

$$ S = \sqrt{|h_g| K_s (\theta_s - \theta_r) S^*}, \quad (10) $$

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:

$$ 
\left\{ 
\begin{aligned}
S_M(h_0, h_1) &= \sqrt{|h_g| K_s (\theta_s - \theta_r) S_M^* (h_0^*, h_1^*)} \\
S_M^*(h_0^*, h_1^*) &= \left( \underbrace{\int_{S_c(h_0^*)} (S_c(h_1^*) + S_c - 2S_c(h_0^*)) D^*(S_c) dS_c}_{\theta_c} + \int_{h_c}^{h_1} (S_c(h_1^*) + S_c(h^*) - 2S_c(h_0^*)) K_r(h^*) dh^* \right)^{1/2}, \quad (11)
\end{aligned}
\right. $$

where $S_M^*(h_0^*, h_1^*)$ is the scaled version of the proposed mixed formulation, $S_M(h_0, h_1)$, with $h_0^* = h_0/|h_g|$ and $h_1^* = h_1/|h_g|$. Equation (11) can be demonstrated by changing the integration variable $\theta \to 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_0^* = -\infty$, we replace the input $h_0$ (that may be infinite) with $S_{e,0} = S_e(h_0^*)$, which always remains bounded in the final expression of the mixed formulation $S_M^*$:

$$ S_M^*(S_{e,0}, h_1^*) = \sqrt{\underbrace{\int_{S_{e,0}}^{S_c(h_0^*)} (S_c(h_1^*) + S_c - 2S_{e,0}) D^*(S_c) dS_c}_{\theta_c} + \int_{h_c}^{h_1} (S_c(h_1^*) + S_c(h^*) - 2S_{c,0}) K_r(h^*) dh^*} \quad (12) $$

Several options exist for the choice of the intermediate water pressure head $h_c^*$ and intermediate saturation degree $S_{e,c} = S_e(h_c^*)$. 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_{e,c} = (S_{e,0} + S_{e,1})/2$.
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Concept of the mixed formulation, $S_M(h_0, h_1)$: the integration of $\int_{h_0}^{h_1} (\theta_1 + \theta(h) - 2\theta_0) K(h) dh (= S_M^2(h_0, h_1))$ is converted into the sum of the integration of two bounded functions over bounded intervals, $\int_{\theta_0}^{\theta_c} (\theta_1 + \theta - 2\theta_0) D(\theta) d\theta$ and $\int_{h_c}^{h_1} (\theta_1 + \theta(h) - 2\theta_0) K(h) dh$. Note that the data are depicted on a log scale for clarity, but the integration is performed directly on the integrands instead of their log-scaled counterparts. Consequently, the integrals do not correspond to the areas below the curves, conversely to the case of linear scales. An illustrative case of the computation of sorptivity, $S_M(h_0, h_1)$, for a synthetic loamy soil with an initial water pressure head of $-10^3$ mm and a final water pressure head of $-10^{-1}$ mm.</figcaption></figure>

However, under certain circumstances (e.g., for soils with gradual water retention functions; see Results section), the value of $h^*(S_{e,c})$ may reach very large values, leading to numerical instabilities. Therefore, we use the following criteria to ensure that $h^*(S_{e,c})$ remains finite:

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

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

In the following section, the mixed formulation $S_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 $\theta_1 = \theta_s$ (or, $S_{e,1} = 1$) for a null water pressure head at the surface, $h_1 = 0$, and for very dry initial conditions $\theta_0 \to \theta_r$ (or, $h_0^* \to -\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_0 \to -\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_0$. Such an approach was used by Di Prima et al. (2020) to estimate $S_K(h_0, h_1)$ for very dry soils. In this case, the maximum sorptivity, $S(-\infty, h_1)$, is approached as follows:

$$ \begin{aligned} S(-\infty, h_1) &= \sqrt{\int_{-\infty}^{h_1} (\theta_1 + \theta(h) - 2\theta_r) K(h) dh} \\ &= \lim_{h_0 \to -\infty} \sqrt{\int_{\underline{h_0}}^{h_1} (\theta_1 + \theta(h) - 2\theta(\underline{h_0})) K(h) dh} \\ &= \lim_{h_0 \to -\infty} S_K(\underline{h_0}, h_1). \end{aligned} \quad (14) $$

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(x)} = \lim \sqrt{f(x)}$, with $x \to \sqrt{x}$ being a continuous function. In practical applications of this method, $S_K(h_0, h_1)$
<figure><img src="image_4.png" /><figcaption>**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(-\infty, h_1)$ using either $S_K(h_0, h_1)$, Eq. (5), (a) or $S_{K-V2}(h_0, h_1)$, Eq. (16), (b). The integration proceeds from the given value of $h_0$ to the right, $h \ge h_0$, corresponding to the solid zones, and the value of $h_0$ is lowered step by step to reach the limit $S(-\infty, h_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_0, h_1)$ or $S_{K-V2}(h_0, h_1)$ for a synthetic loamy soil with a final water pressure head of $h_1 = -1$ mm; successive values are considered for the initial water pressure heads, $h_0 \in \{-10^1, -10^2, -10^3\}$ mm.</figcaption></figure>

is computed for decreasing values of $h_0$ 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(-\infty, h_1)$. 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_0$.

We propose an alternative procedure that employs a specific integrand to compute the same limit (Fig. 2b). In this case, the water content $\theta_0$ is set equal to $\theta_r$ in the integrand so as to correspond to the targeted initial conditions $\theta(h_0 = -\infty) = \theta_r$:

$$ 
\begin{aligned}
S(-\infty, h_1) &= \sqrt{\int_{-\infty}^{h_1} (\theta_1 + \theta(h) - 2\theta_r) K(h) dh} \\
&= \lim_{h_0 \to -\infty} \sqrt{\int_{h_0}^{h_1} (\theta_1 + \theta(h) - 2\theta_r) K(h) dh} \\
&= \lim_{h_0 \to -\infty} S_{K-V2}(h_0, h_1) \qquad (15)
\end{aligned}
$$

with the specific function $S_{K-V2}$ defined as follows:

$$ S_{K-V2}(h_0, h_1) = \sqrt{\int_{h_0}^{h_1} (\theta_1 + \theta(h) - 2\theta_r) K(h) dh}. \qquad (16) $$

In comparison to $S_K$, the water content $\theta_0$ is replaced with $\theta_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., $(\theta_1 + \theta(h) - 2\theta(h_0)) K(h) \neq (\theta_1 + \theta(h) - 2\theta_r) K(h)$, thus involving an additional source of error (Fig. 2a). Briefly, $S_K$ combines the error due to the integral bound $h_0 > -\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(-\infty, h_1)$ and never used for the computation of other cases, i.e., $S(h_0 \neq -\infty, h_1)$, since the water content in the integrand is fixed at $\theta_r$, which corresponds exclusively to the case of $h_0 = -\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{aligned}
S^*(-\infty, h_1^*) &= \lim_{h_0^* \to -\infty} \sqrt{\int_{h_0^*}^{h_1^*} (S_{e,1} + S_e(h^*) - 2S_e(h_0^*)) K_r(h^*) dh^*} \\
&= \lim_{h_0^* \to -\infty} S_K^*(h_0^*, h_1^*) \qquad (17)
\end{aligned}
$$
$$ 
\begin{aligned}
S^*(-\infty, h_1^*) &= \lim_{h_0^* \to -\infty} \sqrt{\int_{h_0^*}^{h_1^*} (S_{e,1} + S_e(h^*)) K_r(h^*) dh^*} \\
&= \lim_{h_0^* \to -\infty} S_{K-V2}^*(\underline{h_0^*, h_1^*})
\end{aligned} 
\quad (18) 
$$

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

$$ 
\left\{ 
\begin{aligned}
S_K^*(h_0^*, h_1^*) &= \sqrt{\int_{h_0^*}^{h_1^*} (S_{e,1} + S_e(h^*) - 2S_e(h_0^*)) K_r(h^*) dh^*} \\
S_{K-V2}^*(h_0^*, h_1^*) &= \sqrt{\int_{h_0^*}^{h_1^*} (S_{e,1} + S_e(h^*)) K_r(h^*) dh^*}
\end{aligned}
\right. 
\quad (19) 
$$

The derivation of these equations involved scaling Eqs. (9) and (10) along with the change of variable $\theta \to S_e$.

### 2.2.2 Computing sorptivity with $S_D$ for null water pressure head at surface, $h_1 = 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(\theta_0, \theta_1)$ with $\theta_1 \to \theta_s$ as suggested, for instance, by Fernández-Gálvez et al. (2019):

$$ 
\begin{aligned}
S_u(\theta_0, \theta_s) &= \sqrt{\int_{\theta_0}^{\theta_s} (\theta_s + \theta - 2\theta_0) D(\theta) d\theta} \\
&= \lim_{\theta_1 \to \theta_s} \sqrt{\int_{\theta_0}^{\theta_1} (\underline{\theta_1} + \theta - 2\theta_0) D(\theta) d\theta} \\
&= \lim_{\theta_1 \to \theta_s} S_D(\theta_0, \underline{\theta_1}).
\end{aligned} 
\quad (20) 
$$

Note that with this method, we can only account for the unsaturated part of sorptivity, $S_u(\theta_0, \theta_s) = \sqrt{\int_{\theta_0}^{\theta_s} (\theta_s + \theta - 2\theta_0) D(\theta) d\theta}$, and we systematically miss the saturated portion of sorptivity $2(\theta_s - \theta_0) K_s(h_1 - h_a)$, as mentioned in Sect. 1. The total sorptivity corresponds to the sum of its two components (see Eq. 6): $S(h_0, h_1 > h_a) = \sqrt{S_u^2(\theta_0, \theta_s) + 2(\theta_s - \theta_0) K_s(h_1 - h_a)}$.

The subscript “u” in $S_u(\theta_0, \theta_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 $(\theta_1 + \theta(h) - 2\theta_0) D(\theta)$. 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 $[\theta_0, \theta_1]$ instead of $[\theta_0, \theta_s]$ (Fig. 3a, $S_D$). To correct this problem, we define a different estimator, $S_{D-V2}$, to integrate directly the targeted integrand,

($\theta_s + \theta(h) - 2\theta_0$) $D(\theta)$ (Fig. 3b, $S_{D-V2}$), and the following developments emerge:

$$ 
\begin{aligned}
S_u(\theta_0, \theta_s) &= \sqrt{\int_{\theta_0}^{\theta_s} (\theta_s + \theta - 2\theta_0) D(\theta) d\theta} \\
&= \lim_{\theta_1 \to \theta_s} \sqrt{\int_{\theta_0}^{\theta_1} (\theta_s + \theta - 2\theta_0) D(\theta) d\theta} \\
&= \lim_{\theta_1 \to \theta_s} S_{D-V2}(\theta_0, \underline{\theta_1})
\end{aligned} 
\quad (21) 
$$

with the function $S_{D-V2}$ defined as follows:

$$ S_{D-V2}(\theta_0, \theta_1) = \sqrt{\int_{\theta_0}^{\theta_1} (\theta_s + \theta - 2\theta_0) D(\theta) d\theta}. \quad (22) $$

As mentioned above, for $S_{K-V2}$, $S_{D-V2}$ should be only used for the determination of $S_u(\theta_0, \theta_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 $\theta_1 = \theta_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{aligned}
S_u^*(S_{e,0}, 1) &= \lim_{S_{e,1} \to 1} \sqrt{\int_{S_{e,0}}^{S_{e,1}} (S_{e,1} + S_e - 2S_{e,0}) D^*(S_e) dS_e} \\
&= \lim_{S_{e,1} \to 1} S_D^*(S_{e,0}, \underline{S_{e,1}})
\end{aligned} 
\quad (23) 
$$

$$ 
\begin{aligned}
S_u^*(S_{e,0}, 1) &= \lim_{S_{e,1} \to 1} \sqrt{\int_{S_{e,0}}^{S_{e,1}} (1 + S_e - 2S_{e,0}) D^*(S_e) dS_e} \\
&= \lim_{S_{e,1} \to 1} S_{D-V2}^*(S_{e,0}, \underline{S_{e,1}}).
\end{aligned} 
\quad (24) 
$$

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

$$ 
\left\{ 
\begin{aligned}
S_D^*(S_{e,0}, S_{e,1}) &= \sqrt{\int_{S_{e,0}}^{S_{e,1}} (S_{e,1} + S_e - 2S_{e,0}) D^*(S_e) dS_e} \\
S_{D-V2}^*(S_{e,0}, S_{e,1}) &= \sqrt{\int_{S_{e,0}}^{S_{e,1}} (1 + S_e - 2S_{e,0}) D^*(S_e) dS_e}
\end{aligned}
\right. 
\quad (25) 
$$

In the following sections, we compare these previously used strategies, based on the use of $S_K^*$, and $S_D^*$, and the improved versions designed for the purpose of this study, $S_{K-V2}^*$ and $S_{D-V2}^*$ with the proposed mixed formulation $S_M^*$, by quantifying the accuracy and efficiency of each.
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Illustration of the regular strategies for estimating the limits for the case of saturation $\theta_1 \to \theta_s$, with the estimation of $S_u(\theta_0, \theta_s)$ using either $S_D(\theta_0, \theta_1)$, Eq. (4), (a) or $S_{D-V2}(\theta_0, \theta_1)$, Eq. (22), (b). The integration proceeds from the given value of $\theta_1$ to the left, $\theta \ge \theta_1$, defining the solid zones. The values of $\theta_1$ are increased step by step to reach $S_u(\theta_0, \theta_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(\theta_0, \theta_1)$ or $S_{D-V2}(\theta_0, \theta_1)$ for a synthetic loamy soil with an initial water content of $\theta_0 = 0.25$; successive values are considered for the final water contents, $\theta_1 \in \{0.39, 0.40, 0.41\}$.</figcaption></figure>

## 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:

$$ \left\{ \begin{array}{l} \theta_{BC}(h) = \begin{cases} \theta_s & h \ge h_{BC} \\ \theta_r + (\theta_s - \theta_r) \left(\frac{h_{BC}}{h}\right)^{\lambda_{BC}} & h < h_{BC} \end{cases} \\ K_{BC}(\theta) = K_s \left(\frac{\theta - \theta_r}{\theta_s - \theta_r}\right)^{\eta_{BC}} \end{array} \right. . \quad (26) $$

- The van Genuchten–Burdine (vGB) model combines the van Genuchten (1980) model with the Burdine condition ($m = 1 - 2/n$) 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:

$$ \left\{ \begin{array}{l} \theta_{vGB}(h) = \theta_r + (\theta_s - \theta_r) \left(1 + \left(\frac{h}{h_{vGB}}\right)^{n_{vGB}}\right)^{-m_{vGB}} \\ m_{vGB} = 1 - \frac{2}{n_{vGB}} \\ K_{vGB}(\theta) = K_s \left(\frac{\theta - \theta_r}{\theta_s - \theta_r}\right)^{\eta_{vGB}} \end{array} \right. . \quad (27) $$

- The van Genuchten–Mualem (vGM) model combines the van Genuchten (1980) model with Mualem's condition ($m = 1 - 1/n$) 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:

$$ \left\{ \begin{array}{l} \theta_{vGM}(h) = \theta_r + (\theta_s - \theta_r) \left(1 + \left(\frac{h}{h_{vGM}}\right)^{n_{vGM}}\right)^{-m_{vGM}} \\ m_{vGM} = 1 - \frac{1}{n_{vGM}} \\ K_{vGM}(\theta) = K_s \left(\frac{\theta - \theta_r}{\theta_s - \theta_r}\right)^{l_{vGM}} \left(1 - \left(1 - \left(\frac{\theta - \theta_r}{\theta_s - \theta_r}\right)^{\frac{1}{n_{vGM}}}\right)^{m_{vGM}}\right)^2 \end{array} \right. . \quad (28) $$

- 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 charac-
terization of soils (Fernández-Gálvez et al., 2019). The *KG model* reads as follows:

$$ egin{cases} 
\theta_{\text{KG}}(h) = \theta_r + \left(\frac{\theta_r - \theta_l}{2}\right) \text{erfc} \left(\frac{\ln\left(\frac{h}{h_{\text{KG}}}\right)}{\sqrt{2}\sigma_{\text{KG}}}\right) \\ 
K_{\text{KG}}(\theta) = K_s \left(\frac{\theta_r - \theta_l}{\theta_r - \theta_l}\right)^{l_{\text{KG}}} \left(\frac{1}{2} \text{erfc} \left(\text{erfc}^{-1}\left(2\left(\frac{\theta_r - \theta_l}{\theta_r - \theta_l}\right) + \frac{\sigma_{\text{KG}}}{\sqrt{2}}\right)\right)\right)^2. 
\end{cases} 
\quad (29) $$

where erfc stands for the complementary error function.

These models involve the following common scale hydraulic parameters: residual water content, $\theta_r$; saturated water content, $\theta_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 $\lambda_{BC}$ and $\eta_{BC}$ for the BC model; $m_{vGB}$, $n_{vGB}$, and $\eta_{vGB}$ for the vGB model; $m_{vGM}$, $n_{vGM}$, and $l_{vGM}$ for the vGM model; and, lastly, $\sigma_{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):

$$ \eta = \frac{2}{\lambda} + 2 + p, \quad (30) $$

where $\lambda = \lambda_{BC}$ for the BC model and $\lambda = mn$ for the vGB models. In addition, the shape parameters $l_{vGM}$ and $l_{KG}$ are usually fixed at $1/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:

$$ egin{cases} 
S_{e,BC}(h^*) = (1 - H(1 + h^*))|h^*|^{-\lambda_{BC}} + H(1 + h^*) \\ 
K_{r,BC}(S_e) = S_e^{\eta_{BC}}. 
\end{cases} 
\quad (31) $$

– vGB model:

$$ egin{cases} 
S_{e,vGB}(h^*) = (1 + |h^*|^{n_{vGB}})^{-m_{vGB}} \\ 
\text{with } m_{vGB} = 1 - \frac{2}{n_{vGB}} \\ 
K_{r,vGB}(S_e) = S_e^{\eta_{vGB}}. 
\end{cases} 
\quad (32) $$

– vGM model:

$$ egin{cases} 
S_{e,vGM}(h^*) = (1 + |h^*|^{n_{vGM}})^{-m_{vGM}} \\ 
\text{with } m_{vGM} = 1 - \frac{1}{n_{vGM}} \\ 
K_{r,vGM}(S_e) = S_e^{l_{vGM}} \left( 1 - \left( 1 - S_e^{\frac{1}{m_{vGM}}} \right)^{m_{vGM}} \right)^2. 
\end{cases} 
\quad (33) $$

– KG model:

$$ egin{cases} 
S_{e,KG}(h^*) = \frac{1}{2} \text{erfc} \left( \frac{\ln(|h^*|)}{\sqrt{2}\sigma_{KG}} \right) \\ 
K_{r,KG}(S_e) = S_e^{l_{KG}} \left( \frac{1}{2} \text{erfc} \left( \text{erfc}^{-1} (2S_e) + \frac{\sigma_{KG}}{\sqrt{2}} \right) \right)^2. 
\end{cases} 
\quad (34) $$

These hydraulic models have the following hydraulic diffusivity functions, $D^*(S_e) = K(S_e) \frac{dh^*}{dS_e}$ (Lassabatere et al., 2021):

$$ egin{cases} 
D_{BC}^*(S_e) = \frac{1}{\lambda_{BC}} S_e^{\eta_{BC}} \left( \frac{1}{\lambda_{BC}} + 1 \right) \\ 
D_{vGB}^*(S_e) = \frac{1 - m_{vGB}}{2m_{vGB}} S_e^{\eta_{vGB}} \left( 1 - S_e^{\frac{1}{m_{vGB}}} \right)^{-\frac{1+m_{vGB}}{2}} \\ 
D_{vGM}^*(S_e) = \frac{1 - m_{vGM}}{m_{vGM}} S_e^{l_{vGM}} \left( \left( 1 - S_e^{\frac{1}{m_{vGM}}} \right)^{-m_{vGM}} + \left( 1 - S_e^{\frac{1}{m_{vGM}}} \right)^{m_{vGM}} - 2 \right) \\ 
D_{KG}^*(S_e) = \frac{1}{2} \sqrt{\frac{\pi}{2}} \sigma_{KG} S_e^{l_{KG}} \left( \text{erfc} \left( \text{erfc}^{-1} (2S_e) + \frac{\sigma_{KG}}{\sqrt{2}} \right) \right)^2 e^{\left( \text{erfc}^{-1} (2S_e) \right)^2 + \sqrt{2} \sigma_{KG} \text{erfc}^{-1} (2S_e)}. 
\end{cases} 
\quad (35) $$

These equations are needed for the computation of the dimensionless sorptivity with the proposed mixed formulation $S_M^{*2}$ and the regular formulations $S_K^{*2}$, $S_{K-v2}^{*2}$, $S_D^{*2}$, and $S_{D-v2}^{*2}$.

The studied hydraulic models exhibit contrasting and challenging features for the computation of sorptivity, including non-null water pressure heads $h_a^* < 0$ and infinite hydraulic diffusivity close to saturation $\lim_{S_e \to 1} D^*(S_e) = +\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_e (h^*) = 1/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 $|h^*|$ (see Fig. 2, in Lassabatere et al., 2021, and also Sect. 3). Conversely, when $x$ gets close to unity, WR functions approach the stepwise function with an abrupt decrease of $S_e$ with the increase of $|h^*|$. Lassabatere et al. (2021) defined the WR shape index $x$ as follows:

$$ egin{cases} 
x_{BC} = \frac{\lambda_{BC}}{2+\lambda_{BC}} \\ 
x_{vGB} = m_{vGB} \\ 
x_{vGM} = m_{vGM} \\ 
x_{KG} = \frac{1}{1+\sigma_{KG}} 
\end{cases} 
\quad (36) $$

Lassabatere et al. (2021) also analytically determined the maximum squared scaled sorptivity $S^{*2}(-\infty, 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:
$$
\left\{
\begin{aligned}
c_{\text{p,BC}}(x) &= 2 + \frac{1-x}{5x+1} + \frac{1-x}{7x+1} \\
c_{\text{p,vGB}}(x) &= \Gamma\left(\frac{3-x}{2}\right) \left[ \frac{\Gamma\left(\frac{1+5x}{2}\right)}{\Gamma(1+2x)} + \frac{\Gamma\left(\frac{1+7x}{2}\right)}{\Gamma(1+3x)} \right] \\
c_{\text{p,vGM}}(x) &= \Gamma(2-x) \left[ \frac{\Gamma\left(\frac{3}{2}x\right)}{\left(\frac{3}{2}x-1\right)\Gamma\left(\frac{1}{2}x\right)} + \frac{\Gamma\left(\frac{5}{2}x\right)}{\left(\frac{5}{2}x-1\right)\Gamma\left(\frac{3}{2}x\right)} \right] \\
&\quad + (1-x) \left[ \frac{\Gamma\left(\frac{3}{2}x\right)\Gamma(1+x)}{\left(\frac{3}{2}x-1\right)\Gamma\left(\frac{5}{2}x\right)} + \frac{\Gamma\left(\frac{5}{2}x\right)\Gamma(1+x)}{\left(\frac{5}{2}x-1\right)\Gamma\left(\frac{7}{2}x\right)} \right] \\
&\quad - 2 \left( \frac{1}{\frac{3}{2}x-1} + \frac{1}{\frac{5}{2}x-1} \right)
\end{aligned}
\right.
\quad (37)
$$

Note that no analytical formulation can be found for the case of Kosugi's hydraulic model, so $c_{\text{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_{\text{vGM}} = l_{\text{KG}} = \frac{1}{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 $\hat{S}$ of the nominal dimensional sorptivity $S$, and the related scaled variables, $\hat{S}^*$ and $S^*$, the following relations emerge between the relative errors, $E_r$, related to dimensional and scaled sorptivity:

$$
\begin{align}
E_r(S) &= \frac{\hat{S} - S}{S} \\
&= \frac{\hat{S}^* \cdot \sqrt{|h_g| K_s (\theta_s - \theta_r)} - S^* \cdot \sqrt{|h_g| K_s (\theta_s - \theta_r)}}{S^* \cdot \sqrt{|h_g| K_s (\theta_s - \theta_r)}} \\
&= \frac{\hat{S}^* - S^*}{S^*} \\
&= E_r(S^*), \tag{38}
\end{align}
$$

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^*(-\infty, 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_0 = -\infty$, and saturated final conditions with null water pressure head, $h_1 = 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 \in \{0.01, 0.02, \dots, 0.99\}$. For the second goal of the study, point (ii), we compared the values provided by the proposed mixed formulation $S_M^*$ with the nominal (error-free) values of sorptivity, $S^* = S^*(-\infty, 0) = \sqrt{c_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_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_M^*$, to evaluate the efficiency of those previously used strategies. We considered several scenarios for the use of $S_K^*$, $S_{K-V2}^*$, $S_D^*$, and $S_{D-V2}^*$, with several values of the lower water pressure head, $h_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_0, h_1)$, 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_f = -150$ or $h_f = 30$ mm) along with an initial water pressure head of $-10$ m ($h_i = -10^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 $|h^*| \le 1$, i.e., $h^* \ge -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_M^*$ compared to the usual use of $S_D^*$ for the cases of problematic computation close to saturation (i.e., $h_1 \to 0$ and $\theta_1 \to \theta_s$). Similarly, the accuracy of the commonly used $S_K^*$ version can be challenged when integrating over infinite intervals $(-\infty, 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_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_D^*$ for large values of $x$, with quasi-infinite values for the hydraulic diffusivity close to saturation, i.e., $S_e \to 1$ (Fig. 4h, l, and p).

## 3.2 Validation of the proposed mixed formulation, $S_M^*$, against the nominal sorptivity, $S^*(-\infty, 0)$

The computations using $S_M^*$ (Eqs. 12 and 13) with $h_c^* = h^*\left(\frac{S_e,0+S_e,1}{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_{vGB} \in \{0.01, 0.02\}$), and the first 17 values for the vGM model ($x_{vGM} \in \{0.01, 0.02, \dots, 0.17\}$). The value of $z=0$ was enough to allow the computation in all cases, apart from the case of the vGM model for which the value of $z=-1$ had to be considered for $x_{vGM} \in \{0.02, \dots, 0.07\}$ and

<table><caption>**Table 1.** Absolute values of relative errors, $|E_r|$, between the proposed mixed formulation, $S_M^*$, and the targeted scaled sorptivity, $S^*(-\infty, 0) = \sqrt{c_p}$, with the mean $|E_r|$, the standard deviation ($\sigma_{|E_r|}$), the minimum and the maximum values for the three hydraulic models whose sorptivity is analytically tractable: Brooks and Corey (BC), van Genuchten-Burdine (vGB), van Genuchten-Mualem (vGM). Note that $10^{-16}$ corresponds to the relative precision of numbers in Scilab.</caption><thead><tr><th>$|E_r|$</th><th>BC</th><th>vGB</th><th>vGM</th></tr></thead><tbody><tr><td>$|E_r|$</td><td>1.445 &times; 10<sup>-13</sup></td><td>5.594 &times; 10<sup>-13</sup></td><td>3.309 &times; 10<sup>-9</sup></td></tr><tr><td>$\sigma_{|E_r|}$</td><td>2.774 &times; 10<sup>-13</sup></td><td>1.014 &times; 10<sup>-12</sup></td><td>1.934 &times; 10<sup>-8</sup></td></tr><tr><td>Min</td><td>&lt; 10<sup>-16</sup></td><td>&lt; 10<sup>-16</sup></td><td>8.255 &times; 10<sup>-16</sup></td></tr><tr><td>Max</td><td>1.201 &times; 10<sup>-12</sup></td><td>6.037 &times; 10<sup>-12</sup></td><td>2.000 &times; 10<sup>-7</sup></td></tr></tbody></table>

$z = -2$ for $x_{vGM} = 0.01$. Note that, as long as convergence is obtained, the values of $S_M^*(-\infty, 0)$ do not depend on $z$. With this strategy involving a threshold, the proposed mixed formulation, $S_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_M^*(-\infty, 0)$, and the targeted sorptivity, $S^*(-\infty, 0) = \sqrt{c_p}$, was analyzed in terms of means, standard deviations, and minimum and maximum values (Table 1):

$$E_r = \frac{S_M^*(-\infty, 0) - \sqrt{c_p}}{\sqrt{c_p}}. \quad (39)$$

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

In this section, we compare the estimates provided by the strategies commonly considered (i.e., $S_D^*$, $S_{D-V2}^*$, $S_K^*$, and $S_{K-V2}^*$) with the reference values of the targeted scaled sorptivity $S^* = S^*(-\infty, 0)$. For these comparisons, we consider the analytical formulations of Eq. (37) for the BC, vGB, and vGM
<figure><img src="image_6.png" /><figcaption>**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 $\lambda_{BC}$, $m_{vGM}$, $m_{vGB}$, and $\sigma_{KG}$ were computed as a function of $x$ using Eq. (36) with $l_{vGM} = l_{KG} = \frac{1}{2}$. Adapted from Lassabatere et al. (2021).</figcaption></figure>

models and use the proposed mixed formulation $S_M^*$ for the KG model. Indeed, no analytical expressions are available for this last model, whereas the estimates provided by $S^* = S_M^*$ are very accurate for the BC, vGB, and vGM models (see Sect. 3.2), and thus $S_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_K^*$ and $S_{K\_V2}^*$ towards $S^*$. These functions, in particular, $S_K^*$, are often used without special attention regarding their accuracy. $S_K^*$ and $S_{K\_V2}^*$ converge to the limit $S^*$ when $|h_0^*|$ becomes large enough (Fig. 5, left column). For instance, for the BC model with $\lambda = 0.56$, the use of $S_K^*$ with $h_0^* = -10$ and $h_0^* = -100$ has respective relative errors

of $E_r = -4.1\%$ and $E_r = -15\%$ (Fig. 5a, dashed red line). The second estimator, $S_{K\_V2}^*$, converges much faster than $S_K^*$. With the same values of $h_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_0^* = -10$, with high accuracy. Such an improvement results from setting the initial saturation degree, $S_{e,0}$, at its target value, $S_{e,0} = 0$, in the integrand (see Eq. 18 versus Eq. 17). The same conclusions hold regardless of the hydraulic model (Fig. 5c, e, and g).

The convergence of the two functions, $S_D^*$ and $S_{D\_V2}^*$, is depicted for the case of BC functions in Fig. 5b. For both estimators, the estimates are far from the target values, $S^*$, with $|E_r|$ close to 50%. This large error results from the omission of the saturated part
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Relative errors of the regular functions $S_D^*(0, S_{e,1})$ and $S_{D-V2}^*(0, S_{e,1})$ (b, d, f, h) and $S_K^*(h_0^*, 0)$ and $S_{K-V2}^*(h_0^*, 0)$ (a, c, e, g) towards $S^*(-\infty, 0)$ for the four hydraulic models, for the following specific cases: BC model with $\lambda_{BC} = 56$ ($x_{BC} = 0.22$), vGB model with $n_{vGB} = 3$ ($x_{vGB} = 1/3$), vGM model with $n_{vGM} = 2$ ($x_{vGM} = 1/2$) with $l_{vGM} = 0.5$, and KG model with $\sigma_{KG} = 1.5$ ($x_{KG} = 0.4$) and $l_{KG} = 0.5$.</figcaption></figure>
of sorptivity, as explained above (Sect. 2.2). By design, $S_D^*$ and $S_{D-V2}^*$ converge towards the unsaturated part of sorptivity, $S_u(0, 1) = \sqrt{\int_0^1 (1 + S_e) D^*(S_e) dS_e}$, and miss the additional saturated part of the scaled sorptivity, $2(S_{e,1} - S_{e,0}) K_f(S_{e,1}) (h_1^* - h_a^*) = 2$, given that $h_1^* = 0$ ($S_{e,1} = 1$), $h_0^* = -\infty$ ($S_{e,0} = 0$), and $h_a^* = 1$. When the saturated portion is added to the estimators, the convergence becomes excellent (Fig. 5b, green line). These results show that the computation of sorptivity using the integration with regards to saturation degree, like Eqs. (4) and (25), leads to erroneous estimations. Accounting for the term $2(\theta_s - \theta_0)K_s|h_a|$ 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_D^*$ and $S_{D-V2}^*$ are significantly better (Fig. 5d, f, and h). Indeed, the air-entry water pressure head is null in these models, which removes the large underestimation due to the omission of the saturated part of sorptivity, as observed for the BC model. However, the underestimation remains substantial, with relative errors on the order of 15%–20% for $S_{e,1} = 0.99$ (Fig. 5d, f and h). Even when $S_{e,1} = 0.999$, the relative errors are larger than 10%. From these results, we conclude that determining the sorptivity, $S^*$, by setting the final saturation degree to near-unity does not provide good estimates, regardless of the estimators considered. These poor estimates result from the fact that the diffusivity and, thus, the integrand are infinite. The convergence of the integrals defined by Eq. (25) and thus Eq. (4) is slow and prevents accurate estimations. Conversely, in the case of $S_K^*$ and $S_{K-V2}^*$, changing the integrand by fixing it to the targeted integrand does not significantly improve the convergence.

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

The previous results were presented for specific values of the shape parameters in Fig. 5. One may wonder if the accuracy of estimators $S_K^*$ and $S_{K-V2}^*$, on the one hand, and $S_D^*$ and $S_{D-V2}^*$, on the other, varies with the shape parameters and indices. Figure 6 depicts the relative errors of the estimators $S_K^*$ and $S_{K-V2}^*$ as a function of the shape parameters for different values of $h_0^*$ (left column of Fig. 6) and of the estimators $S_D^*$ and $S_{D-V2}^*$ 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_D^*$ and $S_{D-V2}^*$ are obtained for intermediate values of WR shape indices (Fig. 6, right column). Discrepancies increase for both small or large values of the WR shape index. As detailed above (Sect. 3.3.1), accurate predictions require the practitioner to fix the upper integration boundary, $S_{e,1}$, to a minimum of 0.999, to get errors less than 10%. However, estimates remain poor, even with $S_{e,1} = 0.9999$ or $S_{e,1} = 0.99999$, when the WR shape indices tend towards unity (Fig. 6d, f, and h).

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

The results obtained in this study also revealed particular behaviors for the different formulations used to compute sorptivity, specifically $S_K^*$, $S_{K-V2}^*$, $S_D^*$, and $S_{D-V2}^*$. The approaches that compute sorptivity by integration with regards to $h^*$ (i.e., $S_K^*$ and $S_{K-V2}^*$) have errors that come from the choice of the lower integration boundary $h_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_D^*$ and $S_{D-V2}^*$) have more error associated with the omission of the saturated part, as well as slow convergence of the integration procedure when dealing with infinite diffusivity values. Among the usual procedures, the $S_K^*$ estimator is better than $S_D^*$, and the $S_{K-V2}^*$ estimate is superior to both. Indeed, $S_{K-V2}^*$ does not miss the saturated part of sorptivity, like $S_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 $|E_r| < 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_M^*$. Consequently, the use of the mixed formulation, $S_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
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 6.** Relative errors of estimators $S_D^*$ and $S_{D-V2}^*$ (**b**, **d**, **f**, **h**) and $S_K^*$ and $S_{K-V2}^*$ (**a**, **c**, **e**, **g**) for the four selected hydraulic models, i.e., BC, vGB, vGM, and KG models, as a function of the WR shape indices.</figcaption></figure>
<table><caption>**Table 2.** Hydraulic parameters, i.e., values of the residual and saturated water contents, $\theta_r$ and $\theta_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).</caption><thead><tr><th></th><th>vGM</th><th>BC</th><th>KG</th></tr></thead><tbody><tr><td>θ<sub>r</sub></td><td>0.078</td><td>0.078</td><td>0.078</td></tr><tr><td>θ<sub>s</sub></td><td>0.43</td><td>0.43</td><td>0.43</td></tr><tr><td>h<sub>g</sub> [mm]</td><td>-277.8</td><td>-111.5</td><td>-1018</td></tr><tr><td>h<sub>a</sub> [mm]</td><td>0</td><td>-111.5</td><td>0</td></tr><tr><td>K<sub>s</sub> [mm s<sup>-1</sup>]</td><td>2.288 × 10<sup>-3</sup></td><td>3.667 × 10<sup>-3</sup></td><td>2.288 × 10<sup>-3</sup></td></tr><tr><td>Shape parameter</td><td>n<sub>vGM</sub> = 1.56</td><td>λ<sub>BC</sub> = 0.34</td><td>σ<sub>KG</sub> = 1.997</td></tr><tr><td>Tortuosity parameter</td><td>l<sub>vGM</sub> = 1/2</td><td>p<sub>BC</sub> = 1</td><td>l<sub>KG</sub> = 1/2</td></tr></tbody></table>

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

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_f = -150$ mm and an initial water pressure head of $h_i = -10^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 di-

mensional 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_{e,i} = S_e(h_i)$ and $S_{e,f} = S_e(h_f)$.

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

3. Compute the intermediate saturation degree from the intermediate water pressure head, $S_{e,c} = S_e(h_c^*)$.

4. Integrate the lower part of the squared sorptivity

$$A = \int_{S_{e,i}}^{S_{e,c}} (S_{e,f} + S_e - 2S_{e,i}) D^*(S_e) dS_e.$$

5. Integrate the upper part of the squared sorptivity

$$B = \int_{h_c^*}^{h_f^*} (S_{e,f} + S_e(h^*) - 2S_{e,i}) K_r(h^*) dh^*.$$

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

7. Upscale by multiplying the scaled sorptivity,

$$S_M = \sqrt{|h_g| K_s (\theta_s - \theta_r) S_M^*}.$$

For the studied case, with $h_i = -10$ m and $h_f = -150$ mm, it leads to the following results:

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

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

3. The corresponding saturation degree takes the value of $S_{e,c}^* = 0.780$.

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

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

6. The combination of the two parts gives the scaled sorptivity: $S_M^* = 0.2931$.

7. The upscaling finalizes the computation of sorptivity: $S_M = 0.156 \text{ mm s}^{-1/2}$.
<figure><img src="image_16.png" /><figcaption>**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$.</figcaption></figure>

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_f = -150$ mm (Fig. 7c) and as a function of final water pressure head, $h_f$, for an initial water pressure head of $h_i = -10^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 \ge 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

<figure><img src="image_17.png" /></figure>

with $h_f$ but remain equal to the unsaturated part of sorptivity: $S_D(\theta_i, (h_f \ge h_a)) = S_D(\theta_i, \theta_s) = S_u(\theta_i, \theta_s)$ (Fig. 7d, plateaus with dashed lines). Note for vGM and KG models, the air-entry water pressure head is zero, then the plateaus begin at $h_f = 0$, whereas the air-entry water pressure head is lower than zero for the BC model, $h_a = -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 con-
ductivity 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_i = -10^4$ mm). To model water infiltration, we used the analytical model developed for $h_f \le 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:

$$ \frac{2\Delta K^2}{S^2} t = \frac{1}{1-\beta} \left( \frac{2\Delta K}{S^2} (I_{1D} - K_i t) - \ln \left( \frac{e^{\frac{2\beta\Delta K}{S^2}} (I_{1D} - K_i t) + \beta - 1}{\beta} \right) \right), \quad (40) $$

where $\Delta K$ stands for the difference between final and initial values of hydraulic conductivity $\Delta K = K_f - K_i$, and $\beta$ 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$:

$$ \left\{ \begin{array}{ll} \gamma_t = \frac{S^2}{2\Delta K} & I = \gamma_t I^* + K_i t \\ \gamma_t = \frac{S^2}{2\Delta K^2} & t = \gamma_t t^* \end{array} \right. . \quad (41) $$

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

$$ I_{\text{QEI}}(t) = \gamma_t I_{\text{QEI}}^* \left(\frac{t}{\gamma_t}\right) + K_i t \quad (42) $$

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

$$ t^* = \frac{1}{1-\beta} \left( I^* - \ln \left( \frac{e^{\beta I^*} + \beta - 1}{\beta} \right) \right). \quad (43) $$

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

$$ I^* = \frac{\sigma}{q^* - 1} + \frac{1 - \sigma}{\beta} \ln \left( 1 + \frac{\beta}{q^* - 1} \right) \quad (44) $$

$$ \begin{aligned} t^* = & \frac{1 - \sigma}{\beta(1 - \beta)} \ln \left( 1 + \frac{\beta}{q^* - 1} \right) + \frac{\sigma}{q^* - 1} \\ & - \frac{1 - \sigma\beta}{1 - \beta} \ln \left( 1 + \frac{1}{q^* - 1} \right), \end{aligned} \quad (45) $$

where $q^* = \frac{q}{\gamma_q}$ is the scaled infiltration rate and $\gamma_q = \Delta 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^*(q^*)$, defined by Eq. (45), and the second relating the cumulative infiltration and the infiltration rate, $I^* = I^*(q^*)$, defined by Eq. (44). To compute the cumulative infiltration as a function of time, the infiltration rate $q^*$ must be retrieved as the root of the Eq. (45) before being inserted into Eq. (44) to define the scaled infiltration as a function of the scaled time, $I_{\text{QEI\_ext}}^*(t^*)$. In the following, this relationship is referred to as the extended version of the QEI model and denoted as the QEI\_ext model. This set of scaled models, Eqs. (44) and (45), must be upscaled considering the same scaling factors, as defined in Eq. (41) with the additional parameter, $\sigma$, which quantifies the relative contribution of the saturated part of sorptivity to the total sorptivity:

$$ \sigma = \frac{2(\theta_s - \theta_i) K_s (h_f - h_a)}{S^2}. \quad (46) $$

The QEI\_ext model addresses all cases. When $h_f \le h_a$, the saturated part of sorptivity is null, leading to $\sigma = 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.

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
<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /><figcaption>**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 \le h_a$ (BC model, $h_i = -10$ m and $h_f = -150$ mm) and (d) the case of saturated conditions, $h_f \ge h_a$ (BC model, $h_i = -10$ m and $h_f = 30$ mm).</figcaption></figure>

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 cu-

mulative 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_i = -10^4$ mm and a final water pressure head of $h_f = -150$ mm (Fig. 8c) and for the case of the same initial water pressure head, $h_i = -10^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 ($< 10^{-7}$). Conversely, the use of regular estimates for sorptivity leads to large underestimates in many instances, in particular when the final water pressure head exceeds the air-entry water pressure head.

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

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

**Data availability.** No data sets were used in this article.

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

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

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

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

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

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

## References

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

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

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

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

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

*Brooks, R. and Corey, A.: Hydraulic Properties of Porous Media, Hydrology Papers, Colorado State University, 1964.
Campbell, S. L., Chancelier, J. P., and Nikoukhah, R.: Modeling and Simulation in SCILAB, Modeling and Simulation in Scilab/Scicos with ScicosLab 4.4, Springer, https://doi.org/10.1007/978-1-4419-5527-2, 2010.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Lassabatere, L., Peyneau, P.-E., Yilmaz, D., Pollacco, J., Fernández-Gálvez, J., Latorre, B., Moret-Fernández, D., Di Prima, S., Rahmati, M., Stewart, R. D., Abou Najm, M., Hammeker, C., and Angulo-Jaramillo, R.: A scaling procedure for straightforward computation of sorptivity, Hydrol. Earth Syst. Sci., 25, 5083–5104, https://doi.org/10.5194/hess-25-5083-2021, 2021.

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

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

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

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

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

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

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

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

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

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

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

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

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

Varado, N., Braud, I., Ross, P. J., and Haverkamp, R.: Assessment of an efficient numerical solution of the 1D Richards' equation on bare soil, J. Hydrol., 323, 244–257, 2006.
White, I. and Sully, M. J.: Macroscopic and microscopic capillary length and time scales from field infiltration, Water Resour. Res., 23, 1514–1522, 1987.

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

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