VZJ Journal of Natural Resources and Life Sciences Education
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (2)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Chan, T. P.
Right arrow Articles by Govindaraju, R. S.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Chan, T. P.
Right arrow Articles by Govindaraju, R. S.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Chan, T. P.
Right arrow Articles by Govindaraju, R. S.
Related Collections
Right arrow Soil Hydrology
Right arrow Vadose Zone Processes and Chemical Transport
Published in Vadose Zone Journal 3:1443-1454 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA

ORIGINAL RESEARCH

Estimating Soil Water Retention Curve from Particle-Size Distribution Data Based on Polydisperse Sphere Systems

T. P. Chan and R. S. Govindaraju*

School of Civil Engineering, Purdue University, 500, Stadium Mall Drive, West Lafayette, IN 47907-2501
* Corresponding author (govind{at}ecn.purdue.edu)

Received 28 January 2004.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 MODEL FORMULATION
 COMPARISON WITH SOIL DATA
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Soil hydraulic properties are necessary for understanding water and solute movement in subsurface vadose zone environments. Previous statistical models of soil hydraulic properties have postulated a pore-size distribution from the particle-size distribution by assuming an empirical relationship between pore sizes and particle sizes. This pore-size distribution has been used to develop the soil water retention curve and/or the relative conductivity curve. In this study, we develop soil water retention models from particle-size distribution data based on the results obtained by Torquato and coworkers on polydisperse particle systems. We model soils as systems of randomly placed spheres of random sizes and derive expressions for the soil water retention curve based on the statistical description of the systems. Two special cases of homogeneous systems of polydisperse spheres are considered: fully penetrable spheres (FPS) and totally impenetrable spheres (TIS). A lognormal particle-size distribution is assumed. More than 100 soils are used to corroborate the models. It is found that the TIS model works well for sand and loamy sand soils.

Abbreviations: cdf, cumulative distribution function • cmf, cumulative mass fraction • FPS, fully penetrable spheres • pdf, probability density function • TIS, totally impenetrable spheres • UNSODA, Unsaturated Soil Hydraulic Database


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
MODEL FORMULATION
 COMPARISON WITH SOIL DATA
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
ONE OF THE FUNDAMENTAL soil hydraulic properties is the water retention characteristic curve that relates capillary pressure to the degree of water saturation in a particular soil. This relationship must be known to model or predict water and solute movement in the subsurface vadose zone where the soil is not fully saturated. Many models of relative hydraulic conductivity (see Mualem, 1986) also rely on the knowledge of the soil water retention curve in their formulation. Unfortunately, it is generally very time-consuming to measure this unsaturated soil property, although progress has been made to develop more efficient methods and apparatus (see Salehzadeh and Demond, 1994; Dane and Topp, 2002, p. 675–720). Consequently, much effort has been directed at relating the water retention relationship to some structural soil properties such as particle-size distribution data (which can be obtained through sieve analysis and/or hydrometry). In this study, we contribute to this ongoing effort by developing a statistically based soil water retention model from soil particle-size distribution based on many-particles systems.

Most existing models of water retention curve are empirical curve-fitting equations that contain two or more fitting parameters. The more widely used ones are the Brooks and Corey (1964), van Genuchten (1980), and Russo (1988) models. Of the three models, the equation developed by van Genuchten provides the most flexible form that can fit very well to a wide range of soils. Although these models offer simple functional expressions that can be easily incorporated into the modeling of unsaturated flow, they lack a sound theoretical basis that allows for physical understanding of the fitting parameters. As a result, alternate approaches have been pursued by many researchers.

Kosugi (1994) developed a three-parameter model using a lognormal pore-size distribution. Kosugi argued that since the particle-size distribution of many soils is approximately lognormal (Shirazi and Boersma, 1984; Buchan, 1989), it is reasonable to assume the same for their pore-size distribution. It was found that this model performed as well as the van Genuchten model and other existing empirical ones. Although Kosugi emphasized physical interpretation of the parameters, this model is nonetheless a purely empirical curve-fitting equation.

There is a need to relate easily measurable physical properties of a soil to its water retention characteristics to lend a physical basis to the model (Perrier et al., 1996). Only if such a model is developed can parameters reflect their true physical significance. One way to bridge this gap between soil data and hydraulic properties is through the use of pedotransfer functions (Wösten et al., 2001). Parameters of well-known empirical water retention functions, such as the van Genuchten equation, are correlated using multiple-linear regression (e.g., Vereecken et al., 1989; Wösten et al., 1999) or artificial neural networks (Schaap et al., 1998) with various soil properties, such as texture, structure, organic content, and bulk density. Wösten et al. (2001) provided a comprehensive review of the pedotransfer function approach, and Cornelis et al. (2001) compared several pedotransfer water retention functions and found that the one proposed by Vereecken et al. (1989) gave the best performance. The pedotransfer function approach seeks to find the relationship between fitting parameters of retention models and soil properties, but it ignores the importance of the physics that governs the drainage and wetting of the soil.

A different approach to mending the missing link between soil properties and water retention curves is to make use of the fact that the water retention curve is governed by the spatial arrangement of soil grains. As a result, the retention function can be related directly to the particle-size distribution of the soil. Arya and Paris (1981) attempted to do so by developing a model that uses particle-size distribution of a soil and its bulk density data to predict the water retention curve. They proposed a linear relationship between the mean pore radius and the mean particle radius with a proportionality constant being a function of some measurable parameters of the soil and an empirical scaling constant to account for tortuosity of the pores. Different methods to estimate this scaling factor were proposed and introduced by Arya et al. (1999). Haverkamp and Parlange (1986) also assumed a linear relationship between the particle and pore-size distribution in their analysis. They adopted the van Genuchten expression for fitting the particle-size distribution. By integrating with the Brooks and Corey relationship, their new model yielded good results when compared with retention data of sandy soils that do not contain organic matter. Assouline et al. (1998) proposed a power-law relationship between the normalized particle volume and the normalized pore volume. Using the Weibull distribution to model the particle volume distribution (as a result of a random fragmentation process), a two-parameter water retention expression was developed and was fitted to soil data. Although this model performs well with various types of soils, it has been reduced to an empirical curve-fitting expression. More recently, Zhuang et al. (2001) employed the nonsimilar media concept proposed by Miyazaki (1996) to estimate the water retention curve from particle-size distribution and bulk density. Essentially, they assumed that the mean particle diameter for each size fraction is proportional to the mean pore diameter with a function that involves a shape factor. The shape factor is also related to the mean particle diameter using a logistic function. Their model applied well to the soils tested, except for sandy soils, and the results were consistent with Arya and Paris' model.

Tyler and Wheatcraft (1989) introduced the fractal concept to reinterpret the relationship between the particle- and pore-size distributions based on analytical power functions such as the Brooks and Corey (1964) equation. Fractal water retention models were subsequently proposed by Tyler and Wheatcraft (1990) (using Sierpinski carpet to model a fractal pore space), Rieu and Sposito (1991a)( 1991b) (based on fragmentation of aggregated soil), and Hunt and Gee (2002) (derived from a fractal particle-size distribution). These results were further generalized by the model proposed by Perrier et al. (1996) under the assumption that the pore-size distribution is fractal in nature. The critical parameter in these models is the fractal dimension D that can be related to the particle-size distribution (Rieu and Sposito, 1991b; Perrier et al., 1996; Hunt and Gee, 2002).

As seen from the brief review above, previous statistical models of soil water retention have obtained an empirical relationship of the pore-size distribution from the particle-size distribution by assuming a linear or power relationship between pore sizes and particle sizes. An attempt to obtain a statistically based relationship was made by Chan and Govindaraju (2003). They derived the pore-size distribution based on a system of overlapping spheres with uniformly distributed sizes. The model parameters, as a result, are related to the soil particle-size distribution; however, the soil particle-size distribution data were not directly utilized in fitting the model. The relationship between the model parameters and the soil particle-size distribution was not pursued.

In this study, we developed a new method to infer the pore-size distribution and further develop the soil water retention curve from particle-size distribution of soils based on the results obtained by Torquato (2002) on polydisperse sphere systems. The new models differ from the earlier model of Chan and Govindaraju (2003) in that they incorporate particle distribution information and porosity of the soil to predict the water retention curve. Moreover, the system of impenetrable spheres was also considered in addition to the overlapping spheres. More than 100 soils were selected from a soil database for corroboration of the proposed models.


   MODEL FORMULATION
TOP
 ABSTRACT
 INTRODUCTION
 MODEL FORMULATION
COMPARISON WITH SOIL DATA
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
A soil is a porous medium that can be conceptualized as a mixture of particles of different sizes and shapes. The packing or arrangement of these particles can be vastly different from soil to soil, and so can be their spatial distribution. Traditionally, characterization of the soil structure has been all but impossible because of this random heterogeneity. However, Torquato (2002) recently extended the formalism of statistical mechanics to quantify the microstructure of many-particle systems, thus allowing the better understanding of, and perhaps explicit definition of, the linkage between microscopic properties of these systems and their macroscopic properties. In this study, we apply some of these results in the context of natural soils. Of particular interest to us are the results obtained by Lu and Torquato (1991)(1992) on the so-called polydisperse particle systems. These systems consist of sphere of many sizes, randomly placed in space. We believe some soil textures, such as sand and loamy sand that are granular in nature, can be effectively modeled based on this conceptualization. In the following discussion we present only those results that help us develop a water retention curve from the statistical description of polydisperse spheres. Readers are referred to the literature references provided in Torquato (2002) for detailed derivations.

Two special cases of homogeneous systems of polydisperse spheres are considered here: fully penetrable spheres and totally impenetrable spheres (see Fig. 1) . In a fully penetrable sphere system, spheres are allowed to overlap; thus, their placement in space follows a Poisson distribution. In the second case, the system consists of spheres that are nonoverlapping, or impenetrable by each other. These two special cases will yield different pore-size distributions and water retention curves, as we will see in the following discussion.



View larger version (23K):
[in this window]
[in a new window]
  		Fig. 1. Schematic of the polydisperse sphere systems in two dimensions: (left) impenetrable hard spheres and (right) overlapping spheres.

 
Particle-Size Distribution
Polydisperse sphere systems consist of spheres of random sizes. Various distributional forms, such as lognormal, Weibull, and Schulz, can be used for the sphere-size distribution. In this study, we assume that sphere sizes (radii), R, are lognormally distributed, since many advantages of a lognormal distribution have been identified by Buchan (1989), and several studies have adopted this assumption for soils (e.g., Shirazi and Boersma, 1984; Buchan et al., 1993). However, it should noted here that in most of the cases, the lognormal distribution was applied as a distribution of particle sizes in terms of mass as opposed to, in this study, the particle number. Let y = lnR, then y is normally distributed. The probability density function (pdf) for the sphere sizes is therefore

[1]
where µy is the mean and {sigma}y is the standard deviation of lnR. The nth moment, mn, of this distribution is given by

[2]
where < > denotes expectation. Equation [2] can be generalized as (Govindaraju et al., 2001)

[3]
where erfc = {int}{infty}texpd{xi} is the complementary error function and {xi} is a dummy integration variable. The cumulative distribution function (cdf) of the lognormal particle sizes is given by Eq. [3] when n = 0.

For use in subsequent derivation, we define a dimensionless reduced density, {eta}, as

[4]
where {rho} is the total number density (number of spheres per unit volume) and <v(R)> is the average sphere volume given as 4{pi}m3/3. The average sphere surface area is given as

[5]

Fully Penetrable Spheres
For a statistically homogeneous fully penetrable sphere system, the porosity, {phi}, is given by its 1-point probability function (Lu and Torquato, 1991) for the void phase:

[6]

The specific surface, s, is given as

[7]

Torquato (2002) defined the pore-size probability density function, p({delta}), for isotropic media as follows:

[8]

In other words, p({delta})d{delta} is the probability that one can insert a "test" sphere of radius {delta} at a randomly chosen point in the void space. Note that at {delta} = 0, the pore-size probability density is s/{phi}, the interfacial area per unit pore volume, and when {delta} approaches infinity, the probability density is zero. As pointed out by Torquato, the pore-size distribution derived from this definition is different from the "effective" pore-size distribution obtained by a mercury intrusion experiment. We will discuss this problem in detail below and propose a scaling ratio to derive the "effective" pore-size distribution based on Torquato's definition.

Without the scaling factor, the pdf and the cdf of the pore-size distribution are, respectively,

[9]

[10]
where hV(·) and eV(·) are the void nearest-surface pdf and its complementary cdf, respectively. The definition of hV(·) is the same as that of p({delta}) with the exception that hV(·) is not conditioned on being in the void space, but rather any arbitrary point in the system. So if r is the distance from a reference point in the system to the nearest particle surface, a negative value of r indicates that the reference point lies within the particle phase. Lu and Torquato (1992) found the expressions for eV(r) and hV(r) as follows:

[11]

[12]
where H(r + R) is the Heaviside step function and R is the size of the reference solid particle. Although r theoretically lies within the interval (–{infty},{infty}), it has to be larger than –R from physical considerations. Substituting Eq. [2] into Eq. [11] gives

[13]
and for [12]

[14]

The expressions for r < 0 are omitted here because they are not required for Eq. [9] and [10].

Totally Impenetrable Spheres
For a statistically homogeneous totally impenetrable sphere system, the reduced density is the volume fraction of the particle phase; therefore, the porosity is

[15]

The pore-size distribution for the totally impenetrable spheres follows the same definition given in Eq. [8]. Equations [9] and [10] are valid as the pore-size pdf and cdf, respectively.

Lu and Torquato (1992) derived approximate expressions for the void nearest-surface complementary cdf and pdf for r ≥ 0 by invoking the Carnahan–Starling approximation:

[16]

[17]
where S is the surface area ratio given by

[18]
and the coefficients a0, a1, and a2 are as follows:

[19]

The expressions for r < 0 are again not given since they are not required to calculate the pore-size distribution functions.

Effective Pore-Size Distribution
As we mentioned above, the pore-size distribution given by Eq. [9] and [10] is not the usual pore-size distribution obtained experimentally by mercury porosimetry. Torquato's definition of pore size is biased toward smaller values because the pore radius is defined as the distance from a randomly selected point in the void space to the nearest pore–solid interface. The use of "nearest" quantities results in this bias. Taking a spherical pore of radius rp as an example, for a point picked at random within the spherical pore, {delta} is the distance from the point to the pore perimeter (equivalent to the nearest solid surface), and the expected value of {delta} is found to be

[20]

Equation [20] shows that the effective pore radius, rp, is related to the expected value of {delta} by a factor of 4.

The relationship between {delta} and rp is given above for a single pore. To upscale this relationship to an ensemble of pores prescribed by the many-particle systems, one must know the joint distribution of {delta} and rp. However, the derivation of this joint distribution is not trivial. One major obstacle lies in the complex geometry of the pore structure. It is, therefore, very difficult (subjective) to define what constitutes a pore or pore-throat (see Dullien, 1992). Chan and Govindaraju (2003) adopted a definition of pore size based on the smallest separation between two spheres given that they do not overlap—this definition is also biased toward smaller pore sizes. They did not propose a way to deal with this bias but argued that their definition provides a measure of pore throats that essentially governs the drainage of a porous medium.

Here we propose a linear relationship between {delta} and rp; that is, rp = {alpha}{delta}, where {alpha} is a scaling factor obtained by equating the first moment shown in Eq. [20] (i.e., {alpha} = 4). Then, the pore-size pdf and cdf can be rewritten, respectively, as

[21]

[22]
where the subscript "e" denotes the effective pore-size distributions since the scaling factor is used.

Soil Water Retention Curve
The effective soil water saturation, {Theta}, is defined as

[23]
where {theta}s and {theta}r are the saturated and residual volumetric water contents, respectively. The effective pore-size pdf was related to this effective saturation by (Kosugi, 1994; Chan and Govindaraju, 2003):

[24]

Equation [24] assumes that all pores described by the pore-size distribution are accessible for drainage. Integrating Eq. [24], we obtain the effective water saturation that is equivalent to the pore-size cdf given by Eq. [10]:

[25]

On the basis of Eq. [25], the effective soil water saturation is interpreted as the volumetric fraction of water saturated pores that have a diameter ≤rp. We can further relate the capillary pressure head, h, to the pore radius for a given saturation by using the Young–Laplace equation:

[26]
where {sigma} is the interfacial tension between air and water, {psi} the contact angle, {rho} the density of water, and g the gravitational acceleration. At room temperature and for a contact angle of 0, c = 0.149 cm2. Transforming rp in Eq. [25] into h using Eq. [26], we have

[27]

This yields the main equation for the soil water retention models. The void nearest-surface function, eV(·), is given by Eq. [13] and [16] for fully penetrable spheres and totally impenetrable spheres, respectively.


   COMPARISON WITH SOIL DATA
TOP
 ABSTRACT
 INTRODUCTION
 MODEL FORMULATION
 COMPARISON WITH SOIL DATA
RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
To verify the proposed models, we need soils that have data not only for the complete range of the water retention curve but also the particle-size distribution. The Unsaturated Soil Hydraulic Database (UNSODA) (Nemes et al., 2001) was the source of the soil data used in this study. One criterion for soil selection was the number of data points available for the particle-size and retention data obtained in the laboratory. The number of data points required for inclusion in this study was set to ≥7 for both the particle-size and the retention data, so a comprehensive data set can be obtained while enough data points are ensured to provide meaningful fits to the particle-size distribution and retention curve. Another criterion is related to the soil texture. The TIS model consists of individual spheres that represent granular particles; therefore, it can best describe sandy soils. The FPS model, on the other hand, allows spheres to overlap, and intuitively one may think that it can be used to model soils with aggregation. However, soil particles in aggregates are not exactly overlapping each other, rather they conglomerate, or stick together. So we believe that the FPS model will be inferior to the TIS model if the soil particle-size data are applied equivalently to both models (without the correction for overlapping spheres in the FPS case). On the basis of these considerations, we limited our soil selection to sands and loamy sands, in which aggregation is less likely to occur. Using these criteria, an initial pool of 137 soils was extracted from the database. Out of the 137 soils, 18 were eliminated because of either or both of the following reasons:

  1. The water retention data did not span the whole range of saturation.
  2. Data were deemed unreliable because of outliers or errant data points.

The particle-size data provided by UNSODA are given in terms of cumulative mass fraction; however, the lognormal distribution (Eq. [1]) assumed for the models describes the distribution not of the sphere mass but the sphere number. The relationship in Eq. [3] turns out to be very useful because it can be used to relate the cumulative mass fraction of particle sizes to Eq. [1]. Assuming soil particles are spheres [v(R) = 4{pi}R3/3], the cumulative mass fraction can be written as

[28]
where {rho}s is the particle density. Therefore, Eq. [28] is used to fit the particle-diameter–mass fraction data pairs of each soil to obtain estimates of µy and {sigma}y. A nonlinear least-squares curve-fitting routine in MATLAB (Mathworks, Natick, MA) was used to perform all curve-fitting analysis in this study. To obtain the retention curve for a soil, Eq. [27] and [23] are used while the void nearest-surface function, eV(·), is given by Eq. [13] and [16] for the FPS and TIS models, respectively. The retention equations require the following parameters: {phi}, {theta}s, {theta}r, and {alpha} besides µy and {sigma}y. For most of the selected soils, neither the porosity, {phi}, nor the saturated water content, {theta}s, is provided. Also, our models do not differentiate between the connected and unconnected pore space in estimating the pore-size distribution and therefore treat all void space as being fully accessible to the displacement of the fluid. Thus, in this study, we assume that the porosity is equivalent to the saturated water content in predicting the retention curve. For soils that do not provide a porosity nor a saturated water content, the first point of the retention data that corresponds to the lowest capillary pressure and highest saturation are used as {theta}s. For cases where the porosity or saturated water content is given but is larger than the first point of the curve, {theta}s is allowed to be fitted within the interval bounded by the two values. The residual water content, {theta}r, is also a fitting parameter. The scaling factor, {alpha}, is prescribed using the method of moments described previously; however, one can regard it as a fitting parameter, and variability of {alpha} can be obtained. The resulting water retention curves from the proposed models are also compared with the fitted two-parameter van Genuchten model (van Genuchten, 1980).


   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 MODEL FORMULATION
 COMPARISON WITH SOIL DATA
 RESULTS AND DISCUSSION
CONCLUSIONS
 REFERENCES
 
Particle-Size Distribution
Figure 2a shows the pdf and cdf of the lognormal particle-size distribution given by Eq. [1] and G0(R) in Eq. [3], respectively, for different values of µy and with {sigma}y = 0.5. The corresponding cumulative mass fraction (cmf) of particle sizes is also plotted using Eq. [28]. From the figure, the cmf appears to have the exact shape of the cdf. In fact by substituting in Eq. [2] and [3], Eq. [28] can be written as

[29]



View larger version (19K):
[in this window]
[in a new window]
  		Fig. 2. Particle-size distributions: (a) with constant {sigma}y and varying µy (thicker lines for µy = 2, thinner line for µy = 4); (b) with constant µy and varying {sigma}y (thicker lines for {sigma}y = 0.5, thinner line for µy = 0.8). The cumulative mass fraction (cmf) is plotted as dash–dot line. cdf, cumulative distribution function; pdf, probability density function.

 
This shows that the cmf, M(R), also follows a lognormal cumulative distribution function where lnR would have a mean of µy + 3{sigma}2y and a standard deviation of {sigma}y. The cmf is therefore shown on the figure as only a shift of the cdf to the right by a distance of exp(3{sigma}2y). Increasing µy shifts the cdf and cmf to the right while lowering the peak of the pdf. On the other hand, increasing {sigma}y provides more spread for the distribution, as shown in Fig. 2b. Since µy is kept constant, an increase in {sigma}y does not affect the median of the cdf but does increase that of the cmf.

Figure 3 shows a fit of the lognormal distribution to the particle-size data of the Beerse podzol II sand (code no. 4061). The particle radius, R, is plotted on the x axis on a log scale while the cumulative mass fraction is plotted on the y axis on a normal probability scale. Therefore, data falling on a straight line imply a lognormal distribution. For most of the sandy soils used in this study, their particle-size distributions resemble the one shown in Fig. 3. The lognormal cdf tends to deviate from the measured data below 10 to 20% of the cumulative mass. This observation is obvious when the measured and estimated cumulative mass fractions for all soils are plotted in Fig. 4 . Data points fall below the diagonal solid line where the cumulative mass fraction is <0.2. This is also reflected by the negative intercept of the linear regression. The point of this deviation usually lies at the boundary between the silt and sand (R {approx} 25 µm) texture classes. This fraction of clay and silt often exists as clay–silt coatings around sand grains, which are very common in sandy soils (FitzPatrick, 1993), and likely has little influence on the pore-size distribution. Nonetheless, the inability of the lognormal distribution to describe the silt and clay fractions can affect the predictability of the proposed water retention models, especially in the dry region where the pore sizes seem to be largely governed by the smallest particles. However, this part of the curve is asymptotic to the residual water content, {theta}r, and is "nailed down" by this parameter. The issue of proper estimation of {theta}r is discussed below in the section Performance of Models with {alpha} = 4. A bimodal lognormal model suggested by Shiozawa and Campbell (1991) should provide a better fit to the sand and loamy sand soils in this study (Hwang et al., 2002). Further studies are required to investigate the use of other distributional forms. Despite the discrepancy that occurs at smaller-size fractions, the overall fit of the lognormal distribution is good. The linear regression performed between the fitted and measured cumulative mass fractions has a very high R2 of 0.99 and a small RMSE of 0.034. The RMSE is calculated as follows:

[30]
where N is the number of data points and the subscripts "est" and "meas" stand for "estimated" and "measured," respectively. The same formula applies to other RSME calculations where M can be substituted with any other variable of interest.



View larger version (23K):
[in this window]
[in a new window]
  		Fig. 3. Particle-size distribution of Beerse podzol II sand (code no. 4061) plotted on a normal probability scale (y axis). Data points falling on a straight line follow a lognormal distribution. The solid line represents the best-fit curve given by Eq. [29].

 


View larger version (35K):
[in this window]
[in a new window]
  		Fig. 4. Measured and estimated cumulative mass fractions for all 119 sandy soils. The dashed line is a linear regression of the data points. The R2 and RMSE imply the goodness of fit of the dashed line.

 
Water Retention Curve
Figure 5a shows the effective pore-size distributions of the polydisperse sphere models. For the same particle-size distribution (identical µy and {sigma}y) and porosity, the pore-size distribution of the FPS model displays a longer tail than that of the TIS model. Spheres are allowed to overlap in the FPS model; therefore, the pore sizes tend to be larger than those in the TIS model. What this translates into is a water retention curve that has a lower air-entry pressure and a general decrease in pressure across the whole range of the water saturation (see Fig. 5b). If {alpha} is set to be variable, increasing its value will rescale pore-size distribution to a larger mean and variance; as a result, the suction heads will decrease due to larger pore sizes. This effect is similar to that of increasing the µy and {sigma}y of the particle-size distribution. As {sigma}y increases, the retention curve shifts downward, as shown in Fig. 6 .



View larger version (14K):
[in this window]
[in a new window]
  		Fig. 5. Comparison of the proposed totally impenetrable spheres (TIS) and fully penetrable spheres (FPS) models' (a) effective pore-size distribution and (b) water retention curve.

 


View larger version (20K):
[in this window]
[in a new window]
  		Fig. 6. The water retention curves of the totally impenetrable sphere models with varying {sigma}y.

 
Scaling Factor {alpha} as a Fitting Parameter
To investigate the variability of the scaling factor {alpha} and the validity of the assumption of a linear relationship between {delta} and rp, {alpha} along with {theta}r in the proposed water retention models are fitted for each soil. The measured and model-estimated normalized water content data pairs are plotted in Fig. 7 for all 119 soils. Note that the normalized water content, {theta}/{theta}s, is used to allow for a more meaningful comparison. Since {alpha} is fitted for each soil and for each proposed model, the overall performances of the proposed models are quite good. All proposed models perform similarly because their particle-size distributions have very similar shapes, as shown in Fig. 5a. The TIS model provides slightly better fit than the FPS model, with R2 closer to one and with smaller RMSE. In general, the proposed models compare favorably with the van Genuchten model, as shown in Fig. 7.



View larger version (29K):
[in this window]
[in a new window]
  		Fig. 7. Estimated and measured normalized water contents for all 119 soils. The scaling factor {alpha} is treated as a fitting parameter in obtaining the totally impenetrable spheres (TIS) and fully penetrable spheres (FPS) model estimates. VG, van Genuchten model.

 
It is also observed from Fig. 7 that quite a few data points near full saturation lie above the diagonal line, indicating that the models are overestimating the measured water content. One example is given in Fig. 8 , with Fig. 8a showing the fit to the cmf for Wagram sand (code no. 1142). Figure 8b shows the corresponding water retention curve for the soil. The retention models fit the measured data very well, except for the few data points near saturation. This effect is, however, observed in many unconsolidated sands (Durner, 1994). Although boundary effects during the initial drainage can give rise to this observation, Durner (1994) attributed this phenomenon to a heterogeneous (secondary) pore system in cases where boundary effects can be ruled out. For this soil, the particle-size distribution deviates from a lognormal distribution in the silt and clay fractions like most of the soils used in this study. However, it is not immediately clear that a bimodal particle-size distribution will be able to account for the considerable change in water content near saturation. The secondary pore system suggested by Durner may not be reflected entirely in the particle-size distribution, but the arrangement of particles seems to be more critical. In this study, we assume random packing of the soil particles and therefore the proposed retention models fail to capture the effects of any heterogeneous pore system (of the kind proposed by Durner) exhibited by the actual soil.



View larger version (17K):
[in this window]
[in a new window]
  		Fig. 8. Wagram sand (code no. 1142): (a) measured particle-size distribution data and the lognormal distribution fit (solid line); (b) measured water retention curve and the proposed models with fitted {alpha}. Considerable change in the water content occurs in the high saturation range, perhaps indicating a secondary pore system. TIS, totally impenetrable spheres model; FPS, fully penetrable spheres model; and VG, van Genuchten model.

 
Caution has to be exercised when applying the proposed models based on random arrangement of particles, or any other unimodal retention model, to soils like the one shown in Fig. 8. It is often observed that a slight discrepancy in the retention curve can translate into a rather substantial difference in predicting the relative hydraulic conductivities (Chan and Govindaraju, 2003). This kind of sensitive dependence between the {theta}–K relationship also exists with empirical models (see van Genuchten, 1980). In this case, it might be justifiable to fit {theta}s and ignore the initial change occurring near full saturation. Otherwise, a different water retention model involving a bimodal pore distribution should be used instead.

Variability of {alpha}
The curve-fitting results show a moderate variability of {alpha} where the coefficients of variation (the standard deviation divided by the mean) are 0.37 and 0.39 for the TIS and FPS model, respectively (Table 1). The theoretical value of {alpha} (Eq. [4]) is close to the median value of fitted {alpha} (4.29) for TIS. The median for the FPS model, however, is considerably smaller than the theoretically derived value. As mentioned, the overlapping of spheres in FPS results in larger pore sizes. The scaling factor therefore accounts for not only the bias but also the correction for the overlapping particles.


View this table:
[in this window]
[in a new window]
  		Table 1. Statistics of the scaling factor {alpha} in the totally impenetrable spheres (TIS) and fully penetrable spheres (FPS) retention models.

 
Figure 9 reveals the possible correlation between {alpha} and the parameters of the particle-size distribution. As seen in the figure, there is a slight increase in {alpha} as µy decreases. When {alpha} is plotted against {sigma}y, a similar trend is observed, except that now an increasing {sigma}y corresponds to a slight increase in {alpha}. For µy < 1.5 and {sigma}y > 1, the estimated values of {alpha} are significantly higher than the theoretical values. A larger {alpha} suggests that the effective pore sizes are larger than those predicted by the theoretical value. These data points correspond to the loamy sand soils, which also have the larger values of {sigma}y. One possible explanation of this trend is that the arrangement of the soil particles can no longer be assumed random with increasing fractions of silt and clay. As discussed above, the silt and clay can form a coating around sand grains that sometimes results in bridges that connect individual sand particles together. Larger pores would be expected as a result. The variability of {alpha} also suggests that the assumed linear relationship between {delta} and rp might not be universal for all texture classes. However, for the TIS model, {alpha} is rather constant with exceptions only for a few loamy sand soils.



View larger version (40K):
[in this window]
[in a new window]
  		Fig. 9. Scaling factor {alpha} vs. the mean parameter µy and the standard deviation parameter {sigma}y of the lognormal distribution. The dashed line indicates the theoretical value of {alpha}. TIS, totally impenetrable spheres model; FPS, fully penetrable spheres model.

 
Performance of Models with {alpha} = 4
If the theoretical value of {alpha} is used for the proposed water retention models, the only remaining unknown parameter is the residual water content, {theta}r, which can be obtained in a number of ways, including using the data point at very high suction (>104 cm), as suggested by van Genuchten (1980). In this study, we use the {theta}r obtained by simultaneously fitting {alpha} and {theta}r of the models to the soil data. The values of {theta}r are therefore different for each soil and for each model. This procedure seems to provide good estimates of {theta}r, as seen from Fig. 7 where the scatter of data pairs at the lower end of the normalized water content is relatively small. Also using the same values of {theta}r enables a reasonable comparison between models with varying {alpha} and with theoretical {alpha}. The treatment of {theta}r deserves some attention here, as it is still a controversial topic (Nimmo, 1991; Luckner et al., 1991; Kosugi, 1994). We conjecture that the residual water is constituted of (i) water that is trapped within small pores or pores that are inaccessible to drainage and (ii) water that is being adsorbed on the soil particles that cannot be removed by hydromechanical forces (Luckner et al., 1989, 1991). In this study, {theta}r is treated as a fitting parameter; therefore, it imposes a cutoff that serves the first conjecture well but violates the second one because the adhesive forces are dependent on the varying water potential. A constitutive yet physically impossible cutoff imposed by {theta}r is therefore questionable. To solve this problem, the adsorption process has to be modeled discretely. A new approach proposed by Tuller et al. (1999) and Or and Tuller (1999) can be adopted to incorporate the adsorptive contribution through the use of the shifted Young–Laplace equation. More work is needed to incorporate this absorptive component into our models.

Using the theoretical {alpha} instead of treating {alpha} as a fitting parameter, one could expect the models not to perform as well. Figure 10 shows the predicted normalized water contents vs. the measured ones. More scatter of the data points is observed in Fig. 10 than in Fig. 7. Less scatter is observed near the dry range of water content because fitted values of {theta}r are used in the model predictions. The goodness of fit by the FPS model is affected the most. The TIS remains the better model, with higher R2 value of 0.91. For each soil, the RMSE of the model prediction was calculated, and the mean and standard deviation of the RMSEs are presented in Table 2. Statistics from Table 2 confirm the observation between Fig. 7 and 10. The increase in error of prediction due to fixing {alpha} is substantial for the FPS model. The TIS model shows smaller increase in the mean and the standard deviation of the RMSEs. The two-parameter van Genuchten model provides the best overall fit to the soils studied with smallest mean RMSE.



View larger version (54K):
[in this window]
[in a new window]
  		Fig. 10. Estimated and measured normalized water contents for all 119 soils. The scaling factor {alpha} is fixed at 4. A wider scatter of data points is expected as a result. TIS, totally impenetrable spheres model; FPS, fully penetrable spheres model.

 

View this table:
[in this window]
[in a new window]
  		Table 2. The mean and standard deviation of RMSEs of the totally impenetrable spheres (TIS), fully penetrable spheres (FPS), and van Genuchten (VG) model predictions.

 
In Fig. 11 , the RMSEs for the TIS model with theoretical {alpha} are plotted against the estimated µy of a soil's particle-size distribution. A group of data points with rather high RMSE is located at the smaller end of µy. These data points correspond to the same ones in Fig. 9 that have high values of {alpha}. This mismatch of {alpha} is demonstrated with an example in Fig. 12 where the soil water retention curve of Troup loamy sand (code no. 1012) is compared with the proposed models for both cases of fitted and fixed {alpha}. The Troup loamy sand has a relatively small µy of 1.843 and a large {sigma}y of 0.9112. When {alpha} is allowed to vary, all models fit very well to the measured data; however, the estimated values of {alpha} are much higher than the theoretical ones for the TIS models. When the model predictions with the theoretical {alpha} are plotted against the measured data, it is obvious that the TIS model overestimates the water contents. The retention curves estimated by these models are well above the measured data. The possible reason for this mismatch between the fitted and theoretical {alpha} was discussed above.



View larger version (24K):
[in this window]
[in a new window]
  		Fig. 11. The RMSE of the totally impenetrable spheres (TIS) model prediction on water retention data vs. the mean parameter µy of the lognormal particle-size distribution.

 


View larger version (23K):
[in this window]
[in a new window]
  		Fig. 12. The measured and estimated soil water retention curves of Troup loamy sand (code no. 1012) (a) with fitted {alpha} and (b) with theoretical {alpha}. The fitted {alpha} value is higher than the theoretical one. TIS, totally impenetrable spheres model; FPS, fully penetrable spheres model; and VG, van Genuchten model.

 
For a number of soils, the high RMSE is not due to the mismatch of {alpha}, but rather, a result of the mismatch in the shape of the model retention curve and the measured data. For the Kootwijk sand (code no. 4520), which has the particle-size characteristics of µy = 3.936 and {sigma}y = 0.5162, the measured water retention curve resembles a step function with a sharp drop-off in water content beyond the air-entry pressure (Fig. 13) . The proposed models cannot capture such step function–like behavior. The TIS model underestimates suction pressure in the wet range of the curve while overestimating a portion of the curve at the dry range. At very high pressures, the water content approaches zero very gradually, whereas the proposed models predicted complete drainage of water at a much lower pressure. This soil case demonstrates two limitations of the proposed theory. First, the step function–like retention curve of the measured data suggests that the soil has a fairly uniform pore system. However, the proposed theory assumes that the soil particles are randomly placed. Therefore, the predicted pore-sizes will not be uniform but follow the specified distribution. Second, the relationship assumed between {delta} and rp may not be a simple linear function. Further research is necessary to find a more appropriate relationship.



View larger version (24K):
[in this window]
[in a new window]
  		Fig. 13. The measured and estimated soil water retention curves of Kootwijk sand (code no. 4520) (a) with fitted {alpha} (b) with theoretical {alpha}. The proposed TIS model does not match the shape of the measured retention curve very well. TIS, totally impenetrable spheres model; FPS, fully penetrable spheres model; and VG, van Genuchten model.

 
The TIS model was found to provide a better estimation of the retention curve than the FPS model because its model assumptions are in better agreement with the sandy and loamy sand soils used in this study. Figure 14 shows the retention curve of Beerze podzol II sand (code no. 4061), and its particle-size distribution is given in Fig. 3. The TIS model is in close agreement with the van Genuchten model and the measured data.



View larger version (23K):
[in this window]
[in a new window]
  		Fig. 14. The measured and estimated soil water retention curves of Beerze podzol II sand (code no. 4061) (a) with fitted {alpha} (b) with theoretical {alpha}. The TIS model provides a very good fit to the data. TIS, totally impenetrable spheres model; FPS, fully penetrable spheres model; and VG, van Genuchten model.

 

   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 MODEL FORMULATION
 COMPARISON WITH SOIL DATA
 RESULTS AND DISCUSSION
 CONCLUSIONS
REFERENCES
 
Models of soil water retention were developed in this study using the theory of polydisperse sphere systems. We modeled the soil as a system of lognormally distributed spheres that are either fully penetrable or totally impenetrable. The void nearest-surface function was used in an attempt to describe the pore structure of these many-particle systems. A linear relationship was assumed to relate this nearest-surface function to the effective pore-size distribution, so the water retention curve could be obtained using the classical Young–Laplace equation. Models of water retention were derived for the totally impenetrable spheres and the fully penetrable spheres.

Since the proposed models are derived from statistical considerations (note, however, that the residual water content, {theta}r, still remains a fitting parameter), the performance is expected to be inferior to the empirically fitted models. In fact, the success of the proposed model is predicated on how closely the soil matches the assumptions of the model. In this study, we found that the TIS model works well in predicting the water retention curve for many sandy and loamy sand soils. However, for soils with µy < 1.5 and {sigma}y > 1, using the theoretical scaling factor ({alpha} = 4) can no longer provide reasonable estimates of the retention curve. We reason that the placement of the soil particles in soils with larger fractions of silt and clay may not be completely random as specified by the theory because the silt and clay can change the packing structure of the sand grains. A polydisperse particle system with "sticky" spheres (Torquato, 2002) may be able to model this phenomenon with better success. A "sticky factor" can be introduced to account for the degree of stickiness of the sphere particles, thus allowing aggregates to form. A sticky-sphere model might serve better in modeling a larger range of soils, as a recent study has found that there is a significant correlation between the median aggregate size and median pore size (Lebron et al., 2002). Another limitation of the proposed model is that the assumed linear relationship between {delta} and rp might not hold for soils other than sand and loamy sand. More research is required to investigate these issues.

The residual water content becomes a relevant issue if we would like the proposed models to be purely predictive without additional calibration or fitting. Further work is necessary to couple the residual water concept with a general physically based adsorption process dominating at low saturation. Finally, future research could be directed toward the prediction of relative hydraulic conductivity and other hydraulic properties of interest.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 MODEL FORMULATION
 COMPARISON WITH SOIL DATA
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES