Analytical expressions for noncapillary soil water retention based on popular capillary retention models

Recently, Weber et al. proposed a new model for noncapillary retention of soil water without introducing new parameters. The model is based on the integral of any given saturation function for the capillary part in log space. In the original paper, the integration was carried out numerically. Here, we derive analytical solutions for the most popular soil water retention models (van Genuchten, Kosugi, and Brooks–Corey).


INTRODUCTION
As part of their new modular framework for modeling unsaturated hydraulic properties of soils, Weber, Durner, Streck, and Diamantopoulos (2019) introduced a new model for noncapillary water retention. The commonly adopted approach for modeling the soil hydraulic property functions (Mualem-van Genuchten) is based on completely filled capillaries and is used over the entire range of soil matric potential. However, it has been shown (Diamantopoulos & Durner, 2015;Tuller & Or, 2001) that it fails beyond pF ≈ 2.3, where pF is the common logarithm of soil water suction (absolute value of soil matric potential) given in head units (cm).
The new model has a number of advantages over existing approaches for noncapillary water retention (e.g., Fredlund & Xing, 1994;Peters, 2013). First, it achieves a nearly linear reduction in water content at high, increasing pF; second, a volumetric water content of zero is attained at predefined pF; third, the function is continuously differentiable; and fourth, no extra parameter is introduced. Essen-Abbreviations: WRC, water retention curve.
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited. © 2020 The Authors. Vadose Zone Journal published by Wiley Periodicals, Inc. on behalf of Soil Science Society of America tially, this is achieved by integrating the capillary saturation function in pF (log) space. In Weber et al. (2019), the integration was carried out numerically. The purpose of this technical note is to provide analytical solutions for the integrals of the most popular capillary saturation functions-that is, those published by Brooks and Corey (1964), Kosugi (1996), and van Genuchten (1980. To this end, we briefly review the mathematical structure of the model of Weber et al. (2019).

THEORY
The soil water retention curve (WRC) gives the volumetric soil water content, θ (m 3 m −3 ), as a function of soil water suction, h. The WRC is subdivided into a capillary and a noncapillary part: where θ c and θ nc are the maximum capillary and noncapillary volumetric water contents, respectively, whereas x = log 10 (h). pF c ( ) and pF nc ( ) are the respective saturation functions. These can be derived from the usual capillary saturation models, Γ pF c ( ) (Table 1). Following a proposal by Iden and Durner (2014), the functions Γ pF c ( ) are forced to vary between zero and unity by rescaling them: where Γ 0 is the saturation at which the soil does not dry out further. Weber et al. (2019) proposed to model the noncapillary part of the WRC by integrating pF c ( ) − 1: * pF where x′ is the dummy variable of integration. The integral equals the area between full saturation (x′ = −∞, corresponding to h = 0 cm) and the capillary saturation at a given pF value. It is slightly different from that given by Weber et al. (2019), who replaced the lower bound of the integral by a finite number to make it numerically tractable. They proposed pF = −3 (h = 0.001 cm), where the soil is very close to full saturation, as the lower bound.

Core Ideas
• The commonly used soil water retention functions fall short in the adsorptive (noncapillary) region.
Upon inserting Equation 5 into Equation 4, the factor 1/(1 − Γ 0 ) cancels out, so that Equation 4 can be rewritten without loss of generality as where we define Equation 7 contains the selected capillary saturation model, Γ pF c ( ). The preceding derivation has shown that rescaling Γ pF c ( ) by Equation 2 has no effect on Γ pF nc ( ).

Note that the denominator in Equation 6
: is the total area above Γ pF nc ( ) − 1 between maximum saturation and maximum dryness.

Table 1 lists the analytical solutions of Equation 7
for the most popular capillary soil water retention models. The derivation is described in the Appendix. The noncapillary saturation function pF nc ( ) can be calculated from the given solution by applying Equation 6.
Taking the van Genuchten equation as an example, Figure 1 compares the new analytical solution with the numerical solution applied by Weber et al. (2019) and Weber and Diamantopoulos (2019). The difference between the solutions is negligible.

C O N F L I C T O F I N T E R E S T
The authors declare no conflict of interest.

A C K N O W L E D G M E N T S
We thank Lauritz Streck for the fruitful discussion of some mathematical aspects in this paper. This work was prepared within SFB 1253 CAMPOS ( Substituting x by y = 1 + (α10 x ) n , which implies that d = 1 ln (10) gives The integral in Equation 11 is of the Chebyshev type and can be solved in terms of beta functions. The incomplete beta function is defined by (NIST, 2019, 8.17.1): In case of z = 1 the incomplete beta function turns into Euler's beta function, that is, B(a, b) = B(1; a, b) (NIST, 2019, 5.12.1).
Combining Equations 11 and 12 yields (letting a = 1/n or a = 1, respectively, and b = 0): The application of Equation 13 is not straightforward because B(1/n, 0) and B(1, 0) are not finite. The difference of the two, however, can be evaluated. To this end, we use the series representation of the beta function (retrieved from http://www.wolframalpha.com/with query [series representation of beta function]): where, by definition, (c) k = c(c + 1)(c + 2) . . . (c + k − 1), which is known as the Pochhammer notation. Rewriting the difference of the two beta functions from Equation 13 gives (1) As can be seen from direct comparison with the convergent series ∑ 1∕ 2 (e.g., Greenspan, Benney, & Turner, 1986, Chapter 4.2), the series in Equation 17 converges. Convergence is rather slow but can be considerably accelerated by rewriting the series as where we used the method of telescoping sums (e.g., Greenspan et al., 1986, Chapter 4.6): ; Resubstitution of y = 1 + (α10 x ) n finally yields Equation C1 in Table 1. It is worth noting that even if only the first term of the series is evaluated, the (absolute) error in Γ pF nc ( ) is below 0.01 (for n ∈ [1,10]). Evaluating the first five (ten) terms reduces the maximum error to less than 0.001 (<0.0003).
Alternatively, and from a physical viewpoint sufficient in most cases, the integral in Equation 9 can be evaluated with a finite lower limit ε (see above). In this case, it is easier to apply the substitution only to the left term in the integral Proceeding along the same lines as above yields (24) and, upon resubstitution of y, finally Equation C2 in Table 1.
The solutions presented in Table 1 have been implemented as functions in the R package spsh (https://CRAN.R-project.org/package=spsh; R Core Team, 2019, version 3.6.0). The incomplete beta function in C1 and C2 (Table 1) was computed with the hypergeo package (Hankin, 2016). This required representing it as hypergeometric function (NIST, 2019, 8.17.8).
Equation C2 can be computed faster than Equation C1 but may be difficult to evaluate for large n. Presuming that ε is chosen such that ε + log 10 (α) = −2, which should in general be sufficiently accurate, the computation of C2 demands a routine that can evaluate B[ (ε); 1 , 0] at y = 1 + 10 −2n .