Published online 8 March 2006
Published in Vadose Zone J 5:184-203 (2006)
DOI: 10.2136/vzj2005.0024
© 2006 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: FROM FIELD- TO LANDSCAPE-SCALE VADOSE ZONE PROCESSES
Stochastic Continuum Transport Equations for Field-Scale Solute Transport
Overview of Theoretical and Experimental Results
Jan Vanderborght*,
Roy Kasteel and
Harry Vereecken
Agrosphere, ICG-IV, Forschungszentrum Jülich GmbH, D-52425 Jülich, Germany
* Corresponding author (j.vanderborght{at}fz-juelich.de)
Received 11 February 2005.
 |
ABSTRACT
|
|---|
One-dimensional transport models that predict field-scale averaged solute fluxes are often used to estimate the risk of nonpoint source groundwater contamination by widespread surface-applied chemicals. However, within-field variability of soil hydraulic properties leads to lateral variation in local solute fluxes. When this smaller scale variability is characterized in a geostatistical sense, stochastic three-dimensional flow and transport equations can be used to predict field-scale averaged transport in terms of geostatistical parameters. We discuss the use of stochastic equations for the parameterization of equivalent one-dimensional models predicting averaged solute fluxes. First, we consider the equivalent one-dimensional convection dispersion model and the equivalent dispersivity, which characterizes the spreading of laterally averaged concentrations or solute fluxes. Second, we discuss the parameterization of a stream tube model to predict local transport variables (i.e., distributions of local concentrations and local arrival times) These local transport variables are shown to be important for predicting nonlinear local transport processes and useful for inversely inferring the spatial structure of soil properties. Stochastic flow and transport equations reveal a dependency of equivalent model parameters on transport distance and flow rate, which reflects the importance of smaller scale heterogeneities on field-scale transport. Approximate solutions of stochastic flow and transport equations are obtained for steady-state and uniform flow. The effect of transient flow conditions on transport is discussed. Throughout the paper we refer to experimental and numerical data that confirm or contradict results from stochastic flow and transport equations.
Abbreviations: BTC, breakthrough curve • CDE, convection–dispersion equation • pdf, probability density function • STM, stream tube model
|
INTRODUCTION
|
|---|
TO MANAGE GROUNDWATER QUALITY in response to land use and the widespread application of chemicals, the field scale can be considered the elementary management unit. Lateral variations in solute fluxes between different fields due to differences in management (e.g., land use), climatic conditions, and soil characteristics are important. Several studies have revealed that soil hydraulic properties can vary considerably within a field (Jury, 1985) as well. This leads to an effective redistribution of uniform water and solute fluxes at the soil surface and considerable variation of subsurface water and solute fluxes within a field (Biggar and Nielsen, 1976). Therefore, upscaling procedures are required to derive effective parameters describing the system's behavior and average fluxes at the field scale. Effective parameters are parameters that lump the system subscale heterogeneity and describe its behavior at a larger scale (Grayson and Blöschl, 2000; Harter and Hopmans, 2004). The upscaling problem, which exists at every scale, can be solved in two ways. In the so-called "scale way" approach (Vogel and Roth, 2003), the spatial arrangement of hydrologically relevant structures is explicitly considered in a small-scale process model. The predicted processes and variables at the smaller scale are then averaged, and effective parameters that predict the spatially averaged processes and variables at the larger scale are derived. Depending on the structure of the medium, the model type that is used to describe the process at the larger scale may be different from the smaller scale process model. In another approach, the average system's behavior is monitored at a larger scale and model parameters are derived using inverse modeling.
Reviews of transport models, their relation to soil structure, and their conceptualization of the transport process structure can be found in, for example, Feyen et al. (1998) and van Dam et al. (2004). The lateral averaging scale is generally assumed to be much larger than the scale of the heterogeneities, so the laterally averaged water flow process is uniform and predominantly vertical in the vadose zone. Since the scale of the solute application area is generally much larger than the scale of the heterogeneities, the field-scale averaged solute transport process can still be represented by a one-dimensional process. Despite the existence of a variety of other model concepts, the one-dimensional convection–dispersion equation (CDE) is often used to predict field-scale leaching of surface-applied chemicals through the vadose zone toward the groundwater, such as for pesticide registration (FOCUS, 2000; Tiktak et al., 2004). The predicted leached mass fraction of substances is sensitive to the dispersion parameter of the CDE, which lumps the effect of smaller scale transport velocity variations on solute spreading (Boesten, 2004). Therefore, the relation between smaller scale spatial variability in soil hydraulic properties and the dispersion parameter is important. Such a relation may be derived by characterizing the three-dimensional smaller scale structure of soil properties in a fully deterministic way and modeling three-dimensional flow and transport processes (Kasteel et al., 2000). As an alternative, the spatial variability of soil properties may be characterized in a stochastic or geostatistical sense, and stochastic three-dimensional flow and transport equations may be solved to upscale local-scale processes and parameterize a field-scale one-dimensional CDE.
Although it is common practice to predict field-scale transport and nonlinear reaction processes using a CDE that is parameterized on the basis of inert tracer experiments, it is nevertheless trivial that the spatial average of a local nonlinear process is different from its prediction using the averaged concentrations in the nonlinear process equation. For nonlinear local transport processes (e.g., nonlinear sorption), the statistical distribution of local-scale concentrations is required to upscale the process. A stream tube model (STM) that represents the field by a set of vertical soil columns in which one-dimensional convective–dispersive transport takes place, and which is parameterized by a distribution of advection velocities and dispersion coefficients (e.g., Bresler and Dagan, 1981; Amoozegarfard et al., 1982; Toride and Leij, 1996a, 1996b), is the simplest tool to predict the distribution of local-scale concentrations. As a consequence, a STM is an attractive model for upscaling nonlinear transport processes. The distribution of stream tube velocities and the stream tube dispersion, which quantifies the dilution of local concentrations, are obviously linked to the spatial variability of the hydraulic parameters and to local-scale dispersion. Solutions of stochastic flow and transport equations offer the possibility to make this link and parameterize the STM on the basis of spatial distributions of hydraulic and local-scale dispersion parameters.
Finally, the application of tomographic geophysical techniques for obtaining three-dimensional spatiotemporal solute concentration data sets is emerging (e.g., Hubbard and Rubin, 2004; Vereecken et al., 2004), which raises the question of how these data sets can be used to infer information about the spatial distribution of hydraulic properties. Again, solutions of stochastic flow and transport and transport equations can be used to link the structure of locally observed transport with the spatial distribution of hydraulic properties.
Overviews of solutions of these stochastic flow and transport equations are given by Dagan, (1989), Gelhar, (1993), Zhang, (2002), and Rubin, (2003). The major focus was on the prediction of expected or averaged concentrations of inert solutes or solutes undergoing linear reactions under steady-state and gravity-driven flow conditions. Our objectives are to give an overview of stochastic modeling of three-dimensional flow and transport in the vadose zone, its application for the parameterization of simpler one-dimensional transport models, its limitations, and how stochastic modeling agrees or disagrees with numerical and experimental observations. We first review the prediction of field-scale averaged concentrations of an inert solute and the parameterization of a "field-scale" or "equivalent" one-dimensional CDE. Next we discuss some recent and new results related to the prediction of local concentration distributions and the parameterization of a STM. The relevance of considering local concentration distributions for upscaling of nonlinear transport and for inferring the structure of local hydraulic properties is emphasized. In the final section, we discuss transport in nonstationary flow fields and the effect of transient flow regimes on transport in a heterogeneous soil profile.
|
PREDICTION OF FIELD-SCALE AVERAGED CONCENTRATIONS
|
|---|
Particle Trajectory and Particle Travel Time Variances
Transport of a nonreactive tracer in the three-dimensional flow field q(x) is described using the CDE:
 | [1] |
where x (L) is the coordinate vector, t (T) is time, 
(x) is the volumetric water content, C(x,t) (M L–3) is the solute resident concentration in the liquid phase, and Dd(x) (L2 T–1) is the local-scale dispersion tensor. Unsaturated flow is described by the Buckingham–Darcy equation:
 | [2] |
where q (L T–1) is the Darcy flux vector, x1 is the vertical coordinate (positive downward), h (L) is the matric head, and K(h;x) (L T–1) is the hydraulic conductivity (assumed isotropic). The matric head is obtained by solving Richards' equation:
 | [3] |
The soil hydraulic functions, (h;x) and K(h;x) vary in space. The deterministic spatial distribution of the hydraulic functions is generally not known, and their parameters are stochastic. The spatial variation of these parameters is expressed in a geostatistical framework often based on two assumptions. First, it is assumed that the hydraulic functions can be linearly scaled to reference functions K*(h*) and *(h*) (Warrick et al., 1977; Hopmans, 1987; Vogel et al., 1991):
 | [4] |
 | [5] |
 | [6] |
where 
K, , and h are the linear scaling factors. The second assumption is that the scaling factors, or their lognormal transforms, are represented by second-order stationary Gaussian random fields (i.e., the expected value and the two-point covariance of the random variable) and are defined and translation invariant.
 | [7] |
 | [8] |
where 
y
represents the expected value of the random variable y in all realizations of the random space function, and Cyy(h) is the two-point spatial covariance for a lag distance h.
Because of the stochastic scaling factors, the dependent variables of the flow and transport equations (, h, q, C) are stochastic as well. By solving the stochastic flow (Eq. [3]) and transport (Eq. [1]) equations, the expected values and spatial (cross)covariances of the dependent variables are derived from those of the hydraulic parameters. Therefore, a functional form of the parameter covariance is assumed, which here is a nonseparable exponential covariance function:
 | [9] |
where 
y2 is the variance and 
i (L) the spatial correlation length in direction i.
To solve the stochastic transport equation, the spatial variability of the water content and flow velocity must first be characterized in geostatistical terms. Approximations of the variance and spatial covariance of and q are obtained using a perturbation of the flow equation. Because of its mathematical convenience, the Gardner exponential hydraulic conductivity function is often used:
 | [10] |
where a(x) (L–1) parameterizes the influence of the capillary forces on the flow (Raats, 1976), aG is the geometric mean of a(x), and KsG is the geometric mean of the saturated hydraulic conductivity, Ks (L T–1). Inserting Eq. [10] in the Richards equation, Eq. [3] can be written as
 | [11] |
For uniform gravity-driven flow in an unbounded medium with constant parameters KG and aG, the mean matric head is constant in the flow domain, and the linearized perturbation equations of matric head and flux are given by (e.g., Harter et al., 1996)
 | [12] |
 | [13] |
where h' (L) and qi' (L T–1) are the perturbations of the matric head and local flux, respectively, and f = ln(K) is the perturbation of the log-transformed saturated hydraulic conductivity, and g = ln(1/h) is the perturbation of the log-transformed a(x) parameter. Multiplying the linearized perturbation equations (Eq. [12] and [13]) by a perturbation of a parameter or a variable and taking the expectation, a set of coupled partial differential equations (i.e., the moment equations) is obtained (Zhang, 2002, p. 249–250). The moment equations relate the (cross)covariance functions of the different stochastic variables and parameters. For second-order stationary variable and parameter fields (i.e., for a steady-state and a gravity-driven flow field in an unbounded domain), spectral methods offer a convenient way to solve the moment equations. Using these closed-form solutions, the spatial covariance of the water flux can be predicted directly from the spatial covariance and cross covariance of lnKs and 1/h or ln(1/h) (Yeh et al. (1985a, 1985b, 1985c; Russo, 1995a, 1995b; Harter et al., 1996).
To solve the stochastic transport equation, two different approaches have been followed. In the first approach, the Eulerian approach, the transport equation is perturbed and solved similarly to the flow equation. In the second approach, the Lagrangian approach, the resident and flux solute concentrations are interpreted as probability distributions of the locations and of the arrival times, respectively, of individual solute particles. The transport equation then describes the probability distribution of the location or arrival time of individual solute particles. The location at time t of a particle that was released at time t = 0 at point a, X(t;a) (L); the travel time of a particle from an injection point a to a reference surface perpendicular to the mean flow direction located at x1, 
(x1;a) (T); and the arrival time at location x of a particle that was injected in a surface at a1, (x;a1) (T) are related to the pore water velocity, v(x) = q(x)/(x) (L T–1) and a dispersive velocity fluctuation vd(t) (L T–1) due to local-scale dispersion as (Dagan, 1989 p. 277; Dagan et al., 1992; Rubin and Ezzedine, 1997; Vanderborght, 2001):
 | [14] |
 | [15] |
 | [16] |
where X–1(t;x) is the backward trajectory or the particle's location at a time t before arriving at x:
 | [17] |
The temporal covariance Cvdivdi(t) (L2 T–2) of vd(t) is related to the pore-scale dispersion tensor Dd as
 | [18] |
where 
(t) (T–1) is the Dirac function.
Since v(x) and vd(t) are stochastic, the particle locations and travel times are stochastic as well. The calculation of the particle locations and arrival times using Eq. [14], [15], [16], and [17] is not straightforward because the particle trajectory along which the velocity needs to be evaluated is not known. To solve this problem, the particle trajectory is approximated by the expected trajectory (e.g., Dagan, 1989, p. 311) and a deviation that accounts for the displacement due to local-scale dispersive velocities, xd(t) = 
vd(t')dt' (Fiori and Dagan, 2000). Although more general expressions for nonstationary transient velocity fields have been developed (see the section below on transport in nonstationary velocity fields), the further derivation largely simplifies when second-order stationary and steady velocity fields are assumed. Then, the expected trajectory and travel or arrival time (using a first-order Taylor expansion of Eq. [15] and [16]) can be written as:
 | [19] |
where v, the expected velocity, is constant in space and time. The first-order approximations of the particle trajectory X' and travel time ' deviations from their expected values are obtained after introducing Eq. [19] into Eq. [14] and [15] and using a first-order Taylor expansion of Eq. [15]:
 | [20] |
 | [21] |
where v'(x) is the deviation of the velocity at point x from its expected value. To simplify the notation, the injection point coordinate a was omitted. Similar equations can be written for the backward trajectories and arrival times. Using Eq. [20], the particle location variances at time t, xixj2(t) (L2), are related to velocity field moments as
 | [22] |
From Eq. [20] and [21], it follows that the particle location and travel time, 2(x1) (L2), variances are related as
 | [23] |
A straightforward derivation of the variances of the backward particle locations and arrival leads to identical relations. Therefore, we will not distinguish further between forward and backward trajectories and between arrival and travel times.
Equation [22] can be written in terms of the spectrum of the velocity perturbations, Svivj(k) (L5 T–2), which is the Fourier transform of the velocity covariance function Cvivj(x) (L2 T–2). After using Eq. [23], this results in:
 | [24] |
where k (L–1) is the wave number vector and i is the imaginary number. To arrive at Eq. [24], the second-order stationary velocity field was represented by its Fourier transform, and the characteristic function of the Gaussian random variable xd(t) – xd(t') with variance–covariance matrix 
xd(t) – xd(t') = 2Dd|(t – t')| was used. Equations [23] and [24] directly link the moments of the particle locations and travel times to the spatial variability of the hydraulic properties through the velocity spectrum and the local-scale dispersion, Dd.
In a second-order stationary velocity field, the moments of the particle locations and arrival times do not depend on the locations of the injection and observation points in the injection and observation surfaces (i.e., a and x do not appear in the right hand sides of Eq. [24]). The moments of the particle location and travel or arrival time probability density functions (pdf) are then related to the spatial and temporal moments of expected concentrations and solute fluxes, respectively:
 | [25] |
 | [26] |
where
 | [27] |
 | [28] |
and xixj2(t = 0) and 2(x1 = a1;a1) are the spatial and temporal moments of the initial solute plume and of the solute injection at the injection surface, Js (M L–2 T–1) is the solute flux, and c and cf the normalized resident and flux concentration.
The expected resident concentration or solute flux at a certain point x represent the average of the concentrations or fluxes at this point in all realizations of the velocity field, that is, the ensemble average. If there exists a set of points x with the same ensemble average of concentrations or solute fluxes, then the ensemble average can be exchanged with the average of the concentrations and fluxes in this set of points in one realization of the velocity field, that is, the spatial average. The equivalence between the ensemble and spatial averages is also called "ergodicity". For a wide initial solute slab or wide injection surface (i.e., much wider than the spatial correlation of the hydraulic conductivity field), the expected concentration and solute flux are independent of the location in a horizontal reference surface (i.e., perpendicular to mean flow direction): C(x,t) = C(x1,t) and Js(x,t) = Js(x1,t). As a consequence, the expected concentrations or fluxes then represent averaged concentrations or fluxes in horizontal reference surfaces. The same holds for the spatial and temporal moments of expected and horizontally averaged concentrations and fluxes.
One-Dimensional Equivalent Convection–Dispersion Equation
The spatial and temporal moments of the ensemble averaged or expected concentrations or fluxes can be used to parameterize an equivalent one-dimensional CDE that predicts distributions and breakthroughs of horizontally averaged concentrations and solute fluxes with the same spatial and temporal moments:
 | [29] |
where 
eq (L) is the equivalent dispersivity. Replacing C(x1,t) by Js(x1,t) or Cf(x1,t) 
Js(x1,t)/q1, solute fluxes or flux concentrations are predicted by Eq. [29]. The equivalent dispersivity is related to the travel distance and travel time variances (Jury and Sposito, 1985) as
 | [30] |
 | [31] |
The equivalent dispersivity eq is a measure of the spreading of a field-scale averaged concentration depth profile or breakthrough curve (BTC) due to local-scale dispersion and spatial variations in local advection velocities.
Through Eq. [23], [24], [30], and [31], which are obtained using first-order solutions of three-dimensional stochastic flow and transport equations, eq is linked to the spatial variability of the velocity field, the hydraulic parameter fields, and the local-scale dispersion. From Eq. [30] and [31] it follows that eq is a function of time or travel distance, whereas the one-dimensional equivalent CDE (Eq. [29]), presumes a constant eq. In Fig. 1
, the change of eq – d1 (d1 = Dd11/v1) (i.e., the part of the equivalent dispersivity that results from spatial variations of local advection velocities) with travel distance predicted by Eq. [24] and [31] is shown for saturated flow conditions and different structures of the hydraulic conductivity. For small travel distances, relative to the correlation scale of lnK in the vertical direction, 1, eq – d1 increases linearly with x1, and this rate of increase becomes smaller with increasing x1. For large travel distances, again relative to 1, eq – d1 and eq reach an asymptotic value. This asymptotic value is reached after a shorter travel distance, relative to 1 for structures that are more elongated in the direction of the mean flow (2,3 < 1), than for structures which are more elongated in the horizontal direction (2,3 > 1) (Russo, 1995a).
Since the travel distance in soils is generally relatively small (1–10 m), the asymptotic regime is rarely reached. The equivalent CDE (Eq. [29]) must then be interpreted as a tool to predict a BTC or concentration depth profile at a certain depth or time assuming a hydrodynamically uniform soil profile and using an equivalent dispersivity for this travel time or distance. If predictions for other travel times or depths are required, the one-dimensional CDE must be reparameterized. Using first-order approximate solutions of stochastic flow and transport equations, this reparameterization is based on the three-dimensional structure of the hydraulic parameter field. However, the assumption of a hydrodynamically uniform profile in the one-dimensional CDE is inconsistent with the large vertical gradient in biological and chemical soil properties in the upper 1 m of the soil profile. To what extend this inconsistency leads to incorrect predictions of the fate of substances undergoing chemical and biological reactions requires further investigation.
The effect of the local-scale dispersion Dd on the solute spreading caused by spatial variability of advection velocities (i.e., eq – d1) is interesting. With increasing Dd, eq – d1 decreases (Fiori, 1996). This decrease is caused by local-scale mixing that reduces the impact of local-scale velocity variations on the spreading of a solute plume. However, the spatial scale of dispersive mixing is considerably smaller than the scale of advection velocity fluctuations, so the impact of Dd on the travel time and particle location moments is small. The ratio of these two scales is defined by a Peclet number:
 | [32] |
In unsaturated media, the variance and spatial correlation of the velocity depend on the average degree of saturation. When lnK and lnh are negatively correlated (e.g., in a Miller–Miller similar medium), the unsaturated conductivity decreases more rapidly with increasing |h| for large than for small Ks. This leads to a crossing of hydraulic conductivity curves at a critical matric head |hcrit|. For 0 < |h| < |hcrit| the variance of the unsaturated conductivity decreases with increasing |h|, whereas the opposite occurs for |hcrit| < |h|. The location of the more conductive regions and the regions with higher flow velocities also switch when the critical matric head is passed (Roth, 1995). For non- or positively correlated lnK and lnh, the variance of the hydraulic conductivity monotonically increases with increasing |h|. Besides an effect on the variance of the unsaturated conductivities, the degree of saturation also influences the structure of the flow field. Under drier conditions, flow is more channeled into "preferential flow paths" that persist for a longer distance than under wetter conditions.
Because of these effects of the mean matric head on the variance of the conductivity and the structure of the flow field, which are related to the mean degree of saturation or the mean flow rate, eq – d1 depends on these mean state variables as well. An important consequence is that eq must not be considered a material constant. First-order approximate solutions of the flow and transport equations are tools to express relations between state variables, the spatial variance and covariance of the hydraulic parameters, and eq, at least in a qualitative sense (e.g., Russo, 1995a, 1995b; Harter et al., 1996). For negatively correlated lnK and lnh, eq – d1 reaches a minimum when |h| reaches |hcrit| and eq – d1 increases with increasing difference between |h| and |hcrit|. For non- or positively correlated lnK and lnh, eq – d1 increases with |h|. The travel times or distances before the asymptotic regime is reached also increase with increasing |h|. Numerical simulations of flow and transport in heterogeneous media (Roth, 1995; Harter and Yeh, 1996, 1998; Roth and Hammel, 1996; Birkholzer and Tsang, 1997; Hammel and Roth, 1998; Cirpka and Kitanidis, 2002; Khaleel et al., 2002) corroborate this behavior of eq – d1, with the degree of saturation predicted by first-order approximations. However, due to the strong nonlinearity of the hydraulic conductivity, the variance and the importance of higher order moments of the unsaturated hydraulic conductivity distributions strongly increase with increasing dryness (de Rooij et al., 2004) so that first-order approximations become invalid. But, the average flow rate is very small for dry conditions, so the spatial scales of diffusive mixing and advection velocity fluctuations become similar; that is, the Peclet number (Eq. [32]) decreases. Therefore, local-scale diffusive fluxes may effectively smooth out concentration differences between regions with different velocities and limit spreading caused by advection velocity variability. As a consequence, the increase of eq – d1 with increasing dryness may be considerably reduced (Cirpka and Kitanidis, 2002).
Although the relation between the spatial variability of hydraulic properties and solute spreading in unsaturated soils has been the topic of several studies, the theoretical and numerical studies largely outnumber studies providing direct experimental support of these relations. Reviews of tracer experiments in soils (Beven et al., 1993; Vanderborght et al., 2001) provide indirect evidence that flow and transport in soils is indeed a spatially variable process that leads to an increase of equivalent dispersivity with travel distance. Also the effect of the flow rate or degree of saturation on eq has been demonstrated (Forrer et al., 1999; Vanderborght et al., 2001). However, except for tracer experiments in unsaturated heterogeneously packed sand containers (Wildenschild and Jensen, 1999; Ursino et al., 2001b), eq was not found to increase with decreasing flow rate or decreasing degree of saturation. This suggests that leaching experiments were mostly performed using relatively high flow rates under relatively wet conditions, that is, in the range where, for negatively correlated lnK and lnh, eq increases with saturation or flow rate.
The difficulty of characterizing the variance and spatial structure of the hydraulic conductivity may be the principal reason for the scarcity of direct experimental links between observed transport and spatial variability of soil hydraulic properties. The spatial variability of the hydraulic properties may vary considerably depending on the experimental method used to determine these properties (Mallants et al., 1997) and the method used to derive the parameters of the spatial covariograms (Ünlu et al., 1990). The linear scaling approach as presented in Eq. [4], [5], and [6] may not completely represent the spatial variability of the unsaturated hydraulic functions (Deurer et al., 2000). The estimation of the spatial correlation length may be problematic when the scale of the spatial correlation and the distance between soil samples is of the same order of magnitude. This is especially so for the correlation length in the vertical direction which is in the order of 10–1 m (Russo and Bouton, 1992; Rockhold et al., 1996; Russo et al., 1997). It should be noted that field sites where spatial correlation lengths of hydraulic properties could be derived are located in more arid regions and often in poorly developed soils with a coarser texture (Entisols) (Ünlü et al., 1990; Russo and Bouton, 1992; White and Sully, 1992; Rockhold et al., 1996). However, a similar number of studies show variograms of hydraulic parameters with a large nugget, indicating that no spatial structure could be identified (Mohanty et al., 1994; Mallants et al., 1996; Kasteel, 1997; Vereecken et al., 1997; Hammel et al., 1999; Deurer et al., 2000).
Since the spatial correlation of the hydraulic conductivity is difficult to determine, especially in the vertical direction of the mean flow, an indirect estimation of this correlation scale from the scaling of eq with travel distance seems interesting. By observing the scaling of eq, derived from, for example, BTCs or concentration profiles at different depths or times, information about the structures of the flow and hydraulic parameters fields may be inferred inversely (Russo, 1997). Vanderborght et al. (1997) compared eq derived from transport simulations in generated heterogeneous conductivity fields with eq derived from measured field-scale averaged BTCs. On the basis of this comparison they identified an experimental method that yields relevant information about the variability of unsaturated hydraulic conductivity.
Besides problems with inferring statistical parameters of soil hydraulic properties, second-order stationary Gaussian fields may not be adequate to characterize the spatial variability of hydraulic properties in soils that are layered systems that developed as a result of a long-term leaching process. Hammel et al. (1999) investigated the effect of soil layering and wedges on transect-scale transport. The effect of larger scale structures was found to be relatively localized and depends on their size and hydraulic properties that define their potential to divert flow laterally. Another aspect is the connectivity of regions with higher hydraulic conductivity. Deurer et al. (2001) attributed the rapid leaching of the tracer in the subsoil to a high connectivity and convergence of high conductivity zones in the subsoil layer. This structure cannot be described by a Gaussian random field, which assumes that high and low conductivity zones are equally connected and that extreme values are not connected (Gomez-Hernandez and Wen, 1998). A high connectivity of high conductivity zones in soils can be understood as the result of long-term leaching, which plays a crucial role in the development of soil profiles. Also in Gaussian random conductivity fields, a network of connected high flow regions develops. When soil properties are changed by long-term leaching along this connected flow pattern, a connected pattern of soil properties develops. The effect of these connected, high-conductivity regions on transport under saturated conditions was shown to be substantial and to lead to an early breakthrough and a long tailing (Wen and Gomez-Hernandez, 1998; Zinn and Harvey, 2003). The effect of the connectivity properties of the medium on unsaturated flow was investigated by (Neuweiler and Cirpka, 2005).
|
PREDICTION OF LOCAL-SCALE CONCENTRATIONS
|
|---|
Local concentration predictions require consideration for basically three reasons. First, for risk assessment, the probability of exceeding a critical threshold concentration may be more relevant than the average or expected concentration. Second, chemical equilibria and reaction rates that depend in a nonlinear way on local concentrations cannot be simply scaled up by using expected or spatially averaged local concentrations in the equilibria and reaction rate equations. Finally, information about the spatial distribution of local concentrations may be used to infer the spatial structure of the hydraulic parameter fields. In this section, we define parameters that characterize BTCs and the spatial covariance of local concentrations and relate these parameters to the spatial covariance of the velocity field and local-scale transport parameters. Most of the results have been obtained for saturated steady-state flow conditions. However, if the spatial covariance of the flow field can be determined for unsaturated flow conditions (i.e., for steady-state and gravity driven flow), these results can be directly transferred to unsaturated flow conditions.
Prediction of local-scale concentrations basically comes down to the prediction of the dilution caused by local dispersive fluxes. How the dilution process can be coupled to the local-scale dispersion and the spatial structure of the hydraulic conductivity field was the topic of several studies on saturated zone transport (Kapoor and Gelhar, 1994; Kitanidis, 1994; Kapoor and Kitanidis, 1996, 1998; Zhang and Neuman, 1996; Dagan and Fiori, 1997; Andricevic, 1998; Fiori and Dagan, 1999, 2000; Pannone and Kitanidis, 1999). In heterogeneous velocity fields, the solute plume is distorted and the internal plume surface, across which local-scale dispersion can take place, is much larger than in a homogeneous flow field. As a consequence, the dilution is considerably larger in a heterogeneous than in a homogeneous flow field. When the local-scale dispersion process cannot keep up with smoothing out lateral concentration differences generated by the heterogeneous advection velocity field, the dilution will be smaller than the dilution that is predicted in a homogeneous flow field using an equivalent dispersion coefficient that predicts the spreading of a tracer plume. Therefore, neither the local scale nor the equivalent dispersion coefficient can be used to predict the dilution process in a heterogeneous flow field.
Stream Tube Model Parameters and Their Prediction
Several parameters have been proposed to quantify dilution; for example, Kitanidis (1994) introduced the dilution coefficient and reactor ratio. The STM, which represents the complex three-dimensional flow field by a bundle of one-dimensional stream tubes, is an attractive alternative tool to predict both spreading and dilution in a heterogeneous flow fields. In the STM, dilution is quantified by a stream tube dispersivity, s (Cirpka and Kitanidis, 2000; Vanderborght and Vereecken, 2002). The locally observed breakthrough in a given realization of the velocity field, which is actually the result of a three-dimensional transport process, is interpreted to be the result of an equivalent one-dimensional convection dispersion process along a one-dimensional stream tube:
 | [33] |
where vs(x) (L T–1) and s(x) (L) are the stream tube velocity and dispersivity, respectively, and 
is the stream tube coordinate. When the stream tube is assumed to be a straight tube, corresponds with x1. Since the locally observed BTCs are spatially variable or stochastic, it is obvious that also the parameters vs(x) and s(x) are stochastic. The parameterization of the STM therefore boils down to defining the probability distributions of vs(x) and s(x). Using first-order solutions of stochastic three-dimensional flow and transport equations, the variance of vs(x) and the expected value of s(x) are derived.
Similar to the equivalent dispersivity, eq, the stream tube dispersivity s is defined in terms of the variance of the solute arrival time at a certain point in a given realization of the velocity field 2(x)|v(x) (T2):
 | [34] |
In contrast to 2(x), which represents the variance of arrival times at point x in all realizations of the velocity field and quantifies the spreading of a spatially averaged BTC, 2(x)|v(x) represents the variance of arrival times at point x in one realization of the velocity field and quantifies the spreading of a local BTC.
The expected variance of arrival times at a certain point in a given realization of the velocity field, 2(x)|v(x) is derived from the arrival time variance at a reference surface, 2(x1) and the covariance of the arrival times of two different particles that cross the reference surface in a given realization of the velocity field at the same location x, (x,x) (T2):
 | [35] |
A derivation of Eq. [35] is given in the Appendix.
The two-particle arrival time covariance function (x,y) represents the covariance of arrival times of particles at locations x and y on a reference plane. It is predicted from the advection velocity perturbations v', and the random displacements of the two particles due to local-scale dispersion, xd(t) and yd(t), in a similar way as the travel time variance 2(x1) using a first-order approximation of the particle arrival time (Fiori and Dagan, 2000; Vanderborght and Vereecken, 2002):
 | [36] |
Similar to the travel time variance, Eq. [36] is written in terms of the velocity perturbation spectrum, Sv1v1(k), and the local-scale dispersion tensor, Dd:
 | [37] |
Here, the characteristic function of the Gaussian random variable xd(t) – yd(t') with a variance covariance matrix, xd(t) – yd(t'), which by virtue of the independence between xd(t) and yd(t'), equals the sum of the covariance matrices of xd(t) and yd(t'), 2Dd(t + t') is used. It should be noted that in the case of purely advective transport, the two-particle arrival time covariance for x = y, (x,x), is equal to the arrival time variance at a reference surface, 2(x1). Inserting Eq. [37] into Eq. [35] and [34] leads to the following first-order approximate prediction of the stream tube dispersivity (Vanderborght et al., 2005):
 | [38] |
From Eq. [37] it follows that the first-order approximation of the two-particle travel time covariance function (x,y) is proportional to the two-particle trajectory covariance (see Eq. [15] in Fiori and Dagan (2000) and Eq. [B5] in Fiori et al. (2002)). The first-order approximation of the stream tube dispersivity, which is defined here in terms of the variance of the arrival time distribution at a given location, is therefore similar to the first-order approximation of an effective dispersion coefficient, which is defined in terms of the derivative with respect to time of the second spatial moment of a point-scale injected solute (Dentz et al., 2000; Cirpka, 2002).
In Fig. 2
, the equivalent dispersivity, eq, and the stream tube dispersivity, s, derived from a flow and transport simulation in a three-dimensional heterogeneous conductivity field are plotted together with first-order predictions of eq and s vs. travel distance from the injection surface. For the considered heterogeneity, the first-order predictions correspond fairly well with the dispersivities derived from the Monte-Carlo simulations. Also stream tube dispersivities that were derived from locally measured BTCs in a sand tank with known heterogeneous structure could be reproduced using first-order approximate predictions (Eq. [38]) (Jose et al., 2004). This implies that the local-scale dispersion coefficient may be inferred by fitting Eq. [38] to experimentally determined s values (Vanderborght and Vereecken, 2002; Jose et al., 2004). Figure 2 also illustrates that s increases with travel distance but considerably more slowly than eq, and is considerably larger than d, except for reference surfaces close to the injection surface. This reflects the difference of the dilution and spreading processes and the larger dilution in heterogeneous flow fields compared with homogeneous flow fields.
View larger version (19K):
[in this window]
[in a new window]
|
Fig. 2. Variance of stream tube velocities, vs2 (red), and dispersivities that characterize: the spreading of an averaged breakthrough curve, equivalent dispersivity, eq (black), the spreading of locally measured breakthrough curve, stream tube dispersivity, s (blue), and the pore-scale dispersivity d (green), in reference planes at different distances from the injection surface in an heterogeneous aquifer with f2 = 1, 1 = 2 = 5 m, 3 = 1 m, dL = 0.1 m, dT = 0.01 m. Thick lines are first-order predictions; symbols refer to parameters derived from a numerical experiment in a realization of the conductivity field.
|
|
Cirpka and Kitanidis (2002) investigated the effect of the degree of saturation on the stream tube dispersivity in heterogeneous fields with a negative correlation between lnK and lnh. They assumed a uniform local-scale dispersion that did not depend on the flow rate or degree of saturation. For small infiltration rates, s and eq converged to small values. Under these conditions, the time scales of local diffusion and advection are similar, so local variations in advection velocities are smoothed out and the dilution and spreading processes occur at a similar rate. The behavior of s with increasing flow rate and degree of saturation is more complex. At higher flow rates, the time for local-scale diffusive lateral mixing decreases, lowering s, whereas the variability of the flow field increases, increasing s. How these effects counterbalance each other depends on the heterogeneity of the hydraulic properties, the correlation between lnK and lnh, and the relative importance of local-scale diffusion compared with hydrodynamic dispersion. Dye tracer patterns in a heterogeneous sand tank for different leaching rates (Ursino et al., 2001a), partly contradicted the numerical and theoretical results of Cirpka and Kitanidis (2002). The dilution indices of the dye patterns, which are similar measures for dilution as s, did not change considerably with decreasing flow rate despite a large increase of the flow variability and plume spreading with decreasing flow rate. This indicates that the pore-scale mixing becomes less efficient (i.e., d decreases) with decreasing degree of saturation, due to an increase in tortuosity for local-scale diffusion.
The variance of the stream tube velocities, vs2 (L2 T–2) is related to the two-particle arrival time covariance:
 | [39] |
From Eq. [37] and [39] follows that vs2(x) is constant in a reference plane perpendicular to the mean flow direction. This implies that vs2(x) represents the variance of stream tube velocities observed in a reference plane in one realization of the three-dimensional velocity field. The stream tube velocity represents the average velocity of the particles in a stream tube along their trajectory from the injection surface to the observation point. Close to the injection surface, vs2 equals the variance of the pore water velocity, v12. Due to convergence and divergence of stream lines and due to diffusive exchange of solutes between different stream lines, the velocity of a particle changes along its trajectory with a concurrent decrease of vs2 with increasing travel distance (Fig. 2). The decrease of vs2 with increasing travel distance is often called "lateral mixing" (Flühler et al., 1996). However, also in the case of purely advective transport (i.e., without diffusive mixing and dilution), the lateral mixing process is active and results from the three-dimensional heterogeneous conductivity field, which leads to lateral flow components and a convergence and divergence of flow lines. Therefore, the notion of lateral mixing should be decoupled from a dilution of local concentrations.
Since both vs2 and s change with travel distance, the STM should not be interpreted as an exact representation of the transport process, which implies the assumption of no lateral exchange of solutes between different stream tubes and a constant velocity in the stream tube. As a consequence, the STM as it is used in this context is not a tool to derive vs2 from the variability of the hydraulic conductivity and the water content (Dagan and Bresler, 1979; Bresler and Dagan, 1981; Destouni, 1992) or to predict the breakthrough of averaged solute concentrations at different depths, assuming that vs2 remains constant with depth (Jury, 1982; Jury et al., 1986). It is rather a tool to describe essential transport phenomena, such as spreading and dilution at a certain observation or prediction depth given that vs2 and s are known for that depth. First-order approximate solutions of the stochastic flow and transport equations offer the possibility to derive vs2 and s and their change with depth from the spatial structure of the hydraulic conductivity field and the local-scale dispersion. By changing vs2 and s with depth or travel distance, the lateral mixing between stream tubes and the change of particle velocities along their trajectories can be implicitly accounted for.
Application of a Stream Tube Model to Predict Field-Scale Nonlinear Transport
In general, the STM can be formulated as
 | [40] |
where C(x1,t;C0,P) is the concentration at distance x1 from the injection surface and time t predicted by the transport and reaction equations in a one-dimensional stream tubes for boundary condition C0 and a set of transport and reaction parameters P (i.e., including the stream tube velocity and dispersivity), and pdf(C0,P) is the pdf of the parameter set P and the boundary condition C0. Because the numerical solution of transport and reaction equations in a set of one-dimensional stream tubes is computationally less demanding than solving the three-dimensional transport and reaction equations, the STM may be used to predict field-scale reactive transport. Often, purely advective transport is assumed. For uniform chemical rate parameters the general STM then simplifies to (e.g., Dagan and Cvetkovic, 1996)
 | [41] |
where = x1/vs is the advective travel time along the stream line, Pr is the parameter set of the reaction equations assumed to be constant in space, and pdf(,x1) is the travel time pdf. For certain reaction types, closed-form solutions of C(x1,t;,C0,Pr) can be used to solve Eq. [41], assuming a certain functional form of pdf(,x1). For purely advective transport, the particle arrival time variance 2 is identical to the two-particle arrival time covariance (x,y) and related to the stream tube velocity variance by Eq. [39]. In another approach, a numerical solution of the one-dimensional advective-reaction equations in Eq. [41] is used. For instance, Finkel et al. (1999) and Malmstrom et al. (2004) implemented geochemical models in the stream tube approach to model multicomponent reactions and transport in heterogeneous flow fields.
However, the effect of local-scale mixing plays an important role for certain reaction types. A typical example is a bimolecular reaction between equally mobile interacting compounds. In that case, the mixing between the two compounds can only take place through local-scale dispersion processes. Cirpka and Kitanidis (2000) and Ginn (2001) showed that a stream tube dispersivity, s, must be included in the STM to model such reaction rates in a heterogeneous flow field.
To illustrate the difference between solute spreading and dilution of local concentrations and the consequence of this difference for the upscaling of nonlinear reaction processes in heterogeneous media, the leaching of substances undergoing a nonlinear sorption characterized by a Freundlich sorption isotherm (S = kfCnf, where S is the sorbed mass concentration [M M–1]) was simulated. The two-dimensional flow and transport equations were numerically solved in 10 realizations of a two-dimensional heterogeneous conductivity field using the TRACE (Vereecken et al., 1994) and PARTRACE (Neuendorf, 1997) codes. Figure 3
shows the concentration patterns of the inert and reactive tracers. The fronts of the nonlinearly sorbing substances show more pronounced fingers than the fronts of the inert tracer, and this fingering gets more pronounced with increasing nonlinearity. This suggests that spatial variability of the water flow has a larger impact on the transport of nonlinearly sorbing than on linearly sorbing tracers. Simulated and predicted BTCs of the reactive substances at two different depths in the heterogeneous flow field are shown in Fig. 4
. Breakthrough curves are predicted by an equivalent field-scale CDE and by an advective STM. In the former, the spreading of an inert tracer is parameterized by eq, and the sorption is calculated from the field-scale averaged concentrations and the nonlinear sorption isotherm:
 | [42] |
where 
b (M L–3) is the bulk density. The equivalent dispersivity eq was derived for two observation depths using a first-order approximation of the stochastic flow and transport equations (Eq. [24] and [31]).
View larger version (18K):
[in this window]
[in a new window]
|
Fig. 4. Breakthrough curves of flux averaged concentrations at two depths, 0.5 and 0.8 m, below the injection surface in a heterogeneous two-dimensional saturated flow field (f2 = 1, 1 = 2 = 10 cm, dL = 0.2 cm, dT = 0.02 cm, s = 0.5) of two different tracers: a nonlinearly sorbing reactive tracer (kf = 1, nf = 0.67) (left panel), and a strongly nonlinearly sorbing tracer (kf = 0.6, nf = 0.3) (right panel). Solid lines are numerically simulated breakthrough curves in realizations of the conductivity field, dashed lines are the predictions by the equivalent CDE (Eq. [42]) and symbols are predictions by a STM (Eq. [41], [43], [44], [45]).
|
|
In the advective STM, the breakthrough in a stream tube at depth x1 is calculated using the method of characteristics:
 | [43] |
where tp (x1;,C0,Pr) is the time when the solute front arrives at depth x1 in a stream tube with travel time , and H is the Heaviside function. The depth of the solute front at a certain time, x1p(t;,C0,Pr), is derived from a mass balance equation:
TOP
ABSTRACT
INTRODUCTION
