A scaling procedure for straightforward computation of sorptivity

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

$$\begin{array}{}\text{(1)}& \left\{\begin{array}{l}I\left(t\right)=S\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{1}}\right)\sqrt{t}\\ S\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{1}}\right)=\underset{{\mathit{\theta }}_{\mathrm{0}}}{\overset{{\mathit{\theta }}_{\mathrm{1}}}{\int }}\mathit{\chi }\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\end{array}\right..\end{array}$$

In the above equations, S(θ₀,θ₁) stands for the sorptivity between θ₀ and θ₁, $\mathit{\chi }\left(\mathit{\theta }\right)=\frac{x(\mathit{\theta },t)}{\sqrt{t}}$ stands for the Boltzmann transformation variable, θ₀ is the initial water content, θ₁ 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:

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

where $D\left(\mathit{\theta }\right)=K\left(\mathit{\theta }\right)\frac{\mathrm{d}h}{\mathrm{d}\mathit{\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(θ(h)):

$$\begin{array}{}\text{(3)}& \begin{array}{rl}{S}_{\mathrm{K}}^{\mathrm{2}}\left(h_{\mathrm{0}},h_{\mathrm{1}}\right)& =\underset{h_{\mathrm{0}}}{\overset{h_{\mathrm{1}}}{\int }}\left(\mathit{\theta }\left(h_{\mathrm{1}}\right)+\mathit{\theta }\left(h\right)-\mathrm{2}\mathit{\theta }\left(h_{\mathrm{0}}\right)\right)K\left(h\right)\mathrm{d}h\\ & =\underset{h_{\mathrm{0}}}{\overset{h_{\mathrm{1}}}{\int }}\left({\mathit{\theta }}_{\mathrm{1}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)K\left(h\right)\mathrm{d}h,\end{array}\end{array}$$

where h is the water pressure head, h₀ and h₁ are respectively the initial and final water pressure heads, and θ₀=θ(h₀) and θ₁=θ(h₁). For the sake of clarity, the functions $S_{\mathrm{D}}^{\mathrm{2}}$ and $S_{\mathrm{K}}^{\mathrm{2}}$ are respectively referred to as the “diffusivity” and “conductivity” forms of sorptivity. $S_{\mathrm{D}}^{\mathrm{2}}$ and $S_{\mathrm{K}}^{\mathrm{2}}$ are equivalent so long as the water retention function θ(h) is bijective over the water pressure head interval [h₀,h₁], which is the case when the water pressure head at the surface is lower than the air-entry water pressure head, i.e., h₁≤h_a. Then, Eq. (3) can be easily deduced from Eq. (2) by a simple change of variable from θ to h:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}{S}_{\mathrm{K}}^{\mathrm{2}}\left(h_{\mathrm{0}},h_{\mathrm{1}}\le h_{\mathrm{a}}\right)& =S_{\mathrm{D}}^{\mathrm{2}}\left(\mathit{\theta }\left(h_{\mathrm{0}}\right),\mathit{\theta }\left(h_{\mathrm{1}}\right)\right)\\ & =S_{\mathrm{D}}^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{1}}\right).\end{array}\end{array}$$

Otherwise, when the surface water pressure head exceeds the soil air-entry pressure, i.e., $h_{\mathrm{1}}>h_{\mathrm{a}}(\le \mathrm{0})$, sorptivity must be computed using Eq. (3) (Ross et al., 1996). Indeed, the function θ(h) is no longer bijective over the interval [h₀,h₁]. 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₀,h_a] ensuring the bijectivity of the function θ(h) and the other part retaining the integration interval [h_a,h₁] corresponding to the saturated part of the integration (Ross et al., 1996). Then, the change of variable from θ to h in the first integral leads to the retrieval of $S_{\mathrm{D}}^{\mathrm{2}}$:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{S}_{\mathrm{K}}^{\mathrm{2}}\left(h_{\mathrm{0}},h_{\mathrm{1}}>h_{\mathrm{a}}\right)& =\underset{h_{\mathrm{0}}}{\overset{h_{\mathrm{a}}}{\int }}\left({\mathit{\theta }}_{\mathrm{s}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)K\left(\mathit{\theta }\right(h\left)\right)\mathrm{d}h\\ & +\underset{h_{\mathrm{a}}}{\overset{h_{\mathrm{1}}}{\int }}\left({\mathit{\theta }}_{\mathrm{s}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)K\left(\mathit{\theta }\right(h\left)\right)\mathrm{d}h\\ & =\underset{h_{\mathrm{0}}}{\overset{h_{\mathrm{a}}}{\int }}\left({\mathit{\theta }}_{\mathrm{s}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)K\left(\mathit{\theta }\right(h\left)\right)\mathrm{d}h\\ & +\mathrm{2}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}\right)K_{\mathrm{s}}\underset{h_{\mathrm{a}}}{\overset{h_{\mathrm{1}}}{\int }}\mathrm{d}h\\ & =\underset{h_{\mathrm{0}}}{\overset{h_{\mathrm{a}}}{\int }}\left({\mathit{\theta }}_{\mathrm{s}}+\mathit{\theta }\left(h\right)-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)K\left(\mathit{\theta }\right(h\left)\right)\mathrm{d}h\\ & +\mathrm{2}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}\right)K_{\mathrm{s}}\left(h_{\mathrm{1}}-h_{\mathrm{a}}\right)\\ & =\underset{{\mathit{\theta }}_{\mathrm{0}}}{\overset{{\mathit{\theta }}_{\mathrm{s}}}{\int }}\left({\mathit{\theta }}_{\mathrm{s}}+\mathit{\theta }-\mathrm{2}{\mathit{\theta }}_{\mathrm{0}}\right)D\left(\mathit{\theta }\right)\mathrm{d}\mathit{\theta }\\ & +\mathrm{2}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}\right)K_{\mathrm{s}}\left(h_{\mathrm{1}}-h_{\mathrm{a}}\right)\\ & =S_{\mathrm{D}}^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{s}}\right)+\mathrm{2}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}\right)K_{\mathrm{s}}\left(h_{\mathrm{1}}-h_{\mathrm{a}}\right).\end{array}\end{array}$$

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 $\left(m=\mathrm{1}-\frac{\mathrm{2}}{n}\right)$ 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 $\left(m=\mathrm{1}-\frac{\mathrm{1}}{n}\right)$ for the WR function and the Mualem (1976) capillary model for the HC function. The vGM 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:

$$\begin{array}{}\text{(6)}& \left\{\begin{array}{l}{\mathit{\theta }}_{\mathrm{d}}\left(h\right)={\mathit{\theta }}_{\mathrm{r}}+\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)H\left(\mathrm{1}+\frac{h}{\left|h_{\mathrm{d}}\right|}\right)\\ K_{\mathrm{d}}\left(\mathit{\theta }\right)=K_{\mathrm{s}}H\left(\frac{\mathit{\theta }-{\mathit{\theta }}_{\mathrm{r}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}-\mathrm{1}\right)\end{array}\right.\end{array}$$

BC model:

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

vGB model:

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

vGM model:

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

KG model:

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

where H stands for the one-sided Heaviside step function: $H(x<\mathrm{0})=\mathrm{0}$, $H(x\ge \mathrm{0})=\mathrm{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, θ_r, saturated water content, θ_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 indetermination 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^{*\mathrm{2}}(-\mathrm{\infty },\mathrm{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., $(-\mathrm{\infty },\mathrm{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_{\mathrm{K}}^{\mathrm{2}}(h_{\mathrm{0}},h_{\mathrm{1}})$ to the maximum square sorptivity $S_{\mathrm{K}}^{\mathrm{2}}(-\mathrm{\infty },\mathrm{0})$ by isolating the effect of initial and final conditions (h₀ and h₁) and (ii) scaling hydraulic functions and sorptivity to split the contributions of shape and scale parameters and facilitate the computation of $S_{\mathrm{K}}^{\mathrm{2}}(h_{\mathrm{0}},h_{\mathrm{1}})$. In our case, we consider that the final conditions always involve a null water pressure head, h₁=0, and saturated conditions, θ₁=θ_s, while initial conditions (h₀,θ₀) 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 (${\mathit{\theta }}_{\mathrm{0}}\le \frac{\mathrm{1}}{\mathrm{4}}{\mathit{\theta }}_{\mathrm{s}}$) and a zero water pressure head at surface (i.e., h₁=0 and θ₁=θ_s), the following approximation applied:

$$\begin{array}{}\text{(11)}& \left\{\begin{array}{l}{S}^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{s}}\right)\approx S^{\mathrm{2}}\left(\mathrm{0},{\mathit{\theta }}_{\mathrm{s}}\right)\frac{K_{\mathrm{s}}-K_{\mathrm{0}}}{K_{\mathrm{s}}}\frac{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}}{{\mathit{\theta }}_{\mathrm{s}}}\\ S^{\mathrm{2}}\left(\mathrm{0},{\mathit{\theta }}_{\mathrm{s}}\right)=c_{\mathrm{p}}\left|h_{\mathrm{g}}\right|K_{\mathrm{s}}{\mathit{\theta }}_{\mathrm{s}}\end{array}\right.,\end{array}$$

where K₀ corresponds to the initial hydraulic conductivity K₀=K(θ₀). The parameter $c_{\mathrm{p}}=\frac{S^{\mathrm{2}}(\mathrm{0},{\mathit{\theta }}_{\mathrm{s}})}{\left|h_{\mathrm{g}}\right|K_{\mathrm{s}}{\mathit{\theta }}_{\mathrm{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²(θ₀,θ_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, $\left|h_{\mathrm{g}}\right|$, K_s, and θ_s, and initial conditions, θ₀:

$$\begin{array}{}\text{(12)}& S^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{s}}\right)\approx c_{\mathrm{p}}\left|h_{\mathrm{g}}\right|K_{\mathrm{s}}{\mathit{\theta }}_{\mathrm{s}}{\displaystyle \frac{K_{\mathrm{s}}-K_{\mathrm{0}}}{K_{\mathrm{s}}}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}}{{\mathit{\theta }}_{\mathrm{s}}}}.\end{array}$$

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₁=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:

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{S^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{s}}\right)}{S^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{r}},{\mathit{\theta }}_{\mathrm{s}}\right)}}\approx {\displaystyle \frac{K_{\mathrm{s}}-K_{\mathrm{0}}}{K_{\mathrm{s}}}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}.\end{array}$$

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

The following derivations are then proposed:

$$\begin{array}{}\text{(14)}& \begin{array}{rl}{S}_{\mathrm{K}}^{\mathrm{2}}\left(h_{\mathrm{0}},h_{\mathrm{1}}=\mathrm{0}\right)& =S_{\mathrm{D}}^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{0}},{\mathit{\theta }}_{\mathrm{s}}\right)+\mathrm{2}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}\right)K_{\mathrm{s}}\left(h_{\mathrm{1}}-h_{\mathrm{a}}\right)\\ & ={\displaystyle \frac{K_{\mathrm{s}}-K_{\mathrm{0}}}{K_{\mathrm{s}}}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}{S}_{\mathrm{D}}^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{r}},{\mathit{\theta }}_{\mathrm{s}}\right)\\ & +\mathrm{2}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}\right)K_{\mathrm{s}}\left|h_{\mathrm{a}}\right|\mathrm{since}{h}_{\mathrm{1}}=\mathrm{0}\\ & ={\displaystyle \frac{K_{\mathrm{s}}-K_{\mathrm{0}}}{K_{\mathrm{s}}}}{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}{S}_{\mathrm{D}}^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{r}},{\mathit{\theta }}_{\mathrm{s}}\right)\\ & +{\displaystyle \frac{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}}\left(\mathrm{2}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)K_{\mathrm{s}}\left|h_{\mathrm{a}}\right|\right),\end{array}\end{array}$$

where $S_{\mathrm{D}}^{\mathrm{2}}({\mathit{\theta }}_{\mathrm{r}},{\mathit{\theta }}_{\mathrm{s}})$ and $\mathrm{2}({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}})K_{\mathrm{s}}\left|h_{\mathrm{a}}\right|$ refer to the unsaturated and saturated parts of the square sorptivity $S_{\mathrm{K}}^{\mathrm{2}}(-\mathrm{\infty },\mathrm{0})$ (see Eq. 5 with $h_{\mathrm{0}}=-\mathrm{\infty }$, i.e., θ₀=θ_r, and h₁=0). Equation (14) can be simplified by introducing correcting factors R_θ and R_K to account for the contribution of the initial conditions:

$$\begin{array}{}\text{(15)}& S_{\mathrm{K}}^{\mathrm{2}}\left(h_{\mathrm{0}},\mathrm{0}\right)=R_{\mathrm{K}}{R}_{\mathit{\theta }}{S}_{\mathrm{D}}^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{r}},{\mathit{\theta }}_{\mathrm{s}}\right)+R_{\mathit{\theta }}\left(\mathrm{2}\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)K_{\mathrm{s}}\left|h_{\mathrm{a}}\right|\right),\end{array}$$

with R_θ and R_K defined as follows:

$$\begin{array}{}\text{(16)}& \left\{\begin{array}{l}{R}_{\mathit{\theta }}=\frac{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}\\ R_{\mathrm{K}}={\displaystyle \frac{K_{\mathrm{s}}-K_{\mathrm{0}}}{K_{\mathrm{s}}}}.\end{array}\right.\end{array}$$

2.1.2 Scaling sorptivity and defining parameters c_p and $c_{\mathrm{p}}^{\prime }$

So far, the computation of the square sorptivity $S_{\mathrm{K}}^{\mathrm{2}}(h_{\mathrm{0}},\mathrm{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_{\mathrm{K}}^{\mathrm{2}}(h_{\mathrm{0}},\mathrm{0})$ can be easily related to the square scaled sorptivity $S_{\mathrm{K}}^{*\mathrm{2}}(h_{\mathrm{0}}^{*},\mathrm{0})$ by scaling variables (water content, water pressure head, and hydraulic conductivity) as follows (Ross et al., 1996):

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

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^{*}\left(S_{\mathrm{e}}\right)=K_{\mathrm{r}}\left(S_{\mathrm{e}}\right)\frac{\mathrm{d}{h}^{*}}{\mathrm{d}{S}_{\mathrm{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, θ_s=1, h_g=1, and K_s=1, except the residual water content that is fixed at zero, θ_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², to the square scaled sorptivity (unit soil), S^*2 (Ross et al., 1996):

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

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_{\mathrm{0}}^{*}$ and $h_{\mathrm{1}}^{*}$, or saturation degrees, $S_{\mathrm{e},\mathrm{0}}=S_{\mathrm{e}}\left(h_{\mathrm{0}}^{*}\right)$ and $S_{\mathrm{e},\mathrm{1}}=S_{\mathrm{e}}\left(h_{\mathrm{1}}^{*}\right)$:

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

$S_{\mathrm{D}}^{*\mathrm{2}}$ and $S_{\mathrm{K}}^{*\mathrm{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

$$\begin{array}{}\text{(20)}& \left\{\begin{array}{l}{S}_{\mathrm{K}}^{*\mathrm{2}}\left(h_{\mathrm{0}}^{*},h_{\mathrm{1}}^{*}\le h_{\mathrm{a}}^{*}\right)=S_{\mathrm{D}}^{*\mathrm{2}}\left(S_{\mathrm{e},\mathrm{0}},S_{\mathrm{e},\mathrm{1}}\right)\\ S_{\mathrm{K}}^{*\mathrm{2}}\left(h_{\mathrm{0}}^{*},h_{\mathrm{1}}^{*}>h_{\mathrm{a}}^{*}\right)=S_{\mathrm{D}}^{*\mathrm{2}}\left(S_{\mathrm{e},\mathrm{0}},\mathrm{1}\right)+\mathrm{2}\left(\mathrm{1}-S_{\mathrm{e},\mathrm{0}}\right)\left(h_{\mathrm{1}}^{*}-h_{\mathrm{a}}^{*}\right)\end{array}\right.\end{array}$$

with $h_{\mathrm{a}}^{*}=\frac{h_{\mathrm{a}}}{\left|h_{\mathrm{g}}\right|}$. The scaled sorptivity functions, $S_{\mathrm{D}}^{*\mathrm{2}}$ and $S_{\mathrm{K}}^{*\mathrm{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_{\mathrm{K}}^{\mathrm{2}}(-\mathrm{\infty },\mathrm{0})$ and its unsaturated part $S_{\mathrm{D}}^{\mathrm{2}}({\mathit{\theta }}_{\mathrm{r}},{\mathit{\theta }}_{\mathrm{s}})$, which are involved in Eq. (15), to define the parameters c_p and $c_{\mathrm{p}}^{\prime }$:

$$\begin{array}{}\text{(21)}& \left\{\begin{array}{l}{c}_{\mathrm{p}}=S_{\mathrm{K}}^{*\mathrm{2}}(-\mathrm{\infty },\mathrm{0})=\frac{S_{\mathrm{K}}^{\mathrm{2}}(-\mathrm{\infty },\mathrm{0})}{\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)K_{\mathrm{s}}\left|h_{\mathrm{g}}\right|}\\ c_{\mathrm{p}}^{\prime }=S_{\mathrm{D}}^{*\mathrm{2}}(\mathrm{0},\mathrm{1})=\frac{S_{\mathrm{D}}^{\mathrm{2}}\left({\mathit{\theta }}_{\mathrm{r}},{\mathit{\theta }}_{\mathrm{s}}\right)}{\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)K_{\mathrm{s}}\left|h_{\mathrm{g}}\right|}.\end{array}\right.\end{array}$$

The application of Eq. (19) with the proper lower and upper initial and final conditions ($S_{\mathrm{e},\mathrm{0}}=\mathrm{0}$, $h_{\mathrm{0}}=-\mathrm{\infty },S_{\mathrm{e},\mathrm{1}}=\mathrm{1}$, h₁=0) leads to the following expression for the two parameters c_p and $c_{\mathrm{p}}^{\prime }$:

$$\begin{array}{}\text{(22)}& \left\{\begin{array}{l}{c}_{\mathrm{p}}=\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}\left(\mathrm{1}+S_{\mathrm{e}}\left(h^{*}\right)\right)K_{\mathrm{r}}\left(h^{*}\right)\mathrm{d}{h}^{*}\\ c_{\mathrm{p}}^{\prime }=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(\mathrm{1}+S_{\mathrm{e}}\right)D^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}.\end{array}\right.\end{array}$$

The parameters c_p and $c_{\mathrm{p}}^{\prime }$ depend exclusively on the soil shape parameters. The application of Eq. (20) provides a very simple equation linking the two parameters:

$$\begin{array}{}\text{(23)}& c_{\mathrm{p}}=c_{\mathrm{p}}^{\prime }+\mathrm{2}\left|h_{\mathrm{a}}^{*}\right|.\end{array}$$

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_{\mathrm{K}}^{*\mathrm{2}}(h_{\mathrm{0}},\mathrm{0})$:

$$\begin{array}{}\text{(24)}& \begin{array}{rl}{S}_{\mathrm{K}}^{*\mathrm{2}}\left(h_{\mathrm{0}}^{*},\mathrm{0}\right)& =c_{\mathrm{p}}^{\prime }{R}_{\mathrm{K}}{R}_{\mathit{\theta }}+\mathrm{2}\left|h_{\mathrm{a}}^{*}\right|R_{\mathit{\theta }}\\ & =\left(c_{\mathrm{p}}-\mathrm{2}\left|h_{\mathrm{a}}^{*}\right|\right)R_{\mathrm{K}}{R}_{\mathit{\theta }}+\mathrm{2}\left|h_{\mathrm{a}}^{*}\right|R_{\mathit{\theta }},\end{array}\end{array}$$

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:

$$\begin{array}{}\text{(25)}& \left\{\begin{array}{l}{R}_{\mathit{\theta }}=\frac{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}=\mathrm{1}-S_{\mathrm{e},\mathrm{0}}\\ R_{\mathrm{K}}={\displaystyle \frac{K_{\mathrm{s}}-K_{\mathrm{0}}}{K_{\mathrm{s}}}}=\mathrm{1}-K_{\mathrm{r}}\left(S_{\mathrm{e},\mathrm{0}}\right).\end{array}\right.\end{array}$$

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_{\mathrm{K}}^{\mathrm{2}}(h_{\mathrm{0}},\mathrm{0})$:

$$\begin{array}{}\text{(26)}& \left\{\begin{array}{l}{S}_{\mathrm{K}}^{\mathrm{2}}\left(h_{\mathrm{0}},\mathrm{0}\right)=S_{\mathrm{K}}^{*\mathrm{2}}\left(h_{\mathrm{0}}^{*},\mathrm{0}\right)\left({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}\right)K_{\mathrm{s}}\left|h_{\mathrm{g}}\right|\\ S_{\mathrm{K}}^{*\mathrm{2}}\left(h_{\mathrm{0}}^{*},\mathrm{0}\right)=\left(c_{\mathrm{p}}-\mathrm{2}\left|h_{\mathrm{a}}^{*}\right|\right)R_{\mathrm{K}}{R}_{\mathit{\theta }}+\mathrm{2}\left|h_{\mathrm{a}}^{*}\right|R_{\mathit{\theta }}\\ R_{\mathit{\theta }}=\mathrm{1}-S_{\mathrm{e},\mathrm{0}}\\ R_{\mathrm{K}}=\mathrm{1}-K_{\mathrm{r}}\left(S_{\mathrm{e},\mathrm{0}}\right).\end{array}\right.\end{array}$$

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:

$$\begin{array}{}\text{(27)}& \left\{\begin{array}{l}{S}_{\mathrm{e},\mathrm{d}}\left(h^{*}\right)=H\left(\mathrm{1}+h^{*}\right)\\ K_{\mathrm{r},\mathrm{d}}\left(S_{\mathrm{e}}\right)=H\left(S_{\mathrm{e}}-\mathrm{1}\right)\end{array}\right.\end{array}$$

BC model:

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

vGB model:

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

vGM model:

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

KG model:

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

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^{*}\left(S_{\mathrm{e}}\right)=K_{\mathrm{r}}\left(S_{\mathrm{e}}\right)\frac{\mathrm{d}{h}^{*}}{\mathrm{d}{S}_{\mathrm{e}}}$ (see Appendix A):

$$\begin{array}{}\text{(32)}& {\displaystyle }{D}_{\mathrm{d}}^{*}\left(S_{\mathrm{e}}\right)=\mathit{\delta }\left(S_{\mathrm{e}}\right)\text{(33)}& {\displaystyle }{D}_{\mathrm{BC}}^{*}\left(S_{\mathrm{e}}\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{BC}}}}{S}_{\mathrm{e}}^{{\mathit{\eta }}_{\mathrm{BC}}-\left(\frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{BC}}}+\mathrm{1}\right)}\text{(34)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{vGB}}^{*}\left(S_{\mathrm{e}}\right)=& {\displaystyle \frac{\mathrm{1}-m_{\mathrm{vGB}}}{\mathrm{2}{m}_{\mathrm{vGB}}}}{S}_{\mathrm{e}}^{{\mathit{\eta }}_{\mathrm{vGB}}-\frac{\mathrm{1}+m_{\mathrm{vGB}}}{\mathrm{2}{m}_{\mathrm{vGB}}}}\\ & {\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m_{\mathrm{vGB}}}}\right)}^{-\frac{\mathrm{1}+m_{\mathrm{vGB}}}{\mathrm{2}}}\end{array}\text{(35)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{vGM}}^{*}\left(S_{\mathrm{e}}\right)=& {\displaystyle \frac{\mathrm{1}-m_{\mathrm{vGM}}}{m_{\mathrm{vGM}}}}{S}_{\mathrm{e}}^{l_{\mathrm{vGM}}-\frac{\mathrm{1}}{m_{\mathrm{vGM}}}}\\ & \left({\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m_{\mathrm{vGM}}}}\right)}^{-m_{\mathrm{vGM}}}\right.\\ & \left.+{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m_{\mathrm{vGM}}}}\right)}^{m_{\mathrm{vGM}}}-\mathrm{2}\right)\end{array}\text{(36)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{KG}}^{*}\left(S_{\mathrm{e}}\right)=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sqrt{{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}}{\mathit{\sigma }}_{\mathrm{KG}}{S}_{\mathrm{e}}^{l_{\mathrm{KG}}}\left(\mathrm{erfc}\left({\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}{S}_{\mathrm{e}}\right)\right.\right.\\ & {\left.\left.+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{KG}}}{\sqrt{\mathrm{2}}}}\right)\right)}^{\mathrm{2}}{e}^{{\left({\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}{S}_{\mathrm{e}}\right)\right)}^{\mathrm{2}}+\sqrt{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{KG}}{\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}{S}_{\mathrm{e}}\right)},\end{array}\end{array}$$

where δ 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, η, with the combination of those of the WR function, λ=mn, as follows:

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

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_{\mathrm{vGM}}=l_{\mathrm{KG}}=\mathrm{1}/\mathrm{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

$$\begin{array}{}\text{(38)}& {\displaystyle }{D}_{\mathrm{BC}}^{*}\left(S_{\mathrm{e}}\right)={\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{BC}}}}{S}_{\mathrm{e}}^{\frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{BC}}}+\mathrm{2}}\text{(39)}& {\displaystyle }{D}_{\mathrm{vGB}}^{*}\left(S_{\mathrm{e}}\right)={\displaystyle \frac{\mathrm{1}-m_{\mathrm{vGB}}}{\mathrm{2}{m}_{\mathrm{vGB}}}}{S}_{\mathrm{e}}^{\frac{\mathrm{1}+\mathrm{3}{m}_{\mathrm{vGB}}}{\mathrm{2}{m}_{\mathrm{vGB}}}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m_{\mathrm{vGB}}}}\right)}^{-\frac{\mathrm{1}+m_{\mathrm{vGB}}}{\mathrm{2}}}\text{(40)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{vGM}}^{*}\left(S_{\mathrm{e}}\right)& ={\displaystyle \frac{\mathrm{1}-m_{\mathrm{vGM}}}{m_{\mathrm{vGM}}}}{S}_{\mathrm{e}}^{\frac{m_{\mathrm{vGM}}-\mathrm{2}}{\mathrm{2}{m}_{\mathrm{vGM}}}}\left({\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m_{\mathrm{vGM}}}}\right)}^{-m_{\mathrm{vGM}}}\right.\\ & \left.+{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m_{\mathrm{vGM}}}}\right)}^{m_{\mathrm{vGM}}}-\mathrm{2}\right)\end{array}\text{(41)}& {\displaystyle }\begin{array}{rl}{D}_{\mathrm{KG}}^{*}\left(S_{\mathrm{e}}\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sqrt{{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}}{\mathit{\sigma }}_{\mathrm{KG}}{S}_{\mathrm{e}}^{\frac{\mathrm{1}}{\mathrm{2}}}\left(\mathrm{erfc}\left({\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}{S}_{\mathrm{e}}\right)\right.\right.\\ & {\left.\left.+{\displaystyle \frac{{\mathit{\sigma }}_{\mathrm{KG}}}{\sqrt{\mathrm{2}}}}\right)\right)}^{\mathrm{2}}{e}^{{\left({\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}{S}_{\mathrm{e}}\right)\right)}^{\mathrm{2}}+\sqrt{\mathrm{2}}{\mathit{\sigma }}_{\mathrm{KG}}{\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}{S}_{\mathrm{e}}\right)}.\end{array}\end{array}$$

The set of equations (Eqs. 38–41) shows that each hydraulic diffusivity function involves only one shape parameter, i.e., λ_BC, m_vGB, m_vGM, and σ_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^{*}<-\mathrm{1}$ and a unit value K_r=1 for $h^{*}\ge -\mathrm{1}$, as featured by Eq. (27). In this case, Eq. (22) yields

$$\begin{array}{}\text{(42)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{d}}& =\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}\left(\mathrm{1}+S_{\mathrm{e},\mathrm{d}}\left(h^{*}\right)\right)K_{\mathrm{r},\mathrm{d}}\left(h^{*}\right)\mathrm{d}{h}^{*}\\ & =\underset{-\mathrm{\infty }}{\overset{-\mathrm{1}}{\int }}(\mathrm{1}+\mathrm{0})\cdot \mathrm{0}\cdot \mathrm{d}{h}^{*}+\underset{-\mathrm{1}}{\overset{\mathrm{0}}{\int }}(\mathrm{1}+\mathrm{1})\cdot \mathrm{1}\cdot \mathrm{d}{h}^{*}\\ & =\mathrm{2}.\end{array}\end{array}$$

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_{\mathrm{a}}^{*}=-\mathrm{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_{\mathrm{BC}}^{*}\left(S_{\mathrm{e}}\right)$ 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):

$$\begin{array}{}\text{(43)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{BC}}\left({\mathit{\lambda }}_{\mathrm{BC}},{\mathit{\eta }}_{\mathrm{BC}}\right)& =\mathrm{2}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{BC}}{\mathit{\eta }}_{\mathrm{BC}}-\mathrm{1}}}\\ & +{\displaystyle \frac{\mathrm{1}}{{\mathit{\lambda }}_{\mathrm{BC}}{\mathit{\eta }}_{\mathrm{BC}}+{\mathit{\lambda }}_{\mathrm{BC}}-\mathrm{1}}}.\end{array}\end{array}$$

When the relation between η_BC and λ_BC is ruled by Eq. (37) with p=1, the analytical expression of c_p turns into

$$\begin{array}{}\text{(44)}& c_{\mathrm{p},\mathrm{BC}}\left({\mathit{\lambda }}_{\mathrm{BC}}\right)=\mathrm{2}+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}{\mathit{\lambda }}_{\mathrm{BC}}+\mathrm{1}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}{\mathit{\lambda }}_{\mathrm{BC}}+\mathrm{1}}}.\end{array}$$

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_{\mathrm{a}}^{*}=\mathrm{0}$. Equation (22) shows that c_p reverts to the diffusivity form of the square scaled sorptivity, $\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}(\mathrm{1}+S_{\mathrm{e}})D^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}$. In addition, the diffusivity function, $D_{\mathrm{vGB}}^{*}\left(S_{\mathrm{e}}\right)$ (Eq. 34), makes the square scaled sorptivity analytically integrable, leading to (see demonstration in Appendix B2)

$$\begin{array}{}\text{(45)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{vGB}}& \left(m_{\mathrm{vGB}},n_{\mathrm{vGB}},{\mathit{\eta }}_{\mathrm{vGB}}\right)=\mathrm{\Gamma }\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{n_{\mathrm{vGB}}}}\right)\\ & \left[{\displaystyle \frac{\mathrm{\Gamma }\left(m_{\mathrm{vGB}}{\mathit{\eta }}_{\mathrm{vGB}}-\frac{\mathrm{1}}{n_{\mathrm{vGB}}}\right)}{\mathrm{\Gamma }\left(m_{\mathrm{vGB}}{\mathit{\eta }}_{\mathrm{vGB}}\right)}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }\left(m_{\mathrm{vGB}}{\mathit{\eta }}_{\mathrm{vGB}}+m_{\mathrm{vGB}}-\frac{\mathrm{1}}{n_{\mathrm{vGB}}}\right)}{\mathrm{\Gamma }\left(m_{\mathrm{vGB}}{\mathit{\eta }}_{\mathrm{vGB}}+m_{\mathrm{vGB}}\right)}}\right],\end{array}\end{array}$$

where Γ is the gamma function:

$$\begin{array}{}\text{(46)}& \mathrm{\Gamma }\left(z\right)=\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}{t}^{z-\mathrm{1}}{e}^{-t}\mathrm{d}t.\end{array}$$

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

$$\begin{array}{}\text{(47)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{vGB}}\left(m_{\mathrm{vGB}}\right)& =\mathrm{\Gamma }\left({\displaystyle \frac{\mathrm{3}-m_{\mathrm{vGB}}}{\mathrm{2}}}\right)\left[{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\mathrm{1}+\mathrm{5}{m}_{\mathrm{vGB}}}{\mathrm{2}}\right)}{\mathrm{\Gamma }\left(\mathrm{1}+\mathrm{2}{m}_{\mathrm{vGB}}\right)}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\mathrm{1}+\mathrm{7}{m}_{\mathrm{vGB}}}{\mathrm{2}}\right)}{\mathrm{\Gamma }\left(\mathrm{1}+\mathrm{3}{m}_{\mathrm{vGB}}\right)}}\right].\end{array}\end{array}$$

The expression corresponding to Eq. (45) was already proposed and discussed by Haverkamp et al. (2005).

2.3.4 c_p 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_p (see demonstration in Appendix B3):

$$\begin{array}{}\text{(48)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{vGM}}& \left(m_{\mathrm{vGM}},l_{\mathrm{vGM}}\right)=\mathrm{\Gamma }\left(\mathrm{2}-m_{\mathrm{vGM}}\right)\\ & \left({\displaystyle \frac{\mathrm{\Gamma }\left(m_{\mathrm{vGM}}\left(\mathrm{1}+l_{\mathrm{vGM}}\right)\right)}{\left(m_{\mathrm{vGM}}\left(\mathrm{1}+l_{\mathrm{vGM}}\right)-\mathrm{1}\right)\mathrm{\Gamma }\left(m_{\mathrm{vGM}}{l}_{\mathrm{vGM}}\right)}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }\left(m_{\mathrm{vGM}}\left(\mathrm{2}+l_{\mathrm{vGM}}\right)\right)}{\left(m_{\mathrm{vGM}}\left(\mathrm{2}+l_{\mathrm{vGM}}\right)-\mathrm{1}\right)\mathrm{\Gamma }\left(m_{\mathrm{vGM}}\left(\mathrm{1}+l_{\mathrm{vGM}}\right)\right)}}\right)\\ & +\left(\mathrm{1}-m_{\mathrm{vGM}}\right)\\ & \left[\left({\displaystyle \frac{\mathrm{\Gamma }\left(m_{\mathrm{vGM}}\left(\mathrm{1}+l_{\mathrm{vGM}}\right)\right)\mathrm{\Gamma }\left(\mathrm{1}+m_{\mathrm{vGM}}\right)}{\left(m_{\mathrm{vGM}}\left(\mathrm{1}+l_{\mathrm{vGM}}\right)-\mathrm{1}\right)\mathrm{\Gamma }\left(m_{\mathrm{vGM}}\left(\mathrm{2}+l_{\mathrm{vGM}}\right)\right)}}\right.\right.\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }\left(m_{\mathrm{vGM}}\left(\mathrm{2}+l_{\mathrm{vGM}}\right)\right)\mathrm{\Gamma }\left(\mathrm{1}+m_{\mathrm{vGM}}\right)}{\left(m_{\mathrm{vGM}}\left(\mathrm{2}+l_{\mathrm{vGM}}\right)-\mathrm{1}\right)\mathrm{\Gamma }\left(m_{\mathrm{vGM}}\left(\mathrm{3}+l_{\mathrm{vGM}}\right)\right)}}\right)\\ & \left.-\mathrm{2}\left({\displaystyle \frac{\mathrm{1}}{m_{\mathrm{vGM}}\left(\mathrm{1}+l_{\mathrm{vGM}}\right)-\mathrm{1}}}+{\displaystyle \frac{\mathrm{1}}{m_{\mathrm{vGM}}\left(\mathrm{2}+l_{\mathrm{vGM}}\right)-\mathrm{1}}}\right)\right].\end{array}\end{array}$$

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

$$\begin{array}{}\text{(49)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{vGM}}& \left(m_{\mathrm{vGM}}\right)=\mathrm{\Gamma }\left(\mathrm{2}-m_{\mathrm{vGM}}\right)\\ & \left({\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\mathrm{3}{m}_{\mathrm{vGM}}}{\mathrm{2}}\right)}{\left(\frac{\mathrm{3}{m}_{\mathrm{vGM}}}{\mathrm{2}}-\mathrm{1}\right)\mathrm{\Gamma }\left(\frac{m_{\mathrm{vGM}}}{\mathrm{2}}\right)}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\mathrm{5}{m}_{\mathrm{vGM}}}{\mathrm{2}}\right)}{\left(\frac{\mathrm{5}{m}_{\mathrm{vGM}}}{\mathrm{2}}-\mathrm{1}\right)\mathrm{\Gamma }\left(\frac{\mathrm{3}{m}_{\mathrm{vGM}}}{\mathrm{2}}\right)}}\right)\\ & +\left(\mathrm{1}-m_{\mathrm{vGM}}\right)\left[\left({\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\mathrm{3}{m}_{\mathrm{vGM}}}{\mathrm{2}}\right)\mathrm{\Gamma }\left(\mathrm{1}+m_{\mathrm{vGM}}\right)}{\left(\frac{\mathrm{3}{m}_{\mathrm{vGM}}}{\mathrm{2}}-\mathrm{1}\right)\mathrm{\Gamma }\left(\frac{\mathrm{5}{m}_{\mathrm{vGM}}}{\mathrm{2}}\right)}}\right.\right.\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\mathrm{5}{m}_{\mathrm{vGM}}}{\mathrm{2}}\right)\mathrm{\Gamma }\left(\mathrm{1}+m_{\mathrm{vGM}}\right)}{\left(\frac{\mathrm{5}{m}_{\mathrm{vGM}}}{\mathrm{2}}-\mathrm{1}\right)\mathrm{\Gamma }\left(\frac{\mathrm{7}{m}_{\mathrm{vGM}}}{\mathrm{2}}\right)}}\right)\\ & \left.-\mathrm{2}\left({\displaystyle \frac{\mathrm{1}}{\frac{\mathrm{3}{m}_{\mathrm{vGM}}}{\mathrm{2}}-\mathrm{1}}}+{\displaystyle \frac{\mathrm{1}}{\frac{\mathrm{5}{m}_{\mathrm{vGM}}}{\mathrm{2}}-\mathrm{1}}}\right)\right].\end{array}\end{array}$$

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_p 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{array}{}\text{(50)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{KG}}\left({\mathit{\sigma }}_{\mathrm{KG}},l_{\mathrm{KG}}\right)& =\underset{-\mathrm{\infty }}{\overset{\mathrm{0}}{\int }}\left(\mathrm{1}+S_{\mathrm{e},\mathrm{KG}}\left(h^{*}\right)\right)K_{\mathrm{r},\mathrm{KG}}\left(h^{*}\right)\mathrm{d}{h}^{*}\\ & =\underset{-\mathrm{\infty }}{\overset{h_{\mathrm{KG}}^{*}\left(\frac{\mathrm{1}}{\mathrm{2}}\right)}{\int }}\left(\mathrm{1}+S_{\mathrm{e},\mathrm{KG}}\left(h^{*}\right)\right)K_{\mathrm{r},\mathrm{KG}}\left(h^{*}\right)\mathrm{d}{h}^{*}\\ & +\underset{h_{\mathrm{KG}}^{*}\left(\frac{\mathrm{1}}{\mathrm{2}}\right)}{\overset{\mathrm{0}}{\int }}\left(\mathrm{1}+S_{\mathrm{e},\mathrm{KG}}\left(h^{*}\right)\right)K_{\mathrm{r},\mathrm{KG}}\left(h^{*}\right)\mathrm{d}{h}^{*}\\ & =\underset{\mathrm{0}}{\overset{\frac{\mathrm{1}}{\mathrm{2}}}{\int }}\left(\mathrm{1}+S_{\mathrm{e}}\right)D_{\mathrm{KG}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}\\ & +\underset{h_{\mathrm{KG}}^{*}\left(\frac{\mathrm{1}}{\mathrm{2}}\right)}{\overset{\mathrm{0}}{\int }}\left(\mathrm{1}+S_{\mathrm{e},\mathrm{KG}}\left(h^{*}\right)\right)K_{\mathrm{r},\mathrm{KG}}\left(h^{*}\right)\mathrm{d}{h}^{*},\end{array}\end{array}$$

where $h_{\mathrm{KG}}^{*}\left(\frac{\mathrm{1}}{\mathrm{2}}\right)$ is the water pressure head corresponding to $S_{\mathrm{e}}=\frac{\mathrm{1}}{\mathrm{2}}$. In the last expression of Eq. (50), c_p 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{\mathrm{1}}{\mathrm{2}}$.

2.4 Shape indexes for comparing c_p between the selected hydraulic models

The approach described below allows c_p to be determined for the selected BC, vGB, vGM, and KG models. In addition to its contribution to simplifying Eq. (26), c_p 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_p should not depend on the choice of the hydraulic models. We then investigate the dependency of c_p upon the selected hydraulic model. Consequently, we designed shape indexes to compare c_p 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_vGB and m_vGM to be the appropriate shape indexes, i.e., x_vGB=m_vGB and x_vGM=m_vGM. Next, we designed the shape index for BC, x_BC, deriving from vGB, m_vGB, considering that vGB and BC functions describe close WR curves. Indeed, for large values of n_vGB, the vGB model converges towards a power function similar to the BC model (Haverkamp et al., 2005):

$$\begin{array}{}\text{(51)}& \begin{array}{rl}\underset{n_{\mathrm{vGB}}\to +\mathrm{\infty }}{\lim }& {\left(\mathrm{1}+|h^{*}{|}^{n_{\mathrm{vGB}}}\right)}^{-m_{\mathrm{vGB}}}\approx |h^{*}{|}^{-n_{\mathrm{vGB}}{m}_{\mathrm{vGB}}}\\ & \approx |h^{*}{|}^{-{\mathit{\lambda }}_{\mathrm{vGB}}}\\ & \mathrm{with}{\mathit{\lambda }}_{\mathrm{vGB}}=n_{\mathrm{vGB}}{m}_{\mathrm{vGB}}=n_{\mathrm{vGB}}-\mathrm{2}.\end{array}\end{array}$$

The equation ${\mathit{\lambda }}_{\mathrm{BC}}={\mathit{\lambda }}_{\mathrm{vGB}}=n_{\mathrm{vGB}}-\mathrm{2}$ defines a relation between λ_BC and n_vGB that ensures a similar state for WR functions. Substituting n_vGB with m_vGB according to $m_{\mathrm{vGB}}=\mathrm{1}-\frac{\mathrm{2}}{n_{\mathrm{vGB}}}$ leads to

$$\begin{array}{}\text{(52)}& m_{\mathrm{vGB}}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{BC}}}{\mathrm{2}+{\mathit{\lambda }}_{\mathrm{BC}}}}.\end{array}$$

Since m_vGB is the appropriate shape index for the vGB model, we consider its equivalent, $\frac{{\mathit{\lambda }}_{\mathrm{BC}}}{\mathrm{2}+{\mathit{\lambda }}_{\mathrm{BC}}}$, to be the appropriate shape index for the BC model, leading to

$$\begin{array}{}\text{(53)}& x_{\mathrm{BC}}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{BC}}}{\mathrm{2}+{\mathit{\lambda }}_{\mathrm{BC}}}}.\end{array}$$

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

$$\begin{array}{}\text{(54)}& x_{\mathrm{KG}}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}+{\mathit{\sigma }}_{\mathrm{KG}}}}.\end{array}$$

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:

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

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:

$$\begin{array}{}\text{(56)}& \left\{\begin{array}{l}{c}_{\mathrm{p},\mathrm{BC}}\left(x\right)=\mathrm{2}+\frac{\mathrm{1}-x}{\mathrm{5}x+\mathrm{1}}+\frac{\mathrm{1}-x}{\mathrm{7}x+\mathrm{1}}\\ c_{\mathrm{p},\mathrm{vGB}}\left(x\right)=\mathrm{\Gamma }\left(\frac{\mathrm{3}-x}{\mathrm{2}}\right)\left[\frac{\mathrm{\Gamma }\left(\frac{\mathrm{1}+\mathrm{5}x}{\mathrm{2}}\right)}{\mathrm{\Gamma }\left(\mathrm{1}+\mathrm{2}x\right)}+\frac{\mathrm{\Gamma }\left(\frac{\mathrm{1}+\mathrm{7}x}{\mathrm{2}}\right)}{\mathrm{\Gamma }\left(\mathrm{1}+\mathrm{3}x\right)}\right]\\ c_{\mathrm{p},\mathrm{vGM}}\left(x\right)=\mathrm{\Gamma }(\mathrm{2}-x)\left[\frac{\mathrm{\Gamma }\left(\frac{\mathrm{3}}{\mathrm{2}}x\right)}{\left(\frac{\mathrm{3}}{\mathrm{2}}x-\mathrm{1}\right)\mathrm{\Gamma }\left(\frac{\mathrm{1}}{\mathrm{2}}x\right)}\right.\\ \left.+\frac{\mathrm{\Gamma }\left(\frac{\mathrm{5}}{\mathrm{2}}x\right)}{\left(\frac{\mathrm{5}}{\mathrm{2}}x-\mathrm{1}\right)\mathrm{\Gamma }\left(\frac{\mathrm{3}}{\mathrm{2}}x\right)}\right]+\left(\mathrm{1}-x\right)\left[\frac{\mathrm{\Gamma }\left(\frac{\mathrm{3}}{\mathrm{2}}x\right)\mathrm{\Gamma }\left(\mathrm{1}+x\right)}{\left(\frac{\mathrm{3}}{\mathrm{2}}x-\mathrm{1}\right)\mathrm{\Gamma }\left(\frac{\mathrm{5}}{\mathrm{2}}x\right)}\right.\\ \left.+\frac{\mathrm{\Gamma }\left(\frac{\mathrm{5}}{\mathrm{2}}x\right)\mathrm{\Gamma }\left(\mathrm{1}+x\right)}{\left(\frac{\mathrm{5}}{\mathrm{2}}x-\mathrm{1}\right)\mathrm{\Gamma }\left(\frac{\mathrm{7}}{\mathrm{2}}x\right)}-\mathrm{2}\left(\frac{\mathrm{1}}{\frac{\mathrm{3}}{\mathrm{2}}x-\mathrm{1}}+\frac{\mathrm{1}}{\frac{\mathrm{5}}{\mathrm{2}}x-\mathrm{1}}\right)\right]\\ c_{\mathrm{p},\mathrm{KG}}\left(x\right)=\underset{\mathrm{0}}{\overset{\frac{\mathrm{1}}{\mathrm{2}}}{\int }}\left(\mathrm{1}+S_{\mathrm{e}}\right)D_{\mathrm{KG}\left(x\right)}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}\\ +\underset{h_{\mathrm{KG}\left(x\right)}^{*}\left(\frac{\mathrm{1}}{\mathrm{2}}\right)}{\overset{\mathrm{0}}{\int }}\left(\mathrm{1}+S_{\mathrm{e},\mathrm{KG}\left(x\right)}\left(h^{*}\right)\right)K_{\mathrm{r},\mathrm{KG}\left(\mathrm{x}\right)}\left(h^{*}\right)\mathrm{d}{h}^{*}\\ \mathrm{with}{\mathit{\sigma }}_{\mathrm{KG}}\left(x\right)=\frac{\mathrm{1}-x}{x}.\end{array}\right.\end{array}$$

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 $\mathit{\eta }=\frac{\mathrm{2}}{\mathit{\lambda }}+\mathrm{2}+p$ and $l_{\mathrm{vGM}}=l_{\mathrm{KG}}=\frac{\mathrm{1}}{\mathrm{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_{\mathrm{p},\mathrm{d}}=\mathrm{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 sorptivity c_p.

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 λ_BC, m_vGM, m_vGB, and σ_KG were computed as a function of x using Eq. (55) with $l_{\mathrm{vGM}}=l_{\mathrm{KG}}=\frac{\mathrm{1}}{\mathrm{2}}$. The dashed line represents the delta model.

The comparison between the hydraulic models (with the same shape index value) reveals some similarities and discrepancies (Fig. 1). Three of the water retention models (vGB, vGM, KG) exhibit an inflection point with a continuous increase in S_e(h^*) over the whole interval ($-\mathrm{\infty },\mathrm{0}$] (Fig. 1a, vGB, vGM, and KG), while the BC model reaches the asymptote S_e=1 at $h^{*}=-\mathrm{1}$ with full saturation for $h^{*}\ge -\mathrm{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 inflection 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^{*}=-\mathrm{1}$ for the KG model. The other models have inflection points positioned 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 saturation compared to intermediate saturation degrees. In particular, these two models exhibit a very large increase close to saturation, whereas BC and vGB models have a gradual increase (Fig. 1b, vGM and KG versus BC and vGB, close to S_e=1). This feature allows the vGM and KG models to simulate large drops in hydraulic conductivity close to saturation that are typical of certain soils. The functions K_r(h^*) are depicted in Fig. 1c and combine the properties of the functions 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 -\mathrm{1}$ followed by a plateau (Fig. 1c, BC and vGB). The vGM and KG models 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 progressive 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 maximum equal to λ_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 infinity close to saturation, at S_e=1 (Fig. 1d, vGB versus BC). The vGB model defines an S shape that reflects the larger increases 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).

Figure 2 Impact of the shape index, x, on the WR and HC functions versus the selected hydraulic models: WR functions, S_e(h^*) (a, e, i, m); HC functions, K_r(S_e) (b, f, j, n) or K_r(h^*) (c, g, k, o); and diffusivity function D^*(S_e) 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 λ_BC, m_vGM, m_vGB, and σ_KG were computed as a function of x using Eq. (55), with $l_{\mathrm{vGM}}=l_{\mathrm{KG}}=\frac{\mathrm{1}}{\mathrm{2}}$.

Varying the shape index changes the WR and HC functions in an expected way (Fig. 2). For the WR functions, increasing the shape index from 0 to 1 makes the shift from a gradual and moderate to an abrupt increase in saturation degree, respectively. Values close to unity make the WR functions close to a stepwise function correspond to the delta model (Fig. 2, first column, arrows). As for the WR functions, 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 -\mathrm{1}$, whereas K_r(h^*) remains equal to unity above (Fig. 2c, third column). Conversely, for the vGB, vGM, and KG models, the increase in the shape index has two antagonist effects: a decrease of K_r(h^*) for $h^{*}\le -\mathrm{1}$, followed by an increase for $h^{*}\ge -\mathrm{1}$ (Fig. 2g, k, o, third column, arrows). Briefly, as expected, the water retention and the relative hydraulic conductivity tend towards stepwise functions when the shape index 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 conductivity K_r(h^*) is not greatly impacted close to $h^{*}=\mathrm{0}$ (Fig. 2c and g). Conversely, for the vGM and KG models, K_r(h^*) tends towards zero in the vicinity of $h^{*}=\mathrm{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).

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

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_{\mathrm{r},\mathrm{BC}}\left(h^{*}\right)$ functions decrease over the whole interval (−∞, 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_{\mathrm{p},\mathrm{BC}}\left(\mathrm{0}\right)=\mathrm{4}$ and $c_{\mathrm{p},\mathrm{BC}}\left(\mathrm{1}\right)=\mathrm{2}$. The air-entry water pressure head is non-null, so that the c_p and $c_{\mathrm{p}}^{\prime }$ are related by $c_{\mathrm{p}}=c_{\mathrm{p}}^{\prime }+\mathrm{2}\left|h_{\mathrm{a}}^{*}\right|=c_{\mathrm{p}}^{\prime }+\mathrm{2}$, since by convention $h_{\mathrm{a}}^{*}=-\mathrm{1}$. Consequently, $c_{\mathrm{p}}^{\prime }$ decreases from 2 to 0, with a simple vertical shift compared to c_p (Fig. 3a, $c_{\mathrm{p}}^{\prime }$). For the other hydraulic models, the air-entry water pressure head is null, so that $c_{\mathrm{p}}=c_{\mathrm{p}}^{\prime }$ (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≥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 -\mathrm{1}$ and an increase for $h^{*}\ge -\mathrm{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_{\mathrm{p},\mathrm{vGB}}\left(\mathrm{0}\right)=\mathrm{2}\mathrm{\Gamma }\left(\frac{\mathrm{3}}{\mathrm{2}}\right)\mathrm{\Gamma }\left(\frac{\mathrm{1}}{\mathrm{2}}\right)=\mathrm{\Gamma }{\left(\frac{\mathrm{1}}{\mathrm{2}}\right)}^{\mathrm{2}}=\mathit{\pi }$ and an upper limit of $c_{\mathrm{p},\mathrm{vGB}}\left(\mathrm{1}\right)=\mathrm{\Gamma }\left(\mathrm{1}\right)\left[\frac{\mathrm{\Gamma }\left(\mathrm{3}\right)}{\mathrm{\Gamma }\left(\mathrm{3}\right)}+\frac{\mathrm{\Gamma }\left(\mathrm{4}\right)}{\mathrm{\Gamma }\left(\mathrm{4}\right)}\right]=\mathrm{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_{\mathrm{p},\mathrm{vGM}}\left(\mathrm{0}\right)=\mathrm{0}$ (numerical determination) and $c_{\mathrm{p},\mathrm{vGM}}\left(\mathrm{1}\right)=\mathrm{\Gamma }\left(\mathrm{1}\right)\left[\frac{\mathrm{\Gamma }\left(\frac{\mathrm{1}}{\mathrm{2}}\right)}{\mathrm{\Gamma }\left(\frac{\mathrm{1}}{\mathrm{2}}\right)}+\frac{\mathrm{\Gamma }\left(\frac{\mathrm{3}}{\mathrm{2}}\right)}{\mathrm{\Gamma }\left(\frac{\mathrm{3}}{\mathrm{2}}\right)}\right]=\mathrm{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 [\mathrm{0},\mathrm{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_{\mathrm{p},\mathrm{d}}=\mathrm{2}$ (see Eq. 42). Such a result indicates that similar results should be obtained for soils with a narrow pore size distribution (like coarse soils with narrow pore size distributions), regardless 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 distributions. 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_p for the KG and vGB models are in line with the large decrease in D^*(S_e) over the whole interval [0, 1] when the shape index, x, tends towards 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_p,BC(x), c_p,vGB(x), c_p,vGM(x), and c_p,KG(x) exhibit contrasting evolution, with c_p,vGM(x) and c_p,KG(x) increasing monotonously, c_p,BC(x) decreasing monotonously, and c_p,vGB(x) exhibiting a two-step behavior of decreasing followed by increasing values.

This contrasting evolution of functions c_p,BC(x), c_p,vGB(x), c_p,vGM(x), and c_p,KG(x) points at different implications 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_{\mathrm{K}}^{\mathrm{2}}(h_{\mathrm{0}},\mathrm{0})$, for a given soil and given initial conditions, h₀. 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 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, $({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}})K_{\mathrm{s}}\left|h_{\mathrm{g}}\right|$, to compensate for the very low values of c_p (see Eq. 26). Given that the scale parameters θ_s, θ_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, $\left|h_{\mathrm{g}}\right|$, 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 $\left|h_{\mathrm{g}}\right|$ must tends towards being infinite when the shape index tends towards zero for these two models. For the KG model, such a relation, $\left|h_{\mathrm{KG}}\right|\sim x_{\mathrm{KG}}$, is related to the relation $\left|h_{\mathrm{KG}}\right|\sim {\mathit{\sigma }}_{\mathrm{KG}}$ since $x_{\mathrm{KG}}=(\mathrm{1}+{\mathit{\sigma }}_{\mathrm{KG}}{)}^{-\mathrm{1}}$. Our statements imply that very large values of σ_KG(x_KG→0) should be associated with very large values of $\left|h_{\mathrm{KG}}\right|$, 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.

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

3.3 Upscaling sorptivity S_K(h₀,0) from c_p

In this section, we elaborate on the use of Eq. (26) for the easy and straightforward computation of $S_{\mathrm{K}}^{\mathrm{2}}(h_{\mathrm{0}},\mathrm{0})$ from the tabulated values of c_p (Table 1). The proposed scaling procedure Eq. (26) allows for the computation of $S_{\mathrm{K}}^{\mathrm{2}}(h_{\mathrm{0}},\mathrm{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 (λ_BC, m_vGB, m_vGM, or σ_KG) to compute the related shape index, x, considering the following definitions: x_vGB=m_vGB, x_vGM=m_vGM, $x_{\mathrm{KG}}=\frac{\mathrm{1}}{\mathrm{1}+{\mathit{\sigma }}_{\mathrm{KG}}}$, and $x_{\mathrm{BC}}=\frac{{\mathit{\lambda }}_{\mathrm{BC}}}{\mathrm{2}+{\mathit{\lambda }}_{\mathrm{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, θ₀ or θ(h₀), depending on the description of the initial condition (either θ₀ or h₀).

4. Compute the related hydraulic conductivity, K₀=K(θ₀), using the HC function.

5. Compute the correcting factors $R_{\mathit{\theta }}=\frac{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}$ ($=\mathrm{1}-S_{\mathrm{e},\mathrm{0}}$) and $R_{\mathrm{K}}=\frac{K_{\mathrm{s}}-K_{\mathrm{0}}}{K_{\mathrm{s}}}$ ($=\mathrm{1}-K_{\mathrm{r}}\left(S_{\mathrm{e},\mathrm{0}}\right))$.

6. Compute the scaled air-entry water pressure head, $\left|h_{\mathrm{a}}^{*}\right|=\left|\frac{h_{\mathrm{a}}}{h_{\mathrm{g}}}\right|$.

7. Compute the square scaled sorptivity, $S_{\mathrm{K}}^{*\mathrm{2}}(h_{\mathrm{0}}^{*},\mathrm{0})=R_{\mathrm{K}}{R}_{\mathit{\theta }}(c_{\mathrm{p}}-\mathrm{2}|h_{\mathrm{a}}^{*}\left|\right)+\mathrm{2}{R}_{\mathit{\theta }}\left|h_{\mathrm{a}}^{*}\right|$.

8. Upscale to derive the square sorptivity, $S_{\mathrm{K}}^{\mathrm{2}}(h_{\mathrm{0}},\mathrm{0})=S_{\mathrm{K}}^{*\mathrm{2}}(h_{\mathrm{0}}^{*},\mathrm{0})({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}})K_{\mathrm{s}}\left|h_{\mathrm{g}}\right|$.

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₁=0) and an initial water pressure head of $h_{\mathrm{0}}=-\mathrm{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: θ_r=0.078, θ_s=0.43, $h_{\mathrm{g}}=-\mathrm{277}$ mm, $K_{\mathrm{s}}=\mathrm{2.88}\times {\mathrm{10}}^{-\mathrm{3}}$ mm s⁻¹, 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_{\mathrm{vGM}}=\mathrm{1}-\mathrm{1}/n_{\mathrm{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): θ₀=0.125;

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

5. corresponding correction factors: $R_{\mathit{\theta }}=\frac{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}$ and $R_K=\frac{K_{\mathrm{s}}-K_{\mathrm{0}}}{K_{\mathrm{0}}}$, leading to R_θ=0.865 and R_K=1.000;

6. air-entry water pressure head: no air-entry water pressure head, leading to $\left|h_{\mathrm{a}}^{*}\right|=\mathrm{0}$;

7. square scaled sorptivity: $S_{\mathrm{K}}^{*\mathrm{2}}(h_{\mathrm{0}}^{*},\mathrm{0})=R_{\mathrm{K}}{R}_{\mathit{\theta }}{c}_{\mathrm{p}}$, leading to $S_{\mathrm{K}}^{*\mathrm{2}}(h_{\mathrm{0}}^{*},\mathrm{0})=\mathrm{0.416}$;

8. sorptivity: $S_{\mathrm{K}}^{\mathrm{2}}(h_{\mathrm{0}},\mathrm{0})=S_{\mathrm{K}}^{*\mathrm{2}}(h_{\mathrm{0}}^{*},\mathrm{0})({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}})K_{\mathrm{s}}\left|h_{\mathrm{g}}\right|$, leading to $S_{\mathrm{K}}^{\mathrm{2}}(h_{\mathrm{0}},\mathrm{0})=\mathrm{0.117}$ mm² s⁻¹ and $S_{\mathrm{K}}(h_{\mathrm{0}},\mathrm{0})=\mathrm{0.342}$ mm s^−½.

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: θ_r=0.078, θ_s=0.43, $h_{\mathrm{g}}=-\mathrm{277}$ mm, and $K_{\mathrm{s}}=\mathrm{2.88}\times {\mathrm{10}}^{-\mathrm{3}}$ mm s⁻¹, with a value of λ_BC=0.56. λ_BC was deduced from the previous value of n=1.56 considering the usual relation λ=mn, as suggested by Haverkamp et al. (2005). The application of the proposed procedure leads to the following computations:

1. shape index: $x_{\mathrm{BC}}={\mathit{\lambda }}_{\mathrm{BC}}/(\mathrm{2}+{\mathit{\lambda }}_{\mathrm{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): θ₀=0.125;

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

5. corresponding correction factors: $R_{\mathit{\theta }}=\frac{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{0}}}{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}$ and $R_{\mathrm{K}}=\frac{K_{\mathrm{s}}-K_{\mathrm{0}}}{K_{\mathrm{0}}}$, leading to R_θ=0.866 and R_K=1.000;

6. air-entry water pressure head: significant air-entry water pressure head, with h_BC=h_a, leading to $\left|h_{\mathrm{a}}^{*}\right|=\mathrm{1}$;

7. square scaled sorptivity: $S_{\mathrm{K}}^{*\mathrm{2}}(h_{\mathrm{0}}^{*},\mathrm{0})=R_{\mathrm{K}}{R}_{\mathit{\theta }}(c_{\mathrm{p}}-\mathrm{2}|h_{\mathrm{a}}^{*}\left|\right)+\mathrm{2}{R}_{\mathit{\theta }}\left|h_{\mathrm{a}}^{*}\right|$, leading to $S_{\mathrm{K}}^{*\mathrm{2}}(h_{\mathrm{0}}^{*},\mathrm{0})=\mathrm{2.318}$;

8. sorptivity: $S_{\mathrm{K}}^{\mathrm{2}}(h_{\mathrm{0}},\mathrm{0})=S_{\mathrm{K}}^{*\mathrm{2}}(h_{\mathrm{0}}^{*},\mathrm{0})({\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}})K_{\mathrm{s}}\left|h_{\mathrm{g}}\right|$, leading to $S_{\mathrm{K}}^{\mathrm{2}}(h_{\mathrm{0}},\mathrm{0})=\mathrm{0.651}$ mm² s⁻¹, and $S_{\mathrm{K}}(h_{\mathrm{0}},\mathrm{0})=\mathrm{0.806}$ mm s^−½.

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_{\mathrm{K}}^{\mathrm{2}}(h_{\mathrm{0}},\mathrm{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 ${\mathit{\theta }}_{\mathrm{0}}\le \frac{\mathrm{1}}{\mathrm{4}}{\mathit{\theta }}_{\mathrm{s}}$. For fine soils, even a small initial water pressure head may cause ${\mathit{\theta }}_{\mathrm{0}}>\frac{\mathrm{1}}{\mathrm{4}}{\mathit{\theta }}_{\mathrm{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: θ_r=0.07, θ_s=0.36, $h_{\mathrm{g}}=-\mathrm{2000}$ mm, $K_{\mathrm{s}}=\mathrm{5.555}\times {\mathrm{10}}^{-\mathrm{5}}$ mm s⁻¹, n=1.09, and l=0.5. Considering the same value for the initial water pressure head, i.e., $h_{\mathrm{0}}=-\mathrm{10}$ m, the initial water content is θ₀=0.318, which exceeds $\frac{\mathrm{1}}{\mathrm{4}}{\mathit{\theta }}_{\mathrm{s}}$. The application of the scaling procedure led to an estimated sorptivity of 0.0475 mm s^−½, whereas the targeted sorptivity computed with Eq. (2) was 0.0127 mm s^−½. Such an error corresponds to an overestimation by a factor of 2.73. Thus, we advise that the user verify that ${\mathit{\theta }}_{\mathrm{0}}\le \frac{\mathrm{1}}{\mathrm{4}}{\mathit{\theta }}_{\mathrm{s}}$ when θ_r=0, or more generally $S_{e,\mathrm{0}}\le \frac{\mathrm{1}}{\mathrm{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_θ), 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 ${\mathit{\theta }}_{\mathrm{0}}\le \frac{\mathrm{1}}{\mathrm{4}}{\mathit{\theta }}_{\mathrm{s}}$ when θ_r=0, or more generally $S_{e,\mathrm{0}}\le \frac{\mathrm{1}}{\mathrm{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, σ_KG (Pollacco et al., 2013). In other words, the scale parameter, h_KG, should be parametrized as a function of the shape parameter, σ_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 θ_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 λ, m, n, σ, in order to avoid heavy equations. The dimensionless diffusivity functions were derived from their definition D^*(S_e), applying $D^{*}\left(S_{\mathrm{e}}\right)=K_{\mathrm{r}}\left(S_{\mathrm{e}}\right)\frac{\mathrm{d}{h}^{*}}{\mathrm{d}{S}_{\mathrm{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:

$$\begin{array}{}\text{(A1)}& {\displaystyle }& {\displaystyle }{h}_{\mathrm{BC}}^{*}\left(S_{\mathrm{e}}\right)=-S_{\mathrm{e}}^{-\frac{\mathrm{1}}{\mathit{\lambda }}}\text{(A2)}& {\displaystyle }& {\displaystyle }{h}_{\mathrm{vGB}}^{*}\left(S_{\mathrm{e}}\right)=-{\left(S_{\mathrm{e}}^{-\frac{\mathrm{1}}{m}}-\mathrm{1}\right)}^{\frac{\mathrm{1}}{n}}\mathrm{with}m=\mathrm{1}-{\displaystyle \frac{\mathrm{2}}{n}}\text{(A3)}& {\displaystyle }& {\displaystyle }{h}_{\mathrm{vGM}}^{*}\left(S_{\mathrm{e}}\right)=-{\left(S_{\mathrm{e}}^{-\frac{\mathrm{1}}{m}}-\mathrm{1}\right)}^{\frac{\mathrm{1}}{n}}\mathrm{with}m=\mathrm{1}-{\displaystyle \frac{\mathrm{1}}{n}}\text{(A4)}& {\displaystyle }& {\displaystyle }{h}_{\mathrm{KG}}^{*}\left(S_{\mathrm{e}}\right)=-e^{\sqrt{\mathrm{2}}\mathit{\sigma }{\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}{S}_{\mathrm{e}}\right)},\end{array}$$

where erfc⁻¹ is the inverse function of the complementary error function. These functions can be differentiated to define their relative derivatives, $\frac{\mathrm{d}{h}^{*}}{\mathrm{d}{S}_{\mathrm{e}}}$:

$$\begin{array}{}\text{(A5)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{h}_{\mathrm{BC}}^{*}}{\mathrm{d}{S}_{\mathrm{e}}}}\left(S_{\mathrm{e}}\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }}}{S}_{\mathrm{e}}^{-\frac{\mathrm{1}}{\mathit{\lambda }}-\mathrm{1}}\text{(A6)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{h}_{\mathrm{vGB}}^{*}}{\mathrm{d}{S}_{\mathrm{e}}}}\left(S_{\mathrm{e}}\right)={\displaystyle \frac{\mathrm{1}-m}{\mathrm{2}m}}{S}_{\mathrm{e}}^{-\frac{\mathrm{1}+m}{\mathrm{2}m}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{-\frac{\mathrm{1}+m}{\mathrm{2}}}\text{(A7)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{h}_{\mathrm{vGM}}^{*}}{\mathrm{d}{S}_{\mathrm{e}}}}\left(S_{\mathrm{e}}\right)={\displaystyle \frac{\mathrm{1}-m}{m}}{S}_{\mathrm{e}}^{-\frac{\mathrm{1}}{m}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{-m}\text{(A8)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{h}_{\mathrm{KG}}^{*}}{\mathrm{d}{S}_{\mathrm{e}}}}\left(S_{\mathrm{e}}\right)=\sqrt{\mathrm{2}\mathit{\pi }}\mathit{\sigma }{e}^{{\left({\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}{S}_{\mathrm{e}}\right)\right)}^{\mathrm{2}}+\sqrt{\mathrm{2}}\mathit{\sigma }{\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}{S}_{\mathrm{e}}\right)}.\end{array}$$

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

The scaled hydraulic conductivity functions are

$$\begin{array}{}\text{(A9)}& {\displaystyle }& {\displaystyle }{K}_{\mathrm{r},\mathrm{BC}}\left(S_{\mathrm{e}}\right)=S_{\mathrm{e}}^{\mathit{\eta }},\text{(A10)}& {\displaystyle }& {\displaystyle }{K}_{\mathrm{r},\mathrm{vGB}}\left(S_{\mathrm{e}}\right)=S_{\mathrm{e}}^{\mathit{\eta }},\text{(A11)}& {\displaystyle }& {\displaystyle }{K}_{\mathrm{r},\mathrm{vGM}}\left(S_{\mathrm{e}}\right)=S_{\mathrm{e}}^l\left(\mathrm{1}-\mathrm{2}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^m+{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{\mathrm{2}m}\right),\text{(A12)}& {\displaystyle }& {\displaystyle }{K}_{\mathrm{r},\mathrm{KG}}\left(S_{\mathrm{e}}\right)=S_{\mathrm{e}}^l{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\mathrm{erfc}\left({\mathrm{erfc}}^{-\mathrm{1}}\left(\mathrm{2}{S}_{\mathrm{e}}\right)+{\displaystyle \frac{\mathit{\sigma }}{\sqrt{\mathrm{2}}}}\right)\right)}^{\mathrm{2}}.\end{array}$$

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

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

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_{\mathrm{a}}^{*}=-\mathrm{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 $\left|h_{\mathrm{a}}^{*}\right|=\mathrm{1}$. Then, the first part $\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}(\mathrm{1}+S_{\mathrm{e}})D_{\mathrm{BC}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}$ can be integrated analytically given that the hydraulic diffusivity $D_{\mathrm{BC}}^{*}\left(S_{\mathrm{e}}\right)$ obeys a power law. The following developments can be done:

$$\begin{array}{}\text{(B1)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{BC}}& =\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(\mathrm{1}+S_{\mathrm{e}}\right)D_{\mathrm{BC}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}+\mathrm{2}\left|h_{\mathrm{a}}^{*}\right|\\ & =\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(\mathrm{1}+S_{\mathrm{e}}\right)D_{\mathrm{BC}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}+\mathrm{2}\end{array}\end{array}$$

$$\begin{array}{}\text{(B2)}& {\displaystyle }\begin{array}{rl}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(\mathrm{1}+S_{\mathrm{e}}\right)& D_{\mathrm{BC}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{D}_{\mathrm{BC}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}\\ & +\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}{D}_{\mathrm{BC}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}\\ & =\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }}}{S}_{\mathrm{e}}^{\mathit{\eta }-\left(\frac{\mathrm{1}}{\mathit{\lambda }}+\mathrm{1}\right)}\mathrm{d}{S}_{\mathrm{e}}+\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }}}{S}_{\mathrm{e}}{S}_{\mathrm{e}}^{\mathit{\eta }-\left(\frac{\mathrm{1}}{\mathit{\lambda }}+\mathrm{1}\right)}\mathrm{d}{S}_{\mathrm{e}}\\ & ={\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }}}\left(\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}^{\mathit{\eta }-\left(\frac{\mathrm{1}}{\mathit{\lambda }}+\mathrm{1}\right)}\mathrm{d}{S}_{\mathrm{e}}+\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}^{\mathit{\eta }-\frac{\mathrm{1}}{\mathit{\lambda }}}\mathrm{d}{S}_{\mathrm{e}}\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathit{\lambda }}}\left({\displaystyle \frac{\mathrm{1}}{\mathit{\eta }-\frac{\mathrm{1}}{\mathit{\lambda }}}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\eta }-\frac{\mathrm{1}}{\mathit{\lambda }}+\mathrm{1}}}\right)\\ & ={\displaystyle \frac{\mathrm{1}}{\mathit{\eta }\mathit{\lambda }-\mathrm{1}}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\eta }\mathit{\lambda }+\mathit{\lambda }-\mathrm{1}}}.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B3)}& c_{\mathrm{p},\mathrm{BC}}(\mathit{\lambda },\mathit{\eta })={\displaystyle \frac{\mathrm{1}}{\mathit{\eta }\mathit{\lambda }-\mathrm{1}}}+{\displaystyle \frac{\mathrm{1}}{\mathit{\eta }\mathit{\lambda }+\mathit{\lambda }-\mathrm{1}}}+\mathrm{2}.\end{array}$$

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

$$\begin{array}{}\text{(B4)}& c_{\mathrm{p},\mathrm{BC}}\left(\mathit{\lambda }\right)=\mathrm{2}+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}\mathit{\lambda }+\mathrm{1}}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\lambda }+\mathrm{1}}}.\end{array}$$

These development demonstrate the equations proposed for c_p for the BC model, i.e., Eqs. (43)–(44). This demonstration is in line with Varado et al. (2006).

B2 Parameter c_p for the vGB model

For the vGB model, and the remaining models, there is no air-entry water pressure head, $h_{\mathrm{a}}^{*}=\mathrm{0}$, leading to

$$\begin{array}{}\text{(B5)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{vGB}}& =\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(\mathrm{1}+S_{\mathrm{e}}\right)D_{\mathrm{vGB}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}+\mathrm{2}\left|h_{\mathrm{a}}^{*}\right|\\ & =\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(\mathrm{1}+S_{\mathrm{e}}\right)D_{\mathrm{vGB}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}.\end{array}\end{array}$$

Then, the integral $\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}(\mathrm{1}+S_{\mathrm{e}})D_{\mathrm{vGB}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}$ can be decomposed into well-known integrals:

$$\begin{array}{}\text{(B6)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{vGB}}& (m,\mathit{\eta })=\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(\mathrm{1}+S_{\mathrm{e}}\right)D_{\mathrm{vGB}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}\\ & =\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{D}_{\mathrm{vGB}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}+\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}{D}_{\mathrm{vGB}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}\\ & =\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\displaystyle \frac{\mathrm{1}-m}{\mathrm{2}m}}{S}_{\mathrm{e}}^{\mathit{\eta }-\frac{\mathrm{1}+m}{\mathrm{2}m}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{-\frac{\mathrm{1}+m}{\mathrm{2}}}\mathrm{d}{S}_{\mathrm{e}}\\ & +\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}{\displaystyle \frac{\mathrm{1}-m}{\mathrm{2}m}}{S}_{\mathrm{e}}^{\mathit{\eta }-\frac{\mathrm{1}+m}{\mathrm{2}m}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{-\frac{\mathrm{1}+m}{\mathrm{2}}}\mathrm{d}{S}_{\mathrm{e}}\\ & ={\displaystyle \frac{\mathrm{1}-m}{\mathrm{2}m}}\left(\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}^{\mathit{\eta }-\frac{\mathrm{1}+m}{\mathrm{2}m}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{-\frac{\mathrm{1}+m}{\mathrm{2}}}\mathrm{d}{S}_{\mathrm{e}}\right.\\ & \left.+\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}{S}_{\mathrm{e}}^{\mathit{\eta }-\frac{\mathrm{1}+m}{\mathrm{2}m}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{-\frac{\mathrm{1}+m}{\mathrm{2}}}\mathrm{d}{S}_{\mathrm{e}}\right)\\ & ={\displaystyle \frac{\mathrm{1}-m}{\mathrm{2}m}}\left(\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}^{\mathit{\eta }-\frac{\mathrm{1}+m}{\mathrm{2}m}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{-\frac{\mathrm{1}+m}{\mathrm{2}}}\mathrm{d}{S}_{\mathrm{e}}\right.\\ & \left.+\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}^{\mathit{\eta }+\mathrm{1}-\frac{\mathrm{1}+m}{\mathrm{2}m}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{-\frac{\mathrm{1}+m}{\mathrm{2}}}\mathrm{d}{S}_{\mathrm{e}}\right)\\ & ={\displaystyle \frac{\mathrm{1}-m}{\mathrm{2}m}}\left(\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\left(S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{m\mathit{\eta }-\frac{\mathrm{1}+m}{\mathrm{2}}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{-\frac{\mathrm{1}+m}{\mathrm{2}}}\mathrm{d}{S}_{\mathrm{e}}\right.\\ & \left.+\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\left(S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{m\mathit{\eta }+m-\frac{\mathrm{1}+m}{\mathrm{2}}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{-\frac{\mathrm{1}+m}{\mathrm{2}}}\mathrm{d}{S}_{\mathrm{e}}\right).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B7)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{vGB}}(m,\mathit{\eta })& ={\displaystyle \frac{\mathrm{1}-m}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{y}^{m\mathit{\eta }+\frac{m}{\mathrm{2}}-\frac{\mathrm{3}}{\mathrm{2}}}(\mathrm{1}-y{)}^{-\frac{\mathrm{1}+m}{\mathrm{2}}}\mathrm{d}y\\ & +{\displaystyle \frac{\mathrm{1}-m}{\mathrm{2}}}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{y}^{m\mathit{\eta }+\frac{\mathrm{3}m}{\mathrm{2}}-\frac{\mathrm{3}}{\mathrm{2}}}(\mathrm{1}-y{)}^{-\frac{\mathrm{1}+m}{\mathrm{2}}}\mathrm{d}y.\end{array}\end{array}$$

We now use the beta function B, which has the following properties (assuming x>0 and y>0):

$$\begin{array}{}\text{(B8)}& \begin{array}{rl}B(x,y)& =\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{t}^{x-\mathrm{1}}(\mathrm{1}-t{)}^{y-\mathrm{1}}\mathrm{d}t\\ & ={\displaystyle \frac{\mathrm{\Gamma }\left(x\right)\mathrm{\Gamma }\left(y\right)}{\mathrm{\Gamma }(x+y)}},\end{array}\end{array}$$

where the Γ function is already defined by Eq. (46): $\mathrm{\Gamma }\left(z\right)=\underset{\mathrm{0}}{\overset{+\mathrm{\infty }}{\int }}{t}^{z-\mathrm{1}}{e}^{-t}\mathrm{d}t$ (z>0). Then, the parameter c_p,vGB can be expressed as follows:

$$\begin{array}{}\text{(B9)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{vGB}}(m,\mathit{\eta })& ={\displaystyle \frac{\mathrm{1}-m}{\mathrm{2}}}B\left(m\mathit{\eta }+{\displaystyle \frac{m}{\mathrm{2}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},{\displaystyle \frac{\mathrm{1}-m}{\mathrm{2}}}\right)\\ & +{\displaystyle \frac{\mathrm{1}-m}{\mathrm{2}}}B\left(m\mathit{\eta }+{\displaystyle \frac{\mathrm{3}m}{\mathrm{2}}}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}},{\displaystyle \frac{\mathrm{1}-m}{\mathrm{2}}}\right)\\ & ={\displaystyle \frac{\mathrm{1}-m}{\mathrm{2}}}\mathrm{\Gamma }\left({\displaystyle \frac{\mathrm{1}-m}{\mathrm{2}}}\right)\left({\displaystyle \frac{\mathrm{\Gamma }\left(m\mathit{\eta }+\frac{m}{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}}\right)}{\mathrm{\Gamma }\left(m\mathit{\eta }\right)}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }\left(m\mathit{\eta }+\frac{\mathrm{3}m}{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}}\right)}{\mathrm{\Gamma }(m\mathit{\eta }+m)}}\right)=\mathrm{\Gamma }\left({\displaystyle \frac{\mathrm{3}-m}{\mathrm{2}}}\right)\\ & \left({\displaystyle \frac{\mathrm{\Gamma }\left(m\mathit{\eta }+\frac{m}{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}}\right)}{\mathrm{\Gamma }\left(m\mathit{\eta }\right)}}+{\displaystyle \frac{\mathrm{\Gamma }\left(m\mathit{\eta }+\frac{\mathrm{3}m}{\mathrm{2}}-\frac{\mathrm{1}}{\mathrm{2}}\right)}{\mathrm{\Gamma }(m\mathit{\eta }+m)}}\right).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B10)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{vGB}}(n,m,\mathit{\eta })& =\mathrm{\Gamma }\left(\mathrm{1}+{\displaystyle \frac{\mathrm{1}}{n}}\right)\left({\displaystyle \frac{\mathrm{\Gamma }\left(m\mathit{\eta }-\frac{\mathrm{1}}{n}\right)}{\mathrm{\Gamma }\left(m\mathit{\eta }\right)}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }\left(m\mathit{\eta }+m-\frac{\mathrm{1}}{n}\right)}{\mathrm{\Gamma }(m\mathit{\eta }+m)}}\right).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B11)}& c_{\mathrm{p},\mathrm{vGB}}\left(m\right)=\mathrm{\Gamma }\left({\displaystyle \frac{\mathrm{3}-m}{\mathrm{2}}}\right)\left({\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\mathrm{1}+\mathrm{5}m}{\mathrm{2}}\right)}{\mathrm{\Gamma }(\mathrm{1}+\mathrm{2}m)}}+{\displaystyle \frac{\mathrm{\Gamma }\left(\frac{\mathrm{1}+\mathrm{7}m}{\mathrm{2}}\right)}{\mathrm{\Gamma }(\mathrm{1}+\mathrm{3}m)}}\right).\end{array}$$

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{array}{}\text{(B12)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{vGM}}(m,l)& =\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}\left(\mathrm{1}+S_{\mathrm{e}}\right)D_{\mathrm{vGM}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}\\ & =\underset{A}{\underbrace{\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{D}_{\mathrm{vGM}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}}}+\underset{B}{\underbrace{\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}{D}_{\mathrm{vGM}}^{*}\left(S_{\mathrm{e}}\right)\mathrm{d}{S}_{\mathrm{e}}}}.\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B13)}& \begin{array}{rl}A& ={\displaystyle \frac{\mathrm{1}-m}{m}}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}^{l-\frac{\mathrm{1}}{m}}\left({\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{-m}+{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^m-\mathrm{2}\right)\mathrm{d}{S}_{\mathrm{e}}\\ & ={\displaystyle \frac{\mathrm{1}-m}{m}}\left(\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}^{l-\frac{\mathrm{1}}{m}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{-m}\mathrm{d}{S}_{\mathrm{e}}\right.\\ & \left.+\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}^{l-\frac{\mathrm{1}}{m}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^m\mathrm{d}{S}_{\mathrm{e}}-\mathrm{2}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}^{l-\frac{\mathrm{1}}{m}}\mathrm{d}{S}_{\mathrm{e}}\right)\\ & ={\displaystyle \frac{\mathrm{1}-m}{m}}\left(\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\left(S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{ml-\mathrm{1}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{-m}\mathrm{d}{S}_{\mathrm{e}}\right.\\ & \left.+\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{\left(S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^{ml-\mathrm{1}}{\left(\mathrm{1}-S_{\mathrm{e}}^{\frac{\mathrm{1}}{m}}\right)}^m\mathrm{d}{S}_{\mathrm{e}}-\mathrm{2}\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}^{l-\frac{\mathrm{1}}{m}}\mathrm{d}{S}_{\mathrm{e}}\right).\end{array}\end{array}$$

The last term of A can be simplified easily:

$$\begin{array}{}\text{(B14)}& \underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{S}_{\mathrm{e}}^{l-\frac{\mathrm{1}}{m}}\mathrm{d}{S}_{\mathrm{e}}={\left[{\displaystyle \frac{S_{\mathrm{e}}^{l-\frac{\mathrm{1}}{m}+\mathrm{1}}}{l-\frac{\mathrm{1}}{m}+\mathrm{1}}}\right]}_{\mathrm{0}}^{\mathrm{1}}={\displaystyle \frac{m}{(l+\mathrm{1})m-\mathrm{1}}}.\end{array}$$

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

$$\begin{array}{}\text{(B15)}& \begin{array}{rl}A& =(\mathrm{1}-m)\left(\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{y}^{ml+m-\mathrm{2}}(\mathrm{1}-y{)}^{-m}\mathrm{d}y\right.\\ & \left.+\underset{\mathrm{0}}{\overset{\mathrm{1}}{\int }}{y}^{ml+m-\mathrm{2}}(\mathrm{1}-y{)}^m\mathrm{d}y-{\displaystyle \frac{\mathrm{2}}{(l+\mathrm{1})m-\mathrm{1}}}\right).\end{array}\end{array}$$

In this case, we need to assume that $ml+m-\mathrm{2}>-\mathrm{1}$, i.e. $m>\frac{\mathrm{1}}{l+\mathrm{1}}$, to use the beta function. Such a condition corresponds to $m>\frac{\mathrm{2}}{\mathrm{3}}$, for the by-default value of $l=\frac{\mathrm{1}}{\mathrm{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_p:

$$\begin{array}{}\text{(B16)}& \begin{array}{rl}A& =(\mathrm{1}-m)\left(B\left(m\right(l+\mathrm{1})-\mathrm{1},\mathrm{1}-m)+B\left(m\right(l+\mathrm{1})-\mathrm{1},\mathrm{1}+m)\right.\\ & \left.-{\displaystyle \frac{\mathrm{2}}{(l+\mathrm{1})m-\mathrm{1}}}\right)=(\mathrm{1}-m)\left({\displaystyle \frac{\mathrm{\Gamma }\left(m\right(l+\mathrm{1})-\mathrm{1})\mathrm{\Gamma }(\mathrm{1}-m)}{\mathrm{\Gamma }\left(ml\right)}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }\left(m\right(l+\mathrm{1})-\mathrm{1})\mathrm{\Gamma }(\mathrm{1}+m)}{\mathrm{\Gamma }\left(m\right(l+\mathrm{2}\left)\right)}}-{\displaystyle \frac{\mathrm{2}}{(l+\mathrm{1})m-\mathrm{1}}}\right).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B17)}& \begin{array}{rl}B& =(\mathrm{1}-m)\left({\displaystyle \frac{\mathrm{\Gamma }\left(m\right(l+\mathrm{2})-\mathrm{1})\mathrm{\Gamma }(\mathrm{1}-m)}{\mathrm{\Gamma }(ml+m)}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }\left(m\right(l+\mathrm{2})-\mathrm{1})\mathrm{\Gamma }(\mathrm{1}+m)}{\mathrm{\Gamma }\left(m\right(l+\mathrm{3}\left)\right)}}-{\displaystyle \frac{\mathrm{2}}{(l+\mathrm{2})m-\mathrm{1}}}\right).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(B18)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{vGM}}(m,l)& =\mathrm{\Gamma }(\mathrm{2}-m)\left({\displaystyle \frac{\mathrm{\Gamma }\left(m\right(l+\mathrm{1})-\mathrm{1})}{\mathrm{\Gamma }\left(ml\right)}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }\left(m\right(l+\mathrm{2})-\mathrm{1})}{\mathrm{\Gamma }(ml+m)}}\right)\\ & +(\mathrm{1}-m)\left[\left({\displaystyle \frac{\mathrm{\Gamma }\left(m\right(l+\mathrm{1})-\mathrm{1})\mathrm{\Gamma }(\mathrm{1}+m)}{\mathrm{\Gamma }\left(m\right(l+\mathrm{2}\left)\right)}}\right.\right.\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }\left(m\right(l+\mathrm{2})-\mathrm{1})\mathrm{\Gamma }(\mathrm{1}+m)}{\mathrm{\Gamma }\left(m\right(l+\mathrm{3}\left)\right)}}\right)\\ & \left.-\mathrm{2}\left({\displaystyle \frac{\mathrm{1}}{(l+\mathrm{1})m-\mathrm{1}}}+{\displaystyle \frac{\mathrm{1}}{(l+\mathrm{2})m-\mathrm{1}}}\right)\right]\end{array}\end{array}$$

$$\begin{array}{}\text{(B19)}& \begin{array}{rl}{c}_{\mathrm{p},\mathrm{vGM}}(m,l)& =\mathrm{\Gamma }(\mathrm{2}-m)\left({\displaystyle \frac{\mathrm{\Gamma }\left(m\right(l+\mathrm{1}\left)\right)}{\left(m\right(l+\mathrm{1})-\mathrm{1})\mathrm{\Gamma }\left(ml\right)}}\right.\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }\left(m\right(l+\mathrm{2}\left)\right)}{\left(m\right(l+\mathrm{2})-\mathrm{1})\mathrm{\Gamma }\left(m\right(\mathrm{1}+l\left)\right)}}\right)\\ & +(\mathrm{1}-m)\left[\left({\displaystyle \frac{\mathrm{\Gamma }\left(m\right(l+\mathrm{1}\left)\right)\mathrm{\Gamma }(\mathrm{1}+m)}{\left(m\right(l+\mathrm{1})-\mathrm{1})\mathrm{\Gamma }\left(m\right(l+\mathrm{2}\left)\right)}}\right.\right.\\ & \left.+{\displaystyle \frac{\mathrm{\Gamma }\left(m\right(l+\mathrm{2}\left)\right)\mathrm{\Gamma }(\mathrm{1}+m)}{\left(m\right(l+\mathrm{2})-\mathrm{1})\mathrm{\Gamma }\left(m\right(l+\mathrm{3}\left)\right)}}\right)\\ & \left.-\mathrm{2}\left({\displaystyle \frac{\mathrm{1}}{(l+\mathrm{1})m-\mathrm{1}}}+{\displaystyle \frac{\mathrm{1}}{(l+\mathrm{2})m-\mathrm{1}}}\right)\right].\end{array}\end{array}$$

Note that, as stated above, the equation theoretically should apply only for the case of $m>\frac{\mathrm{1}}{\mathrm{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{\mathrm{1}}{\mathrm{1}+l}$ and can be considered for any value of $m\in [\mathrm{0},\mathrm{1}]$ provided that $m\ne \frac{\mathrm{1}}{\mathrm{1}+l}$ or $m\ne \frac{\mathrm{1}}{\mathrm{2}+l}$. Equation (B19) demonstrates Eq. (48). The combination of these equations with the capillary model Eq. (37) leads to Eq. (49).