Published online 26 April 2005
Published in Vadose Zone J 4:255-263 (2005)
DOI: 10.2136/vzj2004.0076
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
Unsaturated Flow Through Spherical Inclusions with Contrasting Sorptive Numbers
Alex Furmana,* and
A. W. Warrickb
a Hydrology and Water Resources, University of Arizona, Tucson, AZ 85721. Now at Institute of Soil, Water, and Environmental Sciences, ARO, Bet Dagan 50250, Israel
b Water and Environmental Sciences, University of Arizona, Tucson, AZ 85721
* Corresponding author (alexf{at}agri.gov.il)
Received 11 May 2004.
 |
ABSTRACT
|
|---|
The analytic element method (AEM) is used to model unsaturated flow through a spherical inclusion of contrasting hydraulic properties. The steady state Richards' equation is combined with the Gardner model for unsaturated hydraulic conductivity to form the Helmholtz equation. The later is solved by means of the AEM. The background and inclusion materials are assumed to have different saturated hydraulic conductivities and different sorptive numbers; hence, the conditions are more general than treatments of spherical inclusions. Continuity of the interfacial head boundary condition leads to a nonlinear system of equations, whose solution requires an iterative solution. Analysis includes the effect of the hydraulic properties and of the background flux and evaluation of computational efficiency for contrasting hydraulic properties. To examine water contents, a methodology is presented for matching Gardner and van Genuchten parameters. The new solution is more realistic than the previous solution for a spherical inclusion with a sorptive number the same as the background but is computationally more tedious.
 |
INTRODUCTION
|
|---|
PREVIOUSLY, WE ANALYZED the unsaturated flow through circular and spherical inclusions (Warrick and Knight, 2002, 2004). This expanded the earlier work of Philip et al. (1989) and Knight et al. (1989) dealing with fully impervious inclusions relevant to descriptions of flow in regions proximate to macropores, tunnels, and cavities and around underground obstacles such as stones or structures. This relates to a series of solutions to the refraction (and reflection) problem, originally addressed in other fields by Maxwell and by Rayleigh (e.g., Strutt, 1871; Maxwell, 1873), and later by Philip (1998) for open parabolic shapes. Also, contributions are by Bear (1972), Kacimov (2004)(personal communication, August 2004), Frenkel (1949), Bystrov (1959), and Irmay (1964). These new papers demonstrated the use of the AEM to an unsaturated flow domain. Past applications of the AEM had been limited to steady state saturated flow, described by the Laplace equation. Examples include solutions for large number of inhomogeneities (e.g., Barnes and Jankovi
, 1999; Jankovi
and Barnes, 1999). Recently, the method was extended to solutions of the Helmholtz equation (rather than the Laplace equation) hence allowing development of solutions for multi-aquifers (e.g., Bakker and Strack, 2003), transient flow (through the use of Laplace transform; Furman and Neuman, 2003), and for unsaturated flow through inclusions (Warrick and Knight, 2002, 2004; Bakker and Nieber, 2004).
All the AEM solutions for unsaturated flow made use of the Gardner (1958) function (also known as the Philip model) to describe the relations between the hydraulic conductivity and the pressure head. This allows the transformation between Richards' equation to the Helmholtz equation (using the matric flux potential, or the so-called Kirchhoff potential). Warrick and Knight (2002)(2004) were the first to introduce the AEM to describe unsaturated flow. They considered circular and spherical inclusions with a uniform distribution of the Gardner coefficient
. Bakker and Nieber (2004) extended the solution by considering circular inclusions (and circular drains) and allowing contrasts in
. This is a valuable contribution by adding flexibility and realism. The contrasting
values do lead to an increase in computational difficulty because a dual set of coefficients is used, which doubles the number of equations to solve and, more importantly, leads to nonlinear relations which require an iterative solution.
The objective of the present study is to apply the AEM to the analysis of flow through spherical inclusions in an unsaturated domain with contrasting sorptive numbers (Philip, 1986). The analysis remains for a "Gardner soil" (Gardner, 1958) with unsaturated hydraulic function of the form K = Kiexp(
ih) with the pressure head, h < 0. For spherical inclusions as considered by Warrick and Knight (2004) the sorptive number
i was a single value for the entire domain; we now consider the more difficult problem with both Ki and
i changing between the background and the inclusion regions.
 |
THEORY
|
|---|
General
Consider the spherical inclusion of radius r1 in Fig. 1
. Within the background domain and the inclusion, the saturated hydraulic conductivities are K0 and K1, respectively, and the corresponding sorptive numbers (Philip, 1986) are
0 and
1. Thus, the hydraulic conductivity function is K = Ki exp(
ih) with i = 0 corresponding to the outer region (background) and i = 1 to the interior of the inclusion. The origin is defined in the center of the inclusion with z > 0 for the upward direction. Cartesian coordinates x, y, and z are related to r, the spherical radial coordinate (r2 = x2 + y2 + z2), and
, the zenith angle by
 | [1] |
 | [2] |
 | [3] |
Also, the cylindrical polar radius is related to the Cartesian coordinates by
2 = x2 + y2. The Darcian flow is
 | [4] |
with the matric flux potential equal to
 | [5] |
The first part of Eq. [4] is flow due to pressure difference and the second part is due to gravity. For these definitions, Richards' equation becomes
 | [6] |
Note that Eq. [6] does not contain the saturated hydraulic conductivity, Ki.
Equation [6] can be related to the Helmholtz equation,
 | [7] |
using Moon and Spencer (1961)(p. 27). The resulting solution follows by separation of variables:
 | [8] |
for the exterior of the inclusion (r
r1), and
 | [9] |
for its interior (r
r1). The separation between internal and external solutions results from the requirement for finite potential at infinity. H is the Helmholtz variable,
2 =
/2, superscript "+" indicates the region external to the inclusion interface, and "" the region internal to the inclusion interface. Si is given by Si = (
ir1)/2, P are associated Legendre polynomials, p is a separation variable, and i and k are spherical Bessel functions of the first and third kind, respectively, given by
 | [10] |
and
 | [11] |
with I and K being modified Bessel functions of the first and second kind. To solve the Helmholtz problem, a general solution should be obtained, and boundary conditions should be satisfied at the inclusion interface. The solution, which is in terms of H, can then be converted to terms of pressure head, h.
The normalization of the spherical function is not necessary, but useful for matching the boundary conditions (i.e., for the evaluation of head at r = r1 the ratio of the spherical Bessel functions is one; hence, it does not require the computation).
Note that to avoid controversy we use the term Helmholtz, and not modified Helmholtz, throughout this paper. The only difference between the two is whether
is a real number or imaginary, leading to either modified Bessel functions or ordinary Bessel functions with real arguments.
 |
BOUNDARY CONDITIONS
|
|---|
General
There are two types of boundary conditions to apply: external and internal. External boundary conditions require that the pressure head (and hence H) be finite at r
0 and at r
. This implies the use of only spherical Bessel functions of the first kind inside the inclusion and only spherical Bessel functions of the third kind outside the boundaries.
There are two requirements for internal boundary conditions: continuity of normal flux and continuity of pressure head across the inclusion interface. The condition for the continuity of the pressure head, h+= h becomes
 | [12] |
which for
0
1 is a nonlinear equation.
The other requirement, for continuity of normal flux implies
 | [13] |
which is a linear partial differential equation regardless of the relations between the external and internal sorptive numbers. The application of Eq. [12] and [13] at a number of interface locations is used to solve for the coefficients in the Helmholtz solution.
Note that for a multiple spherical inclusions, in addition to the change in the form of the general solution (see later in text), internal boundary conditions include mutual contributions of all inclusions in the entire domain.
Application
Head boundary
The matching of pressure head across the inclusion interface requires the fulfillment of Eq. [12]. This can also be expressed as
 | [14] |
For contrasting
's, this is a nonlinear equation and the solution must be iterative. We linearized the equation following Bakker and Nieber (2004) in the form
 | [15] |
where
+ is the value of
+ evaluated using an estimate of coefficients from a previous iteration.
The boundary condition is written for M equally spaced boundary nodes along the inclusion interface (m is the index for the nodes, starting at 0). For evaluation of the coefficients, Eq. [15] can be written as
 | [16] |
where
 | [17] |
 | [18] |
 | [19] |
Flux boundary
A second set of equations is written to ensure the continuity of normal flux across the inclusion interface. Expanding Eq. [13] yields
 | [20] |
where (') indicates derivative with respect to the argument (radius). This can be rewritten with the help of Abramowitz and Stegun (1964)(Eq. 10.2.21):
 | [21] |
The inner solution flux is found analogously:
 | [22] |
Written in terms of Eq. [16], we obtain
 | [23] |
 | [24] |
 | [25] |
where C1 = (K1
1)/(K0
0).
For a convergence criterion, we require that the maximal change of any coefficient will be smaller than 105. This does not guarantee accuracy of the solution but was found to be sufficient for most cases as will later be demonstrated. We used N = 20 terms in the series approximations and M = 100 for the discretization of boundary conditions. These values, although detailed error analysis is not reported here, were found to be more than sufficient for this case. Similar investigations were reported previously (e.g., Bakker and Nieber, 2004) for circular inclusions.
 |
EXAMPLES AND DISCUSSION
|
|---|
Effect of Contrast in 
The contrast in
is the major difference between the solution presented here and that of Warrick and Knight (2004). We start by comparing the pressure head distributions around spheres of contrasting conductivity and sorptive numbers. Two combinations of hydraulic conductivities and sorptive numbers are presented in Fig. 2
. Also presented are the equivalent cases with contrast in hydraulic conductivity but no contrast in the sorptive number. Properties for all figures are listed in Table 1. Values are in consistent length and time dimensions.
It is clear that the contrast in
has a significant influence on pressure head distribution within the spherical inclusion. A contrast in
translates to a different range of actual hydraulic conductivities within the inclusion. This may lead to a complete reversal of the hydraulic properties (e.g., soil with higher saturated hydraulic conductivity may perform as a flow resistant body).
Examining the shape of the head contours shows that a higher sorptive number externally leads to a higher curvature of the pressure head contours within the inclusion (compare Fig. 2a with 2c and Fig. 2b with 2d). However, changes in pressure head distribution outside the sphere are much less dramatic. Also note that the range of the pressure heads, especially outside the inclusion, is very similar in both cases (i.e., homogeneous and heterogeneous sorptive numbers). This suggests that for some purposes an approximation that makes use of uniform
may be reasonable, gaining much higher computational efficiency.
Effect of Background Flux
The solutions of Warrick and Knight (2002)(2004) assumed a single value of the sorptive number. Consequently, an inclusion with higher saturated hydraulic conductivity will remain more conductive at any range of water content distribution. In our case, however, since conductivities vary at different rates, it is possible to have a reversal of conductance. One way such a reversal may happen is by changing the background flux so that the range of pressure heads at the inclusion vicinity will change. Figure 3 depicts the pressure head distributions for the same inclusions presented in Fig. 2, but with a fivefold increase in the background flux. Comparing Fig. 2b to Fig. 3f shows that such a reversal has indeed occurred.
It is important to note that some of the figures seem similar (e.g., Fig. 2a and Fig. 3e), suggesting the same distribution of pressure head (at different scales and ranges). But they are not identical. Warrick and Knight (2002)(2004) produced solutions that were conveniently presented by normalizing with respect to the background flux. Here, however, there is apparently no advantage of using a normalized form.
Comparison to Two-Dimensional Inclusion
Generally, it is easier to model (and visualize) flow problems in two dimensions. However, the dimensional reduction may sometimes be physically inappropriate. To examine the influence of the third dimension we compare the pressure head distribution caused by a two-dimensional inclusion, with the same distribution caused by a sphere of the same properties. The solution for the two-dimensional case of Warrick and Knight (2002) is now generalized, with two separate sets of coefficients (i.e., for the interior and exterior of the inclusion) and using the same linearization used here (see Eq. [15]) to accommodate contrasting sorptive numbers.
Pressure head distributions are presented in Fig. 4
. The perturbation due to the inclusion is greater in two dimensions than it is in the three-dimensional case. A circular inclusion translates into an infinite cylinder in three dimensions and accounts for the difference in flow proximate to the inclusion.
Stream Function and Flow Nets
The pressure head contours provide a full solution to the flow problem; however, an understanding of the flow phenomena requires a detailed examination of contour labels. For example, the contours external to the inclusion in Fig. 3e and 3g are very similar in shape, but close examination shows that while in Fig. 3e the pressure head is increased at the top of the inclusion (indicating higher water content), the opposite is observed in Fig. 3g. Furthermore, if uniform
is considered (e.g., compare Fig. 2c to Fig. 3g, and Fig. 2d to Fig. 3h) the contours are identical in shape and the values are scaled by
0.
To visualize better the flow phenomena, we use a dimensionless stream function, similar to that used by Warrick and Knight (2004). We define the dimensionless stream function as
 | [26] |
where Jz is the vertical flux, Jz,B is the background vertical flux, and u is an integration variable. The Jz is evaluated using
 | [27] |
To evaluate Eq. [27], equations 33 and 34 of Warrick and Knight (2004) are useful. Note however that since we use here a slightly different conversion between the Helmholtz function, H, and
(see Eq. [9]), there are corresponding changes [use of (
/
0) instead of (
/
0 1) when evaluating within the inclusion]. Resultant flow nets are presented in Fig. 5
.
Water Content
We choose to combine the Gardner model with the van Genuchten (1980) function to relate water content with pressure head. This requires matching of Gardner's
with the corresponding van Genuchten parameters
vg and nvg (nvg is related to mvg through mvg = 1 1/nvg). One way to perform the matching is through the calculation of the capillary length,
(
= 1/
).
 | [28] |
with Kdry, Kwet, hdry, and hwet being the hydraulic conductivities and pressure heads at the extremes of the matching range. Alternatively, the limits of integration can be taken from negative infinity to zero if a single sorptive number is needed for several applications. To match between the van Genuchten parameters and that of Gardner, we need to numerically evaluate the right hand side of Eq. [28] and directly obtain the sorptive number. This approach is useful as soil property databases (e.g., see Rosetta, Schaap, and Leij, 1998) tend to use the van Genuchten parameters rather than the Gardner one. Going in the other direction there are extra degrees of freedom, which typically requires fixing nvg by using literature values and evaluating
vg using Eq. [28]. Instead we choose here to minimize the sum of squares of differences between the hydraulic conductivity functions at fixed (50) number of points along in the range of pressure heads. The result is a unique set of van Genuchten parameters for every case presented and not a general match which might be biased at different ranges.
Table 2 lists the matched van Genuchten parameters used to calculate the water content. Note that since the matching is done locally (for the actual range of pressure heads), different numeral properties are obtained for each soil at different conditions.
For a complete calculation of the water content we also need extreme values (essentially residual and saturated water content). Obtained water content are dramatically different between the interior and the exterior of the spheres presented above. Figure 6 depicts the saturation distribution for a case with contrasting conductivity and sorptive numbers. Hydraulic conductivities and sorptive numbers for that case are 5 and 1 (conductivity) and 3 and 2 (sorptive numbers) for exterior and interior, respectively. Matched van Genuchten parameters are 1.162 and 0.745 (
vg) and 2.142 and 1.939 (nvg). Matching was performed for a pressure head range of (absolute values) 1.15 to 1.18. Note the jump in saturation across the inclusion interface.

View larger version (22K):
[in this window]
[in a new window]
|
Fig. 6. Distribution of effective saturation. See text for Gardner and matched van Genuchten parameters.
|
|
The Gardner function is somewhat less flexible than the van Genuchten function (and others) in describing the full range of the hydraulic conductivity function. However, if only a local range is considered (e.g., the maximal range of pressure heads in Fig. 3 of about 0.8), a very good match can be achieved. The changes in matched values presented in Table 2 show that increased accuracy can be achieved in localized matching. Compare, for example,
vg,0 for Fig. 5a (0.381) to
vg,1 for Fig. 5b (2.083). These two values correspond to the same soil but under dramatically different conditions.
Numerical Efficiency and Error Analysis
The major conceptual difference between the AEM and ordinary numerical techniques (e.g., finite differences, finite elements) is that the AEM fulfills exactly the partial differential equation while the other two only approximate it. In this case the Helmholtz equation is solved leading to Eq. [8] and [9]. This implies the use of the Gardner model and the Kirchhoff transformation.
The numerical nature of the AEM is primarily through the approximation of boundary conditions. The numerical error in the AEM model then is the highest at the vicinity of the boundaries and is affected primarily by the level to which the boundary is discretized. Another parameter which can affect the results is the truncation of the series. In the case of contrasting sorptive numbers there is also linearization of the boundary conditions, which may introduce additional error to the vicinity of the boundaries.
To quantify the error associated with the use of the AEM we make use of the maximal error at the interface. For alternative numbers of discretization nodes along the sphere interface, we compare both the head and the flux, calculated from opposite sides (i.e., from the inner and outer sides of the sphere). Table 3 lists the maximal error for all the examples listed above. In addition, the number of iterations required for each of the cases to converge is included. This gives an indication to numerical efficiency and the degree of linearization required.
The cases where the ratio
0/
1 is larger than one (i.e., fine inclusion in coarser soil) are significantly more difficult. The linearization for those cases requires large number of iterations, and the achieved error is up to 10 orders of magnitude larger than the error for the reversed cases (coarse soil in a finer one). However, note that even for those cases the error is a few orders of magnitude smaller than the actual values (e.g., compare maximal error of approximately 9 x 104 for head in Fig. 2b, with head values of [absolute value] h = 0.8).
Validity of the Gardner model
The Gardner model for describing the dependence of the soil hydraulic conductivity in the pressure head attracts modelers, primarily as it allows linearization of Richards' equation. However, its flexibility is limited compared to other models. This results typically in the inability to use a single Gardner model to match well for the whole range of the pressure head.
It is also important to reiterate that the Gardner model is valid only in unsaturated conditions. Once a single point in the domain (typically the top or the bottom of an inclusion) exceeds saturation (positive pressure head), the model becomes invalid. In such conditions the model domain should be separated into Laplace and Helmholtz domains. To the best of our knowledge this was not done yet in the framework of the AEM.
However, if only a limited range of pressure head is considered, the Gardner model can be appropriately matched with experimental results or with other models. The implication is that the soil characteristics are not represented by a single constant value, but by a value that depends on the range of pressure heads, roughly
=
(h). Note however that for the solution presented here to be valid the sorptive number needs to be constant in each domain; hence, the matching should be done for a range of pressure heads in a stepwise fashion, and not for a single value of h.
The technicalities for matching between the Gardner and the van Genuchten models through the Gardner capillary length were presented above (see under water content). Similar matching can be performed for other models. Alternatively one can use a form of the Gardner model that includes air entry pressure {i.e., use exp [
(h ha)]}, and use the air entry value ha as a matching parameter. Changes in the AEM algorithm are minor with this addition (e.g., Bakker and Nieber, 2004).
Solution for Multiple Inclusions
For the case of two dimensions, differences between the solution for a single inclusion and the one for multiple inclusions were limited to the mutual effect of inclusions on each other (and as a result, although not necessary, the preference to solve iteratively). There is generally a loss of symmetry with respect to the zenith angle. For the three-dimensional case, differences are analogous but there is a change in the form of the base solution (see Eq. [8] and [9]). The solution for a single inclusion practically reduces one of the dimensions (i.e., of
). This cannot be done for multiple inclusions; hence, a third component needs to be added to base solution (i.e., sine and cosine of
in addition to the spherical Bessel functions with argument of r and to the Legendre function with argument of
). In addition, the Legendre function takes the more complicated form with both degree and order.
 |
SUMMARY AND CONCLUSIONS
|
|---|
An AEM for a spherical inclusion in a uniform vertical flow field is presented. The model allows a contrast in both the hydraulic conductivity and the sorptive numbers. The solution requires the numerical matching of both pressure head and normal flux at sphere interface and is iterative due to nonlinearity in the head boundary condition. Solutions for multiple spherical inclusions of varying conductivities, sorptive numbers, and sizes can also be solved using the same techniques. This is more difficult, due primarily to the lack of symmetry and, thus, remains work in progress.
Several examples are presented for spheres of varying properties subjected to different background fluxes. With contrasting sorptive numbers, reversal of inclusion properties may occur for different background fluxes. The contrast in the sorptive numbers are shown to have dramatic effect at the interior of the sphere. A more limited effect due to the contrast are observed away from the inclusion, suggesting that for some applications the assumption of uniform sorptive number may be reasonable with a gain of computational efficiency.
A comparison is made between the solution for a sphere and the analogous solution for a single circular inclusion (in two dimensions). The solution for the circle has a greater effect on the flow, which is expected, although the pressure head patterns are similar. To compare water contents, we matched conductivity functions to the van Genuchten function and then made use of the water retention part of the latter. By localizing the matching (for the actual range of pressure heads), we obtained a much better agreement between the functions and, hence, increased confidence in the estimation of the water content.
 |
REFERENCES
|
|---|
- Abramowitz, M., and I.A. Stegun. 1964. Handbook of mathematical functions. Vol. 55. U.S. Gov. Print. Office, Washington, DC.
- Bakker, M., and J.L. Nieber. 2004. Analytic element modeling of cylindrical drains and cylindrical inhomogeneities in steady two-dimensional unsaturated flow. Available at www.vadosezonejournal.org. Vadose Zone J. 3:10381049.[Abstract/Free Full Text]
- Bakker, M., and O.D.L. Strack. 2003. Analytic elements for multiaquifer flow. J. Hydrol. (Amsterdam) 271:119129.
- Barnes, R., and I. Jankovi
. 1999. Two-dimensional flow through large numbers of circular inhomogeneities. J. Hydrol. (Amsterdam) 226:204210.
- Bear, J. 1972. Dynamics of fluids in porous media. Dover Publications, New York.
- Bystrov, K.N. 1959. On two-dimensional steady flows in a layer with exponentially varying thickness. (In Russian.) Trans. Moscow District Pedagogical Inst. 75:3159.
- Frenkel, Y.I. 1949. Theory of phenomena of atmospheric electricity. (In Russian.) Gostehizdat, Leningrad.
- Furman, A., and S.P. Neuman. 2003. Laplace-transform analytic element solution of transient flow in porous media. Adv. Water Resour. 26:12291237.[CrossRef]
- Gardner, W.R. 1958. Some steady state solutions of the unsaturated moisture flow equation with application to evaporation from a water table. Soil Sci. 85:228232.
- Irmay, S. 1964. Refraction of flow at the boundary between two different anisotropic porous media. (In French.) CRH Acad. Sci. 259:509511.
- Jankovi
, I., and R. Barnes. 1999. Three-dimensional flow through large numbers of spheroidal inhomogeneities. J. Hydrol. (Amsterdam) 226:224233.
- Kacimov, A.R. 2004. Capillary fringe and unsaturated flow in a porous reservoir bank. J. Irrig. Drain. Div. Am. Soc. Civ. Eng. (in press).
- Knight, J.H., J.R. Philip, and R.T. Waechter. 1989. The seepage exclusion problem for spherical cavities. Water Resour. Res. 25:2937.
- Maxwell, C. 1873. Treatise on electricity and magnetism. Vol. 1. Clarendon Press, Oxford.
- Moon, P., and D.E. Spencer. 1961. Field theory handbook. Springer-Verlag, Berlin.
- Philip, J.R. 1986. Linearized unsteady multidimensional infiltration. Water Resour. Res. 22:17171727.
- Philip, J.R. 1998. Seepage shedding by parabolic capillary barriers and cavities. Water Resour. Res. 43(11):28272835.[CrossRef]
- Philip, J.R., J.H. Knight, and R.T. Waechter. 1989. Unsaturated seepage and subterranean holes: Conspectus, and exclusion problem for circular cylindrical cavities. Water Resour. Res. 25:1628.
- Schaap, M.G., and F.J. Leij. 1998. Database related accuracy and uncertainty of pedotransfer functions. Soil Sci. 163:765779.[CrossRef]
- Strutt, J.W. (Lord Rayleigh). 1871. On the light from the sky, its polarisation and colour. Philos. Mag. 41:107120, 274279.
- 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]
- Warrick, A.W., and J.H. Knight. 2002. Two-dimensional unsaturated flow through a circular inclusion [Online]. Available at www.agu.org/journals/wr/. Water Resour. Res. 38(7). doi:10.1029/2001WR001041.
- Warrick, A.W., and J.H. Knight. 2004. Unsaturated flow through a spherical inclusion [Online]. Available at www.agu.org/journals/wr/. Water Resour. Res. 40. W05101. doi:10.1029/2003WR002890.
This article has been cited by other articles:

|
 |

|
 |
 
A. Furman
Modeling Coupled Surface-Subsurface Flow Processes: A Review
Vadose Zone J.,
May 27, 2008;
7(2):
741 - 756.
[Abstract]
[Full Text]
[PDF]
|
 |
|