The primary objectives of this study were to: (i) elucidate the impacts of nonlinear scale transformations on the shapes and parameter values of soil water release and moisture capacity curves; and (ii) demonstrate how implicit characteristics of some established soil water release and moisture capacity models can impact model-data fits and estimates of model parameters. Nonlinear scale transformations of the tension head (h) axis (e.g., log10h, h1/2) were found to distort release and capacity curve shapes, create fictitious curve inflections and modes, and occasionally erase visual evidence of actual inflections and modes. The popular van Genuchten–Mualem and Assouline–Grant models were shown to always generate a release curve inflection and a capacity curve mode, even when inflections and modes did not exist in the data, and this in turn caused poor model-data fits in the critical near-saturated region. The van Genuchten model with four independently fitted parameters and the Dexter–Weibull model could accurately fit data sets with no inflection or mode, but this resulted in a physically unrealistic zero-angle intersection between the release curve and the water content axis. It was concluded that nonlinear h axis transforms should not be used when determining inflections, modes, pore size distributions, soil structure parameters, or soil quality indexes from soil water release and moisture capacity data-sets. It was also recommended that more flexible release curve models should be developed that do not assume the existence of inflections and modes, and also produce physically realistic angles of intersection between the water content axis and the fitted model.
Les principaux objectifs de cette étude étaient les suivants : i) élucider l’impact de la conversion non linéaire des échelles sur la forme et la valeur des paramètres des courbes de la pression capillaire et de la rétention d’eau du sol, et ii) montrer comment les particularités implicites de certains modèles de la pression capillaire et de la rétention d’eau du sol peuvent modifier l’ajustement modèle-données ainsi que l’estimation des paramètres du modèle. Les auteurs ont constaté qu’en convertissant non linéairement l’échelle de l’axe de la pression interstitielle (h) (à savoir, log10h, h1/2), on fausse la forme des deux courbes, on crée une inflexion et un mode fictifs de la courbe et, parfois, on efface les preuves visuelles d’une inflexion et d’un mode véritables. Les auteurs montrent que les populaires modèles de van Genuchten–Mualem et d’Assouline–Grant confèrent toujours une inflexion à la courbe de la pression capillaire et un mode à celle de la rétention d’eau, même si les données en sont privées, ce qui entraîne un piètre ajustement modèle-données dans la zone critique, proche du point de saturation. Le modèle de van Genuchten à quatre paramètres ajustés indépendamment et celui de Dexter-Weibull pourraient bien s’ajuster aux jeux de données sans inflexion ni mode, mais l’intersection à angle nul qu’on obtient entre la courbe de la pression capillaire et l’axe de la teneur en eau est physiquement irréaliste. Les auteurs en concluent qu’on ne devrait pas se servir des transformées non linéaires de l’axe h pour déterminer les inflexions, les modes, la répartition des pores selon leur taille, les paramètres de la structure du sol ou les indices de la qualité du sol à partir des données sur la pression capillaire et la rétention d’eau. Ils préconisent aussi l’élaboration de modèles plus adaptables de la courbe de la pression capillaire qui ne présument pas l’existence d’inflexions ou de modes et qui produisent des angles physiquement plus réalistes de l’intersection entre l’axe de la teneur en eau et le modèle ajusté. [Traduit par la Rédaction]
1. Introduction
Effective agri-environmental management of field crop production requires detailed knowledge of the storage and transmission of water in soil. This in turn requires rigorous and realistic descriptions of the soil water release curve, and the soil moisture capacity curve.
The soil water release curve (also variously referred to as the “capillary pressure” curve, “characteristic” curve, “desorption or imbibition” curve, and “retention” curve) is used primarily for describing water storage, and consists of the equilibrium functional relationship between the amount and energy status of water in soil pores (Or and Wraith 2002). The amount of water in undisturbed soil is frequently expressed as water volume per unit bulk soil volume and is referred to as pore water content (θ), while water energy status is often expressed on a per unit water weight basis and referred to as pore water tension head (h). A soil's water storage characteristics are therefore often expressed by a release curve in the form of a θ versus h relationship (e.g., Fig. 1a). Important applications of the soil water release curve include determination of soil effective porosity (accounting for entrapped air and nondrainable soil water), field capacity (FC) water content, plant permanent wilting point water content, and root zone capacity for storing plant-available water and air—all of which are crucial soil attributes affecting crop productivity, irrigation scheduling, greenhouse gas generation, and water balance modelling (see e.g., chapter 3 in Kutilek and Nielsen 1994; chapter 21 in Hillel 1998; chapter 2 in Radcliffe and Šimůnek 2010).
Fig. 1.
Hypothetical natural exponential data-set plotted on a reduced linear h scale, h* (circles), and a reduced log10h scale, log10h* (diamonds): (a) normalized water release curve, Θ(h); (b) corresponding moisture capacity curve, −dΘ/dh (on h* scale) and − dΘ/d(log10h) = −ln(10)hdΘ/dh (on log10h* scale). h* = h/hI = h/hM; log10h* = log10h/log10hI = log10h/log10hM; hI = hM = k−1 = 400 cm locates the inflection and mode on the log10h axis. The “data” (circles and diamonds), hI and hM were generated using eq. 28 with scale parameter, k = 0.0025 cm−1. [Colour online.]
The soil moisture capacity curve (also known as the “water or hydraulic capacity” curve, “specific moisture” curve, “specific water capacity” curve, and “differential water capacity” curve) is used in the description of water transmission through soil. The moisture capacity curve is simply the first derivative (slope) of the water release curve, i.e., dθ/dh versus h (e.g., Fig. 1b), and it may be physically interpreted as the volume of water per unit bulk soil volume that is released or imbibed by a soil per unit change in tension head. Equation 1.1 shows how moisture capacity can be incorporated into a one-dimensional form of Richards (1931) equation for describing water transmission in porous media:
where z [L] is the vertical space coordinate and K(θ) [LT−1] is the soil unsaturated hydraulic conductivity versus water content function, i.e., K versus θ.
Water release and moisture capacity curves are also useful for describing soil physio-hydraulic attributes because of the functional relationship between tension head and equivalent pore radius (i.e., capillary rise equation or Jurin's law):
where h [L] is the pore water tension head, σ [MT−2] is the pore air–water interfacial surface tension, γ [°] is the water contact angle with the soil pore wall, r [L] is the soil pore radius (equivalent cylinder), ρW [ML−3] is the soil water density, g [LT−2] is the gravitational acceleration, and β [L] is the soil capillary length (ratio of upward surface tension force to downward gravitational force). Replacing h with r in the release curve (i.e., θ vs. r) gives a soil pore size cumulative distribution, and replacing h with r in the capacity curve (i.e., dθ/dh vs. r) gives a pore size frequency distribution (e.g., pp. 23–25, Kutilek and Nielsen 1994). The cumulative and frequency distributions of pore size can be used to “quantify” soil physical quality or health, as well as assess soil quality/health responses to changes in land management. For example, Reynolds et al. (2009) used a range of soil types and land managements to propose representative pore size frequency distributions and water release curves for soils with good physical quality; Ferreira et al. (2019) used changes in pore size frequency distribution to elucidate soil physio-hydraulic responses to agricultural liming; and Jensen et al. (2020) used both cumulative and frequency pore size distributions to assess the soil physical impacts of converting long-term grassland to annual arable cropping and vice versa.
Soil water release curves are typically constructed from discrete (h,θ) data points measured in the field or from laboratory soil samples. Because the data are often sparse, noisy, and unevenly spaced with h extending over several orders of magnitude, it is usually advantageous (and often necessary) to replace the data with a smooth, continuous θ versus h function that has been “fitted” to the data (Or and Wraith 2002). Ideally, the fitted function is differentiable throughout its entire domain so that the corresponding capacity function (dθ/dh vs. h) is easily calculated, smooth, and internally consistent with the release function (θ vs. h).
Due to the extreme complexity of soil pore water systems, θ versus h functions are almost exclusively empirical, and many parametric forms have been proposed over the past 60+ years (e.g., Gardner 1956; Brooks and Corey 1964; Visser 1966; Campbell 1974; Haverkamp et al. 1977; van Genuchten 1980; Fredlund and Xing 1994; Assouline et al. 1998; Groenevelt and Grant 2004; Dexter et al. 2008; Reynolds et al. 2021). The ability of these empirical functions to fit and predict θ versus h data (and dθ/dh vs. h data) varies substantially, however, as they contain both obvious and implicit characteristics that may or may not be compatible with actual data behaviour. For example, the continuous Assouline et al. (1998) and Groenevelt and Grant (2004) functions assume a release curve with gradual air entry followed by an inflection, while the discontinuous Brooks and Corey (1964) and Campbell (1974) functions assume abrupt air entry with no inflection.
As mentioned above, the relevant h-domain of soil water release and moisture capacity data usually extends over several orders of magnitude. For agronomic applications, the h-domain is typically 0 ≤ h ≤ 15 000 cm (from saturation at h = 0, to plant permanent wilting point at h = 15 000 cm), although the domain can extend to much greater h values for some drought-tolerant crops (Or and Wraith 2002). As a result, the h axis is frequently scaled using nonlinear transformations (usually log10h or logeh) to allow visualization of both data and fitted models over their entire domains (see e.g., figs. 1–3 in Assouline et al. 1998; fig. 1 in Reynolds et al. 2009). Unfortunately, nonlinear scale transformations can introduce artefacts into the presentation, analysis, and interpretation of water release and moisture capacity data which are not explained in soils textbooks, methods manuals, or publications, and often not recognized by practitioners.
In view of the above, the objectives of this study were to: (i) elucidate the impacts of the logbh scale transformation and its artefact effects on the shapes and parameter values of soil water release and moisture capacity curves; (ii) make recommendations regarding the use of nonlinear h axis transformations; and (iii) illustrate some implicit characteristics of several established soil water release and moisture capacity functions, or “models,” and demonstrate how these characteristics can impact model-data fits and estimates of model parameters. It should perhaps also be mentioned that this study makes no attempt to delineate the underlying physical or hydrologic conditions or processes responsible for release and capacity curve shapes and parameter values.
2. Materials and methods
2.1. Soil water release relationships
Some versatile and widely used soil water release relationships include those of van Genuchten (1980), Assouline et al. (1998), Groenevelt and Grant (2004), and Dexter et al. (2008). An adapted Weibull relationship was also proposed recently by Reynolds et al. (2021).
The van Genuchten (1980) model is often expressed in the form
where h [L] is the soil water tension head, θ [L3L−3] is the volumetric soil water content, θR [L3L−3] is the residual water content (i.e., θ → θR as h → ∞), θS [L3L−3] is the saturated water content (i.e., θ at h = 0), α [L−1] is an empirical scale constant, and n [−] and m [−] are empirical shape constants. In most soils research, α, n, m, and θR are treated as curve fitting parameters, and θR is constrained to values ≥ 0. Note also that the m parameter can be fitted independently, or it can be restricted to m = 1 − (1/n) with n > 1 when the Mualem (1976) pore connectivity model is assumed, or to m = 1 − (2/n) with n > 2 when the Burdine (1953) pore connectivity model is assumed (van Genuchten 1980). The van Genuchten function with independently fitted m will be referred to here as the van Genuchten-m-Independent model (vG-I), while van Genuchten with the Mualem (1976) restriction will be referred to as the van Genuchten–Mualem model (vG-M). The Burdine (1953) restriction was not considered in this study, as it is rarely used in soils research.
The Assouline et al. (1998) water release function has the form
where h, θ, and θS are as defined above, kA [Lc] and c [−] are empirical curve fitting constants related to scale and shape, respectively, and (hL, θL) is a specified tension head-water content coordinate. The (hL, θL) coordinate is often set to the plant permanent wilting point values (i.e., θL = θWP at hL = hWP = 15 000 cm; Assouline et al. 1998).
The originally bimodal water release model of Dexter et al. (2008) can be written in the monomodal form
where h and θ are as defined above, and A0 [L3L−3], A1 [L3L−3], and h1 [L] are empirical curve fitting constants with A0 + A1 = θS at h = 0.
A convenient form of the Groenevelt and Grant (2004) water release function is given by (Grant et al. 2010)
where h and θ are as defined above, k1 [L3L−3] is an empirical curve fitting constant related to soil water content, k2 [L] and c [−] are empirical curve fitting constants related to function scale and shape, respectively, and (hA, θA) is a specified tension head-water content coordinate. The (hA, θA) coordinate can take on any value between saturation and air dry (Grant et al. 2010).
The adapted Weibull function of Reynolds et al. (2021) is given by
where h, θ, θR, and θS are as defined above, and kW [L−c] and c [−] are empirical curve fitting constants related to function scale and shape, respectively. The h0 [L] value in eq. 6 may be used to specify the air-entry value for soil water desorption from saturation (i.e., the tension head at which air first enters a saturated soil), or the water-entry value for soil wetting from dryness (i.e., the tension head at which wetting soil spontaneously saturates).
Comparative application of the above functions is most rigorous when they are all anchored at the same data point on the soil water release curve. Since the van Genuchten function (eq. 2) was designed to have no distinct air/water-entry value and to yield θ(h) = θS at h = 0, eqs. 3–6 were adjusted to behave similarly. This was achieved for the Assouline et al. (1998) model (eq. 3) by setting hL = ∞ to obtain
for the Grant et al. (2010) model (eq. 5) by setting (hA,θA) = (0,θS) and k1 = (θS – θR): for the Dexter et al. (2008) model (eq. 4) by setting A0 = θR and A1 = (θS – θR): and for the adapted Weibull model (eq. 6) by setting h0 = 0 and rearranging to produce:
For these conditions, it is easily shown that eqs. 7 and 8 are identical, and that eq. 9 is a special case of eq. 10 with c = 1 and kW = 1/h1. Hence, the four models collapse to just two models when they are anchored at soil saturation with no air/water entry value. That is, eqs. 7 and 8 become
which will be referred to here as the “Assouline–Grant” (A-G) model where q [Lp] and p [−] are empirical curve fitting constants controlling model scale and shape, respectively; and eqs. 9 and 10 become which will be called the “Dexter–Weibull” (D-W) model with scale constant, k [L−c], and shape constant, c [−]. For the most part, release curve range or position on the h axis is determined by the scale constant (α, q, k), and release curve slope is determined by the shape constant (n, m, p, c). The vG-M (eq. 2), vG-I (eq. 2), A-G (eq. 11), and D-W (eq. 12) models were compared using metrics and data-sets described below in Sections 2.4 and 2.6, respectively.
2.2. Soil moisture capacity relationships
As mentioned above, the soil moisture capacity curve is simply the first derivative, or slope, of the soil water release curve, i.e., dθ/dh versus h. Hence, the vG-I, vG-M, A-G, and D-W moisture capacity functions are given, respectively, by
and
It is readily shown using limit theory that dθ/dh → 0 as h → ∞ for all four relationships, and that dθ/dh → 0 as h → 0 for eqs. 13.2 and 14. For eqs. 13.1 and 15, however, the value of dθ/dh as h → 0 depends on the value of n and c, respectively. Specifically, dθ/dh = 0 at h = 0 if n and c > 1; dθ/dh → −∞ as h → 0 if n and c < 1; dθ/dh = −mα(θS − θR) at h = 0 if n = 1 (eq. 13.1); and dθ/dh = −k(θS − θR) at h = 0 if c = 1 (eq. 15). This means as a consequence that vG-M and A-G release curves must always be sigmoid in shape, and intersect the θ axis at a right angle (since dθ/dh → 0 as h → 0). The vG-I and D-W release curves, on the other hand, are convex-monotonic in shape with a zero-angle intersection on the θ axis when n and c < 1 (since dθ/dh → -∞ as h → 0), convex-monotonic in shape with an acute angle intersection when n = c = 1 (since dθ/dh = −mα(θS – θR) or − k(θS − θR) at h = 0), and sigmoid in shape with a right-angle intersection when n and c > 1 (since dθ/dh = 0 at h = 0). The above limits further dictate that vG-M and A-G moisture capacity curves must be bell-shaped with zero endpoints, while vG-I and D-W moisture capacity curves can be bell shaped with zero endpoints (n, c > 1), convex-monotonic with acute angle intersection on the θ axis (n = c = 1), or convex-monotonic with zero-angle intersection on the θ axis (n, c < 1).
Moisture capacity measurements are rare and difficult to obtain directly, but they can be reasonably estimated from the slope of release curve data, providing the data are not too noisy. Several estimation approaches are possible (e.g., slopes calculated from fitted splines); however, simple first-order finite difference slope estimates seem to work reasonably well, i.e.,
where and N is the number of measured (hi, θi) data pairs in the release curve. Equation 16 was consequently used to estimate moisture capacity data for comparison to the corresponding vG-I, vG-M, A-G, and D-W moisture capacity functions (eqs. 13.1, 13.2, 14, and 15, respectively).
2.3. Model fitting to data
The vG-I, vG-M, A-G, and D-W models were fitted to soil water release data (Section 2.6) using nonlinear least squares regression applied via the Solver® algorithm in Excel®. The adjustable curve fitting parameters included α, n, m, q, p, k, c, and θR, with initial guess values obtained by preliminary “trial-and-error” matching of model prediction to release data. Fits were conducted using several initial guess values to ensure that consistent global solutions were obtained. Selected Solver® options included generalized reduced gradient numerical iteration, slope calculation using central derivatives, automatic scaling, and a 10−10 convergence criterion. The tension head (h) data were untransformed to minimize potential curve-fitting bias (e.g., Miller 1984), and fitting parameter values were constrained to be ≥ 0 to ensure physically plausible results (especially for θR). The corresponding moisture capacity functions (eqs. 13.1, 13.2, 14, and 15) were not fitted because the moisture capacity “data” were derived using eq. 16, and therefore both approximate (first-order finite difference estimates) and not independent from the release curve data.
2.4. Model-data fit metrics
The accuracy and plausibility of the model-data fits were quantified using adjusted coefficient of determination (RA2), normalized mean bias error (MBEN), normalized standard error of regression (SERN), and normalized probability ranking (PX).
The RA2, MBEN, and SERN values are calculated using
where
is the regression sum of squared errors, 
is the regression total sum of squares, 
and 
are, respectively, the model-predicted and measured θ values at each h value, 
is the mean of the measured θ values, N is the number of data points (measurements), and L is the number of model fitting parameters (L = 3 for vG-M, A-G, and D-W; L = 4 for vG-I). The RA2 value indicates degree to which the fitted model explains systematic data variation, and it accounts for the number of model fitting parameters (L). MBEN (also known as “normalized mean prediction error,” MPEN) indicates degree of systematic model bias, with positive values indicating net overestimate of the data by the fitted model, and negative values indicating net underestimate. The SERN value indicates the “predictive power” of a fitted model with L fitting parameters, where SERN = 0 indicates perfect prediction (i.e., no discrepancies between predicted and measured values). Generally speaking, larger RA2, smaller MPEN, and smaller SERN indicates better goodness of model-data fit, with “perfect fit” indicated by RA2 = 1, MPEN = 0, and SERN = 0. As described in Reynolds et al. (2021), 0 ≤ SERN < 10% indicates excellent model prediction of the data; 10% ≤ SERN ≤ 20% signifies good prediction; 20% < SERN ≤ 30% indicates fair prediction; and SERN > 30% signifies poor prediction. Further detail on the RA2, MBEN, and SERN calculations can be found in Reynolds et al. (2021).The PX metric ranks the likelihood or plausibility of the models from the perspective of optimized goodness of model-data fit and parsimony (Banks and Joyner 2017). In this study, PX is given by
where PvG-I, PvG-M, PA-G, and PD-W are, respectively, the relative probability (0 ≤ PX ≤ 100%) that vG-I, vG-M, A-G, or D-W is the “most likely” or “most probable” of the four models. The ωvGI, ωvGM, ωAG, and ωDW parameters are relative weights (ωvGI + ωvGM + ωAG + ωDW = 1) based on the corrected Akaike information criterion (AICC) for each model: where SSEX is the regression sum of squared errors for the fitted model in question, and N and L are as defined above. The fitted model with the largest PX value (and smallest AICC value) is deemed the “most likely/probable” of the four models, and thereby the most suitable for representing the release curve data (at least from a statistical perspective). Further detail on the theory and calculation of PX is given in Banks and Joyner (2017) and Reynolds et al. (2021).
The RA2, MBEN, SERN, PX, and AICC metrics were used to determine goodness of model-data fit and relative probability of the vG-I, vG-M, A-G, and D-W release curve models (θ vs. h). RA2, SERN, and AICC were also used to quantify the fit between model-predicted moisture capacity (eqs. 13.1, 13.2, 14, 15) and the estimated moisture capacity data (obtained viaeq. 16).
2.5. Soil water release and moisture capacity parameters
Several characterizing parameters can be derived from fitted soil water release and moisture capacity models (see e.g., Reynolds et al. 2021). Two of the most important include the (hI, θI) coordinate that locates the inflection (maximum slope) of the release curve and the (hM, dθ/dh|M) coordinate that locates the mode (maximum value or peak) of the capacity curve. The inflection and mode are important because they are related in physically fundamental (albeit still largely unknown) ways to changes in soil mechanical behaviour, soil health, and the ability of soil pores to store and transmit water and air. For example, the inflection has been related to soil physical quality and crop rootability (Dexter 2004a); to soil workability, friability, and hard setting (Dexter 2004b); and to a physically based matching point in hydraulic conductivity models (Dexter 2004c; Ghanbarian and Skaggs 2022). Dexter and Bird (2001) relate the inflection to the optimum water content for tillage; Dexter and Richard (2009) use the inflection to demark the boundary between soil structure/tillage pores and soil texture pores; andAssouline and Or (2014) relate the inflection to loss of hydraulic continuity and thereby a sudden decrease in soil drainage rate. In addition, Liakopoulos (1965), Silva et al. (2014), Liang et al. (2016), Turek et al. (2019), and many others have used the inflection to demark the soil water pressure head and water content at FC. The capacity curve mode, on the other hand, gives the dominant or characteristic pore size in the soil (e.g., pp. 22–25 in Kutilek and Nielsen 1994), and has been used (along with capacity curve shape) to quantify the effects of soil structure, land management, crop type, amendment application, etc. on the soil's water-air storage capacity, permeability, and physical quality/health (e.g., Durner 1994; Dexter et al. 2008; Dexter and Richard 2009; Reynolds et al. 2009, 2014; Ghanbarian and Skaggs 2022). Note also that because inflection and mode are both determined by maximum release curve slope, hI must have the same value as hM.
The tension head at the inflection and mode is obtained by setting the curvature of the fitted release function to zero, solving for h to obtain hI and hM, and ensuring that the curvature changes sign on either side of hI and hM (see e.g., Reynolds et al. 2021 for details). The release function curvature is given by the second derivative with respect to h, i.e., 
. As soil water release and moisture capacity curves are often strongly right skewed, they may be plotted on a linear h axis or on a transformed h axis. The type of h axis used changes the h expression for inflection and mode, however. For example, release function curvature on a logbh scale, 
, is equivalent to 
, where b is the logarithm base (van Genuchten 1980; Durner 1994; Dexter et al. 2008). The resulting hI and hM expressions for the vG-I and vG-M functions are
The matching θI and dθ/dh|M values are then obtained by substituting the appropriate hI or hM expression (or numerical value) into eqs. 2, 11, and 12, and into eqs. 13.1, 13.2, 14, and 15.
2.6. Soil water release and moisture capacity data-sets
Six illustrative soil water release and moisture capacity data-sets were used, including two hypothetical data-sets based on analytic functions, and four measured data-sets obtained from porous media. One hypothetical data-set was generated using the normalized natural exponential function:
where k = 0.0025 cm−1 and 0 ≤ h ≤ 5000 cm; and the other hypothetical data-set was generated using the normalized vG-M function: where α = 0.0119 cm−1, n = 2.7 (−), and 0 ≤ h ≤ 15 000 cm. These hypothetical data-sets simulate the two most common release curve and capacity curve shapes, i.e., release and capacity curves both convex-monotonic (eq. 28), and sigmoidal release curve with bell-shaped capacity curve (eq. 29). The four measured data-sets were obtained (using the tension table, tension plate, and pressure plate extractor methods in chapter 71 of Carter and Gregorich 2008) from repacked cylinders of glass beads (Spheriglass® solid glass microspheres, A-Glass, 2024), and from intact cores of Berrien sandy loam soil, Harrow loam soil, and Brookston clay loam soil (Table 1). The hypothetical data-sets were used to illustrate scaling artefacts induced by logbh transformation of the h axis, with moisture capacity calculated using dθ/d (logbh) = ln(b)h(dθ/dh), as mathematically required (e.g., Durner 1994; Nekola et al. 2008). The measured data-sets were used to elucidate implicit model characteristics, and analyzed using linear h, with (h, θ) data points located at h = 0, 5, 10, 30, 50, 100, 225, 350, and 15 000 cm.
Table 1.
Selected porous medium characteristics including mean grain size proportions (sand, silt, and clay), organic carbon content (OC), oven-dry bulk density (BD), and texture classification.
3. Results and discussion
3.1. Scale transformation artefacts and implications
As mentioned above, soil water release and moisture capacity data extend over a large h range (0 ≤ h ≤ 15 000 cm), and are often strongly right skewed. One consequence of this is obscuration of near-saturated θ and dθ/dh information in the small (near-zero) h range when the data are plotted on a linear h axis. To remedy this, the data are often plotted using a base 10 logarithmic h axis (e.g., fig. 3.15 in Or and Wraith 2002), which increases or “stretches out” the spacing between small h data points to better reveal near-saturated (and hydrologically important) θ and dθ/dh characteristics (Durner 1994). Log transformation not only stretches out small h data; however, it also decreases or “compresses” the spacing between large h data points (e.g., Nekola et al. 2008). An unavoidable consequence of these opposing processes is to create an artefact inflection in the release curve, and an artefact mode or peak in the corresponding capacity curve. Note as well that log10h transformation (and many other nonlinear transformations) actually forces the wet end of the capacity curve to zero because dθ/d(log10h) = ln(10)hdθ/dh = 0 at h = 0 regardless of the value of dθ/dh. The impacts of log10h transformation are illustrated below using the two hypothetical data-sets (eqs. 28 and 29).
By definition, natural exponential y versus x data (eq. 28) is convex-monotonic with no inflection or mode (e.g., Nekola et al. 2008); and this is indeed the case for our exponential θ and dθ/dh data-sets when plotted on a linear h axis, where it is seen that both θ (Fig. 1a) and dθ/dh (Fig. 1b) are continuous convex negative-slope monoclines throughout their h-domains. When the same data are plotted using a log10h axis, however, a transformation-induced “artefact” inflection appears causing θ versus log10h to become sigmoidal (Fig. 1a), and an artefact mode appears causing dθ/d(log10h) versus log10h to become bell shaped with zero endpoints (Fig. 1b). Note also that the artefact inflection and mode occur at hI = hM = k−1, where k is the exponential function scale parameter (eq. 28). Hence, the scale parameter determined the “neutral” or “pivot” point on the log10h axis (i.e., hI = hM = k−1), with progressive data stretching where h < hI or hM, and progressive data compression where h > hI or hM. Data stretching and compression occurs to varying degrees with many other nonlinear scale transformations commonly used in the soils and geology literature, as exemplified in Fig. 2 for log10h, lnh, log2h, and h(1/2) transforms applied to the natural exponential data-set. Note also from Fig. 2 that curve shape, inflection coordinates, and mode coordinates change with the h-transform used; and contrary to the conjecture of Dexter and Bird (2001), log10h has no special physical or theoretical significance relative to linear h or to any other nonlinear transform.
Fig. 2.
Hypothetical natural exponential function (eq. 28) plotted using reduced linear h (i.e., h*), log10h, lnh, log2h, and h1/2 scales: (a) soil water release curve; (b) soil moisture capacity curve. h* = h/hk where hk = k−1 = 400 cm. Squares locate inflections on nonlinear axes; diamonds locate modes on nonlinear axes; natural exponential scale constant, k = 0.0025 cm−1. [Colour online.]
Regardless of h axis scale, truly sigmoidal soil water release data always displays an inflection point, and the corresponding moisture capacity curve is always bell-shaped with a peak or mode. The curve shapes and the inflection and mode coordinates change; however, depending on whether the h axis is linear or logarithmic. Using the hypothetical sigmoidal data-set to illustrate (eq. 29), curve shapes morphed from strongly right-skewed on linear h* axis, to near-symmetrical on log10h* axis (Fig. 3); and release curve inflection changed coordinate location from (0.71,0.74) on linear h* to (1,0.55) on log10h* (Fig. 3a), while capacity curve mode changed coordinate location from (0.71,0.0068) on linear h* to (1,1.32) on log10h* (Fig. 3b). As with the hypothetical exponential data, plotting on a logarithmic h axis progressively stretched the release and capacity curve data as h decreased below hI or hM, and it progressively compressed the data as h increased beyond hI or hM (Fig. 3). Unlike the natural exponential data, however, hI and hM were determined by the vG-M shape parameter (n), as well as by the scale parameter (α) (eq. 29).
Fig. 3.
Hypothetical van Genuchten (vG-M) data-set plotted on a reduced linear h scale, h* (circles), and a reduced log h scale, log10h* (diamonds): (a) normalized water release curve, Θ(h); (b) corresponding moisture capacity curve, −dΘ/dh (on h* scale) and − dΘ/d(log10h) = −ln(10)hdΘ/dh (on log10h* scale). h* = h/hI = h/hM; log10h* = log10h/log10hI = log10h/log10hM; hI = hM = 71 cm locates the inflection and mode on the linear h* axis; hI = hM = 100 cm locates the inflection and mode on the log10h axis. The “data” (circles and diamonds), hI and hM were generated using eq. 29 with α = 0.011 869 cm−1, n = 2.7, and m = 1 – (1/n). [Colour online.]
It is clear from Figs. 1–3 that nonlinear h axis transformations can have potentially dramatic impacts on the shapes and interpretations of soil water release and moisture capacity curves. If the release and capacity data are actually convex-monotonic (e.g., Figs. 1 and 2; also Figs. 6–8), curve shape, inflection, and mode obtained using log10h (and many other nonlinear transforms) are pure “artefictions,” and have nothing to do with actual water release and moisture capacity characteristics, soil pore sizes, soil structure, or soil quality indexes (further illustrated below in Fig. 7 where log10h transformation predicted inflections that clearly did not exist in the data). If the data do actually possess an inflection and mode (e.g., Fig. 3), use of log10h (and many other nonlinear transforms) will distort curve shapes and change the locations and magnitudes of inflection and mode coordinates. Note also from Fig. 3 that plotting on a log10h axis erased all visual evidence of the data's actual inflection and mode, which were located at (h,Θ) = (0.71,0.74) and (h,−dΘ/dh) = (0.71,0.0068), respectively. Hence, analyses based on release or capacity curve shape, slope, inflection or mode may be compromised if a nonlinear h axis is used. For example, use of semilog release curves (θ vs. logbh) to locate an inflection point or determine index water contents or release curve shapes (e.g., Dexter and Bird 2001; Dexter and Richard 2009; Dexter et al. 2012; Sillers et al. 2001; Silva et al. 2014; Fu et al. 2021; Ghanbarian and Skaggs 2022) may lead to fictitious or physically unrealistic results. Similarly, use of log-transformed moisture capacity curves (−dθ/d logbh vs. logbh) may yield distorted pore size distributions, and potentially inaccurate (or fictitious) modal, index and optimal pore sizes (e.g., Dexter and Richard 2009; Dexter et al. 2008; Sillers et al. 2001; Ferreira et al. 2019; Jensen et al. 2020; Ghanbarian and Skaggs 2022).
The degree to which nonlinear axis transformations affect model estimates of inflection points and modes can be determined using an hI ratio, i.e.,
R for the vG-I model is therefore given by while R for vG-M has the form R for the A-G model is and R for the D-W model is given by where it is seen that R depends only on the function shape parameters (n, m, p, and c); and it is understood in eqs. 30.2 and 30.5 that the vG-I and D-W models do not have an inflection or mode when n and c are ≤1. The RvG-I, RvG-M, RA-G, and RD-W relationships (Fig. 4) all form concave-monotonic curves with R = 0 and R = 1 endpoints, and it is seen that R changes rapidly when n, p, or c are less than about 3. The RvG-M, RA-G, and RD-W relationships are single value, while RvG-I is a family of relationships (due to dependence on both n and m) that converge to the RD-W relationship as m gets large (because eq. 30.2 collapses to the same form as eq. 30.5 as m → ∞). It can also be shown that hI and R are independent of logarithm base (e.g., de Jong van Lier 2014).
Fig. 4.
Tension head at inflection (hI) or mode (hM) on linear h axis divided by hI or hM on log10h axis (RX) versus model shape constant (n, p, c). RvG-I = van Genuchten model with independently fitted m parameter (vG-I, shape constant n); RvG-M = van Genuchten model with Mualem restriction, m = 1 – (1/n) (vG-M, shape constant n); RA-G = Assouline–Grant model (A-G, shape constant p); RD-W = Dexter–Weibull model (D-W, shape constant c). m = 106.7 applies for vG-I fitted to glass bead data (Fig. 5); m = 0.26 applies for vG-I fitted to Brookston clay loam data (Fig. 8). The ratios for the three soils all fall within the marked rectangle with symbols eliminated to improve clarity. Details on the glass beads, loam, sandy loam, and clay loam are given in Table 1 and Figs. 5–8. [Colour online.]
For the repacked glass bead medium (Table 1), fitted n, p, and c values were greater than 3.9, and the corresponding R values were greater than 0.92, indicating that there was less than 8% difference between hI obtained using a linear h axis or a log10h axis (Fig. 4). For intact Berrien sandy loam, Harrow loam and Brookston clay loam, on the other hand, fitted n, p, and c were all less than 1.23 with R less than 0.063, indicating that hI obtained using log10h was more than 16 times greater than hI obtained using linear h (Fig. 4). In fact, the A-G estimate of hI for Brookston clay loam was more than 14 400 times greater when using log10h (hI = 597 cm) than when using linear h (hI = 0.04 cm). Hence, log10 transformation of the h axis induced minimal artefact effects for the glass beads, but huge effects for the soils. This is presumably related to the underlying pore size distributions of the materials. The glass beads have a narrow pore size distribution because of highly uniform bead shape and size (Table 1), and this generates a release curve that is steep, sigmoidal, and near symmetrical (e.g., Fig. 5a) with large n, p, and c values (Fig. 4). The intact natural soils, on the other hand, have much broader pore size distributions because of cracks, biopores, aggregates, plant residues, and particle size variation; and this in turn generates release curves that tend to be flat and strongly right skewed (e.g., Figs. 6a, 7a, and 8a) with small n, p, and c values (Fig. 4). There is also some evidence that the log10h artefact effect increases with increasing range in pore sizes, as the RvG-M value changed from 0.063 for Berrien sandy loam, to 0.035 for Harrow loam, to 0.013 for Brookston clay loam. Note as well from the R-ratio relationships (eqs. 30.1–30.5) that the two hI values approach equality (i.e., R → 1) only as n, p, and c → ∞, which corresponds to an extremely narrow pore size distribution and a near step-function release curve.
Fig. 5.
van Genuchten model with independently fitted m (vG-I), van Genuchten model with Mualem restriction (vG-M), Assouline–Grant model (A-R) and Dexter–Weibull model (D-W) fitted to glass bead data: (a) water release curve; (b) moisture capacity curve. Circles are measured water content; diamonds are finite difference (FD) estimates of moisture capacity (eq. 16); squares and triangles demark model inflections and modes (see text); vertical “T-bars” are standard error (N = 5). The vG-I fit required the Solver® convergence criterion to be relaxed from 10−10 to 10−3. [Colour online.]
Fig. 6.
van Genuchten model with independently fitted m (vG-I), van Genuchten model with Mualem restriction (vG-M), Assouline–Grant model (A-R) and Dexter–Weibull model (D-W) fitted to Berrien sandy loam data: (a) water release curve; (b) moisture capacity curve. Circles are measured water content; diamonds are finite difference (FD) estimates of moisture capacity (eq. 16); squares are model inflections and modes; vertical “T-bars” are standard error (N = 11). The vG-I and D-W models did not produce inflections or modes on the linear h scale. [Colour online.]
Fig. 7.
van Genuchten model with independently fitted m (vG-I), van Genuchten model with Mualem restriction (vG-M), Assouline–Grant model (A-R) and Dexter–Weibull model (D-W) fitted to Harrow loam data: (a) water release curve; (b) moisture capacity curve. Circles are measured water content; diamonds are finite difference (FD) estimates of moisture capacity (eq. 16); squares are model inflections and modes; triangles are model inflections on log10h scale; vertical “T-bars” are standard error (N = 10). The D-W inflection determined using log10h scale is (309,0.25). The vG-I and D-W models did not produce inflections or modes on the linear h scale. [Colour online.]
Fig. 8.
van Genuchten model with independently fitted m (vG-I), van Genuchten model with Mualem restriction (vG-M), Assouline–Grant model (A-R), and Dexter–Weibull model (D-W) fitted to Brookston clay loam data: (a) water release curve; (b) moisture capacity curve. Circles are measured water content; diamonds are finite difference (FD) estimates of moisture capacity (eq. 16); squares are model inflections and modes; vertical “T-bars” are standard error (N = 10). The A-G mode is off scale at (h, dθ/dh) = (0.04,0.0295). The vG-I and D-W models did not produce inflections or modes on the linear h scale. [Colour online.]
3.2. Implicit function characteristics
3.2.1. Model-data fits
The glass bead data and all four fitted models suggest a sigmoidal release curve and a bell-shaped capacity curve that are both approximately symmetric (Fig. 5). The vG-M and A-G models were effectively equivalent, yielding excellent visual model-data fits and fit metrics (RA2 > 0.998, MBEN ≤ 0.69%, SERN ≤ 2.93%), large and effectively equal probability rankings (PvG-M = 51.31%, PA-G = 46.30%), and identical residual water contents (θR = 0.038 m3 m−3) (Fig. 5a, Table 2). Although visual model-data fit and fit metrics for the vG-I model were actually slightly better than those for vG-M and A-G (Fig. 5a, Table 2), the parsimony penalty induced by one additional fitting parameter (m) caused a much lower probability ranking (PvG-I = 2.39%). The D-W fit to the release data was also visually good; however, the fit metrics were poorer (RA2 = 0.982, MBEN = −1.35%, SERN = 10.92%), a slightly smaller residual water content was obtained (θR = 0.035 m3 m−3), and the model's ranked probability was very low (PD-W ≈ 0) (Fig. 5a, Table 2). The fits of all four models to the corresponding moisture capacity data were relatively poor (RA2 ≤ 0.961, SERN ≥ 34.67%), which was not unexpected as the finite difference estimates of slope (eq. 16) are only first-order accurate and can magnify variability (Fig. 5b, Table 2). Nevertheless, all four models tracked the moisture capacity data reasonably well (Fig. 5b); and interestingly, the D-W model produced the second best fit to the capacity data (RA2 = 0.8332, SERN = 71.36%) after vG-M (RA2 = 0.9606, SERN = 34.67%) despite D-W being the worst fit to the release data (Table 2). This appears to be due largely to D-W fitting the wet-end moisture capacity estimates (h ≈ ≤30 cm) somewhat better than vG-I, vG-M, and A-G (Fig. 5b). All four models produced similar and plausible estimates of the tension head at the glass bead inflection and mode, i.e., hI = hM = 51.9, 55.3, 59.4, and 61.0 cm from A-G, vG-M, D-W, and vG-I, respectively.
Table 2.
Parameter values and associated metrics for model fits to measured water release data, θ versus h, and calculated moisture capacity data, −dθ/dh versus h (eq. 16), from uniform spherical glass bead media (Table 1, Figs. 5a and 5b).
For Berrien sandy loam, the D-W and vG-I models produced the best and second best fits, respectively, to both the water release data (RA2 = 0.9986 and 0.9982, MBEN = 0.03% and −0.02%, SERN = 1.11% and 1.27%) and the moisture capacity data (RA2 = 0.6552 and 0.4940, SERN = 82.32% and 99.72%) (Table 3). D-W had by far the largest ranked probability (PD-W = 97.94%), while vG-I was much lower (PvG-I = 2.05%) due to the parsimony penalty associated with having the extra “m” fitting parameter (Fig. 6, Table 3). These favourable results occur because both the data and the D-W and vG-I fits indicate convex-monotonic curves throughout the h range (0–15 000 cm), while A-G and vG-M produce inflections and modes at hI = hM = 4.4 cm and hI = hM = 10.7 cm, respectively (Fig. 6). The vG-M and A-G models consequently provide relatively poor representations of at least the wet-end data (h ≈ <40 cm), which is evidenced visually in Fig. 6, and quantitatively in Table 3 (e.g., RA2 values for moisture capacity are only − 0.408 and 0.089 for vG-M and A-G, respectively). Further evidence that vG-M and A-G were less appropriate models is provided by the fact that their fitted θR values were physically unrealistic (negative) when they were not constrained by the θR ≥ 0 criterion (data not shown), while D-W and vG-I produced plausible sandy loam values of 0.105 and 0.083 m3 m−3, respectively (Table 3). Although some researchers treat θR as an empirical curve fitting parameter (thereby allowing negative values, e.g.,van Genuchten and Nielsen 1985; Groenevelt and Grant 2001; Haverkamp et al. 2005), others provide strong thermodynamic and physical arguments that θR should be ≥ 0, and represent either the oven-dry (or air-dry) water content (e.g., Groenevelt and Grant 2004; Inforsato et al. 2020; Ross et al. 1991), or the water content at which hydraulic discontinuity first occurs on the main drainage curve (e.g., p. 30, Kutilek and Nielsen 1994; Dexter et al. 2012).
Table 3.
Parameter values and associated metrics for model fits to measured water release data, θ versus h, and calculated moisture capacity data, −dθ/dh versus h (eq. 16), from a Berrien sandy loam soil (Table 1, Fig. 6a, 6b).
All four models provided good fits to the Harrow loam water release and moisture capacity data-sets (Fig. 7, Table 4), with vG-I providing the best fit metrics for water release (RA2 = 0.9980, MBEN = −0.03%, SERN = 1.35%), vG-I and D-W effectively tying for highest ranked probability (PD-W = 41.48%, PvG-I = 31.29%), and vG-M and A-G providing the best and second best fits, respectively, for moisture capacity (RA2 = 0.9115 and 0.8528, SERN = 39.52% and 50.96%). Note, however, that vG-M and A-G achieved their good fits to the convex-monotonic data-set by generating inflection and modal tension heads that were near zero (i.e., hI = hM = 1.61 cm from A-G; hI = hM = 1.97 cm from vG-M; Fig. 7). It was additionally found that vG-M produced a physically unrealistic (negative) θR value when the Solver® fitting routine was not constrained by the θR ≥ 0 criterion (Table 4). Note as well that although distinct data and model inflections were evident on the log10h scale (data not shown), plotting these coordinates on linear h (Fig. 7a) shows that inflections did not exist, and were in fact nothing more than transformation-induced “artefictions” with no connection to the Dexter and Bird (2001) “maximum rate of increase in soil air content.”
Table 4.
Parameter values and associated metrics for model fits to measured water release data, θ versus h, and calculated moisture capacity data, −dθ/dh versus h (eq. 16), from a Harrow loam soil (Table 1, Figs. 7a, 7b).
The A-G model provided the best and most likely fit to the Brookston clay loam release data (RA2 = 0.9878, MBEN = −0.10%, SERN = 2.45%, PA-G = 84.89%), with D-W second (RA2 = 0.9813, MBEN = −0.13%, SERN = 3.03%, PD-W = 12.43%), vG-M third (RA2 = 0.9719, MBEN = 0.16%, SERN = 3.71%, PvG-M = 2.02%), and vG-I fourth (RA2 = 0.9801, MBEN = −0.19%, SERN = 3.12%, PvG-I = 0.66%) (Fig. 8, Table 5). The corresponding moisture capacity data, on the other hand, was fitted best by vG-M (RA2 = 0.9572, SERN = 28.24%), with A-G second (RA2 = 0.9004, SERN = 43.08%), vG-I third (RA2 = 0.7744, SERN = 64.84%) and D-W fourth (RA2 = 0.7228, SERN = 71.87%). It is also seen, however, that vG-I, vG-M, and A-G all produced unrealistic (negative) θR values when the Solver® fitting routine was not constrained to θR ≥ 0 (Table 5), and that vG-M and A-G fitted the strongly convex-monotonic data by predicting inflection and modal tension heads that were virtually zero (hI = hM = 0.04 cm from A-G, hI = hM = 0.62 cm from vG-M; Fig. 8). From a physical perspective, D-W therefore seems be the most appropriate model for representing the Brookston clay loam data-set, despite A-G having a better fit to the release curve data, and vG-I and vG-M having better fits to the capacity curve data.
Table 5.
Parameter values and associated metrics for model fits to measured water release data, θ versus h, and calculated moisture capacity data, −dθ/dh versus h (eq. 16), from a Brookston clay loam soil (Table 1, Fig. 7a, 7b).
3.2.2. Intrinsic model behaviours and potential implications
As discussed above, the vG-M and A-G models always produce sigmoidal water release curves and bell-shaped moisture capacity curves because α > 0 and n > 1 in vG-M (eqs. 2, 13, and 25), and q > 0 and p > 0 in A-G (eqs. 11, 14, and 26). Hence, fitted vG-M and A-G models always have inflections and modes regardless of data trends; and this may as a consequence produce poor and/or unrealistic model-data fits in one or more h-regions, such as demonstrated for h ≈ <40 cm in Fig. 6 and h ≈ <50 cm in Fig. 7. The vG-I and D-W models, on the other hand, allow sigmoid or convex-monotonic release curves (e.g., Figs. 5a–8a), and bell-shaped or convex-monotonic capacity curves (e.g., Figs. 5b–8b).
Poor model-data fits may cause important inaccuracies in estimation of soil physio-hydraulic parameters and in simulation of soil water dynamics, especially if the fitting error occurs in the “wet end” (small h range) where water storage and transmission often change rapidly with changing h (e.g., Durner 1994). Consider, for example, water transmission according to the diffusivity form of the Darcy flux equation for infiltration into uniformly unsaturated soil:
where v [LT−1] is the water flux density, dΘ/dx [L−1] is the water content gradient, and D(Θ) [L2 T−1] is the hydraulic diffusivity function: where −dΘ/dh is the moisture capacity function based on the normalized D-W soil water release function: and K(Θ) [LT−1] is the Brooks and Corey (1964) unsaturated hydraulic conductivity function: where KS [LT−1] is saturated soil hydraulic conductivity, and β [−] is a dimensionless empirical constant. For the parameter values selected, Fig. 9 shows that the hydraulic diffusivity function (Fig. 9b) can differ by more than 3 orders of magnitude as the moisture capacity function (Fig. 9a) varies from bell shaped with dΘ/dh = 0 at h = 0 (c = 1.5), to convex-monotonic with dΘ/dh = finite at h = 0 (c = 1.0), to convex-monotonic with dΘ/dh → -∞ as h → 0 (c = 0.5). For a given water content gradient (dΘ/dx), Darcy flux (v) could therefore also vary by more than 3 orders of magnitude depending on the underlying characteristics of the fitted water release and moisture capacity functions. Note also that the hydraulic diffusivity relationship [D(Θ) vs. Θ] changes shape substantially with changing c value—ranging from sigmoidal when c > 1, to concave-monotonic when c = 1, and to concave-nonmonotonic when c = 0.5 (Fig. 9b). Interestingly, horizontal infiltration into repacked laboratory soil columns (e.g., Bruce and Klute 1956; Ahuja and Swartzendruber 1972) and pressure-induced vertical outflow from intact cores (e.g., Parker et al. 1985) often produce diffusivity relationships that are roughly sigmoidal in shape, suggesting c > 1, while vertical infiltration into intact soil (e.g., Clothier and White 1981) can produce diffusivity relationships that are approximately concave-monotonic or concave-nonmonotonic, suggesting c ≤ 1. Hence, all c > 0 seem possible (or all n > 0 in the case of vG-I). As a result, flexible models capable of producing both sigmoidal or convex-monotonic release curves, such as vG-I and D-W, may be more generally applicable over models that produce only sigmoidal release curves, such as vG-M and A-G. Supporting this is the common observation (e.g., Durner 1994; Jabro et al. 2009) that in situ measurements and intact soil core samples often produce convex-monotonic release curves (e.g., Figs. 6a, 7a, and 8a), rather than sigmoidal release curves (e.g., Fig. 5a). If this is in fact the case, water release and moisture capacity curves from in situ measurements and intact soil core samples would often have fictitious inflections and modes and greatly misleading shapes (cf. Figs. 1 and 2) when presented using a nonlinear h axis, such as log10h.
Fig. 9.
Effect of Dexter–Weibull (D-W) shape constant, c, on: (a) relative moisture capacity function, C(h)* versus h; and (b) relative hydraulic diffusivity function, D(Θ)* versus Θ. Relative moisture capacity is given by C(h)* = dΘ/dh/(dΘ/dh|MAX) where dΘ/dh|MAX is maximum calculated moisture capacity. Relative diffusivity is given by D(Θ)* = D(Θ)/D(Θ)MAX, where D(Θ)MAX is maximum calculated diffusivity. Specified Dexter–Weibull scale constant, k = 0.1 cm−1; specified Brooks–Corey K(Θ) parameters, KS = 0.001 cm s−1 and β = 3. Note that the domain of h in panel (a) differs from the domain of Θ in panel (b), and that the Θ axis is reversed to match the orientation of the h axis. [Colour online.]
3.2.3. Model advantages and limitations
All four water release models appear capable of providing “excellent” or “good” fits to a diverse range of release curve data-sets (i.e., SERN < 11%; Tables 2–5), and they all were either “most probable” model (largest PX), or effectively tied for most probable model (second largest PX almost equal to largest PX), for at least one of the four data-sets (Tables 2–5). Nevertheless, all four models have implicit theoretical and/or practical advantages and limitations that may dictate how and when they are used.
The A-G model (eq. 11) has the theoretical advantage of being consistent with capillary rise theory (Jurin's law; eq. 1.2), which shows that porous medium water content and equivalent pore diameter are inversely related to pore water tension head (see also Assouline et al. 1998). An implicit limitation, however, is that A-G can produce only sigmoidal release curves and bell-shaped capacity curves with zero endpoints, which were shown above to often be inaccurate and unrealistic in the near-saturated range of intact soils (e.g., Figs. 6a, 7a, and 8a). In addition, Jurin's law becomes increasingly inaccurate as pore size increases beyond the soil's capillary length (eq. 1.2; Liu et al. 2018).
The vG-I (eq. 2) and D-W (eq. 12) models have the practical advantage of being able to produce convex-monotonic release and capacity curves in addition to sigmoidal release curves and bell-shaped capacity curves. As shown above, this flexibility often allows for more realistic fits to near-saturated water release and moisture capacity data (e.g., Figs. 6–8). On the other hand, implicit limitations of these models include lack of clear theoretical underpinning, and the need for a fourth fitting parameter in the case of vG-I (i.e., “m”).
The vG-M model (eq. 2) has the dual limitations of lack of theoretical underpinning, and inability to produce convex-monotonic release and capacity curves. Important practical advantages, however, are that it is by far the most widely applied of the four water release models, and it is available as a built-in option in almost all soil water simulation models.
An additional implicit limitation of all four models is their restricted flexibility and unrealistic options for predicting the release curve slope at the θ axis intercept (which is also the moisture capacity value, dθ/dh, at h = 0). The A-G and vG-M models have only the zero slope option (i.e., p, n > 1, dθ/dh = 0 at h = 0), while the vG-I and D-W models allow just three possibilities: zero slope (n, c > 1, dθ/dh = 0 at h = 0), infinite slope (n, c < 1, dθ/dh → −∞ at h = 0), and a fixed finite slope when an exponential release curve is fitted (n = 1, dθ/dh = −mα(θS − θR) at h = 0; c = 1, dθ/dh = −k(θS − θR) at h = 0). Soil water release data in the literature (e.g., Fig. 2 in Durner 1994) and in this study (Figs. 5–8) suggests, however, that the slope of the release curve at the θ axis intercept (h = 0) is always finite (i.e., never infinite), and can take on virtually any value from steep at an acute angle (e.g., Fig. 7a; top left of Fig. 2 in Durner 1994), to zero (e.g., Fig. 5a; bottom left of Fig. 2 in Durner 1994). Given that soil water transmission is maximum and soil hydraulic properties often change rapidly as h → 0 (e.g., Fig. 9; van Genuchten and Nielsen 1985), it seems advisable to develop more flexible models that can accurately fit near-saturated soil water release data including the θ axis intercept.
4. Conclusions
Nonlinear transformations of the tension head (h) axis (e.g., logbh or h1/b, where b = 2, e, 10, etc.) can have dramatic effects on the shapes and interpretations of soil water release and moisture capacity curves. If the release and capacity data are convex-monotonic (e.g., Figs. 1, 2, and 6–8), curve shapes, inflections and modes obtained by plotting on a nonlinear h axis are pure artefacts of the axis transformation, and have nothing (or very little) to do with actual water release and moisture capacity characteristics. If the data do actually possess an inflection and mode, plotting on a nonlinear h axis will distort curve shapes, as well as change the locations, magnitudes and equations of the inflection and mode coordinates (Fig. 3, eqs. 25–27). A nonlinear h axis may also erase all visual evidence of the data's actual inflection and mode (Fig. 3), and the degree of artefact effect on inflection and mode coordinates may increase with increasing (broader) pore size distribution (Fig. 4). Use of a nonlinear h axis is therefore not recommended when determining inflection points, modes, pore size distributions, soil structure parameters, or soil quality indexes from soil water release and moisture capacity curves.
The vG-I, vG-M, A-G, and D-W models all seem viable and useful, as they were all capable of providing “excellent” or “good” fits (SERN < 11%) to the release curve data-sets, as well as highest or second highest ranked probability (largest or second largest PX) for at least one data-set (Tables 2–5). However, all four models have advantages and implicit limitations that may affect how and when they are applied. The vG-M model (eq. 2) has the advantage of widespread incorporation into most soil water storage and transmission models, while A-G (eq. 11) has the advantage of consistency with capillary rise theory (eq. 1.2). Both models are limited, however, by their implicit assumption that water release curves must be sigmoidal with inflections, and moisture capacity curves must be bell shaped with modes. The vG-I (eq. 2) and D-W (eq. 12) models, on the other hand, have the advantage of being able to simulate convex-monotonic release and capacity curves in addition to sigmoidal release curves and bell-shaped capacity curves, but may be limited by lack of theoretical underpinning. All four models are additionally limited by restricted ability to simulate the angle of release curve intersection with the θ axis; and this may in turn force inaccurate model-data fits at near saturation (e.g., Fig. 6) where soil hydraulic properties and water flow change rapidly with h (Fig. 9). The vG-M and A-G models allow just right-angle intersection with the θ axis (equivalent to dθ/dh = 0 at h = 0), while vG-I and D-W allow right-angle intersection, zero-angle intersection (dθ/dh → −∞), or a fixed acute angle intersection when an exponential release curve is fitted. Release curves from the literature and this study (Figs. 5–8) suggest, however, that inflections and modes are frequently absent, and angles of intersection on the θ axis vary continuously from acute-angle to right angle. It is consequently recommended that more flexible release curve models should be developed that do not assume the existence of inflections or modes, and can also simulate a wide range of realistic intersection angles on the θ axis. In the meantime, it may be advisable to use the vG-I and D-W models over the vG-M and A-G models.
Acknowledgements
We gratefully acknowledge J. Gignac, J. Huffman, and various students for collection and analysis of the soil core and glass bead samples.
Data availability
The data-sets presented in this study are available on request to the corresponding author. The data can be used for scholarly research only, and its source must be acknowledged in writing.
Author contributions
WDR conceived the research and conducted the mathematical and computer analyses; WDR supervised data collection; WDR and CFD contributed to data interpretation and writing of the manuscript.
Funding information
Funding for this work was provided by the Science and Technology Branch of Agriculture and Agri-Food Canada, A-Base Study 2380.
