Scaling procedure for straightforward computation of sorptivity

Sorptivity is a parameter of primary importance in the study of unsaturated flow in soils. This integral parameter is often considered for modeling the computation of water infiltration into vertical soil profiles (1D or 3D axisymmetric geometry). Sorptivity can be directly estimated from the knowledge of the soil hydraulic functions (water retention of hydraulic conductivity), using the integral formulation of Parlange (Parlange, 1975). However, it requires the prior determination of the soil 5 hydraulic diffusivity and its numerical integration between the initial and the final saturation degrees, which may be tricky for some instances (e.g., coarse soil with diffusivity functions quasi-infinite close to saturation). In this paper, we present a specific scaling procedure for the computation of sorptivity considering slightly positive water pressure heads at the soil surface and initial dry conditions (corresponding to most water infiltration on the field). The square sorptivity is related to the square dimensionless sorptivity (referred to as cp parameter) corresponding to a unit soil (i.e., unit values of all the scaled 10 parameters and zero residual water content) utterly dry at the initial state and saturated at the final state. The cp parameter was computed numerically and analytically for five current models: delta functions (Green and Ampt model), Brooks and Corey, van Genuchten-Mualem, van Genuchten-Burdine, and Kosugi models as a function of the shape parameters. The values are tabulated and can be easily used to determine any dimensional sorptivity value for any case. We propose brand-new analytical expressions for some of the models and validate previous formulations for the other models. Our numerical results also showed 15 that the relation between the cp parameters and shape parameters strongly depends on the chosen model, with either increasing 1 https://doi.org/10.5194/hess-2021-150 Preprint. Discussion started: 24 March 2021 c © Author(s) 2021. CC BY 4.0 License.


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 modelling 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., 2009Lassabatere et al., , 2014Lassabatere et al., , 2019 and 25 must be properly quantified from the soil hydraulic functions and the initial and surface hydric conditions. 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) for the modelling of 1D gravity-free water infiltration: (1) 30 In the above equations, S(θ 0 , θ 1 ) stands for the sorptivity between θ 0 and θ 1 , χ = χ (θ) = x(θ) √ t stands for the Boltzmann transformation variable, θ 0 is the initial water content, and θ 1 the final water content corresponding also to the water content applied at the soil surface. The main shortcoming is that such a procedure requires to proceed to the numerical modeling of the horizontal infiltration, for given hydraulic functions and initial and final conditions. Then, the Boltzmann transformation must be computed from the modeled water content profile and integrated between the initial and final water contents (see Eq. (1)). 35 Such a procedure may be time-consuming and subject to numerical instabilities spoiled with numerical errors.
To avoid such a complexity, Parlange (1975) proposed a formulation that relates directly sorptivity to the hydraulic functions and the initial and final water contents: where D(θ) = K(θ)dh/dθ is the hydraulic diffusivity function. While the above equation provides the diffusivity form for 40 sorptivity determination, it can be equally defined as a function of the hydraulic conductivity function, K(h) = K (θ (h)):

55
The computation of sorptivity with Eqs.
(2) or (3) requires the choice of a set of hydraulic functions from a wide range of available models for the characterization of soils. Here, we consider five of the most widely used hydraulic functions. Firstly, the Brooks and Corey (1964) model, referred to as the "BC" model, is among the first hydraulic functions of soil physics (Hillel). The BC model involves power functions for both the water retention (WR) and hydraulic conductivity (HC) functions, thus allowing for analytical integration of Eq.
(2) and leading to analytical expressions for sorptivity (e.g. Varado et al., 2006). 60 Secondly, the van Genuchten -Burdine (vGB) model that combines van Genuchten (1980) model with Burdine conditions for water retention (WR) and the Brooks and Corey (1964) model for Hydraulic Conductivity (HC). The vGB model has been used for the development of BEST methods and for the characterization of soil hydraulic properties (Lassabatere et al., 2006;Yilmaz et al., 2010;Bagarello et al., 2014). Thirdly, the van Genuchten -Mualem (vGM) model that combines van Genuchten widely-used models and often used for the numerical modelling of flow in the vadose zone (Šimůnek et al., 2003). Fourthly, the Kosugi (KG) model (Kosugi, 1996), is quite popular since it relates the water retention function to physical characteristics of the soil pore size distribution assuming log-normal distributions. Finally, the Delta model (d), corresponding to the Dirac delta functions, i.e., the use of stepwise functions for the description of WR and HC, was added to the list of studied models.
Indeed, this model is often considered for the analytical resolution of Richards equations and the determination of analytical Eq.
(2) may pose numerical shortcomings for infinite hydraulic diffusivity, which is the case of 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 100 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 theory section details the scaling procedure that relates the square sorptivity to the square scaled sorptivity, c p = S * 2 K (−∞, 0), the product of scale parameters, and correcting factors accounting for the contribution of initial hydric conditions. The square scaled sorptivity corresponds to the sorptivity 105 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 head, 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 was then computed for the hydraulic models defined by Eqs. (6)-(10), either analytically when feasible or numerically, otherwise. For each model, it was computed and tabulated as a function of a shape index that characterizes the spreads of 110 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 was 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 provide easily the sorptivity corresponding to zero water pressure head at surface any relatively small initial water content. The global scaling procedures rely on several steps emanating from previous studies. Scaling has often been used to separate the effects of sorptivity drivers and simplify the computation of sorptivity. Ross et al. (1996) suggested to scale water content, water pressure head and hydraulic conductivity using the following equations: 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 * dSe . These dimensionless hydraulic functions define the hydraulic characteristics of the unit soil that has the same values for the shape parameters and unit value for all the scale parameters, θ s = 1, h g = 1, K s = 1, excepted the residual water content that is fixed at zero. The use of the 125 scaling Eqs. (11) allows us to relate the square sorptivity (dimensional soil), S 2 , to the square scaled sorptivity (unit soil), S * 2 , as follows (Ross et al., 1996): 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 squared scaled sorptivity and are related to each other by Eqs. (4) and (5), leading to: With h * a = 0 for the soils with a null air-entry water pressure head or h * a < 0 otherwise. Eqs. (13) and (14) show that the 135 dimensionless sorptivity functions, S * 2 D and S * 2 K depend only on the shape parameters and the initial and final conditions. Consequently, one advantage to the simplification approach is that it allows to separate the effect of different scale parameters.
To further simplify, (Haverkamp et al., 2005) isolated the contributions of the initial and final conditions to sorptivity. They considered that, at dry initial states (θ 0 ≤ 1 4 θ s ) and zero water pressure head at surface (i.e., h 1 = 0 and θ 1 = θ s ), the following approximation applied: Ks θs−θ0 θs Where K 0 corresponds to the initial hydraulic conductivity (K 0 = K (θ 0 )) and c p is a proportionality constant that depends only on shape hydraulic parameters (Haverkamp et al., 2005). By combining the two Eqs. (15), 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 hydraulic parameter c p ), the scale parameters |h g |, K s , θ s , and initial and final conditions: Equation (16) 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 (16) 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, 150 h a = 0. Besides, it was mostly developed and used for the vGB model.
In this study, we adapt Eq. (16) to any type of hydraulic models, including those with non-null residual water contents and air-entry water pressure heads. First of all, we consider that the residual water content θ r must be accounted for. We suggest to replace Eq. (15) with 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. (15) involves only the unsaturated part of sorptivity, i.e., S 2 D . 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). In that case, the following derivations can be proposed: In the right-hand side of Eq. (18), the terms S 2 D (θ r , θ s ) and 2(θ s − θ r )K s |h a | refers to the unsaturated and saturated parts of the sorptivity S 2 K (−∞, 0), given that the application of Eq. (5) to the case of h 0 = −∞ and h 1 = 0, leads to: Consequently, Eq. (18) can be rewritten as follows: Then, this simplifying approach (Eq. 20), can be combined with the scaling procedure proposed by Ross et al. (1996), i.e., Eq. (11) and Eq. (12) to give the following expressions: or, equivalently, The coefficient c p , as pioneered by Haverkamp et al. (2005), corresponds to the square sorptivity of the unit soil for infinite initial water pressure head and null water pressure head at surface, h * 0 = −∞ and h * 1 = 0. c p can be quantified by using Eq. (13) or Eq. (14) with the right initial and final conditions, i.e., (h * 0 = −∞, h * 1 = 0) and (S e,0 = 0, S e,1 = 1): Eqs. (22) is a new version of the approximation developed by Haverkamp et al. (2005), and can be applied to any type of hydraulic function. It separates the contribution of shape parameters (that lump in the term c p = c p − 2|h * a |), the scale parameters, involving the product (θ s − θ r )K s |h g |, and lastly the air-entry water pressure head, in its scaled version, |h * a |. The set of Eqs. (22) makes it very easy to compute the sorptivity corresponding to zero water pressure head at surface and any initial water pressure head h 0 or water content θ 0 , provided that the values of c p are known. The following part aims at computing 180 and tabulating c p for the hydraulic functions defined in Eqs. (6)-(10).

General expressions
The first steps of the determination of c p requires the computations of the dimensionless functions, S e (h * ), K r (S e ), and D * (S e ) to be injected in Eqs. (23). The application of the scaling Eqs. (11) to the hydraulic functions defined by Eqs. (6)-(10) 185 leads to the following expressions: Where H stands for the one-sided Heaviside step function. Note that the scaling parameter for water pressure head, h g , used in the 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.
Eqs. (24)-(28) 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 * dSe (see Appendix A): Where δ stands for the one-sided Dirac delta function (Triadis and Broadbridge, 2012). These expressions are demonstrated in 210 the appendix A.

Further simplifications
To reduce the number of shape parameters, we here propose several additional simplifications, that are usually considered (Angulo-Jaramillo et al., 2016). Several authors used capillary models to relate the unsaturated hydraulic conductivity to the water retention functions. For the vGB model, Haverkamp et al. (2005) linked the shape parameter related to the hydraulic 215 conductivity, η, with the combination of those of the water retention function, λ = m n as follows: Where the tortuosity parameter, p, takes the values of 1 for the case of the 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 result section). In addition, further simplifications involved the values of the tortuosity parameters, l vGM and l KG in the vGM and KG models. These 220 were fixed at the by-default value: l vGM = l KG = 1/2 (Šimůnek et al., 2003;Kosugi, 1996;Kosugi and Hopmans, 1998).
In practice, these shape parameters are rarely varied (Haverkamp et al., 2005). With those supplementary considerations, the diffusivity functions for BC and vGB models become (see demonstration in appendix A): The set of equations Eqs. (35)-(38) shows that the studied hydraulic functions and hydraulic diffusivity functions involve only one shape parameter, i.e., λ BC , m vGB , m vGM and σ KG for BC, vGB, vGM and KG hydraulic models, respectively.
In the following, we consider the both cases of generic Eqs.

Integral determination of parameter c p
Once the dimensionless diffusivity functions were determined, the use of Eq. (23) allowed the determination of c p , either analytically of numerically, depending on the considered hydraulic models.

c p for the Delta model
This case is the easiest for the determination of c p . Indeed, the hydraulic conductivity function is characterized by a null 240 hydraulic conductivity for h * < −1 and a unit value K r = 1 for h * ≥ −1, as featured by Eq. (24). The determination of c p then follows straight from Eq. (23), with the following steps: 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.

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

255
In this case, there is no air-entry pressure head, i.e., h * a = 0. Eq. (23) shows that c p reverts to the diffusivity form of the squared dimensionless sorptivity, 1 0 (1 + S e ) D * (S e )dS e . Besides, the diffusivity function Eq. (26) makes the dimensionless sorptivity analytically integrable, leading to (see demonstration in appendix B2): where Γ is the gamma function: Considering the relations between m and n Eq. (34) and the relation between η and λ = m n Eq. (34), the following simplification comes out: The expression corresponding to Eq. (42) was already proposed and discussed by (Haverkamp et al., 2005).
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. (45) can be simplified to: These sets of equations have never been proposed and constitute one of the inputs of this study. The complexity of algebraic developments for the derivation of Eqs. (42) and (44)

c p for the Kosugi (KG) model
No analytical formulation was found for the case of Kosugi's hydraulic functions. Therefore, the dimensionless scaled sorptivity was computed numerically with a generic procedure that can be applied to any type of model for WRHC 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 version of Eq. (47), c p is composed of two integrals of continuous functions over closed intervals, and thus integrable at all times. Again, a simplified version is proposed assuming that the shape parameter l KG is fixed at: l KG = 1 2 .
2.4 Shape indexes for comparing c p between the selected hydraulic models 285 The approach described below allows c p determination for the selected BC, vGB, vGM, and KG models. The dependency of c p upon the chosen hydraulic model could be questioned. Indeed, in addition to its contribution to simplifying Eq. (22), c p has real physical meaning: it corresponds to the sorptivity of unit soils for the case of zero water pressure head at the soil surface and utterly dry initial profile. Therefore the dependency of such sorptivity upon the selected hydraulic model has also to be questioned. Consequently, we designed shape indexes that allow comparing c p between models. These shape indexes 290 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 function between two extreme states: (i) a value close to zero for gradual WR function corresponding to soils with broad pore size distributions, (ii) a value of unity for stepwise WR function 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 be-295 tween a gradual shape (m = 0) and a stepwise function for (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 converge towards a power function similar to the BC model (Haverkamp et al., 2005): The equation λ BC = n vGB − 2 defines a relation between λ BC and then n vGB that ensures a similar state for WR functions.
Substituting n vGB by 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 , as 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. (50), i.e. by using a ratio, we propose the following shape index, Finally, the 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]. Equations (52) provide 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 shape parameters and indexes, c p can easily be related to the shape index using 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. (52)  3 Results

Analysis of the hydraulic functions and diffusivity
The hydraulic functions and diffusivity functions are plotted in Fig. 1 for the shape index value of 0.275, and their sensitivity 330 upon each model's 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. (13). 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 squared scaled sorptivity c p .

335
The comparison between the hydraulic models (at 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 340 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, (Fig. 1a).

345
Regarding the relative hydraulic conductivity, the BC and vGB models have similar shapes for K r (S e ), with both typical of power functions (Fig. 1b, "BC and "vGB"). In contrast, the vGM and KG models depict an inflection point, with 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 350 drops in hydraulic conductivity close to saturation that are typical for 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 depicted above. The function K r (h * ) exhibit 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 355 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). were computed as a function of x using Eq. (52) with lvGM = lKG = 1 2 . The dashed line represents the "delta" model.
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. (35) at S e = 1 (Fig. 1d, "BC"). In opposite, the vGB model diverts from the concave shape to tend towards infinity close to saturation, 360 at S e = 1 (Fig. 1d, "vGB" versus "BC"). The vGB model defines a S shape that reflects the larger increases at both low and high saturation degrees with 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"). Varying the shape index changed the WRHC 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 makes the WR functions close to a stepwise function corresponding to the Delta model (Fig. 2, 1 st column, arrows). As for the WR functions, the increase in the shape index put the curves K r (h * ) close to stepwise functions (Fig. 2, 3 rd 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 370 above (Fig. 2c, 3rd column). In opposite, 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,n, 3 rd 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, 1 st and 3 rd columns). This trend is less evident for the diffusivity functions (Fig. 2, 4 th column).
These results show that, regardless of the selected model, increasing the shape index put the hydraulic functions closer to 375 the Delta models that correspond to soils with narrow pore size distribution. In opposite, 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,g). In opposite, for the vGM and KG models, K r (h * ) tends towards zero in the vicinity of h * = 0 (Fig. 2k,o, inverted arrows). Similarly, the dimensionless diffusivity D * (S e ) tends towards zero over the whole interval 380 [0, 1] when the shape index tends towards zero (Fig. 2l,p, inverted arrows). Consequently, the features of K r (h * ) and D * (S e ) functions presume very small values of the square scaled sorptivity for vGM and KG models, when the shape index tends towards zero (see results below).

Squared dimensionless scaled sorptivity as a function of shape indexes
The squared dimensionless scaled sorptivity c p (x) is plotted as its function to the shape index, x (Fig. 3). The results show 385 contrasting evolutions. For the BC model, we note a decrease of c p,BC (x) from 4 down to 2 (Fig. 3). Such a decrease is expected since K r,BC (h * ) functions decrease over the whole interval (−∞, 0] with the shape index (see Fig. 1c). This leads to the decrease of the integral involved in Eq. (23), and thus c p . The upper and lower limits can be easily determined by applying Eq. (53), leading to c p,BC (0) = 4 and c p,BC (1) = 2.
In contrary to the BC model, c p,vGB (x) for the vGB model does not decrease monotonically (Fig. 3c). Instead, c p,vG (x) 390 decreases to 1.5 before increasing up to 2 for x ≥ 0.52. This feature is 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,vG (x) can be easily found by using Eq. (53) with a lower limit of c p,vGB (0) = 2Γ 3 2 Γ 1 2 = Γ 1 2 2 = π and an upper limit of 395 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,d).
This increase is in line with the fact that the shape indexes mainly increase the diffusivity function, D * (S e ) over the whole   remain constant regardless of the choice of hydraulic models 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, the choice of the hydraulic model strongly impacts the estimation of c p , for soils with broad pore size distributions, i.e., when the shape index gets close to zero: 420 the BC and vGB models predicts non-null values for c p , thus ensuring non-null values of the sorptivity (see Eq. (20)), 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. (20)). 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 parameters for water pressure head, |h g |, is expected to compensate for the very 425 low values of c p when the vGM and KG models are used. We conclude that the value of |h g | must tends towards infinite when the shape index tends towards zero for these two models. For KG models, such relation |h KG | x KG is related to the relation 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 vGM model given our findings (and to avoid null sorptivity for soils with broad pore size distributions). More investigations are needed to verify for real soils such link between average pore size and standard deviations. More investigations are also needed to guide on the proper choice of the hydraulic 435 models as a function of the type of soil, 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 Eqs. (53). These values are expected to be perfect without any error since their correspond to the application of exact analytical formulations. Instead, we used the mixed numerical formulation defined by Eqs. (53) for the KG 440 model that relies on the numerical integration of the hydraulic conductivity or diffusivity functions. In that case, the numerical integration may bring some numerical errors. The mixed form Eq. (47) was designed to minimize numerical indetermination and uncertainty. Such formulation was applied to the other models (BC, vGB and vGM) and the resulting values were compared against the analytical formulations (considered as the benchmark). A perfect agreement was obtained (errors < 1%), thus validating the numerical mixed formulation and making the authors very confident on the values tabulated in Table 1. Note that 445 the promotion of the numerical mixed formulation and the study of its uncertainty will be the subject of another study.

Upscaling sorptivity S
In this section, we elaborate on the use of Eq. (22) 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. (22) allows the computation of S 2 K (h 0 , 0) given initial hydric conditions (water contents or the water pressure heads), hydraulic shape and scale parameters, and specific hydraulic 450 models selected among the studied hydraulic models 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 . Table 1 the value of c p corresponding to the shape index, x, and the chosen WRHC functions.

Choose in
3. Consider or compute the initial water content θ 0 or θ (h 0 ) depending on the observation for the description of the initial condition (either θ 0 or h 0 ). 4. Compute the related hydraulic conductivity, K 0 = K (θ 0 ), using the HC function.
6. Compute the scaled air-entry water pressure head: |h * a | = | ha hg | 7. Compute the square scaled sorptivity S 2 * 8. Upscale to derive the square sorptivity: As an illustrative example, let consider the case of a loamy soil submitted to 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  4. Initial hydraulic conductivity: computed from the initial water content using vGM-HC function,i.e., Eqs. (9): K 0 = 1.8710 −9 mm s −1 .
5. Corresponding correction factors: R θ = θs−θ0 θs−θr and R K = Ks−K0 K0 , leading to: R θ = 0.865 and R K = 1.000 6. Air-entry water pressure head: no air-entry water pressure head, consequently, |h * a | = 0 475 7. Square scaled sorptivity: To check the accuracy of the proposed approximation, we computed the nominal value of the sorptivity, using the regular 480 Eq. (2). We found a very close value, with less than 0.5% relative error, demonstrating the accuracy of the proposed scaling procedure Eq. (22). As a second illustrative example, we consider the computation of sorptivity for the case of 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 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 λ = m n, as suggested by Haverkamp et al. (2005). The application of the proposed procedure leads to the following computations: 1. Shape index: x BC = λ BC / (2 + λ BC ) leading to a value of x BC = 0.219 2. Corresponding value of c p : c p = 2.678 (see Table 1, "BC model" column for x BC = 0.22) 490 3. Initial water content: computed from the initial water pressure head of -10 m using BC-WR function, i.e., Eqs. (7): 4. Initial hydraulic conductivity: computed from the initial water content using BC-HC function, i.e., Eqs. (7): K 0 = 3.342 10 −9 mm s −1 .
(3) and lead 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. (20) for the two cases of hydraulic functions with and without air-entry water pressure heads.

505
Equation (20) 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 dry initial state. Indeed, Haverkamp et al. (2005) stated that their approximation Eq. (15) was ensured only when θ 0 ≤ 1 4 θ s . For fine soils, even a small initial water pressure head may cause θ 0 > 1 4 θ s , which may spoil the proposed scaling procedure. To illustrate this point, we investigated the case of the silty clay soil, as defined by Carsel and Parrish (1988). This soil is defined for the following parameters: θ r = 0.07, θ s = 0.36, h g = −2000 mm, and 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 lead 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 error corresponds to an overestimation by a factor of 2.73. Thus, we advise that the user verify that θ 0 ≤ 1 4 θ s before using the proposed scaling procedure.

515
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 zero water pressure head imposed at surface and dry initial state (θ 0 ≤ 1 4 θ s ). We based our approach on the combination and 520 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 famous sets of hydraulic models: Brooks and Corey, van Genuchten -Mualem, van-Genuchten -Burdine 525 and Kosugi models. The values of c p were tabulated as function of specific shape indexes representing similar states of WR functions (well-graded versus stepwise shapes) between hydraulic models. The proposed scaling procedure is very easy of use.
Once a given hydraulic model is selected with related shape and scale parameters, the procedure steps are easy to perform: computation of the shape index from the 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 530 sorptivity from c p and these 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 ). In addition to providing a straightforward method for the determination of sorptivity, this study brings very interesting 535 findings on the square scaled sorptivity c p and its dependency upon the shape index, x and the chosen hydraulic models.
The results show that the function c p (x) strongly depends on the hydraulic model selected for the WRHC 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 540 hydraulic model: one soil submitted to peculiar initial conditions, 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 545 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 WRHC 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 hypothesis has never been suggested and requires further investigations. These results show the need to better understand the 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 analytically the scaled sorptivity 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 if of great interest and could be implemented into soil hydraulic characterization methods. For instance, additional BEST methods could be developed ob the basis of the 555 use of the proposed formulations for square scaled sorptivity to relate sorptivity to shape and scale parameters. In more details, the prior estimation of shape parameters allows the determination of the parameter c p using Eq. (45). Then, the estimation of saturated hydraulic conductivity, and sorptivity allows the determination of scale parameter h g once c p is determined. A similar procedure may be proposed for the vGM model, using the Eq. (46) that defines the parameter c p to the shape parameter m vGM .
The development of BEST method for the specific vGM hydraulic model, that is much more used than the vGB model, will 560 be the subject of further investigations. It would also be interesting to derive somehow the residual water content, and not to assume it to be equal to zero as it might alter the shape of the soil water retention function. The use of the scaled sorptivity for these purposes are the subject of ongoing studies. In the appendices, for the sake of clarity the notations of the shape parameters were 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 * dSe . This task requires first to derive the inverse functions for the dimensionless water retention curves. The following 570 equations can be easily found through usual algebraic developments: where erf c −1 is the inverse function of the complementary error function. These functions can be differentiated to define their relative derivatives, dh * dSe : The differentiation of the function h * KG involves the following usual rules of differentiation (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 function erf function, erf (x) = 2 √ π e −x 2 , and the relation between erf c and erf functions, erf c(x) = 1 − erf (x).

650
However, thanks to the analyticity of the functions involved in the expression, this equality remains valid for m > 1 1+l and can be considered for any value of m ∈ [0, 1] provided that m = 1 1+l and m = 1 2+l . The Eq. (B19) demonstrates the Eqs. (42). The combination of these equations with the capillary model Eq. (34) leads to Eqs. (46).