Published in Vadose Zone Journal 3:513-526 (2004)
© 2004 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: UNCERTAINTY IN VADOSE ZONE FLOW AND TRANSPORT PROPERTIES
A Solute Flux Approach to Transport through Bounded, Unsaturated Heterogeneous Porous Media
Alexander Y. Sun*,a and
Dongxiao Zhangb
a Center for Nuclear Waste and Regulatory Analyses, Southwest Research Institute, 6220 Culebra Road, San Antonio, TX 78238
b EES-6, Los Alamos National Laboratory, Los Alamos, NM 87545
* Corresponding author (asun{at}cnwra.swri.edu).
Received 16 April 2003.
 |
ABSTRACT
|
|---|
In this paper, we present a solute flux approach for analyzing solute transport statistics in statistically nonstationary, unsaturated flow. Flow nonstationarity in the vadose zone may arise from a number of factors. It is useful to develop a systematic approach that incorporates these factors into an uncertainty analysis. We first derive the general forms for solute flux moments. The solute flux moments are associated with one- and two-particle joint probability distribution functions (JPDF). We illustrate our results for certain forms of one-particle and two-particle JPDFs, in which the particle travel time is assumed to be lognormally distributed and the particle transverse displacement normally distributed. In the numerical examples, the Eulerian velocity moments is obtained by solving the head moment equations numerically using a finite-difference method. Our results show that flow nonstationarity has a significant impact on the statistics of solute fluxes and solute breakthrough curves.
Abbreviations: CP, control plane • JPDF, joint probability distribution function • PDF, probability distribution function
|
INTRODUCTION
|
|---|
MODELING OF SOLUTE TRANSPORT in the unsaturated zone has received broad attention during the last two decades. Such analyses can provide valuable information to environmental risk assessment projects. For example, a travel time analysis is required in feasibility studies for potential nuclear waste repository sites. In case of leakage, the contaminant-laden solute flux may eventually reach the underlying groundwater aquifer and subsequently pose a direct threat to humans. It is thus important to assess beforehand the risks associated with such events. Both the U.S. Nuclear Regulatory Commission and USDOE have specified site-selection criteria. For example, a disqualifying condition specified by USDOE is that "a site shall be disqualified if the pre-waste-emplacement groundwater travel time from the disturbed zone to the accessible environment is expected to be less than 1,000 yr along any pathway of likely and significant radionuclide travel" (Altman et al., 1996). A numerical simulation performed by Altman et al. (1996) revealed that particle travel times from the proposed Yucca Mountain repository horizon to the water table vary from approximately 50 yr to more than 1 million years. This large variability in travel time was attributed to nonuniform infiltration rates at the ground surface and also to spatially heterogeneity of geohydrologic parameters.
Current approaches for modeling solute transport in a porous medium adopt either an Eulerian or Lagrangian framework. In Eulerian approaches, measurements of a quantity, say the solute concentration, are taken at fixed locations by preinstalled samplers. At any time the concentration field is characterized through measurements obtained at various locations. For transport in heterogeneous porous media, the concentration field is generally spatially and temporally nonstationary due to the evolving nature of mass transport. In Lagrangian approaches, the contaminant plume is envisioned as being composed of a large number of particles. The modeler tracks the movement of all particles. Thus, the Lagrangian approach offers a flexible and conceptually intuitive way for modeling contaminant transport. The derivation of a Lagrangian framework is independent of any specific velocity field. Consequently, once a Lagrangian framework is formulated, it can be used to analyze both saturated and unsaturated flow. One of the most difficult steps in Lagrangian analysis, however, is relating variables in a Lagrangian framework to those in an Eulerian framework, since by convention the velocity fields are given in the Eulerian form. Several approximations exist for such purposes (see a review by Rubin, 1997). In the past, geohydrologists used both Eulerian and Lagrangian approaches to derive macro-dispersion coefficients. In recent years, as Rubin (1997) pointed out, more recognition has been given to the true nature of the two approaches, with the Eulerian approach being used mainly in static type of analyses and Lagrangian approach in dynamic type of analyses. Examples of the latter are particle travel time and solute flux studies. Particle travel time problems are not new and have been an ongoing topic of research in other disciplines, such as in chemical engineering and quantum physics (van Kampen, 1992).
Shapiro and Cvetkovic (1988) derived moments of solute arrival times in terms of statistics of the hydraulic conductivity field. While the flow field is three-dimensional, the authors considered solute movement in the longitudinal direction only. They concluded that in the near field, the variance of arrival time, 
2
, is a quadratic function of the longitudinal travel distance. At large distances from the injection point (i.e., in the far field), 2 varies linearly with distance, an indication of Fickian transport behavior.
To increase the dimensionality of earlier solute flux analyses (e.g., Jury 1982; Simmons, 1982), Dagan et al. (1992) formulated a general Lagrangian framework for modeling conservative solute transport in three-dimensional, steady-state flow. The starting point of their analysis was to replace the particle trajectory vector by an alternative representation, which consists of the particle arrival time and the intersection coordinates of the particle with the control plane (CP). For a particle crossing the CP at location (x2 = 
, x3 = 
) and time , the particle trajectory equation, x = X(t; t0, a) is then replaced by t – t0 = (x1; a), (x1; a) = X2(; a) and (x1; a) = X3(; a), where t0 is the starting time and a is the particle origin. By doing so, the movement of each solute parcel is related to a macrostreamline (or a streamtube) and to its travel time along that streamline. The ensemble statistics of the solute fluxes, cumulative solute fluxes, and cumulative solute mass across the CP are then expressed with the aid of one- and two-particle travel time probability distribution functions (PDFs). As a special case, Dagan et al. (1992) considered steady uniform mean flow in an unbounded domain. They assumed lognormal distribution for one-particle travel time PDFs, and joint lognormal distribution for two-particle travel time PDFs in their examples (Cvetkovic et al., 1992).
For modeling transport in the vadose zone, Destouni (1992a)(1992b) derived moments of solute travel times and fluxes. She improved the Dagan–Bresler model (Dagan and Bresler, 1979) by considering the effect of vertical soil heterogeneity in her analysis. The mean travel time was found to increase linearly with depth, while the travel time variance exhibits a compression–expansion behavior for cases where both saturated and unsaturated conditions exist in the flow field; for "predominately unsaturated" flow, the system does not manifest this compression–expansion behavior. Vanderborght et al. (1998) conducted numerical studies to compare four different methods for obtaining the moments of travel time and transverse displacement. They also compared the effects of different concentration boundary conditions on travel time and transverse displacement statistics.
Virtually all studies mentioned above are based on the assumption of statistically stationary flow. This assumption is not appropriate in several situations. For example, it is not uncommon to encounter sites with layered formations, or with high or low conductivity lenses. In addition to medium nonstationarity, the mean flow nonuniformity can also be caused by external factors such as pumping or injection, or by internal factors such as the water table. The latter is important when considering flow in a shallow vadose zone.
Consequently, it is our purpose in this work to develop a Lagrangian framework for modeling solute transport in nonstationary flow domains. The results will provide insights into the effect of nonstationarity and allow one to assess, at least qualitatively, solute transport characteristics in nonstationary domains. Zhang et al. (2000) developed a Lagrangian framework for obtaining solute flux statistics in nonstationary saturated flow. Their work is extended in this paper to modeling transport in unsaturated flow. Spatial variability in the soil moisture content and soil hydraulic conductivity is accounted for in this study. In the following theoretical development, we first define the transport problem in a Lagrangian framework, and then derive general forms of the solute flux and travel time statistics. As a special case, we assume lognormal distribution for travel time and normal distribution for traverse displacement and then construct the one-particle and two-particle JPDFs accordingly. Finally, we illustrate our results through several numerical examples.
|
THEORY
|
|---|
The Lagrangian Framework
The problem formulation presented here is along the lines of Zhang et al. (2000). Consider a solute particle that starts moving from x = a at time t = t0. The particle travels along a streamline and intersects the CP at y = [y = (x; a), z = (x; a)] at time (x; a), where x is the longitudinal axis (the vertical axis in the vadose zone) and boldface letters represent vectors. Here the solute travel time, (x; a), is regarded as the first-passage time across the CP. In lieu of the Lagrangian particle trajectory equation, dX(t; a)/dt = V[X(t; a)] (Dagan, 1984), the particle trajectory can alternatively be described with the aid of travel time (x; a) and the streamline coordinates (x; a). The particle travel time (x; a) is obtained by inverting the relationship x = X1(t – t0; a), which yields t – t0 = (x; a). In addition, we have x2 = (x; a) and x3 = (x; a), which are simply the equations of streamlines and particle trajectories in a steady-state velocity field (Dagan et al., 1992). The system of ordinary differential equations that describe a particle trajectory then becomes (Dagan et al., 1992)
 | [1] |
with initial conditions (0; ax) = 0, (0; ax) = ay, (0; a1) = az, respectively, where Vi(x, , ), with i = 1, 2, 3, are components of the pore velocity vector, V(x, , ). Since the pore velocity is a random space function in a heterogeneous medium, the parameters and also become random functions.
Moments of Solute Flux
We denote the mass of the imaginary solute particle by 
m. If we know the initial solute concentration at a, say, C0(a) (M L–3), we then have the expression m = 
(a)C0(a)a, where a = axayaz is a differential volume surrounding a and (a) is the moisture content at a. The starting point of the analysis by Dagan (1982)(1984) is the following equation relating the concentration field to m:
 | [2] |
This equation states that C(x, t; a, t0) is equal to C when x 
X, and equal to zero elsewhere. Using the alternative set of parameters given in the previous subsection, the above equation becomes (Cvetkovic and Dagan, 1994)
 | [3] |
If we define qs = CV1, then Eq. [3] can be rewritten as
 | [4] |
where qs is the advection part of the solute flux [M (L2T)–1] that crosses the CP at time t. Here the contribution of pore-scale dispersion is omitted since the system is assumed to be advection-dominated. Equation [4] is simply a statement of mass continuity. In other words, one can reverse the order of derivation presented here by obtaining Eq. [4] first and then deriving the particle trajectory kinematics.
If the solute is reactive, Eq. [4] has to be revised to include a retention function, 
(t – t0, ), namely (Cvetkovic and Dagan, 1994):
 | [5] |
The above formulation is valid for linear reactions. For conservative solutes or tracers, (t – t0, ) becomes 
(t – t0 – ). In the case of continuous injection, Eq. [5] is modified to the following form (Andricevic and Cvetkovic, 1998; Zhang et al., 2000):
 | [6] |
where 
0(a) is the solute mass density at a, (y – ) is equal to (y – )(z – ), and the function 
(t – t0, ) is defined as
 | [7] |
where 
(t) (T–1) is the injection rate. For instantaneous injection, (t) = (t), and (t – t0, ) = (t – t0, ). Thus, among the three expressions for qs given in Eq. [4], [5], and [6], Eq. [6] is the most general one. We proceed with Eq. [6] in the following development. For simplicity, we assume t0 = 0.
With Eq. [6], the expected value of the solute flux can be obtained as
 | [8] |
where 
0 is the source volume, and f1[ = (x; a), = (x; a)] is the one-particle JPDF between and (i.e., f1 represents the probability of a particle, which originates at x = a, crossing the CP between and + and and + ). By using properties of the Dirac-delta function, we can simplify Eq. [8] to
 | [9] |
The total solute discharge Q(t; x) and the cumulative mass M(t; x) can be expressed in terms of the solute flux as follows:
 | [10] |
 | [11] |
The expected values of Q(t; x) and M(t; x) can then be obtained by performing an additional integration with respect to ,
 | [12] |
 | [13] |
Note that for conservative solutes, Eq. [13] in Dagan et al. (1992) can be retrieved from Eq. [9], [12], and [13] by setting (t, ) to (t – ).
The autovariance function of the solute flux is defined as
 | [14] |
where q's = qs – 
and the functions F1(·) and F2(·) are given by
 | [15] |
 | [16] |
where F2(·) is the two-particle JPDF; that is, f212 12 gives the probability of two particles, one originating at x = a and the other at x = b, crossing the CP between 1 + 1 
1 + 1 2 + 2 2 + 2.
One measure of the uncertainty associated with the solute flux prediction is the variance of the solute flux, which can be retrieved from Eq. [14] by setting t = t' and y = y'; that is,
 | [17] |
The autocovariances of Q(t; x) and M(t; x) can be derived in a similar manner.
One- and Two-Particle JPDFs
To obtain solute flux statistics, explicit forms of the JPDFs f1(, ; x, a) and f2(1, 1, 2, 2; x, a) have to be determined. In practice, this is not a trivial task since usually only the low-order moments of and are accessible, whereas an infinite number of moments might be required to completely characterize a PDF. As a consequence, many studies have opted to fit relatively simple distributions to the data, such as using the normal distributions, which are completely defined through their first two moments. Dagan et al. (1992) showed that for statistically stationary saturated flow, and are independent random variables based on the equations given in Eq. [1]. As a result, the one-particle JPDF, f1(, ; x, a), is separable; that is,
 | [18] |
where 0 = (ay, az) are the coordinates of the starting point in the transverse directions. For similar reasons, Dagan et al. (1992) showed that the two-particle JPDF is separable as well.
In the following development, we shall focus on two-dimensional flow and assume that the travel time is lognomally distributed and the transverse displacement normally distributed. The same assumption has been adopted in previous saturated zone studies (e.g., Zhang et al., 2000; Andricevic and Cvetkovic, 1998; Cvetkovic et al., 1992), and in vadose zone studies (e.g., Destouni, 1992a; Destouni and Graham, 1995; Russo, 1993). In this paper, we adopt the JPDFs presented in Zhang et al. (2000). For the sake of completeness we recap the major results in the next paragraph.
The JPDF between a lognormal random variable and a normal random variable is
 | [19] |
where = (x; a), = (x; a), ' = – <>, (ln)' = ln – <ln>, and r is a factor accounting for the correlation between and ,
 | [20] |
The mean and variance of ln are related to those of through the following pair of formulas:
 | [21] |
where CV = /<> is the coefficient of variation of . Note that the cross correlation between and is generally not zero for nonstationary flows. The two-particle JPDF is obtained similarly:
 | [22] |
where the residual vector X = [ln1 – <ln1>, 1 – <1>, ln 2 – <ln2>, 2 – <2> ]T, and the covariance matrix 
is given as
 | [23] |
The matrix is symmetric and must be positive definite. Some of the elements of are defined as
while other elements can be easily defined using Eq. [21]. In the next section, we derive approximations for all moments needed to evaluate the one- and two-particle PDFs presented here.
Moments of Travel Time and Transverse Displacement
The moments of (x; a) and (x; a) can be related to the statistics of the Eulerian velocity field. The starting point is the particle kinematics Eq. [1], from which we obtain the following equations by integrating with respect to x:
 | [24] |
where L denotes the distance of CP from the source. We now expand the velocity components in a Taylor series around the mean particle path, <>; that is,
 | [25] |
Decomposing the velocity components into the sums of means and zero-mean perturbations [i.e., Vi(x, ) = Ui(x, ) + ui(x, )] and substituting them into Eq. [25], we get
 | [26] |
 | [27] |
where Ui (i = 1, 2) are the means and ui are zero-mean perturbations. Notice that the nonstationarity in the velocity field is accounted for in the above expansions through the dependence of mean velocities on spatial coordinates and by the appearance of the last term on the right-hand side of Eq. [26] and [27]. Substituting Eq. [26] and Eq. [27] into Eq. [24] and making use of Taylor's expansion for 1/V1, we obtain to a first-order approximation (Zhang et al., 2000):
 | [28] |
 | [29] |
After taking expectations on both sides of Eq. [28], we obtain to a first-order approximation, the mean travel time,
 | [30] |
Subtracting Eq. [30] from Eq. [28] then yields the first-order perturbation term, '(L; a),
 | [31] |
where b1(x, 
) = 
U1(x, )/|=. Similarly, we obtain the mean and perturbation of (L; a); that is,
 | [32] |
and
 | [33] |
In its current form, Eq. [33] is an implicit equation since ' appears on both sides of the equation. Fortunately, we are able to simplify Eq. [33] to a more simple form. First, we rewrite Eq. [33] into an equivalent differential form:
 | [34] |
with the initial condition '(x = a1; a) = 0. The last two terms on the right-hand side of Eq. [34] are now treated separately. We can rearrange them as follows:
 | [35] |
As a first-order approximation, the mean pore velocity Ui can be replaced by the ratio of mean Darcy flux, Wi, and the mean moisture content, 
; that is,
 | [36] |
With the above formulation, the derivatives of Ui (i = 1, 2) with respect to are obtained as
 | [37] |
Substituting Eq. [37] into Eq. [35] and also making use of the following relationship obtained from Eq. [32]:
 | [38] |
we are able to rewrite Eq. [35] as
 | [39] |
In arriving at the first equality in the above equation, we have made use of the steady-state continuity equation, W1(x, )/x + W2(x, )/ = 0. In arriving at the second equality, we have made use of the chain rule. Substituting the final results of Eq. [39] into the original Eq. [34], and after some algebraic manipulations, we obtain
 | [40] |
Equation [40] is similar to the saturated zone result of Zhang et al. (2000)(Eq. [A6]) except that U1 is replaced by W1 on the left-hand side of the equation, and an additional term appears on the right-hand side of the equation. It is easily seen that the two forms are identical when is constant, which is the case in the Zhang et al. (2000) saturated zone study. Integrating Eq. [40] with respect to x yields
 | [41] |
The right-hand side of Eq. [41] contains the Darcy flux (W1), pore velocities (Ui, i = 1, 2), as well as the moisture content (). We proceed below using the form of '(L; a) given by Eq. [41].
We now derive the second moments for (L; a) and (L; a). The variance of (L; a) is
 | [42] |
in which 1 = 1(x1; a) and 2 = 2(x2; a). Similarly, for two particles originating at x = a and b, respectively, their travel time cross covariance 12 is given by
 | [43] |
where 1 = 1(x1; a) and 2 = 2(x2; b). Equations [42] and [43] contain two new terms that need to be defined, the cross-covariances 
and 
. We define them for the two-particle case using Eq. [41]:
 | [44] |
 | [45] |
where 
(x,<>) = U2(x,<>)/U1(x,<>). The above expressions again imply that 1 = 1(x1; a) and 2 = 2(x2; b). The variance of , 2, is a special case of Eq. [45]:
 | [46] |
Finally, the two-particle cross covariance between and is
 | [47] |
With the moments derived in this section, the JPDFs f1(, ; x, a) and f2(1, 1, 2, 2; x, a) can be completely defined under the assumptions made there. For the special case of uniform mean flow in the vertical direction, the coefficients b1(x,<>) and (x,<>) vanish, and we can retrieve the results obtained by Dagan et al. (1992) for stationary flow. For example, the variances 2 and 2 become
where Cuii(i = 1,2) are the velocity covariances. For nonstationary flows, the cross covariances given in this section are generally not zero.
|
NUMERICAL EXAMPLES AND DISCUSSION
|
|---|
Before calculating the solute flux statistics, we first need to obtain the statistics of the Eulerian velocity field and the travel time and transverse displacement moments. To obtain the velocity statistics, we shall use the finite difference code developed originally by Zhang (1999)(2002). This code solves the first-order moment equations of flow quantities, including the means and covariances of the capillary head and the moisture content, and the cross covariance between them. The results are then used to compute the statistics of the Darcy flux and the pore velocity. To compute the moments of the travel time and the transverse displacement, we use as inputs the velocity statistics and mean moisture content, and subsequently solve the particle kinematic equations using a fourth-order Runge–Kutta method (Press et al., 1996).
Effect of Water Table on 2 and 2
In this subsection, we illustrate the impact of flow field nonstationarity on travel time and transverse displacement variances. The flow nonstationarity referred to in this section is related solely to the water table boundary effect (Zhang and Winter, 1998; Sun and Zhang, 2000). When the water table is included, the vertical mean flow is no longer uniform, such as in the case of gravity-dominated flows. The latter assumes that the water table is located at infinity. As we will show in this section, the nonstationarity resulting from a water table boundary can significantly affect the flow and transport statistics.
In this study we assume that the vadose zone hydraulic properties can be well described by the Gardner–Russo model (Gardner, 1958; Russo, 1988) as given below
 | [48] |
 | [49] |
where 
(x) is the capillary pressure head, K(x) is the unsaturated hydraulic conductivity, Ks(x) is the saturated hydraulic conductivity, 
(x) is a soil pore-size distribution parameter, and e(x), s, and r are the effective, saturated, and residual moisture contents, respectively. We assume that s and r are deterministic constants due to their small variability. The parameter m in the exponent of Eq. [49] is related to soil tortuosity (Russo, 1988). For simplicity we set m to zero in the present study. The input parameters Ks(x) and (x) are assumed to be second-order stationary random functions with an exponential spatial covariance of the form
 | [50] |
where x' and x'' are two points in space, I represents the integral scale (or correlation length) of the corresponding random function, and 2 is the variance of the random function. A detailed discussion of the models and assumptions mentioned in this paragraph can be found in Sun (2000).
The two-dimensional domain in this example is assumed to be impermeable on both the left- and right-hand sides. For all cases, the domain is 20If in the vertical direction and 10If in the horizontal direction, where If is the integral scale of lnKs. The resolution of the numerical block is 0.5If x 0.5If. Values of the saturated and residual moisture contents are 0.3 and 0.0, respectively. A steady-state infiltration rate is applied at the upper boundary and the water table is set as the lower boundary. Thus, the mean velocity is unidirectional, but varies in the vertical direction due to changes in the mean moisture content.
The upper plot of Fig. 1
shows the travel time variance, 2, as a function of travel distance and for different values of 2a, where 2a is the variance of soil pore size distribution parameter, (x). For comparison purposes, the corresponding results for mean gravity-dominated flow are given in the lower plot in Fig. 1. Table 1 lists the parameter values used in these plots. The solid, dashed and dotted curves in Fig. 1 correspond to Cases 1, 2, and 3 in Table 1, respectively. In the presence of water table, the variance of 2 grows in a nonlinear fashion all the way to the water table. In contrast, with all other parameters fixed, 2 grows linearly at large distance when the mean flow is gravity dominated, and the magnitude is much smaller. Recall from Eq. [42] that 2 depends on the vertical velocity covariance, the mean vertical velocity, and the derivative of the mean velocity. Although the double summation of the longitudinal velocity covariance reaches an asymptotic value near the water table, the mean velocity (in the denominator) decreases gradually at the same time. This phenomenon was discussed in Sun and Zhang (2000). The net result is a rapid increase of 2 near the water table. This finding suggests that the travel time prediction at or near the water table can have a relatively large uncertainty associated with it. Approximating the flow as being gravity-dominated would not capture this trend in 2. In addition, we observe that 2 increases as 2a increases. In the figure, we make 2 dimensionless by multiplying it with 
2, where (U01 is the minimum vertical velocity (in this example, the velocity at the upper boundary) and If is the integral scale of lnKs.
Figure 2
shows 2 as a function of travel distance on three separate plots, each for a different 2a. The behavior of 2 is similar to that of the transverse particle trajectory autocovariance, X22, except that the independent variable in 2 is the travel distance instead of time. We observe that the magnitude of 2 is smaller in the presence of a water table than without a water table. Sun and Zhang (2000) observed the same trend in X22 through Monte Carlo simulations. To further explain this behavior, we see from the first-order approximation of 2 (Eq. [46]) that 2 depends not only on the reciprocal of the mean vertical velocity, but also on the transverse velocity covariance along the vertical direction, both of which increase toward the water table. The transverse velocity covariance is a hole-type function (Dagan, 1984) when the two points involved the covariance calculation do not coincide. The net area under the velocity covariance curve is positive, but its value for the stationary case is much larger than for the nonstationary case. Adding these net areas (one for each node) leads to the deviations shown in Fig. 2. Among the three cases, the largest CV of is 7.9%. To obtain dimensionless values, we divide 2 by (If)2.
Solute Flux Statistics
In this example we illustrate the solute flux statistics. To give readers a better feeling of the results, we use the soil data reported by Russo et al. (1997). These authors obtained in situ measurements of Ks and in a field plot 40 m long, 6 m wide, and 2 m deep. The site, which is located near Bet Dagan, Israel, is characterized by fine-textured soil. The Guelph permeameter method was used to obtain the soil hydraulic parameters. In that paper, Russo et al. estimated statistics of Ks and using data obtained from 130 observation points. The mean and variance of lnKs are –0.636 and 1.242, and the integral scale of lnKs is 22 cm. For this study, we converted ln statistics to the corresponding statistics using Eq. [21]. The resulting mean, variance, and integral scale of are 0.024 cm–1, 3.1 x 10–4 cm–2, and 6 cm, respectively. Since the measured parameter variances at the site were relatively high for a first-order analysis, we use 0.5 for the variance of lnKs and 1 x 10–5 cm–2 for the variance of . These choices were based on our experience gained from the previous first-order studies and Monte Carlo simulations (Sun and Zhang, 2000). Values of the saturated and residual water contents were not reported by Russo et al., and are assumed here to be 0.45 and 0.20, which are typical values for a sandy loam.
The domain is 20If deep by 10If wide. We used the water table as the lower boundary. Figures 3 and 4
show the solute flux statistics for two infiltration rates of 0.01 and 0.001 cm d–1. First, the mean head profiles are plotted in Fig. 3. The horizontal axis of the plot is normalized by the total thickness of the domain. Notice that for q0 = 0.001 cm d–1, the mean pressure head varies over nearly the entire domain, indicating that the water table boundary has a strong effect. For shallow vadose zones, other factors such as the water table fluctuations may also need to be taken into consideration. The term shallow vadose zone is only relative. Deep vadose zones also have nonstationary flow regimes. The difference is that in a deep vadose zone, a large portion of the domain may be gravity dominated, and one could apply statistical theories derived for stationary flow. The existence of a nonstationary flow regime may well alter the solute transport statistics near the water table. This should be a concern when studying the long-term fates of contaminants in the vadose zone.
Figure 4 shows the mean and standard deviation of the solute flux at a fixed time along the transverse direction (i.e., a snapshot of the solute plume in the transverse direction). The control plane location is 4If from the source, and the fixed time is taken to be the mean solute arrival time. The source is located near the surface and its horizontal dimension is 1If. To make the results dimensionless, we use three quantities, namely, U01, the minimum vertical velocity; If, the integral scale of lnKs; and M0, the initial solute mass. In this example, U01 is 0.004 cm d–1 and M0 is 20 g. The plots in Fig. 4 represent the superposed contributions of all streamlines originated from the source. Between the two cases shown in Fig. 4, the higher infiltration rate results in less spreading, more rapid transport, and thus a higher peak value (i.e., the mean travel time) at the observation time. Both the mean and variance of the solute flux have symmetric, Gaussian-like profiles. This is not surprising since the mean flow is uniform vertically and zero horizontally, and the transverse displacement distribution was assumed to be normal. If the mean flow is different from what is considered here, the flux profile may not have a Gaussian shape, even though the transverse displacement distribution is still assumed to be normal.
Figures 5 and 6
show standard deviations of the solute flux and the cumulative solute flux as a function of time and for two infiltration rates. In the lower part of both figures, we also show the CV, which is regarded as a better measure for uncertainty. The minimum value of CV occurs near the mean arrival time. It increases significantly in both directions, where the mean solute flux is generally lower. This result is known in the literature and was reported by, for example, Cvetkovic et al. (1992) in their saturated zone study. The difference is that the magnitudes of CV in the vadose zone are generally larger.
Examples of Medium Nonstationarity
Only the impact of water table–induced nonstationarity has been examined so far. Medium nonstationarity is another common factor that attributes to nonstationary flow. With the flexibility of the numerical code we use, it is relatively easy to devise an example that demonstrates the effect of medium nonstationarity, in addition to the nonstationarity caused by a water table. We consider the case of a high-permeability block being embedded in a low-permeability formation. Such a setting is conceptually analogous to that of a capillary barrier. The simplest capillary barrier consists of a fine soil layer (i.e., liner) overlying a coarse layer. Since the capillary pressure of the fine soil layer is lower than that of the coarser layer, water is kept away from the coarse layer with the help of a water collection system. Ho and Webb (1998) were concerned with the effectiveness of heterogeneous capillary barriers, which may occur naturally in many vadose zones. They conducted three numerical simulations, in which the materials composing the capillary barriers were assumed to be homogeneous, layered, and random heterogeneous, respectively. Ho and Webb (1998) found that the random heterogeneous capillary barrier performed the worst among the three because of excessive channeling within the barrier. Their conclusion implies that the natural capillary barriers often cannot be trusted. Instead, engineered barriers should be used whenever possible.
We would like to illustrate the working mechanism of a capillary barrier by constructing an example where a high-Ks block exists in the domain. In this example, the properties of the background low-permeability material are <> = 0.04 cm–1, 2a = 3 x 10–6 cm–2, s = 0.30, <lnKs> = –3, and 2f = 0.5. Table 2 lists the properties of the coarse-textured, high-permeability blocks used in Case I and II. For comparison purposes, we also present results for a case where no anomaly exists in the domain. Figure 7a and 7b
show plots of <> along the vertical axis for Cases I and II, respectively. The <> profiles indicate that the water is diverted away from the region of high permeability not only because of an increase in the transverse velocity resulting from lower moisture content, but also because of higher pressures in the coarse-textured block. Comparing Case I with Case II, we see that the degree of divergence in Case I is higher because the block in Case I has a higher mean log-conductivity, a larger <>, and a lower s.
View larger version (30K):
[in this window]
[in a new window]
|
Fig. 7. Plot of mean transverse displacement along the vertical axis when a high-conductivity block exists in less permeable background material. The mean lnKs of the background material is –3.0. (a) Case I for <lnKs> = 1.5, (b) Case II for <lnKs> = –1.0. Other parameter values are given in Table 2.
|
|
Figure 8
depicts the mean and standard deviation of solute breakthrough curves for the three cases being considered. The mean flux in Case I (solid line with squares) has the largest time span and the lowest peak. In contrast, the case having no block has the fastest peak arrival time and the narrowest time span. The standard deviation plot, however, shows that Case II has the lowest uncertainty at the peak arrival time. An explanation for this is that when medium variances are similar, the case of higher contrast (Case I) has larger variability than the case of lower contrast (Case II). In Fig. 9
, the solute fluxes are plotted along a control plane that is located at 92% of the total domain length. The plots were obtained by setting the time for each case equal to the mean arrival time along the central streamline. The curves shown on the figure hence correspond to different times. The spread of each breakthrough curve is as expected.
View larger version (20K):
[in this window]
[in a new window]
|
Fig. 8. Comparison of the mean and standard deviations of the solute breakthrough curves for Cases I and II.
|
|
|
CONCLUSIONS
|
|---|
A Lagrangian framework was established for analyzing solute flux statistics in statistically nonstationary, unsaturated flow systems. The starting point was the Zhang et al. (2000) analysis for solute transport in nonstationary, saturated flow systems. In this study, we added nonstationarity caused by variability in the moisture content. Particle movement was represented in a Lagrangian coordinate system that consists of the travel time of the particle to the CP, and the interception coordinates of the particle with the CP. The solute flux statistics, specifically, the first- and second-order moments, were then related to JPDFs of the solute travel time and the transverse coordinates. As a special case, solute transport in a two-dimensional steady flow field was considered. We assumed that the statistical distribution of the particle travel time was lognormal and the distribution of the particle transverse displacement normal. With these assumptions, the JPDFs can be entirely determined by the first- and second-order moments of and , and the cross covariance between the two. Using a first-order perturbation expansion, the moments of and were related to the statistics of the nonstationary Eulerian velocity field. Due to the complex nature of the unsaturated flow equation and the need to incorporate medium nonstationarity in our study, we elected to use a numerical solution.
As compared with mean gravity-dominated flow, the boundary effect induced by a water table causes the variance of the travel time to increase dramatically in the nonstationary regime. The same effect is responsible for a decrease in the transverse displacement variance near the water table. The mean solute breakthrough curve was found to slow down near the water table and have a larger span. The variance of the solute breakthrough curve manifested a similar behavior.
In addition to the capillary boundary effect, we also demonstrated the impact of medium nonstationarity. In the vadose zone, as a result of the spatial variability in the moisture content, the behavior of flow near anonymous blocks differs from that of the saturated flow. Our stochastic analysis can be used to assess solute flux statistics of layered soils, commonly encountered in the field. The first-order analysis followed in our study is limited to cases of mild heterogeneity. In future work we plan to test the robustness of our first-order results with Monte Carlo simulation.
|
ACKNOWLEDGMENTS
|
|---|
The authors wish to thank Dr. van Genuchten and the anonymous reviewers for their careful review of the manuscript.
|
REFERENCES
|
|---|
- Altman, S.J., B.W. Arnold, C.K. Ho, S.A. McKenna, R.W. Barnard, G.E. Barr, and R.R. Eaton. 1996. Flow calculations for Yucca Mountain groundwater travel time (GWTT-95). SAND96–0819. Sandia National Laboratories, Albuquerque, NM.
- Andricevic, R., and V. Cvetkovic. 1998. Reactive dispersion for solute flux in aquifers. J. Fluid Mech. 361:145–174.
- Cvetkovic, V., and G. Dagan. 1994. Transport of kinetically sorbing solute by steady random velocity in heterogeneous porous formations. J. Fluid Mech. 265:189–215.
- Cvetkovic, V., A.M. Shapiro, and G. Dagan. 1992. A solute flux approach to transport in heterogeneous formation. 2. Uncertainty analysis. Water Resour. Res. 28:1377–1388.
- Dagan, G. 1982. Stochastic modeling of groundwater flow by unconditional and conditional probabilities: 2. The solute transport. Water Resour. Res. 18:835–848.
- Dagan, G. 1984. Solute transport in heterogeneous porous formations. J. Fluid Mech. 145:151–177.
- Dagan, G., and E. Bresler. 1979. Solute dispersion in unsaturated heterogeneous soil at field scale. I. Theory. Soil Sci. Soc. Am. J. 43:461–467.[Abstract/Free Full Text]
- Dagan, G., V. Cvetkovic, and A. Shapiro. 1992. A solute flux approach to transport in heterogeneous formations: 1. The general framework. Water Resour. Res. 28:1369–1376.
- Destouni, G. 1992a. Prediction uncertainty in solute flux through heterogeneous soil. Water Resour. Res. 28:793–801.
- Destouni, G. 1992b. The effect of vertical soil heterogeneity on field scale solute flux. Water Resour. Res. 28:1303–1309.
- Destouni, G., and W. Graham. 1995. Solute transport through an integrated heterogeneous soil-groundwater system. Water Resour. Res. 31:1935–1944.
- 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:228–232.
- Ho, C.K., and S.W. Webb. 1998. Capillary barrier performance in heterogeneous porous media. Water Resour. Res. 34:603–609.
- Jury, W.A. 1982. Simulation of solute transport with a transfer function model. Water Resour. Res. 18:363–368.
- Press, W.H., B.P. Flannery, S.A. Teukilsky, and W.T. Vetterling. 1996. Numerical recipes. The art of scientific computing. Cambridge University Press, London.
- Rubin, Y. 1997. Transport of inert solutes by groundwater: Recent developments and current issues. p. 115–132. In G. Dagan and S.P. Neuman (ed.) Subsurface flow and transport: A stochastic approach. Cambridge University Press, Cambridge, UK.
- Russo, D. 1988. Determining soil hydraulic properties by parameter estimation: On the selection of a model for the hydraulic properties. Water Resour. Res. 24:453–459.
- Russo, D. 1993. Stochastic modeling of solute flux in a heterogeneous partially saturated porous formation. Water Resour. Res. 29:1731–1744.
- Russo, D., I. Russo, and A. Laufer. 1997. On the spatial variability of parameters of the unsaturated hydruaulic conductivity.
TOP
ABSTRACT
INTRODUCTION
