Title: Simple consistent models for water retention and hydraulic conductivity in the complete moisture range
Abstract: Water Resources ResearchVolume 49, Issue 10 p. 6765-6780 Regular ArticleFree Access Simple consistent models for water retention and hydraulic conductivity in the complete moisture range A. Peters, Corresponding Author A. Peters Institut für Ökologie, Technische Universität Berlin, GermanyCorresponding author: A. Peters, Institut für Ökologie, Technische Universität Berlin, Ernst-Reuter-Platz 1, DE-10587 Berlin, Germany. ([email protected])Search for more papers by this author A. Peters, Corresponding Author A. Peters Institut für Ökologie, Technische Universität Berlin, GermanyCorresponding author: A. Peters, Institut für Ökologie, Technische Universität Berlin, Ernst-Reuter-Platz 1, DE-10587 Berlin, Germany. ([email protected])Search for more papers by this author First published: 07 October 2013 https://doi.org/10.1002/wrcr.20548Citations: 122AboutSectionsPDF ToolsRequest permissionExport citationAdd to favoritesTrack citation ShareShare Give accessShare full text accessShare full-text accessPlease review our Terms and Conditions of Use and check box below to share full-text version of article.I have read and accept the Wiley Online Library Terms and Conditions of UseShareable LinkUse the link below to share a full-text version of this article with your friends and colleagues. Learn more.Copy URL Share a linkShare onFacebookTwitterLinkedInRedditWechat Abstract [1] The commonly used hydraulic models only account for capillary water retention and conductivity. Adsorptive water retention and film conductivity is neglected. This leads to erroneous description of hydraulic properties in the dry range. The few existing models, which account for film conductivity and adsorptive retention are either difficult to use or physically inconsistent. A new set of empirical hydraulic models for an effective description of water dynamics from full saturation to complete dryness is introduced. The models allow a clear partitioning between capillary and adsorptive water retention as well as between capillary and film conductivity. The number of adjustable parameters for the new retention model is not increased compared to the commonly used models, whereas only one extra parameter for quantifying the contribution of film conductivity is required for the new conductivity model. Both models are mathematically simple and thus easy to use in simulation studies. The new liquid conductivity model is coupled with an existing vapor conductivity model to describe conductivity in the complete moisture range. The new models were successfully applied to literature data, which all reach the dry to very dry range and cannot be well described with the classic capillary models. The investigated soils range from pure sands to clay loams. A simulation study with steady-state water transport scenarios shows that neglecting either film or vapor conductivity or both can lead to significant underestimation of water transport at low water contents. Key Points Empirical models distinguish between capillary, adsorptive and film components New models are easy to use in parameter estimation procedures and simulations Neglecting film flow can lead to a significant underestimation of water fluxes 1. Introduction [2] Modeling fluxes of water and solutes in soils are an essential means to address many problems in applied soil science, such as water, nutrient, and salinity management research. Typically, the Richards equation is used in order to model the behavior of water in unsaturated soils. An accurate knowledge of the soil hydraulic functions is required to solve this equation, i.e., the soil water retention function θ(h) and the hydraulic conductivity function K(h), where θ is the volumetric water content, h (L) is the suction, and K (L T−1) is the hydraulic conductivity. This knowledge implies both the appropriate soil hydraulic models and the correct parameterization of these models. [3] The selection of the correct model combination is of crucial importance. In many cases, the well-established retention functions of Brooks and Corey [1964], van Genuchten [1980], or more recently Kosugi [1996] in combination with the capillary bundle models of Mualem [1976a] or Burdine [1953] for conductivity prediction are used and well suited for specific problems. However, a lot of studies show that other model combinations are of great benefit for specific soils or boundary conditions. Soils with structural elements are often better described by bimodal, multimodal [Ross and Smettem, 1993; Durner, 1994] or, for the most flexible description, free from spline retention functions [Iden and Durner, 2007]. [4] Usually, the water retention functions assume a distinct residual water content, θr, which is asymptotically reached at very high suctions. The residual water content is either interpreted as the water held by adsorptive forces [Corey and Brooks, 1999] or as a mere fitting parameter. In the very dry range, even the adsorptive water content finally reaches a value of 0 and thus the concept of residual water content is inappropriate. Several retention models have been proposed to account for this fact [Campbell and Shiozawa, 1992; Fredlund and Xing, 1994; Rossi and Nimmo, 1994; Fayer and Simmons, 1995; Khlosi et al., 2006]. Since the measured water contents are usually based on oven drying, it is conventional to assign a finite suction, at which the water content becomes zero (herein denoted as h0), to a value corresponding to oven dry conditions at 105°C. This yields a suction of ≈ 6.3 × 106 to 107 cm depending on laboratory conditions. [5] The frequently used retention models of Fayer and Simmons [1995] or Khlosi et al. [2006], which account for water adsorption in the medium to dry range, fail for soils with wide pore-size distributions, because in these models θ does not reach 0 at h0. The same applies to the Zhang [2011] retention model. Fredlund and Xing [1994] developed a retention model, where θ equals 0 at h0, regardless of pore-size distribution. Unfortunately, their model does not allow a partition of capillary and adsorptive water. Moreover, Peters et al. [2011] showed that the models of the type introduced by Fayer and Simmons [1995] (including the model of Khlosi et al. [2006]) can lead to the physically unrealistic case of increasing water content with increasing suction. They solved this problem by introducing a slight modification together with appropriate parameter constraints. [6] The capillary bundle models often fail to express hydraulic conductivity in the medium to dry range, because they neglect film and corner flow [Tuller and Or, 2001]. Extensions of the capillary models accounting for film flow can improve hydraulic conductivity prediction [Peters and Durner, 2008a; Lebeau and Konrad, 2010; Zhang, 2011]. In the very dry range, liquid water flow might completely cease, but water can still be conducted by vapor flow [Philip and de Vries, 1957; Saito et al., 2006]. An accurate description of θ(h) and K(h) in the medium to dry range is of great importance for simulating evaporation or root water uptake processes. [7] Peters and Durner [2008a] described liquid conductivity in the complete moisture range with a weighted sum of capillary and film conductivity. In their model, capillary conductivity is described by a capillary bundle model [e.g., Mualem, 1976a] and film conductivity is given by a simple empirical power function of saturation. They showed that their film conductivity model is well suited to describe the measured conductivity data. However, both capillary and film conductivity depend on the capillary retention characteristics with residual water content. Thus, their film conductivity part is coupled with an inappropriate retention model. This conceptual inconsistency was first overcome by Lebeau and Konrad [2010]. They developed a more consistent model, where the retention model of Khlosi et al. [2006] is partitioned into a capillary and an adsorptive part. The capillary and film conductivities are linked to the capillary and adsorptive retention characteristics, respectively. In their model, capillary conductivity is again given by a capillary bundle model, whereas the film conductivity is described by a more complex hydrodynamic model. The same applies to the Zhang [2011] conductivity model. He also partitioned his slightly modified version of the Fayer and Simmons [1995] retention function into capillary and adsorptive water and linked capillary and film conductivity to the according retention parts. [8] The physically based film conductivity in the model of Lebeau and Konrad [2010] is further partitioned into thick and thin film conductivity, accounting for different viscosities in the vicinity of solid surfaces. However, their model is difficult to implement due to its mathematical complexity. Furthermore, the used physical constants might differ from the literature values because of heterogeneities of mineral and organic matter surfaces as well as the usually unknown composition of the fluid. [9] The film conductivity model of Zhang [2011] also uses a number of physical constants, which might differ from the literature values. Moreover, it is coupled to his retention function, which is difficult to use with regard to the determination procedure of the critical point distinguishing between capillary and adsorptive retention. [10] Summarizing, the models accounting for water adsorption and film conductivity are either difficult to use or physically inconsistent. Up to now no physically consistent and easy-to-use hydraulic models for the complete moisture range exist. [11] In this paper, a new retention model is presented allowing a direct partition of the capillary and adsorptive water and guaranteeing that θ = 0 at h0. A new empirical conductivity model is introduced, which links the capillary and film conductivity to capillary and adsorptive water retention. Both models are easy to use. Thus, the gap between simple straightforward-to-use models with physical inconsistencies on the one hand and physically consistent but more cumbersome-to-use models on the other hand is closed. [12] The new liquid conductivity model is coupled with an existing vapor conductivity model to describe conductivity in the complete moisture range. The new models are tested with literature data reaching dry to very dry conditions. Finally, the contributions of capillary, film, and vapor conductivity to total flow are exemplarily shown in a steady-state simulation study. 2. Materials and Methods Theory 2.1.1. New Retention Function 2.1.1.1. Complete Retention Function [13] Total water retention is described by two different terms accounting for water held by capillary and adsorptive forces: (1)where the superscripts tot, cap, and ad mean total, capillary, and adsorptive, θ is the volumetric water content and h (L) is the suction, defined positive for unsaturated soils. This notation will be used throughout the rest of the paper. Saturation (S) can be defined for both fractions. Thus, the new total saturation model is given by the weighted sum of a capillary and an adsorptive saturation term: (2)where w is the weighting factor subjected to w ≤ 1. The soil water retention model is given by , where θs is the saturated water content of the soil. The capillary and adsorptive bound water of equation 1 are given by and (see Figure 1). The capillary bound water is 0 at a total water content of , which might be interpreted as the residual water content of the capillary part. This is the same water content at which the adsorptive fraction is saturated (see Figure 1 for illustration). Figure 1Open in figure viewerPowerPoint Exemplary illustration of the contributions of capillary and adsorptive components to the proposed retention model (equation 7) using the Kosugi function as basic function. Used parameter values are: hm = 20 cm, σ = 0.7, w = 0.8, θs = 0.5, and ha = hm. The capillary held water water θcap was replaced by to emphasize the contributions of capillary and adsorptive retention to total water retention. Note that is identical to the original Kosugi retention function with residual water content, which is given by θs(1−w). [14] Note that saturation is defined as here, where θs can be but does not have to be equivalent to porosity. Capillary saturation can be interpreted as or (see Figure 1), where and resemble the total and residual water contents in the commonly applied retention models with residual water content. Thus, Scap is equivalent to the effective saturation in the publications of e.g., van Genuchten [1980] or Kosugi [1996]. The adsorptive saturation is given by . 2.1.1.2. Basic Functions [15] In the following, the two parts of equation 2 will be explained in detail. The adsorptive water usually decreases linearly toward 0 on a semilog scale [Campbell and Shiozawa, 1992]. This is approximately given by a slight modification of the correction function of Fredlund and Xing [1994]: (3)with [16] where h0 (L) is the suction at water content of 0 and ha (L) is the suction below which X = 1. Xm is a fictitious parameter slightly greater than 1 adopted from the retention model of Vogel et al. [2001] (see Figure 2 for illustration). The adsorptive saturation is simply given by (4) Figure 2Open in figure viewerPowerPoint X as a function of suction (equation 3) with fictitious parameter Xm. [2001]. Parameters were set to ha = 20 cm and h0 = 6.3 × 106 cm. See text for explanation. [17] Since X is also used as correction for the general capillary saturation function (see below), it is distinguished between Sad and X. The introduction of the fictitious parameter Xm and a distinct value for ha, above which Sad is 1, is necessary to guarantee that with increasing suction water of the adsorptive fraction does not decrease before water of the capillary fraction. [18] The basic saturation functions for the capillary part used here are the constrained function of van Genuchten [1980] and the function of Kosugi [1996]: (5)and (6)where α (L−1) and n in equation 5 are curve shape parameters. In equation 6, hm (L) is the suction corresponding to the median pore radius, σ is the standard deviation of the log-transformed pore-size distribution density function and erfc[] is the complementary error function. Since Γ is used as sole capillary saturation in the simple form of the new retention model and as part of the corrected capillary saturation function (see next sections), this study distinguishes between Scap and Γ. 2.1.1.3. Simple Form of New Retention Model [19] Using the basic functions Γ as functions for capillary water retention, i.e., Scap = Γ and Sad = X, the new simplified retention model can be written as: (7) [20] If Γ is expressed by appropriate functions, Mualem's or Burdine's capillary conductivity models have analytical solutions. Setting w = 1 reduces equation 7 to the original saturation function. [21] The parameter ha marks the suction below which saturation of the adsorptive part is 1. Since adsorptive water will not leave the soil before capillary water drains, ha should have a value above the air entry value. One possibility is to treat it as a free fitting parameter. Since the shape of the complete retention model is not very sensitive with regard to ha when fitting the new models to the data used in this study (not shown), it is set at a certain value. Here, ha is expressed in dependence on the basic capillary retention function. The choice is kept simple by defining ha = hm for the Kosugi function. This means that ha is given by h at Γ = 0.5. For the van Genuchten function, ha at Γ = 0.5 is given by . For simplicity, ha is set here at ha = α−1. Parameter α−1 corresponds to a suction where 1 > Γ > 0.5, hence this choice is justified. Note that ha should not be interpreted in a strict physical sense but rather as a shape parameter of the soil water retention function. [22] Setting h0 = 6.3 × 106 cm [see Schneider and Goss, 2012], only four parameters (either θs, w, hm, and σ or θs, w, α and n) are needed to describe the complete retention function. The contributions of capillary and adsorptive water retention to the complete retention characteristics are exemplarily shown in Figure 1. The course of θtot is dominated by θcap in the wet moisture range and by θad in the dry range. The model has a transition zone between these moisture ranges, where both capillary and adsorptive retention are important. This is in accordance to the concept of Lebeau and Konrad [2010] (see their Figure 2). [23] The water content at h = h0 is given by θsΓ(h0). Thus, the same bias as in the retention models of Fayer and Simmons [1995], Khlosi et al. [2006], or Zhang [2011] is allowed. In case of narrow to medium pore-size distributions, Γ reaches approximately 0 at h0, hence . If θ(h0) is significantly greater than 0, which is usually found for soils with wide pore-size distributions combined with a lack of data in the suction range close to h0, a correction is introduced as shown in the next section. 2.1.1.4. New Retention Function for Pore-Size Distributions of Any Width [24] In case of wide pore-size distributions, the saturations of the basic retention functions (equations 5 and 6), and thus the saturation of equation 7 do not reach a value of 0 at h = h0. In fact, this is true for any pore-size distribution, since θ only asymptotically reaches a value of 0 in equations 5 and 6. This problem is solved without extra adjustable parameters by the correction [Fredlund and Xing, 1994], leading to: (8) [25] The benefit of this correction is shown in Figure 3 (top). For small values of n or high values of σ, Γ does not reach a value close to 0 if h → h0. Thus, Scap(h) is not well represented by Γ. If h → h0 then X(h) → 0. Thus, X(h)Γ(h) approaches also 0, regardless of the shape of Γ(h). For narrow pore-size distributions, X(h)Γ(h) ≈ Γ(h). The corresponding complete new saturation models are shown in Figure 3 (bottom). The disadvantage of this capillary saturation model is that no analytical solution for the Mualem conductivity model exists (see below). Figure 3Open in figure viewerPowerPoint Sensitivity of capillary model with respect to pore-size distribution index. (top left) ΓVG and ΓVGX are the van Genuchten capillary saturation function and the corrected function; (top right) ΓKos and ΓKosX are the Kosugi capillary saturation function and the corrected function. (bottom) Complete new saturation functions (including adsorptive saturation) with correction. Numbers indicate values for parameters n and σ. Other parameters were set to w = 0.8, α = 0.005 cm−1, hm = 103 cm, h0 = 6.3 × 106 cm, and ha = α−1 (left) or ha = hm (right). See text for further explanations. [26] A value of 10−3 for θ(h0), corresponding to one-tenth of the assumed measurement error of the water content data (see below), is accepted before switching from the new simple model (equation 7) to the corrected model (equation 8). [27] Note that the parameters α and n or hm and σ have different values in equations 8 and 7. Therefore, they should be interpreted as mere shape parameters without further meaning when equation 8 is used. 2.1.2. New Conductivity Model 2.1.2.1. Complete Model for Liquid Conductivity [28] The hydraulic conductivity function for liquid flow is expressed by the sum of capillary and film conductivity: (9)where the superscripts liq, cap, and film stand for liquid, capillary, and film. The relative hydraulic conductivity function for liquid flow can be expressed by the weighted sum of capillary and film conductivity [Peters and Durner, 2008a], where the capillary and film parts depend on the capillary and adsorptive water saturation functions: (10) [29] The entire hydraulic conductivity function for the liquid phase is: (11)where Ks (L T−1) is the saturated overall conductivity. Equations 9-11 can also be interpreted as follows: (12)where and are capillary and film conductivity at saturation given by and . Note that both terms in equation 10 were coupled with the complete capillary saturation function in the original Peters and Durner [2008a] model, whereas here the concepts of Lebeau and Konrad [2010] and Zhang [2011] are followed by defining them explicitly for the capillary and adsorptive water fractions, respectively. Figure 4 exemplarily shows the contributions of capillary and film conductivity together with isothermal vapor conductivity (see below). These three parts of the complete conductivity model will be explained in the next three sections. Figure 4Open in figure viewerPowerPoint Exemplary illustration of the contribution of capillary, film, and vapor components to the proposed hydraulic conductivity model. Used parameter values are: τ = 0.5, Ks = 100 cm d−1, ω = 10−4, and a = −1.5. The parameters for water retention are the same as in Figure 1. 2.1.2.2. Model for Capillary Conductivity [30] An expression for relative unsaturated hydraulic conductivity for capillary flow is derived from the capillary water retention characteristics by pore-bundle models. The general form of these models can be expressed as follows [Hoffmann-Riem et al., 1999]: (13)where is the relative hydraulic conductivity and x is a dummy variable of integration. The parameters τ, κ, and β can be varied to get more specific functional expressions. For the Burdine model [Burdine, 1953], τ = 2, κ = 2, and β = 1. In Mualem's model [Mualem, 1976a], τ = 0.5, κ = 1, and β = 2, whereas in the model of Alexander and Skaggs [1986] τ = 1, κ = 1, and β = 1. The parameter τ accounts for tortuosity and connectivity in Mualem's original interpretation, hence in a physical sense τ must be positive. However, its physical meaning must be questioned [Hoffmann-Riem et al., 1999], and τ is often treated as a free fitting parameter that is frequently negative [Schaap and Leij, 2000]. Peters et al. [2011] derived boundaries for the lower allowed value of τ, which guarantee physical consistency of the complete function, whereas τ is interpreted as a mere shape parameter. [31] If the new simple retention model (equation 7) is used and Scap is expressed by equations 5 or 6, the analytical solutions of Mualem's capillary-bundle model are given by (14)and (15)where m = 1−1/n and erfc−1[] is the inverse of the complementary error function. If equation 8 is used instead, Mualem's model has to be solved by numerical integration. 2.1.2.3. New Model for Film Conductivity [32] Most of the data in literature show that conductivity as a function of suction decreases more or less linearly on the log-log scale at low water contents (see data below). This is in accordance with the Langmuir-based film flow model for monodisperse particles derived by Tokunaga [2009]. The film conductivity is proportional to the third power of film thickness in his model, leading to a linear relationship between film conductivity and film thickness on the log-log scale. The logarithm of film thickness in turn depends not exactly but approximately linearly on log h (see Figure 3 in Lebeau and Konrad [2010]). Therefore, the hydrodynamically derived models [e.g., Tuller and Or, 2001; Tokunaga, 2009; Lebeau and Konrad, 2010; Zhang, 2011] usually show the exact or approximate linear relationship between film conductivity and log h. Tokunaga [2009] showed that film conductivity is proportional to h−1.5 for constant viscosity, low ionic strength, high surface electrostatic potentials and at high suctions. In other words, the slope on the log-log scale is −1.5. [33] The exact value for film conductivity in real soils depends on many, usually unknown, soil and fluid properties like particle-size distribution, surface roughness, surface charge, ionic strength of the fluid, and so on. Furthermore, soils usually consist of a heterogeneous mixture of different minerals and organic matter of variable sizes and surface properties. Thus, the exact prediction of film conductivity in natural soils with a physically based model like the models of Tokunaga [2009] or Lebeau and Konrad [2010] might need a correction to account for the unknown properties of a certain soil. Zhang [2011] used the derivation of Tokunaga [2009] and introduced a correction factor for saturated film conductivity. Introducing a correction, which accounts for all the different soil and fluid properties, makes the use of different literature constants of surface and fluid properties questionable. Therefore, the linear relationship shall be accounted for but without further physical interpretation. This linearity at suctions greater than ha can be empirically described with a simple power function: (16)where a is the slope on the log-log scale. Following Tokunaga [2009], a is here set at −1.5. Since the dynamic viscosity of the fluid will increase at very high suctions [Or and Tuller, 2000] and the dominant forces acting on the water molecules switch from ionic electrostatic to molecular forces, the slope will slightly change as suction increases and thus thickness decreases [Lebeau and Konrad, 2010]. In this paper, it is assumed that a sufficient approximation is given by an effective slope, which encompasses the complete course of film conductivity. Note that in cases of increasing ionic strength of the fluid or decreasing surface electrostatic potential, the slope will deviate from −1.5. [34] The capillary conductivity is usually given as a function of capillary saturation. Thus, expressing the film conductivity as a function of Sad is more convenient. Solving the unsaturated part of equation 3 with respect to h yields: (17) [35] Since h0/ha >> 1, Xm ≈ 1 and film flow is important when Sad < 1, equation 17 is approximately given by (18)which can be combined with equation 16 yielding: (19)where the slope on the semilog plot with respect to Sad is now given by . The new film conductivity model as a function of either h or Sad is shown in Figure 5 with different values for ω. Figure 5Open in figure viewerPowerPoint Contribution of relative film conductivity in equation 10 (left) as a function of suction and (right) as function of adsorptive water saturation, where is given by equation 19. Numbers indicate values for parameter ω. Other parameters were set to h0 = 6.3 × 106 cm, ha = 20 cm, and a = −1.5. 2.1.2.4. Model for Vapor Conductivity [36] In the very dry range, liquid conductivity will completely cease and water transport is solely governed by vapor flow [Saito et al., 2006]. Thus, when evaporation takes place long enough to dry out the soil surface toward suctions in the surrounding air, vapor transport might be the main process for water flow. The isothermal vapor conductivity Kvap (L T−1) is calculated as described in Appendix Appendix A. Since the gradient in gravitational potential is negligible when vapor transport is the dominant flow process, the complete isothermal conductivity is approximately given by (20)where the superscript lv stands for liquid + vapor. Figure 4 shows the contributions of capillary, film, and vapor conductivity to total conductivity as a function of suction, combining equations 20, 10, 11, 15, 19, and A1. There are three different suction ranges in which either capillary, film, or vapor conductivity predominates. The transition zones are relatively narrow. Note that most available data are in the zone where either capillary or film conductivity is dominant. Therefore, it is not the object of this paper to test the vapor conductivity model. Since the laboratory temperatures for most of the literature data are unknown, the temperature for vapor conductivity is assumed to be 20°C, which seems reasonable for typical laboratory conditions. 2.1.3. Parameter Estimation [37] The models for θ(h) and K(h) were fitted to the data with a nonlinear regression algorithm by minimizing the sum of weighted squared residuals between model prediction and data pairs: (21)where r and k are the number of data pairs for the retention and the conductivity function, wθ and wK are the weights of water content and conductivity data, , , Ki, and are the measured and model predicted values, respectively, and b is the parameter vector. In case of unknown Ks and θs, b consisted of seven adjustable parameters: for the Kosugi and for the van Genuchten function. With known Ks and θs, the size of b was reduced to 5 (see below). [38] The predicted water contents in equation 21 were either calculated in a standard manner as the point water contents at mean suction (“classic method”) or, if the column height was known, as the mean water content of the whole column (“integral method”) to avoid systematic errors [Peters and Durner, 2006]. [39] For normally distributed and uncorrelated measurement errors with zero mean, the single weights can be set to the reciprocal of the variance of the measurement error. This is in accordance with the maximum likelihood principle for the method of least squares [Omlin and Reichert, 1999]. The errors for the retention and the log10(K) data were assumed to be σθ = 0.01 and , leading to wθ = 10000 and wK = 16. [40] A descriptive measure giving the mean deviation between model and data is the root-mean-square error: (22)where yi and are measured and model predicted quantities. For a sound representation of the data by the model, the values of RMSE should