Published online 12 October 2005
Published in Vadose Zone J 4:959-966 (2005)
DOI: 10.2136/vzj2005.0012
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
Modeling the Primary Drainage Curve of Prefractal Porous Media
E. Perfect*
Dep. of Earth and Planetary Sciences, Univ. of Tennessee, Knoxville, TN 37996-1410
* Corresponding author (eperfect{at}utk.edu)
Received 20 January 2005.
 |
ABSTRACT
|
|---|
Fractal models for the soil water retention curve have largely ignored hysteresis. Such models generally require that all pores be connected to the atmosphere via a network of similar-sized or larger pores. In random mass fractals and natural porous media, however, not all pores of a given size class empty during drainage due to incomplete pore connectivity and/or the presence of surrounding smaller pores. A modification of Rieu and Sposito's (1991a) prefractal water retention equation is proposed to accommodate this hysteresis during monotonic drying. The resulting expression is identical in form to the empirical Brooks and Corey (1964) model and contains three physically-based parameters: the proportion of nondraining pores (Pd), the scale factor (b), and the apparent fractal dimension of the primary drainage curve (Dd), which is related to the underlying mass fractal dimension (D) of the porous medium by Dd = D + [log(Pd)]/[log(b)]. Model testing consisted of fitting the modified equation to previously published water retention data for 30 randomized Sierpinski carpets and seven soils. Error sums of squares ranged from 0.001 to 0.093 for the carpets, and from 0.015 to 0.047 for the soils. In contrast, the unmodified Rieu and Sposito (1991a) equation performed very poorly. Best fit estimates of Dd ranged from 0.295... to 1.640... for the two-dimensional carpets, and from 2.467... to 2.902... for the three-dimensional soils. Prediction of D, based on estimates of Dd derived from the primary drainage curve, will require additional research on how to obtain Pd and b. Based on the analytical approach outlined here, it should also be possible to model the primary wetting curve, thus providing a more complete fractal description of soil water hysteresis.
Abbreviations: ESS, error sums of squares
 |
INTRODUCTION
|
|---|
THE RELATIONSHIP BETWEEN water content and suction differs depending on whether the soil is wetting or drying; this path dependence is known as hysteresis. Hysteresis occurs because of differences in pore shape, size, and interconnectivity (Hillel, 1998). With irregularly-shaped pores, drainage is determined by the smallest opening or neck, while water entry is controlled by the dimensions of the main body. Small pores fill up first and empty last, while large pores empty first and fill up last. However, large pores that are connected to the atmosphere via smaller ones will not empty until the smaller ones have drained. The presence of entrapped air, differences in contact angle during wetting vs. drying, and shrinkage and swelling phenomena also contribute to hysteresis (Hillel, 1998).
As a result of hysteresis, the water retention curve depends on the wetting and drying history of the soil. For any given suction, the soil water content during monotonic drying will be greater than or equal to that during monotonic wetting. Assuming continuous drainage occurs from complete saturation to oven dryness and vice versa, the water retention curve can be characterized by its main drying and wetting branches. Under natural conditions, however, the water content at any given suction may vary between these envelopes depending on the wetting and drying history of the soil. As a result, a family of scanning loops bounded by the main wetting and drying curves is necessary to fully describe the water retention curve (e.g., Poulovassilis, 1962; Poulovassilis and Childs, 1971).
For modeling purposes, experimental water content vs. suction data are frequently parameterized by fitting an appropriate mathematical relationship. A variety of empirical equations have been developed for this task (e.g., Brooks and Corey, 1964; van Genuchten, 1980). These models are generally fitted to the main drying branch of the water retention curve. Based on the resulting parameter estimates it is possible to construct the primary wetting curve and any intermediate scanning loops (e.g., Mualem, 1984a; Parker and Lenhard, 1987; Braddock et al., 2001). Jaynes (1984)(/85) and Viaene et al. (1994) have evaluated the performance of these different modeling approaches in terms of their ability to fit hysteretic water retention data. However, none of them can be used to back out information about the connectivity and size-distribution of pores.
Equations based on fractal geometry are increasingly being used to parameterize the soil water retention curve (e.g., Tyler and Wheatcraft, 1990; Toledo et al., 1990; Rieu and Sposito, 1991a; Perfect, 1999; Bird et al., 2000). Soil aggregates have been shown to exhibit mass fractal scaling over a finite range of lengths scales (Anderson and McBratney, 1995). The theoretical soil water retention function for this type of "prefractal" porous medium is given by (Rieu and Sposito, 1991a):
 | [1] |
where S
/
is the relative saturation,
is the volumetric water content,
is the porosity, h is the suction, hmin is the suction that drains the largest pores, E is the Euclidean dimension of the initiator (see below), and D is the mass fractal dimension. Equation [1] is valid for hmin
h
hmax, where hmax = hmin(1-
)1/(DE). Note that when h = hmin, S = 1, and when h = hmax, S = 0.
Equation [1], like other fractal models, is based on a number of underlying assumptions. Individual pores are assumed to be either empty or full, their status being determined solely by the Young-Laplace capillary relation (de Gennes et al., 2004). The possibility that two or more pores may coalesce, resulting in drainage at a lower suction, is neglected. Furthermore, all pores are assumed to be connected to the atmosphere via a network of progressively larger pores. This research is concerned with the later assumption.
In random fractals and natural porous media, a certain portion of the pore space does not drain during monotonic drying because pores are either completely disconnected or are connected to the atmosphere via smaller pores. Numerical capillary drainage simulations in random fractal structures have demonstrated how a lack of pore connectivity causes the observed water retention curve to deviate from that predicted by analytical models (Perrier et al., 1995; Bird and Dexter, 1997). The magnitude of this discrepancy increases as D approaches the critical fractal dimension, Dc, for pore percolation (Sukop et al., 2001). Pore networks in structures with mass fractal dimensions close to Dc are poorly connected and, as a result, air has only limited access to the interior. In contrast when D < < Dc, pores are sufficiently well connected to facilitate air invasion of the pore network, so that observed relative saturations are much closer to model predictions (Sukop et al., 2001).
Equation [1] has been fitted to the main drying branches of water retention curves determined experimentally on soils from a wide range of textural classes (Rieu and Sposito, 1991a, 1991b). The resulting best fit values for D are generally >2.90, and approach three with increasing clay content. Such estimates are assumed to represent the mass fractal dimension of the entire pore network. Because of hysteresis, however, it is to be expected that "apparent" D values, estimated from the main drying branch of the water retention curve, will deviate systematically from actual D values for fractal porous media.
This paper presents a scale-invariant conceptualization of incomplete pore drainage during monotonic drying of random fractal porous media. This conceptualization is incorporated into the analytical prefractal model used to derive Eq. [1]. The resulting theoretical expression is then evaluated using previously published numerical and experimental drainage curves for random fractals and soils.
 |
THEORY
|
|---|
Equation [1] assumes that soil can be modeled as a mass prefractal porous medium. Such fractals are constructed from a solid initiator by an iterative process of mass removal and re-scaling. As a concrete example, consider the well-known Menger sponge that scales by a constant ratio of b = 3 (Mandelbrot, 1982). The initiator is a solid cube (E = 3) of unit length. A generator is then defined by subdividing the initiator into bE = 27 smaller cubes of length
= 1/b = 1/3, and removing n = 7 of them. Construction continues by repeatedly applying the generator to the remaining solid cubes. Note that
depends on b as
= 1/bi, where i = 0,1,2,3... is the level of iteration of the fractal algorithm. The number of solid cubes of length
,
s(
), at the first iteration is
s(1/3) = 20. At the second iteration
s(1/9) = 400, and so on. In general, we have:
 | [2] |
where D is the mass fractal dimension defined by the ratio log(bEn)/log(b), with n being the number of cubes removed in the generator.
The number of pores of length
, Np(
), in a mass prefractal porous medium is given by:
 | [3] |
where
s(b
) is the number of solid cubes of length b
. At saturation, all of the pores are assumed to be completely filled with water. As drying occurs not all pores of a given size will empty at the appropriate suction because of incomplete pore connectivity and/or the presence of adjacent smaller pores. Let the actual number of pores draining be denoted by Nd(
). It then follows from Eq. [3] that,
 | [4] |
where Dd is the fractal dimension for the drained pore network. Depending on the geometrical arrangement (lacunarity) of the prefractal porous medium, those pores of length
that do not drain at the appropriate suction may remain full, or empty into pores of length
/b,
/b2,
/b3... or
/bi1 as h
. Similarly, nondraining pores of length
/b may never empty, or later drain into pores of length
/b2,
/b3,
/b4... or
/bi1. To model this complex drainage process, it is assumed that the proportion of empty pores at the ith level does not change over time. In other words, the number of pores filling due to water draining from larger pores is offset by an equal number of pores emptying due to water draining into smaller pores.
The probability of pores of length
emptying during drainage, Pd(
), is defined by,
 | [5] |
By assuming the same proportions of pores empty at each iteration level, that is, Pd(
) = Pd(b
) in Eq. [5], it is easy to show that
 | [6] |
where Pd is the scale-invariant drainage probability for the pore network.
Taking the logarithms of both sides of Eq. [6] and rearranging gives:
 | [7] |
If Pd and b are known, Eq. [7] can be used to estimate the mass fractal dimension of a porous medium from the fractal dimension for the drained pore network. Alternatively, if D and Dd are known, Eq. [7] can be used to infer values for Pd and b. For example, Crawford et al. (1995) reported Dd values (estimated from the water retention curve) ranging from 2.90 to 2.97, along with corresponding D values (obtained from thin section analysis) of between 2.94 and 2.98, for eight Japanese soils. Inserting these values into Eq. [7] suggests that 0.8 < Pd
1 for a wide range of possible b values between 2 and 100.
Using the concepts developed above it is possible to modify the Rieu and Sposito (1991a) analytical model to allow for scale-invariant incomplete drainage. The number of boxes of length
needed to fill the pore phase of a mass prefractal porous medium,
p(
), is given by,
 | [8] |
Based on Eq. [8], the cumulative volume of the pore phase, Vp(
), can be calculated from:
 | [9] |
The porosity of the prefractal porous medium is defined when i
j and
min in Eq. [9], where j is an integer <
and
min > 0 is the length of the smallest pores present, that is,
 | [10] |
The number of boxes of length
needed to fill only the drained portion of pore network,
d(
), is given by
 | [11] |
It follows from Eq. [11], that when i = 0,
d(1) = 0 and when i = 1,
d(1/b) = Pd(bE bD) = Pd
p(1/b). The cumulative volume of the drained pore space, Vd(
) can be calculated from Eq. [11] using:
 | [12] |
The volumetric water content of the incompletely drained prefractal porous medium is given by
 | [13] |
Expressed in terms of relative saturation, Eq. [13] becomes
 | [14] |
Invoking the Young-Laplace expression (de Gennes et al., 2004), 1/bi in Eq. [14] can be replaced with hmin/h giving
 | [15] |
Substituting Eq. [7] into Eq. [15] yields
 | [16] |
Equation [16] is valid for hmin
h
hmax. When h = hmin, S = 1, and when h = hmax, S = Sr, the residual relative saturation, which is defined by setting h = hmin(1
)1/(DE) in Eq. [16]. Note that when Pd = 1, Eq. [16] is identical to Eq. [1], and Sr = 0.
The above derivation is similar to that of most other fractal models in its neglect of pore coalescence and assumption of purely capillary behavior. It differs from most others, however, in that it explicitly considers partial pore connectivity.
For the purpose of fitting Eq. [16] to the main drying branch of experimentally determined soil water retention curves, it is convenient to simplify it by collecting the b, Pd, and D parameters into two compound terms, that is
 | [17] |
where
= Pd
and ß = Dd = D +
. Assuming E and
are known, Eq. [17] can be fitted to experimental drainage curves as a nonlinear model with three unknown parameters: hmin,
, and ß. The ß term provides an explicit linkage between the main drying branch of a soil water retention curve, and the mass fractal dimension of the underlying matrix. However, to estimate D from such data it is necessary to know, or be able to independently estimate b and Pd.
 |
MATERIALS AND METHODS
|
|---|
Numerical Drainage Curves
Equation [16] was fitted to numerically simulated monotonic drainage curves for the random two-dimensional prefractal porous media investigated by Sukop et al. (2001). These authors generated 10 realizations of E = 2, b = 3, and j = 5 randomized Sierpinski carpets for each of three different generators (n = 1, 2, and 3, that is, D = 1.892...,1.771..., and 1.630..., respectively) using the homogenous algorithm. The main drying branches of the water retention curves for these prefractal structures were simulated using a numerical model for capillary drainage developed by Bird and Dexter (1997). Each drainage curve consisted of six paired values of S and h/hmin. Basic information on these simulations is summarized in Table 1. The reader is referred to the papers by Sukop et al. (2001) and Bird and Dexter (1997) for further details on the programs employed.
View this table:
[in this window]
[in a new window]
|
Table 1. Physical attributes (mass fractal dimension, D, scale factor, b, maximum iteration level, j, porosity, , and scaled suction, h/hmin) of the randomized Sierpinski carpets and drainage simulations from Sukop et al. (2001).
|
|
Equation [1] was used for the forward prediction of S based on h/hmin and the known values of D, E, and
for each prefractal construction. Equation [16] was fitted to the numerically simulated drainage curves using the Marquardt nonlinear regression method in PROC NLIN of the SAS statistical software package (SAS Institute Inc., 1999). Pd in Eq. [16] was the only unknown parameter estimated. All of the fits converged according to the software default criterion. Goodness of fit was quantified in terms of the error sums of squares (ESS) for the individual fits, and by linear regression analysis of the pooled observed and predicted S values.
Experimental Drainage Curves
Since values of b and D are unknown a priori for natural porous media, Eq. [1] was applied inversely along with Eq. [17]. Both equations were fitted to the main drying branches of experimentally determined water retention curves for the six Washington State soils investigated by Campbell and Shiozawa (1992). Each curve consisted of between 31 and 39 paired measurements of S and h, with h ranging from 3.1 x 100 to 3.3 x 105 kPa. Details about the physical characterization of these soils are given in Table 2. The reader is referred to the original publication for complete information on the different methods used to construct the monotonic drainage curves.
A segmented nonlinear regression procedure was used for the fitting, with S = 1 when h
hmin and either Eq. [1] or Eq. [17] applicable when h > hmin. The Marquardt method in PROC NLIN of the SAS statistical software package (SAS Institute Inc., 1999) was employed, with E = 3,
set to the volumetric water content at the lowest suction, and hmin and D in the case of Eq. [1], or hmin,
, and ß in the case of Eq. [17], treated as unknown parameters. All of the fits converged according to the software default criterion. Goodness of fit was assessed using the ESS for the individual fits, and by linear regression analysis of the pooled observed and predicted S values.
 |
RESULTS
|
|---|
Numerical Drainage Curves
The total porosity of the randomized Sierpinski carpets decreased as the D-value was increased (Table 1). As discussed by Sukop et al. (2001) and illustrated in Fig. 1
, the numerically simulated drainage curves for these prefractal structures deviated from those predicted by the Rieu and Sposito (1991a) analytical model, Eq. [1]. The magnitudes of these deviations were greatest for the high D-value (low
) carpets, and least for the low D-value (high
) carpets. This trend can be explained by increases in the degree of interconnectivity of the pore network with decreasing D or increasing
.

View larger version (10K):
[in this window]
[in a new window]
|
Fig. 1. Example of a simulated drainage curve for a b = 3, j = 5, D = 1.771... randomized Sierpinski carpet (realization #1), where S is the relative saturation and logb(h/hmin) is log to base b of the scaled suction. Also shown are the predicted curve based on Eq. [1] and the best fit curve from Eq. [16] with Pd = 0.786 (error sums of squares, ESS = 0.005).
|
|
Figure 1 shows the numerical drainage curve for a specific realization of a b = 3, D = 1.771..., j = 5 randomized Sierpinski carpet. As noted above, h/hmin = bi in the Young-Laplace equation. Taking logarithms of both sides yields log(h/hmin) = ilog(b), or logb(h/hmin) = i, indicating the scaled log of suction corresponds directly to the iteration level. The critical fractal dimension for pore phase percolation in randomized b = 3, j = 5 Sierpinski carpets is Dc
1.716 (Sukop et al., 2001). When D
Dc (Fig. 1), the pore phase is poorly connected and air has only limited access to the interior of the structure. As a result, relative saturation never approaches zero as predicted by Eq. [1]. In contrast, when D << Dc (e.g., D = 1.630...), the pore phase is sufficiently well connected to permit significant air invasion, and relative saturations approach those predicted by Eq. [1].
When treated as a model with one unknown parameter (i.e., Pd), Eq. [16] fitted the numerical drainage data very well, as can be seen from the example in Fig. 1. Overall, the ESS for the nonlinear fits ranged from 0.001 to 0.093 (Table 3). Since there were no significant differences between the mean ESS values for the different carpets (Table 3), goodness of fit can be considered independent of the degree of pore phase connectivity.
Pooled observed vs. predicted S values for all of the realizations are graphed in Fig. 2
. Linear regression analysis indicated the relationship between the numerical drainage data and analytical model predictions was approximately one to one (Fig. 2). Equation [16] assumes that the number of empty pores of a given size which fill due to drainage from larger pores is offset by an equal number of filled pores that empty due to water draining into smaller pores. This assumption is reasonable for intermediate and high suctions, but breaks down at low suctions since only a finite number of large water-filled pores is available to replenish those pores that drain over time. As a result, the model might be expected to overestimate S at low suctions. While this trend was clearly evident in some simulations (e.g., Fig. 1), there was no general pattern of over prediction at high S values (Fig. 2). In fact, the outliers in Fig. 2 suggest a slight tendency for the best fit predictions to underestimate S as i
0 in the low D-value structures. The reasons for this phenomenon are unknown.

View larger version (14K):
[in this window]
[in a new window]
|
Fig. 2. Comparison of relative saturation values from all of the numerical simulations (Numerical S) with the corresponding values predicted from the best fits of Eq. [16] (Analytical S).
|
|
Best fit estimates of the Pd parameter in Eq. [16] ranged from 0.173 to 0.973 (Table 3). The mean Pd value increased progressively with decreasing D-value (Table 3). This statistically significant trend denotes an increase in the scale-invariant fraction of pores that drain at each iteration level due to increased interconnectivity as the mass fractal dimension decreases. It is also apparent from Table 3 that the spread of best fit Pd values decreased with decreasing D. Randomized Sierpinski carpets with D-values at or above the critical value for percolation exhibited more variability in their drainage characteristics from one realization to another. Based on Eq. [7], the Pd values in Table 3 indicate a range in Dd of between 0.295... and 1.640... for these two-dimensional prefractal structures. Interestingly, this range was produced entirely by the estimates of Pd for the D = 1.892... carpets.
Experimental Drainage Curves
An example of one of the experimental water retention curves from Campbell and Shiozawa (1992) is shown in Fig. 3
. As expected the relative saturation values decreased in a curvilinear fashion with increasing suction.

View larger version (12K):
[in this window]
[in a new window]
|
Fig. 3. Example of an experimentally determined drainage curve for soil (Palouse silt loam), where S is the relative saturation and h is the suction. Also shown are the best fit curves for Eq. [1] and [17]. Parameter estimates for these curves are given in Tables 4 and 5, respectively.
|
|
Also included in Fig. 3 are the predictions from the best fits of Eq. [1] and [17]. It can be seen that the analytical model incorporating hysteresis, Eq. [17], did a much better job of describing the data than the nonhysteretic model, Eq. [1]. Because Eq. [1] has only two fitting parameters it is much less flexible than Eq. [17], which has three unknown parameters. As a result, Eq. [1] tended to underpredict the experimental S data at low and high suctions, and overpredict them at intermediate suctions (Fig. 3). This trend occurred for all of the soils, as can be seen from the pooled predicted vs. observed relative saturations (Fig. 4)
. Linear regression analysis of these data showed that the predictions from Eq. [17] were much closer to a one-to-one relationship with the observed values than those from Eq. [1]. The overall better goodness-of-fit of Eq. [17] as compared to Eq. [1] can also be seen in the ESS values from the individual fits (Tables 4 and 5). The mean EES for Eq. [17] was approximately an order of magnitude lower than that for Eq. [1].

View larger version (20K):
[in this window]
[in a new window]
|
Fig. 4. Comparison of relative saturation values for all of the soils (Observed S) with the corresponding values predicted from the best fits of Eq. [1] and [17] (Predicted S).
|
|
Best estimates of hmin and D obtained by fitting Eq. [1] to the experimental drainage curves are given in Table 4. The suction draining the largest pores present ranged from 0.001 to 2.385 kPa. The lowest values of hmin occurred with the sand and sandy loam soils, while the silt loam and silty clay soils gave the highest values. The mass fractal dimension only ranged from 2.941 to 2.990, with the sand and sandy loam soils giving the highest values, and the silt loam and silty clay soils the lowest values. Rieu and Sposito (1991a)(1991b) and Perfect (1999) reported a much wider range for D, and an opposite trend with soil texture. The reasons for this discrepancy are unclear, but may be related to the inflexibility of Eq. [1] when fitted to experimental drainage data as a two-parameter model. Specifying
(as was done here), instead of treating it as a fitting parameter, appears to have induced a negative correlation between hmin and D.
Application of Eq. [17] to the same experimental data yielded a different suite of parameter estimates (Table 5) that must be interpreted differently from those associated with Eq. [1]. The only parameter common to both models was hmin. Estimates of hmin from the best fits of Eq. [17] were approximately an order of magnitude greater than those from Eq. [1]. In this case, the sandy loam and silt loam soils produced the highest hmin values, while the sand and silty clay soils gave the lowest values. There was no significant correlation between the two sets of hmin estimates. Estimates of the
parameter ranged from 0.170 to 0.729, and generally increased with increasing clay content. The ß parameter ranged from 2.467 to 2.902, with the sand and sandy loam soils giving the lowest values, and the silt loam and silty clay soils the highest values. This is a reversal of the trend observed for the D values from Eq. [1], which is consistent with the results of Rieu and Sposito (1991a)( 1991b) and Perfect (1999). However, keep in mind that ß is a compound parameter, and cannot be directly equated with D unless Pd = 1. To estimate D from ß it is necessary to know the values for b and Pd.
 |
DISCUSSION AND CONCLUSIONS
|
|---|
The Brooks and Corey (1964) empirical water retention equation is widely used for parameterizing main drying branch data. This model takes the form:
 | [18] |
where Se is the effective saturation,
r is the residual water content, and
is a fitting parameter known as the pore-size distribution index. Tyler and Wheatcraft (1990) showed that in the limit i
,
1 and
r
0; in this case Eq. [1] and Eq. [18] are identical with
= E-D. However, a porosity of one is physically unrealistic for natural porous media. As shown previously, such media are best modeled as randomized prefractals with either Eq. [16] or Eq. [17] describing the main drying branch of the water retention curve. Rearranging Eq. [16], [17], and [18] so that
is the sole dependent variable reveals the following equivalencies:
 | [19] |
 | [20] |
Note that hmin appears in all three expressions. Equations [19] and [20] appear to provide the first physically-based interpretation of the Brooks and Corey (1964) parameters that does not require the unrealistic condition
1.
The advantage of using Eq. [16] or [17] over an empirical model for the water retention curve is that the parameters b, Pd, and D have precise physical meanings, and thus can potentially be measured independently. How to obtain D for the underlying matrix from estimates of ß (i.e., Dd) remains a challenge. The forward use of Eq. [7] requires knowledge of the scale-invariant drainage probability, as well as the scale factor. It might be possible to estimate Pd using the renormalization group method (Turcotte, 1997), and b from the sample size divided by the size of the largest pores present, which can be computed directly from hmin (see Perfect, 1997). However, further research is needed to investigate these approaches. In the mean time, it is certainly possible to obtain these parameters inversely. Image analysis of soil thin sections (Crawford et al., 1995; Anderson et al., 1996), small-angle X-ray scattering (Martin and Hurd, 1987; Schmidt, 1991; Beaucage, 1996), and X-ray computed tomography (Peyton et al., 1994; Perret et al., 2003) can all be used to measure D directly. By comparing estimates of D obtained with these methods to Dd values derived from water retention curves measured on the same material would provide an opportunity to more fully test Eq. [7], and estimate Pd and b inversely.
The Rieu and Sposito (1991a) water retention model, Eq. [1], on which the current approach is based, assumes that all pores evacuate completely in accordance with the Young-Laplace equation. In reality, some water will be retained in the form of thin films adsorbed onto pore surfaces and pendular structures in pore corners (Toledo et al., 1990; Tuller et al., 1999). The possibility of such partial pore drainage should be considered as a refinement in future modeling efforts. It is generally accepted that water retention by capillarity gives way to an adsorption controlled regime at high suctions. However, Fig. 3 and Table 5 demonstrate that Eq. [17] can be fitted over an extremely wide range of suctions without any obvious break in scaling that suggests a crossover point. This observation lends support to the simplification of capillary equilibrium that underlies the hysteretic model.
The present work has focused exclusively on the main drying branch of the water retention curve because this is the most commonly measured soil hydraulic property. However, information about the main wetting branch is also needed to more thoroughly describe hysteresis. There is no reason why the analytical equations developed here cannot be applied to monotonic wetting up from residual saturation. A scale-invariant wetting probability, Pw, would need to be defined, and then relationships derived between this parameter and the mass fractal dimension. Assuming constant D and b values for a porous medium, it is to be expected that that Pd < Pw, and thus Dd < Dw, where Dw is the fractal dimension for the water-filled pore network. This is because drainage is controlled mainly by the size of pore necks, whereas rewetting occurs in response to the size of pore bodies. Further theoretical and experimental research is needed to test this hypothesis.
If the main drying and wetting curves can both be described in terms of fractal parameters representing the underlying pore space geometry, this might allow a more physically-based evaluation of competing approaches for predicting one envelope curve based on characterization of the other (e.g., Mualem, 1977, 1984b; Parlange, 1976, 1980; Hogarth et al., 1988). Furthermore, if both primary curves can be established then existing conceptual models, such those proposed by Mualem (1974)(1984a), could be invoked to predict all intermediate scanning loops of interest.
In conclusion, a modification of Rieu and Sposito's (1991a) prefractal water retention equation has been presented, that allows for incomplete drainage of pores dues to hysteresis. The proportion of nondraining pores, Pd, is assumed to be constant and independent of pore size and suction level. The modified equation is equivalent to the full Brooks and Corey (1964) expression. It was fitted to 30 numerically simulated drainage curves for randomized Sierpinski carpets, and seven experimentally determined water retention curves for soils. For the Sierpinski carpets, fractal dimensions estimated from the main drying branch data were lower than the generator mass fractal dimensions. This difference can be quantitatively related to Pd, and is indicative of hysteresis. It means that fractal dimensions estimated from water retention curves do not equate to the mass fractal dimension of the porous medium. Based on these results we recommend adoption of the terminology "apparent fractal dimension" in future studies focused on the fractal modeling of soil hydraulic properties. For the soils investigated, the apparent fractal dimension ranged from 2.467 for a sand to 2.902 for a silty clay. Further research is needed before these values can be used for the forward estimation of mass fractal dimensions. A logical follow up would be to extend the current model to monotonic wetting.
 |
ACKNOWLEDGMENTS
|
|---|
J.F. McCarthy and Jie (Joe) Zhuang provided valuable comments on an earlier version of this paper. The Campbell and Shiozawa (1992) data were kindly supplied in worksheet format by J.R. Nimmo.
 |
REFERENCES
|
|---|
- Anderson, A.N., and A.B. McBratney. 1995. Soil aggregates as mass fractals. Aust. J. Soil Sci. 33:757772.[CrossRef]
- Anderson, A.N., A.B. McBratney, and E.A. FitzPatrick. 1996. Soil, mass, surface, and spectral fractal dimensions estimated from thin section photographs. Soil Sci. Soc. Am. J. 60:962969.[Abstract/Free Full Text]
- Beaucage, G. 1996. Small-angle scattering from polymeric mass fractals of arbitrary mass-fractal dimension. J. Appl. Crystallogr. 29:134146.[CrossRef]
- Bird, N.R.A., and A.R. Dexter. 1997. Simulation of soil water retention using random fractal networks. Eur. J. Soil Sci. 48:633641.[CrossRef]
- Bird, N.R.A., E. Perrier, and M. Rieu. 2000. The water retention function for a model of soil structure with pore and solid fractal distributions. Eur. J. Soil Sci. 51:5563.[CrossRef]
- Braddock, R.D., J.Y. Parlange, and H. Lee. 2001. Application of a soil water hysteresis model to simple water retention curves. Transp. Porous Media 44:407420.[CrossRef]
- Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous media. p. 2227. Hydrology Paper 3. Colorado State Univ., Fort Collins.
- Campbell, G.S., and S. Shiozawa. 1992. Prediction of hydraulic properties of soils using particle size distribution and bulk density data. p. 317328. In Int. Workshop on Indirect Methods for Estimating the Hydraulic Properties of Unsaturated Soils. Univ. of California Press, Berkeley.
- Crawford, J.W., N. Matsui, and I.M. Young. 1995. The relation between the moisture-release curve and the structure of soil. Eur. J. Soil Sci. 46:369375.[CrossRef]
- de Gennes, P.-G., F. Brochard-Wyart, and D. Quéré. 2004. Capillarity and wetting phenomena: Drops, bubbles, pearls, and waves. Springer-Verlag, New York.
- Hillel, D. 1998. Environmental soil physics. Academic Press, San Diego, CA.
- Hogarth, W.L., J. Hopmans, J.-Y. Parlange, and R. Haverkamp. 1988. Application of a simple soil-water hysteresis model. J. Hydrol. (Amsterdam) 98:2129.
- Jaynes, D.B. 1984/85. Comparison of soil-water hysteresis models. J. Hydrol. (Amsterdam) 75:287299.
- Mandelbrot, B.B. 1982. The fractal geometry of nature. W.H. Freeman and Co., New York.
- Martin, J.E., and A.J. Hurd. 1987. Scattering from fractals. J. Appl. Crystallogr. 20:6178.[CrossRef]
- Mualem, Y. 1974. A conceptual model of hysteresis. Water Resour. Res. 10:514520.
- Mualem, Y. 1977. Extension of the similarity hypothesis used for modeling the soil water characteristics. Water Resour. Res. 13:773780.
- Mualem, Y. 1984a. A modified dependent-domain theory of hysteresis. Soil Sci. 137:283291.
- Mualem, Y. 1984b. Prediction of the soil boundary wetting curve. Soil Sci. 137:379390.
- Parker, J.C., and R.J. Lenhard. 1987. A model for hysteretic constitutive relations governing multiphase flow: 1. Saturation-pressure relations. Water Resour. Res. 23:21872196.
- Parlange, J.-Y. 1976. Capillary hysteresis and the relationship between drying and wetting curves. Water Resour. Res. 12:224228.
- Parlange, J.-Y. 1980. Water transport in soils. Ann Rev. Fluid Mech. 12:77102.[CrossRef][ISI]
- Perret, J.S., S.O. Pasher, and A.R. Kacimov. 2003. Mass fractal dimension of soil macropores using computed tomography: From the box counting to the cube-counting algorithm. Eur. J. Soil Sci. 54:569579.[CrossRef]
- Perfect, E. 1997. Fractal models for the fragmentation of rocks and soils: A review. Eng. Geol. (Amsterdam) 48:185198.
- Perfect, E. 1999. Estimating soil mass fractal dimensions from water retention curves. Geoderma 88:221231.
- Perrier, E., C. Mullon, M. Rieu, and G. d Marsily. 1995. Computer construction of fractal soil structures: Simulation of their hydraulic and shrinkage properties. Water Resour. Res. 31:29272943.
- Peyton, R.L., C.J. Gantzer, S.H. Anderson, B.A. Haeffner, and P. Pfeifer. 1994. Fractal dimension to describe soil macropore structure using X-ray computed tomography. Water Resour. Res. 30:691700.[CrossRef]
- Poulovassilis, A. 1962. Hysteresis of pore water, an application of the concept of independent domains. Soil Sci. 93:405412.
- Poulovassilis, A., and E.E. Childs. 1971. The hysteresis of pore water: The non independence of domains. Soil Sci. 112:301312.
- Rieu, M., and G. Sposito. 1991a. Fractal fragmentation, soil porosity, and soil water properties: 1. Theory. Soil Sci. Soc. Am. J. 55:12311238.[Abstract/Free Full Text]
- Rieu, M., and G. Sposito. 1991b. Relation pression capillaire-teneur en eau dans les milieux poreux fragmentés et identification du caractère fractal de la structure des sols. C.R. Acad. Sci. Paris 312(II):14831489.
- SAS Institute Inc. 1999. SAS/STAT user's guide. Version 8. SAS Inst., Cary, NC.
- Schmidt, P.W. 1991. Small-angle scattering studies of disordered porous and fractal systems. J. Appl. Crystallogr. 24:414435.[CrossRef]
- Sukop, M.C., E. Perfect, and N.R.A. Bird. 2001. Water retention of prefractal porous media generated with the homogeneous and heterogeneous algorithms. Water Resour. Res. 37:26312636.[CrossRef]
- Toledo, P.G., R.A. Novy, H.T. Davis, and L.E. Scriven. 1990. Hydraulic conductivity of porous media at low water content. Soil Sci. Soc. Am. J. 54:673679.[Abstract/Free Full Text]
- Tuller, M., D. Or, and L.M. Dudley. 1999. Adsorption and capillary condensation in porous media: Liquid retention and interfacial configurations in angular pores. Water Resour. Res. 35:19491964.[CrossRef]
- Turcotte, D.L. 1997. Fractals and chaos in geology and geophysics. 2nd ed. Cambridge Univ. Press, New York.
- Tyler, S., and S.W. Wheatcraft. 1990. Fractal processes in soil water retention. Water Resour. Res. 26:10471054.[CrossRef]
- van Genuchten, M.Th. 1980. A closed form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892898.[Abstract/Free Full Text]
- Viaene, P., H. Vereecken, J. Diels, and J. Feyen. 1994. A statistical analysis of six hysteresis models for the moisture retention characteristic. Soil Sci. 157:345355.
This article has been cited by other articles:

|
 |

|
 |
 
A. Cihan, E. Perfect, and J. S. Tyner
Water Retention Models for Scale-Variant and Scale-Invariant Drainage of Mass Prefractal Porous Media
Vadose Zone J.,
October 8, 2007;
6(4):
786 - 792.
[Abstract]
[Full Text]
[PDF]
|
 |
|