Soil Science Society of America Journal 66:1104-1114 (2002)
© 2002 Soil Science Society of America
DIVISION S-1—SOIL PHYSICS
Analytical Models for Soil Pore-Size Distribution After Tillage
Feike J. Leij*,a,
Teamrat A. Ghezzeheib and
Dani Orb
a USDA-ARS, George E. Brown Jr. Salinity Lab. and Dep. of Environmental Sciences, Univ. of California, Riverside, CA 92507-4851
b Dep. of Plants, Soils, and Biometeorology, Utah State Univ., Logan, UT 84322
* Corresponding author (fleij{at}ussl.ars.usda.gov)
 |
ABSTRACT
|
|---|
Tillage causes soil fragmentation thereby increasing the proportion of interaggregate (structural) pore space. The resulting tilled layer tends to be structurally unstable as manifested by a gradual decrease in interaggregate porosity until a new equilibrium has been reached between external loads and internal capillary forces at a rate governed by the soil rheological properties. The soil pore-size distribution (PSD) will change accordingly with time. We have previously applied the Fokker-Planck equation (FPE) to describe the evolution of the PSD as the result of drift, dispersion, and degradation processes that affect the pore space in unstable soils. In this study, we provide closed-form solutions for PSD evolution, which can be used to predict temporal behavior of unsaturated soil hydraulic properties. Solutions and moments of the PSD were obtained in case: (i) drift and degradation coefficients depend on time and the dispersivity is constant and (ii) drift and dispersivity are also linearly related to pore size. Both solutions can model the reduction in pore size during the growing season while the second solution can account for a reduction in the dispersion of the PSD. The solutions for PSD were plotted for a mathematically convenient expression for the drift and degradation coefficients and for an expression derived from a model for soil aggregate coalescence. Experimental data on the settlement of a Millville (coarse-silty, carbonatic, mesic Typic Haploxeroll) silt loam during wetting and drying cycles were used to determine time-dependent drift and degradation coefficients according to this coalescence model. The solution for the PSD was used to independently predict the water retention curve, which exhibited a satisfactory agreement with experimental retention data at the end of two drying cycles.
Abbreviations: ADE, Advection-Dispersion equation • FPE, Fokker-Planck equation • PSD, pore-size distribution • SWC, soil water characteristic
|
INTRODUCTION
|
|---|
THE SOIL PSD greatly affects the movement of fluids and dissolved substances and impacts the thermal and mechanical properties of soils and other porous media. Soil tillage is intended to create conditions for gas, water, chemical, and heat movement in agricultural soils that provide an optimal habitat for crop growth. The upper part of the soil will generally be unstable after tillage. The porosity and median pore size will decrease over time because of soil settlement and filling of pore space instigated by mechanical compaction, wetting and drying, and biological activity. Consequently, associated soil hydraulic and transport properties will also vary over time (e.g., Mapa et al., 1986). Ahuja et al. (1998) reviewed the literature and reported that tillage temporarily increases water retention in the wet range. These authors proposed two approaches to predict the change in the retention curve according to the Brooks and Corey (1964) equation. Van Es et al. (1999) reported on the sources of variability of infiltration and retention parameters determined in 3-yr field studies. Temporal variability was more significant than spatial variability for all properties whereas tillage practices contributed the most to the variation of infiltrability.
Quantifying the temporal dynamics of the PSD may have important applications, but to do so in a deterministic manner would be a difficult undertaking because of the complexity of the pore geometry and our incomplete understanding of relevant processes. Or et al. (2000) therefore proposed to describe the evolution of the PSD with the Fokker-Planck Equation (FPE) (Risken, 1989). The FPE is frequently used to quantify natural processes where the underlying processes or conditions cannot be precisely captivated. The coefficients of the FPE encompass our understanding of the mathematical behavior of the PSD in response to physical processes. A somewhat similar approach was followed by Ozkan and Ortoleva (2000), who used a Markov model for the evolution of the particle-size distribution because of breaking and shearing in fault zones of porous media.
Or et al. (2000) developed a conceptual framework for the evolution of the PSD after soil tillage by postulating that the PSD can be partitioned into a time-invariant textural component and a time-dependent structural component (Nimmo, 1997). Such partitioning is supported by experimental evidence that the tillage effect manifests itself only in the wet range for suctions <300 hPa (Ahuja et al., 1998). Following the work by Hara (1984), the drift coefficient in the FPE was defined by Or et al. (2000) from the displacement of the mean pore size and the dispersion coefficient as the product of drift coefficient and dispersivity. A pore-loss term accounts for the collapse (disappearance) of pores. The evolution of the PSD was assessed by numerically solving the FPE for a Millville silt loam (Or et al., 2000). The water retention curve or soil water characteristic (SWC) was determined by calculating the soil water pressure head from the pore size with the Laplace-Young equation while the saturated and unsaturated conductivities were estimated based on the Kozeny-Carman equation (Hillel, 1980) and the approach by Mualem (1976), respectively.
The current work is intended to augment the conceptual approach by Or et al. (2000) in two ways. First, the numerical solution procedure is quite adequate to predict the PSD but analytical solutions are more suitable to elucidate the effect of different coefficients on the evolution of the PSD and they can be used for moment analysis to determine transport parameters and predict trends. For the solution, we can utilize results obtained for the advection-dispersion equation (ADE), which is closely related to the FPE. For example, Su (1995) obtained explicit solutions for transport in heterogeneous media by making simplifying assumptions regarding the FPE. Second, we want to apply our solutions for the PSD to cases where the coefficients can be predicted according to a mechanistic model and compare our predictions with independent measurements. Ghezzehei and Or (2000) developed a model for coalescence of soil aggregates that we will use to predict the reduction in pore size and, hence, to estimate the coefficients in the FPE. Furthermore, we will use retention data obtained during the settling of a Millville silt loam to compare with the retention curve predicted from the solution for the PSD.
In view of the above, the objectives of the current work are: (i) derive closed-form solutions for the PSD and its first three moments with coefficients that depend on time and pore size, (ii) illustrate the behavior of these solutions for both a mathematically convenient and a physically based set of coefficients, and (iii) apply the solution for the PSD with time-dependent coefficients according to a coalescence model to describe experimental settlement and retention data for a Millville silt loam.
|
MATERIALS AND METHODS
|
|---|
Mathematical Models
General
The ADE will be used to describe the relative number of pores for each size fraction (i.e., the PSD) as a function of pore size and time. The ADE is a special form of the Fokker-Planck or second Kolmogorov equation (Bear, 1972). This equation has traditionally been used to model the probability density function associated with stochastic processes. A time-dependent normalization factor, not applied here, would be needed to convert the PSD into a true probability density function. The general mathematical problem is given by:
 | [1] |
where f is the PSD (L-1), t is time (T), r is the pore size or radius (L), V is the drift coefficient (L T-1), D is a dispersion coefficient (L T-2), and M is a decay coefficient (T-1). A source-sink term can readily be added to captivate changes in the soil PSD by external factors such as tillage, biological activity, and solution composition. The mathematical conditions are:
 | [2] |
 | [3] |
 | [4] |
Because the total number of pores can only be reduced because of degradation, the boundary conditions require a zero probability flux. For the initial distribution we select the following lognormal distribution of pore radii, which is already in use to describe soil water retention (compare Kosugi, 1994):
 | [5] |
where r0 is the median pore size, 
is the standard deviation of the log-transformed pore size, while 
s and r are the saturated and residual water contents.
If the transport coefficients do not depend on r and t, a solution to the problem may be readily written down from results for the solute transport equation (e.g., van Genuchten and Alves, 1982; Jury and Roth, 1990). It is reasonable to assume that the magnitude of the transport coefficients depends on time because changes in PSD will likely be most pronounced immediately after tillage and then taper off. Furthermore, pore-size evolution may affect larger pores more than smaller pores. As a first step, we are interested in a scenario where V, D, and M depend only on time. Later on, we consider a scenario where there is also a simple dependency on pore size. For both scenarios we assume that the dispersion coefficient is linearly related to the absolute value of the drift coefficient, i.e.:
 | [6] |
where 
is the dispersivity (L). The drift coefficient will be negative since the size of pores typically decreases over time after soil tillage or with compaction.
Drift Varies with Time
It is convenient to recast the time variable of Eq. [1] into a new independent variable for time (compare Barry and Sposito, 1989; Huang and van Genuchten, 1995):
 | [7] |
To further simplify the analysis, the first-order term, M, can be dropped from the governing equation through the following change of the dependent variable:
 | [8] |
If a constant dispersivity is used for the dispersion coefficient defined by Eq. [6], the problem given by Eq. [1] can be rewritten in terms of the new variables as:
 | [9] |
The mathematical conditions are
 | [10] |
 | [11] |
 | [12] |
The solution of the resulting problem follows from setting D = , v = -1, and t = T in the solution to the similar ADE as obtained with Green's functions for an arbitrary initial distribution (cf. Leij et al., 2000). The solution for the PSD in terms of the original variable, f, is:
The second integral in the solution was evaluated with Gauss-Chebyshev quadrature. The first three moments corresponding to this solution are given in the Appendix.
The time-dependent coefficients V and M still need to be specified. At this stage we will ignore degradation and use the following mathematically convenient and popular expression for V (compare Thornley, 1990):
 | [14] |
where <r> is the mean pore size with initial value <r0> (L) while a (T-1) and b (L) are empirical coefficients that characterize the temporal and the absolute value of the drift term. Note that T = <r> - <r0>. Figure 1
shows the PSD after 30 and 120 d as well as the initial distribution. For the solution, we assumed that there was no decay, a constant dispersivity = 1 µm (Fig. 1a) or = 0.1 µm (Fig. 1b) and a time-dependent drift term V(t). The behavior of V(t) is shown in the insert of Fig. 1a. Values for moments, which are further discussed in the Appendix, are included to characterize mass balance (m0), drift (M1), and variance (µ2). The figures illustrate that the drift of probability from larger to smaller pores gradually diminishes with time as demonstrated by the values for M1 at the three times whereas the variance increases. The distribution moves toward smaller pores. In Fig. 1a, dispersion causes the maximum probability to decrease with time and to occur at a lower pore size, i.e., f(r = 6.8 µm, t = 0 d) = 0.062 µm-1, f(6.3, 30) = 0.051 µm-1, and f(5.4, 120) = 0.042 µm-1. Because a zero probability flux is stipulated at r = 0, probability is reflected for larger times (t = 120 d). The value for m0 suggests that probability was preserved. Figure 1b demonstrates that a smaller dispersivity results in considerably less change of the PSD with maxima f(6.1, 30) = 0.061 µm-1, and f(5.1, 120) = 0.059 µm-1.
Drift Varies with Time and Pore Size
Few closed-form solutions can be found for the FPE (or ADE) where the coefficients depend on both time and pore size (position). In the following, we will consider a special case where both dispersion and drift coefficients depend on pore size and time. The drift term is written as the product of time- and pore-dependent parts, i.e.,
 | [15] |
The components u(r) (L) and w(t) (T-1) are taken positive and negative, respectively. The dispersion coefficient is defined as D(r,t)= -(r) u(r)w(t). For the behavior of the drift and dispersivity terms with pore size, we followed the simple relationship used by Zoppou and Knight (1997) to solve the ADE for a first-type boundary value problem by assuming:
 | [16] |
where 0 and u0 are dimensionless constants. The following governing equation can now be formulated from Eq. [1]:
 | [17] |
With the help of the independent and dependent variables
 | [18] |
 | [19] |
the governing equation can be written as:
 | [20] |
subject to:
 | [21] |
 | [22] |
|
| [23] |
After Laplace transformation with T, an Euler-type ordinary differential equation is obtained that can be readily converted into an equation with constant coefficients by using the transformation z = lnr (this transformation prompted the use of r = 1 rather than r = 0 as lower boundary). This equation is solved with a Laplace transform with respect to z (compare Leij et al., 1991). The solution for the PSD in terms of pore size, r, and cumulative time, T, is:
 | [24] |
Expressions for the associated zero-, first-, and second-order moments are also given in the Appendix.
Figure 2
shows the evolution of the same initial pore-size distribution as sketched in Fig. 1. Again, there is no decay and the time-dependent part of the drift term, w(t), is given by Eq. [14] and the previously defined parameters while u(r) and (r) are given by Eq. [16] with u0 = 0.5 and 0 = 0.25 (Fig. 2a) or 0.05 (Fig. 2b). The inserts in Fig. 2a illustrate the behavior of the functions w(t) and u(r). Considerable drift and dispersion of probability occurs during the initial stage when pore sizes are relatively large. The maximum probabilities in Fig. 2a are f(6.8, 0) = 0.062 µm-1, f(4.2, 30) = 0.046 µm-1, and f(2.3, 120) = 0.053 µm-1. There are slight errors in mass balance for t = 30 and 120 d of 0.1 and -0.6%, respectively. At greater times, probability is reflected at r = 1 µm because of the no-flux condition. The analytical solution can be computed more accurately for smaller 0 (compare Fig. 2b).
View larger version (21K):
[in this window]
[in a new window]
|
Fig. 2. Evolution of pore-size distribution over time according to Eq. [24] with no decay, a drift term v(r,t) = u(r)w(t) where w(t) is given by Eq. [14], with a = 0.01 d-1, b = 5 µm, <r0> = 7.5 µm (upper insert), and u(r) is given by Eq. [16] with uo = 0.5 (lower insert), and a dispersion term D(r,t) = -(r)v(r,t) where (r) is given by Eq. [16] with: (a) o = 0.25 and (b) o = 0.05. The initial distribution is again given by Eq. [5] with s = 0.469, r = 0.191, = 0.253, and r0 = 7.3 µm.
|
|
The values for the moments suggest that probability is preserved while the mean pore size decreases with time. For a higher dispersivity (Fig. 2a) the variance increases, but µ2 becomes smaller with time for the lower 0 (compare Fig. 2b). The earlier solution for the PSD with time-dependent coefficients cannot be used to describe such a reduction in spreading with time. However, experimental PSDs suggest that both mean and variance decrease over the growing season. Laliberte and Brooks (1967) determined PSDs for different porosities of a Columbia sandy loam soil; the mean and variance decrease with porosity (see also Fig. 3a of Or et al., 2000).
Applications
Coalescence Model
Ghezzehei and Or (2000) presented a model for the coalescence of soil aggregates based on an energy balance for capillary and rheological processes. The energy released by the change in the air-water interface is used for the viscous deformation of soil aggregates; aggregates coalesce and pore space is reduced. The soil aggregates are represented as spheres as shown in Fig. 3. Soil material at the contact points of aggregates is displaced to the periphery of the solid neck by viscous flow. The interaggregate strain 
(t) is defined as the ratio h/a, where h is the reduction in distance between aggregate center and plane of contact between two aggregates as a result of coalescence. An analytical solution was derived for the axial strain in case of a constant soil matric potential, 
(compare section 3.2 of Ghezzehei and Or, 2000). The interaggregate strain is given by the following expression:
 | [25] |
where a is the radius of a spherical aggregate in a cubic packing and h as the reduction in radius [L] because of the displacement of soil. The initial strain is given by o = (0) while the constants 
and ß are defined as:
 | [26] |
 | [27] |
Ghezzehei and Or (2000) determined by optimization that A = 10.873, B = 0.3887, and C = -4.891. Furthermore, is the surface tension for the air-water interface (e.g., 0.0729 Pa m for water at 20°C), 
is the yield stress (Pa) and 
is the coefficient of plastic viscosity of the aggregates (Pa s). The latter two were measured with a rotational rheometer using a parallel plate sensor system (Ghezzehei and Or, 2001), and were parameterized with the equations aexp(b) and aexp(b), respectively. For a Millville silt loam, Ghezzehei and Or (2000) reported that a = 35 Pa, b = 3 x 10-4 Pa-1, a = 5.6 x 104 Pa s, and b = 3 x 10-4 Pa-1.
Or et al. (2000) derived approximations to the FPE coefficients that depend on the above aggregate coalescence model. For mathematical simplicity, Or et al. (2000) considered a unit-cell model with cubic packing of spherical aggregates (Fig. 3a–b) to represent the entire structural (interaggregate) pore space. The unit-cell pore radius is equated to the geometric-mean radius (r0) of the initial PSD. The complex void of the unit cell is represented by a volume-equivalent spherical pore with radius, r. The ratio of the pore volume to the total volume of the unit-cell defines the porosity of the unit-cell, 
. Pore radius and porosity are related to the aggregate radius and the time-dependent strain by
 | [28] |
 | [29] |
Based upon these findings, we defined the drift and degradation coefficients as follows:
 | [30] |
 | [31] |
where the proportionality constant 
denotes the ratio of the porosity loss because of complete pore closure to reduction in pore size because of viscous deformation as given by Eq. [29]. The transformation to account for first-order degradation becomes:
 | [32] |
The strain rate should always be nonnegative as stipulated by Eq. [16] of Ghezzehei and Or (2000), who provide further constraints on the maximum and minimum matric potential. The matric potential has considerable bearing on the deformation process. The coalescence model is obviously inapplicable when the pore space becomes completely filled with either water or air. Deformation will often occur over a smaller range of potentials. During dry conditions, the high yield stress of the soil precludes deformation. If the soil is wet, the yield stress will be reduced considerably but may pose a threshold for coalescence because of the lack of an appreciable capillary stress.
We will illustrate the model for a Millville silt loam assuming an aggregate size a = 50 µm and a soil matric potential = -10900 Pa and an initial contact strain o = 0. These parameters are somewhat arbitrary. We assume that there is degradation. The time-dependent drift and degradation terms are given by Eq. [30] and [31], respectively, with = 0.1. The initial pore-size distribution of the interaggregate pore space is given by Eq. [5] with r0 = 48.45 µm, an interaggregate porosity (s - r) = 0.2 as well as = 0.2.
Drift Varies with Time
The evolution of the interaggregate pore-size distribution of a Millville silt loam is calculated according to Eq. [13]. The dispersion coefficient was calculated from the drift term assuming = 0.5 µm. Figure 4
shows the distribution at t = 0, 30, and 90 min as well as the corresponding values for the moments m0, M1, and µ2. The curves and the respective moments suggest that over time: (i) the total number of pores decreases because of degradation, (ii) the mean pore size decreases from the initial 49.4 to 43.1 µm after 90 min, and (iii) pore size extends over a wider range of values. The temporal behavior of several parameters is shown in Fig. 5
. Figure 5a sketches the strain (left axis) and potential (right axis) as a function of time. The strain, as calculated with Eq. [25], increases rapidly to 5.77% at 90 min. After 104 min, the minimum potential for coalescence has reached the prevailing soil water potential of -10900 and deformation will cease. At that time, the upper limit on the potential is -2090 Pa. In view of the lower limit, pore-size evolution was only simulated during the first 90 min. Figure 5b shows the drift term (left axis) and degradation factor (right axis). Both show a gradual decrease as the strain increases from V= -0.105 to -0.044 µm min-1 and from M = 0.162 to 0.073 x 103 min-1 during the 90-min period. Moments of the predicted PSD are given as a function with time in Fig. 5c. They accentuate the trend shown in Fig. 4 that the mean pore radius decreases and that the PSD becomes more dispersed.
Drift Varies with Time and Pore Size
The distribution of the pore sizes is now calculated with drift and dispersion coefficients that also depend on pore size using the solution given by Eq. [24]. Many of the model parameters are identical to those used for the previous example (Fig. 3 and 4). However, the drift coefficient is defined by Eq. [15] with the time-dependent part of the drift term, w(t), still given by Eq. [30] while uo = 0.1 for the pore-size dependent u(r). For the dispersion coefficient we used o = 0.05. Figure 6 shows the PSD at t = 0, 30, and 90 min and values for m0, M1 and µ2. The results suggest that over the total number of pores decreases, the mean pore size decreases from the initial 49.4 to 28.0 µm after 90 min, and that the maximum for the distribution remains fairly high and increases slightly (f[24, 90] = 0.0087 µm-1 and f[46, 0] = 0.0084 µm-1). A comparison of the moments listed in Fig. 4 and 6 reveals that m0 is identical because the degradation factor is independent of pore size, but that the linear dependency of drift on pore size resulted in an accelerated decrease of the mean pore-size. The dispersion coefficient has a quadratic relationship with pore size and initially the variance increases more rapidly (after 30 min. µ2 = 104 µm2 in Fig. 7
versus 102 µm2 in Fig. 5). As probability moves toward the lower boundary, the distribution starts to become compressed and the reduction in pore size implies a rapid decrease in the dispersion coefficient. The final variance of 85.1 µm2 is lower than the initial value for µ2 in Fig. 6. The following values of the skew were calculated from the predicted f(r): 0.595 (t = 0), 0.794 (t = 30 min), and 1.02 (t = 90 min). These values demonstrate that the PSD becomes more asymmetric as was already suggested by Fig. 6; the initial rapid drift at high values for r diminishes over time.
View larger version (15K):
[in this window]
[in a new window]
|
Fig. 7. Temporal behavior of the first-order normalized (M1) and second-order central (µ2) moments for the interaggregate pore-size distribution of a Millville silt loam shown in Fig. 5.
|
|
The temporal behavior of strain, matric potential range, and degradation terms are shown in Fig. 5. The behavior of drift and dispersion coefficients is different because of the added dependency on r. The behavior of the mean and variance, as inferred from pore moments, is given in Fig. 7. The mean pore size decreases rapidly and continuously during the 90-min period. At the start, the variance increases to a maximum of 105 µm2 at 18 min, after which the variance decreases because the reduction in size—and the accompanying decrease in V and D—of prevalent pores.
Wetting and Drying
The dynamic behavior of the PSD may be caused by many other factors than tillage. As an example, we used data by Silva (1995) on the settlement of a Millville silt loam obtained in a study of wetting-induced changes in near-surface soil physical properties. Soil columns with a height and diameter of 200 mm were subjected to two 30-min wetting and 30-min drainage cycles. The settlement data were obtained by monitoring the vertical location of thin layers of colored chalk powder placed at 3-cm intervals in the soil during packing with respect to marks on the plexiglass column. The initial thickness of the thin layers was 30 mm. The data reported here are cumulative settlements of six thin layers for Millville soil aggregates (2–4 mm in diameter). Details of the experimental setup can also be found in Or (1996). In a parallel experiment, the PSD and SWC were determined with Tempe pressure cells where ethanol—a non-polar liquid—served as wetting fluid (Silva, 1995). Alcohol retention curves were determined on soil samples prior to the initial wetting-drying cycle and after each subsequent cycle. The alcohol retention curves were transformed to the SWC according to Aggelides (1987). We applied the coalescence model to quantify the one-dimensional settlement observed in the experiment using the axial strain of a unit element consisting of equally sized (0.06 mm) aggregates in a rhombic packing. This choice of aggregate size and packing are believed to be more realistic. The use of an aggregate size obtained by sieving and a cubic packing of regular spheres will overestimate the pore size of natural porous media. The calculations were carried out as described by Ghezzehei and Or (2000) for coalescence under a constant matric potential = 9475 Pa. Figure 8a shows the modeled relative axial strain function, (t), and values for obtained from settlement data. Most of the settlement took place during the drying cycles (i.e., 30–60 and 90–120 min).
View larger version (17K):
[in this window]
[in a new window]
|
Fig. 8. Parameters inferred from the settlement of a Millville silt loam during wetting (0–30 and 60–90 min) and drying (30–60 and 90–120 min) cycles: (a) observed and optimized axial strain and (b) drift coefficient according to the coalescence model.
|
|
The FPE coefficients were calculated with the coalescence model using a unit cell comprising of spherical aggregates in rhombic packing (compare Eq. [28] and [29]). The pore radius and porosity of unit cell for rhombic packing are given by
 | [33] |
 | [34] |
Note that the experimental information does not allow us to quantify the FPE coefficients as a function of pore size. The drift coefficient (dr/dt) is given in Fig. 8b as a function of time. The model predicts an increased drift with time because dr/d tends to increase with for a rhombic unit cell. Coalescence continues until the model constraints are exceeded. Especially during the second drying cycle (t > 60 min), the conditions for coalescence quickly cease to exist. The drift and degradation coefficients were parameterized with exponential functions for use in the analytical solution Eq. [13]. The transformed time (Eq. [7]) and porosity are given by
 | [35] |
 | [36] |
where T is the optimized cumulative drift term (µm), t is time (s) and R2 denotes the coefficient of determination.
The independent retention measurements () allow us to assess the predicted PSD curves. The initial PSD is obtained from the SWC according to (Morel-Seytoux, 1969):
 | [37] |
with as the volumetric water content and as porosity while "<>" denotes an ensemble average. To evaluate the expression on the left-hand side of Eq. [37], we used the soil water capacity C (= d/d) and the Laplace equation. A bimodal PSD, which is based on Eq. [5], was fitted to this initial distribution for a structural component (r0 = 25 µm, calculated from a for rhombic packing, = 0.3, porosity of 0.227) and a textural component (r0 = 1.5:m, = 0.75, porosity of 0.315) The PSD at the end of the first and second drying cycle was calculated with Eq. [13] using the aforementioned initial parameters and a cumulative drift term as well as a dispersivity = 0.5 µm and a degradation term = 0.1. The values for the latter two variables are chosen based on earlier simulations. Figure 9a
displays the initial and the predicted curves for the PSD. The curves are similar to previous figures, the median and mean pore size decrease with time because of drift towards smaller pores as predicted by the coalescence model.
View larger version (22K):
[in this window]
[in a new window]
|
Fig. 9. Pore-size distribution and retention curves for the Millville silt loam: (a) initial PSD determined from soil water characteristics (SWC) and pore-size distributions (PSDs) predicted at the end of two drying cycles using the coalescence model, and (b) corresponding SWCs and experimental retention data.
|
|
Subsequently, the SWC was estimated from the PSD at the end of the drying cycles, also according to Eq. [37]. Figure 9b shows the predicted and initial SWCs as well as independently measured retention data. The predicted SWCs describe the retention data reasonably well, except that the model does not adequately account for the increased retention in the intermediate range of water contents (0.25 < <0.35). It is also possible to predict the evolution of the hydraulic conductivity. Or et al. (2000) used the changes in the porosity inferred from the soil settlement to calculate the evolution of the saturated hydraulic conductivity, Ks, according to the Kozeny-Carman equation. The unsaturated hydraulic conductivity can then be calculated from the predicted Ks and the SWC using, for example, the model by Mualem (1976).
|
SUMMARY AND CONCLUSIONS
|
|---|
Two analytical solutions have been presented for the evolution of the pore-size distribution (PSD) according to the model postulated by Or et al. (2000). For the first solution, the drift and degradation coefficients V and M were assumed to depend on time. For the second solution, drift and dispersion term were also dependent on pore size. The solutions were illustrated with a time-dependent drift term given by Eq. [14], which can be applied to describe experimentally observed (mean) pore sizes for a wide range of drift profiles. The mean pore size and variance of the PSD respectively became smaller and greater with time (Fig. 1). Distributions resulting from drift and dispersion terms that also depend on pore size are shown in Fig. 2. Especially for a small dispersion term, the PSD may become quite narrow over time with a decrease in the variance. Although we did not derive hydraulic properties from these distributions, it is clear that the changing shape of the PSD has direct implications for the retention and transmission of water and solutes by soils.
Analytical expressions were derived for the first three pore size moments, which convey information about the total number of pores, the mean pore size and the variance of the PSD. The expressions provide a tool to both analyze experimental data and to quantify key features of the PSD.
The solutions for the PSD were also applied with drift and degradation terms that are based on a previously published model for aggregate coalescence. The model has the advantage of being physically based, but it requires several soil mechanical properties and it is subject to constraints related to yield stress and soil matric potential. Both solutions were used to illustrate the possible evolution of the PSD for a Millville silt loam. In the example, the model predicted coalescence for slightly more than 100 min. For a time-dependent drift term, the mean pore size decreased and the variance increased. For the solution where V and depend also on pore size, the mean pore size decreased more rapidly while the variance initially increased but later on became smaller for a "compressed" PSD. We analyzed settlement data during drying and wetting of a Millville silt loam to obtain drift and degradation coefficients according to the coalescence model. The PSD that was predicted from these results was used to infer the SWC, which we feel showed a reasonable good agreement with experimental retention data.
We conclude by noting that there is a lack of detailed data on the evolution of the PSD with time and pore size to truly evaluate and refine our approach. Nevertheless, the presented expressions for the PSD and the moments constitute quantitative tools that may provide the impetus for more detailed experimental studies of the temporal behavior of soil pores and hydraulic properties.
APPENDIX:
MOMENT ANALYSIS
Moments of the PSD can be used to elucidate the evolution of the PSD or to analyze and quantify experimental results. Moments are defined by:
 | [A1] |
A lower boundary of r = 1 is used in conjunction with the solution for a drift term depending on both time and pore size. Normalized moments are obtained through division by the zero-order moment (i.e., Mn = mn/m0). It is customary to characterize the drift term with the first-order normalized moment while the variance is given by the second-centralized moment µ2 = M2 - M21. Expressions for moments were obtained by inserting the solution according to Eq. [13] or [24] into Eq. [A1] and integration with the Mathematica package (Wolfram, 1999).
Drift Varies with Time
The expressions for the zero-order moment of the solution by Eq. [13], which accounts for a drift term that may vary with time, is given by:
The mathematical approximation is motivated by the fact that if no pores are lost, as stipulated by the boundary conditions, the zero-order moment should be equal to the area under initial PSD. It can be shown that the term between braces in Eq. [A2] is approximately equal to unity. The expressions for the first- and second-order moments are:
Drift Varies with Time and Pore Size
The corresponding moments for the solution that also accounts for a linear dependency of drift and dispersivity on pore size (i.e., Eq. [24]) are given by:
|
ACKNOWLEDGMENTS
|
|---|
The partial funding by USDA/NRI through contract 97-00814, by the United States-Israel Binational Agricultural and Research Fund (BARD) through contract 2794-96, by the Army Research Organization through contract 39153-EV, and by the National Science Foundation through subcontract Y542146 (SAHRA, University of Arizona), as well as the partial support of the Utah Agricultural Experimental Station (UAES) are gratefully acknowledged.
Received for publication May 18, 2001.
|
REFERENCES
|
|---|
- Aggelides, S. 1989. Pore size distribution of soils as determined from soil characteristic curves using non-polar liquid. p. 59–68. Proc. of the workshop on soil compaction, consequences and structural regeneration processes. Avignon, Commission of the European Communities. Balkema Publishers, Rotterdam, the Netherlands.
- Ahuja, L.R., F. Fiedler, G.H. Dunn, J.G. Benjamin, and A. Garrison. 1998. Changes in soil water retention curves due to tillage and natural reconsolidation. Soil Sci. Soc. Am. J. 62:1228–1233.[Abstract/Free Full Text]
- Barry, D.A., and G. Sposito. 1989. Analytical solution of a convection-dispersion model with time-dependent transport coefficients. Water Resour. Res. 25:2407–2416.
- Bear, J. 1972. Dynamics of Fluids in Porous Media. Elsevier, New York.
- Brooks, R.H., and A.T. Corey. 1964. Hydraulic properties of porous media. Hydrology Paper 3. Colorado State University, Fort Collins, CO.
- Ghezzehei, T.A., and D. Or. 2000. Dynamics of soil aggregate coalescence governed by capillary and rheological processes. Water Resour. Res. 36:367–379.
- Ghezzehei, T.A., and D. Or. 2001. Rheological properties of wet soils and clays under steady and oscillatory stresses. Soil Sci. Soc. Am. J. 65:624–637.[Abstract/Free Full Text]
- Hara, T. 1984. A stochastic model and the moment dynamics of the growth and size distribution in plant populations. J. Theor. Biol. 109:191–215.[ISI][Medline]
- Hillel, D. 1980. Applications of soil physics. Academic Press, New York.
- Huang, K., and M.Th. van Genuchten. 1995. An analytical solution for predicting solute transport during ponded infiltration. Soil Sci. 159:217–223.
- Jury, W.A., and K. Roth. 1990. Transfer functions and solute movement through soil. Theory and applications. Birkhäuser Verlag, Basel.
- Kosugi, K. 1994. Three-parameter lognormal distribution model for soil water retention. Water Resour. Res. 30:891–901.
- Laliberte, G.E., and R.H. Brooks. 1967. Hydraulic properties of disturbed soil materials affected by porosity. Soil Sci. Soc. Am. Proc. 31:451–454.
- Leij, F.J., E. Priesack, and M.G. Schaap. 2000. Solute transport modeled with Green's functions with application to persistent solute sources. J. Contam. Hydrol. 41:155–173.
- Leij, F.J., T.H. Skaggs, and M.Th. van Genuchten. 1991. Analytical solutions for solute transport in three-dimensional semi-infinite porous media. Water Resour. Res. 27:2719–2733.
- Mapa, R.B., R.E. Green, and L. Santo. 1986. Temporal variability of soil hydraulic properties with wetting and drying subsequent to tillage. Soil Sci. Soc. Am. J. 50:1133–1138.[Abstract/Free Full Text]
- Morel-Seytoux, H.J. 1969. Flow of immiscible liquids in porous media. p. 455–516. In R.J.M. De Wiest (ed.) Flow through porous media. Academic Press, New York.
- Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:593–622.
- Nimmo, J.R. 1997. Modeling structural influences on soil water retention, Soil Sci. Soc. Am. J. 61:712–719.
- Or, D. 1996. Wetting-induced soil structural changes: The theory of liquid phase sintering. Water Resour. Res. 32:3041–3049.
- Or, D., F.J. Leij, V. Snyder, and T.A. Ghezzehei. 2000. Stochastic model for posttillage soil pore space evolution. Water Resour. Res. 36:1641–1652.
- Ozkan, G., and P.J. Ortoleva. 2000. Evolution of the gouge particle size distribution: A Markov model. Pure Appl. Geophys. 157:449–468.
- Risken, H. 1989. The Fokker-Planck Equation. Springer-Verlag, Heidelberg.
- Silva, H.R. 1995. Wetting-induced changes in near surface soil physical properties affecting surface irrigation. Ph.D. dissertation, Utah State University, Logan, UT.
- Su, N. 1995. Development of the Fokker-Planck equation and its solutions for modeling transport of conservative and reactive solutes in physically heterogeneous media. Water Resour. Res. 31:3025–3032.
- Thornley, J.H.M. 1990. A new formulation for the logistic growth equation and its application to leaf area growth. Ann. Bot. (London) 66:309–311.[Abstract/Free Full Text]
- van Es, H.M., C.B. Ogden, R.L. Hill, R.R. Schindelbeck, and T. Tsegaye. 1999. Integrated assessment of space, time, and management-related variability of soil hydraulic properties. Soil Sci. Soc. Am. J. 63:1599–1608.[Abstract/Free Full Text]
- van Genuchten, M.Th., and W.J. Alves. 1982. Analytical Solutions of the one-dimensional convective-dispersive solute transport equation. USDA Tech. Bull. No. 1661. U.S. Salinity Lab. Riverside, CA.
- Wolfram, S. 1999. The Mathematica book, 4th ed. Wolfram Media, Champaign, IL, and Cambridge University Press, New York.
- Zoppou, C., and J.H. Knight. 1997. Analytical solutions for advection and advection-diffusion equations with spatially variable coefficients. J. Hydr. Eng., ASCE 123:144–148.
This article has been cited by other articles:
|
|
|
|
 
L. R. Ahuja, L. Ma, and D. J. Timlin
Trans-Disciplinary Soil Physics Research Critical to Synthesis and Modeling of Agricultural Systems
Soil Sci. Soc. Am. J.,
February 2, 2006;
70(2):
311 - 326.
[Abstract]
[Full Text]
[PDF]
|
|
|
TOP
ABSTRACT
INTRODUCTION
