Abstract
Darcy velocities for lognormal hydraulic conductivity with small variance and finite correlation length were approximated by periodic random fields. Accurate simulations of two‐dimensional advection‐dispersion processes were achieved with the global random walk algorithm, using 1010 particles in every transport realization. Reliable statistical estimations were obtained by averaging over 256 realizations.
The main result is a numerical evidence for the mean square convergence of the actual concentrations to the macrodispersion process predicted by a known limit theorem. For small initial plumes the ergodic behavior can be expected after thousands of advection timescales, when the deviation from the theoretical prediction of the cross‐section space‐averaged concentration monotonously decays and falls under 20%. The increase of the transverse dimension of the initial plume slows down the approach to the quasi‐ergodic state and has a nonlinear effect on the variability of the actual concentrations and dispersivities.
Authors
N. Suciu
Institute of Applied Mathematics, Friedrich-Alexander University of
Erlangen-Nuremberg, Erlangen, Germany
C. Vamos
Tiberiu Popoviciu Institute of Numerical Analysis, Romanian Academy, Cluj-Napoca, Romania
J. Vanderborght
Institute of Agrosphere (ICG-IV), Research Center Julich, Julich, Germany.
H. Hardelauf
Institute of Agrosphere (ICG-IV), Research Center Julich, Julich, Germany.
H. Vereecken
Institute of Agrosphere (ICG-IV), Research Center Julich, Julich, Germany.
Keywords
?
Paper coordinates
N. Suciu, C. Vamoş, J. Vanderborght, H. Hardelauf, H. Vereecken (2006), Numerical investigations on ergodicity of solute transport in heterogeneous aquifers, Water Resour. Res., 42, W04409, doi: 10.1029/2005WR004546
About this paper
Print ISSN
Not available yet.
Online ISSN
Not available yet.
Google Scholar Profile
Numerical investigations on ergodicity of solute
transport in heterogeneous aquifers
N. Suciu,
1
C. Vamos¸,
2
J. Vanderborght,
3
H. Hardelauf,
3
and H. Vereecken
3
Received 31 August 2005; revised 10 December 2005; accepted 17 January 2006; published 18 April 2006.
[1] Darcy velocities for lognormal hydraulic conductivity with small variance and finite
correlation length were approximated by periodic random fields. Accurate simulations of
two-dimensional advection-dispersion processes were achieved with the global random
walk algorithm, using 10
10
particles in every transport realization. Reliable statistical
estimations were obtained by averaging over 256 realizations. The main result is a
numerical evidence for the mean square convergence of the actual concentrations to the
macrodispersion process predicted by a known limit theorem. For small initial plumes the
ergodic behavior can be expected after thousands of advection timescales, when the
deviation from the theoretical prediction of the cross-section space-averaged concentration
monotonously decays and falls under 20%. The increase of the transverse dimension of the
initial plume slows down the approach to the quasi-ergodic state and has a nonlinear effect
on the variability of the actual concentrations and dispersivities.
Citation: Suciu, N., C. Vamos¸, J. Vanderborght, H. Hardelauf, and H. Vereecken (2006), Numerical investigations on ergodicity of
solute transport in heterogeneous aquifers, Water Resour. Res., 42, W04409, doi:10.1029/2005WR004546.
1. Introduction
[2] It is generally admitted that groundwater quality is
mainly affected by the transport of dissolved chemicals
through soils and aquifers. The classical model is based
on a dispersion and advection mechanism which describes
the transport at some ‘‘local scale.’’ Further, one assumes
that the variability of the solute movement in subsurface
water is caused by the heterogeneity of the hydraulic
conductivity which, for a given natural formation, is effi-
ciently described as a realization of a random space
function [Sposito et al., 1986; Hassan et al., 1998]. The
corresponding advection velocity field becomes a random
function also and the transport in natural porous media is
described by a stochastic model [Dagan, 1984] which
follows the approach for turbulent diffusion in atmosphere
[Taylor, 1921], or, in terms of the mathematical theory of
stochastic processes, as a diffusion in a random velocity
field [Matheron and de Marsily, 1980; Avellaneda and
Majda, 1992].
[3] As shown by Sposito et al. [1986], once the advec-
tion-dispersion model for local scale has been inferred, the
stochastic modeling of transport in groundwater consists of
two successive steps. First, the behavior of the ensemble
averaged concentration has to be investigated to look for the
existence of an upscaled diffusive behavior called ‘‘macro-
dispersion.’’ The second issue is to assess the applicability
of the ensemble statistics to predictions made for a single
groundwater system. The latter, which is the central prob-
lem in stochastic modeling, is generally referred to as
‘‘ergodicity’’ in hydrogeological literature [Dagan, 1984,
1987; Kabala and Sposito, 1994; Sposito, 1997; Fiori,
1998; Trefry et al., 2003; Jankovic´ et al., 2003]. It is in
this sense that the term ergodicity will be used in the
following.
[4] Even though abundant literature has been produced in
the last two decades, the investigations on ergodicity were
in most cases limited to the study of the second moments of
the solute plume and the conclusions often disagree. For
instance, it is accepted that the length scale to reach the
asymptotic limit predicted by the stochastic model is im-
practical for real contamination problems [Berkowitz, 2001;
Schwarze et al., 2001; Dentz et al., 2002, 2003; Eberhard,
2004]. However, recent numerical investigations conclude
that, for extended initial plumes, the dispersivities behave
ergodically, at relatively small distances from the injection
domain [Jankovic´ et al., 2003]. These conclusions do not
agree with the results of Trefry et al. [2003], which show
that the dispersivities significantly differ from realization to
realization even after hundreds of heterogeneity scales and
the concentrations do not reach the Gaussian limit predicted
by the stochastic macrodispersion model. The result of
Trefry et al. [2003] is consistent with other numerical
simulations which show that the attainment of a quasi-
ergodic state is more complicated than indicated by some
analytical approaches [Naff et al., 1998b].
[5] To examine whether, when and how accurately the
stochastic model predicts the behavior in actual aquifers, we
propose numerical quantitative estimations of the quasi-
ergodic behavior in the spirit of the ‘‘operational ergodic-
ity’’ of Kabala and Sposito [1994]. We use a rigorous
mathematical proof of the existence of the upscaled macro-
dispersion process to check whether the numerical method
is accurate enough for our purpose. Further, we use the
same theoretical result to investigate the ergodicity. That is
to say, we are looking for those indicators of the contam-
1
Institute of Applied Mathematics, Friedrich-Alexander University of
Erlangen-Nuremberg, Erlangen, Germany.
2
‘‘T. Popoviciu’’ Institute of Numerical Analysis, Romanian Academy,
Cluj-Napoca, Romania.
3
Institute of Agrosphere (ICG-IV), Research Center Ju¨lich, Ju¨lich,
Germany.
Copyright 2006 by the American Geophysical Union.
0043-1397/06/2005WR004546
W04409
WATER RESOURCES RESEARCH, VOL. 42, W04409, doi:10.1029/2005WR004546, 2006
1 of 17
ination in actual plumes that can be predicted by the
theoretical result, in the limits of acceptable errors. The
numerical task is carried out with the ‘‘global random walk’’
algorithm (GRW) [Vamos¸ et al., 2003]. Superseding the
limitations encountered by the classical particle tracking
method, the GRW algorithm performs the simulation of
advection-dispersion displacements over thousands of ad-
vection timescales of tens of billions of particles, initially
distributed over hundreds of heterogeneity scales.
[6] For the sake of clarity and for computational reasons,
we consider only the classical stochastic model, two-dimen-
sional transport problems and small variance of the velocity
field in our numerical investigations. This permits detailed
analyses of transport in velocity fields with finite correlation
lengths, which are not possible in less restrictive conditions.
In the following we discuss these limitations and present
some reasons to choose this methodological frame.
[7] More recently, stochastic models free of the Fickian
hypothesis were proposed as alternatives to the classical
stochastic model. These models generally describe ensem-
ble averaged concentrations and their parameters can, in
principle, be derived from measurable properties of the
medium (breakthrough curves or hydraulic conductivities).
The theory proposed by Cushman and Moroni [2001]
generalizes the statistical mechanical approach for Hamil-
tonian systems to an arbitrary statistical ensemble with
invariant probability measure and it enables the description
of anomalous dispersion induced by velocity fluctuations
evolving over a hierarchy of scales on the scale of obser-
vation. Berkowitz and Scher [2001] presented a general
approach, based on the ensemble average of a master
equation, equivalent to a continuous time random walk over
a range of length scales on which the statistical homogene-
ity can be assumed. This approach was extended to a
Fokker-Planck equation with memory term, which integra-
tes the transport behavior over homogeneous units of the
medium, each of them described by the ensemble averaged
master equation [Cortis et al., 2004]. Another promising
stochastic model for transport in saturated porous media
describes the trajectories of the solute particles, in the
position-velocity phase space, by means of a Langevin
stochastic differential equation which accounts for non-
Fickian features of the transport as well [Kurbanmuradov
et al., 2003]. The ensemble average of the classical advec-
tion-dispersion equations can be obtained under appropriate
limiting conditions as particular cases of the three stochastic
models mentioned above.
[8] The classical model is also able to capture the
complexity of transport in groundwater. Anomalous diffu-
sion in the preasymptotic transport regime presumably
occurs under quite broad conditions [Trefry et al., 2003].
For more restrictive conditions (perfectly stratified aqui-
fers), the ensemble averaged solution of the classical model
can be superdiffusive at all times [Matheron and de Marsily,
1980]. However, what really makes the classical model
attractive for applications is the existence of mathematical
proofs for the diffusive behavior of the ensemble averaged
concentration, for typical groundwater transport problems,
characterized by velocity fields with finite correlations
lengths. Since the classical model describes the transport
in given velocity fields, the quantitative approach to ergo-
dicity is straightforward and consists in analyzing the
sample-to-sample fluctuations and the deviation of the
ensemble averaged solutions from the macrodispersion
process. Therefore, in the following we consider only the
model based on the advection-dispersion mechanism. The
conclusions of such an investigation could be useful in
assessing the predictive power of the stochastic modeling
which does not assume the Fickian behavior at a local scale.
[9] Two-dimensional models can be very useful tools to
predict contamination in natural aquifers [Hassan et al.,
1998]. They may be applied to the case of hydraulic
conductivity which is isotropic in the horizontal plane
but has a much smaller correlation length in the vertical
direction as well as to transport at regional scale [Dagan,
1987]. Two-dimensional numerical simulations in vertical
planes oriented along the mean flow direction were suc-
cessful in reproducing the experimental results obtained in
tracer tests [Moltyaner et al., 1993]. The two-dimensional
simulations also provide insight into the convergence of the
transport process to the Gaussian limit [Trefry et al., 2003].
This is mainly relevant when ergodicity is investigated, as
suggested by Dagan [1984, 1987], through averages over
space domains with large transverse dimensions. Then, the
behavior of the space mean concentration is mainly gov-
erned by the longitudinal dispersivities which tend in
ensemble mean to the same limit and are quite close after
a few tens of heterogeneity scales, for both two- and three-
dimensional transport in aquifers with moderate variability
of the hydraulic conductivity [Dagan, 1987; Fiori et al.,
2003; Dentz et al., 2002, 2003].
[10] The presumption that the hydraulic conductivity has
a lognormal distribution and small variance is an accepted
simplification leading to Gaussian velocity fields which can
capture the essential features of the stochastic model [Cortis
et al., 2004; Eberhard, 2004]. In the present context, this
choice is not a limitation since the existence of the upscaled
Gaussian distribution for the ensemble averaged concentra-
tion, used as reference in our investigations, is ensured for
Gaussian velocity fields when the velocity variance goes to
zero.
[11] The paper is structured as follows. Section 2 contains
some definitions and results concerning the notions of
macrodispersion and ergodicity. Section 3 contains some
general considerations on numerical modeling as well as the
statement of our numerical method. Section 4 presents the
main numerical results. Conclusions are drawn in section 5.
Appendix A presents the GRW algorithm and details on the
numerical computations. Appendix B is dedicated to tech-
nical details concerning the transport problem which ap-
proximately fulfils the theoretical requirements for the
existence of upscaled Gaussian diffusion and to the statistics
of numerically generated velocity fields.
2. Macrodispersion and Ergodicity
2.1. Macrodispersion Stochastic Model
[12] For slowly variable porosity which can be taken as a
constant, and for nonreactive solutes, the mathematical
model of transport in saturated porous media is given by
an advection-dispersion equation for the concentration field
c(x, t),
@
t
c þ Vrc ¼ Dr
2
c: ð1Þ
2 of 17
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
The constant ‘‘local dispersion coefficient’’ D accounts for
both the molecular diffusion and the hydrodynamic mixing
due to the small-scale variability of the velocity field
[Sposito et al., 1986; Kapoor and Gelhar, 1994; Labolle
and Fogg, 2001; Jankovic´ et al., 2003]. The stochastic
approach considers stationary velocities V(x) which are
realizations of a random field (random space function) that
is statistically homogeneous.
[13] A stochastic process has ‘‘diffusive behavior’’ when
the mean squared displacement, or variance, is linear as
function of time. The typical example is the Brownian
motion, i.e., the one-dimensional Gaussian process X(t)
with zero mean and variance s
2
(t)= hX
2
i (t)=2Dt. In this
case, the diffusion coefficient D describes both the shape of
the Gaussian distribution and the width of the diffusion
front hX
2
i. In the case of the process described by the
equation (1), and, generally, in systems with space-time
variable properties, the width of the diffusion fronts is no
longer given by the diffusion coefficients alone. Instead, the
rate of increase of the second centered moment s
2
(t) defines
‘‘effective diffusion coefficients’’
D
eff
t ðÞ¼
s
2
t ðÞ
2t
; ð2Þ
which can be used to check whether the process has a
diffusive behavior [Avellaneda and Majda, 1992]. The
existence of a constant limit for t !1 of the effective
coefficients (2) is ensured, for every realization of the
velocity field, by some weak conditions on the partial
derivatives of V (bounded vector potential) [Avellaneda and
Majda, 1989; Tatarinova et al., 1991].
[14] The definition (2) was mainly used in hydrological
applications for comparisons between theory and field tests
[Dagan, 1987] and to compute the first-order approxima-
tions of dispersivities by traveltime statistics [Vanderborght
and Vereecken, 2002; Ferna`ndez-Garcia et al., 2005a,
2005b]. Sometimes, in analytical approaches [Attinger et
al. , 1999; Dentz et al. , 2000a] or numerical studies
[Tompson and Gelhar, 1990; Trefry et al., 2003] the
dispersion is described by the local slope, ds
2
/dt or by the
mean slope of the variance
e
D t ðÞ¼
s
2
t ðÞ s
2
0 ðÞ
2t
¼ D
eff
t ðÞ
s
2
0 ðÞ
2t
: ð3Þ
At large times
e
D tends to D
eff
and the slope of the variance
can be used as well to define the asymptotic effective
coefficients. However, the time behavior of the solute plume
is properly described by the rate of increase of the second
moment (2), which, unlike the mean slope (3), is positive-
definite at all times.
[15] Since disconnected, bimodal or asymmetric plumes
could have the same second moment as a Gaussian plume,
the existence of the coefficients (2) do not yet prove the
existence of the Gaussian limit process. A limit theorem was
demonstrated by Kesten and Papanicolaou [1979]. Neglect-
ing the local dispersion in (1), considering velocity fields
with nonvanishing mean U and small fluctuations eu, e
1, V = U + eu, and making the assumption that the field has
some ‘‘strong mixing’’ property, as characterized by a
suitable fast decay of the correlation function r(x)/e
2
=
hu(0)u(x)i, the authors proved that the average of the
transport process over the realizations of the velocity field
can be upscaled to a Gaussian diffusion. In the ‘‘weakly
random limit’’ (e ! 0, t !1, e
2
t = constant), the ensemble
averaged concentration verifies an advection-diffusion
equation with the upscaled diffusion coefficients given by
D* ¼
Z
1
0
rUt ð Þdt :
Under these conditions, Taylor’s [1921] statement gets a
rigorous proof.
[16] Using a scaling argument, Winter et al. [1984] show
that when the local coefficients D are of the order e
2
, an
extension of the exact result of Kesten and Papanicolaou
[1979] to advection-diffusion processes is possible and the
upscaled diffusion coefficients are of the form
D* ¼ D þ
Z
1
0
rUt ð Þdt : ð4Þ
For incompressible velocity fields, the upscaled velocity
equals the mean velocity U. We remark here that the
upscaled coefficient (4) has the same form as the first-order
approximation of the ‘‘macrodispersion coefficient’’ derived
by Dagan [1984]. A test for the accuracy of our numerical
simulations which mimic the above conditions is to check
whether the deviations of the computed upscaled coeffi-
cients from the theoretical values given by (4) are one order
of magnitude smaller than the local dispersion coefficient.
2.2. Ergodicity
[17] Strictly speaking, the macrodispersion model is di-
rectly applicable to a single aquifer if the same asymptotic
Gaussian approximation holds for each realization [Sposito
et al. , 1986]. Hereafter we call this strong property
‘‘asymptotic ergodicity.’’ Sufficient conditions for asymp-
totic ergodicity are provided if the ensemble average of the
actual concentrations tend to the solution of the macro-
dispersion model and the sample-to-sample fluctuations tend
to zero. In this case, the effective coefficients in every
realization necessarily also tend to the macrodispersion
coefficients.
[18] The assumption underlying the concept of ergodicity
is that suitable space averages of the actual concentration
can be described by the solution of the macrodispersion
model, when the spatial variability of the concentration
encompasses the variability from realization to realization.
Dagan [1984, 1987] assumed that the space and the
ensemble mean are interchangeable if the variance of the
space averaged concentration tends to zero. It was shown
that the spatially averaged concentration in single realiza-
tions is close to the macrodispersion solution after tens of
correlation lengths of the hydraulic conductivity, if the
initial solute body or the domain of the space average
extends over a few correlation lengths across the mean flow
[Dagan, 1984, p. 165]. This result was derived in a
Lagrangian frame for a normal distribution of the displace-
ments of the solute particles, inferred by a first-order
approximation of the transport equations. A similar ap-
proach, using another approximation technique [Dagan,
1987, equation (3.14)], led however to a different result
indicating a nonergodic behavior, manifested by a finite
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY
3 of 17
W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
limit of the concentration fluctuations at the plume center of
mass, which is practically reached after a few correlation
lengths [Dagan and Fiori, 1997, Figure 3]. The same
approximation led to a finite asymptotic variance of the
longitudinal second moment of the plume [Fiori, 1998,
Figures 4 and 5], which also indicates nonergodic behavior.
These results contradict other Lagrangian approaches
[Pannone and Kitanidis, 1999] and Eulerian theories (which
are limited in turn by closure approximations) [Kapoor and
Gelhar, 1994; Kapoor and Kitanidis, 1998]. A renewed
Lagrangian result of Fiori and Dagan [2000] corrects the
previous one, showing that the concentration coefficient of
variation tends to zero (after much larger traveltimes,
corresponding to thousands of correlation scales) but pre-
dicts a large time behavior which is still different from that
obtained in Eulerian approaches.
[19] The ergodic behavior of individual plumes was
recently investigated numerically. Jankovic´ et al. [2003]
found a good agreement between the behavior of the
individual plumes and the predictions of the macrodisper-
sion model, for both two- and three-dimensional transport
simulations, which suggests an ergodic behavior of the
extended plumes as predicted by Dagan [1984]. On the
contrary, Trefry et al. [2003], who simulated the concentra-
tion field in single realizations of transport for two-dimen-
sional plumes with initial extension of tens of correlation
lengths across the mean flow, found that the ergodic
behavior is not reached even when the plumes have traveled
hundreds of correlations lengths. The last result is appar-
ently in agreement with that found by Fiori and Dagan
[2000]. Yet, the two numerical approaches are not com-
pletely comparable. Trefry et al. [2003] performed a direct
investigation, by comparing the dispersivities in given
realizations, defined by the effective coefficients (3) divided
by the mean velocity, to the macrodispersivities corres-
ponding to the upscaled coefficients (4) of the stochastic
model. They also checked whether the concentration
becomes Gaussian. Jankovic´ et al. [2003] did not investi-
gate the concentration field and in their paper the macro-
dispersivities were computed through a Taylor formula,
similar to (4) [see Jankovic´ et al., 2003, equation (14)]
which gives smoother results and can explain the different
conclusions in the two papers.
[20] In spite of the common belief that extended initial
plumes have an ergodic behavior, there are no general
mathematical proofs for that. The relevance of the macro-
dispersion model for a single realization of the medium is
not an immediate consequence of the existence of the
Gaussian upscaling and is not ensured either by the ergodic
properties of the random velocity field [Sposito et al., 1986;
Kabala and Sposito, 1994]. Moreover, in the limiting case
of vanishing local dispersion, it was proved that groundwa-
ter flows governed by the Darcy law are geometrically
constrained against becoming chaotic and, consequently,
the purely advective transport is in general not ergodic
[Sposito, 1997, 2001]. The stratified aquifer model of
Matheron and de Marsily [1980] provides another counter
example, for the complete advection-dispersion model. The
exact expressions derived by Clincy and Kinzelbach [2001]
show that the fluctuations of the longitudinal effective
coefficient tend to a finite limit, even if the flow is not
aligned with the strata and the transport has asymptotic
diffusive behavior. Since the asymptotic effective coeffi-
cients in given realizations do not tend to their expectation
(a necessary condition for ergodicity), the transport is not
ergodic. Quantitative numerical estimations for the devia-
tion from the macrodispersion model of the transport in
realizations of a Gaussian velocity field with finite correla-
tion lengths (in both transverse and longitudinal directions)
were recently obtained, for the case of point-like initial
plumes [Eberhard, 2004]. The numerical computations of
Eberhard [2004], based on approximate solutions of the
transport equation, show that the fluctuations of the longi-
tudinal effective coefficient tend to zero, indicating that a
necessary condition for ergodicity is fulfilled. To prove the
ergodicity of transport in velocity fields with finite correla-
tion lengths it is necessary to check whether the sufficient
conditions, formulated for actual concentrations, are veri-
fied. It is also useful to investigate the behavior of the
effective coefficients for extended plumes and their relation
with ergodicity.
[21] To achieve our goals, we introduce a definition
which is general enough to include various conceptual
features of the term ergodicity cited above. Let A(t) be the
value of an observable at time t, A*(t) the theoretical
prediction and DA = h(A A*)
2
i
1/2
the root mean square
deviation from theory (the angular brackets denote averag-
ing over the ensemble of realizations of the velocity filed).
The observable A is ergodic within the range h, h > 0, if
DA h. Since the theory usually predicts A* as an asymp-
totic limit and the observable A is known at finite times, in
practice a more flexible definition is necessary. From the
relation (DA)
2
= s
A
2
+(DhAi)
2
, where s
A
= h(A hAi)
2
i
1/2
is
the standard deviation of A and DhAi = hAi A* the
deviation of the mean hAi with respect to A*, one obtains
the equivalent definition
Ergodicity condition e1
D A h i j j h
1
Ergodicity condition e2
s
A
h
2
:
ð5Þ
If the conditions in (5) are fulfilled, then the observable A is
ergodic within the range h =(h
1
2
+ h
2
2
)
1/2
. The asymptotic
ergodicity in a strict sense, i.e., as discussed by Sposito et al.
[1986], is ensured when the observable A is the concentra-
tion and (5) holds in the large time limit for arbitrary small
and positive h
1
and h
2
. In other words, asymptotic
ergodicity requires the existence of an upscaled macro-
dispersion solution for the ensemble averaged concentra-
tion, as expressed by condition e1, and a ‘‘self-averaging’’
property, namely the convergence of the actual concentra-
tions to their ensemble average, described by condition
e2. When the observable A is a space averaged concentra-
tion and the conditions (5) are fulfilled at finite times for
increasing dimensions of the averaging domain, then the
space averaged concentration converges in the mean square
limit to the ensemble averaged concentration and the space
and ensemble averages are interchangeable. Condition e2
for space averaged concentration was used by Dagan [1984,
4 of 17
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
1987] to investigate the ergodicity under the assumption
that the ensemble averaged concentration is already
Gaussian at finite times. Condition e2 also formulates the
self-averaging property of the longitudinal effective coeffi-
cient, investigated by Clincy and Kinzelbach [2001] and
Eberhard [2004]. The condition e1 alone, applied to
effective coefficients or to second moments of the plume,
is often referred to as ‘‘ergodicity condition’’ [Fiori, 1998;
Naff et al., 1998b; Zhang and Seo, 2004; Dagan, 2004]. For
finite times and ranges h, the definition (5) corresponds to
the operational interpretation for ergodicity, proposed by
Kabala and Sposito [1994], which seeks conditions that
lead to acceptably small deviations of the experimentally
observable concentrations from the predictions of the
stochastic model. The observables A used to quantify the
ergodicity in this study are the cross-section space-averaged
concentration and the effective diffusion coefficients. The
corresponding theoretical predictions A* are the solutions of
the upscaled macrodispersion model.
3. Numerical Approach
3.1. Prerequisites for Numerical Investigations
[22] Numerical investigations on the existence of the
diffusion limit and ergodicity require the simulation of large
ensembles of realizations of the transport over large traveled
distances. Because this task generally surpasses the available
computing resources, some compromise was accepted and
the efforts focused by now on one or two of the three
objectives: (1) accurate concentration field in a given reali-
zation, (2) simulations of the transport over large timescales,
and (3) reliable statistics for the ensemble of realizations.
[23] For instance, objective 1 was pursued by Smith and
Schwartz [1980], Tompson and Gelhar [1990], Moltyaner et
al. [1993], Naff et al. [1998b], Kapoor and Kitanidis
[1998], Trefry et al. [2003], and Jankovic´ et al. [2003],
where given realizations of the transport were simulated
over distances still too small to describe the asymptotic
behavior (with the exception of the confined aquifer case
considered by Kapoor and Kitanidis [1998]). Reliable
statistical ensembles (objective 3), aiming to compute the
upscaled coefficients (4), were obtained by tracking one
particle in 1500 realizations of the velocity field [Bellin
et al., 1992], 20 particles in 500 realizations [Salandin
and Fiorotto, 1998], 20 particles in 1600 realizations
[Zhang and Seo, 2004], or 10000 particles in 20 realization
[Ferna`ndez-Garcia et al., 2005a], and averaging over
realizations. The computational constraints did not allow
the simulation over more than a few of velocity correlation
lengths and the asymptotic regime, i.e., constant coefficients
D*, was not reached.
[24] Objectives 2 and 3 were aimed at by Schwarze et al.
[2001] and Dentz et al. [2002, 2003], in the study of the large
time behavior of the effective coefficients in the case of point-
like injection. The velocity fluctuations eu were numerically
approximated using the generator of Kraichnan [1970]. The
coefficients were computed in two ways: as average over the
trajectories of a single particle in thousands of realizations,
which corresponds to D* in (4), and as average over realiza-
tions of the coefficients corresponding to D
eff
(2), obtained by
tracking a number of particles (between tens and hundreds)
released at the same point in each realization. The time
behavior of the two coefficients is different but both tend to
the same constant limit, after travel distances of more than
1000 correlations lengths. The numerically simulated limit
coefficients were close to the approximations of the order e
2
provided by perturbation analyses.
[25] The objective 1 has not yet been attained at the same
time with 2 and 3 due to the limitations of the computational
resources. For instance, with regard to the particle tracking
method which is the most frequently used in large-scale
simulations, it is recognized that the number of particles
should be enormous if we want to obtain accurate concen-
trations in the limit of a few significant figures [Sun, 1996,
p. 95]. Recent investigations show that convergent simu-
lations for single realizations of the transport problem
considered here require a number of particles of the order
of tens of billions, which is prohibitive for the particle
tracking algorithms [Vamos¸ et al., 2003; Suciu et al., 2004].
[26] Our numerical approach uses the GRW algorithm, a
generalized random walk method for which there are no
limitations as to the maximum number of particles. Suciu et
al. [2002, 2004] showed that GRW is appropriate to
simulate large-scale transport in groundwater. In this paper
we show that GRW fulfils the requirements 1, 2, and 3 from
above and makes possible a direct investigation on the
ergodicity issue, based on the definition (5).
3.2. Implementation of the GRW Method
[27] We considered two-dimensional divergence-free ve-
locity fields with constant mean hVi = U =(U, 0), U = 1 m/d,
given by the Darcy law for normal log hydraulic conduc-
tivity y, with isotropic correlation length l
y
= 1 m, and
variance s
y
2
= 0.1. They were generated numerically with
the Kraichnan procedure which is frequently used in large-
scale simulations of transport [Jaekel and Vereecken, 1997;
Schwarze et al., 2001; Dentz et al., 2002, 2003; Eberhard,
2004]. The Kraichnan method ensures the incompressibility
and approximates the realizations of the Gaussian velocity
field, given by the first-order approximation in s
y
of Darcy
and continuity equations, by means of a superposition of
sine (or cosine) periodic modes.
[28] In every realization of the velocity field the simula-
tion of an isotropic diffusion (D
ll
= D = 0.01 m
2
/d, D
lm
=0
for l 6¼ m) was conducted for dimensionless times Ut/l
y
corresponding to thousands of correlation lengths, using the
reduced fluctuations GRW algorithm presented in section
A1. In each realization of the velocity field, N = 10
10
particles were initially uniformly distributed in a vertical
band or located at the origin of the grid. Because the
theoretical result presented in section 2.1 was obtained for
unbounded domains, the grid was chosen to be larger than
the maximum extension of the plume. Implementation
details are discussed in section A2. Some preliminary tests
have shown that the longitudinal local dispersion has little
influence on the resulting concentration fields, in agreement
with the results concerning the longitudinal effective coef-
ficients presented by Fiori [1996, Figure 1] and by Dentz et
al. [2000a, Appendix B]. Thus D mainly describes the
strength of the transverse local dispersion.
[29] We chose the parameters s
y
2
= 0.1 and D = 0.01 m
2
/d
after preliminary investigations, presented in section B1.
These parameters minimize possible numerical artifacts
occurring in two-dimensional simulations using particles
methods and Kraichnan generator. The comparison pre-
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY
5 of 17
W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
sented in section B2 shows that the fluctuations have similar
large time behavior for exponential and Gaussian shape for
the correlation of the log hydraulic conductivity field,
provided that the number of periods Np in Kraichnan
routine and the number of realizations R are large enough.
For the detailed investigations on ergodicity presented in the
following we chose the exponential correlation shape, Np =
6400 and R = 256 realizations. In section B3 it was shown
that the numerically generated velocities approximate a
Gaussian random field and the squared fluctuations e
2
have
roughly the same order of magnitude as the local dispersion
coefficient D (see second row of Table B1), in agreement
with the theoretical assumptions leading to an upscaled
Gaussian process. The corresponding theoretical values of
the longitudinal and transverse effective diffusion coeffi-
cients, computed from the first-order estimation of the
correlations under the integral in (4) [Dagan, 1984], have
the values
D
11
*
¼ 0:11 m
2
=day D
22
*
¼ 0:01 m
2
=d: ð6Þ
[30] The simulations started with a uniform initial distri-
bution of the 10
10
particles inside rectangles l
y
Ll
y
, with
L = 10, 50 and 100, oriented across the mean flow. In the
case of instantaneous point source, all the particles were
released from the origin of the grid at the initial time.
The total simulation time was of 2000 days for L = 100,
2700 days for L = 50, and 4000 days for L = 10 and for the
point source. The cross-section space-averaged concentra-
tions in given realizations C(x
1
, t; L) were computed
according to (A1) by the number of particles in a domain
l
y
(160 + L)l
y
divided by the total number of particles
N. The normalized concentrations were obtained through
division by the initial concentrations: C(x
1
, t)= C(x
1
, t; L)/
C(0, 0; L). The averaging domain, oriented across the mean
flow and centered at x
1
, corresponds to the ideal sampling
across reference planes aimed at in field tracer tests
[Vanderborght and Vereecken, 2002] as well as in laboratory
and numerical experiments [Ferna`ndez-Garcia et al., 2005a,
2005b]. We stored the space averaged concentration at the
plume center of mass, C(hx
1
i, t), as well as the concentration
field C(x
1
, t) at several fixed times. The effective coeffi-
cients (2) were computed according to (A2). Further, we
computed the mean concentrations hCi and the mean
effective coefficients hD
ll
eff
i, as averages over the ensemble
of R realizations of the transport. Since in all the simulations
the mean flow velocity was U = 1 m/d, the numerical values
of the effective coefficients coincide with the dispersivities
D
eff
/U, measured in meters.
4. Numerical Results on Ergodicity
4.1. Cross-Section Averaged Concentrations
[31] Since the width of the averaging domain is larger
than the transverse dimension of the plume, the transport
can be described by a one-dimensional problem for the
spatially averaged concentration C(x
1
, t). The corresponding
theoretical concentration C* is the one-dimensional Gauss-
ian distribution of coefficient D*
11
= 0.11 m
2
/d which
describes the cross section average of the concentration
predicted by the stochastic macrodispersion model. Ergo-
dicity condition e1 in (5) was investigated by means of the
relative deviation of the ensemble average of the cross-
section averaged concentration from the theoretical value,
DhCi/C*=(hCi C*)/C*. The deviations DhCi/C*,
computed at the plume center of mass hx
1
i presented in
Figure 1 show that in the case of point source at t = 4000l
y
/U
condition e1 is fulfilled within a range h
1
’ 0.17C*. The
increase of L reduces the deviations at early times. For
instance, for L 50, h
1
’ 0.13C* at 100 advection
timescales l
y
/U (or, equivalently at one dispersion time-
scale l
y
2
/D).
[32] To check the second ergodicity condition e2 we used
the standard deviation of the cross-section averaged con-
centration divided by the theoretical solution s
C
/C*, where
s
C
=(hC
2
ihCi
2
)
1/2
. The results presented in Figure 2
show that for point source at t = 4000l
y
/U, condition e2
is fulfilled within a range h
2
’ 0.11C*. Thus, according
to (5) the ergodicity range is h =(h
1
2
+ h
2
2
)
1/2
’ 0.2C*.
The monotonous decay of the mean square deviation
Figure 1. Deviations from the macrodisperion model of
the mean concentration at the plume center of mass DhCi/
C*(hx
1
i, t).
Figure 2. Standard deviations of the concentration at
plume center of mass s
C
/C*(hx
1
i, t).
6 of 17
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
(DC)
2
= s
C
2
+(DhCi)
2
indicate the convergence in mean
square limit of the cross-section averaged concentration to
the solution of the macrodispersion process.
[33] For L = 100 at one dispersion timescale one obtains
in a similar way h
2
’ 0.1C* and h ’ 0.16C*. The results for
ensemble averaged concentration as function of the distance
from the center of the initial plume for a fixed time of
100l
y
/U, presented in Figures 3 and 4, indicate that h has a
minimum at x
1
= 100l
y
(which corresponds to the mean
center of mass coordinate). A rough prediction of the
maximum concentration estimated by the cross-section
averaged concentration at the plume center of mass, C*±
3h, is reliable for almost all the realizations of the ensemble
defined by a given statistical structure of the log hydraulic
conductivity. In the somewhat ideal case analyzed here, for
initial plumes extending over 100 correlation lengths across
the mean flow, the prediction based on the macrodispersion
model and the ergodic hypothesis of the space averaged
concentration C in actual aquifers is affected by uncertain-
ties of the order 6h ’ 0.9C*.
[34] The behavior of standard deviations presented in
Figure 2 indicate a nonlinear dependence on the plume
dimension at early times. For large sources (L 50) the
concentration fluctuations s
C
/C* are smaller than in the case
of point source, decrease with L and remain almost constant,
at values about 0.1, over more than thousand advection
timescales. The monotonous decrease of the fluctuations,
i.e., the asymptotic ergodicity, does not occur until the end
of the simulations.
4.2. Effective Coefficients
[35] The relative deviation of the ensemble averaged
effective coefficients from the theoretical values given in
(6) can be expressed using (3) by
DhD
eff
ll
i=D
ll
*
¼ hD
eff
ll
i D
ll
*
=D
ll
*
¼ Dh
e
D
ll
i=D
ll
*
þ s
2
0;ll
= 2tD
ll
*
: ð7Þ
Since the contribution of the initial plume s
0,ll
2
/(2tD
*
ll
) is a
deterministic quantity which tends to zero for large times, to
compare the asymptotic behavior of the effective coefficients
for different L we used the deviation of the mean slope of the
second centered moment of the plume Dh
e
D
ll
i/D
*
ll
. The results
for longitudinal coefficients are presented in Figure 5.
Ergodicity condition e1 for the mean slope is fulfilled within
a range h
1
’ 0.11D
*
11
in the case of point source at 4000
advection times and within a range h
1
’ 0.05D
*
11
in the case
L = 100 at 100 advection times. Because s
0,11
2
= 0.1 m
2
in all
cases, at t = 100l
y
/U its contribution in (7) is already only
0.0045 and therefore the same range h
1
characterizes the
longitudinal effective dispersion coefficients.
[36] Because
e
D
ll
and D
ll
eff
differ by a deterministic quan-
tity, as shown by (3), they have the same fluctuations. The
fluctuations s
D
11
eff
/D*
11
of the longitudinal effective disper-
sion coefficient are given in Figure 6. For a small increase
of the plume dimension, the fluctuations increase (similarly
to the results reported by Naff et al. [1998b, Figure 15]) and
they decrease when the plume dimension is further in-
creased, like the fluctuations of the cross-section averaged
Figure 3. Deviations from the macrodisperion model of
the mean concentration as space function DhCi/C*(x
1
,
100l
y
/U).
Figure 4. Standard deviations of the concentration as
space function s
C
/C*(x
1
, 100l
y
/U).
Figure 5. Deviations from the macrodisperion model of
the mean slope of the longitudinal second moment D h
e
D
11
i/
D
11
*.
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY
7 of 17
W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
concentration in Figure 2. Ergodicity condition e2 is fulfilled
within a range h
2
’ 0.16D*
11
for point source at t = 4000l
y
/U
and within a range h
2
’ 0.14 D*
11
in the case L = 100 at t =
100l
y
/U. It follows that, according to (5), the corresponding
ergodicity ranges are h ’ 0.2 D*
11
(point source) and h ’
0.15D*
11
(L = 100). The large time behavior for small sources
(L 10), shown in Figures 5 and 6, indicates the asymptotic
ergodicity of the longitudinal effective coefficient. The
deviations from the macrodispersion model of the mean
slope of the transverse second moment and the fluctuations
of the transverse effective coefficients, presented in Figures 7
and 8, respectively, show that the transverse mean slope
coefficient
e
D
22
behaves asymptotically ergodic. Since s
0,22
2
/
(2tD
22
* ) goes to zero at large times, from (7) and Figure 7
it follows that the effective transverse coefficient is also
asymptotically ergodic.
[37] The asymptotic effective coefficients and their devi-
ation from the theoretical values D*
ll
were computed in the
case L = 100 as follows. First, the temporal averages of the
mean slopes
e
D
ll
between t = 200l
y
/U and t = 2000l
y
/U,
[
e
D
ll
], were computed in every realization. Then, averages
over realizations were used to estimate the mean asymptotic
coefficients D
ll
1
= h[
e
D
ll
]i and the deviations
D
e
D
ll
hi
¼
D
e
D
ll
hi
D
ll
*
2
E
1=2
with respect to D*
ll
. The results are presented in Figure 9 as
functions of the number of realizations R. The fact that the
deviations of the mean coefficients D
1
ll
from D*
ll
are one
order of magnitude smaller than the local dispersion
coefficient D, in very good agreement with (4), constitutes
a test for the accuracy of the numerical simulations. The
behavior of the deviations D[
e
D
ll
] (thin lines in Figure 9)
indicate that the time averaged coefficients [
e
D
ll
] are ergodic
within a range of the order of the local dispersion coefficient
h ’ D = 0.01 m
2
/d.
[38] Even if the mean slope is ergodic within an accept-
able small range, the macrodispersion coefficients D*
ll
Figure 7. Deviations from the macrodisperion model of
the mean slope of the transverse second moment Dh
e
D
22
i/
D
22
*.
Figure 9. Asymptotic effective coefficients and deviations
from the theoretical values (6) for L = 100, as functions of
the number of realizations R.
Figure 8. Fluctuations of the transverse effective disper-
sion coefficient s
D
22
eff
/D
22
*.
Figure 6. Fluctuations of the longitudinal effective
dispersion coefficient s
D
11
eff /D
11
*.
8 of 17
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
describe the time behavior of the second moments of the
plume with errors smaller than 100% only when, according
to (7), s
0,ll
2
/(2tD*
ll
) 1. For the transverse coefficient the
corresponding time s
0,22
2
/(2D) can be very large (tens of
thousands of advection scales for L 50) and increases
proportionally with L
2
. Figures 2, 6, and 8 indicate that the
fluctuations are also characterized by time scales which
increase with the plume dimension L. This finding agrees
with the exact analytical result of Clincy and Kinzelbach
[2001]. Figure 4 in their paper shows that the fluctuations of
the longitudinal coefficient tend to an asymptotic value
(which is finite in the case of the model of G. Matheron
and G. de Marsily studied there) after times of the order of
L
2
/D, i.e., the time for a dispersive spreading of the solute
over the width L of the initial plume. On this timescale the
transport becomes independent of the extent and shape of
the initial plume. This is also indicated by the behavior of
the cross-section concentrations (Figures 1 and 2) and of the
mean slope coefficients (Figures 5 – 8) for L = 10 and point
source, which after 4000 advection times reach almost the
same ergodicity range h. Similar time behavior was found
by Dentz at al. [2000b, Figure 1] for the ensemble averaged
longitudinal dispersion coefficient of extended plumes (de-
fined as half the derivative of the second moment). The total
time in our simulations is still too small to allow an
estimation of the asymptotic ergodicity timescale. However,
the numerical findings and the two theoretical results
mentioned above indicate that the ergodicity scale increases
like L
2
.
4.3. Time Behavior of the Effective
Coefficients and Ergodicity
[39] As shown by equation (A3) in section A1, the GRW
procedure computes the variance of displacements X in a
single realization as an average over the realizations of the
local dispersion process and over the distribution of the
particles inside the initial plume. The ensemble average of
the variance can be written as
hs
2
ll
i
V
¼ hhX
2
l
i
D;X
0
hX
l
i
2
D;X0
i
V
¼ hX
2
l
i
D;X0 ;V
hX
l
i
2
D;X0;V
hhX
l
i
2
D;X0
i
V
hX
l
i
2
D;X0 ;V
;
ð8Þ
where the subscripts D, X
0
and V denote the average over
the realizations of the local dispersion, the initial distribu-
tion of the particles, and the realizations of the random
velocity field respectively. This obvious relation shows that
hs
ll
2
i
V
is the difference between the variance with respect to
the ensemble average of the center of mass R
l
= hX
l
i
D,X
0
and
the variance of the center of mass R
ll
= hR
l
2
i
V
hR
l
i
V
2
.
Assuming that the averages over initial positions and
velocity realizations have the following permutation prop-
erty h i
D,X
0
,V
= h i
D,V,X
0
, the first term in (8) becomes
hhX
2
l
i
D;V
hX
l
i
2
D;V
i
X0
þ hhX
l
i
2
D;V
i
X0
hX
l
i
2
D;X0;V
:
For statistical homogeneous velocity fields it seems reason-
able to assume the independence of the averages h
e
X
l
i
D,V
and h
e
X
l
2
i
D,V
, where
e
X
l
= X
l
X
0l
, from the initial state X
0
.
Then, the last two terms above give the initial variance
s
0,ll
2
= hX
0l
2
i
X
0
hX
0l
i
X
0
2
and the first term becomes
independent of X
0
and represents the variance around the
ensemble averaged center of mass for the process starting
with a point-like injection at
e
X
l
= 0,
e
X
ll
= h
e
X
l
2
i
D,V
h
e
X
l
i
D,V
2
.
Finally, the ensemble average (8) of the second centered
moment of the plume becomes
hs
2
ll
i
V
¼ s
2
0;ll
þ
e
X
ll
R
ll
: ð9Þ
The variance
e
X
ll
is just the ll component of ‘‘the second
spatial covariance of an ergodic plume’’ frequently used in
investigations on ergodicity [Zhang and Lin, 1998; Zhang
and Seo, 2004]. When the local dispersion is neglected, (9) is
identical to relation (11) of Dagan [1990], derived in the
hypothesis of ‘‘Lagrangian stationarity’’ and using the
permutation of averages over initial states and velocity
realizations. The explicit dependence of
e
X
ll
on local
dispersion and velocity correlations can be derived from
descriptions of the transport process in terms of trajectories
[Rajaram and Gelhar, 1993a; Fiori, 1998] or by using the
advection-dispersion equation [Kitanidis, 1988; Rajaram
and Gelhar, 1993b].
[40] Dividing both sides of (9) by 2t and using (2), the
ensemble averaged effective coefficients can be written as
hD
eff
ll
i
V
t ðÞ
s
2
X0;ll
2t
¼ D
erg
ll
t ðÞ D
cm
ll
t ðÞ: ð10Þ
The ‘‘center of mass coefficient’’ D
ll
cm
= R
ll
/(2t) corresponds
to the ‘‘pseudodispersivity’’ investigated numerically by
Naff et al. [1998b]. The first term in the right side of (10),
D
ll
erg
=
e
X
ll
/(2t) is the so-called ‘‘ergodic coefficient’’ which
is expected to become constant and equal to the upscaled
macrodispersion coefficient in the large time limit [Dagan,
1990]. The condition D
ll
cm
= 0 is referred to as ‘‘ergodicity
condition’’ [Fiori, 1998; Naff et al., 1998b; Zhang and Seo,
2004; Dagan, 2004]. The ergodic coefficients correspond to
the ‘‘ensemble coefficients’’ introduced by Attinger et al.
[1999] (where the ensemble and effective coefficients are
defined by the time derivative of the corresponding
variances related by (8)). Under the assumption of statistical
homogeneity of the log hydraulic conductivity the ensemble
coefficients were also shown to be independent of the shape
and dimension of the initial plume and equal to the
coefficients for the case of a point source [Dentz et al.,
2000b; Clincy and Kinzelbach, 2001]. Zhang and Seo
[2004] have shown that, even for statistically anisotropic
aquifers, the longitudinal and transverse ergodic second
moments given by theory can be retrieved in numerical
simulations by using the relation (9). However, analyzing
the fluctuations of the longitudinal effective coefficient from
realization to realization, Naff et al. [1998b] found a
‘‘nonergodic behavior,’’ indicated by the increase of the
fluctuations with the transverse dimension of the plume, and
suggested that the approach to a quasi-ergodic state is more
complicated than described by the equation (10). In the
following we investigate this issue in the light of our GRW
simulations.
[41] The longitudinal center of mass coefficient decreases
with the plume dimension (Figure 10) and the ergodic
coefficient derived from (10) is practically independent of
the plume dimensions (Figure 11). This agrees with the
results for the ergodic second moments of the plume
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY
9 of 17
W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
obtained by Zhang and Lin [1998, Figure 9] and Zhang and
Seo [2004, Figure 5] and the results for the longitudinal
center of mass coefficient of Naff et al. [1998b, Figure 14].
D
ll
cm
is the difference between the ergodic coefficient D
ll
erg
and the ensemble average of the mean slope (3), h
e
D
ll
i
V
,
given by the left side of (10). A small center of mass
coefficient D
ll
cm
(t) h
1
is equivalent to ergodicity condition
e1 in (5), for the deviation of the ensemble average of the
observable A(t)=
e
D
ll
(t) from the theoretical prediction
A*(t)= D
ll
erg
(t). The increase of fluctuations with the plume
dimension reported by Naff et al. [1998b, Figure 15] simply
means that the vanishing center of mass coefficients in (10)
is not sufficient to ensure the ergodicity of the mean slope of
the second moment of the plume. Even though the center of
mass coefficients become negligible quantities and the
ensemble averaged slope approach the ergodic coefficient,
the sample-to-sample fluctuations can be still large. This
situation is dramatically illustrated in the case of the
transverse coefficients. After 1000 advection times the first
ergodicity condition e1 is fulfilled within a range two orders
of magnitude smaller than the local dispersion coefficient,
as shown by the behavior of D
22
cm
in Figure 12. The
fluctuations of the mean slope
e
D
22
, which are equal to the
fluctuations of the effective coefficients in Figure 8, indicate
that the second condition e2 is not fulfilled within the same
range. Moreover, the larger the transverse dimension of the
plume is, the less ergodic the transverse coefficients are.
[42] We comment here that the coefficient corresponding
to L = 100 in Figure 13 has negative values at early times.
This nonphysical behavior occurs because in the left side of
(10) the contribution of the initial plume was extracted from
the total variance. This ‘‘bad result’’ shows that the defini-
tion of the effective coefficients by the mean slope of the
variance fails to describe the plume at finite times. We also
note that the negative effective transverse dispersivities
obtained, in three-dimensional case, by Zhang et al.
[1996] are due to their definition by the local slope of the
second moment, and are not ‘‘an artifact of the first-order
approximation,’’ as the authors suggest. It is indeed easy to
see that there are no negative values if the dispersivities are
computed by using the positive-definite coefficient (2) and
the variances shown in Figure 2b of the quoted paper. The
same explanation is valid for the negative dispersivities
Figure 10. Time behavior of the longitudinal center of
mass coefficient D
11
cm
.
Figure 11. Time behavior of the longitudinal ergodic
coefficient D
11
erg
.
Figure 12. Time behavior of the transverse center of mass
coefficient D
22
cm
.
Figure 13. Time behavior of the transverse ergodic
coefficient D
22
erg
.
10 of 17
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
obtained for two-dimensional plumes by first-order theoret-
ical analysis [Zhang and Zhang, 1997, Figure 2e] and
numerical simulations [Zhang and Lin, 1998, Figure 10].
[43] The large differences between the coefficient for
point source and the coefficients for L 10 shown in
Figure 13 prove that for plumes with large transverse
dimensions the ergodic coefficient for preasymptotic regime
cannot be defined by (10). Unlike the center of mass
coefficients (Figure 12) which go to zero, the differences
between the ergodic coefficients increase with the plume
dimension. These deviations can be neglected only for small
plumes, as indicated by the three-dimensional simulations
for L 4 of Zhang and Seo [2004, Figures 5b and 5c] and
by the two-dimensional simulations for L 10 of Zhang
and Lin [1998, Figure 9]. (Note that the small differences
reported in these papers are also due to the graphical
representation based on relation (9); this representation also
reduces the differences shown in our Figure 13.) This
dependence on the plume dimension could be a conse-
quence of the inherent nonhomogeneity of the numerically
generated velocity field. However, in practice, the statistical
homogeneity is always an approximation. This is shown,
for instance, by the behavior of the confidence intervals
rendered by the uncertainty of the statistical parameter
estimates in carefully designed laboratory experiments
[Ferna`ndez-Garcia et al., 2005b], which are very similar
to those for the numerical velocity field presented in
Figure B7. However, even if the lack of strict statistical
homogeneity can produce the small differences for the
longitudinal coefficients in Figure 11, it does not explain
the large differences shown in Figure 13. Therefore the very
validity of the relation (10) can be questioned.
[44] The assumption that the average over initial positions
permutes with the average over the realizations of the
velocity field, on which (10) is based, is not always true.
For instance, when local dispersion is neglected and the
transport process is described in a Lagrangian framework,
the above permutation of averages is allowed by the
‘‘simplifying assumption of ergodicity in the plume spatial
moments’’ [Sposito and Dagan, 1994, p. 588] (carrying the
implicit hypothesis of ‘‘dynamically identical’’ solute par-
ticles described by a single statistical ensemble [Sposito,
1997]). It was already shown by Sposito and Dagan [1994]
that in the purely advective case the relation (10) is not
complete if the interdependence between the initial posi-
tions and the velocities of the solute particles is taken into
account. Because for Darcy flows such interdependencies
cannot be ignored (unless some restrictions on the flow
domain are imposed) [Sposito, 1997, 2001], it follows that
in general this relation is not true. Since the dynamical
system approach used in the papers quoted above is
restricted to the advective case, the analysis of the relation
(10) in the case of nonvanishing local dispersion calls for
further investigations based on appropriate methods.
5. Conclusions
[45] The results obtained in the present study, indicate
that the numerical approach was appropriate for a numerical
investigation on the ergodicity of transport in heterogeneous
aquifers. The requirements of accurate concentration in
every realization of the random field, long time simulations
of transport and large statistical ensembles of realizations
were carried out with the GRW algorithm. The large-scale
simulations were also possible owing to the fast generator of
the first order approximated Darcy velocity fields, based on
the Kraichnan routine. The approximation of the Gaussian
random velocity fields with periodic fields was shown to be
reliable for 6400 modes in the Kraichnan algorithm. The
quantitative assessment of ergodicity for space averaged
concentrations and effective diffusion coefficients was done
via comparisons with a rigorous mathematical result on the
existence of the upscaled macrodispersion process which
describes the ensemble averaged concentration asymptoti-
cally. We emphasize that the conclusions drawn below are
valid for small heterogeneity of the aquifer, which allows
comparisons with the mathematical result, and cannot be
simply extrapolated for highly heterogeneous formations.
[46] A numerical evidence for the asymptotic ergodic
behavior of the two-dimensional transport was supplied
by the simulations in the case of point-like and small initial
plumes. The time to reach acceptably small deviations from
the predictions of the macrodispersion model is of
thousands of advection timescales. In rapidly fluctuating
velocity fields with advection times of the order of seconds,
for instance in the case of turbulent transport in atmosphere,
the ergodic behavior manifests after a few hours. Since the
advection times in groundwater are of the order of days,
the ergodic behavior in a strict sense can be expected when
the plume has traveled tens of years. This could be useful in
applications for persistent contaminants, like the long life
radionuclides.
[47] For sufficiently large initial plumes, the macrodis-
persion model can be used to predict the contamination
with errors that could be acceptable in forecasting at smaller
timescales. For instance, when the initial plume extends
over one hundred heterogeneity scales across the direction
of the mean flow, and the solute plume has traveled
hundreds of heterogeneity scales, the ‘‘three sigma’’ rule
indicates that the cross-section averaged concentration and
the longitudinal effective coefficient in given realizations
can be predicted within a range of uncertainty of about
90%. This uncertainty remains almost constant over
thousands of heterogeneity scales. In the same conditions,
the uncertainty of the transverse coefficient is tens of times
larger.
[48] Nevertheless, the common belief that large plumes
have ergodic behavior should be amended. The fluctuations
from realization to realization have an intricate nonlinear
dependence on the transverse dimension of the initial
plume, for both cross-section space averaged concentrations
and effective diffusion coefficients. For concentration and
longitudinal effective coefficient, the fluctuations decrease
when the transverse plume dimension is larger than ten
correlation lengths, but the traveltime to reach the monot-
onous decay toward a quasi-ergodic state is much larger
than for small sources. For the transverse effective coeffi-
cient the fluctuations increase with the dimension of the
initial plume at all times. It is expected that the timescale
which characterizes the ergodicity increases like the square
of the transverse dimension of the initial plume.
[49] The evolution of the ensemble averaged effective
coefficients also shows features not accounted for by the
existing theory. The slope of the displacements variance,
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY
11 of 17
W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
often used to estimate the effective coefficients, yields
nonphysical (negative) estimations of the transverse coef-
ficients at tens of advection timescales. Therefore the rate of
increase of the variance, which describes the spatial exten-
sion of the plume and is positive-definite, is the appropriate
definition for the effective coefficients. The condition of
vanishing center of mass coefficient in the large time limit
was also found to be not sufficient for the assessment of
the ergodic behavior. The failure of the usual approach to
define a transverse ‘‘ergodic coefficient’’ for extended
initial plumes indicates a dependence on initial conditions,
similar to the purely advective case, which deserves further
investigations.
Appendix A: Global Random Walk
A1. Algorithm
[50] The GRW algorithm is a generalization of particle
tracking method which increases the speed of the compu-
tations and considerably improves the accuracy of the
numerical simulations [Vamos¸ et al., 2003]. The solution
of a parabolic equation of form (1) is described using N
particles which move in a grid, undergoing advective
displacements and diffusive jumps according to the random
walk law. The concentration field at a given time t = kdt and
at a grid point (x
1
, x
2
)=(i
1
dx
1
, i
2
dx
2
) is given by
cx
1
; x
2
; t ð Þ¼
1
N D
1
D
2
X
s1
i
0
1
¼s1
X
s2
i
0
2
¼s2
ni
1
þ i
0
1
; i
2
þ i
0
2
; k
; ðA1Þ
where D
l
=2s
l
dx
l
, l = 1, 2, are the lengths of the symmetrical
intervals centered at x
l
and n(i
1
, i
2
, k) is the number of
particles which at time step k lie at the grid point (i
1
, i
2
).
[51] The one-dimensional GRW algorithm describes the
scattering of the n(i,k) particles from (x
i
, t
k
) by
nj; k ð Þ¼ dnj; j þ v
j
; k
þ dnj þ v
j
d; j; k
þ dnj þ v
j
þ d; j; k
;
where v
j
are discrete displacements in a given velocity field
and d describes the diffusive jumps. The quantities dn are
Bernoulli random variables and describe respectively, the
number of particles which remain at the same grid site after
an advective displacement, the number of particles jumping
to the left and those jumping to the right (with respect to the
advected position). The distribution of the particles at the
next time (k + 1)dt is given by
ni; k þ 1 ð Þ¼
X
j
dni; j; k ð Þ:
The average number of particles undergoing diffusive jumps
and the average number of particles remaining at the same
node after the displacement v
j
are given by the relations
dnj þ v
j
d; j; k
¼
1
2
r nj; k ð Þ;
dnj; j þ v
j
; k
¼ 1 r ð Þ nj; k ð Þ;
where 0 r 1. The diffusion coefficient D is related to
the grid steps by the relation
D ¼ r
ddx ð Þ
2
2dt
:
For two and three-dimensional cases, the same procedure is
repeated for all space directions.
[52] Because the total number of particles N contained in the
grid is conserved, the GRW algorithm is stable. The condition
r 1, ensures that there is no numerical diffusion. Vamos¸ et al.
[2003] showed that for Gaussian diffusion the numerical
solution converges as O(dx
2
)+ O(N
1/2
), i.e., for large
numbers of particles the convergence order is O(dx
2
), the same
as for the finite differences scheme. A comparison with a
particle tracking code (diffusion over ten time steps of N
particle starting at the center of a cubic grid) shows that while
in GRW algorithm there is practical no limitation, N > 10
9
particles becomes prohibitive for the particle tracking method.
[53] The ‘‘reduced fluctuations’’ GRW algorithm is de-
fined by
dnj þ v
j
d; j; k
¼
n=2 if n is even
n=2 ½ þ q if n is odd
8
<
:
;
where n = n(j, k) dn(j, j + v
j
, k), [n/2] is the integer part of n/2
and q is a variable taking the values 0 and 1 with probability
1/2. This algorithm is appropriate for large-scale problems,
for two reasons. Firstly, the diffusion front does not extend
beyond the limit concentration defined by one particle at a
grid point, keeping a physical significant shape (unlike in
finite differences where a pure diffusion front has a cubic
shape of side (2Dt)
1/2
). Second, the ‘‘reduced fluctuations’’
algorithm requires only a minimum number of calls of the
random number generator. Figure A1 illustrates the reduced
fluctuations GRW algorithm for r = 1, and Figure A2
presents the resulting concentration field computed by (A1).
[54] In the following we describe the computation of the
diagonal effective coefficients, according to the GRW
algorithm. The variance of particles displacements, s
ll
2
, l =
1,2, in dimensionless form, is given by
1
dx ð Þ
2
s
2
ll
k dt ð Þ¼
1
N
X
i1;i2
i
2
l
ni
1
; i
2
; k ð Þ
1
N
X
i1;i2
i
l
ni
1
; i
2
; k ð Þ
" #
2
:
ðA2Þ
Figure A1. Advective displacement and diffusive jumps
of 10
10
particles starting at (0, 0) for d = 1 and r = 1.
12 of 17
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Using (A2), the effective coefficients are computed as
D
eff
ll
k dt ð Þ¼ s
2
ll
= 2k dt ð Þ:
[55] Let us consider N
X
0
points uniformly distributed
inside the initial plume, N/N
X
0
particles at each initial point
and let n(i
1
, i
2
, k; i
01
, i
02
) be the distribution of particles at
the time step k given by the GRW procedure for a diffusion
process starting at (i
01
dx
1
, i
02
dx
2
). Writing the distribution
for the extended plume as
ni
1
; i
2
; k ð Þ¼
X
i01 ;i02
ni
1
; i
2
; k ; i
01
; i
02
ð Þ;
the averages from (A2) can be rewritten in the form
1
N
X
i1 ;i2
ani
1
; i
2
; k ð Þ¼
1
N
X
0
X
i01 ;i02
N
X
0
N
X
i1;i2
ani
1
; i
2
; k ; i
01
; i
02
ð Þ
!
;
ðA3Þ
where a stands for i
l
and i
l
2
respectively. It follows from
(A3) that the variance (A2) is an average over the
trajectories of the diffusion process starting at given initial
positions and over the distribution of the initial positions.
A2. GRW Parameters
[56] The large-scale computations reported by Schwarze
et al. [2001] and Dentz et al. [2002, 2003] were performed
with the particle tracking procedure. Although it was
possible to obtain estimations of the ensemble averaged
effective coefficients, the number of particles used in these
simulations, limited due to computational reasons at N
100, does not suffice to simulate accurate concentrations.
Even in small-scale one-dimensional problems, more than
one million particles should be used in particles methods to
reach the same precision as the finite difference scheme
[Vamos¸ et al., 2003]. Moreover, Suciu et al. [2004] showed
that for the large scale transport problem considered here, a
too small number of particles induces large numerical errors
in the simulation of the time behavior of the effective
coefficients. The statistical convergence of the simulations
for a given realization of the velocity is ensured only when
billions of particles are used. In the present numerical
investigations the number of particles was fixed at N =
10
10
so that all the simulations of transport realizations were
statistically convergent.
[57] The other GRW parameters used for the computation
of the ensembles of realizations presented in this paper are
the space steps dx
1
= dx
2
= dx = 0.1 m, the time step dt = 0.5
days and the amplitude of the diffusive jumps (see
section A1) d = 2. The accuracy of the numerical velocity
field is governed by the ratio of the log-hydraulic conduc-
tivity correlation length to the space step. In our case, l
y
/dx =
10 and fulfils the condition l
y
/dx 4, generally recom-
mended in literature [Ababou et al., 1989; Hassan et al.,
1998]. To reduce the ‘‘overshooting’’ errors in particles
methods one imposes that the mean displacement in a time
step does not surpass a given threshold [Roth and Hammel,
1996; Jankovic´ et al., 2003]. Suciu et al. [2004] showed that
for the same overshooting, the error in GRW simulations is
mainly influenced by the discretization of the velocity, as
described by the parameter Udt/dx. Our choice Udt/dx =5
means that, in average, the particles overpass 5 space steps,
but it also means that the smallest advective displacement
accounted for in the GRW procedure is dx/dt = U/5 = 0.2 m/d.
The tests for a crude estimation of the discretization errors,
for fixed dt = 0.5 day and increasing Udt/dx from 5 to 10
(dx/dt from 0.2 m/d to 0.1 m/d), show that the simulated
effective coefficients differ with less than 2% (N. Suciu et
al., Numerical investigations on ergodicity of solute trans-
port in heterogeneous aquifers, Internal Report ICG-IV
00204, Forschungszentrum Ju¨lich, 2004). A comparison
for the first 100 days with a GRW algorithm without
overshooting led to error estimations for the effective coef-
ficients in given realizations, calculated with the GRW
procedure and the parameters used in this paper, which were
one order of magnitude smaller than the upscaled coeffi-
cients (6) [Suciu et al., 2005].
Appendix B: Large-Scale Simulations
B1. Transport Problem
[58] The transport depends on both the heterogeneity of the
advection velocity field, described by the variance of the log
hydraulic conductivity s
y
2
, and the local dispersion coeffi-
cient D. To select the parameters to be used in the present
numerical investigations, several combinations of s
y
2
and D
were investigated in the case of point source (N. Suciu et al.,
Internal Report ICG-IV 00204, Forschungszentrum Ju¨lich,
2004). The diffusion fronts, defined by grid points containing
at least one particle, are compared in Figure B1. We note that
large variances of the log hydraulic conductivity (s
y
2
= 1)
yield non-Gaussian asymmetric plumes.
[59] It should be noted that the asymmetry of the diffu-
sion fronts can be due to a numerical artifact occurring in
two-dimensional simulations based on particles methods
and Kraichnan fields. In the absence of the third component
of the velocity, the probability of occurrence of very small
or null advection displacements is high enough to delay
some particles with respect to the plume center of mass.
This effect is compensated by diffusive displacements,
when the local dispersion is large enough [Suciu et al.,
2002]. At the limit of zero local dispersion, trapping zones
occur causing the fragmentation of the plume and the linear
increase of the effective coefficients [Dentz et al., 2003].
Figure A2. Concentration field, according to (A1), at t =
10dt, for dx
1
= dx
2
= 0.1 m and D
x
1
= D
x
2
= 1 m.
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY
13 of 17
W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
[60] The averages of the effective coefficients, over 256
realizations of the velocity, for fixed s
y
2
= 0.1 and three
different values of the local dispersion coefficients (D =
0.01 m
2
/d, D = 0.001 m
2
/d, and D = 0.0001 m
2
/d respec-
tively) are presented in Figures B2 and B3 as functions of
time and the corresponding Pe´clet numbers Pe´ = Ul
y
/D.
The comparisons from Figures B2 and B3 show that besides
the small asymmetry shown in Figure B1, the increase of
Pe´ considerably increases the time necessary to reach the
asymptotic coefficients D*
ll
(horizontal lines in Figures B2
and B3).
B2. Number of Periods Np, the Number of
Realizations R and the Correlation Shape
[61] The periodic fields generated with the Kraichnan
algorithm approximate Gaussian fields for Np !1. While
ensemble averages are well approximated for tens of periods
Np in the Kraichnan routine [Jaekel and Vereecken, 1997;
Schwarze et al., 2001], to approximate fluctuations much
larger Np are necessary [Eberhard, 2004]. To assess the
value of Np we compared the fluctuations of the cross-
section concentration and of the longitudinal effective
coefficient (defined as in section 4). We considered a point
source, a fixed number of realizations R = 1024, exponential
and Gaussian shape of the correlation of the log hydraulic
conductivity, with the same l
y
= 1 m and s
y
2
= 0.1. The
results presented in Figures B4 and B5 suggest that Np must
be at least as large as the total number of time steps in
simulations. (The fluctuations of the transverse coefficient,
not presented here, are already reliable for Np = 64.)
Therefore, for times up to t = 4000l
y
/U, we used Np =
6400 to approximate the behavior of the transport in
Gaussian velocity fields.
[62] Figure B6 presents the fluctuations of the longitudinal
effective coefficient for fixed Np = 6400 and increasing
number of realizations R, in the case of exponential correla-
tion of the log hydraulic conductivity. The increase of R from
256 to 1024 has little influence on the time behavior of the
fluctuations. We also found that R = 256 ensures the statistical
reliability for all the quantities investigated in this study.
B3. Velocity Statistics
[63] In most of numerical studies on the stochastic model
[Chin and Wang, 1992; Bellin et al., 1992; Salandin and
Figure B1. Diffusion fronts at t = 1000 days.
Figure B2. Longitudinal ensemble averaged effective
coefficient hD
11
eff
i for different Pe´clet numbers.
Figure B3. Transverse ensemble averaged effective coef-
ficient hD
22
eff
i for different Pe´clet numbers.
Figure B4. Concentration fluctuations at the plume center
of mass for R = 1024 and different Np.
14 of 17
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Fiorotto, 1998; Naff et al., 1998a; Hassan et al., 1998] the
Eulerian statistics of the numerically generated velocity
fields was estimated, under the implicit assumption of
statistical homogeneity, by space averages in given realiza-
tions followed by averages over realizations. This procedure
usually underestimates the true statistical parameters. To
show that, let us consider the components of the velocity
fluctuation u
l
(x)= V
l
(x) U
l
, l = 1, 2, supposed to be
statistically homogeneous variables, and the arithmetic
mean over P space points hu
l
(x
p
)i
P
=
1
P
P
P
p¼1
u
l
(x
p
). Because
the realizations of the velocity computed at different space
points belong to the same statistical ensemble, the order of
the space and ensemble averages can be interchanged and
due to the statistical homogeneity, we have the relations
hhu
l
i
P
i = hhu
l
ii
P
= hu
l
i = 0 and hhu
l
2
i
P
i = hhu
l
2
ii
P
= hu
l
2
i. The
average over realizations of the variance defined through
space averages,
hhu
2
l
i
P
hu
l
i
2
P
i¼hu
2
l
i hhu
l
i
2
P
i;
underestimates the true variance s
ul
2
= hu
l
2
i of the homoge-
neous variable u
l
with the term hhu
l
i
P
2
i, which is the
variance of the space mean hu
l
i
P
. The last quantity vanishes
only when the space mean equals the ensemble mean hu
l
i
P
=
hu
l
i = 0. The numerical fields have poor ergodic properties
and, as already noted by Bellin et al. [1992], they are
not strictly statistically homogeneous. Therefore the statisti-
cal properties of the numerical velocity fields should be
investigated through ensemble averages followed by space
averages. This procedure allows the estimation of the non-
homogeneity. For instance, the mean velocity hu
l
i is esti-
mated through space averages hhu
l
ii
P
, with the standard error
of the mean given by
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
hhu
l
i
2
i
P
hhu
l
ii
2
P
P 1
s
:
In our present numerical investigations we used such
standard errors to estimate the precision for the velocity
moments.
[64] Using 512 velocity values generated by Kraichnan
routine, at 512 different space points inside a square the side
of which was 10 l
y
, the velocity probability densities were
found to be very close to Gaussian homogeneous distribu-
tions. The first three moments of the longitudinal and
transverse velocity are presented in Table B1. The first
and third moment are close to zero. The longitudinal and
transverse variances (the second row in Table B1) are
respectively ’
3
8
s
y
2
and ’
1
8
s
y
2
, in agreement with the first
order asymptotic expansions of the Darcy and the continuity
equations [Dagan, 1984].
[65] The velocity correlations
r
ll
x ð Þ ¼ hhu
l
x
01
; x
02
ð Þu
l
x
01
þ x; x
02
ð Þii
P
; l ¼ 1; 2;
were computed as averages over 512 realizations of the
velocity and over P = 11011 points (x
01
, x
02
) (all the grid
points in a band of l
y
100l
y
, which corresponds to the
largest initial plume in the present simulations). The inte-
grals of the correlation functions r
ll
give the numerical
estimation of the second terms in (4), which describe the
contribution of the velocity fluctuations to the upscaled
diffusion coefficients,
J
ll
t ðÞ¼
Z
t
0
r
ll
Ut
0
ð Þdt
0
¼
1
U
Z
Ut
0
r
ll
x ðÞdx:
The numerical integrations between 0 and t = 5000l
y
/U of the
correlation functions are presented in Figure B7. The values
of J
ll
are close to the theoretical values (6), J*
11
= D*
11
D =
0.1 m
2
/d and J*
11
= D
22
* D =0m
2
/d. Because the correlations
computed by ensemble averages hu
l
(x
01
, x
02
)u
l
(x
01
+ x, x
02
)i
differ from point to point, i.e., the random field is not
Table B1. First Three Moments of the Longitudinal and
Transverse Velocity Components
l =1 l =2
hhu
l
ii
P
0.00214 ± 0.00033 0.00051 ± 0.00021
hh(u
l
hu
l
i)
2
ii
P
0.03801 ± 0.00010 0.01268 ± 0.00004
hh(u
l
hu
l
i)
3
ii
P
0.00014 ± 0.00003 0.00000 ± 0.00001
Figure B5. Fluctuations of the longitudinal effective
coefficient for R = 1024 and different Np.
Figure B6. Fluctuations of the longitudinal effective
coefficient for Np = 6400 and different R.
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY
15 of 17
W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
strictly homogeneous, at large distances the upper estimate of
r
ll
(i.e., space mean plus standard error) is mainly positive
and the lower estimate mainly negative. As a direct
consequence, the confidence intervals of J
ll
grow linearly
with x (see thin lines in Figure B7). Since the standard errors
decrease as (P 1)
1/2
, reliable estimations of the upscaled
effective dispersion coefficients require averaging over large
space domains. This remark is valid not only in the case of
numerical fields generated by the Kraichnan algorithm but
also for all large-scale numerical simulations, where the
accumulation of the numerical errors can result in large
uncertainty of the numerical estimations.
[66] Acknowledgments. The research reported in this paper was
supported in part by the Deutsche Forschungsgemeinschaft grant SU 415/
1-1, awarded to the first author. Part of the work is a contribution to the
‘‘Interdisciplinary Programme for the Prevention of the Major Risk Phe-
nomena at National Level’’of the Romanian Academy. The computations
on the Jump supercomputer at Research center Ju¨lich were done in
cooperation with the John von Neumann Institute for Computing in the
frame of project JICG41. The authors are grateful to U. Jaekel, M. Dentz,
H. Schwarze, J. P. Eberhard, and A. Englert for their comments on an earlier
version of the manuscript and to three anonymous reviewers for construc-
tive remarks. They also acknowledge the technical assistance during the
computations offered by J. Heidbu¨chel and P. Ustohal.
References
Ababou, R., D. McLaughlin, L. W. Gelhar, and A. F. B. Tompson (1989),
Numerical simulation of three-dimensional saturated flow in randomly
heterogeneous porous media, Transp. Porous Media, 4, 549 – 565.
Attinger, S., M. Dentz, H. Kinzelbach, and W. Kinzelbach (1999), Temporal
behavior of a solute cloud in a chemically heterogeneous porous medium,
J. Fluid Mech., 386, 77 – 104.
Avellaneda, M., and M. Majda (1989), Stieltjes integral representation and
effective diffusivity bounds for turbulent diffusion, Phys. Rev. Lett.,
62(7), 753 – 755.
Avellaneda, M., and M. Majda (1992), Superdiffusion in nearly stratified
flows, J. Stat. Phys., 69(3/4), 689 – 729.
Bellin, A., P. Salandin, and A. Rinaldo (1992), Simulation of dispersion in
heterogeneous porous formations: Statistics, first-order theories, conver-
gence of computations, Water Resour. Res., 28(9), 2211 – 2227.
Berkowitz, B. (2001), Dispersion in heterogeneous geological formations:
Preface, Transp. Porous Media, 42, 1–2.
Berkowitz, B., and H. Scher (2001), Probabilistic approaches to transport in
heterogeneous media, Transp. Porous Media, 42, 241 – 263.
Chin, D. A., and T. Wang (1992), An investigation of the validity of first-
order stochastic dispersion theories in isotropic porous media, Water
Resour. Res, 28(6), 1531 – 1542.
Clincy, M., and H. Kinzelbach (2001), Stratified disordered media: Exact
solutions for transport parameters and their self-averaging properties,
J. Phys. A Math. Gen., 34, 7142 – 7152.
Cortis, A., H. Scher, and B. Berkowitz (2004), Numerical simulation
of non-Fickian transport in geological formations with multiple-scale
heterogeneity, Water Resour. Res. , 40, W04209, doi:10.1029/
2003WR002750.
Cushman, J. H., and M. Moroni (2001), Statistical mechanics with three-
dimensional particle tracking velocimetry experiments in the study of
anomalous dispersion. I. Theory, Phys. Fluids, 13(1), 75 – 80.
Dagan, G. (1984), Solute transport in heterogeneous porous formations,
J. Fluid Mech., 145, 151 – 177.
Dagan, G. (1987), Theory of solute transport by groundwater, Water
Resour. Res., 19, 183 – 215.
Dagan, G. (1990), Transport in heterogeneous porous formations: Spatial
moments, ergodicity, and effective dispersion, Water Resour. Res., 26,
1281 – 1290.
Dagan, G. (2004), Comment on ‘‘Exact averaging of stochastic equations
for transport in random velocity field,’’ Transport in Porous Media, 50,
223– 241, 2003, and ‘‘Probability density functions for solute transport
in random field,’’ Transport in Porous Media, 50, 243 – 266, 2003,
Transp. Porous Media, 55, 113 – 116.
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 – 1607.
Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach (2000a),
Temporal behavior of a solute cloud in a heterogeneous porous medium:
1. Point-like injection, Water Resour. Res., 36, 3591 – 3604.
Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach (2000b),
Temporal behavior of a solute cloud in a heterogeneous porous medium
2. Spatially extended injection, Water Resour. Res., 36, 3605 – 3614.
Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach (2002), Tempor-
al behavior of a solute cloud in a heterogeneous porous medium: 3.
Numerical simulations, Water Resour. Res., 38(7), 1118, doi:10.1029/
2001WR000436.
Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach (2003), Numer-
ical studies of the transport behavior of a passive solute in a two-dimen-
sional incompressible random flow field, Phys. Rev. E, 67, 046306.
Eberhard, J. (2004), Approximations for transport parameters and self-
averaging properties for point-like injections in heterogeneous media,
J. Phys. A Math. Gen., 37, 2549 – 2571.
Ferna`ndez-Garcia, D., T. H. Illangasekare, and H. Rajaram (2005a), Differ-
ences in the scale dependence of dispersivity and retardation factors
estimated from forced-gradient and uniform flow tracer tests in three-
dimensional physically and chemically heterogeneous porous media,
Water Resour. Res., 41, W03012, doi:10.1029/2004WR003125.
Ferna`ndez-Garcia, D., H. Rajaram, and T. H. Illangasekare (2005b), Assess-
ment of the predictive capabilities of the stochastic theories in a three-
dimensional laboratory test aquifer: Effective hydraulic conductivity and
temporal moments of breakthrough curves, Water Resour. Res., 41,
W04002, doi:10.1029/2004WR003523.
Fiori, A. (1996), Finite Peclet extensions of Dagan’s solutions to transport in
anisotropic heterogeneous formations, Water Resour. Res., 32, 193 – 198.
Fiori, A. (1998), On the influence of pore-scale dispersion in nonergodic
transport in heterogeneous formations, Transp. Porous Media, 30, 57 – 73.
Fiori, A., and G. Dagan (2000), Concentration fluctuations in aquifer trans-
port: A rigorous first-order solution and applications, J. Contam. Hydrol.,
45, 139 – 163.
Fiori, A., I. Jankovic´, and G. Dagan (2003), Flow and transport in highly
heterogeneous formations: 2. Semianalytical results for isotropic media,
Water Resour. Res., 39(9), 1269, doi:10.1029/2002WR001719.
Hassan, A. H., J. H. Cushman, and J. W. Delleur (1998), A Monte Carlo
assessment of Eulerian flow and transport perturbation models, Water
Resour. Res., 34(5), 1143 – 1163.
Jaekel, U., and H. Vereecken (1997), Renormalization group analysis of
macrodispersion in a directed random flow, Water Resour. Res., 33,
2287 – 2299.
Jankovic´, I., A. Fiori, and G. Dagan (2003), Flow and transport in highly
heterogeneous formations: 3. Numerical simulations and comparison
with theoretical results, Water Resour. Res., 39(9), 1270, doi:10.1029/
2002WR001721.
Kabala, Z. J., and G. Sposito (1994), Statistical moments of reactive solute
concentration in a heterogeneous aquifer, Water Resour. Res., 30, 759–
768.
Kapoor, C., and L. W. Gelhar (1994), Transport in three-dimensionally
heterogeneous aquifers: 1. Dynamics of concentration fluctuations, Water
Resour. Res., 30(6), 1775 – 1788.
Figure B7. Progress of the time integral of the velocity
correlation functions along with mean values and con-
fidence intervals.
16 of 17
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
Kapoor, C., and P. K. Kitanidis (1998), Concentration fluctuation and dilu-
tion in aquifers, Water Resour. Res., 34, 1181 – 1193.
Kesten, H., and G. C. Papanicolaou (1979), A limit theorem for turbulent
diffusion, Commun. Math. Phys., 65, 97–128.
Kitanidis, P. K. (1988), Prediction by the method of moments of transport in
a heterogeneous formation, J. Hydrol., 102, 453 – 473.
Kraichnan, R. H. (1970), Diffusion by a random velocity field, Phys.
Fluids, 13(1), 22 – 31.
Kurbanmuradov, O., K. Sabelfeld,O. F. Smidts, and H. Vereecken (2003),
A Lagrangian stochastic model for the transport in statistically homoge-
neous porous media, Monte Carlo Methods Appl., 9(4), 341 – 366.
Labolle, E. M., and G. Fogg (2001), Role of molecular diffusion in con-
taminant migration and recovery in an alluvial aquifer system, Transp.
Porous Media, 42, 155 – 179.
Matheron, G., and G. de Marsily (1980), Is transport in porous media al-
ways diffusive?, Water Resour. Res., 16, 901 – 917.
Moltyaner, G. L., M. H. Klukas, C. A. Willis, and R. W. D. Killey (1993),
Numerical simulations of Twin Lake natural-gradient tracer tests: A
comparison of methods, Water Resour. Res., 29(10), 3433 – 3452.
Naff, R. L., D. F. Haley, and E. A. Sudicky (1998a), High-resolution Monte
Carlo simulation of flow and conservative transport in heterogeneous
porous media: 1. Methodology and flow results, Water Resour. Res.,
34(4), 663 – 677.
Naff, R. L., D. F. Haley, and E. A. Sudicky (1998b), High-resolution Monte
Carlo simulation of flow and conservative transport in heterogeneous
porous media: 2. Transport results, Water Resour. Res., 34(4), 679 – 697.
Pannone, M., and P. K. Kitanidis (1999), Large time behavior of concen-
tration variance and dilution in heterogeneous formations, Water Resour.
Res., 35(3), 623 – 634.
Rajaram, H., and L. W. Gelhar (1993a), Plume scale-dependent dispersion
in heterogeneous aquifers: 1. Lagrangian analysis in a stratified aquifer,
Water Resour. Res., 29, 3249 – 3260.
Rajaram, H., and L. W. Gelhar (1993b), Plume scale-dependent dispersion
in heterogeneous aquifers: 2. Eulerian analysis and three-dimensional
aquifers, Water Resour. Res., 29, 3261 – 3276.
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.
Salandin, P., and V. Fiorotto (1998), Solute transport in highly heteroge-
neous aquifers, Water Resour. Res., 34, 949 – 961.
Schwarze, H., U. Jaekel, and H. Vereecken (2001), Estimation of macro-
dispersivity by different approximation methods for flow and transport in
randomly heterogeneous media, Transp. Porous Media, 43, 265 – 287.
Smith, L., and F. W. Schwartz (1980), Mass transport: 1. A stochastic
analysis of macroscopic dispersion, Water Resour. Res., 16, 303 – 313.
Sposito, G. (1997), Ergodicity and the ‘‘scale effect’’, Adv. Water Resour.,
20, 309 – 316.
Sposito, G. (2001), Topological groundwater hydrodynamics, Adv. Water
Resour., 24, 793 – 801.
Sposito, G., and G. Dagan (1994), Predicting solute plume evolution
in heterogeneous porous formations, Water Resour. Res., 30(2), 585 –
589.
Sposito, G., W. A. Jury, and V. K. Gupta (1986), Fundamental problems in
the stochastic convection-dispersion model of solute transport in aquifers
and field soils, Water Resour. Res., 22, 77–88.
Suciu, N., C. Vamos¸, H. Vereecken, and J. Vanderborght (2002), Numerical
modeling of solute transport in heterogeneous aquifers and self-aver-
aging, in Proceedings of the first Workshop on Mathematical Modeling
of Environmental Problems, edited by G. Marinovschi and I. Stelian,
pp. 111– 140, Rom. Acad. Publ. House, Bucharest.
Suciu, N., C. Vamos¸, J. Vanderborght, H. Hardelauf, and H. Vereecken
(2004), Numerical modeling of large scale transport of contaminant so-
lutes using the global random walk algorithm, Monte Carlo Methods
Appl., 10(2), 153 – 177.
Suciu, N., C. Vamos¸, P. Knabner, and U. Ru¨de (2005), Biased global ran-
dom walk, a cellular automaton for diffusion, in Simulationstechnique,
18th Symposium in Erlangen, September 2005, edited by F. Hu¨lsemann,
M. Kowarschik, and U. Ru¨ de, pp. 562 – 567, SCS Publ. House, Erlangen,
Germany.
Sun, N.-Z. (1996), Mathematical Modeling in Groundwater Pollution,
Springer, New York.
Tatarinova, E. B., P. A. Kalugin, and A. V. Sokol (1991), What is the
propagation rate of the passive component in turbulent flow limited
by?, Europhys. Lett., 14(8), 773 – 777.
Taylor, G. I. (1921), Diffusion by continuous movements, Proc. London
Math. Soc., 2(20), 196 – 212.
Tompson, A. F. B., and L. W. Gelhar (1990), Numerical simulation of
solute transport in three-dimensional, randomly heterogeneous porous
media, Water Resour. Res., 26(10), 2541 – 2562.
Trefry, M. G., F. P. Ruan, and D. McLaughlin (2003), Numerical simula-
tions of preasymptotic transport in heterogeneous porous media: Depar-
tures from the Gaussian limit, Water Resour. Res. , 39(3), 1063,
doi:10.1029/2001WR001101.
Vamos¸, C., N. Suciu, and H. Vereecken (2003), Generalized random walk
algorithm for the numerical modeling of complex diffusion processes,
J. Comput. Phys., 186(2), 527 – 544.
Vanderborght, J., and H. Vereecken (2002), Estimation of local scale dis-
persion from local breackthrough curves during a tracer test in a hetero-
genaous aquifer: The Lagrangian approach, J. Contam. Hydol., 54, 141 –
171.
Winter, C. L., C. M. Newman, and S. P. Neuman (1984), A perturbation
expansion for diffusion in random velocity field, SIAM J. Appl. Math.,
44(2), 411 – 424.
Zhang, Y.-K., and J. Lin (1998), Numerical simulations of transport of non-
ergodic solute plumes in heterogeneous aquifers, Stochastic Hydrol. Hy-
draul., 12(2), 117 – 139.
Zhang, Y.-K., and B.-M. Seo (2004), Stochastic analyses and Monte Carlo
simulations of nonergodic solute transport in three-dimensional hetero-
geneous statistically anisotropic aquifers, Water Resour. Res., 40,
W05103, doi:10.1029/2003WR002871.
Zhang, Y.-K., and D. Zhang (1997), Time-dependent dispersion of none-
rgodic plumes in two-dimensional heterogeneous aquifers, J. Hydrol.
Eng., 2(2), 91 – 94.
Zhang, Y.-K., D. Zhang, and J. Lin (1996), Nonergodic solute transport in
three-dimensional heterogeneous isotropic aquifers, Water Resour. Res.,
32(9), 2955 – 2963.
H. Hardelauf, J. Vanderborght, and H. Vereecken, Institute of Agro-
sphere (IGC-IV), Research Center Ju¨lich, D-52425 Ju¨lich, Germany.
N. Suciu, Institute of Applied Mathematics, Friedrich-Alexander
University of Erlangen-Nuremberg, Martensstrasse 3, D-91058 Erlangen,
Germany. (suciu@am.uni-erlangen.de)
C. Vamos¸, ‘‘T. Popoviciu’’ Institute of Numerical Analysis, Romanian
Academy, P.O. Box 68-1, 400320 Cluj-Napoca, Romania.
W04409 SUCIU ET AL.: NUMERICAL INVESTIGATIONS ON ERGODICITY
17 of 17
W04409
19447973, 2006, 4, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2005WR004546 by Cochrane Romania, Wiley Online Library on [28/08/2025]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
[1] Ababou, R., D. McLaughlin, L. W. Gelhar, and A. F. B. Tompson (1989), Numerical simulation of three‐dimensional saturated flow in randomly heterogeneous porous media, Transp. Porous Media, 4, 549–565.
CrossRef (DOI)
[2] Attinger, S., M. Dentz, H. Kinzelbach, and W. Kinzelbach (1999), Temporal behavior of a solute cloud in a chemically heterogeneous porous medium, J. Fluid Mech., 386, 77–104.
CrossRef (DOI)
[3] Avellaneda, M., and M. Majda (1989), Stieltjes integral representation and effective diffusivity bounds for turbulent diffusion, Phys. Rev. Lett., 62(7), 753–755.
CrossRef (DOI)
[4] Avellaneda, M., and M. Majda (1992), Superdiffusion in nearly stratified flows, J. Stat. Phys., 69(3/4), 689–729.
CrossRef (DOI)
[5] Bellin, A., P. Salandin, and A. Rinaldo (1992), Simulation of dispersion in heterogeneous porous formations: Statistics, first‐order theories, convergence of computations, Water Resour. Res., 28(9), 2211–2227.
CrossRef (DOI)
[6] Berkowitz, B. (2001), Dispersion in heterogeneous geological formations: Preface, Transp. Porous Media, 42, 1–2.
CrossRef (DOI)
[7] Berkowitz, B., and H. Scher (2001), Probabilistic approaches to transport in heterogeneous media, Transp. Porous Media, 42, 241–263.
CrossRef (DOI)
[8] Chin, D. A., and T. Wang (1992), An investigation of the validity of first‐order stochastic dispersion theories in isotropic porous media, Water Resour. Res, 28(6), 1531–1542.
CrossRef (DOI)
[9] Clincy, M., and H. Kinzelbach (2001), Stratified disordered media: Exact solutions for transport parameters and their self‐averaging properties, J. Phys. A Math. Gen., 34, 7142–7152.
CrossRef (DOI)
[10] Cortis, A., H. Scher, and B. Berkowitz (2004), Numerical simulation of non‐Fickian transport in geological formations with multiple‐scale heterogeneity, Water Resour. Res., 40, W04209.
CrossRef (DOI)
[11] Cushman, J. H., and M. Moroni (2001), Statistical mechanics with three‐dimensional particle tracking velocimetry experiments in the study of anomalous dispersion. I. Theory, Phys. Fluids, 13(1), 75–80.
CrossRef (DOI)
[12] Dagan, G. (1984), Solute transport in heterogeneous porous formations, J. Fluid Mech., 145, 151–177.
CrossRef (DOI)
[13] Dagan, G. (1987), Theory of solute transport by groundwater, Water Resour. Res., 19, 183–215.
CrossRef (DOI)
[14] Dagan, G. (1990), Transport in heterogeneous porous formations: Spatial moments, ergodicity, and effective dispersion, Water Resour. Res., 26, 1281–1290.
CrossRef (DOI)
[15] Dagan, G. (2004), Comment on “Exact averaging of stochastic equations for transport in random velocity field,” Transport in Porous Media, 50, 223–241, 2003, and “Probability density functions for solute transport in random field,” Transport in Porous Media, 50, 243–266, 2003, Transp. Porous Media, 55, 113–116.
CrossRef (DOI)
[16] 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–1607.
CrossRef (DOI)
[17] Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach (2000a), Temporal behavior of a solute cloud in a heterogeneous porous medium: 1. Point‐like injection, Water Resour. Res., 36, 3591–3604.
CrossRef (DOI)
[18] Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach (2000b), Temporal behavior of a solute cloud in a heterogeneous porous medium 2. Spatially extended injection, Water Resour. Res., 36, 3605–3614.
CrossRef (DOI)
[19] Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach (2002), Temporal behavior of a solute cloud in a heterogeneous porous medium: 3. Numerical simulations, Water Resour. Res., 38(7), 1118.
CrossRef (DOI)
[20] Dentz, M., H. Kinzelbach, S. Attinger, and W. Kinzelbach (2003), Numerical studies of the transport behavior of a passive solute in a two‐dimensional incompressible random flow field, Phys. Rev. E, 67, 046306.
CrossRef (DOI)
[21] Eberhard, J. (2004), Approximations for transport parameters and self‐averaging properties for point‐like injections in heterogeneous media, J. Phys. A Math. Gen., 37, 2549–2571.
CrossRef (DOI)
[22] Fernàndez‐Garcia, D., T. H. Illangasekare, and H. Rajaram (2005a), Differences in the scale dependence of dispersivity and retardation factors estimated from forced‐gradient and uniform flow tracer tests in three‐dimensional physically and chemically heterogeneous porous media, Water Resour. Res., 41, W03012, doi:10.1029/2004WR003125.
CrossRef (DOI)
[23] Fernàndez‐Garcia, D., H. Rajaram, and T. H. Illangasekare (2005b), Assessment of the predictive capabilities of the stochastic theories in a three‐dimensional laboratory test aquifer: Effective hydraulic conductivity and temporal moments of breakthrough curves, Water Resour. Res., 41, W04002.
CrossRef (DOI)
[24] Fiori, A. (1996), Finite Peclet extensions of Dagan’s solutions to transport in anisotropic heterogeneous formations, Water Resour. Res., 32, 193–198.
CrossRef (DOI)
[25] Fiori, A. (1998), On the influence of pore‐scale dispersion in nonergodic transport in heterogeneous formations, Transp. Porous Media, 30, 57–73.
CrossRef (DOI)
[26] 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 (DOI)
[27] Fiori, A., I. Janković, and G. Dagan (2003), Flow and transport in highly heterogeneous formations: 2. Semianalytical results for isotropic media, Water Resour. Res., 39(9), 1269.
CrossRef (DOI)
[28] Hassan, A. H., J. H. Cushman, and J. W. Delleur (1998), A Monte Carlo assessment of Eulerian flow and transport perturbation models, Water Resour. Res., 34(5), 1143–1163.
CrossRef (DOI)
[29] Jaekel, U., and H. Vereecken (1997), Renormalization group analysis of macrodispersion in a directed random flow, Water Resour. Res., 33, 2287–2299.
CrossRef (DOI)
[30] Janković, I., A. Fiori, and G. Dagan (2003), Flow and transport in highly heterogeneous formations: 3. Numerical simulations and comparison with theoretical results, Water Resour. Res., 39(9), 1270.
CrossRef (DOI)
[31] Kabala, Z. J., and G. Sposito (1994), Statistical moments of reactive solute concentration in a heterogeneous aquifer, Water Resour. Res., 30, 759–768.
CrossRef (DOI)
[32] Kapoor, C., and L. W. Gelhar (1994), Transport in three‐dimensionally heterogeneous aquifers: 1. Dynamics of concentration fluctuations, Water Resour. Res., 30(6), 1775–1788.
CrossRef (DOI)
[33] Kapoor, C., and P. K. Kitanidis (1998), Concentration fluctuation and dilution in aquifers, Water Resour. Res., 34, 1181–1193.
CrossRef (DOI)
[34] Kesten, H., and G. C. Papanicolaou (1979), A limit theorem for turbulent diffusion, Commun. Math. Phys., 65, 97–128.
CrossRef (DOI)
[35] Kitanidis, P. K. (1988), Prediction by the method of moments of transport in a heterogeneous formation, J. Hydrol., 102, 453–473.
CrossRef (DOI)
[36] Kraichnan, R. H. (1970), Diffusion by a random velocity field, Phys. Fluids, 13(1), 22–31.
CrossRef (DOI)
[37] Kurbanmuradov, O., K. Sabelfeld, O. F. Smidts, and H. Vereecken (2003), A Lagrangian stochastic model for the transport in statistically homogeneous porous media, Monte Carlo Methods Appl., 9(4), 341–366.
CrossRef (DOI)
[38] Labolle, E. M., and G. Fogg (2001), Role of molecular diffusion in contaminant migration and recovery in an alluvial aquifer system, Transp. Porous Media, 42, 155–179.
CrossRef (DOI)
[39] Matheron, G., and G. de Marsily (1980), Is transport in porous media always diffusive? Water Resour. Res., 16, 901–917.
CrossRef (DOI)
[40] Moltyaner, G. L., M. H. Klukas, C. A. Willis, and R. W. D. Killey (1993), Numerical simulations of Twin Lake natural‐gradient tracer tests: A comparison of methods, Water Resour. Res., 29(10), 3433–3452.
CrossRef (DOI)
[41] Naff, R. L., D. F. Haley, and E. A. Sudicky (1998a), High‐resolution Monte Carlo simulation of flow and conservative transport in heterogeneous porous media: 1. Methodology and flow results, Water Resour. Res., 34(4), 663–677.
CrossRef (DOI)
[42] Naff, R. L., D. F. Haley, and E. A. Sudicky (1998b), High‐resolution Monte Carlo simulation of flow and conservative transport in heterogeneous porous media: 2. Transport results, Water Resour. Res., 34(4), 679–697.
CrossRef (DOI)
[43] Pannone, M., and P. K. Kitanidis (1999), Large time behavior of concentration variance and dilution in heterogeneous formations, Water Resour. Res., 35(3), 623–634.
CrossRef (DOI)
[44] Rajaram, H., and L. W. Gelhar (1993a), Plume scale‐dependent dispersion in heterogeneous aquifers: 1. Lagrangian analysis in a stratified aquifer, Water Resour. Res., 29, 3249–3260.
CrossRef (DOI)
[45] Rajaram, H., and L. W. Gelhar (1993b), Plume scale‐dependent dispersion in heterogeneous aquifers: 2. Eulerian analysis and three‐dimensional aquifers, Water Resour. Res., 29, 3261–3276.
CrossRef (DOI)
[46] 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 (DOI)
[47] Salandin, P., and V. Fiorotto (1998), Solute transport in highly heterogeneous aquifers, Water Resour. Res., 34, 949–961.
CrossRef (DOI)
[48] Schwarze, H., U. Jaekel, and H. Vereecken (2001), Estimation of macrodispersivity by different approximation methods for flow and transport in randomly heterogeneous media, Transp. Porous Media, 43, 265–287.
CrossRef (DOI)
[49] Smith, L., and F. W. Schwartz (1980), Mass transport: 1. A stochastic analysis of macroscopic dispersion, Water Resour. Res., 16, 303–313.
CrossRef (DOI)
[50] Sposito, G. (1997), Ergodicity and the “scale effect”, Adv. Water Resour., 20, 309–316.
CrossRef (DOI)
[51] Sposito, G. (2001), Topological groundwater hydrodynamics, Adv. Water Resour., 24, 793–801.
CrossRef (DOI)
[52] Sposito, G., and G. Dagan (1994), Predicting solute plume evolution in heterogeneous porous formations, Water Resour. Res., 30(2), 585–589.
CrossRef (DOI)
[53] Sposito, G., W. A. Jury, and V. K. Gupta (1986), Fundamental problems in the stochastic convection‐dispersion model of solute transport in aquifers and field soils, Water Resour. Res., 22, 77–88.
CrossRef (DOI)
[54] Suciu, N., C. Vamoş, H. Vereecken, and J. Vanderborght (2002), Numerical modeling of solute transport in heterogeneous aquifers and self‐averaging, in Proceedings of the first Workshop on Mathematical Modeling of Environmental Problems, edited by G. Marinovschi, and I. Stelian, pp. 111–140, Rom. Acad. Publ. House, Bucharest.
[55] Suciu, N., C. Vamoş, J. Vanderborght, H. Hardelauf, and H. Vereecken (2004), Numerical modeling of large scale transport of contaminant solutes using the global random walk algorithm, Monte Carlo Methods Appl., 10(2), 153–177.
CrossRef (DOI)
[56] Suciu, N., C. Vamoş, P. Knabner, and U. Rüde (2005), Biased global random walk, a cellular automaton for diffusion, in Simulationstechnique, 18th Symposium in Erlangen, September 2005, edited by F. Hülsemann, M. Kowarschik, and U. Rüde, pp. 562–567, SCS Publ. House, Erlangen, Germany.
[57] Sun, N.‐Z. (1996), Mathematical Modeling in Groundwater Pollution, Springer, New York.
CrossRef (DOI)
[58] Tatarinova, E. B., P. A. Kalugin, and A. V. Sokol (1991), What is the propagation rate of the passive component in turbulent flow limited by? Europhys. Lett., 14(8), 773–777.
CrossRef (DOI)
[59] Taylor, G. I. (1921), Diffusion by continuous movements, Proc. London Math. Soc., 2(20), 196–212.
CrossRef (DOI)
[60] Tompson, A. F. B., and L. W. Gelhar (1990), Numerical simulation of solute transport in three‐dimensional, randomly heterogeneous porous media, Water Resour. Res., 26(10), 2541–2562.
CrossRef (DOI)
[61] Trefry, M. G., F. P. Ruan, and D. McLaughlin (2003), Numerical simulations of preasymptotic transport in heterogeneous porous media: Departures from the Gaussian limit, Water Resour. Res., 39(3), 1063.
CrossRef (DOI)
[62] Vamoş, C., N. Suciu, and H. Vereecken (2003), Generalized random walk algorithm for the numerical modeling of complex diffusion processes, J. Comput. Phys., 186(2), 527–544.
CrossRef (DOI)
[63] Vanderborght, J., and H. Vereecken (2002), Estimation of local scale dispersion from local breackthrough curves during a tracer test in a heterogenaous aquifer: The Lagrangian approach, J. Contam. Hydol., 54, 141–171.
CrossRef (DOI)
[64] Winter, C. L., C. M. Newman, and S. P. Neuman (1984), A perturbation expansion for diffusion in random velocity field, SIAM J. Appl. Math., 44(2), 411–424.
CrossRef (DOI)
[65] Zhang, Y.‐K., and J. Lin (1998), Numerical simulations of transport of non‐ergodic solute plumes in heterogeneous aquifers, Stochastic Hydrol. Hydraul., 12(2), 117–139.
CrossRef (DOI)
[66] Zhang, Y.‐K., and B.‐M. Seo (2004), Stochastic analyses and Monte Carlo simulations of nonergodic solute transport in three‐dimensional heterogeneous statistically anisotropic aquifers, Water Resour. Res., 40, W05103.
CrossRef (DOI)
[67] Zhang, Y.‐K., and D. Zhang (1997), Time‐dependent dispersion of nonergodic plumes in two‐dimensional heterogeneous aquifers, J. Hydrol. Eng., 2(2), 91–94.
CrossRef (DOI)
[68] Zhang, Y.‐K., D. Zhang, and J. Lin (1996), Nonergodic solute transport in three‐dimensional heterogeneous isotropic aquifers, Water Resour. Res., 32(9), 2955–2963.
CrossRef (DOI)