format: text

text:
# Physical Pictures of Transport in Heterogeneous Media: Advection-Dispersion, Random Walk and Fractional Derivative Formulations

## 1. Introduction

Quantification of chemical transport mediated by flow fields in strongly heterogeneous geological environments has received an inordinate amount of attention over the last three decades, and a vast literature dealing with the subject has developed (see, e.g., the recent reviews in Dagan and Neuman ). Existing modeling approaches are generally based on various deterministic and stochastic forms of the advection-dispersion equation (ADE); the former include conditioning the domain of interest by known heterogeneity structures, while the latter include Monte Carlo, perturbation and spectral analyses. A major feature of transport, particularly in more heterogeneous domains, is the appearance of “scale-dependent dispersion” \[e.g., Gelhar et al., 1992\]. Contrary to the fundamental assumptions underlying use of the classical ADE (which assumes a constant flow field and dispersion coefficients), the very nature of the dispersive transport seems to change as a function of time or distance traveled by the contaminant. Such scale-dependent behavior, also sometimes referred to as “pre-asymptotic”, “anomalous” or “non-Gaussian”, is what we shall refer to as “non-Fickian” transport.

Efforts to quantify non-Fickian transport have focused on more general stochastic ADE’s with, e.g., spatially varying velocity fields. Stochastic analyses have provided substantial insight into the dispersion process. They have been shown, through application to well-documented field experiments, to provide predictions of the temporal variation of the first and second order moments of tracer plumes in geological formations characterized by relatively small degrees of heterogeneity (e.g., the Cape Cod site \[Garabedian et al., 1991\]). Other variations based on the classical ADE have also received attention; these include “patch” solutions which include an empirical time- or space-dependent dispersivity, and mobile-immobile and multirate diffusion type models \[e.g., Haggerty and Gorelick, 1995; Harvey and Gorelick, 2000\]. However, the vast majority of these models assume, either explicitly or implicitly, an underlying Fickian transport behavior at some scale \[e.g., Sposito et al., 1986; Rubin, 1997\]. Also, many of these approaches are based on perturbation theory, and they are therefore limited to porous media in which the variance of the log hydraulic conductivity is small.

Other non-local formulations that do not invoke a Fickian transport assumption have been hypothesized and/or developed from various mathematical formalisms \[e.g., Zhang, 1992; Glimm et al., 1993; Neuman, 1993; Deng et al., 1993; Cushman et al., 1994; Dagan, 1997\]. These formalisms, in general, are founded on a fundamental separation between advective and dispersive mechanisms; they yield solutions (for the concentration) that result in definition of a dispersion tensor that is usually formulated in Fourier-Laplace space, whose inversion is difficult to treat and/or apply.

Practical application of these models, to quantify the full evolution of a migrating contaminant plume, has not yet been achieved. In fact, the overwhelming emphasis of these various studies has been limited to moment characterizations of tracer plume migration, and/or to determination of the “macrodispersion” parameter. The complete solutions are not analytically tractable, and their practical utility remains largely undemonstrated.

The difficulty in capturing the complexities of tracer plume migration patterns suggests that local, small-scale heterogeneities cannot be neglected. Evidently, these unresolvable heterogeneities contribute significantly to the occurrence of non-Fickian transport. The apparent existence of hydraulic conductivity fields with coherence lengths that vary over many scales suggests that temporal, as well as spatial issues must be considered in any mathematical formulation. Coupled to this problem is the lack of clarity of how best to use field observations to reduce the inevitable uncertainties of the model. Frequently, the latter issue involves the interplay between ensemble averaging (probabilistic approaches) and spatial scales of resolution of non-stationary geological features.

In this paper, we re-evaluate the basic conceptual picture of tracer migration in heterogeneous media. We derive the general form of the governing equations for conservative chemical transport in heterogeneous geological formations, for single realizations and for ensemble averages of the domain. We emphasize quantification of non-Fickian transport behavior, and show that a general form of the ensemble-averaged transport equation is a continuous time random walk (CTRW). In this framework, we show that non-Fickian transport results from the inapplicability of the central limit theorem to capture the distribution of particle transitions (detailed in the next section). Fractional derivative formulations of the transport equations, both temporal and spatial, are seen to emerge from another set of conditions, and are therefore special cases of the CTRW. We then focus on quantifying transport in non-stationary media, and discuss how best to deal with the coupled problem of integrating ensemble averaging with information on non-stationarity at various scales of resolution.

## 2. Governing Transport Equations for Heterogeneous Media

### 2.1. Physical Framework of the Transport Equations

Contaminants disperse as they migrate within the flow field of the geological maze we call an aquifer. At the outset one must choose an underlying physical model of this process. Two possible models include Taylor dispersion and multiple transitions. Taylor dispersion is based on molecular diffusion of particles in a flowing fluid (e.g., in a pipe) and is governed by an ADE, to be discussed below. An identical formulation can be obtained by considering particle movement in a random network and applying the central limit theorem. The extensive use of the ADE in the hydrology literature is based essentially on the generic concept of Taylor dispersion and works well for relatively homogeneous systems. The particles are assumed to be transported by the average flowing fluid in the medium while the “diffusion” is the dispersion due to local medium irregularities. Larger scale effects (e.g., permeability changes) are treated as perturbations of this model in conventional stochastic treatments.

The prime interest in this work is in highly heterogeneous systems; in these systems contaminant motion can be envisioned as a migrating cloud of particles, each of which executes a series of steps or transitions between changes in velocity v. The spatial extent of these transitions depends on the criterion used to define changes in v. The classical approach is to consider the system divided into representative elementary volumes (REVs) and determine an average v and dispersion D in each REV. In our approach we dispense with the REV idea, because averages can be unreliable in a system of very wide fluctuations about the mean value. The change of concentration $`\mathrm{\Delta }C`$ at each position in a time increment $`\mathrm{\Delta }t`$ is $`\mathrm{\Delta }t\times `$(the net particle flux). The effective volume contributing the net particle flux in $`\mathrm{\Delta }t`$ can vary considerably at different positions in the system. Thus the length scale over which $`\mathrm{\Delta }C`$ varies slowly in space can change considerably over the system. If one fixes a sampling volume at each position, it is important to retain the full distribution (not an average) of the transition times (determined with a physical model) of flux contributing to $`\mathrm{\Delta }C`$. If this distribution is retained, then in our approach one can still use the limit of a spatial continuum (as shown below).

The distribution of transition times, $`\psi (t)`$, can be determined in principle from an analysis of the streamtubes of the flow field and contains the subtle features that can produce non-Fickian behavior. The physical features necessary for non-Fickian transport are the existence of a wide range of transition times (causing large differences in the flow paths of migrating particles) and sufficient encounter with statistically rare, but rate-limiting slow transitions (e.g., low velocity regions) \[Berkowitz and Scher, 1995\]. These general ideas will be developed schematically in the next sections.

### 2.2. Single Realization Transport Equation

For our point of departure we need a transport equation framework that can enumerate all these possible paths and encompass the motion from continuous to discrete over a range of spatial and temporal scales, for any given realization of the domain. An excellent candidate is the “Master Equation” \[Oppenheim et al., 1977; Shlesinger, 1996\]

$$\frac{C(𝐬,t)}{t}=\underset{𝐬^{}}{}w(𝐬^{},𝐬)C(𝐬,t)+\underset{𝐬^{}}{}w(𝐬,𝐬^{})C(𝐬^{},t)$$
(1)
for $`C(𝐬,t)`$, the particle concentration at point s and time $`t`$, where $`w(𝐬,𝐬^{})`$ is the transition rate from s to s (the dimension of $`\mathrm{\Sigma }_𝐬w`$ is reciprocal time). The transition rates describe the effects of the velocity field on the particle motion; the determination of $`w(𝐬,𝐬^{})`$ involves a detailed knowledge of the system. We assume the average effective range of $`w(𝐬,𝐬^{})`$ is a finite distance. The Master Equation has been applied in the context of electron hopping in random systems \[e.g., Klafter and Silbey, 1980a\], and is discussed widely in the physics and chemistry literature.

The transport equation in (1) does not separate the effects of the varying velocity field into an advective and dispersive part of the motion; this separation is an approximation based on the assumption of relatively homogeneous regions in which $`C(𝐬,t)`$ will be slowly varying over a finite length scale (the range of transition rates),

$`C(𝐬^{},t)C(𝐬,t)+(𝐬^{}𝐬)C(𝐬,t)+\frac{1}{2}(𝐬^{}𝐬)(𝐬^{}𝐬):C(𝐬,t)`$ (2)
(with the dyadic symbol : denoting a tensor product). Substituting (2) into (1) leads to a continuum description (i.e., local diffusion in a pressure field $`\pi `$(s)) and a partial differential equation (pde), for a single realization of the domain:

$`{\displaystyle \frac{C(𝐬,t)}{t}}`$ $`=`$ $`{\displaystyle \underset{𝐬^{}}{}}(w(𝐬,𝐬^{})w(𝐬^{},𝐬))C(𝐬,t)+{\displaystyle \underset{𝐬^{}}{}}w(𝐬,𝐬^{})(𝐬^{}𝐬)C(𝐬,t)`$ (3)
$`+{\displaystyle \underset{𝐬^{}}{}}w(𝐬,𝐬^{})\frac{1}{2}(𝐬^{}𝐬)(𝐬^{}𝐬):C(𝐬,t).`$
We note that (3) is close to the form of an ADE with the exception of the term proportional to $`C(𝐬,t)`$. This term is present due to the asymmetry of the transition rates (due to the bias of the pressure field) and/or the non-stationary medium (due to the explicit position dependence of the rates – cf. (4)). It makes a contribution to the final form of the pde for diffusion in a force field. If the system is stationary this term vanishes (as we show below) and thus reduces to the form of an ADE. One can already observe in (3) generalized velocity and dispersion coefficients (in terms of $`w(𝐬,𝐬^{})`$); however we have not yet separated out the effects of the flow field and determined transport coefficients. In order to fully determine the final pde and separate the advection and diffusion contributions, we must specify the $`w(𝐬,𝐬^{})`$ in terms of $`\pi (𝐬)`$, the pressure field.

A general form for a non-stationary medium is

$$w(𝐬,𝐬^{})W(|𝐬^{}𝐬|;𝐬^{})\mathrm{\Omega }(\pi (𝐬^{})\pi (𝐬))$$
(4)
where the asymmetry in the rates is due to $`\pi (𝐬^{})\pi (𝐬)`$, the pressure difference at $`𝐬^{}`$ and $`𝐬`$, and the explicit dependence of the overall rate W on location ($`\mathrm{\Omega }`$ is a function of the pressure difference only). We specify the $`\mathrm{\Omega }`$-function, so that (4) is written as

$`w(𝐬,𝐬^{})`$ $``$ $`W(|𝐬^{}𝐬|;𝐬^{})\mathrm{\Omega }(\pi (𝐬^{})\pi (𝐬))`$ (5)
$``$ $`F(|𝐬^{}𝐬|;𝐬^{})\left[\lambda +\frac{1}{2}(\pi (𝐬^{})\pi (𝐬))\right]`$
where in (5) non-linear terms in the pressure difference have been neglected (i.e., terms proportional to $`(\pi )^2`$) and a contribution to the transition rates is retained even for vanishing pressure difference. The significance of the latter step can be seen by realizing that $`F(\pi (𝐬^{})\pi (𝐬))`$ is a simple advection contribution (with a permeability proportional to $`F`$) and the term $`F\lambda `$ is proportional to a local diffusion contribution to the rates. The $`\lambda `$ term retains the scattering effects of the medium (i.e., the transfers between “streamtubes”) even in the limit of very small local pressure differences. It is also closely associated with the effect of “local” dispersion.

We now also assume $`F(|𝐬^{}𝐬|;𝐬^{})`$ will be slowly varying over some finite length scale. We expand in a Taylor series to second order in $`𝐬^{}𝐬`$,

$`F(|𝐬^{}𝐬|;𝐬^{})`$ $``$ $`F(|𝐬^{}𝐬|;𝐬)+(𝐬^{}𝐬)F`$ (6)
$`+\frac{1}{2}(𝐬^{}𝐬)(𝐬^{}𝐬):F.`$
In (6), the gradient operates on the second argument, s. Combining (5) and (6), and substituting into the first term on the right side of (3), we have

$`w(𝐬,𝐬^{})w(𝐬^{},𝐬)`$ $``$ $`F(|𝐬^{}𝐬|;𝐬^{})\left[\lambda +\frac{1}{2}(\pi (𝐬^{})\pi (𝐬))\right]F(|𝐬^{}𝐬|;𝐬)\left[\lambda \frac{1}{2}(\pi (𝐬^{})\pi (𝐬))\right]`$ (7)
$``$ $`F(|𝐬^{}𝐬|;𝐬)(\pi (𝐬^{})\pi (𝐬))`$
$`+[(𝐬^{}𝐬)F+\frac{1}{2}(𝐬^{}𝐬)(𝐬^{}𝐬):F]\times [\lambda +\frac{1}{2}(\pi (𝐬^{})\pi (𝐬))].`$
Now using a similar expansion for the pressure difference, we have

$$\pi (𝐬^{})\pi (𝐬)(𝐬^{}𝐬)\pi (𝐬)+\frac{1}{2}(𝐬^{}𝐬)(𝐬^{}𝐬):\pi (𝐬).$$
(8)
Substituting (8) into (7) and using

$$\underset{𝐬^{}}{}F(|𝐬^{}𝐬|;𝐬)(𝐬^{}𝐬)=0$$
(9)
because $`F`$ is an even function of the vector difference, we obtain for the expression in (7), summed over s,

$$\left(\frac{𝐃(𝐬)}{\lambda }\pi (𝐬)\right)+𝐃(𝐬)$$
(10)
where the dispersion tensor is defined as

$$𝐃(𝐬)\frac{1}{2}\underset{𝐬^{}}{}F(|𝐬^{}𝐬|;𝐬)(𝐬^{}𝐬)(𝐬^{}𝐬)\lambda .$$
(11)
We insert (10) into (3) and use (4)–(6), (8), (9) (cf. Appendix A) to obtain

$$\frac{C(𝐬,t)}{t}=\left[\frac{𝐃(𝐬)}{\lambda }\pi (𝐬)C(𝐬,t)+(𝐃(𝐬)C(𝐬,t))\right].$$
(12)
The form of (12) is a continuity equation – the time derivative of the concentration is equal to the divergence of the total concentration flux, the sum of the diffusive concentration flux and the advective concentration flux – with an effective permeability of

$$𝐤(𝐬)\frac{𝐃(𝐬)}{\lambda }.$$
(13)
Equation (12) is a generalization to a non-stationary medium of the well-known Smoluchowski equation \[Chandrasekhar, 1943\] which is the basis for describing diffusion in a force field. In our case the force field is $`\pi (𝐬)`$. In the case of electron transfer in a potential field the $`\lambda `$ in (13) can be shown to be $`\kappa `$T (where T is the temperature and $`\kappa `$ is Boltzmann’s constant) and the relation in (13) is the Einstein relation between mobility and diffusion. We use a convention that a product between a tensor T and a vector V is TV yielding a vector. In our case, the vector $`𝐯(𝐬)=𝐤(𝐬)\pi (𝐬)`$ is the velocity field and for an incompressible fluid $`𝐯(𝐬)=0`$. The only term remaining in (12) proportional to $`C(𝐬,t)`$ is $`𝐃(𝐬)C(𝐬,t)`$. The final form for the pde for an incompressible fluid is

$$\frac{C(𝐬,t)}{t}=𝐯(𝐬)C(𝐬,t)+(𝐃(𝐬)C(𝐬,t)).$$
(14)
Equation (14) is a generalization of the ADE. While many simplifications of the ADE are based on (14) with $`𝐃(𝐬)=𝐃`$ (i.e., a constant), the usual (“general”) form of the ADE includes a s-dependent D in (14) but with the second term replaced by $`(𝐃(𝐬)C(𝐬,t))`$. Thus (14) differs from this usual form of the ADE by the addition of two terms: $`𝐃(𝐬)C(𝐬,t)`$ and $`𝐃(𝐬)C`$. The form of (14) is the same as postulated by Kinzelbach , based on the Ito process.

The difference in the general form of the ADE can be traced to starting the derivation with the pressure field $`\pi (𝐬)`$ and not with $`\pi (𝐬)`$, i.e., the expansion (8) is treated on the same basis as the other expansions (2) and (6). Hence, starting with the Master equation (1) and using a general expression for the transfer rates we obtain, for a specific heterogeneous medium, in a continuum limit (slowly varying $`C(𝐬,t)`$ and $`w(𝐬,𝐬^{}))`$ the generalized equation for diffusion in a force field (Smoluchowski) which for irrotational flow is a generalized ADE. We assert that for a non-stationary medium, i.e., s-dependent v and D, (14) should be the starting point for numerical calculations. The main numerical differences between this equation and the usual ADE (with D(s)) should arise in “boundary” regions of more spatially varying D(s). The importance of accounting for D(s) has been demonstrated by, e.g., Labolle et al. .

We will show that the “standard” ADE emerges as the continuum limit of the ensemble averaged Master equation (the term proportional to $`C(𝐬,t)`$ vanishes for stationary transition rates). In general, the continuum limit presents difficulties in regions of increased heterogeneity, such as tightly interspersed permeability layers. The concentration $`C(𝐬,t)`$ will not necessarily vary slowly on the same length scale throughout the system. The point average of v and D can be very sensitive to small changes in the local volume used to determine the average. Conversely, if one fixes the volume to a practical pixel size (e.g., 10 m3) the use of a local average v and D in each volume can be quite limited, i.e., the spreading effects of unresolved residual heterogeneities are suppressed \[e.g., Dagan, 1997\]. We will return to this issue in a broader context in section 4. It essentially involves the degrees of uncertainty and its associated spatial scales. We start, at first, with an ensemble average of the entire medium and discuss the role of this approach in the broader context.

### 2.3. Ensemble Average Transport Equation

We resume our examination of the Master Equation approach, i.e., before assuming any continuum limit. The ensemble average of (1) can be shown \[Klafter and Silbey, 1980b\] to be of the form

$`{\displaystyle \frac{P(𝐬,t)}{t}}={\displaystyle \underset{𝐬^{}}{}}{\displaystyle _0^t}\varphi (𝐬^{}𝐬,tt^{})P(𝐬,t^{})𝑑t^{}+{\displaystyle \underset{𝐬^{}}{}}{\displaystyle _0^t}\varphi (𝐬𝐬^{},tt^{})P(𝐬^{},t^{})𝑑t^{}`$ (15)
where $`P(𝐬,t)`$ is the normalized concentration, and $`\varphi (𝐬,t)`$ is defined below in (20). The form of (15) is a “Generalized Master Equation” (GME) which, in contrast to (1), is non-local in time and the transition rates are stationary (i.e., depend only on the difference ss) and time-dependent. This equation describes a semi-Markovian process (Markovian in space, but not in time), which accounts for the time correlations (or “memory”) in particle transitions.

It is straightforward to show \[Kenkre et al., 1973; Shlesinger, 1974\], using the Laplace transform, that the GME is completely equivalent to a continuous time random walk (CTRW)

$$R(𝐬,t)=\underset{𝐬^{}}{}_0^t\psi (𝐬𝐬^{},tt^{})R(𝐬^{},t^{})𝑑t^{}$$
(16)
where $`R(𝐬,t)`$ is the probability per time for a walker to just arrive at site s at time $`t`$, and $`\psi (𝐬,t)`$ is the probability rate for a displacement s with a difference of arrival times of $`t`$. The initial condition for $`R(𝐬,t)`$ is $`\delta _{𝐬,0}\delta (t0^+)`$, which can be appended to (16). The correspondence between (15) and (16) is

$$P(𝐬,t)=_0^t\mathrm{\Psi }(tt^{})R(𝐬,t^{})𝑑t^{}$$
(17)
where

$$\mathrm{\Psi }(t)=1_0^t\psi (t^{})𝑑t^{}$$
(18)
is the probability for a walker to remain on a site,

$$\psi (t)\underset{𝐬}{}\psi (𝐬,t)$$
(19)
and

$$\stackrel{~}{\varphi }(𝐬,u)=\frac{u\stackrel{~}{\psi }(𝐬,u)}{1\stackrel{~}{\psi }(u)}$$
(20)
where the Laplace transform ($``$) of a function $`f(t)`$ is denoted by $`\stackrel{~}{f}(u)`$.

Equations (16)–(19) are in the form of a convolution in space and time and can therefore be solved by use of Fourier and Laplace transforms \[Scher and Lax, 1973\]. The general solution is

$$𝒫(𝐤,u)=\frac{1\stackrel{~}{\psi }(u)}{u}\frac{1}{1\mathrm{\Lambda }(𝐤,u)}$$
(21)
where $`𝒫(𝐤,u)`$, $`\mathrm{\Lambda }(𝐤,u)`$ are the Fourier transforms ($``$) of $`\stackrel{~}{P}(𝐬,u)`$, $`\stackrel{~}{\psi }(𝐬,u)`$, respectively.

The CTRW accounts naturally for the cumulative effects of a sequence of transitions. The challenge is to map the important aspects of the particle motion in the medium onto a $`\psi (𝐬,t)`$. The identification of $`\psi (𝐬,t)`$ lies at the heart of the CTRW formulation. The CTRW approach allows a determination of the evolution of the particle distribution (plume), $`P(𝐬,t)`$, for a general $`\psi (𝐬,t)`$; there is no a priori need to consider only the moments of $`P(𝐬,t)`$. As we discuss below, a $`\psi (𝐬,t)`$ with a power law (30) for large time leads to the description of anomalous transport (e.g., non-Fickian plumes). Once $`\psi (𝐬,t)`$ is defined one needs to calculate $`\mathrm{\Lambda }(𝐤,u)`$ and then determine the propagator $`P(𝐬,t)`$ by inverting the Fourier and Laplace transform of (21). The latter can be quite challenging.

As shown previously the separation between advection and dispersion occurs in the continuum (diffusion) limit. In an ensemble averaged system this limit leads to an ADE \[Berkowitz and Scher, 2001\]. For clarity and convenience, we reproduce the argument here. The first step is to make a series expansion of $`P(𝐬,t)`$ similar to (2); inserting this into (15) yields

$`{\displaystyle \frac{P(𝐬,t)}{t}}={\displaystyle \underset{𝐬^{}}{}}{\displaystyle _0^t}dt^{}[\varphi (𝐬𝐬^{},tt^{})(𝐬^{}𝐬)P(𝐬,t^{})+\varphi (𝐬𝐬^{},tt^{})\frac{1}{2}(𝐬^{}𝐬)(𝐬^{}𝐬):P(𝐬,t^{})].`$ (22)
We write (22) in a more compact form

$`{\displaystyle \frac{P(𝐬,t)}{t}}={\displaystyle _0^t}dt^{}[𝐯_\psi (tt^{})P(𝐬,t^{})+𝚽_\psi (tt^{}):P(𝐬,t^{})]`$ (23)
$$𝐯_\psi (t)\underset{𝐬}{}\varphi (𝐬,t)𝐬$$
(24)
$$𝚽_\psi (t)\underset{𝐬}{}\varphi (𝐬,t)\frac{1}{2}\mathrm{𝐬𝐬}$$
(25)
Note the sum (over $`𝐬^{}`$) in (22) is independent of $`𝐬`$ in a stationary system; hence we shift the summation variable to obtain (24)-(25). This particular formulation is convenient because, in (23), we can define terms that are familiar in the context of traditional modeling: the “effective velocity” $`𝐯_\psi `$ and the “dispersion tensor” $`𝚽_\psi `$. Note, however, that both of these terms are time-dependent, and most significantly, depend fundamentally on $`\psi (𝐬,t)`$. This equation has the form of an ADE generalized to non-local time responses as a result of the ensemble average.

The next step is a crucial one in distinguishing between normal and anomalous transport. If $`\psi (𝐬,t)`$ has both a finite first and second moment in $`t`$ the transport is normal and one can expand $`\stackrel{~}{\psi }(𝐬,u)`$ as \[Scher and Montroll, 1975\]

$`\stackrel{~}{\psi }(𝐬,u)p_1(𝐬)p_2(𝐬)u+p_3(𝐬)u^2+\mathrm{}`$
$`\text{and}\stackrel{~}{\psi }(u)={\displaystyle \underset{𝐬}{}}\stackrel{~}{\psi }(𝐬,u)1\overline{t}u+du^2+\mathrm{}`$ (26)
with $`_𝐬p_1(𝐬)=1`$, the normalization of $`\psi (𝐬,t)`$, and $`_𝐬p_2(𝐬)\overline{t}`$ and $`_𝐬p_3(𝐬)d`$, the first and second temporal moments of $`\psi (t)`$, respectively. Note that small $`u`$ corresponds to large time in Laplace space. The functions $`p_i(𝐬)`$ are asymmetric due to the bias in the velocity field; $`p_1(𝐬)`$ is the probability to make a step of displacement s. One now inserts (2.3. Ensemble Average Transport Equation) into (20) and expands in a power series of $`u`$. The leading term is independent of $`u`$, which we retain. The correction to this leading term is proportional to $`u`$ and is small. Substituting this expression into the Laplace transform of (23)-(25), which is (53)-(55) (cf. below), and taking the inverse Laplace transform of the result, yields the ADE

$$\frac{P(𝐬,t)}{t}=𝐯P(𝐬,t)+𝐃:P(𝐬,t)$$
(27)
where the effective velocity v is equal to the first spatial moment of $`p_1(𝐬),\overline{𝐬}`$, the mean displacement for a single transition, divided by the mean transition time $`\overline{t}`$, and the dispersion tensor $`𝐃D_{ij}`$ is the second spatial moment divided by $`\overline{t}`$, which can be written as

$$𝐯=\underset{𝐬}{}p_1(𝐬)𝐬/\overline{t}\overline{𝐬}/\overline{t}$$
(28)
$$D_{ij}=v\frac{1}{2}_𝐬p_1(𝐬)s_is_j/\overline{s}$$
(29)
where $`v=|𝐯|`$ and $`\overline{s}=|\overline{𝐬}|`$. If we retain the term proportional to $`u`$ when inserting (2.3. Ensemble Average Transport Equation) into (20), we obtain terms with both spatial and temporal derivatives of $`P(𝐬,t)`$.

Thus, our underlying physical picture of advective-driven dispersion reduces to the familiar ADE when one can assume smooth spatial variation of $`P(𝐬,t)`$ and finite first and second temporal moments of $`\psi (𝐬,t)`$.

### 2.4. Non-Fickian Dispersion

When the $`\psi (𝐬,t)`$ has a power law (algebraic tail) dependence on time at large $`t`$, i.e.,

$$\psi (𝐬,t)t^{1\beta }$$
(30)
the first and second temporal moments do not exist for $`0<\beta <1`$, while the second temporal moment does not exist for $`1<\beta <2`$. The dependence of $`\psi (𝐬,t)`$ in (30) is a manifestation of a wide distribution of event times as encountered in highly heterogeneous media. The relation between the power law behavior (30) and non-Fickian (anomalous) transport has been well documented \[e.g., Scher and Montroll, 1975; Berkowitz and Scher, 2001\]. We sketch the key points of that relationship: The form of $`\psi (𝐬,t)`$ at large time determines the time dependence of the mean position $`\overline{\mathrm{}}(t)`$ and standard deviation $`\overline{\sigma }(t)`$ of $`P(𝐬,t)`$. In the presence of a pressure gradient (or “bias”), and for (30), it can be shown \[Scher and Montroll, 1975; Shlesinger, 1974\] for $`0<\beta <1`$ that

$$\overline{\mathrm{}}(t)t^\beta $$
(31)
$$\overline{\sigma }(t)t^\beta $$
(32)
while for $`1<\beta <2`$

$$\overline{\mathrm{}}(t)t$$
(33)
$$\overline{\sigma }(t)t^{(3\beta )/2}.$$
(34)
Moreover, it can be shown that Fickian-like transport arises when $`\beta >2`$ \[e.g., Margolin and Berkowitz, 2000\].

The unusual time dependence of $`\overline{\mathrm{}}(t)`$ and $`\overline{\sigma }(t)`$ in (31)-(34), resulting from the infinite temporal moments of $`\psi (𝐬,t)`$ (i.e., the conditions of the central limit theorem are not fulfilled), is the hallmark of the non-Fickian propagation of $`P(𝐬,t)`$. This behavior is in sharp contrast to Fickian models where, $`\overline{\mathrm{}}(t)t`$ and $`\overline{\sigma }(t)t^{1/2}`$ (as an outcome of the central limit theorem) and the position of the peak of the distribution coincides with $`\overline{\mathrm{}}(t)`$. Note that in Fickian transport, $`\overline{\mathrm{}}(t)/\overline{\sigma }(t)t^{1/2}`$; an important distinguishing feature of anomalous transport is that $`\overline{\mathrm{}}(t)/\overline{\sigma }(t)constant`$ for $`0<\beta <1`$, and $`\overline{\mathrm{}}(t)/\overline{\sigma }(t)t^{(\beta 1)/2}`$ for $`1<\beta <2`$. The relative shapes of the anomalous transport curves, and the rate of advance of the peak, vary strongly as a function of $`\beta `$. Thus the parameter $`\beta `$ effectively quantifies the contaminant dispersion; this parameter is discussed in detail by, e.g., Margolin and Berkowitz  and Berkowitz and Scher . Hence, the crucial consideration for the appearance of non-Fickian dispersion in a specified scale of a heterogeneous medium are the physical criteria for the power law (30) and its (time) range of applicability. Non-Fickian transport that displays these characteristics has been documented in several analyses of numerical simulations, and laboratory and field data \[Berkowitz and Scher, 1998; Hatano and Hatano, 1998; Berkowitz et al., 2000; Kosakowski et al., 2001\].

The large time regime of $`\psi (𝐬,t)`$ corresponds to the small $`u`$ regime for its Laplace transform and the expansion in $`u`$ (for (30)) is quite different from (2.3. Ensemble Average Transport Equation) \[Shlesinger, 1974\], i.e.,

$$\stackrel{~}{\psi }(𝐬,u)p_1^{}(𝐬)p_2^{}(𝐬)u^\beta +\mathrm{}$$
(35)
for $`u0`$ for $`0<\beta <1`$. Inserting (35) into (20), parallel to the development following (2.3. Ensemble Average Transport Equation), yields a transport equation from (22) which remains non-local in time and is not the ADE. Our development \[Berkowitz and Scher, 1995\] of non-Fickian transport has been based directly on (15). In other words, solutions for the full evolution of a tracer plume, as well as for breakthrough curves (i.e., spatial and temporal distributions of tracer) can be derived directly from (15) \[e.g., Scher and Montroll, 1975; Berkowitz and Scher, 1997, 1998\]. A (fractional) pde form of the transport equation, derived from (22) and holding only for the power law dependence (30), i.e., a special case of CTRW, is exhibited in section 3.2. We observe also that the $`u0`$ expansion of $`\stackrel{~}{\psi }(𝐬,u)`$ for $`1<\beta <2`$ is similar to (2.3. Ensemble Average Transport Equation), but with the $`u^2`$ term replaced by one proportional to $`u^\beta `$. In this case the correction to the $`u`$-independent term $`p_1(𝐬)/\overline{t}`$ used in (28), (29) is proportional to $`u^{\beta 1}`$ and can be significant (especially for $`\beta 1`$).

Finally, we note that the general CTRW formalism (i.e., not restricted to (30)) can be used to model a large number of physical processes. For example, $`\psi (𝐬,t)`$ has been defined for multiple trapping \[e.g., Scher et al., 1991; Hatano and Hatano, 1998\] and as such can be used for multiple-rate models \[Haggerty and Gorelick, 1995\] and to quantify dispersion in stratified formations \[Matheron and de Marsily, 1980\]. Zumofen et al.  have used the CTRW explicitly to model the latter.

## 3. Fractional Differential Equations

There is growing interest in the development and application of fractional differential formulations of transport equations. In particular, fractional differential equations of the diffusion, diffusion–advection, and Fokker–Planck type have been considered in stochastic modeling in physics \[e.g., Hilfer, 2000; Metzler and Klafter, 2000\]. Here we consider fractional derivative equations (FDE) for transport and show how they are special cases of the CTRW equations developed in the previous section. We emphasize that FDE are not different models from the CTRW; rather, they are seen to emerge as asymptotic limit cases of the CTRW theory.

A word of caution: referring to a transport equation as “fractional” can be with respect to the occurrence of fractional order differentiation in time or space, or both. Moreover, a number of definitions for fractional operators exist. Here, we concentrate on two possibilities: the Riemann–Liouville fractional time derivative $`{}_{0}{}^{}D_{t}^{\beta }`$ (for which we will employ the more suggestive notation $`^\beta /t^\beta `$), and the Riesz spatial derivative $`^\mu `$ \[Oldham and Spanier, 1974; Samko et al., 1993\].

The development of FDE in both the time and space variables necessitates a more general starting equation than (22), which depends on the validity of the expansion of $`P(𝐬,t)`$ similar to (2). We return to the general solution (21). In what follows, in order to obtain FDE’s, we need the product form $`p(𝐬)\psi (t)`$ for the $`\psi (𝐬,t)`$ probability density function, which assumes that the transition length and time are statistically independent quantities. Furthermore we need the asymptotic form (30) of $`\psi (t)`$ and/or $`p(𝐬)`$ (cf. below). The indicated power-law decay for $`0<\beta <1`$ causes the divergence of $`\overline{t}`$, the mean transition time (cf. section 2.4). Corresponding to (30) the Laplace transform of $`\psi (t)`$ is

$$\stackrel{~}{\psi }(u)1(uc_t)^\beta $$
(36)
which is (35) summed over $`𝐬`$, where $`c_t`$ is a dimensional constant determined by the physical model. Along the same line we consider the power–law form $`p(𝐬)c_{s}^{}{}_{}{}^{\mu }/|𝐬|^{1+\mu },0<\mu <2`$ for the transition length, where $`c_s`$, is analogous to $`c_t`$, a dimensional constant. Similar to $`\psi (t)`$, the first and second or second (spatial) moment(s) of $`p(𝐬)`$ are infinite for, respectively, $`0<\mu <1`$ and $`1<\mu <2`$. The border case for $`\mu =2`$ is the Gaussian law $`p(𝐬)(4\pi c_{s}^{}{}_{}{}^{2})^1\mathrm{exp}(𝐬^2/(4c_s))`$. For any symmetric Lévy stable law $`p(𝐬)`$, the asymptotic form of the Fourier transform of $`p(𝐬)`$ is given by

$$\stackrel{~}{p}(𝐤)1c_{s}^{}{}_{}{}^{\mu }|𝐤|^\mu 0<\mu 2.$$
(37)
### 3.1 Time-FDE

We concentrate on the case $`0<\beta <1`$ and $`\mu =2`$, for which the spatial moments are finite, but the temporal moments are infinite. We consider first the case with no spatial bias, $`\overline{\mathrm{}}(t)=0`$ (i.e., no advective transport). Insertion of (36) and the low wavenumber expression $`\stackrel{~}{p}(𝐤)1c_{s}^{}{}_{}{}^{2}𝐤^2`$ into (21) leads to

$$\stackrel{~}{𝒫}(𝐤,u)=\frac{1}{u+K_\beta u^{1\beta }𝐤^2}$$
(38)
(dropping the cross term $`(uc_t)^\beta c_{s}^{}{}_{}{}^{2}𝐤^2`$) where the anomalous diffusion constant is defined as $`K_\beta c_{s}^{}{}_{}{}^{2}/c_{t}^{}{}_{}{}^{\beta }`$. The FDE is determined by multiplying (38) by the denominator of the right side and rearranging to yield

$$u𝒫(𝐤,u)1=K_\beta 𝐤^2u(u^\beta 𝒫(𝐤,u)),$$
(39)
where the dimension of the generalized diffusion constant is $`[K_\beta ]=\mathrm{cm}^2\mathrm{sec}^\beta `$. While the two terms on the left correspond to $`P(𝐬,t)/t`$ in $`(𝐬,t)`$ space, with the initial condition $`P(𝐬,0)=\delta (𝐬)`$ (on both sides of (39) the property $`\{dF(t)/dt\}=u\stackrel{~}{F}(u)F(0)`$ is utilized), the factor $`u^\beta `$ on the right poses the problem of finding the corresponding Laplace inversion. One of the definitive responses goes back to Riemann and Liouville who extended the Cauchy multiple integral, in order to define the fractional integral,

$$\frac{^\beta }{t^\beta }P(𝐬,t)\frac{1}{\mathrm{\Gamma }(\beta )}_0^t𝑑t^{}\frac{P(𝐬,t^{})}{(tt^{})^{1\beta }}$$
(40)
which possesses the important property

$$\left\{\frac{^\beta }{t^\beta }P(𝐬,t)\right\}=u^\beta \stackrel{~}{P}(𝐬,u).$$
(41)
The definition (40) explicitly includes the initial value at time $`t=0`$. Note that for a negative index, $`^\beta /t^\beta `$, the Riemann–Liouville operator denotes fractional integration whereas for a positive index, $`^\beta /t^\beta `$, we have fractional differentiation. In our case fractional differentiation is established as the succession of fractional integration and standard differentiation:

$$\frac{^{1\beta }}{t^{1\beta }}P(𝐬,t)=\frac{}{t}\frac{^\beta }{t^\beta }P(𝐬,t).$$
(42)
With these definitions, we can now invert (39), and obtain the fractional diffusion equation

$$\frac{P}{t}=K_\beta \frac{^{1\beta }}{t^{1\beta }}^2P(𝐬,t).$$
(43)
In the limit $`\beta 1`$ (43) reduces to the standard Brownian version.

The generalization to a fractional ADE for anomalous transport $`(0<\beta <1)`$, which includes a spatial bias (advective transport), follows the same procedure as above \[Compte, 1997; Compte et al., 1997; Compte and Càceres, 1998; Metzler et al., 1998; Metzler and Compte, 2000\],

$$\frac{}{t}P(𝐬,t)=\frac{^{1\beta }}{t^{1\beta }}\left(𝐯_\beta +K_\beta ^2\right)P(𝐬,t)$$
(44)
where $`𝐯_\beta `$ is the “generalized drift velocity”. Note that (43) and (44) involve fractional differentiation in time on the spatial derivative terms of the equations. These equations can be rewritten so they do not involve mixed derivatives, if desired \[Metzler and Klafter, 2000\]. We stress that the form of (43) and (44) relies on using (36), and that (44) is valid only for $`0<\beta <1`$; it is modified significantly for $`1<\beta <2`$. We have thus shown that the probability density $`P(𝐬,t)`$ described by the time–fractional ADE (44), is equivalent to the large–time limit of the CTRW with a bias, with the asymptotic form of $`\psi (t)`$ given by (30) (or $`\stackrel{~}{\psi }(u)`$ given by (36)). For a specific class of $`\psi (t)`$ (which also fulfills the asymptotic form (36)), the equivalence between CTRW and FDE can be shown over the entire range of $`t`$ \[Hilfer and Anton, 1995\].

### 3.2. Space-FDE: Lévy Flights

We now consider the opposite case of a transition time distribution with an existing first moment, $`\beta >1`$, $`\stackrel{~}{\psi }(u)1uc_t`$, and a transition length distribution $`p(𝐬)`$ with a diverging second moment ($`0<\mu <2`$) ($`\{p(𝐬)\}`$ in (37)). This case can be shown to be a Markovian process (in contrast to the semi-Markovian process discussed in section 2.3) called a Lévy flight.

To avoid confusion, we stress that a Lévy flight refers to a random movement in space, where the length of the transitions is considered at discrete steps, but time is not involved. Lévy walks, on the other hand, attach a time “penalty”, by assigning a velocity to each transition in space. In the simplest case, this velocity is constant; relaxation of this condition leads back to the more general CTRW formulation of section 2.3 \[Klafter et al., 1987; Shlesinger et al., 1993\]. In any case, Lévy walks cannot be described in terms of simple fractional transport equations \[Metzler, 2000\].

A Lévy flight is characterized through the Fourier–Laplace transform \[Bouchaud and Georges, 1990; Compte, 1996; Metzler and Klafter, 2000\]

$$\stackrel{~}{𝒫}(𝐤,u)=\frac{1}{u+K^\mu |𝐤|^\mu }$$
(45)
from which, upon Fourier and Laplace inversion, the FDE \[Compte, 1996\]

$$\frac{}{t}P(𝐬,t)=K^\mu ^\mu P(𝐬,t)$$
(46)
is inferred. The Riesz operator $`^\mu `$ is defined through \[Samko et al., 1993\]

$$\left\{^\mu P(𝐬,t)\right\}=|𝐤|^\mu 𝒫(𝐤,t).$$
(47)
Note that we use the definition $`K^\mu c_{s}^{}{}_{}{}^{\mu }/c_t`$ for the diffusion constant. From (45), one recovers the characteristic function

$$𝒫(𝐤,t)=\mathrm{exp}\left(K^\mu t|𝐤|^\mu \right),$$
(48)
which is the characteristic function of a centered and symmetric Lévy distribution with the asymptotic power–law behavior \[Lévy, 1925, 1954; Gnedenko and Kolmogorov, 1954\]

$$P(𝐬,t)|𝐬|^{1\mu }.$$
(49)
Lévy distributions are used to generate Lévy flights \[Bouchaud and Georges, 1990\]. Accordingly, the second moment of a Lévy flight diverges:

$$𝐬(t)^2=\mathrm{}.$$
(50)
Observe that Lévy flights are characterized by a transition time distribution $`\psi (t)`$ with a finite first moment; they are thus fundamentally different from those processes underlying the time–fractional dispersion equation (44). As can be seen both descriptions are included in the CTRW framework.

Including a bias into the transition distribution, one obtains for an asymptotic form of $`p(𝐬)`$ the Lévy flight fractional ADE \[Metzler et al., 1998\]

$$\frac{}{t}P(𝐬,t)+𝐯P(𝐬,t)=K^\mu ^\mu P(𝐬,t)$$
(51)
which exhibits Galilei symmetry, i.e., (51) is solved by the Lévy stable solution (49), to be taken at the point $`𝐬𝐯t`$. This means that the symmetric Lévy stable plume is entirely shifted along the velocity vector $`𝐯`$, a situation which strongly contrasts the growing skewness in the CTRW case for long-tailed transition times. Of course, this solution features the same divergence (50) of the second moment of the plume distribution. The first moment of (51) exists for all $`0<\mu <2`$ and follows the usual Galilei symmetry expression

$$𝐬(t)=𝐯t.$$
(52)
### 3.3. Applications

As discussed above, although both time and space FDE forms are special cases of the CTRW, and both represent generalizations of the Fickian-based ADE, there are clear and critical distinctions between the transport equations that result from these two formulations. Here, we assess the Lévy flight description and argue that its characteristics strongly limit its applicability to describing transport in geological formations.

We consider the underlying physical picture of the Lévy flight, as applied to tracer migration in geological formations: a necessary condition for the Lévy flight description is that the domain clearly contain “streaks” of high and low permeability, arranged so as to lead to particle transitions of high and low velocity. In other words, the physical picture of a Lévy flight requires an encounter with a wide range of lengths of permeability streaks to obtain a non-Fickian distribution of particle transitions. And yet, such non-Fickian distributions arise even without the presence of such a permeability distribution, as clearly demonstrated by, e.g., Silliman and Simpson .

In addition, we observe that in mathematical terms, the first and second moments are often used to characterize plume migration. These quantities describe the spatio-temporal distribution of the tracer particles; the particles carry a finite mass, and therefore have a finite velocity. As noted above, the Lévy flight description leads to a diverging second moment of the migrating plume. Given that the macrodispersion parameter is typically defined in terms of the second moment, this divergence property cannot be ignored. Moreover, we observe that through scaling arguments \[Jespersen et al., 1999\], transport only undergoes a “superdiffusive” (faster than linear) process; in the Lévy flight description, subdiffusive transport can never occur.

With respect to the issue of a diverging second moment, one might attempt to work with a finite number of sampled tracer particles in a finite range, during a finite time window; this leads to a truncated Lévy distribution with finite moments. For truncated Lévy distributions it is known that their scaling behaviors in time pertain up to relatively large times \[Mantegna and Stanley, 1994, 1995\]. The difficulty is that to account for the temporal evolution of the particle cloud, the cutoffs would have to be adjusted to the actual space volume explored by the tracer particles, i.e., the cutoffs would themselves become time-dependent \[Jespersen et al., 1999\]. Put somewhat differently, the spatial-fractional formulation is based on an assumed fractal, scale-free nature of the transport process. Truncating the distributions leads, by definition, to a scale-dependent process which invalidates the use of simple fractional operators.

In contrast to the above arguments, the formulation given by, e.g., (44), or, more generally, by (16)-(19), does not suffer from these limitations or assumptions. In realistic field situations, the distribution of particle velocities is expected to vary widely on the order of magnitude of typical spacing between sampling points. Of course, the velocity distribution is bounded by some maximum velocity. In the long time limit, corresponding to the small $`u`$ limit that is of interest in our modeling, the mean effect of this finite variation of velocities can be approximated by a typical velocity. From this point of view, therefore, anomalies in the plume and the related moments should arise from temporal “sticking” processes (i.e., low velocity particle transitions) which are taken into consideration in the CTRW picture. Depending on the range of $`\beta `$ (recall (31)– (34)), both subdiffusive and superdiffusive behaviors for plume spreading can be characterized. Moreover, explicit spatial structure (well-defined conductivity features) can be incorporated within the CTRW framework.

## 4. The use of CTRW-based ensemble averages in non-stationary media: The relation of field scales and uncertainty

We return now to consider the issue, raised in the Introduction, that the interplay between ensemble averaging and spatial scales of non-stationary geological features strongly affects efforts to model transport. Broadly speaking, there exist two approaches to modeling transport in large, field-scale formations. In the first approach, the formation is treated as a single domain, with heterogeneities characterized and distributed according to a random field, with or without correlation and/or anisotropy. Generally speaking, these characterizations treat the domain as a stationary system, although stochastic models that incorporate a deterministic drift component (in the random field generator) have been considered \[e.g., Li and McLaughlin, 1995\]. In the second approach, a physical picture of the domain is constructed which includes explicitly specified (prescribed or known) heterogeneities, so that the resulting domains are non-stationary \[e.g., LaBolle and Fogg, 2001; Koltermann and Gorelick, 1996; Eggleston and Rojstaczer, 1998; Feehley et al., 2000\].

While the study of ensemble-averaged (stationary) domains has given rise to a sub-literature on stochastic methodologies and limiting behavior (e.g., perturbation techniques, macrodispersion) it has not yielded a practical numerical scheme to deal with the large majority of field sites. Anderson  describes in detail heterogeneity and trending structures evident in natural geological formations, and argues convincingly for the need to use facies modeling (coupled with geostatistical techniques) and/or depositional simulation models. These models can provide the underlying hydraulic conductivity structure and flow field of non-stationary domains, conditioned on field measurements, and be integrated with predictive models of transport.

Within the framework of non-stationary domains, explicitly characterized by structural trends, the question then arises as to how best to model transport (or, more precisely, how to deal with the unresolved heterogeneities (residues)). Clearly, there is a critical interplay between length scales associated with the trends and the residues. This gives rise to the associated uncertainty in both the measured/estimated hydraulic parameters and the measured/predicted concentrations. The generally accepted explanation for non-Fickian transport is that heterogeneities which cannot be ignored are present at all scales. Therefore, accounting for these residues is a central consideration for the quantification of non-Fickian transport.

In efforts to combine non-stationarity with local-scale heterogeneity and uncertainty, several recent studies have attempted to use ADE-based modeling approaches in conjunction with facies modeling \[e.g., Eggleston and Rojstaczer, 1998; Feehley et al., 2000\]. However, these studies, which incorporated even highly discretized systems (e.g., with block sizes of the order of 10 m3 in large aquifers), demonstrated an inability to adequately capture the migration patterns; these results suggests that unresolved heterogeneities also exist at these relatively small scales. We note that non-Fickian transport has been observed even in small-scale, relatively homogeneous, laboratory-scale models \[Berkowitz et al., 2000\]. Other related issues that have been considered recently focus on the relative importance of diffusion and local-scale dispersion and on how to separate diffusive mass transfer processes from slow particle velocities \[e.g., Harvey and Gorelick, 2000; LaBolle and Fogg, 2001\]. These questions may be considered to be somewhat moot, especially given that “dispersion” is an artifact of averaging in mathematical formulations, while a definitive separation between diffusion and very low velocity may be unnecessary.

At all of these smaller scales, i.e., within individual facies or depositional structures, the CTRW-based transport equations are highly effective. We therefore suggest that the CTRW-based approach should be used together with these facies and depositional models. As is usually done, a numerical model can be constructed which accounts explicitly for the heterogeneity structure of a formation, and the usual methods to solve for the flow field can be implemented. A CTRW-based transport equation can then be applied, rather than the ADE, over the entire domain. We observe that while the ADE (and the usual definition of “dispersion”) is simpler to apply than the CTRW-based equation, the preceding discussion (both in this section and the previous ones) demonstrate that it cannot and should not generally be applied in realistic field situations.

In this context, we shall consider the use of a hybrid model: known conductivity structures are accounted for explicitly, and within each block (pixel or voxel) of a numerical model we use the CTRW to account for the residues. Precluding the use of $`\psi (𝐬,t)`$ with (spatial) Lévy forms, because the trends are included explicitly in the numerical model, we can start with (23) as a basis for our numerical treatment. The methods developed with the use of the ADE, can be carried out with the Laplace transforms of (23)-(25),

$`u\stackrel{~}{P}(𝐬,u)P_0(𝐬)`$ $`=`$ $`\stackrel{~}{𝐯}_\psi (u)\stackrel{~}{P}(𝐬,u)`$ (53)
$`+\stackrel{~}{𝚽}_\psi (u):\stackrel{~}{P}(𝐬,u)`$
$$\stackrel{~}{𝐯}_\psi (u)=\frac{u\mathrm{\Sigma }_𝐬\stackrel{~}{\psi }(𝐬,u)𝐬}{1\stackrel{~}{\psi }(u)}$$
(54)
$$\stackrel{~}{𝚽}_\psi (u)=\frac{u\mathrm{\Sigma }_𝐬\stackrel{~}{\psi }(𝐬,u)\frac{1}{2}\mathrm{𝐬𝐬}}{1\stackrel{~}{\psi }(u)}$$
(55)
where $`P_0(𝐬)`$ is the initial condition.

The transport equation (53) is very similar to the Laplace transform of the ADE, but with the important exception that $`\stackrel{~}{𝐯}_\psi `$ and $`\stackrel{~}{𝚽}_\psi `$ are $`u`$-dependent. A spatial grid can be employed to numerically solve (53), exactly as can be done with the ADE applied to a non-stationary system. At each grid point, the velocity value determined from the solution to the steady flow problem is used in (53)-(55), along with the corresponding estimate of $`\beta `$, to change the parameters of $`\stackrel{~}{\psi }(𝐬,u)`$ and $`\stackrel{~}{\psi }(u)`$.

In this methodology the interpretation of $`\stackrel{~}{\psi }(𝐬,u)`$ changes somewhat. Instead of single transitions, we consider $`\stackrel{~}{\psi }(𝐬,u)`$ as playing the role of accounting for the transition across an entire element of the spatial grid. This interpretation has been justified by Margolin and Berkowitz .

If we insert (recall (35))

$$\stackrel{~}{\psi }(u)1c_\beta u^\beta ,\mathrm{for}\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}0}<\beta <1$$
(56)
into (53)-(55), we generate non-Fickian transport across each block element (with $`c_\beta `$ proportional to the velocity value at the grid point, divided by a characteristic length, all raised to the $`\beta `$ power). The non-Fickian behavior is due to the unresolved heterogeneities below the scale of the spatial grid. Estimates of $`\beta `$ and $`c_\beta `$ can be obtained for each facies from a standard tracer breakthrough test and subsequent comparison and fitting with analytical solutions (as done, e.g., in Berkowitz et al.  and Kosakowski et al. ); this procedure is exactly analogous to the usual determination of the dispersivity parameter $`\alpha `$ in the ADE.

Using a more complete expression for $`\stackrel{~}{\psi }(𝐬,u)`$ we can also evolve the dynamics of the plume at very long time into a Gaussian (i.e., in a time regime in which $`\psi (𝐬,t)`$ possesses a finite first and second temporal moment). The change in $`\stackrel{~}{\psi }(𝐬,u)`$ across the boundaries can be handled by using suitable averages similar to the ADE-based numerical treatments. Hence one can numerically solve for $`\stackrel{~}{P}(𝐬,u)`$ at each grid point and obtain the normalized concentration $`P(𝐬,t)`$ by calculating $`^1[\stackrel{~}{P}(𝐬,u)]`$. However, the inversion of a Laplace transform can be challenging, and remains a key issue for future research.

Finally, if we include pumping wells at some of the grid points $`𝐬_p`$ (where $`\stackrel{~}{\psi }(𝐬_p,u)=0`$, because the particles enter the well but do not emerge), then we can obtain the accumulated concentration directly from $`\stackrel{~}{P}(𝐬_p,u0)`$. In other words, $`\stackrel{~}{P}(𝐬_p,0)=_0^{\mathrm{}}𝑑tP(𝐬_p,t)`$, and because mass is conserved, each pumping well acts as a sink extracting a fraction of the migrating particles.

## 5. Summary and Conclusions

The application of stochastic approaches to quantification of transport in heterogeneous geological media rests inevitably on the underlying conceptual picture of dispersive mechanisms. The fundamental significance of this picture was pointed out long ago. As noted by Bear , in his discussion of the work of Scheidegger , “…the application of the statistical approach requires…a choice of the type of statistics to be employed, i.e., the probability of occurrence of events during small time intervals within the chosen ensemble. This may take the form of correlation functions between velocities at different points or different times, or joint-probability densities of the local velocity components of the particle as functions of time and space or a probability of an elementary particle displacment. The chosen correlation function determines the type of dispersion equation derived.”

We have developed this early insight into a full, quantitative theory where the joint probability density is the $`\psi (𝐬,t)`$. This joint spatial-temporal distribution allows us to account for the behavior of migrating particles which can encounter a wide range of velocity regions in heterogeneity lenses of different spatial dimensions. This approach is in contrast to most others which have, historically, emphasized spatial formulations of transport equations, motivated by the clear spatial heterogeneity of geological formations.

The overarching framework for our physical picture of transport, and the assumptions (as detailed above) on particle transitions, is the Master Equation. This equation represents a general, yet highly applicable, quantification of transport which recognizes the broad spectrum of particle motions in space and time. We show, under a general assumption of the form of $`w(𝐬,𝐬^{})`$, that the Master Equation can be specialized in any single realization of the geological domain to a generalized form of the ADE.

The ensemble average of the unrestricted Master Equation leads to a Generalized Master Equation, which is exactly equivalent to the CTRW. As a limiting form, under highly restrictive conditions regarding the character of the transport (and therefore of the degree of structural heterogeneity), the conventional ADE can be recovered from this formulation.

Aquifers are inherently heterogeneous over a wide range of scales, and Fickian transport (embodied in the ADE) does not generally occur on practical scales of interest. We therefore suggest that the overwhelming focus on defining “effective” dispersion, or “macrodispersion” coefficients, in Fickian or pseudo-Fickian formulations of the transport problem, is misplaced for field-scale problems. The CTRW theory, which is the basis for our transport equation, quantifies naturally the non-Fickian behavior observed at laboratory and field scales, as well as in numerical simulations. The essential character of the transport can be embodied in an asymptotic form of the $`\psi (𝐬,t)`$, specifically by an exponent $`\beta `$. This exponent, which can be determined from the velocity distribution (based on solution of flow for a given conductivity field) or from a tracer test, parameterizes an entire class of non-Fickian plume evolutions, on scales larger than the size of the heterogeneities. Detailed discussions on the practical identification of $`\psi (𝐬,t)`$ and parameter values is given in Berkowitz and Scher , Kosakowski et al. , and Berkowitz et al. .

We have also shown how fractional derivative formulations of transport equations are special, asymptotic (limit) cases, (30) for $`\psi (𝐬,t)`$, of the CTRW theory. Inserting this limiting form (35) into the Laplace transform of (15), one arrives at the same step necessarily encountered at the outset of the solution of the FDE. Retention of the more general equation (15) has important advantages for a more complete modeling of the transport process. The limiting forms characterized by the exponent $`\beta `$ (which is the fractional order of the derivative in the FDE) apply for a certain time range only. Beyond this range, the $`\psi (𝐬,t)`$ changes in a manner that allows the plume to eventually assume a Gaussian shape (defined by “macrodispersion”) as is reasonable for most physical systems.

Finally, we consider how best to quantify contaminant transport in non-stationary geological formations. We delineate a hybrid approach in which known structural properties are included explicitly, and unresolved (unknown) heterogeneities at smaller scales are accounted for within the CTRW theory. Practical application of this approach is achieved by replacing the usual ADE equation that is integrated into numerical simulation codes by a CTRW-based transport equation. This transport model can be integrated with existing numerical modeling techniques to determine the underlying flow field.

We are currently focusing efforts on implementation of the solution technique suggested here, as well as on deriving analytical solutions for CTRW-based transport equations for forms of $`\psi (𝐬,t)`$ generalized in both space and time.

Appendix A

We showed how the use of (4)–(8) leads to the expression (10) for the first term of the right side of (3). We outline the derivation here for the second and third terms of the right side of (3), using these same equations. We have for the second term

$`{\displaystyle \underset{𝐬^{}}{}}w(𝐬,𝐬^{})(𝐬^{}𝐬)C(𝐬,t)`$ $``$ $`{\displaystyle \underset{𝐬^{}}{}}(F(|𝐬^{}𝐬|;𝐬)+(𝐬^{}𝐬)F)[\lambda +\frac{1}{2}(𝐬^{}𝐬)\pi )](𝐬^{}𝐬)C(𝐬,t)`$ (57)
$``$ $`{\displaystyle \underset{𝐬^{}}{}}\lambda F(|𝐬^{}𝐬|;𝐬)\frac{1}{2}(𝐬^{}𝐬){\displaystyle \frac{\pi }{\lambda }}(𝐬^{}𝐬)C(𝐬,t)`$
$`+{\displaystyle \underset{𝐬^{}}{}}(𝐬^{}𝐬)F\lambda (𝐬^{}𝐬)C(𝐬,t)`$
$`=`$ $`𝐃(𝐬){\displaystyle \frac{\pi }{\lambda }}C(𝐬,t)+2𝐃(𝐬)C(𝐬,t)`$
and for the third term,

$`{\displaystyle \underset{𝐬^{}}{}}w(𝐬,𝐬^{})\frac{1}{2}(𝐬^{}𝐬)(𝐬^{}𝐬):C(𝐬,t)_𝐬^{}F(|𝐬^{}𝐬|;𝐬)\lambda \frac{1}{2}(𝐬^{}𝐬)(𝐬^{}𝐬):C(𝐬,t)=𝐃(𝐬):C(𝐬,t)`$ (58)
We add the results of (A1), (A2) and (10) to obtain (12).

## Acknowledgments

The authors thank Marco Dentz and Gennady Margolin for useful discussions, and two anonymous reviewers for constructive comments. BB thanks the European Commission (Contract No. EVK1-CT-2000-00062) for financial support. RM acknowledges financial support from the Deutsche Forschungsgemeinschaft within the Emmy Noether program.

## References

Anderson, M. P., Characterization of geological heterogeneity, in Subsurface Flow and Transport: A Stochastic Approach, Dagan, G. and S. P. Neuman (Eds.), Cambridge University Press, New York, 23-43, 1997.

Barkai, E., and R. Silbey, Fractional Kramers equation, J. Phys. Chem., 104(16), 3866-3874, 2000.

Bear, J., Dynamics of Fluids in Porous Media, Elsevier, New York, 1972.

Berkowitz, B., and H. Scher, On characterization of anomalous dispersion in porous and fractured media, Water Resour. Res., 31(6), 1461-1466, 1995.

Berkowitz, B., and H. Scher, Anomalous transport in random fracture networks, Phys. Rev. Lett., 79(20), 4038-4041, 1997.

Berkowitz, B., and H. Scher, Theory of anomalous chemical transport in fracture networks, Phys. Rev. E, 57(5), 5858-5869, 1998.

Berkowitz, B., and H. Scher, The role of probabilistic approaches to transport theory in heterogeneous media, Transport Porous Media, 42(1-2), 241-263, 2001.

Berkowitz, B., H. Scher, and S. E. Silliman, Anomalous transport in laboratory-scale, heterogeneous porous media, Water Resour. Res., 36(1), 149-158, 2000. \[Minor correction: 36(5), 1371, 2000.\]

Berkowitz, B., G. Kosakowski, G. Margolin, and H. Scher, Application of continuous time random walk theory to tracer test measurements in fractured and heterogeneous porous media, Ground Water, 39(4), 593-604, 2001.

Bouchaud, J.-P., and A. Georges, Anomalous diffusion in disordered media - statistical mechanisms, models and physical applications, Phys. Rep., 195(4-5), 127-293, 1990.

Chandrasekhar, S., Stochastic problems in physics and astronomy, Rev. Mod. Phys., 15(1), 1-89, 1943.

Compte, A., Stochastic foundations of fractional dynamics, Phys. Rev. E, 53(4), 4191-4193, 1996.

Compte, A., Continuous time random walks on moving fluids, Phys. Rev. E, 55(6), 6821-6831, 1997.

Compte, A., and M. O. Càceres, Fractional dynamics in random velocity fields, Phys. Rev. Lett., 81(15), 3140-3143, 1998.

Compte, A., R. Metzler, and J. Camacho, Biased continuous time random walks between parallel plates, Phys. Rev. E, 56(2), 1445-1454, 1997.

Cushman, J. H., X. Hu, and T. R. Ginn, Nonequilibrium statistical mechanics of preasymptotic dispersion, J. Stat. Phys., 75(5/6), 859-878, 1994.

Dagan, G., and S. P. Neuman (eds.), Subsurface Flow and Transport: A Stochastic Approach, Cambridge University Press, New York, 1997.

Dagan, G., Stochastic modeling of flow and transport: The broad perspective, in Dagan, G. and S. P. Neuman (Eds.), Subsurface Flow and Transport: A Stochastic Approach, 3-19, Cambridge University Press, New York, 1997.

Deng, F.-W., J. H. Cushman, and J. W. Delleur, A fast Fourier transform stochastic analysis of the contaminant transport problem, Water Resour. Res., 29(9), 3241-3247, 1993.

Eggleston, J., and S. Rojstaczer, Identification of large-scale hydraulic conductivity trends and the influence of trends on contaminant transport, Water Resour. Res., 34(9), 2155-2168, 1998.

Feehley, C. E., C. Zheng, and F. J. Molz, A dual-domain mass transfer approach for modeling solute transport in heterogeneous aquifers: Application to the Macrodispersion Experiment (MADE) site, Water Resour. Res., 36(9), 2501-2515, 2000.

Garabedian, S. P., D. R. LeBlanc, L. W. Gelhar, and M. A. Celia, Large-scale natural gradient tracer test in sand and gravel, Cape Cod, Massachusetts 2. Analysis of spatial moments for a nonreactive tracer, Water Resour. Res., 27(5), 911-924, 1991.

Gelhar, L. W., C. Welty, and K. R. Rehfeldt, A critical review of data on field-scale dispersion in aquifers, Water Resour. Res., 28(7), 1955-1974, 1992.

Glimm, J., W. B. Lindquist, F. Pereira, and Q. Zhang, A theory of macrodispersion of the scale-up problem, Transp. Porous Media, 13, 97-122, 1993.

Gnedenko, B. V., and A. N. Kolmogorov, Limit Distributions for Sums of Random Variables, Addison–Wesley, Reading, 1954.

Haggerty, R., and S. M. Gorelick, Multiple-rate mass transfer for modeling diffusion and surface reactions in media with pore-scale heterogeneity, Water Resour. Res., 31(10), 2383-2400, 1995.

Harvey, C., and S. M. Gorelick, Rate-limited mass transfer or macrodispersion: Which dominates plume evolution at the Macrodispersion Experiment (MADE) site?, Water Resour. Res., 36(3), 637-650, 2000.

Hatano, Y., and N. Hatano, Dispersive transport of ions in column experiments: An explanation of long-tailed profiles, Water Resour. Res., 34(5), 1027-1033, 1998.

Hilfer, R., and L. Anton, Fractional master equations and fractal time random walks, Phys Rev. E, 51(2), R848-R851, 1995.

Hilfer, R. (Ed.), Applications of Fractional Calculus in Physics, 472 pp., World Scientific, Singapore, 2000.

Jespersen, S., R. Metzler, and H. C. Fogedby, Lévy flights in external force fields: Langevin and fractional Fokker-Planck equations and their solutions, Phys. Rev. E, 59(3), 2736-2745, 1999.

Kenkre, V. M., E. W. Montroll, and M. F. Shlesinger, Generalized master equations for continuous-time random walks, J. Stat. Phys., 9(1), 45-50, 1973.

Kinzelbach, W., Groundwater Modelling, Developments in Water Science, 25, Elsevier, Amsterdam, 1986.

Klafter, J., and R. Silbey, On electronic energy transfer in disordered systems, J. Chem. Phys., 72(2), 843-848, 1980a.

Klafter, J., and R. Silbey, Derivation of continuous-time random-walk equations, Phys. Rev. Lett., 44(2), 55-58, 1980b.

Klafter, J., A. Blumen, and M. F. Shlesinger, Stochastic pathway to anomalous diffusion, Phys. Rev. A, 35(7), 3081-3085, 1987.

Koltermann, C. E., and S. M. Gorelick, Heterogeneity in sedimentary deposits: A review of structure-imitating, process-imitating and descriptive approaches, Water Resour. Res., 32(9), 2617-2658, 1996.

Kosakowski, G., B. Berkowitz, and H. Scher, Analysis of field observations of tracer transport in a fractured till, J. Cont. Hydrol., 47, 29-51, 2001.

Labolle, E. M., and G. E. Fogg, Role of molecular diffusion in contaminant migration and recovery in an alluvial aquifer system, Transp. Porous Media, 42, 155-179, 2001.

Labolle, E. M., G. E. Fogg, and A. F. B. Tompson, Random-walk simulation of transport in heterogeneous porous media: Local mass-conservation problem and implementation methods, Water Resour. Res., 32(3), 583-593, 1996.

Lévy, P., Calcul des probabilités, Gauthier–Villars, Paris, 1925.

Lévy, P., Théorie de l’addition des variables aléatoires, Gauthier–Villars, Paris, 1954.

Li, S., and D. McLaughlin, Using the nonstationary spectral method to analyze flow through heterogeneous trending media, Water Resour. Res., 31(3), 541-552, 1995.

Mantegna, R. N., and H. E. Stanley, Stochastic process with ultraslow convergence to a Gaussian - The truncated Lévy flight, Phys. Rev. Lett., 73(22), 2946-2949, 1994.

Mantegna, R. N., and H. E. Stanley, Ultra-slow convergence to a Gaussian: The truncated Lévy flight, in Lévy Flights and Related Topics in Physics, edited by M. F. Shlesinger, G. M. Zaslavsky, and U. Frisch, pp. 301-312, Springer-Verlag, New York, 1995.

Margolin, G., and B. Berkowitz, Application of continuous time random walks to transport in porous media, J. Phys. Chem. B, 104(16), 3942-3947, 2000. \[Minor correction: 104(36), 8762, 2000.\]

Margolin, G., and B. Berkowitz, Spatial behavior of anomalous transport, Phys. Rev. E, 65, in press, 2002.

Matheron, G., and G. de Marsily, Is transport in porous media always diffusive? A counter example, Water Resour. Res., 16(5), 901-917, 1980.

Metzler, R., Generalized Chapman-Kolmogorov equation: A unifying approach to the description of anomalous transport in external fields, Phys. Rev. E, 62(5), 6233-6245, 2000.

Metzler, R., and A. Compte, Generalized diffusion-advection schemes and dispersive sedimentation: A fractional approach, J. Phys. Chem. B, 104(16), 3858-3865, 2000.

Metzler, R., and J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep., 339(1), 1-77, 2000.

Metzler, R., J. Klafter, and I. M. Sokolov, Anomalous transport in external fields: Continuous time random walks and fractional diffusion equations extended, Phys. Rev. E, 58(2), 1621-1633, 1998.

Neuman, S. P., Eulerian-Lagrangian theory of transport in space-time nonstationary velocity fields: Exact nonlocal formalism by conditional moments and weak approximations, Water Resour. Res., 29(3), 633-645, 1993.

Oldham, K. B., and J. Spanier, The Fractional Calculus, Academic Press, New York, 1974.

Oppenheim, I., K. E. Shuler, and G. H. Weiss, The Master Equation, MIT Press, Cambridge, 1977.

Rubin, Y., Transport of inert solutes by groundwater: recent developments and current issues, in Subsurface Flow and Transport: A Stochastic Approach, Dagan, G. and S. P. Neuman (eds.), 115-132, Cambridge University Press, Cambridge, 1997.

Samko, S. G., A. A. Kilbas, and O. I. Marichev, Fractional Integrals and Derivatives - Theory and Applications, Gordon and Breach, New York, 1993.

Scheidegger, A. E., Statistical hydrodynamics in porous media, J. Appl. Phys., 25, 994-1001, 1954.

Scheidegger, A. E., Statistical approach to miscible displacement in porous in porous media, Bull. Canad. Inst. Min. Met., 26-30, 1958.

Scher, H., and M. Lax, Stochastic transport in a disordered solid. I. Theory, Phys. Rev. B, 7(10), 4491-4502, 1973.

Scher, H., and E. W. Montroll, Anomalous transit-time dispersion in amorphous solids, Phys. Rev. B, 12(6), 2455-2477, 1975.

Scher, H., M. F. Shlesinger, and J. T. Bendler, Time-scale invariance in transport and relaxation, Physics Today, January, 26-34, 1991.

Shlesinger, M. F., Asymptotic solutions of continuous-time random walks, J. Stat. Phys., 10(5), 421-434, 1974.

Shlesinger, M. F., Random Processes, in Encyclopedia of Applied Physics, Vol. 16, VCH Publishers, Inc., New York, 1996.

Shlesinger, M. F., G. M. Zaslavsky, and J. Klafter, Strange kinetics, Nature, 363(6424), 31-37, 1993.

Silliman, S. E., and E. S. Simpson, Laboratory evidence of the scale effect in dispersion of solutes in porous media, Water Resour. Res., 23(8), 1667-1673, 1987.

Sposito, G., W. A. Jury, and V. Gupta, Fundamental problems in the stochastic convection-dispersion model of solute transport in aquifers and field soils, Water Resour. Res., 22(1), 77-88, 1986.

Zhang, Q., A multi-length-scale theory of the anomalous mixing-length growth for tracer flow in heterogeneous porous media,J. Stat. Phys., 66(1/2), 485-501, 1992.

Zumofen, G., and J. Klafter, Spectral random-walk of a single molecule, Chem. Phys. Lett., 219(3-4), 303-309, 1994.

Zumofen, G., J. Klafter, and A. Blumen, Trapping aspects in enhanced diffusion, J. Stat. Phys., 65(5/6), 991-1013, 1991.


text:
# Metal-Insulator Transition and Ferromagnetism in Diluted Magnetic Semiconductors

## Abstract

We have investigated the interplay between the metal-insulator transition and ferromagnetism in $`\mathrm{III}_{1x}\mathrm{Mn}_x\mathrm{V}`$ semiconductors. Our study is based on a model in which $`S=5/2`$ Mn local moments are exchange-coupled to valence band holes that interact via Coulomb interactions with each other, with ionized Mn acceptors, and with the antisite defects present in these materials. We find quasiparticle participation ratios consistent with a metal-insulator transition that occurs in the ferromagnetic state near $`x0.01`$. By evaluating the distribution of mean exchange fields at Mn moment sites, we provide evidence in favor of the applicability of impurity-band magnetic-polaron and hole-fluid models on insulating and metallic sides of the phase transition respectively.

Introduction.- Recent progress in the growth of diluted magnetic semiconductors semiconductors that exhibit ferromagnetism at relatively high temperatures has suggested exciting new possibilities for devices that combine information processing and storage functionalities in a single material. It is generally accepted that Mn acts as an acceptor in these semiconductors, that its half-filled d-shell contributes a $`S=5/2`$ local moment to the system’s low-energy degrees of freedom, and that ferromagnetism is due to interactions between local moments that are mediated by valence-band holes. Recent experiments demonstrate that ferromagnetism occurs in both metallic and insulating states, and that both magnetic and transport properties are sensitive to the Mn fraction $`x`$ and to the density of compensating antisite and other defects in the material. The highest ferromagnetic critical temperatures ($`T_c`$) appear to occur in the most metallic samples. The role of Coulomb interactions in the ferromagnetism in these materials is subtle. At high carrier densities, well on the metallic side of the metal-insulator transition, exchange and correlation in the hole system is expected to enhance ferromagnetism. Well on the insulator side of the transition, however, Coulomb interactions create a Mott gap that increases the importance of randomness in Mn ion positions, suppresses carrier-hopping between Mn sides, and eventually turns off the coupling between local moments that can lead to ferromagnetism.

In this Letter we report on a numerical Hartree-Fock study of model (III,Mn)V ferromagnets in both metallic and insulating regimes. By evaluating the distribution of $`T=0`$ Mn exchange mean-fields, defined precisely below, we find evidence that is generally supportive of the impurity-band magnetic-polaron picture that has been used to describe ferromagnetism in the insulating regime, and of the hole-fluid picture that has been used in the metallic regime. Our principle results are summarized in Figs. 1 and 2 where we plot $`\mathrm{𝑣𝑠}.T`$ the fraction of Mn sites that have exchange mean-fields larger than $`k_BT`$.

In the metallic case ($`x=0.05`$ and $`p=3.3\times 10^{20}\mathrm{cm}^3`$) the distribution of coupling strengths is peaked at a finite value close to its uniform hole-fluid value and all Mn ions are strongly coupled to holes. In the insulating case ($`x=0.01`$ and hole density $`p=2.8\times 10^{19}\mathrm{cm}^3`$), on the other hand, we find that that vast majority of Mn spins are nearly free the few strongly coupled Mn ions will form magnetic polarons that grow slowly in size and interact more strongly as the temperature is lowered. In the following sections we explain how these results were obtained and elaborate on their significance for models of diluted magnetic semiconductor ferromagnetism.

Model Hamiltonian.- Most theoretical work on (III,Mn)V ferromagnets has started from one of two idealized limits. Impurity-band models achieve simplification by assuming that the holes that mediate interactions between Mn ion moments are strongly localized, whereas hole-fluid models achieve simplification by neglecting disorder due to Mn acceptors and other defects and treating it perturbatively in estimating transport coefficients. An insulator to metal phase transition occurs in these ferromagnets as Mn and hole densities increase; experimentally it appears that the metal-insulator transition usually occurs near $`x0.01`$, likely depending on the density of compensating antisite defects, a quantity that is expected to be sensitive the details of sample-growth and annealing protocols. Although it is generally agreed that an impurity band picture should apply far on the insulator side of the transition and a hole-fluid model far on the metallic side of the transition, it has not been clear which approach is a better starting point in the experimentally relevant parameter ranges. We address this issue by examining the ferromagnetic ground state using a model, and an approximation scheme, that captures the physics of both limits correctly. Our study is built on $`\stackrel{}{k}\stackrel{}{p}`$ envelope-function approach description of the valence band and a Hartree-Fock description of interactions.

The four terms in the single-particle part of our band electron model Hamiltonian, $`H^0=H^K+H^{K.ex.}+H^{Mnh}+H^{ash}`$, require some discussion. $`H^{K.ex.}`$ represents the kinetic exchange interaction between Mn local moments, assumed to be aligned in the ferromagnetic system ground state, $`H^{Mnh}`$ represents the attractive Coulomb interaction between an ionized $`Mn^{2+}`$ acceptor and a valence band hole. In an envelope function formalism, central cell corrections to this interaction are necessary to achieve an accurate description of the isolated bound-acceptor limit. $`H^{ash}`$ describes the repulsive interaction between holes and antisite defects (represented as sites with charge +2). These defects act as deep donors partially compensating the Mn acceptors and reducing the overall hole density, and provide an important additional scattering center. $`H^K`$ is the usual kinetic energy term. In this study we ignore mixing between heavy and light hole bands by using a simple parabolic band approximation.

$`H^0`$ $`={\displaystyle \frac{\mathrm{}^2}{2m}}\stackrel{}{}+{\displaystyle \frac{1}{2}}S{\displaystyle \underset{I}{}}\widehat{\mathrm{\Omega }}_I\stackrel{}{\tau }J(\stackrel{}{r}\stackrel{}{R_I})`$ (1)
$`+`$ $`{\displaystyle \underset{I}{}}({\displaystyle \frac{e^2}{ϵ|\stackrel{}{r}\stackrel{}{R_I}|}}V_0e^{|\stackrel{}{r}\stackrel{}{R_I}|^2/r_0^2})+{\displaystyle \underset{K}{}}{\displaystyle \frac{2e^2}{ϵ|\stackrel{}{r}\stackrel{}{R_K}|}}.`$ (2)
Here $`J(\stackrel{}{r})=(J_{pd}/(2\pi a_0^2)^{3/2})\mathrm{exp}(r^2/2a_{0}^{}{}_{}{}^{2})`$, $`\stackrel{}{\tau }=(\tau _x,\tau _y,\tau _z)`$ is the Pauli spin matrix vector, I labels Mn sites, K labels antisite defect sites, and $`\widehat{\mathrm{\Omega }}_I`$ is the orientation of the $`I^{th}`$ Mn spin. The term in the potential proportional to $`V_0`$ is the central cell correction. Both Mn ions and antisites were distributed randomly in a cube of side $`L`$. The long range of the Coulomb interaction requires overall charge neutrality so that $`n_hN_{Mn}+2N_{as}=0`$, where $`n_h`$ is the density of holes, $`N_{Mn}`$ the density of Mn ions and $`N_{as}`$ the density of antisites. In the ground state, all Mn spins are oriented along the $`\widehat{z}`$ direction and $`H^0`$ is block diagonal in its spin indices. In the spin wave configurations discussed below, $`\widehat{x}`$ and $`\widehat{y}`$ components of the spins are present, doubling the dimension of the matrix that must be diagonalized numerically. The wavevector cutoff in the quasiparticle wavefunction expansion was tested by computing the binding energy of a hole to an isolated Mn, comparing with the results of Bhattacharjee and Benoit a la Guillaume who find that a binding energy of 112meV, 86 meV when the kinetic exchange term is neglected, and 68meV when the central cell correction is also neglected. Our results are 124meV, 88 meV, and 48meV respectively, demonstrating an adequate description of the completely isolated Mn limit.

In order to capture the correct physics of both metallic and impurity band limits, hole-hole interactions must be described using an approximation that accounts for screening effects in the metallic regime and avoids artificial self-interaction effects in the localized regime, motivating our use of Hartree-Fock (HF) theory . The HF quasiparticle Hamiltonian is $`H^{HF}=H^0+(V^H+V^X)`$ where the Hartree and Fock matrix elements, $`V_{\stackrel{}{k}\sigma ,\stackrel{}{k^{}}\sigma ^{}}^H`$ and $`V_{\stackrel{}{k}\sigma ,\stackrel{}{k^{}}\sigma }^X`$ , can be expressed in terms of the density matrix. We evaluate the density matrix in the $`k`$-space representation used to expand our envelope function quasiparticle wavefunctions: $`\rho _{\stackrel{}{k}\sigma ,\stackrel{}{k^{}}\sigma ^{}}=_\alpha n_\alpha c_{\stackrel{}{k}\sigma }^{(\alpha )}c_{\stackrel{}{k^{}}\sigma ^{}}^{(\alpha )}`$, where $`|\alpha =_{\stackrel{}{k}\sigma }c_{\stackrel{}{k}\sigma }^{(\alpha )}|\stackrel{}{k}\sigma `$, and $`n_\alpha `$ is a quasiparticle occupation number. Our HF scheme becomes exact in the strongly localized limit, since a localized quasiparticle does not interact with itself, but is expected to overestimate inhomogeneity in the metallic limit because it neglects quantum fluctuations in the many-hole state.

Participation ratios.- In order to verify that our model correctly describes the metal-insulator transition in these materials, we measured the localization properties of our Hartree-Fock quasiparticles by evaluating their participation ratios, $`P_\alpha =1/(L^3d^3r|\mathrm{\Psi }_\alpha (\stackrel{}{r})|^4)`$, where $`\mathrm{\Psi }_\alpha (\stackrel{}{r})=\stackrel{}{r}|\alpha `$ is a normalized quasiparticle wavefunction. (For a localized state $`P_\alpha (\xi _\alpha /L)^3`$, where $`\xi _\alpha `$ is the localization length.) Fig. 3 compares participation ratios evaluated at $`x=0.05`$ and $`0.0125`$, assuming compensated carrier densities $`n_h=N_{Mn}/3`$ and $`n_h=N_{Mn}/10`$, respectively. The participation fractions extrapolated to infinite volume are clearly finite and therefore consistent with metallic transport at $`x=0.05`$ whereas they are consistent with zero and localized quasiparticles at $`x=0.0125`$, in agreement with experiment.

Impurity-Band Picture vs.Hole-Fluid Picture.- The exchange mean-fields used to construct Fig. 1 and Fig. 2 are given by $`H_{eff,I}=𝑑\stackrel{}{r}J(\stackrel{}{r}\stackrel{}{R}_I)(n_{}(\stackrel{}{r})n_{}(\stackrel{}{r}))/2J_{pd}(n_{}(\stackrel{}{R}_I)n_{}(\stackrel{}{R}_I))/2`$ where the partial densities are determined by solving the Hartree-Fock equations self-consistently. In the extreme impurity band limit, $`H_{eff,I}`$ would have the value $`J_{pd}n_{max}/225\mathrm{m}\mathrm{e}\mathrm{V}`$ for a fraction $`n_h/N_{Mn}`$ of the local moments and much smaller values for all others. Our findings are in qualitative agreement, although the peak $`H_{eff}`$ value is somewhat smaller and there is no sharp peak in the distribution function near the peak value, presumably because of the variable Mn and antisite defect environment experienced by localized holes. (The shallow impurity Bohr radius $`a_B^{}`$ calculated using the heavy-hole mass is $`1.0\mathrm{nm}`$; the peak hole density, $`4.9\times 10^{20}\mathrm{cm}^3`$ is slightly larger than the shallow impurity value because of central cell corrections.) Kaminski and Das Sarma have recently estimated the impurity-band magnetic polaron model $`T_c`$, relating it to $`N_{Mn}`$, $`n_h`$ and the maximum mean-field exchange coupling by $`T_c0.5(N_{Mn}/n_h)^{1/2}(n_ha_{B}^{}{}_{}{}^{3})^{1/6}H_{max}\mathrm{exp}(0.86/n_{h}^{}{}_{}{}^{1/3}a_B^{})/k_B`$. Naively applying this formula to our $`x=0.01`$ case with the shallow impurity Bohr radius yields $`T_c15K`$. Note that $`T_c`$ is expected to decrease rapidly at lower hole densities because of the weakening magnetic polaron coupling. These estimated $`T_c`$’s are somewhat larger than what has typically been observed at $`x0.01`$ in (Ga,Mn)As, possibly because the theory neglects Pauli exclusion principle effects, superexchange interactions, that favor opposite orientations of nearby impurity-band holes and oppose ferromagnetism. Nevertheless, it appears that this picture provides an excellent qualitative description of ferromagnetism at $`x0.01`$ and smaller.

In the metallic case, the exchange mean-field distribution differs qualitatively, with strong couplings common and few weakly coupled moments. Our results likely overestimate the degree of density variation in samples with $`x0.05`$, both because we have neglected Mn acceptor-antisite correlations and because of our use of the Hartree- Fock approximation. (The corresponding self-consistent Hartree calculations lead to more sharply peaked Mn site majority-spin density distributions.) We note that the most likely hole density at a Mn site is approximately 1.5 times larger than the average hole density and that the most likely mean-field coupling strengths are again somewhat larger than the uniform hole-fluid value $`H_{fl}=nJ_{pd}/2190`$K.

There are important differences between metallic and insulating cases in the physics that controls $`T_c`$. The metallic $`T_c`$ can be limited either by the system’s stiffness against collective magnetization orientation variation, or if the stiffness is large enough by the competition between local moment entropy, exchange interactions, and the band energy cost of hole-spin polarization that is captured by the hole-fluid mean-field-theory $`T_c`$ expression. We have estimated the magnetic stiffness of both insulating and metallic states by evaluating the HF electronic energy cost of an imposed spin-wave with wavevector $`Q`$. If disorder were neglected, the procedure we follow in evaluating spin spin-wave energies differs from the theory of König et al. only through retardation effects.

Figure 4 displays average spin reversal energies for metallic and insulating states. Note that the energy cost of slow (long wavelength) spin-direction variations is extremely small in the insulating case because these are sensitive mainly to the weak interactions between magnetic polarons that are necessary for long-range magnetic order. It follows that mean-field will completely fail in estimating the critical temperature of insulating DMS ferromagnets. The energy cost for spin-direction variations that we find in the metallic case is in qualitative agreement with hole-fluid model results, although it is larger at short wavelengths. In the hole-fluid mean-field theory the energy cost of spin-direction variations would equal $`H_{fl}`$ at all wavevectors. It follows that for the parabolic band model we have studied substantial corrections would apply to its mean-field $`T_c`$. It is important to note, however, that valence band spin-orbit interactions stiffen collective spin-variations in these ferromagnets, reducing corrections to mean-field-theory $`T_c`$ estimates.

SREY is supported in part by the KOSEF Quantum-functional Semiconductor Research Center at Dongkuk University. Work at the University of Texas was supported by by the Welch foundation and by DARPA/ONR Award No. N00014-00-1-095.