A scaling procedure for straightforward computation of sorptivity

. Sorptivity is a parameter of primary importance in the study of unsaturated ﬂow in soils. This hydraulic parameter is required to model water inﬁltration into vertical soil proﬁles. 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 ﬁnal saturation degrees, which may be difﬁcult in some situations (e.g., coarse soil with diffusivity functions that are quasi-inﬁnite 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 ﬁnal state. The c p parameter was computed numerically and analytically for ﬁve 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 sorp-tivity.


5084
L. Lassabatere et al.: Scaling sorptivity for easy computation et al., 2006). Sorptivity is incorporated in a wide range of infiltration models (Angulo-Jaramillo et al., 2016;Lassabatere et al., 2009Lassabatere et al., , 2014Lassabatere et al., , 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 In the above equations, S(θ 0 , θ 1 ) stands for the sorptivity between θ 0 and θ 1 , χ (θ ) = x(θ,t) √ t stands for the Boltzmann transformation variable, θ 0 is the initial water content, θ 1 is the final water content corresponding also to the water content applied at the soil surface, and t is the elapsed time. In practical applications, the user must perform numerical modeling of horizontal infiltration using a given set of hydraulic functions and initial and final conditions. Then, the Boltzmann transformation must be computed by integrating the modeled water content profile. This procedure is often timeconsuming 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: where D(θ ) = K(θ ) dh dθ 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)): where h is the water pressure head, h 0 and h 1 are respectively the initial and final water pressure heads, and θ 0 = θ (h 0 ) and θ 1 = θ (h 1 ). For the sake of clarity, the functions S 2 D and S 2 K are respectively referred to as the "diffusivity" and "conductivity" forms of sorptivity. S 2 D and S 2 K are equivalent so long as the water retention function θ (h) is bijective over the water pressure head interval [h 0 , h 1 ], which is the case when the water pressure head at the surface is lower than the airentry water pressure head, i.e., h 1 ≤ h a . Then, Eq. (3) can be easily deduced from Eq. (2) by a simple change of variable from θ to h: Otherwise, when the surface water pressure head exceeds the soil air-entry pressure, i.e., h 1 > h a (≤ 0), sorptivity must be computed using Eq. (3) (Ross et al., 1996). Indeed, the function θ (h) is no longer bijective over the interval [h 0 , h 1 ]. Consequently, Eqs. (2) and (3) are no longer equivalent. In this case, the integration involved in Eq. (3) must be divided into two parts, one part integrating over the interval [h 0 , h a ] ensuring the bijectivity of the function θ (h) and the other part retaining the integration interval [h a , h 1 ] corresponding to the saturated part of the integration (Ross et al., 1996). Then, the change of variable from θ to h in the first integral leads to the retrieval of S 2 D : The computation of sorptivity with Eq.
(2) or Eq. (3) requires a set of hydraulic functions to be chosen from the wide range of available models. Here, we considered five of the most widely used hydraulic models. Firstly, we considered the delta model (d, delta) , which involves Dirac delta functions (stepwise functions) for the description of both the water retention (WR) and hydraulic conductivity (HC) functions. Indeed, this model is often considered for analytical resolutions to the Richards equation and the determination of analytical expressions for water infiltration, like the Green and Ampt approach (Triadis and Broadbridge, 2012). Secondly, we considered the Brooks and Corey (1964) model, referred to as the BC model, since it is among the first hydraulic models of soil physics (Hillel, 1998). The BC model involves power functions for both the WR and HC functions, thus allowing for analytical integration of Eq. (2) and leading to analytical expressions for sorptivity (e.g., Varado et al., 2006). Thirdly, the van Genuchten-Burdine (vGB) model was studied since it has been used for the development of the BEST methods (Beerkan Estimation of Soil Transfer functions) for the characterization of soil hydraulic properties (Lassabatere et al., 2006;Yilmaz et al., 2010;Bagarello et al., 2014). The vGB model combines the van Genuchten (1980) model with Burdine's condition m = 1 − 2 n for the WR function and the Brooks and Corey (1964) model for the HC function. Fourthly, we considered the van Genuchten- Mualem (vGM) model that combines the van Genuchten (1980) model with Mualem's condition m = 1 − 1 n for the WR function and the Mualem (1976) capillary model for the HC function. The vGM is among the most widely used models and often used for the numerical modeling of flow in the vadose zone (Šimůnek et al., 2003). Lastly, the Kosugi (KG) model (Kosugi, 1996), was also considered since it relates the water retention function to physical characteristics of the soil pore size distribution assuming log-normal distributions.
These five models have the following mathematical expressions (Triadis and Broadbridge, 2012;Brooks and Corey, 1964;van Genuchten, 1980;Mualem, 1976;Kosugi, 1996;Nasta et al., 2013): delta model: BC model: vGB model: vGM model: where H stands for the one-sided Heaviside step function: H (x < 0) = 0, H (x ≥ 0) = 1 (Triadis and Broadbridge, 2012); and "erfc" stands for the complementary error function. These models involve several specific hydraulic shape parameters and the following common scale hydraulic parameters: residual water content, θ 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 * 2 K (−∞, 0), which we show is equal to the parameter c p , and (ii) the product of scale parameters and correcting factors accounting for the contribution of initial water contents. The square scaled sorptivity corresponds to the sorptivity of a unit soil (unit value for all the scale parameters, except the residual water content fixed at zero) and for the whole range of water pressure heads, i.e., (−∞, 0]. It depends only on the soil hydraulic shape parameters, and its determination features the main algebraic complexity of the whole scaling procedure; the rest relies on simple algebraic operations (multiplication and sums). It is then computed for the hydraulic models defined by Eqs. (6)-(10), either analytically when feasible or numerically, otherwise. For each model, it is computed and tabulated as a function of a shape index that characterizes the spreads of the water retention function (from gradual to stepwise shapes, corresponding to soils with a broad or a very narrow pore size distribution, respectively). The evolution of the square scaled sorptivity versus the shape index is compared and discussed between models. In the last section, we illustrate the application of the proposed scaling procedure. We show how the tabulated values of the square scaled sorptivity c p can be used to upscale sorptivity and easily provide the sorptivity corresponding to a zero water pressure head at the surface for relatively small initial water contents.

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 2 K (h 0 , h 1 ) to the maximum square sorptivity S 2 K (−∞, 0) by isolating the effect of initial and final conditions (h 0 and h 1 ) and (ii) scaling hydraulic functions and sorptivity to split the contributions of shape and scale parameters and facilitate the computation of S 2 K (h 0 , h 1 ). In our case, we consider that the final conditions always involve a null water pressure head, h 1 = 0, and saturated conditions, θ 1 = θ s , while initial conditions (h 0 , θ 0 ) may vary.

Isolating the effect of initial conditions
The first step of the proposed procedure involves the work of Haverkamp et al. (2005), who isolated the contributions of the initial and final conditions to sorptivity. These authors considered that, at dry initial states (θ 0 ≤ 1 4 θ s ) and a zero water pressure head at surface (i.e., h 1 = 0 and θ 1 = θ s ), the fol-lowing approximation applied: where K 0 corresponds to the initial hydraulic conductivity K 0 = K(θ 0 ). The parameter c p = S 2 (0,θ s ) |h g |K s θ s is a proportionality constant that depends only on the hydraulic shape parameters (Haverkamp et al., 2005). By combining both expressions in Eq. (11), the sorptivity, S 2 (θ 0 , θ s ), can be defined as the product of three different terms that account for the respective contributions of the shape parameters (lumped into the proportionality constant c p ), the scale parameters, |h g |, K s , and θ s , and initial conditions, θ 0 : Equation (12) offers an accurate and practical approximation for the computation of sorptivity in the case of Beerkan runs, i.e., for a zero water pressure head imposed at the soil surface, h 1 = 0. Equation (12) was frequently used for the treatment of Beerkan data and in particular in all BEST methods (Lassabatere et al., 2006;Yilmaz et al., 2010;Angulo-Jaramillo et al., 2019). However, it addresses the case of soils with null residual water content, θ r = 0, and without any airentry 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: Indeed, the denominator θ s must be replaced by (θ s −θ r ) when θ r = 0 to ensure that the ratio S 2 (θ 0 ,θ s ) S 2 (θ r ,θ s ) tends towards unity when θ 0 → θ r . Secondly, we consider that the approximation behind Eq. (11) involves only the unsaturated part of sorptivity, i.e., S 2 D (θ 0 , θ s ). As mentioned above, when the air-entry water pressure head is non-null, the computation of sorptivity S 2 K (h 0 , 0) must be split into its unsaturated and saturated parts, as illustrated by Eq. (5).
The following derivations are then proposed: where S 2 D (θ r , θ s ) and 2(θ s − θ r )K s |h a | refer to the unsaturated and saturated parts of the square sorptivity S 2 K (−∞, 0) (see Eq. 5 with h 0 = −∞, i.e., θ 0 = θ r , and h 1 = 0). Equation (14) can be simplified by introducing correcting factors R θ and R K to account for the contribution of the initial conditions: with R θ and R K defined as follows: 2.1.2 Scaling sorptivity and defining parameters c p and c p So far, the computation of the square sorptivity S 2 K (h 0 , 0) has been treated considering dimensional equations. The second step of the proposed procedure scales the sorptivity in order to alleviate any numerical difficulty in relation to the values of dimensional variables. The square dimensional sorptivity S 2 K (h 0 , 0) can be easily related to the square scaled sorptivity S * 2 K (h * 0 , 0) by scaling variables (water content, water pressure head, and hydraulic conductivity) as follows (Ross et al., 1996): This scaling procedure defines the dimensionless water retention, S e (h * ), the dimensionless (or relative) hydraulic conductivity, K r (S e ), and the dimensionless hydraulic diffusivity function D * (S e ) = K r (S e ) dh * dS e . These dimensionless hydraulic functions define the hydraulic characteristics of a unit soil that has the same values for the shape parameters and the unit value for all the scale parameters, θ 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 2 , to the square scaled sorptivity (unit soil), S * 2 (Ross et al., 1996): Then, S * 2 can be computed by applying Eqs.
(2) or (3) to the dimensionless hydraulic functions as a function of the initial and final water pressure heads, h * 0 and h * 1 , or saturation degrees, S e,0 = S e (h * 0 ) and S e,1 = S e (h * 1 ): S * 2 D and S * 2 K define the diffusivity and conductivity forms of the square scaled sorptivity and are related to each other by Eqs. (4) and (5), leading to with h * a = h a |h g | . The scaled sorptivity functions, S * 2 D and S * 2 K , 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 2 K (−∞, 0) and its unsaturated part S 2 D (θ r , θ s ), which are involved in Eq. (15), to define the parameters c p and c p : The application of Eq. (19) with the proper lower and upper initial and final conditions (S e,0 = 0, h 0 = −∞, S e,1 = 1, h 1 = 0) leads to the following expression for the two parameters c p and c p : The parameters c p and c p depend exclusively on the soil shape parameters. The application of Eq. (20) provides a very simple equation linking the two parameters:

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 * 2 K (h 0 , 0): 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: 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 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.

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: BC model: vGM model: Note that the scaling parameter for the water pressure head, h g , used in Eqs. (6)-(10) was set equal to the air-entry pressure head for the delta and BC WR functions; i.e., h d and h BC were set equal to h a , as mentioned above.
Equations (27)-(31) were then used to derive the following formulations for the dimensionless diffusivity functions, D * (S e ), applying D * (S e ) = K r (S e ) dh * dS e (see Appendix A): where δ stands for the one-sided Dirac delta function (Triadis and Broadbridge, 2012).

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

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.

c p for the delta model
This case is the easiest one. Indeed, the HC function is characterized by a null hydraulic conductivity for h * < −1 and a unit value K r = 1 for h * ≥ −1, as featured by Eq. (27). In this case, Eq. (22) yields 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.

c p for the Brooks and Corey (BC) model
The BC model involves an air-entry water pressure head h * a = −1. We used Eq. (22) while accounting for the saturated part of sorptivity and used the diffusivity form for the determination of the unsaturated sorptivity. Indeed, the diffusivity function D * BC (S e ) shown in Eq. (33) obeys a power law, which makes it possible to integrate the diffusivity form of sorptivity analytically. Simple algebraic operations and integrations of Eq. (22) lead to the following equation (see demonstration in Appendix B1): When the relation between η BC and λ BC is ruled by Eq. (37) with p = 1, the analytical expression of c p turns into Our results are in line with previous studies (Varado et al., 2006).

c p for the van Genuchten-Burdine (vGB) model
In the case of the vGB model, there is no air-entry pressure head, i.e., h * a = 0. Equation (22) shows that c p reverts to the diffusivity form of the square scaled sorptivity, 1 0 (1 + S e )D * (S e )dS e . In addition, the diffusivity function, D * vGB (S e ) (Eq. 34), makes the square scaled sorptivity analytically integrable, leading to (see demonstration in Appendix B2) where is the gamma function: Considering the relations between m and n, i.e., m = 1− 2 n , and the relation between η and λ = mn in Eq. (37) with p = 1, the following simplification emerges: The expression corresponding to Eq. (45) was already proposed and discussed by Haverkamp et al. (2005).

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): Note that this equation requires that m vGM = 1 1+l vGM and m vGM = 1 2+l vGM . Considering that shape parameter l vGM = 1 2 , as usually considered, Eq. (48) can be simplified to 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.

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: where h * KG 1 2 is the water pressure head corresponding to S e = 1 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 1 2 .

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): The equation λ BC = λ vGB = n vGB − 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 vGB = 1 − 2 n vGB leads to Since m vGB is the appropriate shape index for the vGB model, we consider its equivalent, λ BC 2+λ BC , to be the appropriate shape index for the BC model, leading to 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 : 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: 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) 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 η = 2 λ + 2 + p and l vGM = l KG = 1 2 and diffusivity functions Eqs. (38)-(41). Lastly, we computed the square scaled sorptivity c p as a function of the shape index, Eq. (56), and discussed the function c p (x) regarding the choice of the hydraulic model. The values of c p (x) are also compared to those of the delta model (Dirac delta functions), i.e. c p,d = 2.

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 . 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 (−∞, 0] (Fig. 1a, vGB, vGM, and KG), while the BC model reaches the asymptote S e = 1 at h * = −1 with full saturation for h * ≥ −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 * = −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 * ≤ −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 in-creases at both low and high saturation degrees, with a lower increase at intermediate saturation degrees (Fig. 1d, vGB). The two other models, vGM and KG, exhibit the same type of S shape, with an infinite limit close to saturation (Fig. 1d, vGM and KG). Such an infinite limit spoils the numerical integration of Eq. (22) for the determination of c p , requiring the use of the mixed formulation for the KG models defined by Eqs. (50) and (56).
Varying the shape index changes the WR and HC 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 * ≤ −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 * ≤ −1, followed by an increase for h * ≥ −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 * = 0 ( Fig. 2c  and g). Conversely, for the vGM and KG models, K r (h * ) tends towards zero in the vicinity of h * = 0 ( Fig. 2k and o, inverted arrows). Similarly, the dimensionless diffusivity, D * (S e ), tends towards zero over the whole interval [0,1] when the shape index tends towards zero ( Fig. 2l and p, inverted arrows). Consequently, the features of K r (h * ) and D * (S e ) functions suggest that the square scaled sorptivity is very small for vGM and KG models, when the shape index tends towards zero (see results below).

Square scaled sorptivity c p as a function of shape indexes
The square scaled sorptivity, c p (x), is plotted as a function of the shape index, x (Fig. 3). The results show contrasting evolution. For the BC model, we note a decrease of c p,BC (x) from 4 down to 2 (Fig. 3a, c p ). Such a decrease is expected since K r,BC (h * ) functions decrease over the whole inter- The air-entry water pressure head is non-null, so that the c p and c p are related by c p = c p + 2 h * a = c p + 2, since by convention h * a = −1. Consequently, c p decreases from 2 to 0, with a simple vertical shift compared to c p (Fig. 3a, c p ). For the other hydraulic models, the air-entry water pressure head is null, so that c p = c p (Fig. 3b-d). In the following, we compare the evolution of c p as a function of the shape index for each model.
In contrary to the BC model, c p,vGB (x) for the vGB model does not decrease monotonically (Fig. 3c). Instead, c p,vGB (x) decreases to 1.5 before increasing up to 2 for x ≥ 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 * ≤ −1 and an increase for h * ≥ −1 (Fig. 2g, arrows). The numerical computation sorts out these contrasting effects and demonstrates the two-step variation, i.e., a decrease followed by an increase. The two boundaries of the function c p,vGB (x) can be easily found using Eq. (56) with a lower limit of c p,vGB (0) = 2 3 2 1 2 = 1 2 2 = π and an upper limit of c p,vGB (1) = (1) (3) (3) + (4) (4) = 2. For the KG and vGM models, the trend is opposite, and the functions c p,KG (x) and c p,vGM (x) both increase (Fig 3b and d). This increase is in line with the fact that the shape index mainly increases the diffusivity function, = 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 ∈ [0, 0.3] (Fig. 3b).
The four functions c p,BC (x), c p,vGB (x), c p,vGM (x), and c p,KG (x) all reach the value of 2 when the shape index approaches unity, i.e., when the WR and HC functions tend towards stepwise functions. In fact, the value of c p (x) converges to the value obtained for the delta model, c p,d = 2 (see Eq. 42). Such a result indicates that similar results should 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 twostep 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 2 K (h 0 , 0), for a given soil and given initial conditions, h 0 . 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, (θ s − θ r )K s |h g |, 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, |h g |, is expected to compensate for the very low values of c p when the vGM and KG models are used. We conclude that the value of |h g | must tends towards being infinite when the shape index tends towards zero for these two models. For the KG model, such a relation, |h KG | ∼ x KG , is related to the relation |h KG | ∼ σ KG since x KG = (1 + σ KG ) −1 . Our statements imply that very large values of σ KG (x KG → 0) should be associated with very large values of |h KG |, i.e., with a very small pore radius. In other words, fine (coarse) soils with small (large) pores should have broad (narrow) pore size distributions. These considerations are in line the previous studies on the KG model by Pollacco et al. (2013) (see Fig. 1 of their paper) and Fernández-Gálvez et al. (2019). These authors even related the pore radius of soil to the standard deviation of pore size distribution with a strongly decreasing function. Similar trends relating to some extent a relation between pore mean and pore standard deviation should also apply to the vGM model given our findings (and to avoid null sorptivity for soils with broad pore size distributions). More investigations are needed to verify such a link between average pore size and standard deviation for real soils. More investigations are also needed to guide on the proper choice of the hydraulic models as a function of the soil type, as already suggested by Fuentes et al. (1992). These aspects will be the subject of a specific study.
Regarding the numerical accuracy of computed values of c p , we used analytical formulations for the BC, vGB, and vGM models, as detailed in Eq. (56). These values are expected to be perfect without any error since they correspond to the application of exact analytical formulations. Instead, we used the mixed numerical formulation defined by Eq. (50) for the KG model that relies on the numerical integration of the HC or diffusivity functions. In that case, the numerical integration may bring some numerical errors. The mixed form, Eq. (50), was designed to minimize numerical indetermination and uncertainty. Such a formulation was applied to the other models (BC, vGB, and vGM), and the resulting values were compared against the analytical formulations (considered to be the benchmark). A perfect agreement was obtained (errors < 1 %), thus validating the numerical mixed formulation and reinforcing confidence in the values in Table 1. Note that the promotion of the numerical mixed formulation, Eq. (50), and its uncertainty will be the subject of another study.
3.3 Upscaling sorptivity S K (h 0 , 0) from c p In this section, we elaborate on the use of Eq. (26) for the easy and straightforward computation of S 2 K (h 0 , 0) from the tabulated values of c p ( Table 1). The proposed scaling procedure Eq. (26) allows for the computation of S 2 K (h 0 , 0) given initial conditions (water contents or water pressure heads), hydraulic shape and scale parameters, and specific hydraulic models selected among Eqs. (7)-(10): 1. Use the shape parameter (λ 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 KG = 1 1+σ KG , and x BC = λ BC 2+λ 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, θ 0 or θ (h 0 ), depending on the description of the initial condition (either θ 0 or h 0 ).
6. Compute the scaled air-entry water pressure head, |h * a | = | h a h g |.
7. Compute the square scaled sorptivity, S * 2 8. Upscale to derive the square sorptivity, As an illustrative example, let us consider the case of a loamy soil at water saturation with a slightly positive water pressure head at the surface (h 1 = 0) and an initial water pressure head of h 0 = −10 m (dry conditions). The loamy soil has the features of loam as defined in the database of Carsel and Parrish (1988). Its WR and HC functions are described by the vGM model, with the following shape and scale parameters: θ r = 0.078, θ s = 0.43, h g = −277 mm, K s = 2.88 × 10 −3 mm s −1 , n vGM = 1.56, and l vGM = 0.5. The application of the step-by-step procedure gives the following results: 1. shape index: x = m vGM = 1 − 1/n vGM , leading to x = 0.359; 2. corresponding value of c p in Table 1 (vGM model col 3. initial water content: computed from the initial water pressure head of −10 m using the vGM-WR function, i.e., Eq. (9): θ 0 = 0.125; 4. initial hydraulic conductivity: computed from the initial water content using vGM-HC function, i.e., Eq. (9): K 0 = 1.87 × 10 −9 mm s −1 ; 5. corresponding correction factors: R θ = θ s −θ 0 θ s −θ r and R K = K s −K 0 K 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 |h * a | = 0; 7. square scaled sorptivity: 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 g = −277 mm, and K s = 2.88 × 10 −3 mm s −1 , 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: 3. initial water content: computed from the initial water pressure head of −10 m using BC-WR function, i.e., Eq. (7): θ 0 = 0.125; 4. initial hydraulic conductivity: computed from the initial water content using BC-HC function, i.e., Eq. (7): K 0 = 3.342 × 10 −9 mm s −1 ; 5. corresponding correction factors: R θ = θ s −θ 0 θ s −θ r and R K = K s −K 0 K 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 |h * a | = 1; 7. square scaled sorptivity: 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 2 K (h 0 , 0). It must be noted that the proposed scaling procedure applies only for a dry initial state. Indeed, Haverkamp et al. (2005) stated that their approximation, Eq. (11), was valid only when θ 0 ≤ 1 4 θ s . For fine soils, even a small initial wa-ter pressure head may cause θ 0 > 1 4 θ 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 g = −2000 mm, K s = 5.555 × 10 −5 mm s −1 , n = 1.09, and l = 0.5. Considering the same value for the initial water pressure head, i.e., h 0 = −10 m, the initial water content is θ 0 = 0.318, which exceeds 1 4 θ s . The application of the scaling procedure led to an estimated sorptivity of 0.0475 mm s − 1 2 , whereas the targeted sorptivity computed with Eq. (2) was 0.0127 mm s − 1 2 . Such an error corresponds to an overestimation by a factor of 2.73. Thus, we advise that the user verify that θ 0 ≤ 1 4 θ s when θ r = 0, or more generally S e,0 ≤ 1 4 before using the proposed scaling procedure.

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 θ 0 ≤ 1 4 θ s when θ r = 0, or more generally S e,0 ≤ 1 4 . In addition to providing a straightforward method for the determination of sorptivity, this study brings very interesting findings on the square scaled sorptivity, c p , and its dependency upon the shape index, x, and the chosen hydraulic model. The results show that the function c p (x) strongly depends on the hydraulic model selected for the WR and HC curves. If all the functions c p (x) converge for the same value, i.e., 2, close to x = 1 (stepwise WR functions -narrow pore size distribution), they strongly divert close to x = 0 (graded WR functions -broad pore size distribution), with values of 0 for vGM and KG models versus 3-4 for the vGB and BC models. However, the sorptivity should remain the same regardless of the selected hydraulic model: one soil under particular initial conditions means one single sorptivity. Consequently, the contrast of scaled sorptivity must be compensated by a contrast in scale parameters. However, among scale parameters, the residual and saturated water contents and the saturated hydraulic conductivity cannot be changed between models, since they characterize the dry and saturated states of the same soil. Consequently, the value of the scale parameter, h g , must be the one to compensate. Previous studies on the Kosugi model have already hypothesized a strong relation between the scale parameter, h KG , and the standard deviation, σ 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 .
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 * (S e ) = K r (S e ) dh * dS e . This task requires us to first derive the inverse functions for the dimensionless water retention curves. The following equations can be easily found through usual algebraic developments: where erfc −1 is the inverse function of the complementary error function. These functions can be differentiated to define their relative derivatives, dh * dS e : The differentiation of the function (f o g) = f o g · g and f −1 = 1 f o f −1 , considering bijective functions. We also use the usual derivative of the "erf" function, erf (x) = 2 √ π e −x 2 , and the relation between erfc and erf functions, The scaled hydraulic conductivity functions are Note that for the hydraulic conductivity function of the vGM model, K r,vGM , we distributed the terms according to (a + b) 2 = a 2 + 2ab + b 2 . The multiplication of Eqs. (A5)-(A8) by the expressions of relative conductivity, Eqs. (A9)-(A12), leads to the following expressions of dimensionless diffusivity: For the BC model, we need to account for the air-entry pressure, h * a = −1. We remind the reader that by convention, the scale parameter for the water pressure head, h BC , is equal to the air-entry water pressure head, h a . We then use Eqs. (22) and (23) The final expression can be easily computed by combining Eqs. (B1) and (B2): The concatenation of Eq. (B3) with the capillary model, Eq.( 37), leads to the following final expression, assuming p = 1: 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 * a = 0, leading to We now use the beta function B, which has the following properties (assuming x > 0 and y > 0): where the function is already defined by Eq. (46): (z) = +∞ 0 t z−1 e −t dt (z > 0). Then, the parameter c p,vGB can be expressed as follows: The last equation uses the fact that z (z) = (z + 1). Equation ( Note that the proposed simplification using the beta function, Eq. (B9), requires that mη + m 2 − 1 2 = 1 2 + 3 2 + p m > 0, which is quite evident since m > 0. When Eq. (B9) is combined with the capillary model, i.e., Eq. (37), the expression becomes (1 + 2m) + The last term of A can be simplified easily: In this case, we need to assume that ml + m − 2 > −1, i.e. m > 1 l+1 , to use the beta function. Such a condition corresponds to m > 2 3 , for the by-default value of l = 1 2 . The two first integrals can then be replaced using the beta and the gamma functions, leading to the final expression for part A Code availability. Note all computations were done using Scilab free software. The scripts for the computation of Eqs. (28)-(31) for the computation of WR and HC functions, Eqs. (38)-(41) for the computation of the dimensionless diffusivity, and Eqs. (56) for the computation of the c p parameter can be downloaded online at https://doi.org/10.5281/zenodo.4587160 (Lassabatere, 2021). Note also that the computation of the proposed sorptivity was implemented in the SoilWater-ToolBox software that interrelates specific modules to derive the soil hydraulic parameters using a wide range of cost-effective methods, accessible online at https://github. com/manaakiwhenua/SoilWater_ToolBox/ (last access: 7 September 2021, Pollacco and Fernández-Gálvez, 2021) (open source under the GP-3.0 License).
Data availability. No data sets were used in this article.
Author contributions. LL established the question, performed the analytical developments, computed the numerical results, and provided the first draft of the manuscript. PEP performed the analytical developments and verified the numerical results. DY, BL, DMF, SdP, and MR verified parts of the developments and numerical computations. JP and JFG helped with the use of the Kosugi model. RDS and MAN helped with the editing and the presentation of the results, in particular the writing of the application of the proposed approach for end users. All the authors contributed to the editing of the manuscript.
Competing interests. The authors declare that they have no conflict of interest.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Financial support. This research has been supported by the Agence Nationale de la Recherche (grant no. ANR-17-CE04-010).
Review statement. This paper was edited by Nunzio Romano and reviewed by two anonymous referees.