Published online 16 December 2005
Published in Vadose Zone J 5:27-34 (2005)
DOI: 10.2136/vzj2005.0005
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
A Modified Mualem–van Genuchten Formulation for Improved Description of the Hydraulic Conductivity Near Saturation
Marcel G. Schaap* and
Martinus Th. van Genuchten
USDA-ARS, George E. Brown, Jr. Salinity Lab., 450 W. Big Springs Rd., Riverside, CA 92507
* Corresponding author (mschaap{at}ussl.ars.usda.gov)
Received 11 January 2005.
 |
ABSTRACT
|
|---|
The unsaturated soil hydraulic properties are often described using Mualem–van Genuchten (MVG) type analytical functions. Recent studies suggest several shortcomings of these functions near saturation, notably the lack of second-order continuity of the soil water retention function at saturation and the inability of the hydraulic conductivity function to account for macroporosity. We present a modified MVG formulation that improves the description of the hydraulic conductivity near saturation. The modified model introduces a small but constant air-entry pressure (hs) into the water retention curve. Analysis of the UNSODA soil hydraulic database revealed an optimal value of –4 cm for hs, more or less independent of soil texture. The modified model uses a pressure dependent piece-wise linear correction to ensure that deviations between measured and fitted conductivities between pressure heads of 0 and –40 cm were eliminated. A small correction was found necessary between –4 and –40 cm, and a much larger correction was needed between 0 and –4 cm. An average RMSE in logK of only 0.26 remained for a data set of 235 samples. The resulting modified MVG model was found to have small systematic errors across the entire pressure range. The modified model appears well suited for large-scale vadose zone flow and transport simulations, including inverse modeling studies.
Abbreviations: FMVG, Mualem–van Genuchten model with fitted values for Ko and L • MMVG, modified Mualem–van Genuchten • MVVG-A, modified Mualem–van Genuchten refitted with standard retention parameters (hs = 0) • MVG, Mualem–van Genuchten • TMVG, traditional Mualem–van Genuchten model with fixed values for Ko and L
|
INTRODUCTION
|
|---|
MANY VADOSE ZONE flow and transport studies require description of unsaturated soil hydraulic properties over a wide range of pressure heads, h. The hydraulic properties are often described using the pore-size distribution model of Mualem (1976) for the hydraulic conductivity in combination with a water retention function introduced by van Genuchten (1980). The soil water retention equation, 
(h), is given by
 | [1] |
where is the volumetric water content (cm3 cm–3) at pressure head h (cm); r and s are the residual and saturated water contents, respectively (cm3 cm–3); 
(>0, in cm–1) is related to the inverse of the air-entry pressure; n (>1) is a measure of the pore-size distribution (van Genuchten, 1980); and m = 1 – 1/n. The corresponding MVG hydraulic conductivity function, K(h), is
 | [2] |
in which Ko is a matching point at saturation (cm d–1), L is an empirical pore-connectivity parameter, and Se is effective saturation given by
 | [3] |
Accurate measurements of the unsaturated hydraulic conductivity are generally difficult and time-consuming. When the water retention parameters are known, it is therefore common practice to use Eq. [2] to estimate K(Se) or K(h) while making assumptions about the values of Ko and L. Most commonly, a measured saturated hydraulic conductivity (Ks) is used for Ko, while L is fixed at 0.5 (Mualem, 1976). When water retention and unsaturated hydraulic conductivity data are available, an optimization may be performed for several or all of the parameters (Ko, L, r, s, , and n) contained in Eq. [1] and [2] (e.g., van Genuchten et al., 1991). This may be done by optimizing the parameters in Eq. [1] and [2] separately, sequentially, or with a simultaneous procedure in which the objective functions for Eq. [1] and [2] are combined (Yates et al., 1992). Some or all of the parameters in Eq. [1] and [2] are sometimes estimated similarly by means of inverse procedures from measured laboratory and/or field measurements of the pressure head, the water content, and/or fluid fluxes (e.g., Hopmans et al., 2002; Butters and Duchateau, 2002).
While Eq. [1] and [2] have found widespread use, they also have several limitations caused by the particular mathematical properties of Eq. [1] or by the use of default values for Ko and L. Luckner et al. (1989) and Vogel et al. (2000), among others, showed that the retention curve lacks second-order continuity at saturation when n < 2 (i.e., d2/dh2 becomes discontinuous at h = 0). If combined with Eq. [2], the resulting conductivity function, K(h) will then exhibit a sharp drop at pressures just below saturation, especially when n approaches its lower limit of 1.0 (Vogel et al., 2000). For 1.0 < n < 1.3, the reduction in conductivity may become so dramatic that it may cause numerical instabilities in simulations of near-saturated infiltration. Vogel et al. (1985, 2000) solved this problem by forcing the retention function to be constant at s between saturation (h = 0) and some very small negative value of the pressure head, hs. While only minimally affecting the soil water retention curve macroscopically, the modification significantly reduced the nonlinearity in the hydraulic conductivity function near saturation. Vogel et al. (2000) showed that hs values as small as –1 or –2 cm produced much more stable numerical simulations. They obtained substantial improvements in the predicted hydraulic conductivity function of one clay soil when a value of –2 cm was used. However, they did not investigate the optimal value for hs in terms of fitting errors of the retention or the unsaturated hydraulic conductivity properties.
Another limitation of Eq. [2] involves the shape of the hydraulic conductivity function near saturation, especially of structured media (i.e., macroporous soils or unsaturated fractured rock), and how best to estimate the matching point, Ko. van Genuchten and Nielsen (1985) and Luckner et al. (1989) pointed out that the matching point, Ko, ideally should be located at a small negative pressure head to avoid the effects of macropore flow that cannot be captured with Eq. [2]. Schaap and Leij (2000) and Schaap et al. (2001) confirmed that fixing Ko at Ks leads to overpredicted hydraulic conductivity values at most pressure heads. For a range of textures they found that the fitted Ko values were generally about one order of magnitude smaller than the measured Ks values. In addition, Schaap and Leij (2000) found that the fitted L values were often negative, with an optimal value of –1.
Schaap et al. (2001) demonstrated that Eq. [2] with fitted Ko systematically underestimated hydraulic conductivity between 0 and –10 cm. However, as compared to the traditional MVG model (with Ko = Ks and L = 0.5), their approach gave a much better description of the conductivity at more negative pressures. This indicates that Eq. [2] needs to be modified further to describe the unsaturated hydraulic conductivity accurately in both the wet and dry ranges. The main objective of this study hence was to integrate the results found by Schaap and Leij (2000) and Schaap et al. (2001) with the model of Vogel et al. (2000). We investigate the optimal value of hs for four broad soil textural groups in terms of the fitting accuracy for both water retention and the unsaturated hydraulic conductivity. While preserving the advantages of the Vogel et al. (2000) model, the new model has the potential to provide a much better description of the unsaturated hydraulic conductivity near saturation. The performance of our modification (MMVG) will be compared to the traditional Mualem–van Genuchten model with fixed (TMVG) and fitted (FMVG) values for Ko and L.
|
THEORY
|
|---|
Modified Mualem–van Genuchten Model
Vogel et al. (2000) solved the conceptual and numerical limitations of Eq. [1] and [2] by introducing a small minimum capillary height, hs, and a fictitious (extrapolated) parameter m > s in Eq. [1]. Their approach maintains the physical meaning of s as a measurable quantity, while the definition of effective saturation Se is also not affected. The modified water retention equation is given by
 | [4] |
where
 | [5] |
Combination of Eq. [4] with Eq. [2] leads to the modified Mualem–van Genuchten (MMVG) model
 | [6] |
where
 | [7] |
in which
 | [8] |
This model reduces to Eq. [2] for hs = 0. We refer to Vogel et al. (2000) for a detailed discussion of the above equations. They suggested that hs be about –1 or –2 cm, but did not indicate whether other values would lead to better descriptions of observed retention and especially conductivity data. In this study we investigated whether hs can be fixed to a single optimal value applicable to different soil textures. The advantage of having a fixed hs is that it would avoid nonuniqueness problems by not adding another free hydraulic parameter to Ko, L, r, s, , and n, while keeping applications of Eq. [6] as simple as possible. Vogel et al. (2000) showed that the above modification affects the shape of the retention curve only minimally relative to Eq. [1], but that the effects on the unsaturated hydraulic conductivity can be very significant (their Fig. 2 and 3), especially for relatively fine-textured soils when n approaches 1.0.
View larger version (28K):
[in this window]
[in a new window]
|
Fig. 2. Dependence of relative RMSE on values of hs for the four textural groups used in this study: (a) water retention, (b) hydraulic conductivity, and (c) the mean of water retention and conductivity.
|
|
View larger version (8K):
[in this window]
[in a new window]
|
Fig. 3. Average scaled conductivities, R (Eq. [16]), for pressure heads between 0 and –100 cm. The averages were calculated for pressure intervals of 1 cm between 0 and –10 cm, 2 cm between –10 and –50 cm, and intervals of –5 cm between –50 and –100 cm. The piece-wise linear curve indicates the most optimal path through the data between –40 and 0 cm (Eq. [14]).
|
|
Macroporosity Modifications
Schaap and Leij (2000) showed that the fitted values of Ko in Eq. [2] were generally about one order of magnitude smaller than Ks, thus causing an apparent discontinuity in the hydraulic conductivity when Ks would be used for h = 0 and Eq. [2] for h < 0 cm. Conceptually, the discrepancy between Ks and Ko can be explained by the presence of macropores or fractures that dominate the flow regime near saturation and micropores that control matrix flow. Matrix flow would thus be active at all pressures, while macropore flow would be dominant only near saturation and become negligible at some relatively small negative pressure. This conceptualization suggests that Eq. [2] with fitted Ko and L values should be used only for matrix flow. The challenge to be addressed below is how best to account for the effects of macroporosity by further modifying Eq. [6].
Previous attempts to incorporate the effects of macropore flow in the hydraulic conductivity generally focused on composite, dual-porosity type functions in which separate equations are used for the macropore and micropore contributions to the hydraulic conductivity function (e.g., Peters and Klavetter, 1988; Durner, 1994; Mohanty et al., 1997; Köhne et al., 2001). Consistent with this approach, we will use Eq. [6] to describe matrix flow, but adapt this equation for macropore flow near saturation. The modification we propose uses two conductivity terms, Km(h) given by Eq. [6] for matrix flow and a dimensionless pressure-dependent scale factor A(h) to account for the increased conductivity near saturation between a certain pressure head hm and 0. The conductivity function thus becomes
 | [9] |
Ideally, hm would be equal to hs because this would simplify the resulting model. However, no compelling physical reason exists for this since hs was introduced to improve in an empirical way the numerical behavior of Eq. [2] near saturation. We therefore will not assume a priori that hs and hm are equal.
In essence, A(h) introduces a pressure-dependent matching point in Eq. [6]. In this study we assumed that A(h) has the same general functional form for all conductivity data sets, irrespective of soil texture. This avoids the necessity to parameterize A(h) for each conductivity data set separately, which is difficult or impossible given the scarcity of available near-saturated conductivity data. The functional shape of A(h) will be studied by first fitting Eq. [6] to conductivity data using the first-order approximation that hm = hs. Because we explicitly assume that Eq. [6] does not hold near saturation, it is imperative that this equation be fitted only to pressure heads less than hs. Contrary to Eq. [2], pressure heads between hs and 0 cm are thus not used to fit Eq. [6]. We will subsequently derive the shape of A(h) by studying the residuals between the fitted Km(h) and the measured K(h) data. These residuals are expected to exhibit considerable scatter since measured hydraulic conductivity data near saturation are subject to experimental difficulties, while different data sets will lead to different values of Ks, Ko, n, and L in Eq. [6]. To better scrutinize the shape of A(h) we computed scaled residuals according to
 | [10] |
This equation forces R(h) to be between 0 and 1 in the region where Km(h) underestimates measured K(h) data, irrespective of the actual values for Km(h) and Ks. A plot of R vs. h should indicate also whether or not hm can be assumed equal to hs.
|
MATERIALS AND METHODS
|
|---|
The data for this study were taken from the UNSODA database (Leij et al., 1996; Nemes et al., 2001) and are the same as those used by Schaap and Leij (2000). We used 235 laboratory data sets that had at least six –h pairs and at least five K–h pairs. Samples with chaotic data or with limited retention or conductivity ranges were omitted. Figure 1
provides the textural distribution of the samples and their classification in only four textural groups designated as Sands (100 data sets), Loams (41), Silts (58), and Clays (36).
View larger version (34K):
[in this window]
[in a new window]
|
Fig. 1. Textural distribution of samples used in this study and the delineation of the four textural groups.
|
|
The parameters in Eq. [1], [2], [4], and [6] were fitted to water retention and unsaturated hydraulic conductivity data with the simplex or amoeba algorithm (Nelder and Mead, 1965; Press et al., 1988). The objective function for water retention was taken as
 | [11] |
where i and i' are the measured and estimated water contents, respectively, NW is the number of measured water retention data points for each sample, and p is the parameter vector {r, s, , n}. The parameter hs is not fitted but fixed for the entire data set. For optimization of the unsaturated hydraulic conductivity parameters we minimized
 | [12] |
where Ki and Ki' are the measured and estimated hydraulic conductivity, respectively, Nk is the number of measured K(h) data points and p = {Ko, L}. Logarithmic values of Ki were used in Eq. [12] to avoid a bias toward high conductivities in the wet range. For Eq. [6] we used only conductivity data at pressure heads 
hs cm, unless mentioned otherwise. Fitted Ko values were constrained to be smaller than or equal to Ks. In the case of Eq. [2], this produced 62 samples with Ko = Ks (out of 235), while for Eq. [6] the number was only 19. We believe this constraint was hit more often for Eq. [2] because of its shape problems near saturation. We did not carry out simultaneous fits of the water retention and conductivity data.
The goodness of fit was quantified with the RMSE
 | [13] |
where NW,K is the number of water retention or hydraulic conductivity measurements and np is the number of parameters that were optimized. Results will be presented as averages for each textural group as well as for the complete data set. Because logarithmic values were used in the conductivity optimizations, RMSEK results are dimensionless. One RMSEK unit can be interpreted as a factor 10 error in K(h). To more effectively compare errors for different hs values we computed relative RMSE values by dividing RMSEW and RMSEK for hs< 0 with the corresponding RMSE values for hs = 0 (i.e., the MVG model).
To study whether Eq. [2] and [6] underestimate or overestimate observed hydraulic conductivities at particular pressure heads, we computed mean errors for 10 consecutive pressure head ranges between 0, –3.2, –10, –32, –100, –320, –1000, –3200, –10 000, –32 000, and –100 000 cm according to
 | [14] |
where N is the number of observations derived from multiple samples in each pressure head range. As with RMSEK, MEK values are dimensionless because logarithmic values of the conductivities are used.
Parameter Optimization Procedure
The optimal value of hs for water retention and conductivity was studied by comparing fitting errors of Eq. [4] and [6] for hs values between 0 and –20 cm. The most optimal value of hs will be used to compute R given by Eq. [10]. The results should yield information about the functional shape of A(h) and whether or not hm can be assumed equal to hs. We subsequently study whether Eq. [9] (the MMVG model) will provide a better description, in terms of RMSE and ME, than Eq. [2]. The latter equation will be used in two versions, one with fixed parameters (i.e., Ko = Ks and L = 0.5, after Mualem, 1976) and one where Ko and L are optimized (using results from Schaap and Leij, 2000). We also will assess whether the modified model can be simplified by directly using parameters from Eq. [1], rather than parameters from the more complicated Eq. [4].
|
RESULTS AND DISCUSSION
|
|---|
Optimal Value for hs
Figure 2A
shows that relative RMSEW values for water retention are almost constant for –5 < hs < 0 cm. This indicates that the water retention curve within this pressure head range is macroscopically not affected by hs. For hs less than –5 cm, the RMSEW values increased for all soil textural classes, but especially for the Sands. The differences between textural group averages of the fitted r, s, and n parameters in Eq. [1] and [4] (Table 1) are small for an hs value of –4 cm (the value ultimately chosen). Fitted values for Eq. [4] were somewhat larger than those for Eq. [1]. This is probably due to correlation between hs and , since both act as an air-entry pressure. The RMSEW values for Eq. [4] (hs = –4 cm) were virtually the same as those for Eq. [1], indicating that both equations fitted the retention data equally well.
The relative RMSEK values for the hydraulic conductivity (Fig. 2b) decreased sharply between hs = 0 and hs = –1 cm, especially for the Clays and Loams. The reductions are caused by the exclusion of conductivities with h > hs. This approach leads to more realistic Ko and L values and produces better descriptions of the dry parts of the hydraulic conductivity function. The improvement for the Sand group was not expected because this group is normally associated with n values >2 (Table 1). However, further inspection of the database showed that 50 of the 100 samples of the Sand group had n values<2. Thus, the Sand group also benefited from the improved description of the conductivity by using Eq. [6] rather than Eq. [2]. For hs less than –1 cm, the relative RMSEK continued to decrease with lower hs within the entire range between 0 and –20 cm (Fig. 2b).
The RMSEW and RMSEK graphs in Fig. 2a and 2b show different optimal values for hs (i.e., minima in the curves) for the water retention and conductivity data. To arrive at a optimal value for hs we computed the average of the relative RMSE values as (RMSEW + RMSEK)/2. The optimum value of hs is now located between –3 and –5 cm for the Sands, Loams, and Silts, while the Clays have a relatively broad minimum with an optimum value at approximately –8 cm (Fig. 2c). Since the data set as a whole had an optimum at –4 cm, we decided to use this value for the remainder of our study. Use of Laplace's capillary law indicates that this pressure head of –4 cm corresponds with a circular pore of about 0.04-cm radius.
Correction Near Saturation
Figure 3
shows scaled conductivity residuals, R, vs. pressure head, h, for data from all 235 hydraulic conductivity data sets. Values for R are given as averages for 1-cm pressure head intervals between 0 and –10 cm, 2-cm intervals between –10 and –50 cm, and 5-cm intervals for h less than or equal to –50 cm. The graph clearly shows that R decreases from near 1 at h = 0 cm to approximately 0 near h = –40 cm. This indicates that Eq. [6] underestimates observed unsaturated hydraulic conductivities within this range, while below –40 cm no significant under- or overestimation was apparent. The trend in these data suggests a sharp decrease in R near saturation and a slower decrease further away from saturation. A two-element piecewise linear function was used to describe this pattern, with R = 0 at h = –40, a change in slope at h = –4 cm, and R = 1 at h = 0 cm. The change in slope was deliberately located at –4 cm to keep the definition of Eq. [9] as simple as possible, and consistent with the earlier derived value of –4 for hs.
With a least-squares minimization procedure we determined that R at –4 cm should be 0.25 (i.e., the break in the fitted curve of Fig. 3), leading to the following expression for R(h):
 | [15] |
Solving Eq. [10] for K(h) gives
 | [16] |
where Ks is now the measured saturated hydraulic conductivity, R is calculated according to Eq. [15], and Km(h) is calculated using Eq. [6]. Note that this approach does not require an extra parameter, as did the traditional TMVG model. However, in addition to unsaturated hydraulic conductivity data, a measured Ks is now necessary to fit Ko and L with Eq. [16]. Although K(h) is not constant between hs and 0 cm but increases to account for macropore flow effects, Eq. [16] preserves the advantage over Eq. [2] in that it does not have a very steep slope at h = 0 cm. If Km(h) were calculated according to Eq. [2], the numerical problems found by Vogel et al. (2000) would most likely remain.
The need for a two-element piecewise linear correction in Fig. 3 is interesting because it may reflect the contributions of three different sets of pores (0 to –4, –4 to –40, and less than –40 cm) to the hydraulic conductivity, rather than only macropores and micropores. These three pore regions would comprise macropores, mesopores, and micropores (e.g., Luxmoore, 1981), or primary fractures, secondary fractures, and the soil or rock matrix (cf. Simunek et al., 2003). Having three overlapping pore regions would also give credence to the use of multiporosity or multipermeability models for preferential flow in structured media (e.g., Wilson et al., 1992; Gwo et al., 1995; Hutson and Wagenet, 1995), rather than only two regions as is common in most dual-porosity and dual-permeability models (van Genuchten and Sudicky, 1999; Simunek et al., 2003). Still, we note that the deviations between –4 and –40 cm in Fig. 3 could have been caused in part also by mathematical or other limitations in Eq. [4] or [6] not specifically accounted for in our study.
Contrary to Eq. [6], which only holds for hydraulic conductivity data with h hs, Eq. [16] was able to describe all conductivity data, including the near-saturated data. We therefore decided to fit also the Ko and L parameters in the Km term of Eq. [16] directly to all available conductivity data for the entire range of available pressure heads. Table 2 lists fitted parameters and RMSEK values for this model (MMVG) and three other models (TMVG, FMVG, and MVVG-A, modified Mualem–van Genuchten refitted with standard retention parameters [hs = 0] [MMVG-A], the latter to be discussed further below) for all textural groups. We note that the TMVG model has no degrees of freedom whereas the FMVG and MMVG variants have two each (i.e., Ko and L). Fitted log(Ko) values for the MVVG model were considerably lower than those found by Schaap and Leij (2000) for the FMVG model (Eq. [2]) assuming variable Ko and L. The smaller Ko values for MMVG are due to the conceptual split between matrix and macropore flow. Average fitted L values were somewhat larger for the MMVG model than for the FMVG formulation. However, the L values were still negative, which indicates that the new model does not resolve problems with the physical interpretation of the pore interaction term SeL in Eq. [2] and [6] as discussed previously by Kosugi (1999) and Schaap and Leij (2000). The average RMSEK of the FMVG model was about 0.9 order of magnitude lower than the RMSEK of the traditional TMVG model, which uses Ko = Ks and L = 0.5 (Table 2). The RMSEK of the MMVG model is still lower than for the FMVG model (0.261 vs. 0.410 for FMVG), thus indicating that Eq. [16] provides a substantially better description of the unsaturated conductivity than Eq. [2].
View this table:
[in this window]
[in a new window]
|
Table 2. Fitting results for different versions of the Mualem–van Genuchten (MVG) and modified Mualem–van Genuchten (MMVG) models. The results for the traditional Mualem–van Genuchten model with fixed values for Ko and L (TMVG) and Mualem–van Genuchten model with fitted values for Ko and L (FMVG) were derived from Schaap and Leij (2000). The MVVG model includes correction near saturation, and the MVVG-A represents the same model but refitted with standard retention parameters (hs = 0).
|
|
We noted above that the RMSEW and the fitted parameters values obtained with Eq. [1] and [4] differed only marginally. This suggests that the MMVG model can be used directly with retention parameters estimated with Eq. [1] (e.g., using the RETC code; van Genuchten et al., 1991), rather than having to fit the more complicated Eq. [4] to observed retention data. Mean Ko and L parameters for model MMVG-A are listed in Table 2. This model uses Eq. [16] but assumes that all retention parameters (i.e., r, s, , and n) were obtained with Eq. [1]. The differences in averages of fitted parameter values are very small, especially for logKo. Likewise, RMSEK values are only slightly larger for MMVG-A (0.263 vs. 0.261) for all 235 samples. This confirms that Eq. [16] can be used reliably with retention parameters estimated with Eq. [1]. However, it is still necessary to use Eq. [5] through [8] and Eq. [15] and [16] to compute the unsaturated hydraulic conductivity function.
Systematic Deviations vs. Pressure
Mean errors obtained with the TMVG, FMVG, and MMVG models for nine pressure head classes are shown in Fig. 4
. Negative ME values indicate underestimation of the conductivity, while the gray bars denote the number of samples in each pressure head class. Results for MMVG-A are not shown since they were indistinguishable from those for MMVG. Notice that the TMVG model overestimates conductivities between 0 and –100 cm for the Sands, between 0 and –10 000 cm for the Silts, and between 0 and –3000 cm for the Clays. For the Loams, the TMVG model describes conductivity well until a pressure of –100 cm. Conductivities of the Sands and Loams in particular are severely underestimated with the TMVG model at more negative pressures. The FMVG model further underpredicts the conductivity between saturation and –30 cm, thus confirming the findings shown in Fig. 3. Below –30 cm the mean errors are nearly zero. The Sand and the Silts show a small underestimation at –10 000 cm with the FMVG and MMVG models, which likely is a result of the small number of samples in this pressure range (gray bars). Because of the correction inherent in Eq. [15] and [16], the MMVG model as expected performed better than the FMVG model near saturation. A slight underestimation exists for the Sands and Silts at pressures at about –3 cm. For the Loams, the MMVG model overestimated the conductivity between 0 and –10 cm. For all samples, the MMVG model performed much better than either the TMVG model or the FMVG models, thus indicating that this model is able to describe unsaturated conductivities correctly for the entire range of pressure heads between 0 and –10 000 cm.
View larger version (28K):
[in this window]
[in a new window]
|
Fig. 4. Mean errors (ME) for nine pressure head classes for the soil textural groups. The number of observations in each pressure class is indicated with the gray bars (right axis). Note that the scale on the right axis is not the same for all graphs.
|
|
|
SUMMARY AND CONCLUSIONS
|
|---|
A Mualem–van Genuchten model with improved description of conductivity near saturation was presented. As previously proposed by Vogel et al. (1985, 2000), the model introduces a small but non-zero air-entry pressure, hs, in the water retention curve. The model was further modified to account for hydraulic conductivity matching points (Ko) that are much smaller than the measured saturated conductivities (Ks). In this study we found that the optimal value for hs was –4 cm. We also found that the modified model still underpredicted K(h) between 0 and –40 cm. We were able to correct for this underestimation by introducing a piecewise linear function that accounts for the effects of macropores and fractures near saturation. A small correction was needed between –4 and –40 cm (as caused by mesopores), while a stronger correction was necessary between 0 and –4 cm (as caused by macropores). After refitting the Ko and L parameters, an average RMSE of 0.261 remained for the database containing 235 samples. This error is lower than when the traditional Mualem-van Genuchten model was fitted (0.410) with variably Ko and L (model FMVG), and much lower than when default parameters were assumed for the TMVG model (1.309). A plot of mean errors vs. pressure head showed that the MMVG has small systematic deviations from measured hydraulic conductivities. The modifications were found to be beneficial for all four textural groups used in this study.
It is attractive to interpret the modified model in terms of matrix and macropore flow, possibly augmented with flow through intermediate sized mesopores. However, we emphasize here that the corrections were completely empirical and did not rely on statistical pore size distribution concepts such as those inherent in the model of Mualem (1976). The modified model should be used with some caution since hydraulic conductivity measurements near saturation are generally quite unreliable. Other mathematical and conceptual limitations in the Mualem–van Genuchten equations may well have found their way into the derived conductivity corrections as well. Still, we believe that the new MMVG formulation provides a substantial improvement in the macroscopic description of the hydraulic conductivity function. The model seems to be especially suited for large-scale studies that require realistic simulations of saturated or near-saturated infiltration into soils. Since both the retention and conductivity data could be described well, the model may also be well suited for inverse modeling studies. However, because the model in its present form does not explicitly define separate matrix and macropore domains, it may not be well suited for simulations of preferential flow.
|
ACKNOWLEDGMENTS
|
|---|
The authors gratefully acknowledge the support by NSF and NASA (EAR-9804902, EAR-0440024), and the ARO (39153–EV). This paper is based in part also on work supported by SAHRA (Sustainability of semi-Arid Hydrology and Riparian Areas) under the STC Program of the National Science Foundation, Agreement no. EAR-9876800.
|
REFERENCES
|
|---|
- Butters, G.L., and P. Duchateau. 2002. Continuous flow method for rapid measurement of soil hydraulic properties. I. Experimental considerations. Available at www.vadosezonejournal.org. Vadose Zone J. 1:239–251.[Abstract/Free Full Text]
- Durner, W. 1994. Hydraulic conductivity estimation for soils with heterogeneous pore structure. Water Resour. Res. 30:211–223.[CrossRef]
- Gwo, J.P., P.M. Jardine, G.V. Wilson, and G.T. Yeh. 1995. A multiple-pore-region concept to modeling mass transfer in subsurface media. J. Hydrol. (Amsterdam) 164:217–237.
- Hopmans, J.W., J. Simunek, N. Romano, and W. Durner. 2002. Inverse methods. p. 963–1008. In J.J. Dane and G.C. Topp (ed.) Methods of soil analysis. Part 4. SSSA Book Ser. 5. SSSA, Madison, WI.
- Hutson, J.L., and R.J. Wagenet. 1995. A multiregion model describing water flow and solute transport in heterogeneous soils. Soil Sci. Soc. Am. J. 59:743–751.[Abstract/Free Full Text]
- Köhne, J.M., S. Köhne, and H.H. Gerke. 2001. Estimating the hydraulic functions of dual-permeability models from bulk soil data. Water Resour. Res. 38(7)1121. doi:10.1029/2001WR000492.
- Kosugi, K. 1999. General model for unsaturated hydraulic conductivity for soils with lognormal pore-size distribution. Soil Sci. Soc. Am. J. 63:270–277.[Abstract/Free Full Text]
- Leij, F.J., W.J. Alves, M.Th. van Genuchten, and J.R. Williams. 1996. The UNSODA Unsaturated Soil Hydraulic Database; User's manual, Version 1.0. EPA/600/R-96/095, National Risk Management Laboratory, Office of Research and Development, USEPA, Cincinnati, OH..
- Luckner, L., M.Th. van Genuchten, and D.R. Nielsen. 1989. A consistent set of parametric models for the two-phase flow of immiscible fluids in the subsurface. Water Resour. Res. 25:2187–2193.
- Luxmoore, R.J. 1981. Micro-, meso-, and macroporosity of soil. Soil Sci. Soc. Am. J. 45:671–672.[Free Full Text]
- Mohanty, B.P., R.S. Bowman, J.M.H. Hendrickx, and M.Th. van Genuchten. 1997. New piecewise-continuous hydraulic functions for modeling preferential flow in an intermittent flood-irrigated field. Water Resour. Res. 33:2049–2063.[CrossRef]
- Mualem, Y. 1976. A new model predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.[CrossRef]
- Nelder, J.A., and R. Mead. 1965. A simplex method for function minimization. Comput. J. 7:308–313.
- Nemes, A., M.G. Schaap, F.J. Leij, and J.H.M. Wösten. 2001. Description of the unsaturated soil hydraulic database UNSODA version 2.0. J. Hydrol. (Amsterdam) 251:151–162.
- Peters, R.R., and E.A. Klavetter. 1988. A continuum model for water movement in an unsaturated fractured rock mass. Water Resour. Res. 24:416–430.
- Press, W.H., B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling. 1988. Numerical recipes in C. 1st ed. Cambridge Univ. Press, New York.
- Schaap, M.G., and F.J. Leij. 2000. Improved prediction of unsaturated hydraulic conductivity with the Mualem-van Genuchten model. Soil Sci. Soc. Am. J. 64:843–851.[Abstract/Free Full Text]
- Schaap, M.G., F.J. Leij, and M.Th. Van Genuchten. 2001. Rosetta: A computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. J. Hydrol. (Amsterdam) 251:163–176.
- Simunek, J., N.J. Jarvis, M.Th. van Genuchten, and A. Gärdenäs. 2003. Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol. (Amsterdam) 272(1–4):14–35.
- van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.[Abstract/Free Full Text]
- van Genuchten, M.Th., F.J. Leij, and S.R. Yates. 1991. The RETC code for quantifying the hydraulic functions of unsaturated soils. Rep. EPA/600/2–91/065. R.S. Kerr Environmental Research Laboratory, USEPA, Ada, OK.
- van Genuchten, M.Th., and D.R. Nielsen. 1985. On describing and predicting the hydraulic properties of unsaturated soils. Ann. Geophys. 3:615–628.
- van Genuchten, M.Th., and E.A. Sudicky. 1999. Recent advances in vadose zone flow and transport modeling. p. 155–193. In M.B. Parlange and J.W. Hopmans (ed.) Vadose zone hydrology: Cutting across disciplines. Oxford Univ. Press, New York.
- Vogel, T., M. Cislerova, and M. Sir. 1985. Reliability of indirect estimation of soil hydraulic conductivity. Vodohosp Cas 33(2):204–224.
- Vogel, T., M.Th. van Genuchten, and M. Cislerova. 2000. Effect of the shape of the soil hydraulic functions near saturation on variably-saturated flow predictions. Adv. Water Resour. 24(2):133–144.[CrossRef]
- Wilson, G.V., P.M. Jardine, and J.P. Gwo. 1992. Modeling the hydraulic properties of a multiregion soil. Soil Sci. Soc. Am. J. 56:1731–1737.[Abstract/Free Full Text]
- Yates, S.R., M.Th. van Genuchten, A.W. Warrick, and F.J. Leij. 1992. Analysis of measured, predicted, and estimated hydraulic conductivity using the RETC computer program. Soil Sci. Soc. Am. J. 56:347–354.[Abstract/Free Full Text]
This article has been cited by other articles:
|
|
|
|
 
N. K. C. Twarakavi, H. Saito, J. Simunek, and M. Th. van Genuchten
A New Approach to Estimate Soil Hydraulic Parameters Using Only Soil Water Retention Data
Soil Sci. Soc. Am. J.,
February 15, 2008;
72(2):
471 - 479.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
Z. F. Zhang, M. Oostrom, and A.L. Ward
Saturation-Dependent Hydraulic Conductivity Anisotropy for Multifluid Systems in Porous Media
Vadose Zone J.,
November 20, 2007;
6(4):
925 - 934.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
