HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Web of Science (4) Citing Articles via Google Scholar Google Scholar Articles by Vanderborght, J. Articles by Vereecken, H. Search for Related Content PubMed Articles by Vanderborght, J. Articles by Vereecken, H. GeoRef GeoRef Citation Agricola Articles by Vanderborght, J. Articles by Vereecken, H. Related Collections Solute Transport Models Stochastic Processes Vadose Zone Processes and Chemical Transport

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

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

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION PREDICTION OF FIELD-SCALE... PREDICTION OF LOCAL-SCALE... TRANSPORT IN NONSTATIONARY... SUMMARY APPENDIX REFERENCES

 
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

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION PREDICTION OF FIELD-SCALE... PREDICTION OF LOCAL-SCALE... TRANSPORT IN NONSTATIONARY... SUMMARY APPENDIX REFERENCES

 
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

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION PREDICTION OF FIELD-SCALE... PREDICTION OF LOCAL-SCALE... TRANSPORT IN NONSTATIONARY... SUMMARY APPENDIX REFERENCES

 
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 D_d(x) (L² 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, x₁ 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 C_yy(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 _y² 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), a_G is the geometric mean of a(x), and K_sG is the geometric mean of the saturated hydraulic conductivity, K_s (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 K_G and a_G, 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 q_i' (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 lnK_s 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 x₁, (x₁;a) (T); and the arrival time at location x of a particle that was injected in a surface at a₁, (x;a₁) (T) are related to the pore water velocity, v(x) = q(x)/(x) (L T^–1) and a dispersive velocity fluctuation v_d(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 C_vdivdi(t) (L² T^–2) of v_d(t) is related to the pore-scale dispersion tensor D_d as

[18]

where (t) (T^–1) is the Dirac function.

Since v(x) and v_d(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, x_d(t) = v_d(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, _xixj²(t) (L²), are related to velocity field moments as

[22]

From Eq. [20] and [21], it follows that the particle location and travel time, ²(x₁) (L²), 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, S_vivj(k) (L⁵ T^–2), which is the Fourier transform of the velocity covariance function C_vivj(x) (L² 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 x_d(t) – x_d(t') with variance–covariance matrix _{xd(t) – xd(t')} = 2D_d|(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, D_d.

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 _xixj²(t = 0) and ²(x₁ = a₁;a₁) are the spatial and temporal moments of the initial solute plume and of the solute injection at the injection surface, J_s (M L^–2 T^–1) is the solute flux, and c and c^f 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(x₁,t) and J_s(x,t) = J_s(x₁,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(x₁,t) by J_s(x₁,t) or C^f(x₁,t) J_s(x₁,t)/q₁, 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 = D_d11/v₁) (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 ln_K in the vertical direction, ₁, _eq – _d1 increases linearly with x₁, and this rate of increase becomes smaller with increasing x₁. For large travel distances, again relative to ₁, _eq – _d1 and _eq reach an asymptotic value. This asymptotic value is reached after a shorter travel distance, relative to ₁ for structures that are more elongated in the direction of the mean flow (_2,3 < ₁), than for structures which are more elongated in the horizontal direction (_2,3 > ₁) (Russo, 1995a).

View larger version (17K): [in this window] [in a new window]   Fig. 1. Normalized equivalent dispersivity due to spatial variability of local flow velocity, (eq – dL)/(2f1), as a function of normalized travel distance, x1/1, in media with isotropic structures (1 = 2 = 3), elongated structures parallel (1/2,3 = 5,) and perpendicular (1/2,3 = 0.2) to the direction of the mean flow. Full lines: 1/dL = 500. Dashed line: 1/dL = 50.

The effect of the local-scale dispersion D_d on the solute spreading caused by spatial variability of advection velocities (i.e., _eq – _d1) is interesting. With increasing D_d, _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 D_d 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 ln_K and ln_h 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 K_s. This leads to a crossing of hydraulic conductivity curves at a critical matric head |h_crit|. For 0 < |h| < |h_crit| the variance of the unsaturated conductivity decreases with increasing |h|, whereas the opposite occurs for |h_crit| < |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 ln_K and ln_h, 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 ln_K and ln_h, _eq – _d1 reaches a minimum when |h| reaches |h_crit| and _eq – _d1 increases with increasing difference between |h| and |h_crit|. For non- or positively correlated ln_K and ln_h, _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 ln_K and ln_h, _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

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION PREDICTION OF FIELD-SCALE... PREDICTION OF LOCAL-SCALE... TRANSPORT IN NONSTATIONARY... SUMMARY APPENDIX REFERENCES

 
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 v_s(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 x₁. Since the locally observed BTCs are spatially variable or stochastic, it is obvious that also the parameters v_s(x) and _s(x) are stochastic. The parameterization of the STM therefore boils down to defining the probability distributions of v_s(x) and _s(x). Using first-order solutions of stochastic three-dimensional flow and transport equations, the variance of v_s(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 ²(x)|_v(x) (T²):

[34]

In contrast to ²(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, ²(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, ²(x)|_v(x) is derived from the arrival time variance at a reference surface, ²(x₁) 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) (T²):

[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, x_d(t) and y_d(t), in a similar way as the travel time variance ²(x₁) 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, S_v1v1(k), and the local-scale dispersion tensor, D_d:

[37]

Here, the characteristic function of the Gaussian random variable x_d(t) – y_d(t') with a variance covariance matrix, _{xd(t) – yd(t'}), which by virtue of the independence between x_d(t) and y_d(t'), equals the sum of the covariance matrices of x_d(t) and y_d(t'), 2D_d(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, ²(x₁). 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.

The variance of the stream tube velocities, _vs² (L² T^–2) is related to the two-particle arrival time covariance:

[39]

From Eq. [37] and [39] follows that _vs²(x) is constant in a reference plane perpendicular to the mean flow direction. This implies that _vs²(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, _vs² equals the variance of the pore water velocity, _v1². 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 _vs² with increasing travel distance (Fig. 2). The decrease of _vs² 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 _vs² 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 _vs² 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 _vs² 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 _vs² and _s are known for that depth. First-order approximate solutions of the stochastic flow and transport equations offer the possibility to derive _vs² and _s and their change with depth from the spatial structure of the hydraulic conductivity field and the local-scale dispersion. By changing _vs² 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(x₁,t;C₀,P) is the concentration at distance x₁ from the injection surface and time t predicted by the transport and reaction equations in a one-dimensional stream tubes for boundary condition C₀ and a set of transport and reaction parameters P (i.e., including the stream tube velocity and dispersivity), and pdf(C₀,P) is the pdf of the parameter set P and the boundary condition C₀. 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 = x₁/v_s is the advective travel time along the stream line, P_r is the parameter set of the reaction equations assumed to be constant in space, and pdf(,x₁) is the travel time pdf. For certain reaction types, closed-form solutions of C(x₁,t;,C₀,P_r) can be used to solve Eq. [41], assuming a certain functional form of pdf(,x₁). For purely advective transport, the particle arrival time variance ² 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 = k^fC^nf, 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 (52K): [in this window] [in a new window]   Fig. 3. Maps of simulated concentrations distributions 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 three different tracers: inert tracer (top panel), a nonlinearly sorbing reactive tracer (kf = 1, nf = 0.67), and a strongly nonlinearly sorbing tracer (kf = 0.6, nf = 0.3). Concentration maps are shown when the center of mass of the plume is at 40 cm below the injection surface (x1 = 60 cm in the plots).

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]).

[43]

where t_p (x₁;,C₀,P_r) is the time when the solute front arrives at depth x₁ in a stream tube with travel time , and H is the Heaviside function. The depth of the solute front at a certain time, x_1p(t;,C₀,P_r), is derived from a mass balance equation:

[44]

where M₀() (M L^–2) is the solute mass per unit area injected in a stream tube with travel time and

is the concentration depth profile in the stream tube at time t. For a uniform boundary condition (i.e., a uniform input concentration C₀ and application time t₀ for all stream tubes), the mass applied to a stream tube is proportional to the water flux in the stream tube. Assuming that the water content is the same in all stream tubes, the mass applied to a stream tube is related to the stream tube travel time as

[45]

where q₁ and are the average water flux and travel time, respectively. The depth of the solute front at time t in a stream tube with travel time is obtained by inserting Eq. [45] into Eq. [44] and solving it for x_1p(t;,C₀,P_r). By inverting the function x_1p(t;,C₀,P_r), the arrival time of the solute front, t_p(x₁;,C₀,P_r) is obtained. The breakthrough in a stream tube is calculated using Eq. [43] and the field-scale averaged breakthrough using Eq. [41]. A lognormal pdf was used to represent pdf(;x₁) and the stream tube travel time variance at a certain depth was derived from the first-order approximate two-particle travel time covariance.

The equivalent field-scale CDE clearly underestimates the arrival of the solute front (Fig. 4), especially when the degree of nonlinearity of the sorption isotherm is large (small n^f). In a heterogeneous flow field, the propagation of the averaged concentration at a certain depth represents the average propagation of local concentrations. For a nonlinear and concave sorption isotherm (n^f < 1) the local concentration propagation velocity increases with increasing concentration. The propagation of the average concentration is therefore considerably larger in a heterogeneous than in a homogeneous flow field. The advective STM performs much better since the propagation velocity is first calculated in the individual stream tubes before the field-scale averaged concentration is calculated. However, since local-scale mixing is not considered in the STM, the propagation velocities may be overestimated when local concentrations in the stream tubes are higher than in the heterogeneous flow field. Including _s in the STM may therefore improve the prediction. It remains to be investigated whether _s derived from inert tracer transport can be used to predict the dilution of a nonlinearly sorbing tracer in a heterogeneous flow field. Attinger et al. (2003) found that effective dispersion coefficients predicting the spatial moments of a point-scale injected solute plume in a single realization of the velocity field show the same behavior with travel distance for linear and nonlinear transport. This indicates that stream tube dispersivities derived from local linear transport could be applied to predict transport of nonlinearly sorbing tracers. But, also for purely advective transport, the tailing due to nonlinear sorption leads to a dilution of the concentration at the solute front. Therefore, _s may be neglected when this dilution is more effective than the dilution due to local-scale dispersion.

Similar to the breakthrough of the field-scale averaged concentrations, the equivalent field-scale CDE may fail to predict concentration depth profiles of nonlinearly sorbing substances in heterogeneous flow fields. Kasteel et al. (2002) demonstrated that the concentration depth profile of a nonlinearly sorbing dye tracer showed a larger spreading and a more pronounced leading tail than the predicted depth profile. The self sharpening effect of the invading dye tracer front that is expected in a homogeneous flow field due to nonlinear sorption was not observed in a heterogeneous flow field.

Spatial Structure of Local Concentrations and Stream Tube Velocities
The spatial distribution of local concentrations in a heterogeneous flow field reflects the structure of the flow field. It may therefore be used to infer information about the spatial structure of the hydraulic conductivity field. Forrer et al. (2000) determined spatial distributions of dye tracer concentrations with a 1-mm spatial resolution using quantitative image analysis of photographed dye patterns in soils. Using this technique, Kasteel et al. (2005) derived spatial correlations of dye concentration distributions and found that concentrations were correlated over only small distances (1–10 cm). This corroborates the small and hardly determinable spatial correlation scale of hydraulic properties (see discussion above) and points at much smaller advective than local-scale dispersive mixing time scale. The spatial correlation of the concentrations is predicted in terms of the geostatistical parameters of the hydraulic conductivity field and the local-scale dispersion (Pannone and Kitanidis, 2001; Vanderborght, 2001). For a wide injection surface and a pulse application of solute, the concentration covariance is approximated as

[46]

where M₀ (M L^–2) is the mass applied per unit area, and _x1y1(t,x – y) (L²) is the covariance of the displacement in the direction of the mean flow of two particles that are released in the flow field with a separation distance x – y. When x and y lay on the same distance from the injection surface (x₁ = y₁), then first-order approximations of _x1y1(t,x – y) and (x,y) are related as

[47]

The two-particle arrival time covariance, (x,y), is related to the spatial covariance of stream tube velocities as (Vanderborght et al., 2005)

[48]

Figure 5 shows first-order predictions of the spatial correlation perpendicular to the mean flow direction of local concentrations in the center of a solute plume (x₁ – v₁t = y₁ – v₁t₁ = 0) and of the stream tube velocities in a reference plane for two travel distances and two local-scale dispersion coefficients (Pe = ₁/_dL). The spatial correlations of local concentrations and stream tube velocities increase with travel distance or travel time and increasing local dispersion. This reflects the local dispersion process that smoothes local concentration gradients with time. The local concentrations are correlated over a smaller distance than stream tube velocities. This implies that a higher spatial resolution is required to determine the structure of local concentration patterns than the structure of the stream tube velocities. It also implies that the local-scale dispersion and the time available for mixing (i.e., the travel time) have a relatively larger impact on the correlation of local concentrations than on the correlation of stream tube velocities. For the considered travel distances and Peclet numbers, the spatial correlations of stream tube velocities are similar; that is, they are only slightly influenced or "contaminated" by the local-scale mixing. As a consequence, spatial correlations of stream tube velocities may be used to infer the structure of hydraulic properties, even when the local-scale dispersion is not known or assumed to be 0. The two-particle arrival time covariance and stream tube velocity covariance functions have been used to infer the spatial characteristics of the hydraulic conductivity in aquifers from BTCs measured at different locations on a reference surface (Rubin and Ezzedine, 1997; Vanderborght and Vereecken, 2001; Bellin and Rubin, 2004; Vanderborght et al., 2005) and in a multilevel pumping well (Bellin and Rubin, 2004).

View larger version (18K): [in this window] [in a new window]   Fig. 5. Effect of (a) the travel time/distance and (b) the local-scale dispersion (Pe = 1/dL) on the spatial correlation in the direction transverse to the mean flow direction of the local concentration, CC (open circles), and stream tube velocities, vsvs (solid circles). First-order predictions in a heterogeneous aquifer with f2 = 1, 1 = 2 = 5 m, 3 = 1 m, dL = 0.1 m, dT = 0.01 m are shown.

View larger version (74K): [in this window] [in a new window]   Fig. 6. Transport in a layered soil profile. Top panel represents the spatial distribution of the hydraulic conductivity in the profile [f2 = 1, dL = 0.2 cm, dT = 0.02 cm, s = 0.5, 1 = 2 = 10 cm (top layer), 1 = 20 cm, 2 = 2 cm (bottom layer)], middle panel the simulated vertical component of the pore water velocity, vx1, and the bottom panel the concentration profile when the center of mass of the plume is at 50 cm below the injection surface. In the areas between the two top and two bottom horizontal lines, the spatial correlation in the horizontal direction of the concentrations was calculated (see Fig. 7).

View larger version (15K): [in this window] [in a new window]   Fig. 7. Spatial correlations of (a) concentrations, CC, and (b) stream tube velocities, vsvs, in the direction perpendicular to the mean flow direction in a layered soil profile (see Fig. 6). Symbols refer to correlations derived from numerical simulations and lines from first-order predictions. Correlations of concentrations were derived at 30 cm below the injection (top layer; solid circles, dashed lines) and at 70 cm below the injection (bottom layer; open circles, solid lines) when the center of mass of the plume is at the boundary of the two layers. Correlations of stream tube velocities were derived at the 40-cm depth (top layer; solid circles, dashed lines) and at the 80-cm depth (bottom layer; open circles, solid lines).

	TRANSPORT IN NONSTATIONARY VELOCITY FIELDS

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION PREDICTION OF FIELD-SCALE... PREDICTION OF LOCAL-SCALE... TRANSPORT IN NONSTATIONARY... SUMMARY APPENDIX REFERENCES

 
The previous sections were based on the assumption of stationary velocity fields, which exist for steady-state and gravity-driven (h = 0) flow. Close to the soil surface, the flow regime is determined by transient climatic boundary conditions, whereas at the bottom boundary of the vadose zone, the hydraulic gradient in the capillary fringe above the groundwater table becomes <1. As a consequence, the velocity fields at those boundaries are nonstationary both in time and space. The solution of the stochastic flow and transport equations for transient flow conditions and in bounded domains is far more complex. Mantoglou and Gelhar (1987a, 1987b, 1987c) extended the spectral perturbation approach of Yeh et al. (1985a, 1985b) for transient flow conditions and found that local-scale heterogeneities induce anisotropy and hysteresis of the larger scale effective hydraulic conductivity, which depends on the prevailing flow conditions. These findings were confirmed by numerical flow simulations in a generated heterogeneous conductivity field (Polmann et al., 1991). However, the spectral approach imposes some restrictions (unbounded domains, smooth variations of flow characteristics in space and time, or quasi-stationarity) that limit its application for realistic cases of transient flow in soils. For general nonstationary conditions, linearized moment equations were solved numerically to obtain effective hydraulic parameters and two-point (cross)covariances of dependent variables (Mantoglou, 1992; Zhang, 1999). The effective hydraulic conductivity and the variance of water content and matric head derived from the stochastic theory were in good agreement with field observations (Jensen and Mantoglou, 1992).

The two-point covariance function of the pore water velocity is defined as

[49]

For nonstationary velocity fields, a first-order approximation of the particle trajectory equation (Eq. [14]) assuming purely advective transport (v_d = 0) (Indelman and Rubin, 1996) is

[50]

where

[51]

so that

[52]

For the case of unidirectional flow ( = 0 for i 1), the particle location variance at a certain time is derived from two-point velocity covariance as (Sun and Zhang, 2000):

[53]

Particle location variances for more general flow conditions are given by Lu and Zhang (2003). The nonstationarity of C_vivj(x,t;,t') implies that this function needs to be evaluated for each pair of velocity observations in the space–time domain. If the space domain is discretized using N_i and the time domain N_j nodes, C_vivj(x,t;,t') needs to be evaluated for (N_iN_j)² pairs. For three-dimensional flow and nonstationary boundary conditions, this evaluation becomes a daunting numerical task. This may explain that, although the approach can be applied for transient flow conditions, illustrative examples were given by Sun and Zhang (2000) and Lu and Zhang (2003) for steady-state two-dimensional flow fields. The effect of the capillary fringe above the water table on the compression of a solute plume (i.e., a decrease of the particle location variance with travel distance) could be reproduced using Eq. [53] (Sun and Zhang, 2000). Numerical simulations by Abdou and Flury (2004) illustrate the effect of the capillary fringe on BTCs. Especially in soils with vertical structures, the capillary fringe leads to a retardation of the breakthrough and a smaller coefficient of variation of travel times. However, in soils with horizontal and isotropic structures, the effect of the capillary fringe on the travel time moments was considerably smaller.

Russo et al. (1994a, 1994b) performed numerical experiments of solute transport in generated Gaussian heterogeneous conductivity fields for transient boundary conditions. They found that the importance of the horizontal flow components increases under transient flow conditions, a phenomenon which is typical for unsaturated but not for saturated flow. This leads to a larger lateral redistribution of solutes in the surface soil layer and a reduction of the spatial covariance of the vertical solute particle velocity in the direction of the mean flow. As a consequence, the vertical spreading of a solute plume or the spreading of an averaged BTC is smaller than under steady-state flow conditions. Similar numerical simulations considering water uptake by plants (Russo et al., 1998) suggested that crop water uptake intensifies lateral water redistribution in the root zone even more and therefore further reduces the vertical solute spreading. Using the first-order approximate velocity covariance for a steady-state flow rate equal to the time averaged flow rate in the transient flow scenario combined with the Lagrangian first-order transport approximation led to an overestimation of the solute spreading in the vertical flow direction for the transient flow scenario (Russo et al., 1994b). This is explained by the assumption of straight particle trajectories in the Lagrangian first-order analysis, so the effect of lateral redistribution during the redistribution phase after a rainfall or irrigation event is not considered. A second explanation is that with increasing variability of the unsaturated hydraulic conductivity (i.e., for drier conditions), the first-order approximations underestimate the variability of the flow more for the transverse flow components than for the longitudinal flow components. A third explanation is that close to the soil surface the effective flow rate during periods when solutes are actually leached is considerably larger than the time averaged flow rate (Vanderborght et al., 2000b) that corresponds to the steady-state flux deeper in the soil profile for which the first-order approximate velocity covariance is calculated. For noncorrelated ln_K and ln_h parameter fields that were considered by Russo et al. (1994a, 1998) the equivalent dispersivity decreases monotonically with increasing average water content or flow rate (Russo, 1995b). As a consequence, the smaller vertical spreading close to the surface boundary during transient leaching, as compared with the spreading for steady-state time-averaged flow, can be explained also by the higher water contents and lower relative flow variability during periods when leaching effectively takes place. However, following the same line of reasoning, it is anticipated that transient flow in media that show increased dispersivity with increasing degree of saturation or flow rate (i.e., media with a negative correlation between ln_K and ln_h, Miller–Miller media) above a critical saturation threshold, may actually lead to a larger relative flow variability and dispersion close to the soil surface than under steady-state conditions. Therefore, the conclusion that transient flow conditions and plant water uptake lead to a lower vertical solute spreading or equivalent dispersivity than in steady-state flow conditions cannot be generalized, since experimental data do no confirm this conclusion, or even contradict it (Vanderborght et al., 2000a).

Transient flow conditions also affect leaching rate. Figure 8 shows the simulated average solute flux at two depths in a heterogeneous soil profile for steady-state and for transient flow resulting from climatic boundary conditions plotted vs. the cumulative amount of leachate at the corresponding depths. Of note is the retardation of the solute breakthrough for the transient flow case, which is larger close to the surface boundary. Apparently more water is required to leach the tracer for transient than for steady-state conditions. The larger lateral redistribution of solute in the surface soil layer due to transient flow boundary conditions leads to a homogenization of the solute mass in the surface layer and a redistribution of solute mass from high flow rate regions to those with a lower flow rate. The effect of this redistribution or homogenization on the average leaching rate is analogous to the effect of a uniform initial concentration distribution in a heterogeneous flow field on the leaching rate. For a uniform initial concentration distribution in a heterogeneous flow field, the average solute arrival is inversely proportional to the harmonic average of the local pore water velocities close to the initial concentration profile and to the arithmetic average at larger distances from the injection (Vanderborght et al., 1998; Demmy et al., 1999). Since the harmonic average of local pore water velocities is smaller than the arithmetic one, the effective average leaching rate is smaller close to the initial uniform concentration profile than farther from it, where it converges to the time averaged flow rate divided by the time averaged volumetric water content. For steady-state flow rates, this lateral redistribution and homogenization does not take place in the surface layer and the amount of solute that leaches at a certain location is proportional to the local flow rate. This corresponds with a uniform boundary value problem, where the locally injected solute mass is proportional to the local flow rate. For a uniform boundary value problem, the leaching rate does not change with distance from the injection surface and is equal to the arithmetic average of the local pore water velocities, or the average flow rate divided by the average volumetric water content.

View larger version (12K): [in this window] [in a new window]   Fig. 8. Simulated breakthrough curves of spatially averaged solute fluxes at two depths in a realization of a two-dimensional heterogeneous soil profile. Averaged fluxes are plotted vs. the cumulative amount of leachate. Solid circles represent simulated breakthrough in a steady-state flow field; open circles are the results for a transient flow field simulated for climatic boundary conditions.

	SUMMARY

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION PREDICTION OF FIELD-SCALE... PREDICTION OF LOCAL-SCALE... TRANSPORT IN NONSTATIONARY... SUMMARY APPENDIX REFERENCES

 
Starting with the general limitation of first-order solutions of stochastic continuum equations, it is trivial that solutions that neglect higher order perturbations only apply for relatively small perturbations of the flow and transport equations. When higher order perturbations are neglected, the effect of spatially variable hydraulic properties on higher order moments (skewness and kurtosis) quantifying "tailing" of field-scale averaged BTCs or concentration profiles cannot be predicted. Tailing caused by local transport processes (e.g., nonlinear and non-equilibrium sorption) can be included in a STM and coupled with a stochastic advection velocity. Neither can higher order moments or deviations of second-order stationarity of hydraulic parameter fields (e.g., statistics that quantify connectivity of hydraulic properties) be considered and accounted for. It is important to note the difficulties already encountered in attempts to derive second-order moments of hydraulic parameters, especially the spatial (cross)covariance in the vertical direction (i.e., the direction of the mean flow) based on a discrete or grid-based sampling of the vadose zone. To determine higher order moments of the parameter fields, especially connectivities, with sufficient precision, a continuous three-dimensional mapping of the parameter fields is required. The application of noninvasive tomographic geophysical methods to characterize structures that are relevant for vadose zone flow and transport is very promising.

A second limitation imposed to obtain closed-form solutions of the moment equations is that of uniform gravity-driven steady-state flow. Despite these limitations, the direct link between larger scale transport, the structure of smaller scale soil hydraulic properties, and the state variables of the soil that is obtained from these closed-form solutions is of great importance to understanding and predicting field-scale transport. It explains and predicts differences between field-scale transport and transport in small soil columns and differences between transport in saturated media and in the vadose zone where the structure and variability of the flow field are functions of the degree of saturation.

Field-scale averaged transport of a widespread applied solute is often described using a one-dimensional equivalent CDE model. First-order approximate solutions of the stochastic flow and transport equation predict the variance of particle travel times and particle locations from which the equivalent dispersivity length, _eq, can be derived. A first important result is that first-order solutions predict an increase of _eq, with travel distance until an asymptotic value is reached, which is in accordance with results from several field-scale leaching experiments. The asymptotic state _eq may not be reached in the upper part of the soil profile (first 1 m), where a strong vertical gradient in chemical and biological soil properties also exists. Whether this vertical gradient can be combined with an equivalent CDE requires further investigation.

A second important result is that unlike in saturated media, in the vadose zone _eq depends on the flow rate or degree of saturation. The relation between _eq and flow rate is determined by the correlation between the hydraulic conductivity and the capillary length (i.e., the correlation between ln_K and ln_h). For positively correlated ln_K and ln_h, _eq tends to increase monotonically with decreasing degree of saturation whereas for a negative correlation, the relation between _eq and saturation is nonmonotonic and _eq reaches a minimum at a critical degree of saturation. It should be noted that at this moment, no experimental studies exist that show an increase of _eq with decreasing degree of saturation in real soils. This suggests that ln_K and ln_h are mostly negatively correlated and that leaching experiments have been performed using relatively large flow rates.

Since the use of an equivalent CDE implicitly implies an averaging of local concentrations, this model cannot be applied directly to upscale nonlinear local transport processes, such as nonlinear sorption and nonlinear transformations. For such problems, the use of a STM that reproduces the statistics of local concentration distributions is more adequate. The STM parameters (i.e., the variance of the stream tube velocities and the stream tube dispersivity) are related to the two-particle travel time covariance, (x,y), which can be predicted based on the structure of the hydraulic parameter fields and the local-scale dispersion by first-order solutions of the stochastic flow and transport equations. Similar to the equivalent dispersion model, the first-order solutions predict that the STM is an equivalent model with parameters that scale with travel distance and flow rate. The failure of the equivalent CDE to predict field-scale nonlinear transport and the dependence of the equivalent CDE and STM parameters on transport distance and average flow rate demonstrate the important influence of smaller scale structures and three-dimensional flow and transport processes on larger scale transport. Without information about these structural properties and without its proper implementation in larger scale equivalent models, experimental data on larger scale transport processes cannot be extrapolated to other transport distances, boundary conditions, and other substances.

From the two-particle travel time covariance, (x,y), the spatial covariance of local concentrations and stream tube velocities can be derived. The spatial covariance of stream tube velocities is less influenced by local-scale dispersion and travel time and extends over a larger range of separation distances than the spatial covariance of local concentrations. Therefore, the former reflects the structure of the flow field more clearly and may be used to inversely infer the structure of the hydraulic conductivities. This implies that a noninvasive monitoring of the three-dimensional transport process will provide a better insight (Kemna et al., 2002; Vanderborght et al., 2005) than spatially highly resolved snapshots of concentration distributions.

For transient and nonuniform flow, the moment equations may be solved numerically to obtain the nonstationary velocity covariance function. But it remains a challenge for three-dimensional non-stationary flow and transport and may be even more computationally demanding than Monte-Carlo simulations. The development of computationally efficient methods exploiting parallelisms in the moment equations is a field of further research (Zhang, 1999). However, the first-order approximation of the particle trajectory along which the velocity perturbations are evaluated becomes more problematic for transient flow conditions. The larger lateral components of the flow velocity during the redistribution stages led to a larger lateral mixing and a corresponding decrease in vertical spreading. Another aspect, which received much less attention, is the lateral solute redistribution close to the soil surface into zones that are bypassed during downward flow events. As a consequence, the amount of water required to leach a solute plume below a certain depth may be considerably larger under transient than under steady-state flow conditions. However, these conclusions are based on Monte Carlo simulations in heterogeneous fields without a correlation between ln_K and ln_h.. Whether these results can be extrapolated to fields with a negative ln_K and ln_h correlation requires further investigation. The effect of transient flow conditions on the two-particle travel time covariance, (x,y), and consequently on stream tube dispersivities and velocities has so far not been investigated. Upscaling of nonlinear transport processes in nonstationary velocity fields using a STM therefore remains a challenge.

Finally, it should be noted that we only considered the spatial variability of soil hydraulic properties. Spatial variability of chemical and biological processes has an equally important impact on transport and transformation of substances in the subsurface. However, experimental data are scarce or nonexistent (Albrecht et al., 2002), especially for the spatial variability of biological processes (e.g., root water and nutrient uptake) and their correlation with hydraulic and chemical properties.

	APPENDIX

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION PREDICTION OF FIELD-SCALE... PREDICTION OF LOCAL-SCALE... TRANSPORT IN NONSTATIONARY... SUMMARY APPENDIX REFERENCES

 
Relation between ²(x₁), (x,x), and ²(x)|_v(x)
The unconditional arrival time variance, that is, the variance of arrival times in all realizations of the velocity field, v(x), at a certain point x of particles that are injected in an injection plane perpendicular to the average flow direction and at distance x₁ from the observation point is defined in terms of the unconditional travel time moments as

[A1]

where represents the expected value of a random variable. In a second-order stationary velocity field v(x), the unconditional arrival time variance ²(x) is independent of the location of x on a reference plane perpendicular to the average flow direction, so the unconditional variance of arrival times at a certain point x is identical to the arrival times on the reference surface: ²(x) = ²(x₁). For ergodic velocity fields and large injection and reference surfaces, the unconditional arrival time variance at a reference plane is identical to the arrival time variance at a reference plane conditioned on a certain realization of the velocity field: ²(x₁) = ²(x₁)|_v(x). The variance of arrival times at point x conditioned on a realization of the velocity field (i.e., in a specific realization of the velocity field, ²(x)|_v(x)) is defined in terms of the conditional travel time moments as

[A2]

Since the expected value of the conditional expectation of a random variable is equal to the expected value of the unconditioned random variable, the unconditional travel time variance Eq. [A1] can be written after adding and subtracting (x)²|_v(x) as:

[A3]

The first and last terms of the right-hand side of Eq. [A3] represent the expected value of the variance of arrival times at point x conditioned on the velocity field: ²(x)|_v(x). The second and third terms represent the variance of the expected arrival times at point x conditioned on the velocity field: ²|_v(x)(x). In other words, the unconditional arrival time variance equals the sum of the expected conditioned arrival time variance and the variance about the conditioned mean arrival time.

_|v(x)²(x) can be rewritten as

[A4]

where '⁽¹⁾|_v(x) and '⁽²⁾|_v(x) are independent perturbations of the arrival times ⁽¹⁾|_v(x), ⁽²⁾|_v(x) of two different particles around the expected arrival time in a given realization of the velocity field. Since (x,x) equals ²|_v(x)(x), Eq. [A3] can be written as

[A5]

	ACKNOWLEDGMENTS

 
We thank the reviewers and the associate editor, Gerrit de Rooij, for their helpful comments.

	REFERENCES

TOP EXECUTIVE SUMMARY ABSTRACT INTRODUCTION PREDICTION OF FIELD-SCALE... PREDICTION OF LOCAL-SCALE... TRANSPORT IN NONSTATIONARY... SUMMARY APPENDIX REFERENCES

 

• Abdou, H.M., and M. Flury. 2004. Simulation of water flow and solute transport in free-drainage lysimeters and field soils with heterogeneous structures. Eur. J. Soil Sci. 55:229–241.[CrossRef]
• Albrecht, A., U. Schultze, M. Liedgens, H. Fluhler, and E. Frossard. 2002. Incorporating soil structure and root distribution into plant uptake models for radionuclides: Toward a more physically based transfer model. J. Environ. Radioact. 59:329–350.[Medline]
• Amoozegarfard, A., D.R. Nielsen, and A.W. Warrick. 1982. Soil solute concentration distributions for spatially varying pore water velocities and apparent diffusion-coefficients. Soil Sci. Soc. Am. J. 46:3–9.
• Andricevic, R. 1998. Effects of local dispersion and sampling volume on the evolution of concentration fluctuations in aquifers. Water Resour. Res. 34:1115–1129.
• Attinger, S., J.D. Micha, and W. Kinzelbach. 2003. Multiscale modeling of nonlinearly adsorbing solute transport. Multiscale Model. Simul. 1:408–431.[CrossRef]
• Bellin, A., and Y. Rubin. 2004. On the use of peak concentration arrival times for the inference of hydrogeological parameters. Water Resour. Res. 40:W07401. doi:10.1029/2003WR002179.[CrossRef]
• Beven, K.J., D.E. Henderson, and A.D. Reeves. 1993. Dispersion parameters for undisturbed partially saturated soil. J. Hydrol. (Amsterdam) 143:19–43.
• Biggar, J.W., and D.R. Nielsen. 1976. Spatial variability of leaching characteristics of a field soil. Water Resour. Res. 12:78–84.
• Birkholzer, J., and C.F. Tsang. 1997. Solute channeling in unsaturated heterogeneous porous media. Water Resour. Res. 33:2221–2238.[CrossRef]
• Boesten, J. 2004. Influence of dispersion length on leaching calculated with pearl, pelmo and przm for focus groundwater scenarios. Pest Manage. Sci. 60:971–980.[CrossRef]
• Bresler, E., and G. Dagan. 1981. Convective and pore scale dispersive solute transport in unsaturated heterogeneous fields. Water Resour. Res. 17:1683–1693.
• Cirpka, O.A. 2002. Choice of dispersion coefficients in reactive transport calculations on smoothed fields. J. Contam. Hydrol. 58:261–282.[Medline]
• Cirpka, O.A., and P.K. Kitanidis. 2000. An advective-dispersive stream tube approach for the transfer of conservative-tracer data to reactive transport. Water Resour. Res. 36:1209–1220.[CrossRef]
• Cirpka, O.A., and P.K. Kitanidis. 2002. Numerical evaluation of solute dispersion and dilution in unsaturated heterogeneous media. Water Resour. Res. 38:1220. doi:10.1029/2001WR001262.[CrossRef]
• Dagan, G. 1989. Flow and transport in porous formations. Springer-Verlag, Berlin.
• Dagan, G., and E. Bresler. 1979. Solute-dispersion in unsaturated heterogeneous soil at field scale. 1. Theory. Soil Sci. Soc. Am. J. 43:461–467.[Abstract/Free Full Text]
• Dagan, G., and V. Cvetkovic. 1996. Reactive transport and immiscible flow in geological media. 1. General theory. Proc. R. Soc. London. Ser. A. 452:285–301.[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.[CrossRef]
• Dagan, G., and A. Fiori. 1997. The influence of pore-scale dispersion on concentration statistical moments in transport through heterogeneous aquifers. Water Resour. Res. 33:1595–1605.[CrossRef]
• de Rooij, G.H., R.T.A. Kasteel, A. Papritz, and H. Fluhler. 2004. Joint distributions of the unsaturated soil hydraulic parameters and their effect on other variates. Available at www.vadosezonejournal.org. Vadose Zone J. 3:947–955.[Abstract/Free Full Text]
• Demmy, G., S. Berglund, and W. Graham. 1999. Injection mode implications for solute transport in porous media: Analysis in a stochastic lagrangian framework. Water Resour. Res. 35:1965–1973.[CrossRef]
• Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach. 2000. Temporal behavior of a solute cloud in a heterogeneous porous medium. 1. Point-like injection. Water Resour. Res. 36:3591–3604.[CrossRef]
• Destouni, G. 1992. The effect of vertical soil heterogeneity on field scale solute flux. Water Resour. Res. 28:1303–1309.[CrossRef]
• Deurer, M., W.H.M. Duijnisveld, and J. Bottcher. 2000. Spatial analysis of water characteristic functions in a sandy podzol under pine forest. Water Resour. Res. 36:2925–2935.[CrossRef]
• Deurer, M., W.H.M. Duijnisveld, J. Bottcher, and G. Klump. 2001. Heterogeneous solute flow in a sandy soil under a pine forest: Evaluation of a modeling concept. J. Plant Nutr. Soil Sci. 164:601–610.[CrossRef]
• Feyen, J., D. Jacques, A. Timmerman, and J. Vanderborght. 1998. Modelling water flow and solute transport in heterogeneous soils: A review of recent approaches. J. Agric. Eng. Res. 70:231–256.[CrossRef]
• Finkel, M., R. Liedl, and G. Teutsch. 1999. Modelling surfactant-enhanced remediation of polycyclic aromatic hydrocarbons. Environ. Model. Softw. 14:203–211.
• Fiori, A. 1996. Finite Peclet extensions of Dagan's solutions to transport in anisotropic heterogeneous formations. Water Resour. Res. 32:193–198.
• Fiori, A., S. Berglund, V. Cvetkovic, and G. Dagan. 2002. A first-order analysis of solute flux statistics in aquifers: The combined effect of pore-scale dispersion, sampling, and linear sorption kinetics. Water Resour. Res. 38(8):1137. doi:10.1029/2001WR000678.[CrossRef]
• Fiori, A., and G. Dagan. 1999. Concentration fluctuations in transport by groundwater: Comparison between theory and field experiments. Water Resour. Res. 35:105–112.[CrossRef]
• Fiori, A., and G. Dagan. 2000. Concentration fluctuations in aquifer transport: A rigorous first-order solution and applications. J. Contam. Hydrol. 45:139–163.[CrossRef]
• Flühler, H., W. Durner, and M. Flury. 1996. Lateral solute mixing processes- a key for understanding field-scale transport of water and solutes. Geoderma 70:165–183.[CrossRef][Web of Science]
• FOCUS. 2000. FOCUS groundwater scenarios in the EU review of active substances. EC Document Reference SANCO/321/2000 rev2.
• Forrer, I., R. Kasteel, M. Flury, and H. Flühler. 1999. Longitudinal and lateral dispersion in an unsaturated field soil. Water Resour. Res. 35:3049–3060.[CrossRef]
• Forrer, I., A. Papritz, R. Kasteel, H. Flühler, and D. Luca. 2000. Quantifying dye tracers in soil profiles by image processing. Eur. J. Soil Sci. 51:313–322.
• Foussereau, X., W.D. Graham, G.A. Akpoji, G. Destouni, and P.S.C. Rao. 2000b. Stochastic analysis of transport in unsaturated heterogeneous soils under transient flow regimes. Water Resour. Res. 36:911–921.[CrossRef]
• Foussereau, X., W.D. Graham, and P.S.C. Rao. 2000a. Stochastic analysis of transient flow in unsaturated heterogeneous soils. Water Resour. Res. 36:891–910.[CrossRef]
• Gelhar, L.W. 1993. Stochastic subsurface hydrology. Prentice Hall Inc., Englewood Cliffs, NJ.
• Ginn, T.R. 2001. Stochastic-convective transport with nonlinear reactions and mixing: Finite streamtube ensemble formulation for multicomponent reaction systems with intra-streamtube dispersion. J. Contam. Hydrol. 47:1–28.[Medline]
• Gomez-Hernandez, J.J., and X.-H. Wen. 1998. To be or not to be multi-gaussian? A reflection on stochastic hydrogeology. Adv. Water Resour. 21:47–61.
• Grayson, R., and G. Blöschl. 2000. Spatial patterns in catchment hydrology: Observations and modelling Cambridge Univ. Press, U.K.
• Gutjahr, A.L., S. Hatch, and B. Bullard. 1994. Random field generation using the fast fourier transform. New Mexico Institute of Mining and Technology, Socorro, NM.
• Hammel, K., J. Gross, G. Wessolek, and K. Roth. 1999. Two-dimensional simulation of bromide transport in a heterogeneous field soil with transient unsaturated flow. Eur. J. Soil Sci. 50:633–647.[CrossRef]
• Hammel, K., and K. Roth. 1998. Approximation of asymptotic dispersivity of conservative solute in unsaturated heterogeneous media with steady state flow. Water Resour. Res. 34:709–715.[CrossRef]
• Harter, T., and J.W. Hopmans. 2004. Role of vadose-zone flow processes in regional-scale hydrology: Review, opportunities, and challenges, p. 179–208. In R.A. Feddes et al (ed.) International symposium on unsaturated zone modelling: Progress, challenges and applications. 3–5 Oct. 2004. Kluwer, Wageningen, The Netherlands.
• Harter, T., A.L. Gutjahr, and T.C.J. Yeh. 1996. Linearized cosimulation of hydraulic conductivity, pressure head, and flux in saturated and unsaturated, heterogeneous porous media. J. Hydrol. (Amsterdam) 183:169–190.
• Harter, T., and T.C.J. Yeh. 1996. Stochastic analysis of solute transport in heterogeneous, variably saturated soils. Water Resour. Res. 32:1585–1595.[CrossRef]
• Harter, T., and T.C.J. Yeh. 1998. Flow in unsaturated random porous media, nonlinear numerical analysis and comparison to analytical stochastic models. Adv. Water Resour. 22:257–272.
• Hopmans, J.W. 1987. A comparison of various methods to scale soil hydraulic-properties. J. Hydrol. (Amsterdam) 93:241–256.
• Hubbard, S., and Y. Rubin. 2004. Hydrogeophysics overview. p. 3–21. In Y. Rubin and S. Hubbard (ed.) Hydrogeophysics. Springer, New York.
• Indelman, P., and Y. Rubin. 1996. Solute transport in nonstationary velocity fields. Water Resour. Res. 32:1259–1267.[CrossRef]
• Jensen, K.H., and A. Mantoglou. 1992. Application of stochastic unsaturated flow theory, numerical simulations, and comparisons to field observations. Water Resour. Res. 28:269–284.[CrossRef]
• Jose, S.C., M.A. Rahman, and O.A. Cirpka. 2004. Large-scale sandbox experiment on longitudinal effective dispersion in heterogeneous porous media. Water Resour. Res. 40:W12415. doi:10.1029/2004WR003363.[CrossRef]
• Jury, W.A. 1982. Simulation of solute transport using a transfer-function model. Water Resour. Res. 18:363–368.
• Jury, W.A. 1985. Spatial variability of soil physical parameters in solute migration: A critical literature review. Project 2485-6. EPRI EA-4228. Electric Power Res. Inst., Palo Alto, CA.
• Jury, W.A., and J. Gruber. 1989. A stochastic-analysis of the influence of soil and climatic variability on the estimate of pesticide groundwater pollution potential. Water Resour. Res. 25:2465–2474.
• Jury, W.A., and G. Sposito. 1985. Field calibration and validation of solute transport models for the unsaturated zone. Soil Sci. Soc. Am. J. 49:1331–1341.[Abstract/Free Full Text]
• Jury, W.A., G. Sposito, and R.E. White. 1986. A transfer-function model of solute transport through soil. 1. Fundamental concepts. Water Resour. Res. 22:243–247.
• Kapoor, V., and L.W. Gelhar. 1994. Transport in three-dimensionally heterogeneous aquifers. 1. Dynamics of concentration fluctuations. Water Resour. Res. 30:1775–1788.[CrossRef]
• Kapoor, V., and P.K. Kitanidis. 1996. Concentration fluctuations and dilution in two-dimensionally periodic heterogeneous porous media. Transp. Porous Media 22:91–119.[CrossRef]
• Kapoor, V., and P.K. Kitanidis. 1998. Concentration fluctuations and dilution in aquifers. Water Resour. Res. 34:1181–1193.[CrossRef]
• Kasteel, R. 1997. Solute transport in an unsaturated field soil: Describing heterogeneous flow fields using spatial distribution of hydraulic properties. ETH, Zurich, Switzerland.
• Kasteel, R., M. Burkhardt, S. Giesa, and H. Vereecken. 2005. Characterization of field tracer transport using high-resolution images. Available at www.vadosezonejournal.org. Vadose Zone J. 4:101–111.[Abstract/Free Full Text]
• Kasteel, R., H.J. Vogel, and K. Roth. 2000. From local hydraulic properties to effective transport in soil. Eur. J. Soil Sci. 51:81–91.
• Kasteel, R., H.J. Vogel, and K. Roth. 2002. Effect of non-linear adsorption on the transport behaviour of brilliant blue in a field soil. Eur. J. Soil Sci. 53:231–240.
• Kemna, A., J. Vanderborght, B. Kulessa, and H. Vereecken. 2002. Imaging and characterisation of subsurface solute transport using electrical resistivity tomography (ert) and equivalent transport models. J. Hydrol. (Amsterdam) 267:125–146.
• Khaleel, R., T.C.J. Yeh, and Z.M. Lu. 2002. Upscaled flow and transport properties for heterogeneous unsaturated media. Water Resour. Res. 38(5):1053. doi:10.1029/2000WR000072.[CrossRef]
• Kitanidis, P.K. 1994. The concept of the dilution index. Water Resour. Res. 30:2011–2026.[CrossRef]
• Lu, Z.M., and D.X. Zhang. 2003. Solute spreading in nonstationary flows in bounded, heterogeneous unsaturated-saturated media. Water Resour. Res. 39:1049. doi:10.1029/2001WR000908.[CrossRef]
• Mallants, D., D. Jacques, P.H. Tseng, M.Th. van Genuchten, and J. Feyen. 1997. Comparison of three hydraulic property measurement methods. J. Hydrol. (Amsterdam) 199:295–318.
• Mallants, D., B.P. Mohanty, D. Jacques, and J. Feyen. 1996. Spatial variability of hydraulic properties in a multi-layered soil profile. Soil Sci. 161:167–181.[CrossRef]
• Malmstrom, M.E., G. Destouni, and P. Martinet. 2004. Modeling expected solute concentration in randomly heterogeneous flow systems with multicomponent reactions. Environ. Sci. Technol. 38:2673–2679.[Medline]
• Mantoglou, A. 1992. A theoretical approach for modeling unsaturated flow in spatially-variable soils—Effective flow models in finite domains and nonstationarity. Water Resour. Res. 28:251–267.[CrossRef]
• Mantoglou, A., and L.W. Gelhar. 1987a. Stochastic modeling of large-scale transient unsaturated flow systems. Water Resour. Res. 23:37–46.
• Mantoglou, A., and L.W. Gelhar. 1987b. Capillary tension head variance, mean soil-moisture content, and effective specific soil-moisture capacity of transient unsaturated flow in stratified soils. Water Resour. Res. 23:47–56.
• Mantoglou, A., and L.W. Gelhar. 1987c. Effective hydraulic conductivities of transient unsaturated flow in stratified soils. Water Resour. Res. 23:57–67.
• Mohanty, B.P., M.D. Ankeny, R. Horton, and R.S. Kanwar. 1994. Spatial-analysis of hydraulic conductivity measured using disc infiltrometers. Water Resour. Res. 30:2489–2498.[CrossRef]
• Neuendorf, O. 1997. Numerische three-dimensional simulation des stofftransports in einem heterogenen aquifer Berichte des Forschungszentrums Jülich 3421. Forschungszentrum Jülich GmbH, Jülich, Germany.
• Neuweiler, I., and O.A. Cirpka. 2005. Homogenization of Richards equation in permeability fields with different connectivities. Water Resour. Res. 41:W02009. doi:10.1029/2004WR003329.[CrossRef]
• Pannone, M., and P.K. Kitanidis. 1999. Large-time behavior of concentration variance and dilution in heterogeneous formations. Water Resour. Res. 35:623–634.[CrossRef]
• Pannone, M., and P.K. Kitanidis. 2001. Large-time spatial covariance of concentration of conservative solute and application to the cape cod tracer test. Transp. Porous Media 42:109–132.[CrossRef]
• Polmann, D.J., D. McLaughlin, S. Luis, L.W. Gelhar, and R. Ababou. 1991. Stochastic modeling of large-scale flow in heterogeneous unsaturated soils. Water Resour. Res. 27:1447–1458.[CrossRef]
• Raats, P.A.C. 1976. Analytical solutions of a simplified flow equation. Trans. ASAE 19:683–689.
• Rockhold, M.L., R.E. Rossi, and R.G. Hills. 1996. Application of similar media scaling and conditional simulation for modeling water flow and tritium transport at the Las Cruces trench site. Water Resour. Res. 32:595–609.[CrossRef]
• Roth, K. 1995. Steady-state flow in an unsaturated, two-dimensional, macroscopically homogeneous, Miller-similar medium. Water Resour. Res. 31:2127–2140.[CrossRef]
• Roth, K., and K. Hammel. 1996. Transport of conservative chemical through an unsaturated two-dimensional Miller-similar medium with steady state flow. Water Resour. Res. 32:1653–1663.[CrossRef]
• Rubin, Y. 2003. Applied stochastic hydrogeology. Oxford Univ. Press, New York.
• Rubin, Y., and S. Ezzedine. 1997. The travel times of solutes at the Cape Cod tracer experiment: Data analysis, modeling, and structural parameters inference. Water Resour. Res. 33:1537–1547.[CrossRef]
• Russo, D. 1995a. On the velocity covariance and transport modeling in heterogeneous anisotropic porous formations. 1. Saturated flow. Water Resour. Res. 31:129–137.[CrossRef]
• Russo, D. 1995b. On the velocity covariance and transport modeling in heterogeneous anisotropic porous formations. 2. Unsaturated flow. Water Resour. Res. 31:139–145.[CrossRef]
• Russo, D. 1997. On the estimation of parameters of log-unsaturated conductivity covariance from solute transport data. Adv. Water Resour. 20:191–205.[CrossRef]
• Russo, D., and M. Bouton. 1992. Statistical-analysis of spatial variability in unsaturated flow parameters. Water Resour. Res. 28:1911–1925.[CrossRef]
• Russo, D., I. Russo, and A. Laufer. 1997. On the spatial variability of parameters of the unsaturated hydraulic conductivity. Water Resour. Res. 33:947–956.[CrossRef]
• Russo, D., J. Zaidel, and A. Laufer. 1994a. Stochastic-analysis of solute transport in partially saturated heterogeneous soil. 1. Numerical experiments. Water Resour. Res. 30:769–779.[CrossRef]
• Russo, D., J. Zaidel, and A. Laufer. 1994b. Stochastic-analysis of solute transport in partially saturated heterogeneous soil. 2. Prediction of solute spreading and breakthrough. Water Resour. Res. 30:781–790.[CrossRef]
• Russo, D., J. Zaidel, and A. Laufer. 1998. Numerical analysis of flow and transport in a three-dimensional partially saturated heterogeneous soil. Water Resour. Res. 34:1451–1468.[CrossRef]
• Sun, A.Y., and D.X. Zhang. 2000. Prediction of solute spreading during vertical infiltration in unsaturated, bounded heterogeneous porous media. Water Resour. Res. 36:715–723.[CrossRef]
• Tiktak, A., D.S. de Nie, J.D. Pineros Garcet, A. Jones, and M. Vanclooster. 2004. Assessment of the pesticide leaching risk at the Pan-European level. The Europearl approach. J. Hydrol. (Amsterdam) 289:222–238.
• Toride, N., and F.J. Leij. 1996a. Convective dispersive stream tube model for field-scale solute transport. 1. Moment analysis. Soil Sci. Soc. Am. J. 60:342–352.[Abstract/Free Full Text]
• Toride, N., and F.J. Leij. 1996b. Convective dispersive stream tube model for field-scale solute transport. 2. Examples and calibration. Soil Sci. Soc. Am. J. 60:352–361.[Abstract/Free Full Text]
• Ünlü, K., D.R. Nielsen, J.W. Biggar, and F. Morkoc. 1990. Statistical parameters characterizing the spatial variability of selected soil hydraulic properties. Soil Sci. Soc. Am. J. 54:1537–1547.[Abstract/Free Full Text]
• Ursino, N., T. Gimmi, and H. Flühler. 2001a. Dilution of non-reactive tracers in variably saturated sandy structures. Adv. Water Resour. 24:877–885.[CrossRef]
• Ursino, N., T. Gimmi, and H. Flühler. 2001b. Combined effects of heterogeneity, anisotropy, and saturation on steady state flow and transport: A laboratory sand tank experiment. Water Resour. Res. 37:201–208.
• van Dam, J.C., G.H. de Rooij, M. Heinen, and F. Stagnitti. 2004. Concepts and dimensionality in modeling unsaturated water flow and solute transport. p. 1–36. In R.A. Feddes et al. (ed.) International symposium on unsaturated zone modelling: Progress, challenges and applications. 3–5 Oct. 2004. Kluwer, Wageningen, The Netherlands.
• Vanderborght, J. 2001. Concentration variance and spatial covariance in second-order stationary heterogeneous conductivity fields. Water Resour. Res. 37:1893–1912.[CrossRef]
• Vanderborght, J., D. Jacques, and J. Feyen. 2000b. Deriving transport parameters from transient flow leaching experiments by approximate steady-state flow convection-dispersion models. Soil Sci. Soc. Am. J. 64:1317–1327.[Abstract/Free Full Text]
• Vanderborght, J., D. Jacques, D. Mallants, P.H. Tseng, and J. Feyen. 1997. Comparison between field measurements and numerical simulation of steady-state solute transport in a heterogeneous soil profile. Hydrol. Earth Syst. Sci. 4:853–871.
• Vanderborght, J., A. Kemna, H. Hardelauf, and H. Vereecken. 2005. Potential of electrical resistivity tomography to infer aquifer transport characteristics from tracer studies. A synthetic case study. Water Resour. Res. 41:W06013. doi:10.1029/2004WR003774.[CrossRef]
• Vanderborght, J., D. Mallants, and J. Feyen. 1998. Solute transport in a heterogeneous soil for boundary and initial conditions: Evaluation of first-order approximations. Water Resour. Res. 34:3255–3270.[CrossRef]
• Vanderborght, J., A. Timmerman, and J. Feyen. 2000a. Solute transport for steady-state and transient flow in soils with and without macropores. Soil Sci. Soc. Am. J. 64:1305–1317.[Abstract/Free Full Text]
• Vanderborght, J., M. Vanclooster, A. Timmerman, P. Seuntjens, D. Mallants, D.J. Kim, D. Jacques, L. Hubrechts, C. Gonzalez, J. Feyen, J. Diels, and J. Deckers. 2001. Overview of inert tracer experiments in key Belgian soil types: Relation between transport and soil morphological and hydraulic properties. Water Resour. Res. 37:2873–2888.[CrossRef]
• Vanderborght, J., and H. Vereecken. 2001. Analyses of locally measured bromide breakthrough curves from a natural gradient tracer experiment at Krauthausen. J. Contam. Hydrol. 48:23–43.[Medline]
• Vanderborght, J., and H. Vereecken. 2002. Estimation of local scale dispersion from local breakthrough curves during a tracer test in a heterogeneous aquifer: The Lagrangian approach. J. Contam. Hydrol. 54:141–171.[Medline]
• Vereecken, H., S. Hubbard, A. Binley, and T. Ferre. 2004. Hydrogeophysics: An introduction from the guest editors. Available at www.vadosezonejournal.org. Vadose Zone J. 3:1060–1062.[Free Full Text]
• Vereecken, H., R. Kaiser, M. Dust, and T. Putz. 1997. Evaluations of the multistep outflow method for the determination of unsaturated hydraulic properties of soils. Soil Sci. 162:618–631.[CrossRef]
• Vereecken, H., G. Lindenmayer, O. Neuendorf, U. Döring, and R. Seidemann. 1994. Trace, a mathematical model for reactive transport in 3D variable saturated porous media. KFA, ICG-4 Internal Rep. 501494. Forschungszentrum Jülich GmbH, Jülich, Germany.
• Vogel, H.J., and K. Roth. 2003. Moving through scales of flow and transport in soil. J. Hydrol. (Amsterdam) 272:95–106.
• Vogel, T., M. Cislerova, and J.W. Hopmans. 1991. Porous media with linearly variable hydraulic properties. Water Resour. Res. 27:2735–2741.[CrossRef]
• Warrick, A.W., G.J. Mullen, and D.R. Nielsen. 1977. Scaling field-measured soil hydraulic properties using a similar media concept. Water Resour. Res. 13:355–362.[CrossRef]
• Wen, X.H., and J. Gomez-Hernandez. 1998. Numerical modeling of macrodispersion in heterogeneous media: A comparison of multi-gaussian and non-multi-gaussian models. J. Contam. Hydrol. 30:129–156.[CrossRef]
• White, I., and M.J. Sully. 1992. On the variability and use of the hydraulic conductivity alpha-parameter in stochastic treatments of unsaturated flow. Water Resour. Res. 28:209–213.[CrossRef]
• Wildenschild, D., and K.H. Jensen. 1999. Laboratory investigations of effective flow behavior in unsaturated heterogeneous sands. Water Resour. Res. 35:17–27.[CrossRef]
• Yeh, T.C.J., L.W. Gelhar, and A.L. Gutjahr. 1985a. Stochastic analysis of unsaturated flow in heterogeneous soils. 1. Statistically isotropic media. Water Resour. Res. 21:447–456.
• Yeh, T.C.J., L.W. Gelhar, and A.L. Gutjahr. 1985b. Stochastic-analysis of unsaturated flow in heterogeneous soils. 2. Statistically anisotropic media with variable-alpha. Water Resour. Res. 21:457–464.[CrossRef]
• Yeh, T.C.J., L.W. Gelhar, and A.L. Gutjahr. 1985c. Stochastic analysis of unsaturated flow in heterogeneous soils. 3. Observations and applications. Water Resour. Res. 21:465–471.
• Zhang, D.X. 1999. Nonstationary stochastic analysis of transient unsaturated flow in randomly heterogeneous media. Water Resour. Res. 35:1127–1141.[CrossRef]
• Zhang, D.X. 2002. Stochastic methods for flow in porous media. Academic Press, London.
• Zhang, D.X., and S.P. Neuman. 1996. Effect of local dispersion on solute transport in randomly heterogeneous media. Water Resour. Res. 32:2715–2723.[CrossRef]
• Zinn, B., and C.F. Harvey. 2003. When good statistical models of aquifer heterogeneity go bad: A comparison of flow, dispersion, and mass transfer in connected and multivariate gaussian hydraulic conductivity fields. Water Resour. Res. 39:1051. doi:10.1029/2001WR001146.[CrossRef]

This article has been cited by other articles:

				J. Koestel, R. Kasteel, A. Kemna, O. Esser, M. Javaux, A. Binley, and H. Vereecken Imaging Brilliant Blue Stained Soil by Means of Electrical Resistivity Tomography Vadose Zone J., November 17, 2009; 8(4): 963 - 975. [Abstract] [Full Text] [PDF]

				M. F. Dyck and R. G. Kachanoski Measurement of Transient Soil Water Flux Across a Soil Horizon Interface Soil Sci. Soc. Am. J., August 19, 2009; 73(5): 1604 - 1613. [Abstract] [Full Text] [PDF]

				H. Vereecken, P. Burauel, J. Groeneweg, E. Klumpp, W. Mittelstaedt, H.-D. Narres, T. Putz, J. van der Kruk, J. Vanderborght, and F. Wendland Research at the Agrosphere Institute: From the Process Scale to the Catchment Scale Vadose Zone J., August 11, 2009; 8(3): 664 - 669. [Abstract] [Full Text] [PDF]

				J. Koestel, J. Vanderborght, M. Javaux, A. Kemna, A. Binley, and H. Vereecken Noninvasive 3-D Transport Characterization in a Sandy Soil Using ERT: 1. Investigating the Validity of ERT-derived Transport Parameters Vadose Zone J., August 11, 2009; 8(3): 711 - 722. [Abstract] [Full Text] [PDF]

				M. D. Madsen, D. G. Chandler, and W. D. Reynolds Accounting for Bias and Boundary Condition Effects on Measurements of Saturated Core Hydraulic Conductivity Soil Sci. Soc. Am. J., May 1, 2008; 72(3): 750 - 757. [Abstract] [Full Text] [PDF]

				J. Vanderborght and H. Vereecken Review of Dispersivities for Transport Modeling in Soils Vadose Zone J., January 24, 2007; 6(1): 29 - 52. [Abstract] [Full Text] [PDF]

				J. Vanderborght and H. Vereecken One-Dimensional Modeling of Transport in Soils with Depth-Dependent Dispersion, Sorption and Decay Vadose Zone J., January 24, 2007; 6(1): 140 - 148. [Abstract] [Full Text] [PDF]

				D. L. Corwin, J. Hopmans, and G. H. de Rooij From Field- to Landscape-Scale Vadose Zone Processes: Scale Issues, Modeling, and Monitoring Vadose Zone J., March 8, 2006; 5(1): 129 - 139. [Abstract] [Full Text] [PDF]

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Similar articles in this journal Similar articles in Web of Science Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via Web of Science (4) Citing Articles via Google Scholar Google Scholar Articles by Vanderborght, J. Articles by Vereecken, H. Search for Related Content PubMed Articles by Vanderborght, J. Articles by Vereecken, H. GeoRef GeoRef Citation Agricola Articles by Vanderborght, J. Articles by Vereecken, H. Related Collections Solute Transport Models Stochastic Processes Vadose Zone Processes and Chemical Transport