# A scaling procedure for straightforward computation of sorptivity

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

► To cite this version:

Laurent Lassabatere, Pierre-Emmanuel Peyneau, Deniz Yilmaz, Joseph Pollacco, Jesús Fernández-Gálvez, et al.. A scaling procedure for straightforward computation of sorptivity. Hydrology and Earth System Sciences, 2021, 25 (9), pp.5083-5104. 10.5194/hess-25-5083-2021. hal-03370546

HAL Id: hal-03370546

https://hal.science/hal-03370546

Submitted on 19 Oct 2021

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
# A scaling procedure for straightforward computation of sorptivity

Laurent Lassabatere¹, Pierre-Emmanuel Peyneau², Deniz Yilmaz³, Joseph Pollacco⁴, Jesús Fernández-Gálvez⁵, Borja Latorre⁶, David Moret-Fernández⁶, Simone Di Prima⁷, Mehdi Rahmati⁸,⁹, Ryan D. Stewart¹⁰, Majdi Abou Najm¹¹, Claude Hammecker¹², and Rafael Angulo-Jaramillo¹

¹Univ. Lyon, Université Claude Bernard Lyon 1, CNRS, ENTPE, UMR5023 LEHNA, 69518 Vaulx-en-Velin, France

²GERS-LEE, Univ. Gustave Eiffel, IFSTTAR, 44344 Bouguenais, France

³Civil Engineering Department, Engineering Faculty, Munzur University, Tunceli, Turkey

⁴Manaaki Whenua - Landcare Research, 7640 Lincoln, New Zealand

⁵Department of Regional Geographic Analysis and Physical Geography, University of Granada, 18071 Granada, Spain

⁶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

⁷Agricultural Department, University of Sassari, Viale Italia, 39, 07100 Sassari, Italy

⁸Department of Soil Science and Engineering, Faculty of Agriculture, University of Maragheh, Maragheh, Iran

⁹Forschungszentrum Jülich GmbH, Institute of Bio- and Geosciences: Agrosphere (IBG-3), Jülich, Germany

¹⁰School of Plant and Environmental Sciences, Virginia Polytechnic Institute and State University, Blacksburg, VA, United States

¹¹Department of Land, Air and Water Resources, University of California, Davis, CA 95616, United States

¹²University of Montpellier, UMR LISAH, IRD, Montpellier, France

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

Received: 13 March 2021 – Discussion started: 24 March 2021

Revised: 17 July 2021 – Accepted: 3 August 2021 – Published: 22 September 2021

**Abstract.** Sorptivity is a parameter of primary importance in the study of unsaturated flow in soils. This hydraulic parameter is required to model water infiltration into vertical soil profiles. Sorptivity can be directly estimated from the soil hydraulic functions (water retention and hydraulic conductivity curves), using the integral formulation of Parlange (1975). However, calculating sorptivity in this manner requires the prior determination of the soil hydraulic diffusivity and its numerical integration between initial and final saturation degrees, which may be difficult in some situations (e.g., coarse soil with diffusivity functions that are quasi-infinite close to saturation). In this paper, we present a procedure to compute sorptivity using a scaling parameter, $c_p$, that corresponds to the sorptivity of a unit soil (i.e., unit values for all parameters and zero residual water content) that is utterly dry at the initial state and saturated at the final state. The $c_p$ parameter was computed numerically and analytically for five hydraulic models: delta (i.e., Green and Ampt), Brooks and Corey, van Genuchten-Mualem, van Genuchten-Burdine, and Kosugi. Based on the results, we proposed brand new analytical expressions

for some of the models and validated previous formulations for the other models. We also tabulated the output values so that they can easily be used to determine the actual sorptivity value for any case. At the same time, our numerical results showed that the relation between $c_p$ and the hydraulic shape parameters strongly depends on the chosen model. These results highlight the need for careful selection of the proper model for the description of the water retention and hydraulic conductivity functions when estimating sorptivity.

## 1 Introduction

Soil sorptivity represents the capacity of a soil to absorb or desorb liquid by capillarity and is therefore one of the key factors for modeling water infiltration into soil (Cook and Minasny, 2011). Knowledge of soil sorptivity is also important when deciphering soil physical properties, such as hydraulic conductivity, from infiltration tests (e.g., Lassabatere
et al., 2006). Sorptivity is incorporated in a wide range of infiltration models (Angulo-Jaramillo et al., 2016; Lassabatere et al., 2009, 2014, 2019). Sorptivity varies depending on initial and boundary conditions (e.g., water contents). However, calculated sorptivity values can also vary depending on the chosen soil hydraulic model, making it important to assess the impact of such a choice on the computation of sorptivity. In this study, we address this issue and propose a new scaling procedure to simplify its computation.

One of the first equations proposed for the computation of sorptivity was developed by Philip (1957), who modeled the 1D gravity-free water infiltration as

$$ egin{cases} I(t) = S(\theta_0, \theta_1) \sqrt{t} \\ S(\theta_0, \theta_1) = \int_{\theta_0}^{\theta_1} \chi(\theta) d\theta \end{cases} \quad (1) $$

In the above equations, $S(\theta_0, \theta_1)$ stands for the sorptivity between $\theta_0$ and $\theta_1$, $\chi(\theta) = \frac{\chi(\theta, t)}{\sqrt{t}}$ stands for the Boltzmann transformation variable, $\theta_0$ is the initial water content, $\theta_1$ is the final water content corresponding also to the water content applied at the soil surface, and $t$ is the elapsed time. In practical applications, the user must perform numerical modeling of horizontal infiltration using a given set of hydraulic functions and initial and final conditions. Then, the Boltzmann transformation must be computed by integrating the modeled water content profile. This procedure is often time-consuming and subject to numerical instabilities that lead to substantial errors.

To avoid such complexity, Parlange (1975) proposed a formulation that directly relates sorptivity to the hydraulic functions and the initial and final water contents:

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

where $D(\theta) = K(\theta) \frac{dh}{d\theta}$ is the hydraulic diffusivity function. Note that several integral expansions were proposed for the computation of sorptivity (Angulo-Jaramillo et al., 2016). This point is beyond the framework of this study and will be the subject of another study. While the above equation provides the diffusivity form for sorptivity determination, it can be equally defined as a function of the hydraulic conductivity function, $K(h) = K(\theta(h))$:

$$ egin{aligned} S_K^2(h_0, h_1) &= \int_{h_0}^{h_1} (\theta(h_1) + \theta(h) - 2\theta(h_0)) K(h) dh \\ &= \int_{h_0}^{h_1} (\theta_1 + \theta(h) - 2\theta_0) K(h) dh, \end{aligned} \quad (3) $$

where $h$ is the water pressure head, $h_0$ and $h_1$ are respectively the initial and final water pressure heads, and $\theta_0 = \theta(h_0)$

and $\theta_1 = \theta(h_1)$. For the sake of clarity, the functions $S_D^2$ and $S_K^2$ are respectively referred to as the “diffusivity” and “conductivity” forms of sorptivity. $S_D^2$ and $S_K^2$ are equivalent so long as the water retention function $\theta(h)$ is bijective over the water pressure head interval $[h_0, h_1]$, which is the case when the water pressure head at the surface is lower than the air-entry water pressure head, i.e., $h_1 \le h_a$. Then, Eq. (3) can be easily deduced from Eq. (2) by a simple change of variable from $\theta$ to $h$:

$$ egin{aligned} S_K^2(h_0, h_1 \le h_a) &= S_D^2(\theta(h_0), \theta(h_1)) \\ &= S_D^2(\theta_0, \theta_1). \end{aligned} \quad (4) $$

Otherwise, when the surface water pressure head exceeds the soil air-entry pressure, i.e., $h_1 > h_a(\le 0)$, sorptivity must be computed using Eq. (3) (Ross et al., 1996). Indeed, the function $\theta(h)$ is no longer bijective over the interval $[h_0, h_1]$. Consequently, Eqs. (2) and (3) are no longer equivalent. In this case, the integration involved in Eq. (3) must be divided into two parts, one part integrating over the interval $[h_0, h_a]$ ensuring the bijectivity of the function $\theta(h)$ and the other part retaining the integration interval $[h_a, h_1]$ corresponding to the saturated part of the integration (Ross et al., 1996). Then, the change of variable from $\theta$ to $h$ in the first integral leads to the retrieval of $S_D^2$:

$$ egin{aligned} S_K^2(h_0, h_1 > h_a) &= \int_{h_0}^{h_a} (\theta_s + \theta(h) - 2\theta_0) K(\theta(h)) dh \\ &\quad + \int_{h_a}^{h_1} (\theta_s + \theta(h) - 2\theta_0) K(\theta(h)) dh \\ &= \int_{h_0}^{h_a} (\theta_s + \theta(h) - 2\theta_0) K(\theta(h)) dh \\ &\quad + 2(\theta_s - \theta_0) K_s \int_{h_a}^{h_1} dh \\ &= \int_{h_0}^{h_a} (\theta_s + \theta(h) - 2\theta_0) K(\theta(h)) dh \\ &\quad + 2(\theta_s - \theta_0) K_s (h_1 - h_a) \\ &= \int_{\theta_0}^{\theta_s} (\theta_s + \theta - 2\theta_0) D(\theta) d\theta \\ &\quad + 2(\theta_s - \theta_0) K_s (h_1 - h_a) \\ &= S_D^2(\theta_0, \theta_s) + 2(\theta_s - \theta_0) K_s (h_1 - h_a). \end{aligned} \quad (5) $$

The computation of sorptivity with Eq. (2) or Eq. (3) requires a set of hydraulic functions to be chosen from the wide range of available models. Here, we considered five of the most widely used hydraulic models. Firstly, we considered the delta model (d, delta), which involves Dirac delta
functions (stepwise functions) for the description of both the water retention (WR) and hydraulic conductivity (HC) functions. Indeed, this model is often considered for analytical resolutions to the Richards equation and the determination of analytical expressions for water infiltration, like the Green and Ampt approach (Triadis and Broadbridge, 2012). Secondly, we considered the Brooks and Corey (1964) model, referred to as the BC model, since it is among the first hydraulic models of soil physics (Hillel, 1998). The BC model involves power functions for both the WR and HC functions, thus allowing for analytical integration of Eq. (2) and leading to analytical expressions for sorptivity (e.g., Varado et al., 2006). Thirdly, the van Genuchten–Burdine (vGB) model was studied since it has been used for the development of the BEST methods (Beerkan Estimation of Soil Transfer functions) for the characterization of soil hydraulic properties (Lassabatere et al., 2006; Yilmaz et al., 2010; Bagarello et al., 2014). The vGB model combines the van Genuchten (1980) model with Burdine's condition ($m = 1 - \frac{2}{n}$) for the WR function and the Brooks and Corey (1964) model for the HC function. Fourthly, we considered the van Genuchten–Mualem (vGM) model that combines the van Genuchten (1980) model with Mualem's condition ($m = 1 - \frac{1}{n}$) for the WR function and the Mualem (1976) capillary model for the HC function. The vGM is among the most widely used models and often used for the numerical modeling of flow in the vadose zone (Şimûnek et al., 2003). Lastly, the Kosugi (KG) model (Kosugi, 1996), was also considered since it relates the water retention function to physical characteristics of the soil pore size distribution assuming log-normal distributions.

These five models have the following mathematical expressions (Triadis and Broadbridge, 2012; Brooks and Corey, 1964; van Genuchten, 1980; Mualem, 1976; Kosugi, 1996; Nasta et al., 2013):

delta model:

$$ \left\{ \begin{array}{l} \theta_d(h) = \theta_r + (\theta_s - \theta_r) H \left(1 + \frac{h}{|h_d|}\right) \\ K_d(\theta) = K_s H \left(\frac{\theta - \theta_r}{\theta_s - \theta_r} - 1\right) \end{array} \right. \qquad (6) $$

BC model:

$$ \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. \qquad (7) $$

vGB model:

$$ \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)^{n_{vGB}} \end{array} \right. \qquad (8) $$

vGM model:

$$ \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. \qquad (9) $$

KG model:

$$ \left\{ \begin{array}{l} \theta_{KG}(h) = \theta_r + \frac{(\theta_s - \theta_r)}{2} \operatorname{erfc} \left( \frac{\ln\left(\frac{h}{h_{KG}}\right)}{\sqrt{2\sigma_{KG}}} \right) \\ K_{KG}(\theta) = K_s \left( \frac{\theta - \theta_r}{\theta_s - \theta_r} \right)^{l_{KG}} \left( \frac{1}{2} \operatorname{erfc} \left( \operatorname{erfc}^{-1} \left( 2 \frac{\theta - \theta_r}{\theta_s - \theta_r} \right) + \frac{\sigma_{KG}}{\sqrt{2}} \right) \right)^2, \end{array} \right. \qquad (10) $$

where *H* stands for the one-sided Heaviside step function: $H(x < 0) = 0$, $H(x \ge 0) = 1$ (Triadis and Broadbridge, 2012); and “erfc” stands for the complementary error function. These models involve several specific hydraulic shape parameters and 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$ (i.e., $h_d$, $h_{BC}$, $h_{vGB}$, $h_{vGM}$, or $h_{KG}$), and saturated hydraulic conductivity, $K_s$. The delta and BC models involve a non-null air-entry water pressure head, $h_d$ and $h_{BC}$, meaning that air needs a nonzero suction to enter into the soil and start to desaturate it. For the sake of simplicity, the scale parameter for the water pressure head is fixed at the air-entry pressure head, i.e., $h_g = h_d$ and $h_g = h_{BC}$, respectively.

The computation of sorptivity by applying Eq. (2) or (3) to hydraulic models, and in particular to those selected for this study, i.e., Eqs. (6)–(10), is quite tricky, given the complexity of the hydraulic functions. Such computation might exhibit the following shortcomings. First of all, the diffusivity functions must be determined analytically, by multiplying the hydraulic conductivity by the derivative of the water pressure head with regards to water content, which may involve complex algebraic operations. Then, the integration involved in the right-hand side of Eq. (3) may lead to numerical determination for very low initial water pressure heads, in the case of very dry initial conditions. Meanwhile, the integration involved in the right-hand side of Eq. (2) may pose numerical shortcomings for infinite hydraulic diffusivity, which is the case for some of the hydraulic functions detailed above, Eqs. (6)–(10).

In this study, we propose a specific scaling procedure to avoid all these shortcomings and to simplify the computation of sorptivity for the hydraulic models described in Eqs. (6)–(10) under the boundary conditions of a slightly positive water pressure head at surface and relatively dry initial conditions. We focus on these conditions since they constitute the most common experimental conditions for most water infiltration experiments and related procedures for characterizing soil hydraulic properties (Angulo-Jaramillo et al., 2016). In particular, these conditions feature the Beerkan method
that involves pouring water into a ring placed on the ground (Braud, 2005; Lassabatere et al., 2006).

The rest of the paper is organized as follows. The theory section details the scaling procedure that relates the square sorptivity to the (i) square scaled sorptivity, $S_K^*(−∞, 0)$, which we show is equal to the parameter $c_p$, and (ii) the product of scale parameters and correcting factors accounting for the contribution of initial water contents. The square scaled sorptivity corresponds to the sorptivity of a unit soil (unit value for all the scale parameters, except the residual water content fixed at zero) and for the whole range of water pressure heads, i.e., $(-∞, 0]$. It depends only on the soil hydraulic shape parameters, and its determination features the main algebraic complexity of the whole scaling procedure; the rest relies on simple algebraic operations (multiplication and sums). It is then computed for the hydraulic models defined by Eqs. (6)–(10), either analytically when feasible or numerically, otherwise. For each model, it is computed and tabulated as a function of a shape index that characterizes the spreads of the water retention function (from gradual to stepwise shapes, corresponding to soils with a broad or a very narrow pore size distribution, respectively). The evolution of the square scaled sorptivity versus the shape index is compared and discussed between models. In the last section, we illustrate the application of the proposed scaling procedure. We show how the tabulated values of the square scaled sorptivity $c_p$ can be used to upscale sorptivity and easily provide the sorptivity corresponding to a zero water pressure head at the surface for relatively small initial water contents.

# 2 Theory

## 2.1 Global scaling procedure

The proposed scaling procedure relies on two main steps already used in some previous studies: (i) relating the actual square sorptivity $S_K^2(h_0, h_1)$ to the maximum square sorptivity $S_K^2(-∞, 0)$ by isolating the effect of initial and final conditions $(h_0$ and $h_1)$ and (ii) scaling hydraulic functions and sorptivity to split the contributions of shape and scale parameters and facilitate the computation of $S_K^2(h_0, h_1)$. In our case, we consider that the final conditions always involve a null water pressure head, $h_1 = 0$, and saturated conditions, $θ_1 = θ_s$, while initial conditions $(h_0, θ_0)$ may vary.

### 2.1.1 Isolating the effect of initial conditions

The first step of the proposed procedure involves the work of Haverkamp et al. (2005), who isolated the contributions of the initial and final conditions to sorptivity. These authors considered that, at dry initial states ($θ_0 ≤ 1/4 θ_s$) and a zero water pressure head at surface (i.e., $h_1 = 0$ and $θ_1 = θ_s$), the following approximation applied:

$$ \left\{ \begin{array}{l} S^2(\theta_0, \theta_s) \approx S^2(0, \theta_s) \frac{K_s - K_0}{K_s} \frac{\theta_s - \theta_0}{\theta_s} , \\ S^2(0, \theta_s) = c_p |h_g| K_s \theta_s \end{array} \right. \qquad (11) $$

where $K_0$ corresponds to the initial hydraulic conductivity $K_0 = K(θ_0)$. The parameter $c_p = \frac{S^2(0, θ_s)}{|h_g| K_s θ_s}$ is a proportionality constant that depends only on the hydraulic shape parameters (Haverkamp et al., 2005). By combining both expressions in Eq. (11), the sorptivity, $S^2(θ_0, θ_s)$, can be defined as the product of three different terms that account for the respective contributions of the shape parameters (lumped into the proportionality constant $c_p$), the scale parameters, $|h_g|$, $K_s$, and $θ_s$, and initial conditions, $θ_0$:

$$ S^2(\theta_0, \theta_s) \approx c_p |h_g| K_s \theta_s \frac{K_s - K_0}{K_s} \frac{\theta_s - \theta_0}{\theta_s}. \qquad (12) $$

Equation (12) offers an accurate and practical approximation for the computation of sorptivity in the case of Beerkan runs, i.e., for a zero water pressure head imposed at the soil surface, $h_1 = 0$. Equation (12) was frequently used for the treatment of Beerkan data and in particular in all BEST methods (Lassabatere et al., 2006; Yilmaz et al., 2010; Angulo-Jaramillo et al., 2019). However, it addresses the case of soils with null residual water content, $θ_r = 0$, and without any air-entry pressure head, $h_a = 0$. In addition, it was developed and used for the vGB model.

In this study, we adapt Eq. (12) to any type of hydraulic model, including those with non-null residual water contents, $θ_r > 0$, and air-entry water pressure heads, $h_a < 0$. First of all, we consider that $θ_r$ must be accounted for, and we replace Eq. (11) by the following equation:

$$ \frac{S^2(\theta_0, \theta_s)}{S^2(\theta_r, \theta_s)} \approx \frac{K_s - K_0}{K_s} \frac{\theta_s - \theta_0}{\theta_s - \theta_r}. \qquad (13) $$

Indeed, the denominator $θ_s$ must be replaced by $(θ_s - θ_r)$ when $θ_r ≠ 0$ to ensure that the ratio $\frac{S^2(θ_0, θ_s)}{S^2(θ_r, θ_s)}$ tends towards unity when $θ_0 → θ_r$. Secondly, we consider that the approximation behind Eq. (11) involves only the unsaturated part of sorptivity, i.e., $S_D^2(θ_0, θ_s)$. As mentioned above, when the air-entry water pressure head is non-null, the computation of sorptivity $S_K^2(h_0, 0)$ must be split into its unsaturated and saturated parts, as illustrated by Eq. (5).

The following derivations are then proposed:

$$ \begin{aligned} S_K^2(h_0, h_1 = 0) &= S_D^2(θ_0, θ_s) + 2(θ_s - θ_0) K_s (h_1 - h_a) \\ &= \frac{K_s - K_0}{K_s} \frac{θ_s - θ_0}{θ_s - θ_r} S_D^2(θ_r, θ_s) \\ &\quad + 2(θ_s - θ_0) K_s |h_a| \quad \text{since } h_1 = 0 \\ &= \frac{K_s - K_0}{K_s} \frac{θ_s - θ_0}{θ_s - θ_r} S_D^2(θ_r, θ_s) \\ &\quad + \frac{θ_s - θ_0}{θ_s - θ_r} (2(θ_s - θ_r) K_s |h_a|), \end{aligned} \qquad (14) $$
where $S_D^2(\theta_r, \theta_s)$ and $2(\theta_s - \theta_r)K_s|h_a|$ refer to the unsaturated and saturated parts of the square sorptivity $S_K^2(-\infty, 0)$ (see Eq. 5 with $h_0 = -\infty$, i.e., $\theta_0 = \theta_r$, and $h_1 = 0$). Equation (14) can be simplified by introducing correcting factors $R_\theta$ and $R_K$ to account for the contribution of the initial conditions:

$$S_K^2(h_0, 0) = R_K R_\theta S_D^2(\theta_r, \theta_s) + R_\theta (2(\theta_s - \theta_r) K_s|h_a|), \quad (15)$$

with $R_\theta$ and $R_K$ defined as follows:

$$\left\{ \begin{array}{l} R_\theta = \frac{\theta_s - \theta_r}{\theta_s - \theta_r} \\ R_K = \frac{K_s - K_0}{K_s}. \end{array} \right. \qquad (16)$$

### 2.1.2 Scaling sorptivity and defining parameters $c_p$ and $c'_p$

So far, the computation of the square sorptivity $S_K^2(h_0, 0)$ has been treated considering dimensional equations. The second step of the proposed procedure scales the sorptivity in order to alleviate any numerical difficulty in relation to the values of dimensional variables. The square dimensional sorptivity $S_K^2(h_0, 0)$ can be easily related to the square scaled sorptivity $S_K^{*2}(h_0^*, 0)$ by scaling variables (water content, water pressure head, and hydraulic conductivity) as follows (Ross et al., 1996):

$$\left\{ \begin{array}{l} S_e = \frac{\theta_s - \theta_r}{\theta_s - \theta_r} \\ h^* = \frac{h}{|h_g|} \\ K_r = \frac{K}{K_s}. \end{array} \right. \qquad (17)$$

This scaling procedure defines the dimensionless water retention, $S_e(h^*)$, the dimensionless (or relative) hydraulic conductivity, $K_r(S_e)$, and the dimensionless hydraulic diffusivity function $D^*(S_e) = K_r(S_e) \frac{dh^*}{dS_e}$. These dimensionless hydraulic functions define the hydraulic characteristics of a unit soil that has the same values for the shape parameters and the unit value for all the scale parameters, $\theta_s = 1$, $h_g = 1$, and $K_s = 1$, except the residual water content that is fixed at zero, $\theta_r = 0$. In other words, these hydraulic parameters define the so-called “unit soil”. The use of the scaled expressions in Eq. (17) allows us to relate the square sorptivity (dimensional soil), $S^2$, to the square scaled sorptivity (unit soil), $S^{*2}$ (Ross et al., 1996):

$$S^2 = S^{*2} |h_g| K_s (\theta_s - \theta_r). \qquad (18)$$

Then, $S^{*2}$ can be computed by applying Eqs. (2) or (3) to the dimensionless hydraulic functions as a function of the initial and final water pressure heads, $h_0^*$ and $h_1^*$, or saturation degrees, $S_{e,0} = S_e(h_0^*)$ and $S_{e,1} = S_e(h_1^*)$:

$$\left\{ \begin{array}{l} S_K^{*2}(h_0^*, h_1^*) = \int_{h_0^*}^{h_1^*} (S_{e,1} + S_e(h^*) - 2S_{e,0}) K_r(h^*) dh^* \\ S_D^{*2}(S_{e,0}, S_{e,1}) = \int_{S_{e,0}}^{S_{e,1}} (S_{e,1} + S_e - 2S_{e,0}) D^*(S_e) dS_e. \end{array} \right. \qquad (19)$$

$S_D^{*2}$ and $S_K^{*2}$ define the diffusivity and conductivity forms of the square scaled sorptivity and are related to each other by Eqs. (4) and (5), leading to

$$\left\{ \begin{array}{l} S_K^{*2}(h_0^*, h_1^* \le h_a^*) = S_D^{*2}(S_{e,0}, S_{e,1}) \\ S_K^{*2}(h_0^*, h_1^* > h_a^*) = S_D^{*2}(S_{e,0}, 1) + 2(1 - S_{e,0})(h_1^* - h_a^*) \end{array} \right. \qquad (20)$$

with $h_a^* = \frac{h_a}{|h_g|}$. The scaled sorptivity functions, $S_D^{*2}$ and $S_K^{*2}$, depend only on the shape parameters and the initial and final saturation degrees.

By analogy with the approach of Haverkamp et al. (2005) (see Eq. 11), we scaled the maximum dimensional square sorptivity $S_K^2(-\infty, 0)$ and its unsaturated part $S_D^2(\theta_r, \theta_s)$, which are involved in Eq. (15), to define the parameters $c_p$ and $c'_p$:

$$\left\{ \begin{array}{l} c_p = S_K^{*2}(-\infty, 0) = \frac{S_K^2(-\infty, 0)}{(\theta_s - \theta_r)K_s|h_g|} \\ c'_p = S_D^{*2}(0, 1) = \frac{S_D^2(\theta_r, \theta_s)}{(\theta_s - \theta_r)K_s|h_g|}. \end{array} \right. \qquad (21)$$

The application of Eq. (19) with the proper lower and upper initial and final conditions ($S_{e,0} = 0$, $h_0 = -\infty$, $S_{e,1} = 1$, $h_1 = 0$) leads to the following expression for the two parameters $c_p$ and $c'_p$:

$$\left\{ \begin{array}{l} c_p = \int_{-\infty}^{0} (1 + S_e(h^*)) K_r(h^*) dh^* \\ c'_p = \int_{0}^{1} (1 + S_e) D^*(S_e) dS_e. \end{array} \right. \qquad (22)$$

The parameters $c_p$ and $c'_p$ depend exclusively on the soil shape parameters. The application of Eq. (20) provides a very simple equation linking the two parameters:

$$c_p = c'_p + 2|h_a^*|. \qquad (23)$$

### 2.1.3 Final expansion for computing sorptivity

The combination of the two previous steps (Sect. 2.1.1 and 2.1.2) leads to the final expression for sorptivity. As mentioned above, in order to avoid numerical problems, we first compute the scaled sorptivity before upscaling. The application of the scaling procedure, Eq. (18), to Eq. (15) leads to the following equivalent equations for the square scaled sorptivity $S_K^{*2}(h_0, 0)$:

$$\begin{aligned} S_K^{*2}(h_0^*, 0) &= c'_p R_K R_\theta + 2|h_a^*| R_\theta \\ &= (c_p - 2|h_a^*|) R_K R_\theta + 2|h_a^*| R_\theta, \end{aligned} \qquad (24)$$

where the correcting factors defined in Eq. (16) can be scaled and expressed as a function of the initial saturation degree and relative hydraulic conductivity:

$$\left\{ \begin{array}{l} R_\theta = \frac{\theta_s - \theta_r}{\theta_s - \theta_r} = 1 - S_{e,0} \\ R_K = \frac{K_s - K_0}{K_s} = 1 - K_r(S_{e,0}). \end{array} \right. \qquad (25)$$
These developments, based on the combination of the equation proposed by Haverkamp et al. (2005) to isolate the effect of initial and final conditions and the scaling procedure proposed by Ross et al. (1996), provide a new simple equation for a straightforward computation of sorptivity $S_K^2(h_0, 0)$:

$$ egin{cases} S_K^2(h_0, 0) = S_K^{*2}(h_0^*, 0) (\theta_s - \theta_r) K_s |h_g| \\ S_K^{*2}(h_0^*, 0) = (c_p - 2|h_a^*|) R_K R_\theta + 2|h_a^*| R_\theta \\ R_\theta = 1 - S_{e,0} \\ R_K = 1 - K_r(S_{e,0}). \end{cases} \quad (26) $$

The application of Eq. (26) requires the prior determination of the parameter $c_p$. In the following, we compute the value of the parameter $c_p$ for different hydraulic models.

## 2.2 Scaling hydraulic functions

### 2.2.1 General expressions

The first step of the determination of $c_p$ requires the computation of the dimensionless functions, $S_e(h^*)$, $K_r(S_e)$, and $D^*(S_e)$, to be inserted into Eq. (22). The application of the scaling variables, Eq. (17), to the hydraulic functions defined by Eqs. (6)–(10) leads to the following expressions for the dimensionless models:

delta model:

$$ egin{cases} S_{e,d}(h^*) = H(1+h^*) \\ K_{r,d}(S_e) = H(S_e - 1) \end{cases} \quad (27) $$

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 (28) $$

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 (29) $$

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 (30) $$

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 (31) $$

Note that the scaling parameter for the water pressure head, $h_g$, used in Eqs. (6)–(10) was set equal to the air-entry pressure head for the delta and BC WR functions; i.e., $h_d$ and $h_{BC}$ were set equal to $h_a$, as mentioned above.

Equations (27)–(31) were then used to derive the following formulations for the dimensionless diffusivity functions, $D^*(S_e)$, applying $D^*(S_e) = K_r(S_e) \frac{dh^*}{dS_e}$ (see Appendix A):

$$ D_d^*(S_e) = \delta(S_e) \quad (32) $$

$$ D_{BC}^*(S_e) = \frac{1}{\lambda_{BC}} S_e^{\eta_{BC} - \left(\frac{1}{\lambda_{BC}} + 1\right)} \quad (33) $$

$$ D_{vGB}^*(S_e) = \frac{1 - m_{vGB}}{2m_{vGB}} S_e^{\eta_{vGB} - \frac{1+m_{vGB}}{2m_{vGB}}} \left(1 - S_e^{\frac{1}{m_{vGB}}}\right)^{-\frac{1+m_{vGB}}{2}} \quad (34) $$

$$ D_{vGM}^*(S_e) = \frac{1 - m_{vGM}}{m_{vGM}} S_e^{l_{vGM} - \frac{1}{m_{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) \quad (35) $$

$$ 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)}, \quad (36) $$

where $\delta$ stands for the one-sided Dirac delta function (Triadis and Broadbridge, 2012).

### 2.2.2 Further simplifications

In the following section, several simplifications are proposed based on previous studies and the literature (Angulo-Jaramillo et al., 2016; Haverkamp et al., 2005). Several authors used capillary models to relate the HC to the WR functions. For the vGB model, Haverkamp et al. (2005) linked the shape parameter related to the HC function, $\eta$, with the combination of those of the WR function, $\lambda = mn$, as follows:

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

where the tortuosity parameter, $p$, takes the value of 1 for the case of Burdine's condition. We also consider the same equation for the BC model given its similarity with the vGB model (as demonstrated below, in the Results section). In addition, further simplifications involved the values of the tortuosity parameters, $l_{vGM}$ and $l_{KG}$, in the vGM and KG models. These parameters were fixed at the default values: $l_{vGM} = l_{KG} = 1/2$ (Şimünek et al., 2003; Kosugi, 1996; Kosugi and Hopmans, 1998). In practice, these shape parameters rarely vary (Haverkamp et al., 2005). With these supplementary considerations, the diffusivity functions for the BC, vGB, vGM, and KG models become

$$ D_{BC}^*(S_e) = \frac{1}{\lambda_{BC}} S_e^{\frac{1}{\lambda_{BC}} + 2} \quad (38) $$
$$D_{v\text{GB}}^*(S_e) = \frac{1 - m_{v\text{GB}}}{2m_{v\text{GB}}} S_e^{\frac{1+3m_{v\text{GB}}}{2m_{v\text{GB}}}} \left(1 - S_e^{\frac{1}{m_{v\text{GB}}}}\right)^{-\frac{1+m_{v\text{GB}}}{2}} \quad (39)$$

$$D_{v\text{GM}}^*(S_e) = \frac{1 - m_{v\text{GM}}}{m_{v\text{GM}}} S_e^{\frac{m_{v\text{GM}}-2}{2m_{v\text{GM}}}} \left( \left(1 - S_e^{\frac{1}{m_{v\text{GM}}}}\right)^{-m_{v\text{GM}}} + \left(1 - S_e^{\frac{1}{m_{v\text{GM}}}}\right)^{m_{v\text{GM}}} - 2 \right) \quad (40)$$

$$D_{\text{KG}}^*(S_e) = \frac{1}{2} \sqrt{\frac{\pi}{2}} \sigma_{\text{KG}} S_e^{\frac{1}{2}} \left( \text{erfc} \left( \text{erfc}^{-1}(2S_e) + \frac{\sigma_{\text{KG}}}{\sqrt{2}} \right) \right)^2 e^{\left( \text{erfc}^{-1}(2S_e) \right)^2 + \sqrt{2}\sigma_{\text{KG}} \text{erfc}^{-1}(2S_e)}. \quad (41)$$

The set of equations (Eqs. 38–41) shows that each hydraulic diffusivity function involves only one shape parameter, i.e., $\lambda_{BC}$, $m_{vGB}$, $m_{vGM}$, and $\sigma_{KG}$ for the respective BC, vGB, vGM, and KG hydraulic models, respectively. In the following, we consider both the general (i.e., Eqs. 33–36) and simplified (i.e., Eqs. 38–41) versions of the hydraulic diffusivity functions for the analytical determination of $c_p$, then we perform numerical applications only for the simplified versions in Eqs. (38)–(41).

## 2.3 Integral determination of parameter $c_p$

Once the dimensionless diffusivity functions are determined, the use of Eq. (22) allows for the determination of $c_p$, either analytically or numerically, depending on the considered hydraulic models.

### 2.3.1 $c_p$ for the delta model

This case is the easiest one. Indeed, the HC function is characterized by a null hydraulic conductivity for $h^* < -1$ and a unit value $K_r = 1$ for $h^* \ge -1$, as featured by Eq. (27). In this case, Eq. (22) yields

$$\begin{aligned} c_{p,d} &= \int_{-\infty}^{0} (1 + S_{e,d}(h^*)) K_{r,d}(h^*) dh^* \\ &= \int_{-\infty}^{-1} (1 + 0) \cdot 0 \cdot dh^* + \int_{-1}^{0} (1 + 1) \cdot 1 \cdot dh^* \\ &= 2. \end{aligned} \quad (42)$$

Note that this value of 2 was already proposed by Haverkamp et al. (2005) for the Green and Ampt soils, as defined by these authors, which corresponds to the delta model for the description of WR and HC functions.

### 2.3.2 $c_p$ for the Brooks and Corey (BC) model

The BC model involves an air-entry water pressure head $h_a^* = -1$. We used Eq. (22) while accounting for the saturated part of sorptivity and used the diffusivity form for the

determination of the unsaturated sorptivity. Indeed, the diffusivity function $D_{BC}^*(S_e)$ shown in Eq. (33) obeys a power law, which makes it possible to integrate the diffusivity form of sorptivity analytically. Simple algebraic operations and integrations of Eq. (22) lead to the following equation (see demonstration in Appendix B1):

$$c_{p,BC}(\lambda_{BC}, \eta_{BC}) = 2 + \frac{1}{\lambda_{BC}\eta_{BC} - 1} + \frac{1}{\lambda_{BC}\eta_{BC} + \lambda_{BC} - 1}. \quad (43)$$

When the relation between $\eta_{BC}$ and $\lambda_{BC}$ is ruled by Eq. (37) with $p = 1$, the analytical expression of $c_p$ turns into

$$c_{p,BC}(\lambda_{BC}) = 2 + \frac{1}{3\lambda_{BC} + 1} + \frac{1}{4\lambda_{BC} + 1}. \quad (44)$$

Our results are in line with previous studies (Varado et al., 2006).

### 2.3.3 $c_p$ for the van Genuchten–Burdine (vGB) model

In the case of the vGB model, there is no air-entry pressure head, i.e., $h_a^* = 0$. Equation (22) shows that $c_p$ reverts to the diffusivity form of the square scaled sorptivity, $\int_0^1 (1 + S_e) D^*(S_e) dS_e$. In addition, the diffusivity function, $D_{v\text{GB}}^*(S_e)$ (Eq. 34), makes the square scaled sorptivity analytically integrable, leading to (see demonstration in Appendix B2)

$$c_{p,v\text{GB}}(m_{v\text{GB}}, n_{v\text{GB}}, \eta_{v\text{GB}}) = \Gamma\left(1 + \frac{1}{n_{v\text{GB}}}\right) \left[ \frac{\Gamma\left(m_{v\text{GB}}\eta_{v\text{GB}} - \frac{1}{n_{v\text{GB}}}\right)}{\Gamma(m_{v\text{GB}}\eta_{v\text{GB}})} + \frac{\Gamma\left(m_{v\text{GB}}\eta_{v\text{GB}} + m_{v\text{GB}} - \frac{1}{n_{v\text{GB}}}\right)}{\Gamma(m_{v\text{GB}}\eta_{v\text{GB}} + m_{v\text{GB}})} \right], \quad (45)$$

where $\Gamma$ is the gamma function:

$$\Gamma(z) = \int_{0}^{+\infty} t^{z-1} e^{-t} dt. \quad (46)$$

Considering the relations between $m$ and $n$, i.e., $m = 1 - \frac{2}{n}$, and the relation between $\eta$ and $\lambda = mn$ in Eq. (37) with $p = 1$, the following simplification emerges:

$$c_{p,v\text{GB}}(m_{v\text{GB}}) = \Gamma\left(\frac{3 - m_{v\text{GB}}}{2}\right) \left[ \frac{\Gamma\left(\frac{1+5m_{v\text{GB}}}{2}\right)}{\Gamma(1 + 2m_{v\text{GB}})} + \frac{\Gamma\left(\frac{1+7m_{v\text{GB}}}{2}\right)}{\Gamma(1 + 3m_{v\text{GB}})} \right]. \quad (47)$$
The expression corresponding to Eq. (45) was already proposed and discussed by Haverkamp et al. (2005).

### 2.3.4 c<sub>p</sub> for the van Genuchten–Mualem (vGM) model

In contrast with the vGB model, no analytical expressions have been reported so far in the literature for this model. By analogy with the case of vGB model, analytical developments were proposed to analytically integrate the diffusivity form of the scaled sorptivity, leading to the following analytical expression for c<sub>p</sub> (see demonstration in Appendix B3):

$$c_{p,vGM}(m_{vGM}, l_{vGM}) = \Gamma(2 - m_{vGM}) \left( \frac{\Gamma(m_{vGM}(1+l_{vGM}))}{(m_{vGM}(1+l_{vGM}) - 1)\Gamma(m_{vGM}l_{vGM})} + \frac{\Gamma(m_{vGM}(2+l_{vGM}))}{(m_{vGM}(2+l_{vGM}) - 1)\Gamma(m_{vGM}(1+l_{vGM}))} + (1 - m_{vGM}) \left[ \frac{\Gamma(m_{vGM}(1+l_{vGM}))\Gamma(1+m_{vGM})}{(m_{vGM}(1+l_{vGM}) - 1)\Gamma(m_{vGM}(2+l_{vGM}))} + \frac{\Gamma(m_{vGM}(2+l_{vGM}))\Gamma(1+m_{vGM})}{(m_{vGM}(2+l_{vGM}) - 1)\Gamma(m_{vGM}(3+l_{vGM}))} \right] - 2 \left( \frac{1}{m_{vGM}(1+l_{vGM}) - 1} + \frac{1}{m_{vGM}(2+l_{vGM}) - 1} \right) \right]. \quad (48)$$

Note that this equation requires that $m_{vGM} \neq \frac{1}{1+l_{vGM}}$ and $m_{vGM} \neq \frac{1}{2+l_{vGM}}$. Considering that shape parameter $l_{vGM} = \frac{1}{2}$, as usually considered, Eq. (48) can be simplified to

$$c_{p,vGM}(m_{vGM}) = \Gamma(2 - m_{vGM}) \left( \frac{\Gamma\left(\frac{3m_{vGM}}{2}\right)}{\left(\frac{3m_{vGM}}{2} - 1\right)\Gamma\left(\frac{m_{vGM}}{2}\right)} + \frac{\Gamma\left(\frac{5m_{vGM}}{2}\right)}{\left(\frac{5m_{vGM}}{2} - 1\right)\Gamma\left(\frac{3m_{vGM}}{2}\right)} + (1 - m_{vGM}) \left[ \frac{\Gamma\left(\frac{3m_{vGM}}{2}\right)\Gamma(1 + m_{vGM})}{\left(\frac{3m_{vGM}}{2} - 1\right)\Gamma\left(\frac{5m_{vGM}}{2}\right)} + \frac{\Gamma\left(\frac{5m_{vGM}}{2}\right)\Gamma(1 + m_{vGM})}{\left(\frac{5m_{vGM}}{2} - 1\right)\Gamma\left(\frac{7m_{vGM}}{2}\right)} \right] - 2 \left( \frac{1}{\frac{3m_{vGM}}{2} - 1} + \frac{1}{\frac{5m_{vGM}}{2} - 1} \right) \right]. \quad (49)$$

These sets of equations have never been proposed and constitute one of the novel outputs of this study. The complexity of algebraic developments for the derivation of Eqs. (48) and (49) makes them valuable.

### 2.3.5 c<sub>p</sub> for the Kosugi (KG) model

No analytical formulation was found for the case of Kosugi's hydraulic functions. Therefore, the square scaled sorptivity was computed numerically with a generic procedure that can be applied to any type of hydraulic model (i.e., any set of HC and WR functions). To avoid integration over infinite intervals with respect to $h^*$ and integration of an infinite diffusivity close to saturation, the integral was split into two parts, leading to the following developments:

$$\begin{align*} c_{p,KG}(\sigma_{KG}, l_{KG}) &= \int_{-\infty}^{0} (1 + S_{e,KG}(h^*)) K_{r,KG}(h^*) dh^* \\ &= \int_{-\infty}^{h_{KG}^*(\frac{1}{2})} (1 + S_{e,KG}(h^*)) K_{r,KG}(h^*) dh^* \\ &\quad + \int_{h_{KG}^*(\frac{1}{2})}^{0} (1 + S_{e,KG}(h^*)) K_{r,KG}(h^*) dh^* \\ &= \int_{0}^{\frac{1}{2}} (1 + S_{e}) D_{KG}^*(S_{e}) dS_{e} \\ &\quad + \int_{h_{KG}^*(\frac{1}{2})}^{0} (1 + S_{e,KG}(h^*)) K_{r,KG}(h^*) dh^*, \quad (50) \end{align*}$$

where $h_{KG}^*(\frac{1}{2})$ is the water pressure head corresponding to $S_e = \frac{1}{2}$. In the last expression of Eq. (50), c<sub>p</sub> is composed of two integrals of continuous functions over closed intervals, which are thus well defined and easily computable. Again, a simplified version is proposed, assuming that the shape parameter $l_{KG}$ is fixed to $\frac{1}{2}$.

## 2.4 Shape indexes for comparing c<sub>p</sub> between the selected hydraulic models

The approach described below allows c<sub>p</sub> to be determined for the selected BC, vGB, vGM, and KG models. In addition to its contribution to simplifying Eq. (26), c<sub>p</sub> has a real physical meaning: it corresponds to the square sorptivity of unit soils with a zero water pressure head at the soil surface and utterly dry initial profile. Therefore, c<sub>p</sub> should not depend on the choice of the hydraulic models. We then investigate the dependency of c<sub>p</sub> upon the selected hydraulic model. Consequently, we designed shape indexes to compare c<sub>p</sub> between models. These shape indexes were built to describe the same state of the WR functions regardless of the chosen hydraulic model. We designed these shape indexes to vary WR functions between two extreme states: (i) values close to zero for gradual WR functions, corresponding to soils with broad pore size distributions, and (ii) values of unity for stepwise
WR functions mimicking an abrupt change of water content, corresponding to soils with very narrow pore size distributions. This section presents the design of these shape indexes.

A sensitivity analysis of vGB and vGM models showed that the parameter *m* is adequate for varying the WR functions between a gradual shape (*m* = 0) and a stepwise function (*m* = 1). We thus consider *m*<sub>vGB</sub> and *m*<sub>vGM</sub> to be the appropriate shape indexes, i.e., *x*<sub>vGB</sub> = *m*<sub>vGB</sub> and *x*<sub>vGM</sub> = *m*<sub>vGM</sub>. Next, we designed the shape index for BC, *x*<sub>BC</sub>, deriving from vGB, *m*<sub>vGB</sub>, considering that vGB and BC functions describe close WR curves. Indeed, for large values of *n*<sub>vGB</sub>, the vGB model converges towards a power function similar to the BC model (Haverkamp et al., 2005):

$$
\begin{align}
\lim_{n_{v\text{GB}} \to +\infty} (1 + |h^*|^{n_{v\text{GB}}})^{-m_{v\text{GB}}} &\approx |h^*|^{-n_{v\text{GB}}} m_{v\text{GB}} \nonumber \\
&\approx |h^*|^{-\lambda_{v\text{GB}}} \tag{51}
\end{align}
$$

The equation $\lambda_{BC} = \lambda_{vGB} = n_{vGB} - 2$ defines a relation between $\lambda_{BC}$ and $n_{vGB}$ that ensures a similar state for WR functions. Substituting $n_{vGB}$ with $m_{vGB}$ according to $m_{vGB} = 1 - \frac{2}{n_{vGB}}$ leads to

$$m_{vGB} = \frac{\lambda_{BC}}{2 + \lambda_{BC}}. \qquad (52)$$

Since $m_{vGB}$ is the appropriate shape index for the vGB model, we consider its equivalent, $\frac{\lambda_{BC}}{2+\lambda_{BC}}$, to be the appropriate shape index for the BC model, leading to

$$x_{BC} = \frac{\lambda_{BC}}{2 + \lambda_{BC}}. \qquad (53)$$

For KG functions, we considered that stepwise WR functions are associated with a narrow pore size distribution, i.e., a null standard deviation, $\sigma_{KG}$. In contrast, gradual WR functions correspond to a spread distribution of pore size distribution, i.e. very large values of $\sigma_{KG}$. Consequently, by analogy with Eq. (53), i.e. using a ratio, we propose the following shape index, $\sigma_{KG}$:

$$x_{KG} = \frac{1}{1 + \sigma_{KG}}. \qquad (54)$$

Finally, the hydraulic shape parameters for each model can be expressed as a function of the shape index, by inverting the previous equations. For the sake of simplicity, we use the same letter, *x*, to denote the different shape indexes, $x_{BC}$, $x_{vGB}$, $x_{vGM}$, or $x_{KG}$. We obtain the following relations:

$$
\left\{
\begin{array}{l}
\lambda_{BC} = \frac{2x}{1-x} \\
m_{vGB} = x \\
m_{vGM} = x \\
\sigma_{KG} = \frac{1-x}{x},
\end{array}
\right.
\qquad (55)
$$

where *x* takes values in the interval [0, 1]. Equation (55) provides the shape parameters of the studied models for a given

value of shape index *x*, i.e., for a similar state of WR function between gradual (*x* = 0) and stepwise functions (*x* = 1).

On the basis on these relations between hydraulic shape parameters and indexes, $c_p$ can easily be related to the shape index using Eqs. (44), (47), and (49), obtained for the simplified diffusivity functions, Eqs. (38)–(41), based on the capillary model, Eq. (37). For the KG model, the computations of $c_p$ remained numerical. The following set of equations were obtained:

$$
\left\{
\begin{array}{l}
c_{p,BC}(x) = 2 + \frac{1-x}{5x+1} + \frac{1-x}{7x+1} \\
c_{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_{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)} \right. \\
\qquad + \left(\frac{\Gamma\left(\frac{5}{2}x\right)}{\left(\frac{5}{2}x-1\right)\Gamma\left(\frac{3}{2}x\right)}\right) + (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)} \right. \\
\qquad \left. + \left(\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) - 2\left(\frac{1}{\frac{3}{2}x-1} + \frac{1}{\frac{5}{2}x-1}\right) \right] \\
c_{p,KG}(x) = \int_{0}^{\frac{1}{2}} (1 + S_e) D_{KG}(x) (S_e) \, dS_e \\
+ \int_{0}^{0} (1 + S_e, KG(x) (h^*)) K_{r,KG}(x) (h^*) \, dh^* \\
\text{with } \sigma_{KG}(x) = \frac{1-x}{x}.
\end{array}
\right.
\qquad (56)
$$

We performed a sensitivity analysis by varying the shape index *x* for each model between 0 and 1 by increments of 0.025. We then computed the different shape parameters for BC, vGB, vGM, and KG models using Eq. (55) and then plotted the related hydraulic functions Eqs. (28)–(31) with $\eta = \frac{2}{\lambda} + 2 + p$ and $l_{vGM} = l_{KG} = \frac{1}{2}$ and diffusivity functions Eqs. (38)–(41). Lastly, we computed the square scaled sorptivity $c_p$ as a function of the shape index, Eq. (56), and discussed the function $c_p(x)$ regarding the choice of the hydraulic model. The values of $c_p(x)$ are also compared to those of the delta model (Dirac delta functions), i.e. $c_{p,d} = 2$.

# 3 Results

## 3.1 Analysis of the hydraulic functions and hydraulic diffusivity functions

The hydraulic and diffusivity functions are plotted in Fig. 1 for the shape index value of 0.275, and their sensitivity upon each model shape index is shown in Fig. 2. For the sake of clarity, we plotted the relative hydraulic conductivity both as functions of saturation degree, $K_r(S_e)$, and water pressure head, $K_r(h^*)$, noting that these functions have distinct uses: $K_r(S_e)$ defines the HC functions as a property of the soils, whereas $K_r(h^*)$ is mostly used to compute sorptivity, e.g., as in Eq. (19). Thus, $K_r(h^*)$ has a similar role as the diffusivity
function $D^*(S_e)$, and the shapes and properties of these two
functions determine the values of the square scaled sorptiv-
ity $c_p$.

The comparison between the hydraulic models (with the
same shape index value) reveals some similarities and dis-
crepancies (Fig. 1). Three of the water retention models
(vGB, vGM, KG) exhibit an inflection point with a con-
tinuous increase in $S_e(h^*)$ over the whole interval ($-\infty, 0$]
(Fig. 1a, vGB, vGM, and KG), while the BC model reaches
the asymptote $S_e = 1$ at $h^* = -1$ with full saturation for $h^* \ge$
$-1$ (Fig. 1a, “BC”). Despite that difference, the BC and vGB
models exhibit similar shapes (Fig. 1a, BC versus vGB). The
vGM model exhibits a more progressive increase in $S_e(h^*)$
while remaining asymmetrically distributed across the inflec-
tion point (Fig. 1a, vGM). Lastly, the KG model exhibits
an even more progressive increase and a perfect symmetry
around the inflection point (Fig. 1a, KG). The position of the
inflection points depends on the chosen hydraulic model. By
construction, the inflection point is positioned at $h^* = -1$ for
the KG model. The other models have inflection points po-
sitioned at larger abscissas (in absolute values), with similar
intermediate values for the BC and the vGB models and the
largest abscissas for the vGM model (Fig. 1a).

Regarding the relative hydraulic conductivity, the BC and
vGB models have similar shapes for $K_r(S_e)$, both typical
of power functions (Fig. 1b, BC and vGB). In contrast, the
vGM and KG models have an inflection point, with a larger
increase both at low saturation degrees and close to satura-
tion compared to intermediate saturation degrees. In partic-
ular, these two models exhibit a very large increase close to
saturation, whereas BC and vGB models have a gradual in-
crease (Fig. 1b, vGM and KG versus BC and vGB, close to
$S_e = 1$). This feature allows the vGM and KG models to sim-
ulate large drops in hydraulic conductivity close to satura-
tion that are typical of certain soils. The functions $K_r(h^*)$ are
depicted in Fig. 1c and combine the properties of the func-
tions $K_r(S_e)$ and $S_e(h^*)$, as described above. The function
$K_r(h^*)$ exhibits similar shapes for BC and vGB models, with
a quasi-linear and sharp increase for $h^* \le -1$ followed by
a plateau (Fig. 1c, BC and vGB). The vGM and KG mod-
els exhibit a much more progressive increase, involving a
much larger range of water pressure heads. This feature is
the most pronounced for the KG model. This more progres-
sive increase reflects the more gradual WR functions, $S_e(h^*)$,
as described in Fig. 1a combined with the drop in $K_r(S_e)$ that
reaches unity only for saturation degrees extremely close to
unity (Fig. 1b).

As for the WR and HC functions, the diffusivity functions
exhibit close shapes for the BC and vGB models (Fig. 1d). The
BC model defines a concave shape with a finite maxi-
mum equal to $\lambda_{BC}$ obtained at $S_e = 1$, in line with the use
of Eq. (38) at $S_e = 1$ (Fig. 1d, BC). Conversely, the vGB
model diverts from the concave shape to tend towards infin-
ity close to saturation, at $S_e = 1$ (Fig. 1d, vGB versus BC).
The vGB model defines an S shape that reflects the larger in-

creases at both low and high saturation degrees, with a lower
increase at intermediate saturation degrees (Fig. 1d, vGB).
The two other models, vGM and KG, exhibit the same type
of S shape, with an infinite limit close to saturation (Fig. 1d,
vGM and KG). Such an infinite limit spoils the numerical
integration of Eq. (22) for the determination of $c_p$, requiring
the use of the mixed formulation for the KG models defined
by Eqs. (50) and (56).

Varying the shape index changes the WR and HC func-
tions in an expected way (Fig. 2). For the WR functions, in-
creasing the shape index from 0 to 1 makes the shift from a
gradual and moderate to an abrupt increase in saturation de-
gree, respectively. Values close to unity make the WR func-
tions close to a stepwise function correspond to the delta
model (Fig. 2, first column, arrows). As for the WR func-
tions, the increase in the shape index puts the curves $K_r(h^*)$
close to stepwise functions (Fig. 2, third column, arrows).
For the BC model, we notice a decrease in $K_r(h^*)$ for $h^* \le
-1$, whereas $K_r(h^*)$ remains equal to unity above (Fig. 2c,
third column). Conversely, for the vGB, vGM, and KG mod-
els, the increase in the shape index has two antagonist effects:
a decrease of $K_r(h^*)$ for $h^* \le -1$, followed by an increase
for $h^* \ge -1$ (Fig. 2g, k, o, third column, arrows). Briefly, as
expected, the water retention and the relative hydraulic con-
ductivity tend towards stepwise functions when the shape in-
dex tends towards unity (Fig. 2, first and third columns).
This trend is less evident for the diffusivity functions (Fig. 2,
4th column).

These results show that, regardless of the model selected,
increasing the shape index put the hydraulic functions closer
to the delta model that corresponds to soils with narrow pore
size distribution. Conversely, very small values of the shape
index ensure very gradual shapes for WR and HC functions.
However, the results point at contrasting trends when the
shape index is decreased towards zero. It is clear that for
the vGB and BC models, the relative hydraulic conductiv-
ity $K_r(h^*)$ is not greatly impacted close to $h^* = 0$ (Fig. 2c
and g). Conversely, for the vGM and KG models, $K_r(h^*)$
tends towards zero in the vicinity of $h^* = 0$ (Fig. 2k and o,
inverted arrows). Similarly, the dimensionless diffusivity,
$D^*(S_e)$, tends towards zero over the whole interval [0,1]
when the shape index tends towards zero (Fig. 2l and p,
inverted arrows). Consequently, the features of $K_r(h^*)$ and
$D^*(S_e)$ functions suggest that the square scaled sorptivity is
very small for vGM and KG models, when the shape index
tends towards zero (see results below).

## 3.2 Square scaled sorptivity $c_p$ as a function of shape indexes

The square scaled sorptivity, $c_p(x)$, is plotted as a function
of the shape index, $x$ (Fig. 3). The results show contrasting
evolution. For the BC model, we note a decrease of $c_{p,BC}(x)$
from 4 down to 2 (Fig. 3a, $c_p$). Such a decrease is expected
since $K_{r,BC}(h^*)$ functions decrease over the whole inter-
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Examples of water retention, $S_e(h^*)$ (a), relative unsaturated hydraulic, $K_r(S_e)$ (b) and $K_r(h^*)$ (c), and diffusivity $D^*(S_e)$ functions (d) for the four hydraulic models, Brooks and Corey (BC), van Genuchten–Burdine (vGB), van Genuchten–Mualem (vGM), and Kosugi (KG); the curves were plotted for the same value of the shape index, $x = 0.275$. The hydraulic parameters $\lambda_{BC}$, $m_{vGM}$, $m_{vGB}$, and $\sigma_{KG}$ were computed as a function of $x$ using Eq. (55) with $l_{vGM} = l_{KG} = \frac{1}{2}$. The dashed line represents the delta model.</figcaption></figure>

val $(-\infty, 0]$ with the shape index (see Fig. 2c). This leads to the decrease of the integral involved in Eq. (22) and thus $c_p$. The upper and lower limits can be easily determined by applying Eq. (56), leading to $c_{p,BC}(0) = 4$ and $c_{p,BC}(1) = 2$. The air-entry water pressure head is non-null, so that the $c_p$ and $c'_p$ are related by $c_p = c'_p + 2|h_a^*| = c'_p + 2$, since by convention $h_a^* = -1$. Consequently, $c'_p$ decreases from 2 to 0, with a simple vertical shift compared to $c_p$ (Fig. 3a, $c'_p$). For the other hydraulic models, the air-entry water pressure head is null, so that $c_p = c'_p$ (Fig. 3b–d). In the following, we compare the evolution of $c_p$ as a function of the shape index for each model.

In contrary to the BC model, $c_{p,vGB}(x)$ for the vGB model does not decrease monotonically (Fig. 3c). Instead, $c_{p,vGB}(x)$ decreases to 1.5 before increasing up to 2 for $x \ge 0.52$. This feature is in line with the effect of the shape index on the relative hydraulic conductivity, $K_r(h^*)$, as described above. The shape index has two antagonist effects: a decrease of $K_r(h^*)$ for $h^* \le -1$ and an increase for $h^* \ge -1$ (Fig. 2g, arrows). The numerical computation sorts out these contrasting effects and demonstrates the two-step variation, i.e., a decrease followed by an increase. The two boundaries of the function $c_{p,vGB}(x)$ can be easily found using Eq. (56) with a

lower limit of $c_{p,vGB}(0) = 2\Gamma(\frac{3}{2})\Gamma(\frac{1}{2}) = \Gamma(\frac{1}{2})^2 = \pi$ and an upper limit of $c_{p,vGB}(1) = \Gamma(1)[\frac{\Gamma(3)}{\Gamma(3)} + \frac{\Gamma(4)}{\Gamma(4)}] = 2$.

For the KG and vGM models, the trend is opposite, and the functions $c_{p,KG}(x)$ and $c_{p,vGM}(x)$ both increase (Fig 3b and d). This increase is in line with the fact that the shape index mainly increases the diffusivity function, $D^*(S_e)$, over the whole interval $[0, 1]$ and thus the integral of Eq. (22). The two functions $c_{p,KG}(x)$ and $c_{p,vGM}(x)$ increase from 0 up to 2. For the vGM model, the lower and upper limits can be demonstrated using Eq. (56), leading to $c_{p,vGM}(0) = 0$ (numerical determination) and $c_{p,vGM}(1) = \Gamma(1)[\frac{\Gamma(\frac{1}{2})}{\Gamma(\frac{1}{2})} + \frac{\Gamma(\frac{3}{2})}{\Gamma(\frac{3}{2})}] = 2$. For KG hydraulic functions, the lower and upper limits were determined numerically, leading also to 0 and 2, with almost zero values over a large interval of shape index, i.e., $x \in [0, 0.3]$ (Fig. 3b).

The four functions $c_{p,BC}(x)$, $c_{p,vGB}(x)$, $c_{p,vGM}(x)$, and $c_{p,KG}(x)$ all reach the value of 2 when the shape index approaches unity, i.e., when the WR and HC functions tend towards stepwise functions. In fact, the value of $c_p(x)$ converges to the value obtained for the delta model, $c_{p,d} = 2$ (see Eq. 42). Such a result indicates that similar results should
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Impact of the shape index, *x*, on the WR and HC functions versus the selected hydraulic models: WR functions, *S*<sub>e</sub>(*h*<sup>*</sup>) (a, e, i, m); HC functions, *K*<sub>r</sub>(*S*<sub>e</sub>) (b, f, j, n) or *K*<sub>r</sub>(*h*<sup>*</sup>) (c, g, k, o); and diffusivity function *D*<sup>*</sup>(*S*<sub>e</sub>) functions (d, h, l, p) for the 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 increasing values of the shape index *x*. The hydraulic parameters λ<sub>BC</sub>, *m*<sub>vGM</sub>, *m*<sub>vGB</sub>, and σ<sub>KG</sub> were computed as a function of *x* using Eq. (55), with *l*<sub>vGM</sub> = *l*<sub>KG</sub> = 1/2.</figcaption></figure>

be obtained for soils with a narrow pore size distribution
(like coarse soils with narrow pore size distributions), re-
gardless of the model selected for describing their WR and
HC functions. In other words, the choice of the hydraulic
model should not matter for soils with narrow pore size dis-
tributions. Conversely, contrasting trends are obtained when
the shape index tends towards zero: quasi-null values for
the vGM and KG models, π for the vGB model, and 4 for
the BC model. The null values of c<sub>p</sub> for the KG and vGB
models are in line with the large decrease in D*(S<sub>e</sub>) over
the whole interval [0, 1] when the shape index, x, tends to-
wards zero (see Fig. 2l and p and comments above). Such
results show that the choice of the hydraulic model matters
for soils with graded pore size distributions. Between these
two extreme states, the four functions c<sub>p,BC</sub>(x), c<sub>p,vGB</sub>(x),

c<sub>p,vGM</sub>(x), and c<sub>p,KG</sub>(x) exhibit contrasting evolution, with
c<sub>p,vGM</sub>(x) and c<sub>p,KG</sub>(x) increasing monotonously, c<sub>p,BC</sub>(x)
decreasing monotonously, and c<sub>p,vGB</sub>(x) exhibiting a two-
step behavior of decreasing followed by increasing values.

This contrasting evolution of functions c<sub>p,BC</sub>(x),
c<sub>p,vGB</sub>(x), c<sub>p,vGM</sub>(x), and c<sub>p,KG</sub>(x) points at different impli-
cations regarding the physics of flow and infiltration in soils.
It should be borne in mind that the choice of the hydraulic
models should not impact the value of the square sorptivity,
S<sup>2</sup><sub>K</sub>(h<sub>0</sub>, 0), for a given soil and given initial conditions, h<sub>0</sub>.
Indeed, sorptivity should be independent of the choice of
hydraulic model since it always equals the ratio between
the cumulative infiltration and the square root of time for
gravity-free infiltration, as illustrated by Eq. (1). However,
this work shows that the choice of the hydraulic model
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Square scaled sorptivity, $c_p$, as a function of shape index, $x$, for the four hydraulic models, (a) Brooks and Corey (BC), (b) Ko-sugi (KG), (c) van Genuchten–Burdine (vGB), and (d) van Genuchten–Mualem (vGM); the models were sorted as a function of their monotonic features, with monotonically increasing functions on the right.</figcaption></figure>

strongly impacts the estimation of $c_p$ for soils with broad pore size distributions. When the shape index gets close to zero, the BC and vGB models predict non-null values for $c_p$, thus ensuring non-null values of the sorptivity (see Eq. 26), whereas the KG or vGM models predict quasi-null values of $c_p$. We then expect the product of scale parameters, $(\theta_s - \theta_r)K_s|h_g|$, to compensate for the very low values of $c_p$ (see Eq. 26). Given that the scale parameters $\theta_s$, $\theta_r$, and $K_s$ characterize dry (residual) or saturated states of the soil, these parameters are not expected to vary between the hydraulic models and are supposed fixed. Consequently, only the scale parameter for the water pressure head, $|h_g|$, is expected to compensate for the very low values of $c_p$ when the vGM and KG models are used. We conclude that the value of $|h_g|$ must tends towards being infinite when the shape index tends towards zero for these two models. For the KG model, such a relation, $|h_{KG}| \sim x_{KG}$, is related to the relation $|h_{KG}| \sim \sigma_{KG}$ since $x_{KG} = (1 + \sigma_{KG})^{-1}$. Our statements imply that very large values of $\sigma_{KG}(x_{KG} \to 0)$ should be associated with very large values of $|h_{KG}|$, i.e., with a very small pore radius. In other words, fine (coarse) soils with small (large) pores should have broad (narrow) pore size distributions. These considerations are in line the previous studies on the KG model by Pollacco et al. (2013) (see Fig. 1 of their paper) and Fernández-Gálvez et al. (2019). These authors even related the pore radius of

soil to the standard deviation of pore size distribution with a strongly decreasing function. Similar trends relating to some extent a relation between pore mean and pore standard deviation should also apply to the vGM model given our findings (and to avoid null sorptivity for soils with broad pore size distributions). More investigations are needed to verify such a link between average pore size and standard deviation for real soils. More investigations are also needed to guide on the proper choice of the hydraulic models as a function of the soil type, as already suggested by Fuentes et al. (1992). These aspects will be the subject of a specific study.

Regarding the numerical accuracy of computed values of $c_p$, we used analytical formulations for the BC, vGB, and vGM models, as detailed in Eq. (56). These values are expected to be perfect without any error since they correspond to the application of exact analytical formulations. Instead, we used the mixed numerical formulation defined by Eq. (50) for the KG model that relies on the numerical integration of the HC or diffusivity functions. In that case, the numerical integration may bring some numerical errors. The mixed form, Eq. (50), was designed to minimize numerical indetermination and uncertainty. Such a formulation was applied to the other models (BC, vGB, and vGM), and the resulting values were compared against the analytical formulations (considered to be the benchmark). A perfect agreement was obtained
(errors < 1 %), thus validating the numerical mixed formulation and reinforcing confidence in the values in Table 1. Note that the promotion of the numerical mixed formulation, Eq. (50), and its uncertainty will be the subject of another study.

### 3.3 Upscaling sorptivity $S_K(h_0, 0)$ from $c_p$

In this section, we elaborate on the use of Eq. (26) for the easy and straightforward computation of $S_K^2(h_0, 0)$ from the tabulated values of $c_p$ (Table 1). The proposed scaling procedure Eq. (26) allows for the computation of $S_K^2(h_0, 0)$ given initial conditions (water contents or water pressure heads), hydraulic shape and scale parameters, and specific hydraulic models selected among Eqs. (7)–(10):

1. Use the shape parameter ($\lambda_{BC}$, $m_{vGB}$, $m_{vGM}$, or $\sigma_{KG}$) to compute the related shape index, $x$, considering the following definitions: $x_{vGB} = m_{vGB}$, $x_{vGM} = m_{vGM}$, $x_{KG} = \frac{1}{1+\sigma_{KG}}$, and $x_{BC} = \frac{\lambda_{BC}}{2+\lambda_{BC}}$.

2. Choose the value of $c_p$ in Table 1 corresponding to the shape index, $x$, and the chosen WR and HC functions.

3. Consider or compute the initial water content, $\theta_0$ or $\theta(h_0)$, depending on the description of the initial condition (either $\theta_0$ or $h_0$).

4. Compute the related hydraulic conductivity, $K_0 = K(\theta_0)$, using the HC function.

5. Compute the correcting factors $R_{\theta} = \frac{\theta_s - \theta_0}{\theta_s - \theta_r} (= 1 - S_{e,0})$ and $R_K = \frac{K_s - K_0}{K_s} (= 1 - K_r(S_{e,0}))$.

6. Compute the scaled air-entry water pressure head, $|h_a^*| = |\frac{h_a}{h_g}|$.

7. Compute the square scaled sorptivity, $S_K^{*2}(h_0^*, 0) = R_K R_{\theta}(c_p - 2|h_a^*|) + 2R_{\theta}|h_a^*|$.

8. Upscale to derive the square sorptivity, $S_K^2(h_0, 0) = S_K^{*2}(h_0^*, 0)(\theta_s - \theta_r)K_s|h_g|$.

As an illustrative example, let us consider the case of a loamy soil at water saturation with a slightly positive water pressure head at the surface ($h_1 = 0$) and an initial water pressure head of $h_0 = -10$ m (dry conditions). The loamy soil has the features of loam as defined in the database of Carsel and Parrish (1988). Its WR and HC functions are described by the vGM model, with the following shape and scale parameters: $\theta_r = 0.078$, $\theta_s = 0.43$, $h_g = -277$ mm, $K_s = 2.88 \times 10^{-3}$ mm s$^{-1}$, $n_{vGM} = 1.56$, and $l_{vGM} = 0.5$. The application of the step-by-step procedure gives the following results:

1. shape index: $x = m_{vGM} = 1 - 1/n_{vGM}$, leading to $x = 0.359$;

2. corresponding value of $c_p$ in Table 1 (vGM model column, $x = 0.36$): $c_p = 0.480$;

3. initial water content: computed from the initial water pressure head of $-10$ m using the vGM-WR function, i.e., Eq. (9): $\theta_0 = 0.125$;

4. initial hydraulic conductivity: computed from the initial water content using vGM-HC function, i.e., Eq. (9): $K_0 = 1.87 \times 10^{-9}$ mm s$^{-1}$;

5. corresponding correction factors: $R_{\theta} = \frac{\theta_s - \theta_0}{\theta_s - \theta_r}$ and $R_K = \frac{K_s - K_0}{K_0}$, leading to $R_{\theta} = 0.865$ and $R_K = 1.000$;

6. air-entry water pressure head: no air-entry water pressure head, leading to $|h_a^*| = 0$;

7. square scaled sorptivity: $S_K^{*2}(h_0^*, 0) = R_K R_{\theta} c_p$, leading to $S_K^{*2}(h_0^*, 0) = 0.416$;

8. sorptivity: $S_K^2(h_0, 0) = S_K^{*2}(h_0^*, 0)(\theta_s - \theta_r)K_s|h_g|$, leading to $S_K^2(h_0, 0) = 0.117 \text{ mm}^2 \text{ s}^{-1}$ and $S_K(h_0, 0) = 0.342 \text{ mm s}^{-\frac{1}{2}}$.

To check the accuracy of the proposed approximation, we computed the nominal value of sorptivity, using Eq. (2). We found a very close value, with less than 0.5 % relative error, demonstrating the accuracy of the proposed scaling procedure in Eq. (26).

As a second illustrative example, we consider the computation of sorptivity for the case of the BC model, for the same conditions. The difference with the previous case is that the BC model has a non-null air-entry water pressure head, inducing a non-null saturated sorptivity. We consider the same loamy soil with the following parameters for the BC model: $\theta_r = 0.078$, $\theta_s = 0.43$, $h_g = -277$ mm, and $K_s = 2.88 \times 10^{-3}$ mm s$^{-1}$, with a value of $\lambda_{BC} = 0.56$. $\lambda_{BC}$ was deduced from the previous value of $n = 1.56$ considering the usual relation $\lambda = mn$, as suggested by Haverkamp et al. (2005). The application of the proposed procedure leads to the following computations:

1. shape index: $x_{BC} = \lambda_{BC}/(2 + \lambda_{BC})$, leading to $x_{BC} = 0.219$;

2. corresponding value of $c_p$ in Table 1 (BC model column, $x_{BC} = 0.22$): $c_p = 2.678$;

3. initial water content: computed from the initial water pressure head of $-10$ m using BC-WR function, i.e., Eq. (7): $\theta_0 = 0.125$;

4. initial hydraulic conductivity: computed from the initial water content using BC-HC function, i.e., Eq. (7): $K_0 = 3.342 \times 10^{-9}$ mm s$^{-1}$;

5. corresponding correction factors: $R_{\theta} = \frac{\theta_s - \theta_0}{\theta_s - \theta_r}$ and $R_K = \frac{K_s - K_0}{K_0}$, leading to $R_{\theta} = 0.866$ and $R_K = 1.000$;
<table><caption>**Table 1.** Values of the square scaled sorptivity $c_p$ as a function of the shape index $x$ for the studied hydraulic models: Brooks and Corey (BC), van Genuchten–Burdine (vGB), van Genuchten–Mualem (vGM), and Kosugi (KG).</caption><thead><tr><th>x</th><th>BC</th><th>vGB</th><th>vGM</th><th>KG</th><th>x</th><th>BC</th><th>vGB</th><th>vGM</th><th>KG</th></tr></thead><tbody><tr><td>0.00</td><td>4.000</td><td>3.142 (π)</td><td>0.000</td><td>0.000</td><td>...</td><td>...</td><td>...</td><td>...</td><td>...</td></tr><tr><td>0.02</td><td>3.749</td><td>2.891</td><td>2.559 × 10<sup>-3</sup></td><td>3.38 × 10<sup>-776</sup></td><td>0.52</td><td>2.237</td><td>1.605</td><td>0.820</td><td>0.604</td></tr><tr><td>0.04</td><td>3.549</td><td>2.692</td><td>9.814 × 10<sup>-3</sup></td><td>2.14 × 10<sup>-188</sup></td><td>0.54</td><td>2.220</td><td>1.605</td><td>0.864</td><td>0.686</td></tr><tr><td>0.06</td><td>3.384</td><td>2.529</td><td>2.126 × 10<sup>-2</sup></td><td>3.106 × 10<sup>-81</sup></td><td>0.56</td><td>2.205</td><td>1.608</td><td>0.908</td><td>0.768</td></tr><tr><td>0.08</td><td>3.246</td><td>2.393</td><td>3.646 × 10<sup>-2</sup></td><td>2.125 × 10<sup>-44</sup></td><td>0.58</td><td>2.191</td><td>1.611</td><td>0.953</td><td>0.848</td></tr><tr><td>0.10</td><td>3.129</td><td>2.279</td><td>5.500 × 10<sup>-2</sup></td><td>1.161 × 10<sup>-27</sup></td><td>0.60</td><td>2.177</td><td>1.616</td><td>0.998</td><td>0.927</td></tr><tr><td>0.12</td><td>3.028</td><td>2.183</td><td>7.654 × 10<sup>-2</sup></td><td>9.772 × 10<sup>-19</sup></td><td>0.62</td><td>2.164</td><td>1.623</td><td>1.043</td><td>1.004</td></tr><tr><td>0.14</td><td>2.940</td><td>2.099</td><td>0.101</td><td>1.875 × 10<sup>-13</sup></td><td>0.64</td><td>2.151</td><td>1.631</td><td>1.089</td><td>1.079</td></tr><tr><td>0.16</td><td>2.863</td><td>2.028</td><td>0.127</td><td>4.350 × 10<sup>-10</sup></td><td>0.66</td><td>2.140</td><td>1.640</td><td>1.135</td><td>1.151</td></tr><tr><td>0.18</td><td>2.794</td><td>1.966</td><td>0.156</td><td>7.998 × 10<sup>-8</sup></td><td>0.68</td><td>2.128</td><td>1.651</td><td>1.181</td><td>1.220</td></tr><tr><td>0.20</td><td>2.733</td><td>1.911</td><td>0.187</td><td>3.107 × 10<sup>-6</sup></td><td>0.70</td><td>2.117</td><td>1.662</td><td>1.228</td><td>1.287</td></tr><tr><td>0.22</td><td>2.678</td><td>1.864</td><td>0.219</td><td>4.429 × 10<sup>-5</sup></td><td>0.72</td><td>2.107</td><td>1.676</td><td>1.275</td><td>1.351</td></tr><tr><td>0.24</td><td>2.629</td><td>1.823</td><td>0.253</td><td>3.218 × 10<sup>-4</sup></td><td>0.74</td><td>2.097</td><td>1.690</td><td>1.322</td><td>1.412</td></tr><tr><td>0.26</td><td>2.584</td><td>1.787</td><td>0.288</td><td>1.463 × 10<sup>-3</sup></td><td>0.76</td><td>2.088</td><td>1.706</td><td>1.370</td><td>1.471</td></tr><tr><td>0.28</td><td>2.543</td><td>1.755</td><td>0.325</td><td>4.760 × 10<sup>-3</sup></td><td>0.78</td><td>2.079</td><td>1.723</td><td>1.418</td><td>1.527</td></tr><tr><td>0.30</td><td>2.506</td><td>1.727</td><td>0.362</td><td>1.212 × 10<sup>-2</sup></td><td>0.80</td><td>2.070</td><td>1.741</td><td>1.467</td><td>1.580</td></tr><tr><td>0.32</td><td>2.471</td><td>1.704</td><td>0.401</td><td>2.568 × 10<sup>-2</sup></td><td>0.82</td><td>2.062</td><td>1.761</td><td>1.516</td><td>1.630</td></tr><tr><td>0.34</td><td>2.440</td><td>1.683</td><td>0.440</td><td>4.735 × 10<sup>-2</sup></td><td>0.84</td><td>2.054</td><td>1.782</td><td>1.566</td><td>1.680</td></tr><tr><td>0.36</td><td>2.410</td><td>1.665</td><td>0.480</td><td>7.837 × 10<sup>-2</sup></td><td>0.86</td><td>2.046</td><td>1.804</td><td>1.617</td><td>1.727</td></tr><tr><td>0.38</td><td>2.383</td><td>1.650</td><td>0.521</td><td>0.119</td><td>0.88</td><td>2.039</td><td>1.828</td><td>1.669</td><td>1.771</td></tr><tr><td>0.40</td><td>2.358</td><td>1.637</td><td>0.562</td><td>0.170</td><td>0.90</td><td>2.032</td><td>1.853</td><td>1.721</td><td>1.814</td></tr><tr><td>0.42</td><td>2.334</td><td>1.627</td><td>0.604</td><td>0.229</td><td>0.92</td><td>2.025</td><td>1.879</td><td>1.775</td><td>1.854</td></tr><tr><td>0.44</td><td>2.312</td><td>1.619</td><td>0.646</td><td>0.295</td><td>0.94</td><td>2.018</td><td>1.907</td><td>1.829</td><td>1.893</td></tr><tr><td>0.46</td><td>2.291</td><td>1.613</td><td>0.689</td><td>0.367</td><td>0.96</td><td>2.012</td><td>1.937</td><td>1.885</td><td>1.930</td></tr><tr><td>0.48</td><td>2.272</td><td>1.608</td><td>0.732</td><td>0.444</td><td>0.98</td><td>2.006</td><td>1.967</td><td>1.942</td><td>1.966</td></tr><tr><td>0.50</td><td>2.254</td><td>1.606</td><td>0.776</td><td>0.523</td><td>1.00</td><td>2.000</td><td>2.000</td><td>2.000</td><td>2.000</td></tr></tbody></table>

6. air-entry water pressure head: significant air-entry water pressure head, with $h_{BC} = h_a$, leading to $|h_a^*| = 1$;

7. square scaled sorptivity: $S_K^*(h_0^*, 0) = R_K R_\theta(c_p - 2|h_a^*|) + 2R_\theta|h_a^*|$, leading to $S_K^*(h_0^*, 0) = 2.318$;

8. sorptivity: $S_K^2(h_0, 0) = S_K^2(h_0^*, 0)(\theta_s - \theta_r)K_s|h_g|$, leading to $S_K^2(h_0, 0) = 0.651 \text{ mm}^2 \text{ s}^{-1}$, and $S_K(h_0, 0) = 0.806 \text{ mm s}^{-\frac{1}{2}}$.

Again, the exact value of sorptivity was estimated using Eq. (3) and led to a similar value with a relative error of 1‰. Note that, in this case, due to the non-null air-entry water pressure head, Eq. (3) must be employed instead of Eq. (2) for the determination of the targeted value of sorptivity.

The two preceding applications illustrated the accuracy of the proposed scaling procedure, Eq. (26), for both hydraulic functions, with and without air-entry water pressure heads. Equation (26) proved appropriate and very accurate for the determination of the sorptivity, $S_K^2(h_0, 0)$.

It must be noted that the proposed scaling procedure applies only for a dry initial state. Indeed, Haverkamp et al. (2005) stated that their approximation, Eq. (11), was valid only when $\theta_0 \le \frac{1}{4}\theta_s$. For fine soils, even a small initial wa-

ter pressure head may cause $\theta_0 > \frac{1}{4}\theta_s$, which may spoil the proposed scaling procedure. To illustrate this point, we investigated the case of silty clay soil, as defined by Carsel and Parrish (1988). This soil is defined with the following parameters: $\theta_r = 0.07$, $\theta_s = 0.36$, $h_g = -2000 \text{ mm}$, $K_s = 5.555 \times 10^{-5} \text{ mm s}^{-1}$, $n = 1.09$, and $l = 0.5$. Considering the same value for the initial water pressure head, i.e., $h_0 = -10 \text{ m}$, the initial water content is $\theta_0 = 0.318$, which exceeds $\frac{1}{4}\theta_s$. The application of the scaling procedure led to an estimated sorptivity of $0.0475 \text{ mm s}^{-\frac{1}{2}}$, whereas the targeted sorptivity computed with Eq. (2) was $0.0127 \text{ mm s}^{-\frac{1}{2}}$. Such an error corresponds to an overestimation by a factor of 2.73. Thus, we advise that the user verify that $\theta_0 \le \frac{1}{4}\theta_s$ when $\theta_r = 0$, or more generally $S_{e,0} \le \frac{1}{4}$ before using the proposed scaling procedure.

## 4 Conclusions

The proper estimation of sorptivity is crucial to understand and model water infiltration into soils. However, its estimation may be complicated, requiring complicated algebraic derivations and exhibiting potential numerical shortcomings
when using Eq. (2) or Eq. (3). In this study, we present a new scaling procedure for simplifying the computation of sorptivity for the case of a zero water pressure head imposed at the surface and dry initial state. We based our approach on the combination and adaptation of the scaling procedure proposed by Ross et al. (1996) and the approximation proposed by Haverkamp et al. (2005). We then obtain a simple relation that relates the square sorptivity to the product of the square scaled sorptivity, referred to as $c_p$, the product of scale parameters, and two correction factors that account for the initial conditions (i.e., initial water content and hydraulic conductivity). The value of the square scaled sorptivity, $c_p$, was computed either analytically, when feasible, or numerically, for four hydraulic models: Brooks and Corey, van Genuchten–Mualem, van-Genuchten–Burdine, and Kosugi models. The values of $c_p$ were tabulated as a function of specific shape indexes, representing similar states of WR functions (well-graded versus stepwise shapes) between hydraulic models. The proposed scaling procedure is easy of use, and, given a selected hydraulic model with related hydraulic shape and scale parameters, the following steps are conducted: computation of the shape index from the hydraulic shape parameters, reading of the corresponding value of $c_p$ in Table 1, computation of the correction factors (ratios in hydraulic conductivity and water contents, $R_K$ and $R_\theta$), computation of the square scaled sorptivity from $c_p$ and the correction factors, and, lastly, upscaling by multiplying with the scale parameters. All these steps are easy to conduct and straightforward. Illustrative examples are proposed at the end of this study, and the accuracy of the proposed scaling procedure is clearly demonstrated (with errors less than 1%), provided that the initial water content fulfills the conditions $R_\theta \le \frac{1}{4}R_K$ when $\theta_r = 0$, or more generally $S_{e,0} \le \frac{1}{4}$.

In addition to providing a straightforward method for the determination of sorptivity, this study brings very interesting findings on the square scaled sorptivity, $c_p$, and its dependency upon the shape index, $x$, and the chosen hydraulic model. The results show that the function $c_p(x)$ strongly depends on the hydraulic model selected for the WR and HC curves. If all the functions $c_p(x)$ converge for the same value, i.e., 2, close to $x = 1$ (stepwise WR functions – narrow pore size distribution), they strongly divert close to $x = 0$ (graded WR functions –broad pore size distribution), with values of 0 for vGM and KG models versus 3–4 for the vGB and BC models. However, the sorptivity should remain the same regardless of the selected hydraulic model: one soil under particular initial conditions means one single sorptivity. Consequently, the contrast of scaled sorptivity must be compensated by a contrast in scale parameters. However, among scale parameters, the residual and saturated water contents and the saturated hydraulic conductivity cannot be changed between models, since they characterize the dry and saturated states of the same soil. Consequently, the value of the scale parameter, $h_g$, must be the one to compensate. Previous studies on the Kosugi model have already

hypothesized a strong relation between the scale parameter, $h_{KG}$, and the standard deviation, $\sigma_{KG}$ (Pollacco et al., 2013). In other words, the scale parameter, $h_{KG}$, should be parametrized as a function of the shape parameter, $\sigma_{KG}$, to get plausible WR and HC functions and estimates of sorptivity. We may also expect the same link between the scale parameter, $h_{g,vGM}$, and the shape parameter, $m_{vGM}$, to avoid unphysical scenarios and null sorptivity. However, such a hypothesis has never been suggested and requires further investigations. These results show the need to better understand the mathematical properties of the hydraulic models, including the links between the hydraulic shape and scale parameters, and to better relate these properties to the physical processes of water infiltration into soils (Fuentes et al., 1992).

In addition to the proposed scaling procedure, this study gave the opportunity to derive the scaled sorptivity analytically for the three models, BC, vGB, and vGM, thus confirming the expressions provided by previous studies. For the vGM model, the analytical derivation is brand new and had never been proposed before. Its use is of great interest and could be implemented into soil hydraulic characterization methods. For instance, additional BEST methods could be developed on the basis of the use of the proposed formulations for the parameter $c_p$ to relate sorptivity to shape and scale parameters. In current BEST methods, only the vGB model is considered. The prior estimation of shape parameters allows for the determination of the parameter $c_p$ using Eq. (45). Then, the estimation of saturated hydraulic conductivity and sorptivity allows for the determination of scale parameter $h_g$ once $c_p$ is determined (Lassabatere et al., 2006). A similar procedure may be proposed for the vGM model, using Eq. (49), that relates the parameter $c_p$ to the shape parameter $m_{vGM}$. The development of the BEST method for the specific vGM hydraulic model, that is much more used than the vGB model, will be the subject of further investigations. Another improvement concerns the consideration of the residual water content $\theta_r$ in the proposed scaling procedure. It would also be interesting to somehow derive the residual water content and not to assume it to be equal to zero as it might also alter the shape of the soil water retention function. The use of the scaled sorptivity and the proposed scaling procedure for these purposes is the subject of ongoing studies and will lead to implementations in many methods and tools for the characterization of single- or dual-permeability soils, including the BEST methods (Fernández-Gálvez et al., 2019; Lassabatere et al., 2019) and the SoilWater-ToolBox software developed for the characterization of homogeneous and heterogeneous soils (Fernández-Gálvez et al., 2021).
## Appendix A: Dimensionless hydraulic diffusivity functions, $D^*(S_e)$

In the appendices, for the sake of clarity, the notation of the shape parameters was simplified to $\lambda$, $m$, $n$, $\sigma$, in order to avoid heavy equations. The dimensionless diffusivity functions were derived from their definition $D^*(S_e)$, applying $D^*(S_e) = K_r(S_e) \frac{dh^*}{dS_e}$. This task requires us to first derive the inverse functions for the dimensionless water retention curves. The following equations can be easily found through usual algebraic developments:

$$h_{BC}^*(S_e) = -S_e^{-\frac{1}{\lambda}} \qquad (A1)$$

$$h_{vGB}^*(S_e) = -\left(S_e^{-\frac{1}{m}} - 1\right)^{\frac{1}{n}} \quad \text{with } m = 1 - \frac{2}{n} \qquad (A2)$$

$$h_{vGM}^*(S_e) = -\left(S_e^{-\frac{1}{m}} - 1\right)^{\frac{1}{n}} \quad \text{with } m = 1 - \frac{1}{n} \qquad (A3)$$

$$h_{KG}^*(S_e) = -e^{\sqrt{2}\sigma\text{erfc}^{-1}(2S_e)}, \qquad (A4)$$

where $\text{erfc}^{-1}$ is the inverse function of the complementary error function. These functions can be differentiated to define their relative derivatives, $\frac{dh^*}{dS_e}$:

$$\frac{dh_{BC}^*}{dS_e}(S_e) = \frac{1}{\lambda} S_e^{-\frac{1}{\lambda}-1} \qquad (A5)$$

$$\frac{dh_{vGB}^*}{dS_e}(S_e) = \frac{1-m}{2m} S_e^{-\frac{1+m}{2m}} \left(1 - S_e^{\frac{1}{m}}\right)^{-\frac{1+m}{2}} \qquad (A6)$$

$$\frac{dh_{vGM}^*}{dS_e}(S_e) = \frac{1-m}{m} S_e^{-\frac{1}{m}} \left(1 - S_e^{\frac{1}{m}}\right)^{-m} \qquad (A7)$$

$$\frac{dh_{KG}^*}{dS_e}(S_e) = \sqrt{2\pi}\sigma e^{(\text{erfc}^{-1}(2S_e))^2 + \sqrt{2\sigma}\text{erfc}^{-1}(2S_e)}. \qquad (A8)$$

The differentiation of the function $(f \circ g)' = f' \circ g \cdot g'$ and $(f^{-1})' = \frac{1}{f' \circ f^{-1}}$, considering bijective functions. We also use the usual derivative of the “erf” function, $\text{erf}(x) = \frac{2}{\sqrt{\pi}}e^{-x^2}$, and the relation between erfc and erf functions, $\text{erfc}(x) = 1 - \text{erf}(x)$.

The scaled hydraulic conductivity functions are

$$K_{r,BC}(S_e) = S_e^{\eta}, \qquad (A9)$$

$$K_{r,vGB}(S_e) = S_e^{\eta}, \qquad (A10)$$

$$K_{r,vGM}(S_e) = S_e^l \left( 1 - 2 \left( 1 - S_e^{\frac{1}{m}} \right)^m + \left( 1 - S_e^{\frac{1}{m}} \right)^{2m} \right), \qquad (A11)$$

$$K_{r,KG}(S_e) = S_e^l \left( \frac{1}{2} \text{erfc} \left( \text{erfc}^{-1} (2S_e) + \frac{\sigma}{\sqrt{2}} \right) \right)^2. \qquad (A12)$$

Note that for the hydraulic conductivity function of the vGM model, $K_{r,vGM}$, we distributed the terms according to $(a+b)^2 = a^2+2ab+b^2$. The multiplication of Eqs. (A5)–(A8) by the expressions of relative conductivity, Eqs. (A9)–(A12), leads to the following expressions of dimensionless

diffusivity:

$$\left\{
\begin{aligned}
D_{BC}^*(S_e) &= \frac{1}{\lambda} S_e^{\eta - \left(\frac{1}{\lambda} + 1\right)} \\
D_{vGB}^*(S_e) &= \frac{1-m}{2m} S_e^{\eta - \frac{1+m}{2m}} \left(1 - S_e^{\frac{1}{m}}\right)^{-\frac{1+m}{2}} \\
D_{vGM}^*(S_e) &= \frac{1-m}{m} S_e^{l-\frac{1}{m}} \left(\left(1 - S_e^{\frac{1}{m}}\right)^{-m} + \left(1 - S_e^{\frac{1}{m}}\right)^m - 2\right) \\
D_{KG}^*(S_e) &= \frac{1}{2} \sqrt{\frac{\pi}{2}} \sigma S_e^l \left(\text{erfc}\left(\text{erfc}^{-1}(2S_e) + \frac{\sigma}{\sqrt{2}}\right)\right)^2 \\
e^{(\text{erfc}^{-1}(2S_e))^2 + \sqrt{2\sigma}\text{erfc}^{-1}(2S_e)}
\end{aligned}
\right.
\qquad (A13)$$

The expressions in Eq. (A13) correspond to the expressions of Eqs. (33)–(36). Afterwards, the combination with the capillarity model, Eq. (37), leads to Eqs. (38)–(41).

## Appendix B: Analytical developments for the $c_p$ parameter

### B1 Parameter $c_p$ for the BC model

For the BC model, we need to account for the air-entry pressure, $h_a^* = -1$. We remind the reader that by convention, the scale parameter for the water pressure head, $h_{BC}$, is equal to the air-entry water pressure head, $h_a$. We then use Eqs. (22) and (23) with $|h_a^*| = 1$. Then, the first part $\int_0^1 (1 + S_e) D_{BC}^*(S_e) dS_e$ can be integrated analytically given that the hydraulic diffusivity $D_{BC}^*(S_e)$ obeys a power law. The following developments can be done:

$$\begin{aligned}
c_{p,BC} &= \int_{0}^{1} (1 + S_e) D_{BC}^{*}(S_e) dS_e + 2|h_a^{*}| \\
&= \int_{0}^{1} (1 + S_e) D_{BC}^{*}(S_e) dS_e + 2 \qquad (B1)
\end{aligned}$$
$$
\begin{align*}
\int_0^1 (1 + S_e) D_{BC}^*(S_e) dS_e &= \int_0^1 D_{BC}^*(S_e) dS_e \\
&\quad + \int_0^1 S_e D_{BC}^*(S_e) dS_e \\
&= \int_0^1 \frac{1}{\lambda} S_e^{\eta - (\frac{1}{\lambda} + 1)} dS_e + \int_0^1 \frac{1}{\lambda} S_e S_e^{\eta - (\frac{1}{\lambda} + 1)} dS_e \\
&= \frac{1}{\lambda} \left( \int_0^1 S_e^{\eta - (\frac{1}{\lambda} + 1)} dS_e + \int_0^1 S_e^{\eta - \frac{1}{\lambda}} dS_e \right) \\
&= \frac{1}{\lambda} \left( \frac{1}{\eta - \frac{1}{\lambda}} + \frac{1}{\eta - \frac{1}{\lambda} + 1} \right) \\
&= \frac{1}{\eta\lambda - 1} + \frac{1}{\eta\lambda + \lambda - 1}. \tag{B2}
\end{align*}
$$

The final expression can be easily computed by combining
Eqs. (B1) and (B2):

$$
c_{p,BC}(\lambda, \eta) = \frac{1}{\eta\lambda - 1} + \frac{1}{\eta\lambda + \lambda - 1} + 2. \quad (\text{B3})
$$

The concatenation of Eq. (B3) with the capillary model,
Eq.( 37), leads to the following final expression, assuming
$p = 1$:

$$
c_{p,BC}(\lambda) = 2 + \frac{1}{3\lambda + 1} + \frac{1}{4\lambda + 1}. \quad (B4)
$$

These development demonstrate the equations proposed
for c<sub>p</sub> for the BC model, i.e., Eqs. (43)–(44). This demon-
stration is in line with Varado et al. (2006).

B2 Parameter c<sub>p</sub> for the vGB model

For the vGB model, and the remaining models, there is no
air-entry water pressure head, $h_a^* = 0$, leading to

$$
\begin{align}
c_{p,v\text{GB}} &= \int_0^1 (1 + S_e) D_{v\text{GB}}^*(S_e) dS_e + 2|h_a^*| \nonumber \\
&= \int_0^1 (1 + S_e) D_{v\text{GB}}^*(S_e) dS_e. \tag{B5}
\end{align}
$$

Then, the integral $\int_0^1 (1+S_e) D_{v\text{GB}}^*(S_e) dS_e$ can be decomposed
into well-known integrals:

$$
\begin{align*}
c_{p,v\text{GB}}(m, \eta) &= \int_0^1 (1 + S_e) D_{v\text{GB}}^*(S_e) dS_e \\
&= \int_0^1 D_{v\text{GB}}^*(S_e) dS_e + \int_0^1 S_e D_{v\text{GB}}^*(S_e) dS_e \\
&= \int_0^1 \frac{1-m}{2m} S_e^{\eta - \frac{1+m}{2m}} \left(1 - S_e^{\frac{1}{m}}\right)^{-\frac{1+m}{2}} dS_e \\
&\quad + \int_0^1 S_e \frac{1-m}{2m} S_e^{\eta - \frac{1+m}{2m}} \left(1 - S_e^{\frac{1}{m}}\right)^{-\frac{1+m}{2}} dS_e \\
&= \frac{1-m}{2m} \left( \int_0^1 S_e^{\eta - \frac{1+m}{2m}} \left(1 - S_e^{\frac{1}{m}}\right)^{-\frac{1+m}{2}} dS_e \right. \\
&\qquad \left. + \int_0^1 S_e S_e^{\eta - \frac{1+m}{2m}} \left(1 - S_e^{\frac{1}{m}}\right)^{-\frac{1+m}{2}} dS_e \right) \\
&= \frac{1-m}{2m} \left( \int_0^1 S_e^{\eta - \frac{1+m}{2m}} \left(1 - S_e^{\frac{1}{m}}\right)^{-\frac{1+m}{2}} dS_e \right. \\
&\qquad \left. + \int_0^1 S_e^{\eta+1-\frac{1+m}{2m}} \left(1 - S_e^{\frac{1}{m}}\right)^{-\frac{1+m}{2}} dS_e \right) \\
&= \frac{1-m}{2m} \left( \int_0^1 \left(S_e^{\frac{1}{m}}\right)^{m\eta-\frac{1+m}{2}} \left(1 - S_e^{\frac{1}{m}}\right)^{-\frac{1+m}{2}} dS_e \right. \\
&\qquad \left. + \int_0^1 \left(S_e^{\frac{1}{m}}\right)^{m\eta+m-\frac{1+m}{2}} \left(1 - S_e^{\frac{1}{m}}\right)^{-\frac{1+m}{2}} dS_e \right). \tag{B6}
\end{align*}
$$

The change of variable $y = S_e^{\frac{1}{m}}$ provides the following expressions:

$$
\begin{align}
c_{p,v\text{GB}}(m, \eta) &= \frac{1-m}{2} \int_0^1 y^{m\eta + \frac{m}{2} - \frac{3}{2}} (1-y)^{-\frac{1+m}{2}} \, dy \nonumber \\
&\quad + \frac{1-m}{2} \int_0^1 y^{m\eta + \frac{3m}{2} - \frac{3}{2}} (1-y)^{-\frac{1+m}{2}} \, dy. \tag{B7}
\end{align}
$$
We now use the beta function B, which has the following properties (assuming $x > 0$ and $y > 0$):

$$ 
\begin{aligned}
B(x, y) &= \int_{0}^{1} t^{x-1}(1-t)^{y-1}dt \\
&= \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)},
\end{aligned} 
\quad (B8) $$

where the $\Gamma$ function is already defined by Eq. (46): $\Gamma(z) = \int_{0}^{+\infty} t^{z-1}e^{-t}dt$ ($z > 0$). Then, the parameter $c_{p,vGB}$ can be expressed as follows:

$$ 
\begin{aligned}
c_{p,vGB}(m, \eta) &= \frac{1-m}{2} B\left(m\eta + \frac{m}{2} - \frac{1}{2}, \frac{1-m}{2}\right) \\
&\quad + \frac{1-m}{2} B\left(m\eta + \frac{3m}{2} - \frac{1}{2}, \frac{1-m}{2}\right) \\
&= \frac{1-m}{2} \Gamma\left(\frac{1-m}{2}\right) \left( \frac{\Gamma\left(m\eta + \frac{m}{2} - \frac{1}{2}\right)}{\Gamma(m\eta)} \right. \\
&\qquad \left. + \frac{\Gamma\left(m\eta + \frac{3m}{2} - \frac{1}{2}\right)}{\Gamma(m\eta + m)} \right) = \Gamma\left(\frac{3-m}{2}\right) \\
&\qquad \left( \frac{\Gamma\left(m\eta + \frac{m}{2} - \frac{1}{2}\right)}{\Gamma(m\eta)} + \frac{\Gamma\left(m\eta + \frac{3m}{2} - \frac{1}{2}\right)}{\Gamma(m\eta + m)} \right). 
\end{aligned} 
\quad (B9) $$

The last equation uses the fact that $z\Gamma(z) = \Gamma(z+1)$. Equation (B9) corresponds to the equation suggested by Haverkamp et al. (2005), considering $m = 1 - \frac{2}{n}$:

$$ 
\begin{aligned}
c_{p,vGB}(n, m, \eta) &= \Gamma\left(1 + \frac{1}{n}\right) \left( \frac{\Gamma\left(m\eta - \frac{1}{n}\right)}{\Gamma(m\eta)} \right. \\
&\qquad \left. + \frac{\Gamma\left(m\eta + m - \frac{1}{n}\right)}{\Gamma(m\eta + m)} \right). 
\end{aligned} 
\quad (B10) $$

Note that the proposed simplification using the beta function, Eq. (B9), requires that $m\eta + \frac{m}{2} - \frac{1}{2} = \frac{1}{2} + (\frac{3}{2} + p)m > 0$, which is quite evident since $m > 0$. When Eq. (B9) is combined with the capillary model, i.e., Eq. (37), the expression becomes

$$ c_{p,vGB}(m) = \Gamma\left(\frac{3-m}{2}\right) \left( \frac{\Gamma\left(\frac{1+5m}{2}\right)}{\Gamma(1+2m)} + \frac{\Gamma\left(\frac{1+7m}{2}\right)}{\Gamma(1+3m)} \right). \quad (B11) $$

The expressions of the parameter $c_p$ for the vGB model are accurately demonstrated, leading to Eqs. (45) and (47).

## B3 Parameter $c_p$ for the vGM model

For the vGM model, the same equation, Eq. (B5), applies and can be cut into two parts:

$$ 
\begin{aligned}
c_{p,vGM}(m, l) &= \int_{0}^{1} (1 + S_e) D_{vGM}^{*}(S_e) dS_e \\
&= \underbrace{\int_{0}^{1} D_{vGM}^{*}(S_e) dS_e}_{A} + \underbrace{\int_{0}^{1} S_e D_{vGM}^{*}(S_e) dS_e}_{B}.
\end{aligned} 
\quad (B12) $$

For the sake of clarity, we demonstrate the simplifications of the two terms A and B separately:

$$ 
\begin{aligned}
A &= \frac{1-m}{m} \int_{0}^{1} S_{e}^{l-\frac{1}{m}} \left( \left(1-S_{e}^{\frac{1}{m}}\right)^{-m} + \left(1-S_{e}^{\frac{1}{m}}\right)^{m} - 2 \right) dS_{e} \\
&= \frac{1-m}{m} \left( \int_{0}^{1} S_{e}^{l-\frac{1}{m}} \left(1-S_{e}^{\frac{1}{m}}\right)^{-m} dS_{e} \right. \\
&\qquad \left. + \int_{0}^{1} S_{e}^{l-\frac{1}{m}} \left(1-S_{e}^{\frac{1}{m}}\right)^{m} dS_{e} - 2 \int_{0}^{1} S_{e}^{l-\frac{1}{m}} dS_{e} \right) \\
&= \frac{1-m}{m} \left( \int_{0}^{1} \left(S_{e}^{\frac{1}{m}}\right)^{ml-1} \left(1-S_{e}^{\frac{1}{m}}\right)^{-m} dS_{e} \right. \\
&\qquad \left. + \int_{0}^{1} \left(S_{e}^{\frac{1}{m}}\right)^{ml-1} \left(1-S_{e}^{\frac{1}{m}}\right)^{m} dS_{e} - 2 \int_{0}^{1} S_{e}^{l-\frac{1}{m}} dS_{e} \right). 
\end{aligned} 
\quad (B13) $$

The last term of A can be simplified easily:

$$ \int_{0}^{1} S_{e}^{l-\frac{1}{m}} dS_{e} = \left[ \frac{S_{e}^{l-\frac{1}{m}+1}}{l-\frac{1}{m}+1} \right]_{0}^{1} = \frac{m}{(l+1)m-1}. \quad (B14) $$

For the two first terms, we use the same change of variable $y = S_e^{\frac{1}{m}}$ as above to transform the integrals, leading to

$$ A = (1-m) \left( \int_{0}^{1} y^{ml+m-2} (1-y)^{-m} dy + \int_{0}^{1} y^{ml+m-2} (1-y)^{m} dy - \frac{2}{(l+1)m-1} \right). \quad (B15) $$

In this case, we need to assume that $ml + m - 2 > -1$, i.e. $m > \frac{1}{l+1}$, to use the beta function. Such a condition corresponds to $m > \frac{2}{3}$, for the by-default value of $l = \frac{1}{2}$. The two first integrals can then be replaced using the beta and the gamma functions, leading to the final expression for part A
of c<sub>p</sub>:

$$
\begin{equation}
\begin{aligned}
A = (1-m)(B(m(l+1)-1, 1-m) + B(m(l+1)-1, 1+m)) & \\
\qquad - \frac{2}{(l+1)m-1} = (1-m)\left(\frac{\Gamma(m(l+1)-1)\Gamma(1-m)}{\Gamma(ml)}\right) & \\
\qquad + \frac{\Gamma(m(l+1)-1)\Gamma(1+m)}{\Gamma(m(l+2))} - \frac{2}{(l+1)m-1}.
\end{aligned}
\tag{B16}
\end{equation}
$$

By analogy, the following developments come out for the
term B:

$$
\begin{equation}
\begin{aligned}
B = (1-m) \left( \frac{\Gamma(m(l+2)-1)\Gamma(1-m)}{\Gamma(ml+m)} \right. & \\
& \left. + \frac{\Gamma(m(l+2)-1)\Gamma(1+m)}{\Gamma(m(l+3))} - \frac{2}{(l+2)m-1} \right).
\end{aligned}
\tag{B17}
\end{equation}
$$

The simplification for B is valid as soon as $m > \frac{1}{2+l}$, which
is the case since we suppose that $m > \frac{1}{1+l}$. After rearrang-
ing terms, the following expressions come out for the scaled
sorptivity:

$$
\begin{equation}
\begin{aligned}
c_{p,v\text{GM}}(m,l) ={}& \Gamma(2-m) \left( \frac{\Gamma(m(l+1)-1)}{\Gamma(ml)} \right. \\
& \left. + \frac{\Gamma(m(l+2)-1)}{\Gamma(ml+m)} \right) \\
& + (1-m) \left[ \left( \frac{\Gamma(m(l+1)-1)\Gamma(1+m)}{\Gamma(m(l+2))} \right. \right. \\
& \left. \left. + \frac{\Gamma(m(l+2)-1)\Gamma(1+m)}{\Gamma(m(l+3))} \right) \right. \\
& \left. - 2 \left( \frac{1}{(l+1)m-1} + \frac{1}{(l+2)m-1} \right) \right]
\end{aligned}
\tag{B18}
\end{equation}
$$

$$
\begin{equation}
\begin{aligned}
c_{p,v\text{GM}}(m,l) ={}& \Gamma(2-m) \left( \frac{\Gamma(m(l+1))}{(m(l+1)-1)\Gamma(ml)} \right. \\
& \left. + \frac{\Gamma(m(l+2))}{(m(l+2)-1)\Gamma(m(l+1))} \right) \\
& + (1-m) \left[ \left( \frac{\Gamma(m(l+1))\Gamma(1+m)}{(m(l+1)-1)\Gamma(m(l+2))} \right. \right. \\
& \left. \left. + \frac{\Gamma(m(l+2))\Gamma(1+m)}{(m(l+2)-1)\Gamma(m(l+3))} \right) \right] \\
& -2 \left( \frac{1}{(l+1)m-1} + \frac{1}{(l+2)m-1} \right).
\end{aligned}
\tag{B19}
\end{equation}
$$

Note that, as stated above, the equation theoretically should
apply only for the case of $m > \frac{1}{1+l}$, which is quite restrictive.
However, thanks to the analyticity of the functions involved
in the expression, this equality remains valid for $m < \frac{1}{1+l}$ and
can be considered for any value of $m \in [0, 1]$ provided that
$m \neq \frac{1}{1+l}$ or $m \neq \frac{1}{2+l}$. Equation (B19) demonstrates Eq. (48).
The combination of these equations with the capillary model
Eq. (37) leads to Eq. (49).
*Code availability.* Note all computations were done using Scilab free software. The scripts for the computation of Eqs. (28)–(31) for the computation of WR and HC functions, Eqs. (38)–(41) for the computation of the dimensionless diffusivity, and Eqs. (56) for the computation of the $c_p$ parameter can be downloaded online at https://doi.org/10.5281/zenodo.4587160 (Lassabatere, 2021). Note also that the computation of the proposed sorptivity was implemented in the SoilWater-ToolBox software that interrelates specific modules to derive the soil hydraulic parameters using a wide range of cost-effective methods, accessible online at https://github.com/manaakiwhenua/SoilWater ToolBox/ (last access: 7 September 2021, Pollacco and Fernández-Gálvez, 2021) (open source under the GP-3.0 License).

*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 performed the analytical developments and verified the numerical results. DY, BL, DMF, SdP, and MR verified parts of the developments and numerical computations. JP and JFG helped with the use of the Kosugi model. RDS and MAN helped with the editing and the presentation of the results, in particular the writing of the application of the proposed approach for end users. All the authors contributed to the editing of the manuscript.

*Competing interests.* The authors declare that they have no conflict of interest.

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

*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. Hydrol., 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/sssaj2013.06.0231, 2014.

Braud, I.: Use of scaled forms of the infiltration equation for the estimation of unsaturated soil hydraulic properties (the Beerkan method), Eur. J. Soil Sci., 56, 361–374, https://doi.org/10.1111/j.1365-2389.2004.00660.x, 2005.

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

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: Encyclopaedia of Earth Sciences Series, Springer, Dordrecht, 824–826, https://doi.org/10.1007/978-90-481-3585-1_161, 2011.

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.

Fernández-Gálvez, J., Pollacco, J., Lilburne, L., McNeill, S., Carrick, S., Lassabatere, L., and Angulo-Jaramillo, R.: Deriving physical and unique bimodal soil Kosugi hydraulic parameters from inverse modelling, Adv. Water Resour., 153, 103933, https://doi.org/10.1016/j.advwatres.2021.103933, 2021.

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

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, San Diego, California, 1998.

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

Kosugi, K. and Hopmans, J.: Scaling water retention curves for soils with lognormal pore-size distribution, Soil Sci. Soc. Am. J., 62, 1496–1505, 1998.

Lassabatere, L.: Scilab script for sorptivity, Zenodo [code], https://doi.org/10.5281/zenodo.4587160, 2021.

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–523, 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.

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. R.: The theory of infiltration: 2. the profile at infinity, Soil Science, 83, 257–264, 1957.

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, 2013.

Pollacco, J. A. P. and Fernández-Gálvez, J.: SoilWater-ToolBox, available at: https://github.com/manaakiwhenua/SoilWater-ToolBox.jl, last access: 7 September 2021.

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.

Triadis, D. and Broadbridge, P.: The Green–Ampt limit with reference to infiltration coefficients, Water Resour. Res., 48, WR011747, https://doi.org/10.1029/2011WR011747, 2012.

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.

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.