Fitting the junction model and other parameterizations for the unsaturated soil hydraulic conductivity curve: KRIAfitter version 1.0

Several current models for the unsaturated soil hydraulic conductivity curve consider the conductivity of the domains of capillary water in water-filled pores and adsorbed water in films on soil grains, as well as an equivalent conductivity for water vapour diffusion. These models rely on unrealistic configuration of the domains. A junction model is introduced that sidesteps this problem by assigning all liquid water to films (dry range) or to capillaries (wet range). Combined with a sigmoidal junction model for the soil water retention curve, it has up to six fitting parameters, one fewer than the other multi-domain models. Tests on data for 13 soils show that the junction model and an additive model (that adds all domain conductivities) often produce good fits. Models with six or more parameters may be overparameterized for many soils, giving the more parsimonious junction model an advantage, but for some soils, the extra parameter of the additive model is needed to achieve a good fit. This paper and a user manual document a Fortran code (KRIAfitter) that uses the shuffled complex evolution algorithm to fit the junction, additive, and four other conductivity models for any combination of fixed and fitting parameters or their log transforms. KRIAfitter either maps the root mean square error in the entire parameter space in order to then constrain the parameter space around the likely global minimum or generates many fits and uses those to calculate statistics for individual parameters, as well as the covariance and correlation matrices.

1 Introduction

Madi et al. (2018) revealed that many parameterizations of the soil water retention curve (SWRC) gave physically unrealistic near-saturated behaviour of the soil hydraulic conductivity when applied in combination with Kosugi's (1999) generalized soil hydraulic conductivity parameterization. This motivated de Rooij et al. (2021) and de Rooij (2022) to propose a closed-form expression for the SWRC with a distinct air-entry value, such as in 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 at which the soil was oven dry, with the water content essentially zero. The volumetric water content and derivatives of the sigmoid and logarithmic branches were matched at the matric potential of their junction according to Rossi and Nimmo (1994). This SWRC model (termed RIA, for Rossi–Ippisch adaptation) had a finite slope near saturation, which remedied the problems of most existing parameterizations and also eliminated the asymptotic residual water content at the dry end.

De Rooij et al. (2021) presented an analytical expression for a soil hydraulic conductivity curve that could be used with this new SWRC parameterization. This conductivity curve was a special case of Kosugi's (1999) parameterization. The public discussion of de Rooij et al.'s paper (accessible on-line) revealed the desirability of an alternative formulation that could be fitted separately and had the capability to account for non-capillary flow.

Peters and Durner (2008) and Peters (2013) included non-capillary flow in their model for the unsaturated hydraulic conductivity curve (UHCC) by separating the total liquid water in a domain with adsorbed water (present as films on the surface of the solid phase) and a domain with capillary bound water. They combined parametric models for soil water retention and unsaturated hydraulic conductivity for both. By assuming instantaneous equilibrium between water vapour pressure and the matric potential of liquid water, Peters (2013) could also formulate a model for isothermal vapour flow driven by the gradient in the matric potential. Weber et al. (2019) formalized this approach in a modular setup facilitating different choices for the parameterizations chosen to represent the SWRC and the UHCC of the various domains. In all three papers, the adsorbed and capillary-bound water content were added to arrive at the total water content. The hydraulic conductivities were also added to find the total conductivity at a given matric potential. De Rooij (2024a) posited that this additivity attribute requires that all flow domains are arranged in parallel and developed averaging models using the arithmetic and harmonic means of the domain conductivities for domains in parallel and in series, respectively. He also offered a model using the geometric mean as an intermediate between the other two as well as a non-weighted additive model akin to that of Peters (2013). From his analysis it appears fundamentally impossible to derive the bulk soil hydraulic conductivity from domain conductivities based on domain volumes and configurations.

Given these complications, this paper introduces a new model in which the film domain and the capillary domain do not exist in parallel but instead are joined at a critical matric potential below which all liquid water is in films and above which all water is capillary-bound. This creates a junction model for the soil hydraulic conductivity analogously to the RIA parameterization for the SWRC. The vapour domain is assumed to be parallel to the liquid-water domain. Together with this new model, the paper documents a Fortran program (KRIAfitter) that is able to fit the parameters of the junction model and all models introduced by de Rooij (2024a) to unsaturated hydraulic conductivities observed for either different water content or different matric potentials (de Rooij, 2024b). All models operate in conjunction with the RIA parameterization of de Rooij (2022). The fitting code (RIAfitter) for that model was thoroughly overhauled for this study, which resulted in version 2.0 (de Rooij, 2024c). While KRIAfitter is the main focus of this paper, it is expected that RIAfitter and KRIAfitter will normally be used in tandem.

2 The junction models for the soil water retention and hydraulic conductivity curves

2.1 The model for the soil water retention curve

This section presents the main equations from de Rooij (2022) for clarity. He presented a unimodal model for the SWRC.

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

where 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 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{(1b)}& \begin{array}{rl}{\mathit{\beta }}_{\mathit{\theta }}=& (n_{\mathit{\theta }}-\mathrm{1}){\left|{\mathit{\alpha }}_{\mathit{\theta }}{h}_{\mathrm{j}\mathit{\theta }}\right|}^{n_{\mathit{\theta }}}{\left(\mathrm{1}+{\left|{\mathit{\alpha }}_{\mathit{\theta }}{h}_{\text{ae}\mathit{\theta }}\right|}^{n_{\mathit{\theta }}}\right)}^{\mathrm{1}-\frac{\mathrm{1}}{n_{\mathit{\theta }}}}\\ & \times {\left(\mathrm{1}+{\left|{\mathit{\alpha }}_{\mathit{\theta }}{h}_{\mathrm{j}\mathit{\theta }}\right|}^{n_{\mathit{\theta }}}\right)}^{\frac{\mathrm{1}}{n_{\mathit{\theta }}}-\mathrm{2}}\end{array}\end{array}$$

De Rooij (2022) made h_jθ a derived parameter with the following expression:

$$\begin{array}{}\text{(1c)}& h_{\mathrm{j}\mathit{\theta }}=h_{\mathrm{d}}\text{exp}\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}-n_{\mathit{\theta }}}}\right).\end{array}$$

The five fitting parameters then are h_aeθ, h_dθ, θ_s, α_θ, and n_θ. All of them appear in the expressions for the UHCC as well. Those that can have different values for the UHCC have θ in their subscripts to avoid ambiguity in the remainder of the text. Equation (1c) does not guarantee continuity at h_jθ, so de Rooij (2022, 2024a) introduced a correction factor that is very small for most soils.

$$\begin{array}{}\text{(1d)}& c=\text{exp}\left\{{\displaystyle \frac{\mathrm{1}}{(n_{\mathit{\theta }}-\mathrm{1}){\left[\left|{\mathit{\alpha }}_{\mathit{\theta }}{h}_{\mathrm{d}\mathit{\theta }}\right|\text{exp}\left(\frac{\mathrm{1}}{\mathrm{1}-n_{\mathit{\theta }}}\right)\right]}^{n_{\mathit{\theta }}}}}\right\}-\mathrm{1}\end{array}$$

This correction needs to be applied to h_dθ in the logarithmic branch of the SWRC as shown in Eq. (1a). This correction needs to be carried over to the expression for the soil hydraulic conductivity to liquid water in the UHCC.

2.2 The junction model for the unsaturated hydraulic conductivity curve

A junction model analogous to that used for the SWRC above gives a monotonically increasing UHCC with increasing matric potential. It requires that liquid soil water is entirely allocated to either capillary-bound water or water adsorbed in films, with the water abruptly changing its allocation when the matric potential passes through the matric potential at the junction point h_j (L), which will not necessarily be equal to h_jθ. Because the liquid phase in this model only occupies a single domain at any given matric potential, the issue of the correct averaging of the domain conductivities with its complications (de Rooij, 2024a) becomes moot. The fact that the conductivities (L T⁻¹) for the film domain (K_a) and the capillary domain (K_c) are intrinsic conductivities strictly valid for their respective domains instead of bulk conductivities (de Rooij, 2024a) is resolved implicitly through the fitted values of the parameters, made possible because averaging of domain conductivities is not needed.

The clear distinction between adsorbed and capillary-bound water of Eqs. (1a–d) appears well suited for use with the multi-domain conceptualization of Peters and Durner (2008) and Peters (2013) and with the modular framework presented by Weber et al. (2019). Weber et al.'s (2019) approach is tailored to additive formulations of the SWRC in which capillary-bound and adsorbed water co-exist over the full moisture range, but the SWRC of Eqs. (1a–d) has no capillary-bound water for h<h_jθ. Furthermore, a water film on a soil particle is bounded on one side by a solid–liquid interface and on the other by a liquid–gas interface. The velocity profile in the water film is such that its gradient at the liquid–gas interface equals zero (Eq. 8) of Or and Tuller, 2000). When a pore with a water film on its grain surfaces takes in additional water and becomes fully saturated, the region previously occupied by the water film only has a solid–liquid interface, while at the location of the former liquid–gas interface, there is moving water and a non-zero gradient in the velocity profile. This will increase the flow rate in the region previously occupied by the film. Simple addition of film and capillary conductivities may therefore not be accurate. De Rooij (2024a) showed that it is fundamentally impossible to develop a physically based model for the UHCC based on domain conductivities. The junction model introduced here does not claim to have more solid physical bases than other multi-domain UHCC models, but it is simpler than the existing models. To implement the junction model, the conductivity expressions were adapted as outlined below.

The intrinsic hydraulic conductivity of water in films is modelled according to Peters (2013).

$$\begin{array}{}\text{(2a)}& K_{\mathrm{a}}\left(h\right)=K_{\mathrm{s},\mathrm{a}}\cdot \left\{\begin{array}{ll}\mathrm{0},& h\le h_{\mathrm{d}}\\ {\left(\frac{h}{h_{\mathrm{a}}}\right)}^{-\mathrm{1.5}},& h_{\mathrm{d}}<h\le h_{\mathrm{a}}\\ \mathrm{1},& h>h_{\mathrm{a}}\end{array}\right.\end{array}$$

K_s,a (L T⁻¹) is the value of K_a when the domain with adsorbed water is completely filled, and h_a (L) is the matric potential at which this occurs. The value of the exponent is adopted from Peters (2013). Note that K_a(h) abruptly drops to zero at h_d (the matric potential at oven dryness, L), but K_a at that matric potential is so small that this will generally be insignificant for practical use. The need to carry over the correction of h_dθ in Eq. (1a) results in the following equality:

$$\begin{array}{}\text{(2b)}& h_{\mathrm{d}}=(\mathrm{1}+c)h_{\mathrm{d}\mathit{\theta }}.\end{array}$$

In the junction model, h_a is set equal to h_j. For the range in which capillary-bound water is present, Kosugi's (1999) model is used with κ=1 (see de Rooij, 2024a) for the soil hydraulic conductivity as if all water is capillary-bound, irrespective of the matric potential.

$$\begin{array}{}\text{(3a)}& \begin{array}{rl}& K_{\mathrm{c},\text{jun}}\left(h\right)=K_{\mathrm{s},\mathrm{c},\text{jun}}\\ & \times \left\{\begin{array}{ll}{\left(\frac{G\left(h\right)}{G\left(h_{\text{ae}}\right)}\right)}^{\mathit{\tau }}{\left(\frac{\mathrm{1}-{\left|\mathit{\alpha }h\right|}^{n-\mathrm{1}}G\left(h\right)}{\mathrm{1}-{\left|\mathit{\alpha }{h}_{\text{ae}}\right|}^{n-\mathrm{1}}G\left(h_{\text{ae}}\right)}\right)}^{\mathit{\gamma }},& h_{\mathrm{j}}<h\le h_{\text{ae}},\\ \mathrm{1},& h>h_{\text{ae}},\end{array}\right.\end{array}\end{array}$$

where

$$\begin{array}{}\text{(3b)}& G\left(x\right)={\left(\mathrm{1}+{\left|\mathit{\alpha }x\right|}^n\right)}^{\frac{\mathrm{1}}{n}-\mathrm{1}}.\end{array}$$

Here, K_c,jun is the unsaturated hydraulic conductivity of the capillary domain in the junction model (L T⁻¹) and $K_{\mathrm{s},\mathrm{c},\text{jun}}$ its value at saturation. The parameters h_ae (L), α (L⁻¹), and n are the equivalent of h_aeθ, α_θ, and n_θ, respectively, in the equations for the SWRC. Mualem (1976) proposed the value of 0.5 for τ and 2.0 for γ. Assouline (2001) introduced a simpler expression by setting τ to 0.0.

Equations (1b) and (3a, b) establish a continuous conductivity curve if their values are matched at h_j. This requires that the following equality holds.

$$\begin{array}{}\text{(4)}& K_{\mathrm{s},\mathrm{a}}=K_{\mathrm{s},\mathrm{c},\text{jun}}{\left({\displaystyle \frac{G\left(h_{\mathrm{j}}\right)}{G\left(h_{\text{ae}}\right)}}\right)}^{\mathit{\tau }}{\left({\displaystyle \frac{\mathrm{1}-{\left|\mathit{\alpha }{h}_{\mathrm{j}}\right|}^{n-\mathrm{1}}G\left(h_{\mathrm{j}}\right)}{\mathrm{1}-{\left|\mathit{\alpha }{h}_{\text{ae}}\right|}^{n-\mathrm{1}}G\left(h_{\text{ae}}\right)}}\right)}^{\mathit{\gamma }}.\end{array}$$

Vapour flow is assumed to be diffusive under isothermal conditions. The equilibrium between local vapour pressure and matric potential is assumed to be instantaneous. The model used is de Rooij's (2024a) modification of that of Peters (2013). The equivalent water vapour bulk hydraulic conductivity $K_{\mathrm{v}}^{\mathrm{B}}$ is

$$\begin{array}{}\text{(5)}& K_{\mathrm{v}}^{\mathrm{B}}\left(h\right)={\displaystyle \frac{{\mathit{\rho }}_{\text{sv}}Mg{D}_{\mathrm{a}}}{{\mathit{\rho }}_{\mathrm{w}}R(T+\mathrm{273.15}){\mathit{\theta }}_{\mathrm{s}}^{\mathrm{2}}}}({\mathit{\theta }}_{\mathrm{s}}-\mathit{\theta }{)}^{\frac{\mathrm{10}}{\mathrm{3}}}{e}^{\frac{Mg}{R(T+\mathrm{273.15})}h}.\end{array}$$

Because of the specific nature of the expressions and their constants, units replace dimensions in the explanations of the variables and parameters. The unit of $K_{\mathrm{v}}^{\mathrm{B}}$ is cm d⁻¹, consistent with those for the liquid-water conductivity. The diffusion coefficient of water vapour in air at the ambient temperature T (°C) is D_a (cm² d⁻¹). The densities of saturated water vapour and liquid water are ρ_sv and ρ_w (kg m⁻³), respectively; M (kg mol⁻¹) is the molar mass of water; g (cm d⁻²) is the gravitational acceleration; and R is the universal gas constant ($\mathrm{kg}{\mathrm{cm}}^{\mathrm{2}}{\mathrm{d}}^{-\mathrm{2}}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{mol}}^{-\mathrm{1}}$).

The densities ρ_w and ρ_sv are temperature dependent. Peters (2013) provides an expression for the temperature dependence of ρ_sv:

$$\begin{array}{}\text{(6)}& {\mathit{\rho }}_{\text{sv}}={\displaystyle \frac{\mathrm{0.0010}}{T+\mathrm{273.15}}}\text{exp}\left[\begin{array}{c}\mathrm{31.3716}-\frac{\mathrm{6014.79}}{T+\mathrm{273.15}}\\ -\mathrm{0.00792495}(T+\mathrm{273.15})\end{array}\right].\end{array}$$

The density of liquid water in the range between 0 and 40 °C was approximated according to Brutsaert (2005, p. 17).

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathit{\rho }}_{\mathrm{w}}=& \mathrm{999.8505}+\mathrm{0.06001}T-\mathrm{0.007917}{T}^{\mathrm{2}}\\ & +\mathrm{4.1256}\times {\mathrm{10}}^{-\mathrm{5}}{T}^{\mathrm{3}}\end{array}\end{array}$$

For the temperature-dependency of D_a, the expression of Dorsey (1940, p. 73) was converted to units of cm and day and adapted to a reference temperature of 15 °C.

$$\begin{array}{}\text{(8)}& D_{\mathrm{a}}=\mathrm{2.09}\times {\mathrm{10}}^{\mathrm{4}}{\left({\displaystyle \frac{T+\mathrm{273.15}}{\mathrm{288.15}}}\right)}^{\mathrm{1.75}}\end{array}$$

The bulk vapour conductivity model does not have any fitting parameters.

Combining Eqs. (1b), (3a), (4), and (5) gives the expression for the bulk hydraulic conductivity to water according to the junction model, denoted as $K_{\text{jun}}^{\mathrm{B}}$ (L T⁻¹):

$$\begin{array}{}\text{(9)}& \begin{array}{rl}& K_{\text{jun}}^{\mathrm{B}}\left(h\right)=K_{\mathrm{v}}^{\mathrm{B}}\left(h\right)+K_{\mathrm{s},\mathrm{c},\text{jun}}\\ & \times \left\{\begin{array}{ll}\mathrm{0},& h\le h_{\mathrm{d}},\\ {\left(\frac{G\left(h_{\mathrm{j}}\right)}{G\left(h_{\text{ae}}\right)}\right)}^{\mathit{\tau }}{\left(\frac{\mathrm{1}-{\left|\mathit{\alpha }{h}_{\mathrm{j}}\right|}^{n-\mathrm{1}}G\left(h_{\mathrm{j}}\right)}{\mathrm{1}-{\left|\mathit{\alpha }{h}_{\text{ae}}\right|}^{n-\mathrm{1}}G\left(h_{\text{ae}}\right)}\right)}^{\mathit{\gamma }}{\left(\frac{h}{h_{\mathrm{j}}}\right)}^{-\mathrm{1.5}},& h_{\mathrm{d}}<h\le h_{\mathrm{j}},\\ {\left(\frac{G\left(h\right)}{G\left(h_{\text{ae}}\right)}\right)}^{\mathit{\tau }}{\left(\frac{\mathrm{1}-{\left|\mathit{\alpha }h\right|}^{n-\mathrm{1}}G\left(h\right)}{\mathrm{1}-{\left|\mathit{\alpha }{h}_{\text{ae}}\right|}^{n-\mathrm{1}}G\left(h_{\text{ae}}\right)}\right)}^{\mathit{\gamma }},& h_{\mathrm{j}}<h\le h_{\text{ae}},\\ \mathrm{1},& h>h_{\text{ae}}.\end{array}\right.\end{array}\end{array}$$

Obviously, the last term on the right-hand side represents the hydraulic conductivity of the soil for liquid water. Here, $K_{\mathrm{s},\mathrm{c},\text{jun}}$, γ, and τ are fitting parameters that appear only in the expression for the UHCC, additionally to h_ae, α, and n, which also feature in the RIA model for the SWRC. If desired, some or all of the latter can be fitted independently from those of the SWRC. Because K_s,a is not a fitting parameter, the junction model has one parameter fewer than the models that average or add domain conductivities. In fact, it has the same number of parameters as the Kosugi model (de Rooij, 2024a) that does not consider adsorbed water. The value of h_j follows from that of h_d and n according to Eq. (1c). Hence, only if h_d and/or n are fitted separately for the UHCC can h_j differ from h_jθ. The curve described by Eq. (9) has a discontinuous derivative at h_j, allowing it to reproduce the changing slope of the UHCC observed for some soils.

It is worth noticing that in the limit for h≪0, the following simplification holds for the wet branch of Eq. (9), where $h_{\mathrm{j}}<h\le h_{\text{ae}}$.

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\left({\displaystyle \frac{G\left(h\right)}{G\left(h_{\text{ae}}\right)}}\right)}^{\mathit{\tau }}{\left({\displaystyle \frac{\mathrm{1}-{\left|\mathit{\alpha }h\right|}^{n-\mathrm{1}}G\left(h\right)}{\mathrm{1}-{\left|\mathit{\alpha }{h}_{\text{ae}}\right|}^{n-\mathrm{1}}G\left(h_{\text{ae}}\right)}}\right)}^{\mathit{\gamma }}\\ & \approx {\left({\displaystyle \frac{h}{h_{\text{ae}}}}\right)}^{\mathit{\tau }(\mathrm{1}-n)},h\ll \mathrm{0}\end{array}\end{array}$$

This requires τ to be positive, which is stricter than the limits set by Peters et al. (2011) and Peters (2014). It is physically plausible that the second factor on the left-hand side decreases as h decreases. This is the case when γ is positive.

For coding purposes, Eq. (9) and the corresponding expression for Kosugi's model as formulated by de Rooij (2024a) can be cast in a form that separates terms that need to be calculated once in the case that SWRC parameters are fixed, terms that need to be calculated once for every iteration, and terms that need to be calculated for every iteration and every matric potential corresponding to an observed point on the UHCC. The reformulated equations are given in Appendix A.

3 Fitting the model parameters

The junction model (or any of the other UHCC model accommodated by KRIAfitter) does not require that the parameter values fitted for the SWRC are assumed to be valid for the UHCC as well. Nevertheless, physical consistency between the SWRC and the UHCC requires that θ_s and the matric potential at which liquid water is no longer present in the soil (calculated as (1+c)h_dθ for the SWRC) are the same for both curves. Hence, before KRIAfitter 1.0 (de Rooij, 2024b, where the code, user manual, and example input and output files can be downloaded) can be run to determine the values of the parameters of the chosen UHCC model for a particular soil, the parameters of the SWRC of Eqs. (1a–d) need to be fitted using RIAfitter 2.0 (de Rooij, 2024c, which also hosts its user manual) or higher versions, once they become available. The water content at saturation has to be prescribed in the input for KRIAfitter because θ_s does not appear in the equations for any of the hydraulic conductivity models. The value of (1+c)h_dθ should not be fitted separately for the SWRC and the UHCC because that would create a physically impossible situation where either a non-zero liquid-water conductivity exists although no liquid water is present or liquid water is present but is rendered immobile because the liquid-water hydraulic conductivity is zero. Hence, θ_s is fixed at the SWRC value, and h_d is calculated according to Eq. (2b) from SWRC parameters that are provided on input.

Both KRIAfitter and RIAfitter are coded in Fortran2008, compiled on CygWin's Fortran compiler for 64-bit Windows computers, and tested on such a computer. All input and output consist of ASCII files. Typical run times vary from seconds to several minutes.

All hard-coded parameter values are placed in modules. This allows users to modify them as needed by changing a single value without having to inspect the entire code. The user manuals list all hard-coded values and their respective module.

3.1 The parameter fitting algorithm

KRIAfitter fits up to six or seven parameters for one of six models for the unsaturated soil hydraulic conductivity, depending on the model chosen by the user. For each of these models, the equivalent vapour conductivity can be included or excluded, yielding 12 different conductivity models.

Each of k observed data points provided on input needs to have an estimate of the standard deviation of the measurement error of the matric potential or the volumetric water content and of the associated soil hydraulic conductivity. These are used to calculate the weighted squared error term w_i assigned to the fit for data point i as follows:

$$\begin{array}{}\text{(11)}& w_i={\left({\displaystyle \frac{K_{\mathrm{o},i}-K_{\mathrm{f},i}}{\frac{\mathrm{d}K}{\mathrm{d}x}{\mathit{\sigma }}_{x,i}+{\mathit{\sigma }}_{\mathrm{K},i}}}\right)}^{\mathrm{2}},\end{array}$$

where K denotes the soil hydraulic conductivity (cm d⁻¹); x denotes either the matric potential h (cm H₂O) or the volumetric water content θ; and σ_x and σ_K are the standard deviation of the measurement errors in the independent variable (h or θ) and in K, respectively. Subscript “o” denotes an observed, and subscript “f” a fitted value. Subscript i denotes the number of the observation; hence, $i\in \mathit{\{}\mathrm{1},\mathrm{\dots },k\mathit{\}}$. The slope $\mathrm{d}K/\mathrm{d}x$ is evaluated at the observed value x_o,i of the independent variable and updated for every set of fitted parameter values.

The individual terms w_i are used to compute the root mean square error (RMSE) of the fitted vs. the observed values:

$$\begin{array}{}\text{(12)}& \text{RMSE}={\left[{\displaystyle \frac{\sum _{i=\mathrm{1}}^k{w}_i}{\sum _{i=\mathrm{1}}^k{\left(\frac{\mathrm{d}K}{\mathrm{d}x}{\mathit{\sigma }}_{x,i}+{\mathit{\sigma }}_{\mathrm{K},i}\right)}^{-\mathrm{1}}}}\right]}^{\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

In the code, the error standard deviations are scaled in such a way that the arithmetic mean of σ_K,i equals 1 % of the largest conductivity value in the input file with the observation data. The required scale factor is applied to σ_x,i as well to conserve the relative weight of all observation errors.

The code minimizes the RMSE (the objective function) by the shuffled complex evolution algorithm (SCE). A detailed description of the algorithm and discussions of its parameters are provided by Duan et al. (1992, 1993, 1994). A brief summary of the SCE algorithm is provided in the user manual (de Rooij, 2024b).

The number of dimensions of the parameter space (denoted as NrOfDimensions in the code and below) is the number of parameters whose values are actually fitted as opposed to being set to a fixed value by the user. The number of complexes (sets of points in the parameter space) is determined in one of two ways that can be selected by the user. It can adhere to Duan et al. (1994) for any number of fitting parameters but with a minimum of two. This results in two complexes for any of the conductivity models. Alternatively, the number of complexes is 2 × NrOfDimensions. This option significantly increases the execution time. Test calculations showed the quality of the fit occasionally improves. Per the guidelines of Duan et al. (1994), the size of the complexes and the number of evolution steps between shuffles equal (2 × NrOfDimensions) + 1. The size of the subcomplexes is NrOfDimensions + 1. The number of offspring in each evolution step is hard-coded as 1 (parameter SCEa in module SCEBisParameters).

Thus, the algorithm replaces the single poorest-performing point of each subcomplex during an evolution step by a new point with a lower RMSE or a random point if a better-performing point could not be found. It does so by

• first checking if a reflection point (further removed from the centroid of the subcomplex than the point that will be replaced) has a lower RMSE.

• If a parameter value of the reflection point is outside the permitted range, it is made equal to the minimum or maximum allowed value, whichever is the closest. If the value of h_j lies between h_ae and zero, and its value of h_ae is increased to make h_ae slightly larger (less negative) than h_j.

• If the reflection point does not perform better than the original point, a contraction point (closer to the centroid) is checked.

• If neither point improves the RMSE, a random point in the parameter space replaces the worst “parent” even if its RMSE is larger.

After all evolution steps have been completed, all points in the complexes (including those that evolved in the subcomplexes) are reshuffled and assigned to new complexes. The process of selecting subcomplexes and improving those through evolution steps is repeated for the new complexes. After each shuffle, the code checks for convergence. The code terminates when the parameter fits converged or when the user-set maximum number of evaluations of the objective function has been exceeded. To check for convergence, the code evaluates the following convergence criteria for each fitting parameter after each shuffle. If none of the parameters fails more than the user-specified maximum permitted number of criteria, convergence has been achieved.

1. In the best fits from the last set of shuffles, the range of a parameter does not exceed both the absolute and the relative user-specified tolerance. The number of shuffles in the set is the maximum of 2 times the number of non-fixed fitting parameters and the minimum required number of shuffles (hard-coded parameter MinimumStoredShuffles).

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 RMSE of the fit does not exceed a user-specified tolerance.

Criteria 2, 4, 6, 9, and 10 apply to all parameters simultaneously. For large numbers of parameters, a hypercube that occupies only a small fraction of the hypercube defined by the ranges of the fitting parameters can still allow for an excessively wide range for the individual parameters. The permissible volume therefore becomes smaller with an increasing number of fitting parameters. It is determined by a preset value (hard-coded parameter SCEHyperVolumeTolerance1D) raised to a power equal to the number of fitting parameters.

3.2 Validity of parameter values and their combinations

The SCE algorithm requires a finite parameter space in all dimensions, so permitted ranges of all parameters are required on input. The code checks the sign of all minimum and maximum values, corrects them if needed, and places the minimum and maximum value in correct order if needed. For h_ae, α, n, K_s,c, K_s,a, γ, and τ, physical or mathematical limits are implemented in subroutine RangeChecker, and the parameter range is forced to be within these bounds if needed. The routine also checks limits on θ_s and h_d, but since these parameters are fixed, these should be redundant. They are retained in case users wish to modify the code or use the routine for other purposes.

During the fitting process, for those models for which both K_s,a and K_s,c are fitting parameters, the code makes sure the former never exceeds the latter. Similarly, h_j should not be closer to zero than h_ae. If it is, the value of h_ae is modified. If that parameter is fixed, the value of n is modified. Roundoff errors can create a combination of values for h_d, h_j, and n that violate Eq. (1c). If that is the case, n is modified as well.

Table 1 Boolean variables/arrays that guide the parameter fitting process.

3.3 Configuring the fitting process

In addition to the Boolean input parameter “FewComplexes” that governs the number of complexes (Sect. 3.1), there are three others and a Boolean array that steer the fitting process (Table 1). This section explains their use.

Variable “UseMap” determines if the code generates a regular grid of map points with their associated RMSE that covers the parameter space. The density of the grid is calculated from the maximum number of points supplied by the user, but a minimum number of points is generated that overrides a smaller user-supplied number. For parameters that are log-transformed, the map points are equidistant on the logarithmic scale. The map is used to modify the parameter space in order to focus on a hypercube in which the lowest values of the RMSE are concentrated. If the code finds that these values occur close to a boundary of the parameter space, the range of a parameter can be expanded somewhat.

A small number of parameter optimization runs is performed (hard-coded as 3) within the updated parameter space. In optimization run 1, the first complexes are filled with points that have the lowest RMSE values. The second run adds some random noise to the parameter values of these points to have slightly different initial complexes. The remaining run fills the complexes with randomly generated points.

If no map is desired, a larger number of runs is performed (hard-coded as 50; see below), all of them with randomly generated points that fill the first complexes. This larger number of runs allows for the calculation of statistics of the parameter value populations.

Array “LogTransformedParameters” specifies for each parameter whether untransformed values or log transforms of their absolute values (base 10) should be fitted. In the latter case, the corresponding dimension of the parameter space will be log-transformed as well, and the reflection and contraction steps are performed on log-transformed values. The fitted values that are written to output are back-transformed, but the statistics (mean, standard deviation, correlations and covariances) are for the log-transformed absolute values.

Any parameter can be fixed at a desired value by setting the upper and lower limit equal to that value. In that case, dummy values need to be provided for the absolute and relative tolerances and for the corresponding element of LogTransformedParameters. For the specific case that h_ae and h_j should be equal to those of the SWRC, h_ae and n will have to be fixed to the SWRC values. This can be achieved through input parameter AlphaOnly. In that case, all input for h_ae and n is treated as dummy values, except for their identifying names that will be used in the output.

Finally, the user can specify if the untransformed or log-transformed unsaturated soil hydraulic conductivity should be fitted through input parameter LogK. In the latter case, both the observed and fitted conductivities are log-transformed before the RMSE is computed. In that case, the error standard deviations for the observed conductivities apply to the log-transformed values. The user manual offers some guidance on how to estimate these.

Generating a map of the RMSE in the parameter space of course increases computation time, especially if the number of dimensions is high. It can help constrain (or expand) the parameter space and give the initial complexes some guidance about the location of the global minimum. This can be advantageous if one is not sure about the range of the parameter values.

Executing many runs (in the case no map is generated) also increases processor time. If the output shows a large number of random points compared to the number of reflection and contraction points, it is possible that the minimum was found multiple times. If that is the case, it is recommended to increase the probability of convergence by allowing more criteria to fail (the test runs gave good results for setting this number to 7). When doing so, the maximum number of iterations can also be reduced. If it is limited to 3000, most fits will probably still converge.

When a map is generated, the output includes the mean and standard deviation of the RMSE values of a random sample of the map points, as well as the cumulative probability density function of the RMSE in the parameter space and its histogram (also based on the sample). This information shows if there is a large portion of the parameter space in which the RMSE hardly varies, if there is a well-defined minimum, etc. For each run, correlation matrices of the parameter values are provided based on a sampling of the final iterations. These samples will include randomly generated points, so the reported correlations will underestimate the true correlations.

If no map is generated, the mean, population standard deviation, and histogram of all fitting parameters (or their log-transformed values) are written to output, as well as covariance and correlation matrices. These statistics are based on the fits of all runs, so they are expected to reflect the true values well. This information can be used, for instance, to sample the multivariate distribution of the parameter values in order to carry out a Monte Carlo study to examine the effect of parameter uncertainty on the uncertainty of the output of numerical models for unsaturated flow (de Jong van Lier et al., 2024).

It may be useful to run KRIAfitter with a map of the RMSE, check the updated parameter ranges in the output (see the user manual for details), copy these to the input file, and run KRIAfitter again without a map. In this way, the parameter space is better defined, and one obtains reliable parameter statistics.

4 Model evaluation

4.1 Fitting procedure

The first 13 soils in Table 1 of Weber et al. (2019) were used to evaluate the conductivity models. Instead of Ver P36 C1g in that table (P26 is a typo), we used UNSODA soil 4142, which refers to the same pair of samples. For all these soils, the data points span a wide range of matric potentials, and both retention and conductivity data are available.

Table 2 The sets of fitting parameters used in the model evaluation. N.B. θ_s, and h_d are fixed for all sets.

First, the SWRC parameters were fitted, once with free h_d and once with (1+c)h_dθ fixed at −10^6.8 cm H₂O (Schneider and Goss, 2012). Parameters α_θ and h_dθ were log-transformed because they can vary over several orders of magnitude (de Rooij, 2022). The error standard deviations were based on guesstimates of the accuracy of observed water content and matric potentials. These resulted in error bands r with units depending on the variable. The corresponding error standard deviation of $r/\sqrt{\mathrm{12}}$ was found by assuming the error to be uniformly distributed within that band (Abramowitz and Stegun, 1964, p. 930). Because soil water retention data points obtained with pressure plate extractors (Dane and Hopmans, 2002) were found to be unreliable after the data sets were created (Bittelli and Flury, 2009), matric potentials below −1000 cm H₂O were assigned an error band of ± 0.1  h_o (with the subscript o denoting an observed value).

The parameter space bounds were improved by fits that used the map of the RMSE, and then the improved ranges were used to generate a fit without map. In a few cases (Pachappa, Rehovot sand, SM−1005, UNSODA 4650), the data points were distributed so unevenly that the weights of the data points needed to be modified to compensate for that. For all soils, the fit without a map and with fixed (1+c)h_dθ gave the lowest AICc value. Therefore, these parameter values were used in the fits of the UHCC models. The fits also provided two estimates of the bulk soil hydraulic conductivity at saturation based on Eqs. (1) and (15) of Timlin et al. (1999). The latter estimate could only be calculated if h_ae<0 and n<3 (de Rooij, 2024a). If the data range of conductivity data required it, these estimates were included in the file with observation data and given the same error standard deviations as the other points.

Then, the five conductivity models of de Rooij (2024a) and the junction model (Eq. 2b) were fitted, with α, K_s,c, and K_s,a log-transformed and the fits performed on log-transformed conductivity values. For each model, five sets of parameters were fitted (Table 2).

In set 2, all critical values of the matric potential (h_d, h_ae, and h_j) are set to the values determined for the SWRC. By fixing both h_j and h_d, n is determined as well through Eq. (1c). Of the parameters that also occur in the SWRC, only α remains as a fitting parameter (hence “alpha only”).

As was the case for the SWRC, these fits were first performed using the map of the RMSE to better define parameter ranges. Then the fit was repeated with the improved ranges. If the upper limit of K_s,a exceeded the lower limit of K_s,c, it was capped at the value of the latter. Up to 7 of the 10 convergence criteria were allowed to be missed, and the maximum number of evaluations of the objective function was 3000. The error standard deviations were estimated for the largest observed conductivity only (see the user manual for details) and then assigned to all data points so that the data were all weighted equally. The only exception was soil UNSODA 4010, where the first three observations and the two estimates of the saturated bulk conductivity $K_{\mathrm{s}}^{\mathrm{B}}$ had to be given more weight to obtain good fits. Vapour flow was always included. All data were from laboratory measurements, so the temperature was set at 20 °C.

This required 780 runs of KRIAfitter, in which 20 670 fits were performed in total (3 fits per run with a map, 50 without). In the vast majority of cases, fits converged after several hundred evaluations of the objective function. The best fits of the runs without map are discussed below.

4.2 Fitting results

The RIA parameterization (Eqs. 1a–d) fits the intermediate and wet range quite well for a wide range of shapes of the curves, as the graphs in Appendix B show. Because the retention data points in the dry range are unreliable (see above), their weighting factors are smaller than those of the other data points. With the pF (defined as log(−h) with h in cm equivalent water column) of oven dryness fixed at 6.8, the fitted dry branches are drier than the observations to varying degrees, consistent with the expectation that pressure plate data tend to underestimate |h|. For Pachappa, the fitted value of h_aeθ is −15.6 cm H₂O, which seems a bit arbitrary, but with data points missing for matric potentials between 0.0 and −47.0 cm H₂O, the code has insufficient information for a reliable estimate. The desirability of additional points to guide the estimation of h_aeθ also seems apparent to a more moderate degree for SM–1005 and UNSODA 2571.

Appendix B also shows the 30 fits of the UHCC for each soil separated into one graph each for every of the five sets of fitting parameters listed above. The curves are arranged in the same order as that list when read from left to right and top to bottom. The main findings from a visual analysis of the collection of graphs are presented below. For consistency between this text, Appendix B, and the user manual, the various conductivity models are referred to by their identifying three-character labels as follows:

• ADV: the unweighted additive model, with vapour flow (de Rooij, 2024a);

• AMV: the arithmetic mean model, with vapour flow (de Rooij, 2024a);

• GMV: the geometric mean model, with vapour flow (de Rooij, 2024a);

• HMV: the harmonic mean model, with vapour flow (de Rooij, 2024a);

• JUV: the junction model, with vapour flow (Eq. 2b);

• KGV: the Kosugi model with capillary water only, with vapour flow (de Rooij, 2024a).

The main points to observe from a qualitative analysis of the collection of graphs in Appendix B follow. Pachappa (Fig. B2) is the only soil for which K(θ) was fitted instead of K(h). Despite the number of fitting parameters ranging from three to seven, the fits show only minor variations limited to the dry range, which has a single data point.

Pachepsky (Fig. B3) has K(h) data over a wide range but not at saturation. In view of the small amount of noise in the data, K_s,c and K_s,a are both fitted without providing estimates of saturated bulk soil hydraulic conductivity $K_{\mathrm{s}}^{\mathrm{B}}$ (L T⁻¹). The estimates based on the fitted SWRC are 182 and 110 cm d⁻¹. For the fit with all SWRC parameters fixed (set no. 5 in Table 2), all fitted saturated conductivities are well below these and often even below observed values below saturation. For fitting parameter sets 1 through 4, only JUV and KGV have plausible values of $K_{\mathrm{s}}^B$. GMV and HMV underperform overall. JUV gives the best overall fits when all parameters are fitted (set 1) or τ is fixed at zero (set 3).

The UHCC for Rehovot sand (Fig. B4) has a steep descent that transitions abruptly into a gentler slope. This shape can only be reproduced by ADV and AMV and when all seven parameters are fitted (set 1) or τ is fixed at zero (set 3).

The shape of SM–6–62 (Fig. B5) is accurately fitted by ADV, AMV, and JUV if 6 to 8 parameters can be fitted. All do well for set 1 (all parameters free) and set 3 (one parameter fixed). ADV and AMV also give good fits for set 4 (with γ and τ fixed to Mualem's (1976) values, leaving only four fitting parameters for JUV), although the conductivity may stay high for too long as the soil desaturates.

The fits for SM–22–88 (Fig. B6) are the best for sets 1 and 3. JUV does particularly well if all parameters are free (set 1), and AMV and ADV have an extra parameter and still do well if τ is fixed (set 3).

The UHCC of SM–35–119 (Fig. B7) is fitted well by all conductivity models for all sets except set 4 (all SWRC parameters fixed). For set 1, the spread between the fitted curves is remarkably high in the data-free wet range.

For SM–41–127 (Fig. B8), ADV, AMV, JUV, and KGV all do well if all parameters are free (set 1). For the other sets, the somewhat isolated driest point is only fitted well by ADV and AMV. The observation range covers only a narrow pF–range, resulting in considerable spread outside this range. The spread is contained at saturation by the values estimated from the corresponding SWRC.

The data points for SM–1005 (Fig. B9) are grouped in two clusters. HMV and GMV consistently underperform. ADV and JUV give good fits for up to one fixed parameter (sets 1 and 3), and because of its extra parameter, ADV still does well with two fixed parameters (set 4). Interestingly, AMV gives a curvy fit for set 4 that matches the data well.

The UHCC of UNSODA 2571 (Fig. B10) is linear on the double-log scale and fitted well for all sets by all models except HMV. The spread in the dry range is large, with ADV and AMV separated from the other models.

The peculiar shape of the UHCC of UNSODA 4010 (Fig. B11) cannot be fitted properly by any model. ADV, JUV, and KGV give somewhat acceptable fits for set 2 (with α the only free SWRC parameter).

Table 3 Combinations of conductivity models and fitted parameter sets (see main text for the parameter sets corresponding to the listed numbers) that give the lowest root mean square errors (RMSEs) and values of the corrected Akaike's information criterion (AICc) for soils with acceptable fits of both the soil water retention curve and the unsaturated hydraulic conductivity curve.

The wet range of the UHCC of UNSODA 4031 (Fig. B12) is fitted poorly by HMV and especially GMV. For sets 1 and 3, all other models give good fits over the entire data range, with limited spread outside that range.

Table 4 Root mean square errors (RMSEs) and values of the corrected Akaike's information criterion (AICc) for additive (ADV), junction (JUV), and capillary water only (KGV) conductivity models fitted to data for Pachappa. Best values for a particular conductivity model are in italics. The overall best values are bold.

Table 5 Similar to Table 3 but for Pachepsky's soil.

For UNSODA 4142 (Fig. B13), most models give good and comparable fits for sets 1–3, except AMV, GMV, and HMV in the wet range for set 1 and AMV and GMV for set 3. ADV tends to give higher conductivities in the dry range than the other models. The fits of JUV and KGV are nearly identical for all sets.

None of the fits for UNSODA 4650 (Fig. B14) are very good because the data points appear in two disjointed clusters with a discontinuity between them that can only be somewhat captured by ADV and AMV.

In many cases, GMV and HMV underestimate the conductivity in the wet range. This is related to the effect of the value of K_a (much lower than K_c) on the geometric and harmonic means of the two. Even though the upper limit of K_s,c was permitted to be about 20 times the measured saturated conductivity during the fitting process, this often was still insufficient to let the mean reach the observed values.

Table 6 Similar to Table 3 but for SM–35–119.

Table 7 Similar to Table 3 but for UNSODA 2571.

In 75 % of the plots, KGV closely tracks JUV. For a model that lacks representation of water adsorbed in films, KGV performs remarkably well. AMV overlaps ADV in the wet range at times, but the curves tend to separate in the dry end. ADV often has the highest conductivity in the dry range of all models until they all coalesce on the vapour conductivity, which is the only non-zero component of the bulk conductivity for $h\le -{\mathrm{10}}^{\mathrm{6.8}}$ cm H₂O, where liquid water ceases to exist.

JUV and KGV have one parameter fewer than the other models, yet JUV in particular tends to be among the best overall fits for a given soil based on the visual inspection of the full set of graphs. If it does so even when sharing the top spot with other models in the same set (e.g. Pachepsky, SM–35–119, SM–1005), this means it can match the performance of other models with fewer parameters. This suggests that ADV and the averaging models may be overparameterized, especially if all their seven parameters are fitted. On the other hand, for the Rehovot sand, ADV and AMV are the only models that are versatile enough to produce good fits, so the extra parameter may be needed at times. Although AMV and ADV often diverge in the dry range, AMV does not appear to substantially increase the fitting prowess provided by the combination of JUV and ADV. From the fits it appears that with only ADV and JUV, the best possible result (based on visual inspection) can be found for all soils used in the test.

This can be explored further through the goodness-of-fit measures and parameter correlation matrices. The time needed do so for all 390 fits presented in Appendix B is prohibitive. Instead, of the 30 combinations of conductivity model and fitting parameter set for each soil, the ones that give the smallest values of the RMSE and of Akaike's information criterion corrected for small sample sizes (Hurvich and Tsai, 1989; denoted as AICc) are identified. Only soils for which the visual analysis showed that both the SWRC and the UHCC were fitted well at least once are included in Table 3. The soils not considered are Pachappa (questionable value of h_aeθ in the SWRC due to the lack of data in the wet range), UNSODA 4010, and UNSODA 4650 (both because none of the combinations could fit the conductivity curve well).

Based on the minimal RMSE, ADV performs the best for four soils and AMV and JUV for three. The most flexible parameter set (no. 1) prevails nine times and the second-most flexible one (no. 3) one time. But if the required number of fitting parameters is accounted for using the AICc instead of the RMSE, ADV and JUV both achieve the minimum AICc four times and AMV twice. Set 4, with γ and τ fixed to Mualem's (1976) values, performs the best four times; sets 1 and 3 twice; and sets 2 and 5 once.

A more detailed analysis is carried out for a subset of soils. Based on Table 3, only ADV and JUV are considered as they are the overall best-performing models, with KGV included to see to what degree the more advanced models outperform the classical model. This analysis serves two purposes: an evaluation of conductivity models and fitting parameters sets and a demonstration of ways to interrogate the output of KRIAfitter.

The members of this subset of soils must have a good fit of the SWRC so that spillover effects from that fit do not cloud the assessment of the conductivity models. The subset must also represent different shapes of the SWRC and the UHCC and different qualities of conductivity data. Pachepsky's soil has a sigmoidal SWRC and conductivity data that cover a wide range. SM–35–119 has an SWRC with three linear branches and conductivity data over a narrow range, with estimated saturated conductivity values. UNSODA 2571 has a power-law type SWRC (Brooks and Corey, 1964) and conductivity data over an intermediate range. The subset is completed by Pachappa because it is the only soil for which the conductivity data have θ as the independent variable.

The values of the RMSE and AICc of each of the fits discussed above are presented in Tables 4–7. If only one of two identical values is designated optimal in these tables, they are different when more significant digits are considered.

Values of AICc should never be compared between soils. Because the error standard deviations are scaled by the code using the maximum conductivity value in the input file with the measurements (see the user manual), RMSE values for different soil cannot be compared either.

Because ADV has one fitting parameter more than both JUV and KGV, it achieves the lowest RMSE for three out of four soils, albeit with a shared first place for Pachappa. JUV has the lowest RMSE for SM–35–119. The three RMSE values for Pachappa are nearly identical, and all are achieved by set 3, with τ fixed at 0.0 (Table 4). For each soil, the lowest RMSE is achieved by fitting either set 1 (one case: Pachepsky; Table 5) or set 3 (three cases: Tables 4, 5, and 7).

The more parsimonious sets of fitting parameters not surprisingly tend to achieve the smallest AICc values: sets 2, 3, 4, and 5 all achieve the lowest AICc for one soil. Note that UNSODA 2571 (Table 7) has only eight data points, which means that for six fitting parameters, AICc becomes infinite, and beyond that incalculable (see the user manual and Hurvich and Tsai, 1989, for details). JUV is the most successful conductivity model in terms of minimizing AICc, with KGV close behind. JUV has the lowest AICc for Pachappa with set 4, for SM–35–119 with set 3, and for UNSODA 2571 with set 5 (Tables 4, 6, and 7). In all cases, KGV performs equally well or slightly worse. On average, 6 parameters are needed to achieve the lowest RMSE and 4.25 to achieve the lowest AICc.

Figure 1 The fits that had either the smallest value of the root mean square error (RMSE) or Akaike's information criterion corrected for small samples (AICc) of all fits, i.e. those fits that are identified by the bold values in Tables 4–7. For Pachappa, the best-performing UHCC models were ADV (RMSE) and JUV (AICc), for Pachepsky ADV (both RMSE and AICc), for SM–35–119 JUV (both RMSE and AICc), and for UNSODA 2571 ADV (RMSE) and JUV (AICc).

Figure 1 shows the pairs of best fits based on RMSE and AICc values. Only the case for which the difference in the number of fitting parameters between RMSE and AICc-based fits is larger than two results in significantly different fits. The other three are indistinguishable or, in one case, identical.

The output of KRIAfitter when no map of the RMSE is generated permits an inspection of the correlation structure. The breadth of parameter correlations can be assessed by determining from the correlation matrix in file SCEfittingK.OUT the number of data pairs for which the absolute value of the correlation coefficients exceeds a threshold. Table 8 shows the results for this threshold set to 0.8. There is no consistency for any of the conductivity models or fitted parameter sets between soils. Apparently, the properties of the input data set have far more influence.

Table 8 Number of data pairs with |R|>0.8 for selected soils, conductivity models, and fitted parameter sets (see main text).

Table 9 lists the correlation coefficients (R) with the maximum absolute value. Intriguingly, the only consistency here appears to be that the sets with the largest (set 1) and smallest (set 5) number of fitting parameters have the fewest data pairs with very high correlation coefficients (R²>0.9, corresponding to |R|>0.95). For set 5, this is understandable: with three (JUV, KGV) or four (ADV) parameters resulting in three or six possible parameter pairs, there is limited opportunity for different parameter combinations that give the same RMSE across a range of values. And these will only give high R values if the resulting relationship between the parameter values is linear.

Table 9 The maximum value of |R| for any parameter pair for selected soils, conductivity models, and fitted parameter sets (see main text). Values that give R²>0.9 are in italics.

Table 10 The parameter pairs that give the R²>0.9 pair for selected soils, conductivity models, and fitted parameter sets (see main text).

Finally, Table 10 gives the parameter pairs for which R²>0.9. Here too, the dominant effect of the input data is apparent, but there also appears to be an effect of the set of fitting parameters. Parameter pair (α,τ) appears six times, five of which are for Pachappa. The three appearances of (α,n) are also for Pachappa and for fitting parameter set 4. Seven of the eight appearances of (γ,τ) are for set 2, and all four appearances of (n,γ) are for set 3 for models JUV or KGV. All things considered, no general guidelines can be given. Parameter correlations therefore need to be examined on a case-by-case basis. If highly correlated parameters are found, the means and standard deviations of the individual parameters that are also provided in file SCEfittingK.OUT can help determine a suitable value at which to fix one of the parameters, if so desired.

5 Conclusions

The graphs in Appendix B, the combinations of models and fitted parameter sets that minimize AICc in Table 3 and the quality of the fits in Fig. 1 all indicate that all UHCC curves that can be fitted by at least one of the full set of six models can probably also be fitted by ADV and/or JUV without the need for implausibly high values of K_s,c that arise from averaging hydraulic conductivities.

Many of the fitted curves depicted in Appendix B are remarkably similar despite different numbers of fitted parameters. Only 3 out of 10 soils in Table 3 have the same combination minimizing both the RMSE and the AICc. In only three cases in Table 3 is the AICc minimal for a model with six or more parameters. The two top panels of Fig. 1 have indistinguishable plots based on fits with different numbers of parameters. All this suggests that conductivity models with six or seven fitting parameters are overparameterized for the UHCC of many soils. The junction model and the classical Kosugi model have one parameter fewer than the other models tested in this study and cannot have more than six. Of the two, only the junction model has the ability to reduce its slope abruptly in the dry range (best visible for the UHCC plots for soils SM–6–62, SM–22–88, and SM–1005 in Appendix B). Of the models that incorporate different domains of water in the soil, the junction model is the only one that sidesteps the fundamental impossibility of correctly averaging or adding domain conductivities (de Rooij, 2024a). Together with its parsimony, this makes the junction model an attractive alternative to existing conductivity models.

Appendix A

The expressions for the unsaturated hydraulic conductivity have terms and factors that only depend on the parameters of the SWRC. Other factors depend on the conductivity parameters. Finally, there are terms and factors that depend on the matric potential.

The efficiency of the parameter fitting algorithm for the case when all SWRC parameters are fixed was improved by taking this into account. The equations of the conductivity models were rearranged into separate terms and factors. Below, the rearranged equations are given. The symbols for the new factors (two characters and a number) reflect their names in the parameter estimation code, except for KG5, which is termed SWRCGhae in the code.

For Kosugi's (1999) model (de Rooij, 2024a) and the additive model (de Rooij, 2024a), the rearranged equations are as follows:

$$\begin{array}{}\text{(A1a)}& {\displaystyle }\begin{array}{rl}& K\left(h\right)=K_{\mathrm{s}}\\ & \times \left\{\begin{array}{ll}\mathrm{0},& h\le h_{\mathrm{d}},\\ {\left\{\mathit{\beta }\text{ln}\left[\frac{h_{\mathrm{d}}}{h}\right]\right\}}^{\mathit{\tau }}{\left(\frac{\text{KG1}-\frac{\mathit{\beta }}{\left|h\right|}}{\text{KG4}}\right)}^{\mathit{\gamma }},& h_{\mathrm{d}}<h\le h_{\mathrm{j}},\\ {\left(\frac{G\left(h\right)}{\text{KG5}}\right)}^{\mathit{\tau }}{\left(\frac{\text{KG2}+\frac{F\left(h\right)-\text{KG3}}{\text{KG5}}}{\text{KG4}}\right)}^{\mathit{\gamma }},& h_{\mathrm{j}}<h\le h_{\text{ae}},\\ \mathrm{1},& h>h_{\text{ae}},\end{array}\right.\end{array}\text{(A1b)}& {\displaystyle }F\left(h\right)={\displaystyle \frac{{\left|\mathit{\alpha }h\right|}^{n}G\left(h\right)}{\left|h\right|}},\text{(A1c)}& {\displaystyle }\text{KG1}={\displaystyle \frac{\mathit{\beta }}{\left|h_{\mathrm{d}}\right|}},\text{(A1d)}& {\displaystyle }\text{KG2}=\text{KG1}-{\displaystyle \frac{\mathit{\beta }}{\left|h_{\mathrm{j}}\right|}},\text{(A1e)}& {\displaystyle }\text{KG3}={\displaystyle \frac{{\left|\mathit{\alpha }{h}_{\mathrm{j}}\right|}^{n}G\left(h_{\mathrm{j}}\right)}{\left|h_{\mathrm{j}}\right|}},\text{(A1f)}& {\displaystyle }\text{KG4}=\text{KG2}+{\displaystyle \frac{{\left|\mathit{\alpha }{h}_{\text{ae}}\right|}^n}{\left|h_{\text{ae}}\right|}}-{\displaystyle \frac{\text{KG3}}{G\left(h_{\text{ae}}\right)}},\text{(A1g)}& {\displaystyle }\text{KG5}=G\left(h_{\text{ae}}\right).\end{array}$$

In Eqs. (A1a) and (A1c), derived parameter β is calculated according to Eq. (1b) but using the parameter values of the UHCC instead of those of the SWRC. During parameter fitting, KG1 through KG5 only need to be calculated once if the SWRC parameters are fixed.

For the junction model of Eq. (2b) of the main text, the rearranged equations are as follows. The hydraulic conductivity for liquid water is denoted as K_l,jun.

$$\begin{array}{}\text{(A2a)}& {\displaystyle }\begin{array}{rl}& K_{\mathrm{l},\text{jun}}\left(h\right)=\\ & \left\{\begin{array}{ll}\mathrm{0},& h\le h_{\mathrm{d}}\\ \text{JU4}{\left(\frac{h}{h_{\mathrm{j}}}\right)}^{-\mathrm{1.5}},& h_{\mathrm{d}}<h\le h_{\mathrm{j}}\\ K_{\mathrm{s},\mathrm{c}}{\left(\frac{G\left(h\right)}{G\left(h_{\text{ae}}\right)}\right)}^{\mathit{\tau }}{\left(\frac{\mathrm{1}-{\left|\mathit{\alpha }h\right|}^{n-\mathrm{1}}G\left(h\right)}{\text{JU3}}\right)}^{\mathit{\gamma }},& h_{\mathrm{j}}<h\le h_{\text{ae}}\\ K_{\mathrm{s},\mathrm{c}},& h>h_{\text{ae}}\end{array}\right.\end{array}\text{(A2b)}& {\displaystyle }\text{JU1}=K_{\mathrm{s},\mathrm{c}}{\left({\displaystyle \frac{G\left(h_{\mathrm{j}}\right)}{G\left(h_{\text{ae}}\right)}}\right)}^{\mathit{\tau }}\text{(A2c)}& {\displaystyle }\text{JU2}=\mathrm{1}-{\left|\mathit{\alpha }{h}_{\mathrm{j}}\right|}^{n-\mathrm{1}}G\left(h_{\mathrm{j}}\right)\text{(A2d)}& {\displaystyle }\text{JU3}=\mathrm{1}-{\left|\mathit{\alpha }{h}_{\text{ae}}\right|}^{n-\mathrm{1}}G\left(h_{\text{ae}}\right)\text{(A2e)}& {\displaystyle }\text{JU4}=\text{JU1}{\left({\displaystyle \frac{\text{JU2}}{\text{JU3}}}\right)}^{\mathit{\gamma }}\end{array}$$

JU2 and JU3 need only be calculated once if the SWRC parameters are fixed and JU1 and JU4 for every iteration of the fitting process. If at least one of the SWRC parameters is fitted, all four terms have to be recalculated once per iteration.

Appendix B

This appendix shows the graphs of all fitted SWRCs and UHCCs. The SWRC has a single fit for each of the test soils (Fig. B1). The UHCCs of six conductivity models were fitted for five sets of fitting parameters, resulting in 30 fits for each soil (Fig. B2–B14). In the figure panels, the fitting parameter sets are labelled “all free”, “alpha only”, “Assouline”, “Mualem”, and “SWRC fixed” according to Table 2. For Pachappa, the unsaturated hydraulic conductivity was observed as a function of the volumetric water content and for all other soils as a function of the matric potential.

Figure B1 The observed and fitted soil water retention curves for 13 soils listed by Weber et al. (2019). The fits all had (1+c)h_d fixed at −10^6.8 cm H₂O. The data points for pF > 3 were given less weight to reflect the low reliability of pressure plate retention data (Bittelli and Flury, 2009).

Figure B2 Fitted unsaturated hydraulic conductivity curves (UHCCs) for six models and five sets of fitting parameters (both explained in the main text) for the Pachappa soil of Weber et al. (2019).

Figure B3 As Fig. B2 but for Pachepsky's soil.

Figure B4 As Fig. B2 but for Rehovot sand.

Figure B5 As Fig. B2 but for SM−6−62.

Figure B6 As Fig. B2 but for SM−22−88.

Figure B7 As Fig. B2 but for SM−35−119.

Figure B8 As Fig. B2 but for SM−41−127.

Figure B9 As Fig. B2 but for SM−1005.

Figure B10 As Fig. B2 but for UNSODA 2571.

Figure B11 As Fig. B2 but for UNSODA 4010.

Figure B12 As Fig. B2 but for UNSODA 4031.

Figure B13 As Fig. B2 but for UNSODA 4142.

Figure B14 As Fig. B2 but for UNSODA 4650.