Technical note: A sigmoidal soil water retention curve without asymptote that is robust when dry-range data are unreliable

In a recently introduced parameterization for the soil water retention curve (SWRC) with a sigmoid wet branch and a logarithmic dry branch, the matric potential at the junction point of the sigmoid and the logarithmic branch (h_j) was a fitting parameter, while that at oven dryness (h_d) was derived from the fitting parameters. The latter is undesirable, especially if reliable data in the dry range are limited. Therefore, an alternative is presented in which shape parameter α instead of h_d is a derived parameter, and h_d can be fitted or fixed. The resulting relationship between α and h_j is such that it prevents correct fits for h_j. Fortunately, an expression for h_j is found that allows it to be replaced by α as a fitting parameter. The corresponding parameter space is well behaved and has fewer internal bounds defined by restraining relationships between parameters than the space for h_j as a fitting parameter. The few available values of h_j in the literature are in line with those according to the new expression. The reformulated SWRC is fitted to data of 21 soils by shuffled complex evolution. The paper gives the main features of an accompanying open source fitting code. The curves fit the data well, except for some clayey soils. A theoretical value of h_d performs well for a wide range of soils. The new SWRC simplifies to an earlier junction model of the SWRC based on a well known power-law SWRC if α is very large.

1 Introduction

Recently, de Rooij et al. (2021) proposed a closed-form expression for the SWRC (soil water retention curve) with a distinct air-entry value, like the SWRC proposed by Ippisch et al. (2006), a sigmoid shape in the intermediate range according to van Genuchten (1980), and a logarithmic dry branch terminating at a finite matric potential h_d (L) at which the soil was oven dry, with the water content essentially zero. The volumetric water contents and derivatives of the sigmoid and logarithmic branches were matched at the matric potential of their junction according to Rossi and Nimmo (1994). The rationale for developing the function was to preserve the desirable sigmoid shape of van Genuchten's (1980) curve while removing the physically unrealistic asymptote at some non-zero residual water content (Du, 2020), eliminating the non-converging integral of the SWRC for commonly occurring parameter values (Fuentes et al., 1991), and avoiding the detrimental effect of the non-zero slope at saturation on hydraulic conductivity near saturation (Durner, 1994; Assouline and Or, 2013; Wang et al., 2022).

Erroneous measurements in the dry range can lead to unrealistically low values of the matric potential at oven dryness, h_d (de Rooij et al., 2021). Unfortunately, h_d was not fitted independently but expressed as a function of other parameters that were fitted. Data in the dry range can be unreliable due to lack of equilibrium or other causes (Bittelli and Flury, 2009; Solone et al., 2012). If they are, it would be helpful to fix h_d at a reasonable value, e.g., −10^6.8 cm H₂O (Schneider and Goss, 2012), and either give the unreliable data points a lower weight during the fitting process or remove them altogether. An SWRC with improved behavior in the dry range even if dry-range data are scant can help improve the behavior of Richards' solvers (e.g., SWAP Soil Water Atmosphere Plant, 2022; Šimůnek et al., 2016), be useful to improve conceptualizations of the soil reservoirs in large-scale hydrological models (e.g., Lawrence et al., 2019), and for investigating dielectric properties of dry soil and associated soil backscatter (Ferré and Topp, 2002; Davis and Annan, 2002).

This note presents an alternative to de Rooij et al.'s (2021) model in which h_d is a fitting parameter. In doing so, it uncovers the peculiar behavior of shape parameter α (L⁻¹), which makes it essentially impossible for any fitting algorithm to avoid a local minimum with very inaccurate parameter values. The main objective is therefore to elucidate the behavior of α and formulate a version of de Rooij's (2021) SWRC that has h_d as a fitting parameter but avoids the difficulties caused by the behavior of shape parameter α.

In the testing phase, it was found that a commonly used convergence criterion used in parameter optimization did not necessarily give the best parameter values if the objective function was challenging. The second objective is therefore to present a parameter fitting algorithm that employs multiple convergence criteria, and optionally explores the parameter space prior to the fitting operation to reduce the search area during fitting. The corresponding open source code for fitting the improved SWRC is provided.

2 Theory

De Rooij et al. (2021) introduced a unimodal model for the SWRC by combining those of Rossi and Nimmo (1994) and Ippisch et al. (2006), dubbed “RIA”.

$$\begin{array}{}\text{(1)}& \mathit{\theta }\left(h\right)=\left\{\begin{array}{ll}\mathrm{0},& h\le h_{\mathrm{d}}\\ {\mathit{\theta }}_{\mathrm{s}}\mathit{\beta }\ln \left(\frac{h_{\mathrm{d}}}{h}\right),& h_{\mathrm{d}}<h\le h_{\text{j}}\\ {\mathit{\theta }}_{\mathrm{s}}{\left(\frac{\mathrm{1}+{\left|\mathit{\alpha }h\right|}^n}{\mathrm{1}+{\left|\mathit{\alpha }{h}_{\mathrm{ae}}\right|}^n}\right)}^{\frac{\mathrm{1}}{n}-\mathrm{1}},& h_{\mathrm{j}}<h\le h_{\mathrm{ae}}\\ {\mathit{\theta }}_{\mathrm{s}},& h>h_{\mathrm{ae}}\end{array}.\right.\end{array}$$

Here h denotes the matric potential in equivalent water column (L), subscripts “d” and “ae” denote the value at which the water content reaches zero and the air-entry value, respectively, and subscript j indicates the value of h at which the logarithmic and sigmoid branch are joined. The volumetric water content is denoted by θ, with the subscript “s” denoting its value at saturation. Parameters α (L⁻¹) and n determine the shape of the sigmoid branch (van Genuchten, 1980), while parameter β does so for the logarithmic branch. By assuming that for h≤h_j, all water is adsorbed and for h>h_j, the adsorbed water content is equal to θ(h_j), the total water content can be partitioned in a capillary water content θ_c and an adsorbed water content θ_a.

By requiring the derivatives of the sigmoidal and logarithmic branches to match at h_j, parameter β can be expressed in terms of the other parameters (de Rooij et al., 2021).

$$\begin{array}{}\text{(2)}& \mathit{\beta }=\left(n-\mathrm{1}\right){\left|\mathit{\alpha }{h}_{\text{j}}\right|}^n{\left(\mathrm{1}+{\left|\mathit{\alpha }{h}_{\mathrm{ae}}\right|}^n\right)}^{\mathrm{1}-\frac{\mathrm{1}}{n}}{\left(\mathrm{1}+{\left|\mathit{\alpha }{h}_{\text{j}}\right|}^n\right)}^{\frac{\mathrm{1}}{n}-\mathrm{2}}.\end{array}$$

Using this expression to eliminate β from the equality that arises when the values of both branches are matched at h_j, the resulting expression can be solved for α to establish h_d as one of the fitting parameters. The expression can also be found by rearranging Eq. (9) of de Rooij et al. (2021).

$$\begin{array}{}\text{(3)}& \mathit{\alpha }={\left|h_{\text{j}}\right|}^{-\mathrm{1}}{\left[\left(n-\mathrm{1}\right)\ln \left({\displaystyle \frac{h_{\mathrm{d}}}{h_{\text{j}}}}\right)-\mathrm{1}\right]}^{-\frac{\mathrm{1}}{n}}.\end{array}$$

The five fitting parameters are h_ae, h_j, h_d, θ_s, and n. Equation (3) is only valid if the bracketed term is positive. This is the case if the following criterion is met:

$$\begin{array}{}\text{(4)}& \left|h_{\text{j}}\right|<\left|h_{\mathrm{d}}\right|e^{\frac{\mathrm{1}}{\mathrm{1}-n}}.\end{array}$$

Many soils for which SWRCs of van Genuchten (1980) or de Rooij et al. (2021) are fitted have values for α between roughly 0.001 and 0.3, with sandy soils generally having higher values than fine-textured soils (e.g., de Rooij et al., 2021). When h_d is fixed at −10^6.8 cm H₂O, α is a function of n and h_j only. Its contour map is depicted in Fig. 1. De Rooij et al. (2021) reported four soils with data sets without suspect data points above pF 3 for soils with n>1.2 (soils 1142, 1143, 2110, and 2126, all sands or loamy sands). De Rooij et al.'s (2021) values of α and n for these soils give values of h_j that are all larger (closer to zero) than −150 cm, which is unrealistic. For soil 2126 the value even exceeds h_ae, which is not physically acceptable.

Figure 1 The logarithm of shape factor α (cm⁻¹) as a function of shape factor n and the matric potential at the junction point h_j (cm H₂O) according to Eq. (3), with h_d fixed at −10^6.8 cm. The labels of the contour lines represent log(α). The transparent yellow curve is the limit of the valid domain according to Eq. (4). The black area is the invalid part of the domain.

It therefore appears from Fig. 1 that plausible combinations of α and h_j are not feasible, but Eq. (3) reveals that the relationship α(h_j) is non-monotonous. A combination of large values of both α and h_j is possible in a band too narrow to be visible in Fig. 1, located immediately below the maximum allowed value of h_j, marked by the transparent yellow curve in Fig. 1. In that band, α goes to infinity when h_j approximates its limiting value defined in Eq. (4) (Fig. 2). The partial derivative of Eq. (3) with respect to h_j is as follows:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\displaystyle \frac{\partial \mathit{\alpha }}{\partial h_{\text{j}}}}& ={\displaystyle \frac{\mathrm{1}}{h_{\text{j}}^{\mathrm{2}}}}{\left[\left(n-\mathrm{1}\right)\ln \left({\displaystyle \frac{h_{\mathrm{d}}}{h_{\text{j}}}}\right)-\mathrm{1}\right]}^{-\frac{\mathrm{1}}{n}}\\ & +{\displaystyle \frac{\mathrm{1}-n}{n}}{\displaystyle \frac{\mathrm{1}}{h_{\text{j}}^{\mathrm{2}}}}{\left[\left(n-\mathrm{1}\right)\ln \left({\displaystyle \frac{h_{\mathrm{d}}}{h_{\text{j}}}}\right)-\mathrm{1}\right]}^{-\frac{n+\mathrm{1}}{n}}.\end{array}\end{array}$$

The value of α is at its minimum where its derivative is zero. According to Eq. (5), this condition is met when the following equality holds:

$$\begin{array}{}\text{(6)}& \ln \left({\displaystyle \frac{h_{\mathrm{d}}}{h_{\text{j}}}}\right)={\displaystyle \frac{\mathrm{2}n-\mathrm{1}}{n^{\mathrm{2}}-n}}.\end{array}$$

Because Fig. 1 shows that realistic values of α require excessively low values of h_j in much of the parameter space, it may be better to use Eq. (6) to set a lower limit on the permissible values of h_j as follows:

$$\begin{array}{}\text{(7)}& \left|h_{\text{j}}\right|>\left|h_{\mathrm{d}}\right|e^{\frac{\mathrm{2}n-\mathrm{1}}{n-n^{\mathrm{2}}}}.\end{array}$$

Figure 3 shows the band of valid values of h_j enveloped by the limits set by Eqs. (4) and (7). A large part of the parameter space between these limits of h_j has an excessively large $\partial \mathit{\alpha }/\partial h_{\text{j}}$. Finding the optimum will therefore be very difficult for any parameter fitting algorithm. But if the lower limit is not enforced, trial fits showed that the shuffled complex evolution algorithm (SCE, Duan et al., 1992, 1993) consistently ended up in the region of Fig. 1 corresponding to the area below the lower limit in Fig. 3.

Figure 2 A transect of Fig. 1 for n=1.4 that shows the sharp increase of shape factor α (cm⁻¹) as the matric potential at the junction point h_j (cm H₂O) reaches its physical limit defined in Eq. (4).

Figure 3 The limits imposed on the matric potential at the junction point h_j (scaled by the matric potential at oven dryness h_d) by the requirement that shape factor α be positive (upper limit) and by the minimum value of α for any specific value of n (lower limit). The shading indicates the area with plausible parameter combinations. The dots represent de Rooij et al.'s (2021) fits of Eq. (1) with α instead of h_d as a fitting parameter for 21 soils.

De Rooij et al. (2021) already fitted Eq. (1) with α instead of h_d as a fitting parameter and without any restrictions other than minimum and maximum values imposed on any of the fitting parameters. By calculating h_d from the fitted values of α, n, and h_j according to their Eq. (10), it was possible to see if their fits fell within the limits defined above. Figure 3 shows that the values of all 21 soils fell on the upper limit. This opens the possibility to eliminate h_j as a fitting parameter and replace it by its upper limit. This leads to the following additional equation augmenting Eq. (1):

$$\begin{array}{}\text{(8)}& h_{\text{j}}=h_{\mathrm{d}}{e}^{\frac{\mathrm{1}}{\mathrm{1}-n}}.\end{array}$$

Combining Eq. (8) with Eq. (3) leads to an infinite α, consistent with Fig. 2. But Fig. 2 also shows that a minute change in h_j (much smaller than a realistic number of significant digits would be able to represent) allows α to vary beyond the range of values reported in the literature. In other words, if, for given values of h_d and n, h_j is determined from Eq. (8), α can vary over its entire range. It is therefore better to treat α as a fitting parameter and h_j as a derived parameter, by replacing Eq. (3) by Eq. (8). This has the added advantage that the entire parameter space defined by the minimum and maximum values of the fitting parameters is valid, provided the physical and mathematical limits of each parameter are respected, and the fitted or fixed value of h_d is smaller (more negative) than h_ae.

In the limit as α approximates infinity, the expression for the sigmoid branch of the SWRC simplifies to the power law proposed by Brooks and Corey (1964).

$$\begin{array}{}\text{(9)}& \begin{array}{c}\text{lim}\\ \mathit{\alpha }\to \mathrm{\infty }\end{array}{\mathit{\theta }}_{\mathrm{s}}{\left({\displaystyle \frac{\mathrm{1}+{\left|\mathit{\alpha }h\right|}^n}{\mathrm{1}+{\left|\mathit{\alpha }{h}_{\mathrm{ae}}\right|}^n}}\right)}^{\frac{\mathrm{1}}{n}-\mathrm{1}}={\mathit{\theta }}_{\mathrm{s}}{\left({\displaystyle \frac{h}{h_{\mathrm{ae}}}}\right)}^{\mathrm{1}-n}.\end{array}$$

It can be easily shown that both the values and the derivatives of the dry and the wet branch match at h_j if Eq. (9) is used for the latter. Equation (9) establishes that Rossi and Nimmo's (1994) junction model (without the parabolic smoothing near saturation that Madi et al. (2018) showed to be detrimental to the hydraulic conductivity function) is a special case of the RIA parameterization. Incidentally, this implies that Brooks and Corey's model (1964) is a special case of that of Ippisch et al. (2006).

If Eq. (8) replaces Eq. (3), the fitted value of α no longer ensures continuity at the junction point. The ensuing continuity gap can be closed by applying a correction factor c to the value of h_d used in the logarithmic branch as follows:

$$\begin{array}{}\text{(10)}& \mathit{\theta }={\mathit{\theta }}_{\mathrm{s}}\mathit{\beta }\ln \left[{\displaystyle \frac{\left(\mathrm{1}+c\right)h_{\mathrm{d}}}{h}}\right].\end{array}$$

The correction factor c is defined by the following expression, which is found by replacing h_d in Eq. (1) by (1+c)h_d, requiring the logarithmic and sigmoid branches of Eq. (1) to be equal at h_j, and replacing the ratio $h_{\mathrm{d}}/h_{\text{j}}$ in the resulting equality by exp[1/(n–1)] according to Eq. (8).

$$\begin{array}{}\text{(11)}& c=\exp \left[{\displaystyle \frac{\mathrm{1}}{\mathit{\beta }}}{\left({\displaystyle \frac{\mathrm{1}+{\left|\mathit{\alpha }{h}_{\text{j}}\right|}^n}{\mathrm{1}+{\left|\mathit{\alpha }{h}_{\mathrm{ae}}\right|}^n}}\right)}^{\frac{\mathrm{1}}{n}-\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-n}}\right]-\mathrm{1}.\end{array}$$

In the expression for β, h_d appears indirectly. Its value should not be corrected there because β ensures continuity of the derivatives only if its expression is not modified. Trial calculations showed that c is negligible, except when both α and n are small.

Rossi and Nimmo (1994) fitted h_j values for their parameterization of the SWRC between $-\mathrm{2.6}\times {\mathrm{10}}^{\mathrm{6}}$ and $-\mathrm{2.0}\times {\mathrm{10}}^{\mathrm{4}}$ cm H₂O for seven soils, with only the one clayey soil having a value more negative than $-\mathrm{1.9}\times {\mathrm{10}}^{\mathrm{5}}$ cm H₂O. Tuller and Or (2005) tentatively set the matric potential at which the capillary-bound water content becomes negligible −10⁵ cm H₂O, based on data from a sand mixture and six soils that mostly overlap with those of Rossi and Nimmo (1994) (Or and Tuller, 1999). Beyond this there is little guidance on the value of h_j in the literature. When Eq. (8) is used, such guidance is not necessary. The map of h_j as defined by Eq. (8) in Fig. 4 shows that for n below 1.4, h_j is very sensitive to the value of n. For n>2, log(−h_j) is roughly proportional to log(−h_d) for a given value of n. When h_d is close to −10^6.8 cm H₂O (Schneider and Goss, 2012), values of h_j in the range of those reported in the literature occur for n≥1.25.

Figure 4 The logarithm of the absolute value of the matric potential at the junction point (h_j, cm H₂O) as a function of the matric potential at oven dryness (h_d, cm H₂O) and shape factor n, according to Eq. (8). The labels of the contour lines denote log(−h_j). In the black region, $h_{\text{j}}>-\mathrm{100}$ cm.

Table 1 Fitted parameter values and the root mean square error (RMSE) of the best fits for 21 soils. The corresponding values of the derived parameter h_j are given as well. If parameter h_d was fixed during the fitting operation, its value is denoted in italic font.

For completeness, the multimodal version of Eq. (1) is provided as well. The multimodality is limited to the sigmoid branch, so that the multimodal SWRC has only one value for h_j. Because Eq. (8) allows only a single value of n in that case, only α can be varied between the constituting sigmoid curve sections.

$$\begin{array}{}\text{(12)}& \mathit{\theta }\left(h\right)=\left\{\begin{array}{ll}\mathrm{0},& h\le h_{\mathrm{d}}\\ {\mathit{\theta }}_{\mathrm{s}}{\mathit{\beta }}_{\mathrm{m}}\ln \left(\frac{h_{\mathrm{d}}}{h}\right),& h_{\mathrm{d}}<h\le h_{\text{j}}\\ {\mathit{\theta }}_{\mathrm{s}}\sum _{i=\mathrm{1}}^k{w}_i{\left(\frac{F_i\left(h\right)}{F_i\left(h_{\mathrm{ae}}\right)}\right)}^{\frac{\mathrm{1}}{n}-\mathrm{1}},& h_{\mathrm{j}}<h\le h_{\mathrm{ae}}\\ {\mathit{\theta }}_{\mathrm{s}},& h>h_{\mathrm{ae}}\end{array}.\right.\end{array}$$

For brevity, the following function was introduced in Eq. (12):

$$\begin{array}{}\text{(13)}& F_i\left(x\right)=\mathrm{1}+{\left|{\mathit{\alpha }}_{i}x\right|}^n.\end{array}$$

Here, k denotes the modality, w_i is the weighting factor (adding up to one) of the ith constituting curve, and α_i (L⁻¹) its shape factor. The expression for the multimodal β_m is found be setting the derivatives of the logarithmic and the sigmoid branch equal at h_j and invoking Eq. (8).

$$\begin{array}{}\text{(14)}& {\mathit{\beta }}_{\mathrm{m}}={\displaystyle \frac{n-\mathrm{1}}{e}}{\left|h_{\mathrm{d}}\right|}^{n-\mathrm{1}}\sum _{i=\mathrm{1}}^k{w}_i{\mathit{\alpha }}_i^n{F}_i\left(h_{\mathrm{ae}}{)}^{\mathrm{1}-\frac{\mathrm{1}}{n}}{G}_i\right(h_{\mathrm{d}}{)}^{\frac{\mathrm{1}}{n}-\mathrm{2}}.\end{array}$$

Function G is defined as follows:

$$\begin{array}{}\text{(15)}& G_i\left(x\right)=\mathrm{1}+{\left|{\mathit{\alpha }}_{i}x\right|}^n{e}^{\frac{n}{\mathrm{1}-n}}.\end{array}$$

The continuity correction factor c_m can be found by requiring that the logarithmic and sigmoid branch join at h_j.

$$\begin{array}{}\text{(16)}& c_{\mathrm{m}}=\exp \left[{\displaystyle \frac{\mathrm{1}}{{\mathit{\beta }}_{\mathrm{m}}}}\sum _{i=\mathrm{1}}^k{w}_i{\left({\displaystyle \frac{G_i\left(h_{\mathrm{d}}\right)}{F_i\left(h_{\mathrm{ae}}\right)}}\right)}^{\frac{\mathrm{1}}{n}-\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{1}-n}}\right]-\mathrm{1}.\end{array}$$

Above, the subscript “m” is used to distinguish multimodal versions of β and c from their single-mode equivalents.

3 Materials and methods

3.1 Selected soils

The SWRC has five fitting parameters: θ_s, h_ae, h_d, α, and n. These were fitted to data of the 21 soils selected from the UNSODA database (Nemes et al., 2001; National Agricultural Library, 2015) by de Rooij et al. (2021), with h_d either fixed at −10^6.8 cm H₂O, or its lower limit set at that value. These soils cover a wide range of textures (see Madi et al., 2018 and de Rooij et al., 2021 for details). At the time the database was created, the pressure plate apparatus was widely used in the dry range. Therefore, the standard deviation (SD) of the matric potential of any data point with $h\le -\mathrm{1000}$ cm H₂O was set to half its value, thereby drastically reducing its weight in the objective function. The SD of h for $\mathrm{0}>h>-\mathrm{1000}$ cm H₂O was set to 1.0 cm, and that for h=0 at 0.05 cm. The SD for the volumetric water content was set at 0.01 when h=0, and to 0.02 otherwise. The sample height was set to zero for $h\le -\mathrm{1000}$ cm because the variation in the water content can be neglected at low matric potentials. For h=0 cm, it was assumed that the water pressure in the sample was positive everywhere, or that the porosity was measured. In both cases, there is no effect of the matric potential on the local water content, hence the sample height was set to zero to reflect this. For the intermediate values of h, the sample height was set to the value specified in the UNSODA database. If a sample height was not reported, it was set to 3.0 cm.

3.2 Parameter fitting

The parameters were fitted using the SCE algorithm (Duan et al., 1992, 1993), implemented in Fortran in a code that accompanies this paper. The most important features of the code are summarized here. The code itself, further details of the code and the algorithm, and a user manual can be downloaded (de Rooij, 2022). For each case, the code performs three optimization runs by minimizing the objective function: the root mean square error (RMSE) of the differences between fitted and observed volumetric water contents, weighted according to the error standard deviations of the observed matric potentials and corresponding water contents provided on input, as detailed by de Rooij (2022).

Ten convergence criteria are evaluated. Criteria 1 and 2 take into account the results of the last few shuffles. The number of shuffles considered is twice the number of fitting parameters or an internally set number (5), whichever is larger.

1. In the best fits from the last set of shuffles, the range of a parameter exceeds neither the absolute nor the relative user-specified tolerance.

2. In the best fits from the last set of shuffles, the range of the objective function does not exceed its absolute user-specified tolerance.

3. The parameter range in the final complexes does not exceed the maximum internally set permissible value.

4. The volume of the hypercube enveloping the final complexes does not exceed the maximum internally set permissible value.

5. The parameter range in the most successful complex (minus the point with the highest RMSE) does not exceed the internally set maximum permissible value.

6. The volume of the hypercube enveloping the most successful complex (minus the point with the highest RMSE) does not exceed the internally set maximum permissible value.

7. A parameter does not exceed both the absolute and the relative user-specified tolerance in the final complexes.

8. A parameter does not exceed both the absolute and the relative user-specified tolerance in the most successful complex (without the point with the highest RMSE).

9. The change of the objective function between consecutive shuffles does not exceed the user-specified tolerance.

10. The root mean square error of the fit does not exceed a user-specified tolerance.

A relative improvement criterion similar to criterion 9 is often used as the sole criterion (e.g., in the R-package SoilHyp, Dettmann, 2021). Criteria 1, 3, 5, 7, and 8 are evaluated separately for each parameter. Convergence is achieved when no more than a user-prescribed number of these criteria failed for any of the parameters. The code keeps evolving and shuffling complexes until convergence is achieved or the user-specified maximum allowed number of evaluations of the objective function is exceeded. If not all criteria are considered relevant, the user can either set their thresholds unrealistically strict and increase the number of criteria that are allowed to fail, or set them excessively loose and decrease the number of criteria that are allowed to fail accordingly.

Table 2 The continuity correction factor c (Eq. 11) and the corresponding shift on the pF scale of the dry-branch correction for 21 soils.

Figure 5 Soil water retention data and fitted curves for soils of classes A1 and A2 of Twarakavi et al. (2010). Curves fitted with h_d fixed at −10^6.8 cm H₂O, and with h_d fitted with a cap at that value are shown. The one with the lowest root mean square error is shown as a red solid line. The volume fractions of capillary-bound water and water adsorbed in films is shown for this curve. The other curve is shown as black dash–dot line. This curve has not been corrected for continuity at the junction point. The vertical axis denotes the logarithm (base 10) of the absolute value of the matric potential in cm H₂O.

Figure 6 As Fig. 5, for soils of classes A3 and A4 of Twarakavi et al. (2010).

The algorithm generates large numbers of sets of fitted parameter values. A random sample of these is used to determine the correlation matrix of the parameters. The best fit, its RMSE, and its correlation matrix are reported by the code for each of the three runs, and the run with the overall lowest RMSE is identified. The code returns a table of the fitted curve based on the best run, and reports the correction factor c used to compile this table. These tables are the basis for the plots shown below. If desired, the code also calculates the objective function on points of a regular grid covering the parameter space (map points) and writes a random sample of these to output. Even if this is not desired, a map is calculated based on a three point grid along each principal axis of the parameter space. This resulting output is helpful if a user wishes to verify if the objective function is correctly calculated.

Normally, the first complexes of each run are filled with randomly selected points in the valid regions of the parameter space. Optionally, only the complexes of run 3 are filled with randomly selected points, while the first complexes of run 1 are filled with the map points with the lowest RMSE. The first complexes of run 2 then contain these map points perturbed by adding random noise to the parameter values.

For each fitting parameter, a maximum and minimum value need to be provided. If these values are equal, the parameter is treated as a fixed value, and the dimensionality of the parameter space is reduced accordingly. The number of complexes is two (for eight of fewer fitting parameters) or four (see Duan et al., 1994). If this leads to frequent convergence at local minima, the number of complexes can be set to twice the number of fitting parameters. The number of individuals in a complex and the number of evolution steps are twice the number of fitting parameters plus one. The number of individuals in the subcomplexes is the number of fitting parameter plus one. The number of offspring in each evolution step is one. These settings are all in accordance with Duan et al. (1994).

When the user chooses to use the map of regularly spaced points in the parameter space to set the initial guesses of the first two runs, the code adapts the parameter ranges based on their ranges among the map points selected to fill the initial complexes.

The permitted parameter range for θ_s enveloped the range of observed saturated water contents, with a limited buffer on either side of the range. The range of h_ae was determined based on the wettest unsaturated data point and the driest saturated point, with a generous buffer. The range limit at the wet end was often set to zero. The value of h_d was mostly fixed and occasionally varied over a relatively wide range depending on visual examination of the data points. The permitted range for α was 0.001 to 0.5, except for 1142, where alpha could go as high as 100.0. After some trial and error, the range of n was set relatively narrow (between 1.05 and 2) because even with wider ranges the fitted values fell within this range. Any time a fitted value was close to one of its limits, the fit was repeated with an expanded range.

No more than four convergence criteria were allowed to fail for convergence to be achieved. If convergence was not achieved, up to 20 000 evaluations of the objective function evaluations would be performed. The actual number could be slightly higher because it was checked each time a shuffle had been completed. A map was not created, and therefore all three optimization runs started with random parameter combinations filling the complexes. The maximum allowed relative change in the RMSE between consecutive shuffles was 10⁻⁶. The maximum allowed value of the RMSE was 0.1. For the parameters, the absolute tolerances for θ_s, h_ae, h_d, α, and n were 0.001, 0.1, 1000.0, 0.1, and 0.01, respectively. The relative tolerances were 0.1 for α and 0.01 for the others. This choice reflects the limited sensitivity to α.

The internally set relative tolerances for parameter variations for all complexes combined and for the most successful complex were both 0.01. The internally set required maximum size of the hypercube enveloping the range of fitted parameter values (again for all complexes combined and the most successful one), scaled by the volume of the hypercube defined by the minimum and maximum allowed parameter values, equals 0.01^d, with d the number of fitting parameters that were not set at fixed values by the user.

Figure 7 As Fig. 5, for soils of classes B2 and B4 of Twarakavi et al. (2010).

Figure 8 As Fig. 5, for soils of classes C2 and C4 of Twarakavi et al. (2010).

4 Results and discussion

In the tables that follow, the soil texture classification according to Twarakavi et al. (2010) is provided. The advantage of this classification over the conventional USDA soil texture classification is that it better recognizes differences in hydraulic behavior between its textural classes. This makes it more relevant for soil water flow studies for which the SWRCs are typically used.

Table 1 shows the fitted parameter values of the best fits, the resulting RMSE, and the corresponding value of h_j. The value of α for soil 1142 stands out. As explained above, high values of α lead to a power-law type of SWRC between h_ae and h_j. An additional fit with α capped at 1.0 confirmed that the shape of the SWRC is not very sensitive to large values of α. The increase in the RMSE caused by limiting the range of α was 10⁻⁴, and the changes in θ_s, h_ae, and n did not occur before the fourth significant digit. One can also switch from a sigmoidal model to a power-law model by invoking Eq. (9) if α is very large, keeping the fitted values of the other parameters.

Fuentes et al. (1991) showed that for values of n smaller than two, the asymptotic dry branch of the original parameterization of van Genuchten (1980) would lead to physically unacceptable behavior. All soils in Table 1 have values of n in this range. This highlights the importance of avoiding a dry branch with an asymptote at a residual water content.

The range of values of h_j in Table 1 is only slightly beyond the range reported by Rossi and Nimmo (1994) for a smaller number of soils. This lends credibility to Eq. (8). In the few cases where h_d was fitted, the resulting values in Table 1 are close to the value proposed by Schneider and Goss (2012). Table 2 shows the correction factor c of Eqs. (10) and (11), which ensures continuity of the SWRC. Seven of the 21 soils need a correction of h_d that exceeds 1 %. The resulting shift of the dry-branch pF is also shown. For most soils, the shift is negligible. Only for soils 1122 and 1123 (both fine-textured soils with small values for both α and n), the shift exceeds 0.1 pF unit, but never more than 0.2 unit.

Only the optimizations for soils 1142 and 2104 converged, with convergence criteria 4, 6, 8, 9, and 10 satisfied for all parameters for soil 1142, and criteria 4 through 10 for soil 2104. None of the correlation coefficients of the parameter pairs for either soil exceeded 0.31. The other optimizations ran until the maximum number of objective function evaluations was exceeded. For soil 1120, criteria 9 and 10 were met for all parameters. For the remaining soils, criteria 1, 2, and 9 were satisfied in all cases. For 14 soils, criterion 10 was met as well. For soil 3250, criterion 8 was also satisfied. The lack of convergence forced the code to keep exploring the parameter space, leading to a large proportion of randomly selected points because the reflection and contraction points determined by the SCE algorithm did not improve the fit. If the majority of points is randomly selected, there is no correlation between the parameters, and the correlation matrix does not provide any information.

In all cases, the fitted parameter values for the runs with h_d fixed and h_d free and the three individual runs for each optimization were essentially in agreement, and the parameter values did not look suspect. Therefore, it was considered unnecessary to run the optimizations with modified convergence requirements in order to obtain more meaningful correlation matrices.

The reduced weights assigned to data points with pF>3 are reflected in the plots of Figs. 5 through 8, which show the fitted curves with h_d fixed and h_d fitted. In these plots, the fit with the lowest RMSE is plotted in red, and the corresponding curves for θ_c and θ_a are included. The other curve is shown in black. To illustrate how small the continuity correction c is, this curve is shown without this correction. The discontinuity at the junction point is only visible for soils 1133 (Fig. 6) and 1122 (Fig. 8).

The fraction of adsorbed water increases when moving from sands (Figs. 5 and 6) through loams (Figs. 6 and 7) to clays (Fig. 8). Because the separation between capillary and adsorbed water is abrupt and binary at h_j, this should not be interpreted as representative for the more smooth transition in natural soils. Nevertheless, the direction of the trend is physically plausible.

Most soils (2126 and 1142 in Fig. 5; 1120, 1121, 1143, 2110, and 2132 in Fig. 6; 1142 and 2126 in Fig. 7; 1122, 1180, and 1182 in Fig. 8) have observed saturated water contents that seem to be too large compared to the other data points. The causes (e.g., macropores or air inclusion) are not known. Data points at saturation were assumed to be very accurate and therefore had a high weight, which the plots reflect. It sometimes results in relatively low (more negative) air-entry values in coarse soils (Figs. 5 and 6 and Table 1, most notably soil 1142).

All B2 soils (silt loams) and two out of three B4 soils (both silty clay loams) have high values of h_ae, indicating that the maximum pore size is large (Table 1). Although one would suspect that such fine-textured soils would have a low (more negative) air-entry value, the results are consistent with the data, as Fig. 7 shows. Some of the C2 and C4 soils (1180–1182) have high RMSE values (Table 1). Their plots in Fig. 8 reveal that the multimodal shape of the curves was not captured well by Eq. (1). The remaining soils in Fig. 8 had very few points in the dry range, and fixing h_d was very effective in guiding the dry branch of the SWRC.

This paper focuses on mineral soils with unimodal SWRCs. If an extension to multimodality is desired, Eqs. (12)–(16) provide a starting point. Further testing on soils with specific characteristics, such as volcanic or organic soils, may be worthwhile. In the case of organic soils, the risk of irreversible changes will require some caution when measuring points on the SWRC.