- Email: [email protected]
Towards a unified framework for anomalous transport in heterogeneous media
Chemical Physics 284 (2002) 349–359 www.elsevier.com/locate/chemphys
Towards a unified framework for anomalous transport in heterogeneous media Harvey Scher, Gennady Margolin, Brian Berkowitz * Department of Environmental Sciences, Weizmann Institute of Science, 76100 Rehovot, Israel Received 30 January 2002
Abstract We develop a unified framework to model the anomalous transport of tracers in highly heterogeneous media. While the framework is general, our working media in this study are geological formations. The basis of our approach takes into account the different levels of uncertainty, often associated with spatial scale, in characterizing these formations. The effects on the transport of smaller spatial scale heterogeneities are treated probabilistically with a model based on a continuous time random walk (CTRW), while the larger scale variations are included deterministically. The CTRW formulation derives from the ensemble average of a disordered system, in which the transport in each realization is described by a Master Equation. A generic example of such a system – a 3D discrete fracture network (DFN) – is treated in detail with the CTRW formalism. The key step in our approach is the derivation of a physically based wðs; tÞ, the joint probability density for a displacement s with an event-time t. We relate the wðs; tÞ to the velocity spectrum UðnÞ ðjnj ¼ 1=v; n^ ¼ ^vÞ of the steady flow-field in a fluid-saturated DFN. Heterogeneous porous media are often characterized by a log-normal permeability distribution; the UðnÞ we use in this case is an analytic form approximating the velocity spectrum derived from this distribution. The common approximation of wðs; tÞ pðsÞt1b with a constant b, is evaluated in these cases. For the former case it is necessary to include s t coupling while the latter case points to the presence of an effective t-dependent b. The full range of these features can be included in the CTRW solution but, as is shown, not in the fractional-time derivative equation (FDE) formulation of CTRW. Finally, the methods used for the unified framework are critically examined. Ó 2002 Elsevier Science B.V. All rights reserved.
1. Introduction Nature abounds in the phenomenon of anomalous transport, wherein the standard transport coefficients, based on mean velocity and rate of * Corresponding author. Tel.: +972-8-934-2098; fax: +972-8934-4124. E-mail address: [email protected] (B. Berkowitz).
change of variance, are dependent on the time and/ or space scale of the measurement. The scales, e.g., range from nanoseconds (electronic transport/relaxation) to years (migration in groundwater/ catchment basins). The basic physical origin of this phenomenon is the sufficient encounter of the transporting species with a broad distribution of rate-limiting steps (e.g., rates due to: molecular charge transfers, absorption–desorption, relaxation, trap-release, multi-diffusion, movement in
0301-0104/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 1 - 0 1 0 4 ( 0 2 ) 0 0 5 5 8 - X
350
H. Scher et al. / Chemical Physics 284 (2002) 349–359
fracture fragments and/or local heterogeneities). In particular, it is the encounter with statistically rare event-times (reciprocal rates) that compete with the total accumulation of median event-times. Each type of physical medium presents a challenge to develop transport equations that are based on approximations retaining the response to the full spectrum of these rates (i.e., they do not ‘‘average-out’’ the important rare event-times). Electronic transport in disordered semiconductors and polymers/molecular-solutions has been based on the Master Equation (ME) [1,2], which is in widespread use in statistical mechanics (generally) and chemical kinetics (specifically). The ME conserves mass and incorporates all the transition rates between, e.g., sites, energy levels. The transition rates are computed by the local dynamics, while the ME accounts for the overall kinetics. In the case of hopping motion the ensemble average (over realizations of the disordered system) of the ME has been shown [3] to be a Generalized Master Equation (GME) which is fully equivalent to a continuous time random walk (CTRW). In the CTRW the key quantity embodying the rate spectrum is wðs; tÞ, the joint probability density for a displacement s with an event-time t. The CTRW has been successfully applied to anomalous electronic transport [4] where the physical models of hopping and/or multiple-trapping have been used to calculate wðs; tÞ. It was shown that in many of the materials studied a wðs; tÞ with a power law (algebraic tail) dependence on time at large t accounted for the experimental results over a considerable range of t. Moreover, it was also shown that the approximation wðs; tÞ t1b with a constant b, between 0 and 1, is an approximation over a time range that depends on the material system (e.g., site or molecular density [5,6], width of electron density-of-states [7], chemical-bond structure [8]) and the conditions of the measurement such as temperature and electric field [7]. In recent years the same approach has been applied to scale-dependent transport in geological formations [9–13]. Here the observations are the movement of passive or reacting tracers through aquifers, fracture networks, catchment basins, laboratory models of these environments and computer simulations of the processes. The field
systems are characterized by a very high degree of heterogeneity. The heterogeneities give rise to permeability variations on many spatial scales. The pertinent spatial and permeability scales often encompass over five orders of magnitude. The larger (spatial) scale permeability variations (trends) can be measured or approximated by geostatistical methods, but there are significant uncertainties at the smaller scales (residues). In a thoroughly analyzed tracer field experiment [14] it was concluded that these residues have a considerable effect on the transport (i.e., the residues cannot be averaged-out). The complexity of these geological formations has necessitated further development of our transport equation approach. We propose a unified framework in which the spatial regions comprising the residues are treated probabilistically (i.e., ensemble-averaged) by the CTRW and the larger scale permeability and/or head variations deterministically. First we discuss the separation at various scales of the kinetics (with the ME) and the local dynamics. We see how this works for the wðs; tÞ for a representative case like a discrete fracture network (DFN, where we need to use the coupled space-time w-function). We include a discussion of the relation of CTRW to the fractional derivative equation (FDE) formulation of this problem as a way of illustrating the need to use the more general CTRW when one has to consider non-constant b. Finally we outline some methods for the implementation of the unified framework.
2. Dynamics and kinetics The formations are considered to be saturated with water and the tracers move through the velocity field due to the steady-state flow of the water under a natural or applied head. There are two generic types of geological formations that are of general interest in the transport of tracers: (i) a network of intersecting fractures in a low permeability matrix; (ii) a heterogeneous porous medium. In (i) it is not difficult to envision the transport as transitions through fracture fragments (between fracture intersections). The velocity field is determined by the global solution of the
H. Scher et al. / Chemical Physics 284 (2002) 349–359
constant flow in the network and the transitions are determined by the dynamics of the tracer response to a local flow and the geometric properties of the fragment. This will be considered in more detail in the next section. The overall kinetics generated by these transitions are obtained by solving the ME, X X oCðs; tÞ ¼ wðs0 ; sÞCðs; tÞ þ wðs; s0 ÞCðs0 ; tÞ; ot 0 0 s s ð1Þ where Cðs; tÞ is the particle concentration at point s and time t, wðs; s0 Þ is the transition rate from s0 to s. The application of (1) to (ii) is more subtle. We can consider a range of length scales from the pore level to the residue level. On the pore scale the local dynamics determining the aggregate of rates fwðs; s0 Þg are flow and diffusion between pore bodies. The distribution of the transition rates is due to the (quenched) disorder of the pore/throat distribution [15]. The pore scale is simple to model but completely unrealistic to implement for a field site. Upscaling to the next higher reasonable length scale of fixed disorder is the residue scale. On the residue scale (the mean size of relevant heterogeneities) the rates are determined by constitutive relations between the local head and flow (i.e., effective permeabilities). The change of concentration DC at each point in a time increment is Dt (the net particle flux). This flux contributing to DC can vary considerably at different points in the medium due to the fluctuations of permeabilities. If we pursue the usual approach and assume smooth variation of all the relevant quantities (Cðs; tÞ; wðs; s0 Þ and the pressure field) over some local region, we can generate a generalized ADE (advection–dispersion equation) from (1) [16] for a single realization of the medium (cf. discussion of the ADE in Appendix A.1). However, averaging of this type can be unreliable in a system of very wide fluctuations about the mean. Instead, we perform an ensemble-average of (1) wherein in each local region we retain a full distribution of event-times. Assuming that the size of the system is much greater than the residue scale and that the spatial distribution of the heterogeneities approximates statistical stationarity, it can
351
be shown [3] that the result is of the form of a GME XZ t oP ðs; tÞ ¼ /ðs0 s; t t0 ÞP ðs; t0 Þ dt0 ot 0 s0 XZ t þ /ðs s0 ; t t0 ÞP ðs0 ; t0 Þ dt0 ; ð2Þ 0
s0
where P ðs; tÞ is the normalized concentration, and /ðs; tÞ is defined below in (7). (The generalization to distinctly non-stationary domains is considered in Section 5.) The form of (2) is in contrast to (1); it is non-local in time and the transition rates are stationary (i.e., depend only on the difference s s0 ) and time-dependent. It is straightforward to show [17,18] using the Laplace transform, that the GME is completely equivalent to a CTRW XZ t wðs s0 ; t t0 ÞRðs0 ; t0 Þ dt0 ; ð3Þ Rðs; tÞ ¼ 0
s0
where Rðs; tÞ is the probability per time for a walker to just arrive at site s at time t, and wðs; tÞ has been defined in the previous section. The correspondence between (2) and (3) is Z t P ðs; tÞ ¼ Wðt t0 ÞRðs; t0 Þ dt0 ; ð4Þ 0
where WðtÞ ¼ 1
Z
t
wðt0 Þ dt0
ð5Þ
0
is the probability for a walker to remain on a site, X wðtÞ ¼ wðs; tÞ ð6Þ s
and uw~ðs; uÞ /~ðs; uÞ ¼ ; 1 w~ðuÞ
ð7Þ
where the Laplace transform (L) of a function f ðtÞ is denoted by f~ðuÞ. Eqs. (3)–(6) are in the form of a convolution in space and time and can therefore be solved by use of Fourier and Laplace transforms [5]. The general solution is Pðk; uÞ ¼
1 w~ðuÞ 1 ; u 1 Kðk; uÞ
ð8Þ
352
H. Scher et al. / Chemical Physics 284 (2002) 349–359
where Pðk; uÞ, Kðk; uÞ are the Fourier transform (F) of P~ðs; uÞ, w~ðs; uÞ, respectively. Hence, once wðs; tÞ is determined one can, in principle, determine the entire propagator (particle plume) P ðs; tÞ; there is no a priori need to consider only the moments of P ðs; tÞ. The technical difficulty is the inverse Fourier and Laplace transform of (8). The identification of wðs; tÞ lies at the heart of the CTRW formulation. One procedure is to calculate the dynamics of the particle motion in terms of variables such as the local flow and geometry, as one would define the wðs; s0 Þ, and use a global distribution of the variables (e.g., fracture fragment length, velocity spectrum) to correctly weight them. As is often the case, for the results of these calculations (as shown below), in highly disordered media, the wðs; tÞ can be approximated by a power law (algebraic tail) dependence on time for a considerable range of t, i.e., wðs; tÞ t1b :
ð9Þ
tions, e.g., which are of particular interest as sites for underground repositories to store radioactive and toxic industrial wastes, consist of networks of interconnected fractures embedded in a low permeability host rock (Fig. 1). Fracture network connectivity and permeability are governed by fracture density, orientation, length (or areal extent), and aperture, as well as by clustering and anisotropy properties of natural fractures [21]. In the studies referred to here, the rock mass is considered to be completely saturated by groundwater, with a steady-state flow regime established (subject to fixed macroscopic, regional-scale, boundary conditions). DFN analyses are based on available field data involving geological and geophysical mapping of fracture structure, and experimental testing of hydraulic parameters. From this information, distributions are derived to describe fracture properties. These distributions, together with available information on particular fracture features, are used as input to numerically
The first and second temporal moments do not exist for 0 < b < 1, while the second temporal moment does not exist for 1 < b < 2. The propagators for these values of b have been calculated [19,20] and exhibit anomalous transport, i.e., the coefficients of the mean velocity and variance are scale dependent. One can directly determine the mean position ‘ðtÞ and standard deviation rðtÞ of P ðs; tÞ. In the presence of a pressure gradient (or ‘‘bias’’), and for (9) ‘ðtÞ tb ; ‘ðtÞ t;
rðtÞ tb ;
0 < b < 1;
rðtÞ tð3bÞ=2 ;
1 < b < 2:
ð10Þ ð11Þ
This is in contrast to the normal transport of, e.g., a Gaussian propagator where ‘ðtÞ t and rðtÞ t1=2 . In the next two sections we consider in some detail case studies of the determination of wðs; tÞ.
3. Discrete fracture networks Discrete fracture network models are applied frequently to characterize and quantify fluid flow and tracer (contaminant) transport in fractured geological formations. Crystalline rock forma-
Fig. 1. Illustration of a numerically generated, 3D fracture network. Fractures are represented as 2D polygonal ‘‘platelets’’, with (constant) apertures sampled from a prescribed statistical distribution. Shown here is a network containing fractures oriented and centered randomly (uniformly) in space, with a power-law length distribution.
H. Scher et al. / Chemical Physics 284 (2002) 349–359
generate fracture networks that are conditioned to mimic possible realizations of the study site. With the assignment of representative boundary conditions, fluid flow through the network can be solved for numerically. Then, contaminant ‘‘particles’’ are introduced to the fracture network and their migration through network, and arrival times at a predefined compliance plane, are monitored. Such a fracture network reconstruction and flow simulation was carried out to analyze a € sp€ crystalline rock formation at A o, Sweden (for the so-called TRUE-Block Scale Experiment). Over the last two decades, extensive geophysical and hydraulic testing has been carried out to characterize the formation. The ‘‘Fracman’’ software package has been used to implement the DFN model (Golder Associates). Disc-shaped fractures were generated in a 3D domain, according to the distributions that characterize the fracture statistics (Fig. 1). Intersections among fractures were identified, and all non-conducting fractures and fracture segments were removed, thus yielding the network ‘‘backbone’’. One-dimensional ‘‘pipes’’ (or ‘‘channels’’), with an equivalent permeability, were then defined to topologically replace fracture segments connecting each pair of intersections. Flow through the pipe network was determined for prescribed (fixed pressure) boundary conditions, and particle tracking through the network was then used to determine distribution of particle velocities through the network. Hydrostructural models and resulting velocity histograms were derived from calibrated simulations of a series of 10–100 m scale tracer tests in fracture networks, as performed by Golder Associates. As shown in Fig. 2, the log–log plot of the histogram of nð¼ 1=vÞ exhibits a slowly varying region for almost two decades of large n. Although the noise increases at the highest n, the slope of the ‘‘tail’’ over this region is about 1:5 0:2. Note also a similar slope for the small n range; this does not likely affect the long-time breakthrough tailing of particles, because there is not necessarily a direct correlation between high velocities and long travel distances (as shown for some 2D networks [9]). In order to calculate wðs; tÞ for this generic formation we need to first consider the tracer dy-
353
namics across a typical disc-like fracture fragment with, e.g., a given mean velocity and aperture variation. The simplest approximation for our present purpose is to assume a plug flow in the effective ‘‘channel.’’ The geometry distribution (length s and orientation h) of the ‘‘channels’’, f ðsÞ, and the mean, flux-weighted, (inverse) velocity histogram, UðnÞ, are mapped onto the probability rate for a transition to have a displacement s occurring at a time t, wðs; tÞ / UðnÞf ðsÞ;
ð12Þ
where n ¼ 1=v, t ¼ sn. The variation of s is usually fit with an exponential or log-normal distribution. Field experiments [22–25] frequently demonstrate the occurrence of anomalous transport. DFN particle tracking simulations mimic the observed anomalous transport; this transport is due to the power-law tail of UðnÞ ! n1b at large n (e.g., Fig. 2 with 1 þ b ¼ 1:5 0:2) which results in wðs; tÞ ! f ðsÞs1þb t1b . The s t coupling in this case is seen as an effective change in f ðsÞ to f ðsÞs1þb . It is also manifested (not shown) in the hdependence of UðnÞ, e.g., bðhÞ. The data are not sufficient yet to exhibit these features. The wðs; tÞ can now be used to calculate the propagator and breakthrough curves for the DFN (Fig. 1) and any reasonable extension of its size. Our main point here is to demonstrate the physical basis for using a wðs; tÞ with a power-law tail to describe tracer transport through an important form of geological formation. The nature of the spectrum UðnÞ in a heterogeneous medium plays a key role. We explore another case of this spectrum in the following section.
4. Heterogeneous media, CTRW and FDE Heterogeneous porous media are often characterized by a log-normal permeability distribution. Flow field simulations using this distribution [26,27] exhibit a skewed velocity distribution, especially for large permeability variance. An analytic form that approximates the central feature of the simulations (for scalar n) is a UðnÞ equal to a log-normal distribution, which can be written directly as wðtÞ after the substitution t0 ¼ sn,
354
H. Scher et al. / Chemical Physics 284 (2002) 349–359
Fig. 2. Plot of the velocity histogram (based on over 4000 data points) for particles transported through a fracture network, recon€ sp€ structed from the rock formation at A o, Sweden. Note the log–log scale. Shown here is the velocity frequency, UðnÞ, versus n ¼ 1=v (in units of d/m). The large n (low velocity) regime is characterized by a power-law behavior, n1b , where 1 þ b ¼ 1:5 0:2.
wðtÞ ¼
rffiffiffi c1 expðc ln2 tÞ; pt
ð13Þ
where we absorbed the s into t (i.e., t ¼ t0 =s), c 1=2b2 and b is the log standard deviation. As in (12) one has wðs; t0 Þ / wðt0 =sÞf ðsÞ. Our focus is wðtÞ, which we can consider the probability density for event-times with a single transition length in one direction. In [28] Monte Carlo simulations were performed using the log-normal wðtÞ (13) with different values of b ranging from 2 to 6. The first-passage-time and propagator (distributions) F ðt; xÞ and P ðx; tÞ, respectively, were generated at different distances x and times t, and fitted by the distributions obtained using the CTRW with wðtÞ Ab t1b . The quality of the fits is very good. It has been shown [28] that the effective b as a function of time scales as 2c ln t (in contrast to the naive prediction c ln t one would make from inspection of the form of wðtÞ (13)). It has also been conclusively demonstrated that if the ‘‘effective’’ fitted value of b is close to 1 (e.g., 0:8 ! 1), one must add terms in the small u expansion of w~ðuÞ (e.g., 1 cb ub þ c1 u þ ), which modifies the CTRW distributions calculated only with the leading term ub [20]. Otherwise, the quality of the
fits, with 0:8 < b < 1, deteriorates compared to the fits of F ðt; xÞ and P ðx; tÞ with smaller effective b and without the additional terms. Fig. 3 shows example best fits for Qðx; tÞ ðP ðx; tÞ ¼ oQðx; tÞ= oxÞ clearly demonstrating the fact that for b close to 1 the additional term, represented by a non-zero dimensionless parameter j, should be taken into account (cf. details below). Hence, this case deriving from a log-normal velocity distribution is another example (in addition to the fracture network leading to (12)) of the physical basis of a power law wðtÞ. An expression similar to (13) for wðtÞ was derived for electron hopping between localized states [5] (Fig. 4) with rates W ðrÞ ¼ Wo er=Ro , r ¼ the intersite distance, Ro one-half the effective Bohr radius, wðtÞ ¼ gððln csÞ2 þ a2 Þ Wo h i 3 exp 13 g ðln csÞ þ 3a2 ln cs þ 2a3 ; s ð14Þ where s Wo t, g 4pNR3o , N is the hopping site density and c, a2 , a3 are numbers. As one can
H. Scher et al. / Chemical Physics 284 (2002) 349–359
355
Fig. 3. Plot of Qðx; sÞ vs x for s ¼ 1000 (where s is in units of the first moment of log-normal wðtÞ) and b ¼ 4:25. The propagator P ðx; sÞ ¼ oQðx; sÞ=ox. The inset provides the fitting parameters b, j.
However, it is usually understood that in a physical system the effective value of b is not constant but is a function (although slow) of characteristic times and distances at which the transport process is modeled (as shown above with b 2c ln t and in Fig. 4). Thus, e.g., in Fig. 4 at large enough times the effective b becomes larger than 2 (for all g values) and the particle motion evolves to normal transport. In the following we discuss the details and some consequences of including both the u and the ub terms. It has been shown in [20,28] that in the case of no backflow, the mean and the variance of the 1D propagator for 0 < b < 2, b 6¼ 1 are given by m X jj ð15Þ ‘ðtÞ ’ R C½1 þ cð1 ejÞ j¼0 and Fig. 4. Plot of wðtÞ=Wo vs Wo t (modified from [6]). The numerical labels of each curve refer to values of g=2:17. The curve expðWo tÞ has been added for comparison.
r2‘ ðtÞ ’ R2
m X j¼0
observe in Fig. 4, for the lower values of g, there are large ranges of s where wðtÞ can be approximated by s1b .
j X n¼0
jj
2ðj þ 1Þ C½1 þ cð2 ejÞ
! 1 ; C½1 þ cð1 enÞ C½1 þ cð1 eðj nÞÞ ð16Þ
356
H. Scher et al. / Chemical Physics 284 (2002) 349–359
respectively, where c minðb; 1Þ, e maxðb; 1=bÞ 1 and m ½1=e. In the case 0 < b < 1; R tb =C and j C1 tb1 =C while for 1 < b < 2R t=C1 and j Ct1b =C1 with C and C1 constants. We stress here that b is the constant responsible for the functional behavior, i.e., the type of the process, while C and C1 mainly quantify it. We note that only the terms increasing with time are given for the mean and only the terms growing faster than tc are given for the variance. In approximating the mean and the variance, the summation can be truncated at any m1 6 m in (15) and (16) (keeping in mind that for c ¼ 1 the term with j ¼ 0 in (16) is zero) if this is justified by the smallness of j (which is dimensionless); for b 6 0:5 it is proper to use only the rougher approximation with j ¼ 0. If b is far from 1 then j quickly decreases with time and one can safely use only the first leading terms for the mean and the variance. This is not the case if b becomes close to 1. In this case j hardly changes with time and if it is of order of 1, taking into account more terms is then necessary, as has been demonstrated in [28]. It is important to stress here that one cannot simply resort to increasing t in order for j to decrease significantly, since then the ‘‘effective value’’ of b itself will change. The modification with the additional term(s) also affects the construction of FDE (cf. Appendix A.2) [29], based on algebraic tails with a constant power (and thus FDE are a subset of the more general CTRW). In considering a biased transport, it turns out that most of the FDE presented in the literature (e.g., (A.2), see [29]) may fail for quite large intervals of possible b values, as can be seen from the spatial moments behavior predicted by the CTRW (15) and (16). The reason is that all these FDE include only the terms 1 cb ub in the expansion of w~ðuÞ. However, certain correction terms of similar magnitude might appear. This feature becomes prominent for b being close to 1 (and presumably also close to 2 – this case has not been studied in detail, as yet). Moreover, because of an unfortunate ‘‘propagating error’’ it is generally erroneously believed that the case 1 < b < 2 corresponds already to normal transport; this is true, however, only for the unbiased walks as has been argued in [20].
Even in the biased case, the FDE, which are equivalent to the CTRW case with 0 < b < 1 and retention only of the first leading term, are not applicable at all in the region 1 < b < 2. For the FDE to be able to produce (15) and (16) m þ 1 terms must be included, each one having a fractional derivative of different order, leading to quite complicated equations (similar to the proposed Cattaneo-type equations [30]). For a specific class of wðtÞ (which also fulfills the asymptotic form (9)), the equivalence between CTRW and FDE can be shown over the entire range of t [31]. We finally mention that although we were talking about algebraic tails in time, the same is true in general for any algebraic tails – when the tail exponent is close to an integer more terms should be analyzed.
5. Unified framework We have described the application of CTRW to ensemble-averaged geological formations including the physical basis for (9) (which leads to anomalous transport for 0 < b < 2). 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 realistically model transport. We start with a physical picture of the domain, which is constructed to include specified (prescribed or known) heterogeneities, so that the resulting domains are non-stationary, e.g., [14,32–34]. 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)). These unresolved residues are present at all scales and are understood to be the source of anomalous transport [14]. We note that anomalous transport has been observed even in small-scale, relatively homogeneous, laboratory-scale models [11]. Therefore, accounting for the residues is a central consideration for optimizing the transport model of the domain. For each domain characterization there are length scales that separate levels of uncertainty. A
H. Scher et al. / Chemical Physics 284 (2002) 349–359
scale is chosen for the size of a block (pixel or voxel) and within this block we use the CTRW to account for the residues. The larger scale flow field (calculated with a known permeability structure) is resolved on the block size explicitly. The ensembleaverage over the realizations of a distribution of heterogeneities is only within the block while the trends are included explicitly in a numerical treatment over the entire domain. We can implement this approach in at least two ways. The first relates closely to ADE-based extensive numerical computations [32] of a prescribed structure. As discussed in the Introduction one can assume a smooth variation in space if the full spectrum of the rates has been retained. We insert a series expansion of P ðs; tÞ, P ðs0 ; tÞ P ðs; tÞ þ ðs0 sÞ rP ðs; tÞ 1 þ ðs0 sÞðs0 sÞ : rrP ðs; tÞ 2
ð17Þ
(with the dyadic symbol : denoting a tensor product) into (2), yielding Z oP ðs; tÞ X t 0 ¼ dt /ðs s0 ; t t0 Þðs0 sÞ ot 0 s0 1 rP ðs; t0 Þ þ /ðs s0 ; t t0 Þ ðs0 sÞ 2 ðs0 sÞ : rrP ðs; t0 Þ : ð18Þ We write (18) in a more compact form Z t
oP ðs; tÞ ¼ dt0 vw ðt t0 Þ rP ðs; t0 Þ ot 0 þ Uw ðt t0 Þ : rrP ðs; t0 Þ ; vw ðtÞ
X
/ðs; tÞs;
ð19Þ ð20Þ
s
Uw ðtÞ
X s
1 /ðs; tÞ ss: 2
ð21Þ
This particular formulation in (19) aligns us closer to terms that are familiar in the context of traditional modeling: the ‘‘effective velocity’’ vw and the ‘‘dispersion tensor’’ Uw . Note, however, that both of these terms are time-dependent, and most significantly, depend fundamentally on wðs; tÞ. This
357
equation has the form of an ADE generalized to non-local time responses as a result of the ensemble average. We can make the association clearer; the methods developed with the use of the ADE, can be carried out with the Laplace transforms of (19)– (21), uP~ðs; uÞ P0 ðsÞ ¼ ~vw ðuÞ rP~ðs; uÞ ~ w ðuÞ : rrP~ðs; uÞ; þU ~vw ðuÞ ¼
uRs w~ðs; uÞs ; 1 w~ðuÞ
1 ~ ~ w ðuÞ ¼ uRs wðs; uÞ 2 ss ; U 1 w~ðuÞ
ð22Þ ð23Þ
ð24Þ
where P0 ðsÞ is the initial condition. The transport equation (22) is very similar to the Laplace transform of the ADE (the 3D version of (A.1)), but with the important exception that ~vw , ~ w , are u-dependent. A spatial grid can be and U employed to numerically solve (22), exactly as can be done with the ADE applied to a non-stationary system. At each grid point, the velocity value v determined from the solution to the steady flow problem is used in (22)–(24), along with the corresponding estimate of b, to change the parameters of w~ðs; uÞ and w~ðuÞ. In this methodology w~ðs; uÞ corresponds to single transitions over elements within the spatial grid. In light of the analysis in the previous section a complete expression for w~ðs; uÞ should be used which allows for evolution of the power of the algebraic tail. Hence one can numerically solve for P~ðs; uÞ at each grid point and obtain the normalized concentration P ðs; tÞ by calculating L1 ½P~ðs; uÞ. The second method of implementing our unified framework does not involve the assumption of the expansion in (17). We retain the structure of a CTRW on a lattice where the cell size is the block. As above the wðs; tÞ at each block is different; it is modified by the different dynamics (e.g., b) and flow field vðsb Þ at each block position sb . Problems in non-stationary systems where the wðs; tÞ depends explicitly on lattice position have been solved [35]. The solution depends on the inversion of a matrix whose size is the number of sb , sites.
358
H. Scher et al. / Chemical Physics 284 (2002) 349–359
Each element of the matrix is a lattice Green’s function connecting a pair of the sb sites. The resulting matrix is quite stable with respect to numerical inversion. The technique has been used to solve boundary value problems [36] and transport in a potential field [35]. It is a powerful technique and can be used practically in field applications because the block size can be much bigger than one usually associated with a large numerical computation because the wðs; tÞ already includes much of the smaller scale dynamics. For example, a good accounting of the data in an extensive fieldscale tracer test with the CTRW was carried out successfully up to a scale of 100 m [9]. We have discussed the physical basis of using wðs; tÞ with an algebraic tail in generic geological formations. We have analyzed the limitations of representing the w~ðuÞ by only the leading term in the small u expansion. Finally we assembled this approach into a unified framework for heterogeneous media by using an ensemble average in a component (block) in a larger explicit structural model.
Acknowledgements € sp€ A o test data were analyzed in cooperation with Golder Associates, Inc. for the Japan Nuclear Cycle Development Institute (JNC) participation € sp€ in the International A o Task Force. The authors thank SKB, Stockholm for permission to analyze these data, and Bill Dershowitz of Golder Associates Inc. for access to the data. The financial support of the Israel Science Foundation (Grant No. 147/01-1) is gratefully acknowledged.
Appendix A A.1. Advection–dispersion equation (ADE) The textbook approach to transport of passive tracers in both surface and subsurface systems usually focuses on the ADE [37] oWB =ot þ voWB =ox ¼ Ko2 WB =ox2 ;
ðA:1Þ
e.g., for 1D with constant v the average fluid velocity and K the dispersion coefficient. This equa-
tion governs the temporal evolution of the probability density function (pdf) WB ðx; tÞ of finding the tracer particle, which undergoes dispersive (or Brownian) motion, at a certain position x at time t after release at t ¼ 0. The pdf WB is the normalized concentration profile. A.2. Fractional derivative equations Mathematically, the CTRW (in the special case using wðtÞ (for 0 < b < 1Þ, and the small u expansion 1 cb ub for w~ðuÞ) is interchangeable with the time-fractional ADE [29] oW o o2 1b ¼ 0 Dt vb þ Kb 2 W ðx; tÞ; ðA:2Þ ot ox ox where the Riemann–Liouville operator is defined in terms of the convolution [29] Z t 1 d W ðx; t0 Þ 1b W ðx; tÞ ¼ dt0 ; 0 Dt 1b CðbÞ dt 0 ðt t0 Þ 0 < b < 1:
ðA:3Þ
The form of (A.2) is a natural extension of the ADE (A.1) for power-law forms of w. In this special case, the solutions of (A.2) reproduce those of the CTRW [29] in the spatial continuum limit. References [1] I. Oppenheim, K.E. Shuler, G.H. Weiss, The Master Equation, MIT Press, Cambridge, 1977. [2] M.F. Shlesinger, in: Encyclopedia of Applied Physics, vol. 16, VCH Publishers Inc., New York, 1996. [3] J. Klafter, R. Silbey, Phys. Rev. Lett. 44 (2) (1980) 55. [4] H. Scher, M.F. Shlesinger, J.T. Bendler, Phys. Today (January) (1991) 26, and references therein. [5] H. Scher, M. Lax, Phys. Rev. B 7 (10) (1973) 4491. [6] H. Scher, M. Lax, Phys. Rev. B 7 (10) (1973) 4502. [7] E. Muller-Horsche, D. Haarer, H. Scher, Phys. Rev. B 55 (1987) 1273. [8] M. Abkowitz, M.J. Rice, M. Stolka, Philos. Mag. B 61 (1990) 25. [9] B. Berkowitz, H. Scher, Phys. Rev. E 57 (5) (1998) 5858. [10] Y. Hatano, N. Hatano, Water Resour. Res. 34 (5) (1998) 1027. [11] B. Berkowitz, H. Scher, S.E. Silliman, Water Resour. Res. 36 (1) (2000) 149, Minor correction: 36 (5) 1371. [12] G. Kosakowski, B. Berkowitz, H. Scher, J. Cont. Hydrol. 47 (2001) 29.
H. Scher et al. / Chemical Physics 284 (2002) 349–359 [13] H. Scher, G. Margolin, R. Metzler, J. Klafter, B. Berkowitz, Geophys. Res. Lett. 29 (5) (2002), 10.1029/ 2001GL014123. [14] J. Eggleston, S. Rojstaczer, Water Resour. Res. 34 (9) (1998) 2155. [15] M. Blunt, H. Scher, Phys. Rev. E 52 (1995) 6387. [16] B. Berkowitz, J. Klafter, R. Metzler, H. Scher, Water Resour. Res., in press. [17] V.M. Kenkre, E.W. Montroll, M.F. Shlesinger, J. Stat. Phys. 9 (1) (1973) 45. [18] M.F. Shlesinger, J. Stat. Phys. 10 (5) (1974) 421. [19] G. Margolin, B. Berkowitz, J. Phys. Chem. B 104 (16) (2000) 3942, Minor correction: 104 (36) 8762. [20] G. Margolin, B. Berkowitz, Phys. Rev. E 65 (2002) 031101. [21] E. Bonnet, O. Bour, N.E. Odling, P. Davy, I. Main, P. Cowie, B. Berkowitz, Rev. Geophys. 39 (3) (2001) 347. [22] C.F. Tsang, Y.W. Tsang, F.V. Hale, Water Resour. Res. 27 (12) (1991) 3095. [23] H. Abelin, L. Bigersson, L. Moreno, H. Widen, T. Agren, I. Neretnieks, Water Resour. Res. 27 (12) (1991) 3119. [24] J. Hadermann, W. Heer, J. Cont. Hydrol. 21 (1–4) (1996) 87. [25] C.R. Sidle, B. Nilsson, M. Hansen, J. Fredericia, Water Resour. Res. 34 (10) (1998) 2515.
359
[26] A. Bellin, P. Salandin, A. Rinaldo, Water Resour. Res. 28 (9) (1992) 2211. [27] G. Margolin, B. Berkowitz, H. Scher, Water Resour. Res. 34 (9) (1998) 2103. [28] G. Margolin, B. Berkowitz, Phys. Rev. E, submitted. [29] R. Metzler, J. Klafter, Phys. Rep. 339 (1) (2000) 1, and references therein. [30] R. Metzler, T.F. Nonnenmacher, Phys. Rev. E 57 (6) (1998) 6409. [31] R. Hilfer, L. Anton, Phys. Rev. E 51 (2) (1995) R848. [32] E.M. Labolle, G.E. Fogg, Transp. Porous Media 42 (2001) 155. [33] C.E. Koltermann, S.M. Gorelick, Water Resour. Res. 32 (9) (1996) 2617. [34] C.E. Feehley, C. Zheng, F.J. Molz, Water Resour. Res. 36 (9) (2000) 2501. [35] H. Scher, S. Rackovsky, J. Phys. Chem. 81 (4) (1984) 1994. [36] H. Scher, E.W. Montroll, Phys. Rev. B 12 (6) (1975) 2455. [37] G. Dagan, S.P. Neuman (Eds.), Subsurface Flow and Transport: A Stochastic Approach, Cambridge University Press, Cambridge, UK, 1997.
Towards a unified framework for anomalous transport in heterogeneous media Harvey Scher, Gennady Margolin, Brian Berkowitz * Department of Environmental Sciences, Weizmann Institute of Science, 76100 Rehovot, Israel Received 30 January 2002
Abstract We develop a unified framework to model the anomalous transport of tracers in highly heterogeneous media. While the framework is general, our working media in this study are geological formations. The basis of our approach takes into account the different levels of uncertainty, often associated with spatial scale, in characterizing these formations. The effects on the transport of smaller spatial scale heterogeneities are treated probabilistically with a model based on a continuous time random walk (CTRW), while the larger scale variations are included deterministically. The CTRW formulation derives from the ensemble average of a disordered system, in which the transport in each realization is described by a Master Equation. A generic example of such a system – a 3D discrete fracture network (DFN) – is treated in detail with the CTRW formalism. The key step in our approach is the derivation of a physically based wðs; tÞ, the joint probability density for a displacement s with an event-time t. We relate the wðs; tÞ to the velocity spectrum UðnÞ ðjnj ¼ 1=v; n^ ¼ ^vÞ of the steady flow-field in a fluid-saturated DFN. Heterogeneous porous media are often characterized by a log-normal permeability distribution; the UðnÞ we use in this case is an analytic form approximating the velocity spectrum derived from this distribution. The common approximation of wðs; tÞ pðsÞt1b with a constant b, is evaluated in these cases. For the former case it is necessary to include s t coupling while the latter case points to the presence of an effective t-dependent b. The full range of these features can be included in the CTRW solution but, as is shown, not in the fractional-time derivative equation (FDE) formulation of CTRW. Finally, the methods used for the unified framework are critically examined. Ó 2002 Elsevier Science B.V. All rights reserved.
1. Introduction Nature abounds in the phenomenon of anomalous transport, wherein the standard transport coefficients, based on mean velocity and rate of * Corresponding author. Tel.: +972-8-934-2098; fax: +972-8934-4124. E-mail address: [email protected] (B. Berkowitz).
change of variance, are dependent on the time and/ or space scale of the measurement. The scales, e.g., range from nanoseconds (electronic transport/relaxation) to years (migration in groundwater/ catchment basins). The basic physical origin of this phenomenon is the sufficient encounter of the transporting species with a broad distribution of rate-limiting steps (e.g., rates due to: molecular charge transfers, absorption–desorption, relaxation, trap-release, multi-diffusion, movement in
0301-0104/02/$ - see front matter Ó 2002 Elsevier Science B.V. All rights reserved. PII: S 0 3 0 1 - 0 1 0 4 ( 0 2 ) 0 0 5 5 8 - X
350
H. Scher et al. / Chemical Physics 284 (2002) 349–359
fracture fragments and/or local heterogeneities). In particular, it is the encounter with statistically rare event-times (reciprocal rates) that compete with the total accumulation of median event-times. Each type of physical medium presents a challenge to develop transport equations that are based on approximations retaining the response to the full spectrum of these rates (i.e., they do not ‘‘average-out’’ the important rare event-times). Electronic transport in disordered semiconductors and polymers/molecular-solutions has been based on the Master Equation (ME) [1,2], which is in widespread use in statistical mechanics (generally) and chemical kinetics (specifically). The ME conserves mass and incorporates all the transition rates between, e.g., sites, energy levels. The transition rates are computed by the local dynamics, while the ME accounts for the overall kinetics. In the case of hopping motion the ensemble average (over realizations of the disordered system) of the ME has been shown [3] to be a Generalized Master Equation (GME) which is fully equivalent to a continuous time random walk (CTRW). In the CTRW the key quantity embodying the rate spectrum is wðs; tÞ, the joint probability density for a displacement s with an event-time t. The CTRW has been successfully applied to anomalous electronic transport [4] where the physical models of hopping and/or multiple-trapping have been used to calculate wðs; tÞ. It was shown that in many of the materials studied a wðs; tÞ with a power law (algebraic tail) dependence on time at large t accounted for the experimental results over a considerable range of t. Moreover, it was also shown that the approximation wðs; tÞ t1b with a constant b, between 0 and 1, is an approximation over a time range that depends on the material system (e.g., site or molecular density [5,6], width of electron density-of-states [7], chemical-bond structure [8]) and the conditions of the measurement such as temperature and electric field [7]. In recent years the same approach has been applied to scale-dependent transport in geological formations [9–13]. Here the observations are the movement of passive or reacting tracers through aquifers, fracture networks, catchment basins, laboratory models of these environments and computer simulations of the processes. The field
systems are characterized by a very high degree of heterogeneity. The heterogeneities give rise to permeability variations on many spatial scales. The pertinent spatial and permeability scales often encompass over five orders of magnitude. The larger (spatial) scale permeability variations (trends) can be measured or approximated by geostatistical methods, but there are significant uncertainties at the smaller scales (residues). In a thoroughly analyzed tracer field experiment [14] it was concluded that these residues have a considerable effect on the transport (i.e., the residues cannot be averaged-out). The complexity of these geological formations has necessitated further development of our transport equation approach. We propose a unified framework in which the spatial regions comprising the residues are treated probabilistically (i.e., ensemble-averaged) by the CTRW and the larger scale permeability and/or head variations deterministically. First we discuss the separation at various scales of the kinetics (with the ME) and the local dynamics. We see how this works for the wðs; tÞ for a representative case like a discrete fracture network (DFN, where we need to use the coupled space-time w-function). We include a discussion of the relation of CTRW to the fractional derivative equation (FDE) formulation of this problem as a way of illustrating the need to use the more general CTRW when one has to consider non-constant b. Finally we outline some methods for the implementation of the unified framework.
2. Dynamics and kinetics The formations are considered to be saturated with water and the tracers move through the velocity field due to the steady-state flow of the water under a natural or applied head. There are two generic types of geological formations that are of general interest in the transport of tracers: (i) a network of intersecting fractures in a low permeability matrix; (ii) a heterogeneous porous medium. In (i) it is not difficult to envision the transport as transitions through fracture fragments (between fracture intersections). The velocity field is determined by the global solution of the
H. Scher et al. / Chemical Physics 284 (2002) 349–359
constant flow in the network and the transitions are determined by the dynamics of the tracer response to a local flow and the geometric properties of the fragment. This will be considered in more detail in the next section. The overall kinetics generated by these transitions are obtained by solving the ME, X X oCðs; tÞ ¼ wðs0 ; sÞCðs; tÞ þ wðs; s0 ÞCðs0 ; tÞ; ot 0 0 s s ð1Þ where Cðs; tÞ is the particle concentration at point s and time t, wðs; s0 Þ is the transition rate from s0 to s. The application of (1) to (ii) is more subtle. We can consider a range of length scales from the pore level to the residue level. On the pore scale the local dynamics determining the aggregate of rates fwðs; s0 Þg are flow and diffusion between pore bodies. The distribution of the transition rates is due to the (quenched) disorder of the pore/throat distribution [15]. The pore scale is simple to model but completely unrealistic to implement for a field site. Upscaling to the next higher reasonable length scale of fixed disorder is the residue scale. On the residue scale (the mean size of relevant heterogeneities) the rates are determined by constitutive relations between the local head and flow (i.e., effective permeabilities). The change of concentration DC at each point in a time increment is Dt (the net particle flux). This flux contributing to DC can vary considerably at different points in the medium due to the fluctuations of permeabilities. If we pursue the usual approach and assume smooth variation of all the relevant quantities (Cðs; tÞ; wðs; s0 Þ and the pressure field) over some local region, we can generate a generalized ADE (advection–dispersion equation) from (1) [16] for a single realization of the medium (cf. discussion of the ADE in Appendix A.1). However, averaging of this type can be unreliable in a system of very wide fluctuations about the mean. Instead, we perform an ensemble-average of (1) wherein in each local region we retain a full distribution of event-times. Assuming that the size of the system is much greater than the residue scale and that the spatial distribution of the heterogeneities approximates statistical stationarity, it can
351
be shown [3] that the result is of the form of a GME XZ t oP ðs; tÞ ¼ /ðs0 s; t t0 ÞP ðs; t0 Þ dt0 ot 0 s0 XZ t þ /ðs s0 ; t t0 ÞP ðs0 ; t0 Þ dt0 ; ð2Þ 0
s0
where P ðs; tÞ is the normalized concentration, and /ðs; tÞ is defined below in (7). (The generalization to distinctly non-stationary domains is considered in Section 5.) The form of (2) is in contrast to (1); it is non-local in time and the transition rates are stationary (i.e., depend only on the difference s s0 ) and time-dependent. It is straightforward to show [17,18] using the Laplace transform, that the GME is completely equivalent to a CTRW XZ t wðs s0 ; t t0 ÞRðs0 ; t0 Þ dt0 ; ð3Þ Rðs; tÞ ¼ 0
s0
where Rðs; tÞ is the probability per time for a walker to just arrive at site s at time t, and wðs; tÞ has been defined in the previous section. The correspondence between (2) and (3) is Z t P ðs; tÞ ¼ Wðt t0 ÞRðs; t0 Þ dt0 ; ð4Þ 0
where WðtÞ ¼ 1
Z
t
wðt0 Þ dt0
ð5Þ
0
is the probability for a walker to remain on a site, X wðtÞ ¼ wðs; tÞ ð6Þ s
and uw~ðs; uÞ /~ðs; uÞ ¼ ; 1 w~ðuÞ
ð7Þ
where the Laplace transform (L) of a function f ðtÞ is denoted by f~ðuÞ. Eqs. (3)–(6) are in the form of a convolution in space and time and can therefore be solved by use of Fourier and Laplace transforms [5]. The general solution is Pðk; uÞ ¼
1 w~ðuÞ 1 ; u 1 Kðk; uÞ
ð8Þ
352
H. Scher et al. / Chemical Physics 284 (2002) 349–359
where Pðk; uÞ, Kðk; uÞ are the Fourier transform (F) of P~ðs; uÞ, w~ðs; uÞ, respectively. Hence, once wðs; tÞ is determined one can, in principle, determine the entire propagator (particle plume) P ðs; tÞ; there is no a priori need to consider only the moments of P ðs; tÞ. The technical difficulty is the inverse Fourier and Laplace transform of (8). The identification of wðs; tÞ lies at the heart of the CTRW formulation. One procedure is to calculate the dynamics of the particle motion in terms of variables such as the local flow and geometry, as one would define the wðs; s0 Þ, and use a global distribution of the variables (e.g., fracture fragment length, velocity spectrum) to correctly weight them. As is often the case, for the results of these calculations (as shown below), in highly disordered media, the wðs; tÞ can be approximated by a power law (algebraic tail) dependence on time for a considerable range of t, i.e., wðs; tÞ t1b :
ð9Þ
tions, e.g., which are of particular interest as sites for underground repositories to store radioactive and toxic industrial wastes, consist of networks of interconnected fractures embedded in a low permeability host rock (Fig. 1). Fracture network connectivity and permeability are governed by fracture density, orientation, length (or areal extent), and aperture, as well as by clustering and anisotropy properties of natural fractures [21]. In the studies referred to here, the rock mass is considered to be completely saturated by groundwater, with a steady-state flow regime established (subject to fixed macroscopic, regional-scale, boundary conditions). DFN analyses are based on available field data involving geological and geophysical mapping of fracture structure, and experimental testing of hydraulic parameters. From this information, distributions are derived to describe fracture properties. These distributions, together with available information on particular fracture features, are used as input to numerically
The first and second temporal moments do not exist for 0 < b < 1, while the second temporal moment does not exist for 1 < b < 2. The propagators for these values of b have been calculated [19,20] and exhibit anomalous transport, i.e., the coefficients of the mean velocity and variance are scale dependent. One can directly determine the mean position ‘ðtÞ and standard deviation rðtÞ of P ðs; tÞ. In the presence of a pressure gradient (or ‘‘bias’’), and for (9) ‘ðtÞ tb ; ‘ðtÞ t;
rðtÞ tb ;
0 < b < 1;
rðtÞ tð3bÞ=2 ;
1 < b < 2:
ð10Þ ð11Þ
This is in contrast to the normal transport of, e.g., a Gaussian propagator where ‘ðtÞ t and rðtÞ t1=2 . In the next two sections we consider in some detail case studies of the determination of wðs; tÞ.
3. Discrete fracture networks Discrete fracture network models are applied frequently to characterize and quantify fluid flow and tracer (contaminant) transport in fractured geological formations. Crystalline rock forma-
Fig. 1. Illustration of a numerically generated, 3D fracture network. Fractures are represented as 2D polygonal ‘‘platelets’’, with (constant) apertures sampled from a prescribed statistical distribution. Shown here is a network containing fractures oriented and centered randomly (uniformly) in space, with a power-law length distribution.
H. Scher et al. / Chemical Physics 284 (2002) 349–359
generate fracture networks that are conditioned to mimic possible realizations of the study site. With the assignment of representative boundary conditions, fluid flow through the network can be solved for numerically. Then, contaminant ‘‘particles’’ are introduced to the fracture network and their migration through network, and arrival times at a predefined compliance plane, are monitored. Such a fracture network reconstruction and flow simulation was carried out to analyze a € sp€ crystalline rock formation at A o, Sweden (for the so-called TRUE-Block Scale Experiment). Over the last two decades, extensive geophysical and hydraulic testing has been carried out to characterize the formation. The ‘‘Fracman’’ software package has been used to implement the DFN model (Golder Associates). Disc-shaped fractures were generated in a 3D domain, according to the distributions that characterize the fracture statistics (Fig. 1). Intersections among fractures were identified, and all non-conducting fractures and fracture segments were removed, thus yielding the network ‘‘backbone’’. One-dimensional ‘‘pipes’’ (or ‘‘channels’’), with an equivalent permeability, were then defined to topologically replace fracture segments connecting each pair of intersections. Flow through the pipe network was determined for prescribed (fixed pressure) boundary conditions, and particle tracking through the network was then used to determine distribution of particle velocities through the network. Hydrostructural models and resulting velocity histograms were derived from calibrated simulations of a series of 10–100 m scale tracer tests in fracture networks, as performed by Golder Associates. As shown in Fig. 2, the log–log plot of the histogram of nð¼ 1=vÞ exhibits a slowly varying region for almost two decades of large n. Although the noise increases at the highest n, the slope of the ‘‘tail’’ over this region is about 1:5 0:2. Note also a similar slope for the small n range; this does not likely affect the long-time breakthrough tailing of particles, because there is not necessarily a direct correlation between high velocities and long travel distances (as shown for some 2D networks [9]). In order to calculate wðs; tÞ for this generic formation we need to first consider the tracer dy-
353
namics across a typical disc-like fracture fragment with, e.g., a given mean velocity and aperture variation. The simplest approximation for our present purpose is to assume a plug flow in the effective ‘‘channel.’’ The geometry distribution (length s and orientation h) of the ‘‘channels’’, f ðsÞ, and the mean, flux-weighted, (inverse) velocity histogram, UðnÞ, are mapped onto the probability rate for a transition to have a displacement s occurring at a time t, wðs; tÞ / UðnÞf ðsÞ;
ð12Þ
where n ¼ 1=v, t ¼ sn. The variation of s is usually fit with an exponential or log-normal distribution. Field experiments [22–25] frequently demonstrate the occurrence of anomalous transport. DFN particle tracking simulations mimic the observed anomalous transport; this transport is due to the power-law tail of UðnÞ ! n1b at large n (e.g., Fig. 2 with 1 þ b ¼ 1:5 0:2) which results in wðs; tÞ ! f ðsÞs1þb t1b . The s t coupling in this case is seen as an effective change in f ðsÞ to f ðsÞs1þb . It is also manifested (not shown) in the hdependence of UðnÞ, e.g., bðhÞ. The data are not sufficient yet to exhibit these features. The wðs; tÞ can now be used to calculate the propagator and breakthrough curves for the DFN (Fig. 1) and any reasonable extension of its size. Our main point here is to demonstrate the physical basis for using a wðs; tÞ with a power-law tail to describe tracer transport through an important form of geological formation. The nature of the spectrum UðnÞ in a heterogeneous medium plays a key role. We explore another case of this spectrum in the following section.
4. Heterogeneous media, CTRW and FDE Heterogeneous porous media are often characterized by a log-normal permeability distribution. Flow field simulations using this distribution [26,27] exhibit a skewed velocity distribution, especially for large permeability variance. An analytic form that approximates the central feature of the simulations (for scalar n) is a UðnÞ equal to a log-normal distribution, which can be written directly as wðtÞ after the substitution t0 ¼ sn,
354
H. Scher et al. / Chemical Physics 284 (2002) 349–359
Fig. 2. Plot of the velocity histogram (based on over 4000 data points) for particles transported through a fracture network, recon€ sp€ structed from the rock formation at A o, Sweden. Note the log–log scale. Shown here is the velocity frequency, UðnÞ, versus n ¼ 1=v (in units of d/m). The large n (low velocity) regime is characterized by a power-law behavior, n1b , where 1 þ b ¼ 1:5 0:2.
wðtÞ ¼
rffiffiffi c1 expðc ln2 tÞ; pt
ð13Þ
where we absorbed the s into t (i.e., t ¼ t0 =s), c 1=2b2 and b is the log standard deviation. As in (12) one has wðs; t0 Þ / wðt0 =sÞf ðsÞ. Our focus is wðtÞ, which we can consider the probability density for event-times with a single transition length in one direction. In [28] Monte Carlo simulations were performed using the log-normal wðtÞ (13) with different values of b ranging from 2 to 6. The first-passage-time and propagator (distributions) F ðt; xÞ and P ðx; tÞ, respectively, were generated at different distances x and times t, and fitted by the distributions obtained using the CTRW with wðtÞ Ab t1b . The quality of the fits is very good. It has been shown [28] that the effective b as a function of time scales as 2c ln t (in contrast to the naive prediction c ln t one would make from inspection of the form of wðtÞ (13)). It has also been conclusively demonstrated that if the ‘‘effective’’ fitted value of b is close to 1 (e.g., 0:8 ! 1), one must add terms in the small u expansion of w~ðuÞ (e.g., 1 cb ub þ c1 u þ ), which modifies the CTRW distributions calculated only with the leading term ub [20]. Otherwise, the quality of the
fits, with 0:8 < b < 1, deteriorates compared to the fits of F ðt; xÞ and P ðx; tÞ with smaller effective b and without the additional terms. Fig. 3 shows example best fits for Qðx; tÞ ðP ðx; tÞ ¼ oQðx; tÞ= oxÞ clearly demonstrating the fact that for b close to 1 the additional term, represented by a non-zero dimensionless parameter j, should be taken into account (cf. details below). Hence, this case deriving from a log-normal velocity distribution is another example (in addition to the fracture network leading to (12)) of the physical basis of a power law wðtÞ. An expression similar to (13) for wðtÞ was derived for electron hopping between localized states [5] (Fig. 4) with rates W ðrÞ ¼ Wo er=Ro , r ¼ the intersite distance, Ro one-half the effective Bohr radius, wðtÞ ¼ gððln csÞ2 þ a2 Þ Wo h i 3 exp 13 g ðln csÞ þ 3a2 ln cs þ 2a3 ; s ð14Þ where s Wo t, g 4pNR3o , N is the hopping site density and c, a2 , a3 are numbers. As one can
H. Scher et al. / Chemical Physics 284 (2002) 349–359
355
Fig. 3. Plot of Qðx; sÞ vs x for s ¼ 1000 (where s is in units of the first moment of log-normal wðtÞ) and b ¼ 4:25. The propagator P ðx; sÞ ¼ oQðx; sÞ=ox. The inset provides the fitting parameters b, j.
However, it is usually understood that in a physical system the effective value of b is not constant but is a function (although slow) of characteristic times and distances at which the transport process is modeled (as shown above with b 2c ln t and in Fig. 4). Thus, e.g., in Fig. 4 at large enough times the effective b becomes larger than 2 (for all g values) and the particle motion evolves to normal transport. In the following we discuss the details and some consequences of including both the u and the ub terms. It has been shown in [20,28] that in the case of no backflow, the mean and the variance of the 1D propagator for 0 < b < 2, b 6¼ 1 are given by m X jj ð15Þ ‘ðtÞ ’ R C½1 þ cð1 ejÞ j¼0 and Fig. 4. Plot of wðtÞ=Wo vs Wo t (modified from [6]). The numerical labels of each curve refer to values of g=2:17. The curve expðWo tÞ has been added for comparison.
r2‘ ðtÞ ’ R2
m X j¼0
observe in Fig. 4, for the lower values of g, there are large ranges of s where wðtÞ can be approximated by s1b .
j X n¼0
jj
2ðj þ 1Þ C½1 þ cð2 ejÞ
! 1 ; C½1 þ cð1 enÞ C½1 þ cð1 eðj nÞÞ ð16Þ
356
H. Scher et al. / Chemical Physics 284 (2002) 349–359
respectively, where c minðb; 1Þ, e maxðb; 1=bÞ 1 and m ½1=e. In the case 0 < b < 1; R tb =C and j C1 tb1 =C while for 1 < b < 2R t=C1 and j Ct1b =C1 with C and C1 constants. We stress here that b is the constant responsible for the functional behavior, i.e., the type of the process, while C and C1 mainly quantify it. We note that only the terms increasing with time are given for the mean and only the terms growing faster than tc are given for the variance. In approximating the mean and the variance, the summation can be truncated at any m1 6 m in (15) and (16) (keeping in mind that for c ¼ 1 the term with j ¼ 0 in (16) is zero) if this is justified by the smallness of j (which is dimensionless); for b 6 0:5 it is proper to use only the rougher approximation with j ¼ 0. If b is far from 1 then j quickly decreases with time and one can safely use only the first leading terms for the mean and the variance. This is not the case if b becomes close to 1. In this case j hardly changes with time and if it is of order of 1, taking into account more terms is then necessary, as has been demonstrated in [28]. It is important to stress here that one cannot simply resort to increasing t in order for j to decrease significantly, since then the ‘‘effective value’’ of b itself will change. The modification with the additional term(s) also affects the construction of FDE (cf. Appendix A.2) [29], based on algebraic tails with a constant power (and thus FDE are a subset of the more general CTRW). In considering a biased transport, it turns out that most of the FDE presented in the literature (e.g., (A.2), see [29]) may fail for quite large intervals of possible b values, as can be seen from the spatial moments behavior predicted by the CTRW (15) and (16). The reason is that all these FDE include only the terms 1 cb ub in the expansion of w~ðuÞ. However, certain correction terms of similar magnitude might appear. This feature becomes prominent for b being close to 1 (and presumably also close to 2 – this case has not been studied in detail, as yet). Moreover, because of an unfortunate ‘‘propagating error’’ it is generally erroneously believed that the case 1 < b < 2 corresponds already to normal transport; this is true, however, only for the unbiased walks as has been argued in [20].
Even in the biased case, the FDE, which are equivalent to the CTRW case with 0 < b < 1 and retention only of the first leading term, are not applicable at all in the region 1 < b < 2. For the FDE to be able to produce (15) and (16) m þ 1 terms must be included, each one having a fractional derivative of different order, leading to quite complicated equations (similar to the proposed Cattaneo-type equations [30]). For a specific class of wðtÞ (which also fulfills the asymptotic form (9)), the equivalence between CTRW and FDE can be shown over the entire range of t [31]. We finally mention that although we were talking about algebraic tails in time, the same is true in general for any algebraic tails – when the tail exponent is close to an integer more terms should be analyzed.
5. Unified framework We have described the application of CTRW to ensemble-averaged geological formations including the physical basis for (9) (which leads to anomalous transport for 0 < b < 2). 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 realistically model transport. We start with a physical picture of the domain, which is constructed to include specified (prescribed or known) heterogeneities, so that the resulting domains are non-stationary, e.g., [14,32–34]. 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)). These unresolved residues are present at all scales and are understood to be the source of anomalous transport [14]. We note that anomalous transport has been observed even in small-scale, relatively homogeneous, laboratory-scale models [11]. Therefore, accounting for the residues is a central consideration for optimizing the transport model of the domain. For each domain characterization there are length scales that separate levels of uncertainty. A
H. Scher et al. / Chemical Physics 284 (2002) 349–359
scale is chosen for the size of a block (pixel or voxel) and within this block we use the CTRW to account for the residues. The larger scale flow field (calculated with a known permeability structure) is resolved on the block size explicitly. The ensembleaverage over the realizations of a distribution of heterogeneities is only within the block while the trends are included explicitly in a numerical treatment over the entire domain. We can implement this approach in at least two ways. The first relates closely to ADE-based extensive numerical computations [32] of a prescribed structure. As discussed in the Introduction one can assume a smooth variation in space if the full spectrum of the rates has been retained. We insert a series expansion of P ðs; tÞ, P ðs0 ; tÞ P ðs; tÞ þ ðs0 sÞ rP ðs; tÞ 1 þ ðs0 sÞðs0 sÞ : rrP ðs; tÞ 2
ð17Þ
(with the dyadic symbol : denoting a tensor product) into (2), yielding Z oP ðs; tÞ X t 0 ¼ dt /ðs s0 ; t t0 Þðs0 sÞ ot 0 s0 1 rP ðs; t0 Þ þ /ðs s0 ; t t0 Þ ðs0 sÞ 2 ðs0 sÞ : rrP ðs; t0 Þ : ð18Þ We write (18) in a more compact form Z t
oP ðs; tÞ ¼ dt0 vw ðt t0 Þ rP ðs; t0 Þ ot 0 þ Uw ðt t0 Þ : rrP ðs; t0 Þ ; vw ðtÞ
X
/ðs; tÞs;
ð19Þ ð20Þ
s
Uw ðtÞ
X s
1 /ðs; tÞ ss: 2
ð21Þ
This particular formulation in (19) aligns us closer to terms that are familiar in the context of traditional modeling: the ‘‘effective velocity’’ vw and the ‘‘dispersion tensor’’ Uw . Note, however, that both of these terms are time-dependent, and most significantly, depend fundamentally on wðs; tÞ. This
357
equation has the form of an ADE generalized to non-local time responses as a result of the ensemble average. We can make the association clearer; the methods developed with the use of the ADE, can be carried out with the Laplace transforms of (19)– (21), uP~ðs; uÞ P0 ðsÞ ¼ ~vw ðuÞ rP~ðs; uÞ ~ w ðuÞ : rrP~ðs; uÞ; þU ~vw ðuÞ ¼
uRs w~ðs; uÞs ; 1 w~ðuÞ
1 ~ ~ w ðuÞ ¼ uRs wðs; uÞ 2 ss ; U 1 w~ðuÞ
ð22Þ ð23Þ
ð24Þ
where P0 ðsÞ is the initial condition. The transport equation (22) is very similar to the Laplace transform of the ADE (the 3D version of (A.1)), but with the important exception that ~vw , ~ w , are u-dependent. A spatial grid can be and U employed to numerically solve (22), exactly as can be done with the ADE applied to a non-stationary system. At each grid point, the velocity value v determined from the solution to the steady flow problem is used in (22)–(24), along with the corresponding estimate of b, to change the parameters of w~ðs; uÞ and w~ðuÞ. In this methodology w~ðs; uÞ corresponds to single transitions over elements within the spatial grid. In light of the analysis in the previous section a complete expression for w~ðs; uÞ should be used which allows for evolution of the power of the algebraic tail. Hence one can numerically solve for P~ðs; uÞ at each grid point and obtain the normalized concentration P ðs; tÞ by calculating L1 ½P~ðs; uÞ. The second method of implementing our unified framework does not involve the assumption of the expansion in (17). We retain the structure of a CTRW on a lattice where the cell size is the block. As above the wðs; tÞ at each block is different; it is modified by the different dynamics (e.g., b) and flow field vðsb Þ at each block position sb . Problems in non-stationary systems where the wðs; tÞ depends explicitly on lattice position have been solved [35]. The solution depends on the inversion of a matrix whose size is the number of sb , sites.
358
H. Scher et al. / Chemical Physics 284 (2002) 349–359
Each element of the matrix is a lattice Green’s function connecting a pair of the sb sites. The resulting matrix is quite stable with respect to numerical inversion. The technique has been used to solve boundary value problems [36] and transport in a potential field [35]. It is a powerful technique and can be used practically in field applications because the block size can be much bigger than one usually associated with a large numerical computation because the wðs; tÞ already includes much of the smaller scale dynamics. For example, a good accounting of the data in an extensive fieldscale tracer test with the CTRW was carried out successfully up to a scale of 100 m [9]. We have discussed the physical basis of using wðs; tÞ with an algebraic tail in generic geological formations. We have analyzed the limitations of representing the w~ðuÞ by only the leading term in the small u expansion. Finally we assembled this approach into a unified framework for heterogeneous media by using an ensemble average in a component (block) in a larger explicit structural model.
Acknowledgements € sp€ A o test data were analyzed in cooperation with Golder Associates, Inc. for the Japan Nuclear Cycle Development Institute (JNC) participation € sp€ in the International A o Task Force. The authors thank SKB, Stockholm for permission to analyze these data, and Bill Dershowitz of Golder Associates Inc. for access to the data. The financial support of the Israel Science Foundation (Grant No. 147/01-1) is gratefully acknowledged.
Appendix A A.1. Advection–dispersion equation (ADE) The textbook approach to transport of passive tracers in both surface and subsurface systems usually focuses on the ADE [37] oWB =ot þ voWB =ox ¼ Ko2 WB =ox2 ;
ðA:1Þ
e.g., for 1D with constant v the average fluid velocity and K the dispersion coefficient. This equa-
tion governs the temporal evolution of the probability density function (pdf) WB ðx; tÞ of finding the tracer particle, which undergoes dispersive (or Brownian) motion, at a certain position x at time t after release at t ¼ 0. The pdf WB is the normalized concentration profile. A.2. Fractional derivative equations Mathematically, the CTRW (in the special case using wðtÞ (for 0 < b < 1Þ, and the small u expansion 1 cb ub for w~ðuÞ) is interchangeable with the time-fractional ADE [29] oW o o2 1b ¼ 0 Dt vb þ Kb 2 W ðx; tÞ; ðA:2Þ ot ox ox where the Riemann–Liouville operator is defined in terms of the convolution [29] Z t 1 d W ðx; t0 Þ 1b W ðx; tÞ ¼ dt0 ; 0 Dt 1b CðbÞ dt 0 ðt t0 Þ 0 < b < 1:
ðA:3Þ
The form of (A.2) is a natural extension of the ADE (A.1) for power-law forms of w. In this special case, the solutions of (A.2) reproduce those of the CTRW [29] in the spatial continuum limit. References [1] I. Oppenheim, K.E. Shuler, G.H. Weiss, The Master Equation, MIT Press, Cambridge, 1977. [2] M.F. Shlesinger, in: Encyclopedia of Applied Physics, vol. 16, VCH Publishers Inc., New York, 1996. [3] J. Klafter, R. Silbey, Phys. Rev. Lett. 44 (2) (1980) 55. [4] H. Scher, M.F. Shlesinger, J.T. Bendler, Phys. Today (January) (1991) 26, and references therein. [5] H. Scher, M. Lax, Phys. Rev. B 7 (10) (1973) 4491. [6] H. Scher, M. Lax, Phys. Rev. B 7 (10) (1973) 4502. [7] E. Muller-Horsche, D. Haarer, H. Scher, Phys. Rev. B 55 (1987) 1273. [8] M. Abkowitz, M.J. Rice, M. Stolka, Philos. Mag. B 61 (1990) 25. [9] B. Berkowitz, H. Scher, Phys. Rev. E 57 (5) (1998) 5858. [10] Y. Hatano, N. Hatano, Water Resour. Res. 34 (5) (1998) 1027. [11] B. Berkowitz, H. Scher, S.E. Silliman, Water Resour. Res. 36 (1) (2000) 149, Minor correction: 36 (5) 1371. [12] G. Kosakowski, B. Berkowitz, H. Scher, J. Cont. Hydrol. 47 (2001) 29.
H. Scher et al. / Chemical Physics 284 (2002) 349–359 [13] H. Scher, G. Margolin, R. Metzler, J. Klafter, B. Berkowitz, Geophys. Res. Lett. 29 (5) (2002), 10.1029/ 2001GL014123. [14] J. Eggleston, S. Rojstaczer, Water Resour. Res. 34 (9) (1998) 2155. [15] M. Blunt, H. Scher, Phys. Rev. E 52 (1995) 6387. [16] B. Berkowitz, J. Klafter, R. Metzler, H. Scher, Water Resour. Res., in press. [17] V.M. Kenkre, E.W. Montroll, M.F. Shlesinger, J. Stat. Phys. 9 (1) (1973) 45. [18] M.F. Shlesinger, J. Stat. Phys. 10 (5) (1974) 421. [19] G. Margolin, B. Berkowitz, J. Phys. Chem. B 104 (16) (2000) 3942, Minor correction: 104 (36) 8762. [20] G. Margolin, B. Berkowitz, Phys. Rev. E 65 (2002) 031101. [21] E. Bonnet, O. Bour, N.E. Odling, P. Davy, I. Main, P. Cowie, B. Berkowitz, Rev. Geophys. 39 (3) (2001) 347. [22] C.F. Tsang, Y.W. Tsang, F.V. Hale, Water Resour. Res. 27 (12) (1991) 3095. [23] H. Abelin, L. Bigersson, L. Moreno, H. Widen, T. Agren, I. Neretnieks, Water Resour. Res. 27 (12) (1991) 3119. [24] J. Hadermann, W. Heer, J. Cont. Hydrol. 21 (1–4) (1996) 87. [25] C.R. Sidle, B. Nilsson, M. Hansen, J. Fredericia, Water Resour. Res. 34 (10) (1998) 2515.
359
[26] A. Bellin, P. Salandin, A. Rinaldo, Water Resour. Res. 28 (9) (1992) 2211. [27] G. Margolin, B. Berkowitz, H. Scher, Water Resour. Res. 34 (9) (1998) 2103. [28] G. Margolin, B. Berkowitz, Phys. Rev. E, submitted. [29] R. Metzler, J. Klafter, Phys. Rep. 339 (1) (2000) 1, and references therein. [30] R. Metzler, T.F. Nonnenmacher, Phys. Rev. E 57 (6) (1998) 6409. [31] R. Hilfer, L. Anton, Phys. Rev. E 51 (2) (1995) R848. [32] E.M. Labolle, G.E. Fogg, Transp. Porous Media 42 (2001) 155. [33] C.E. Koltermann, S.M. Gorelick, Water Resour. Res. 32 (9) (1996) 2617. [34] C.E. Feehley, C. Zheng, F.J. Molz, Water Resour. Res. 36 (9) (2000) 2501. [35] H. Scher, S. Rackovsky, J. Phys. Chem. 81 (4) (1984) 1994. [36] H. Scher, E.W. Montroll, Phys. Rev. B 12 (6) (1975) 2455. [37] G. Dagan, S.P. Neuman (Eds.), Subsurface Flow and Transport: A Stochastic Approach, Cambridge University Press, Cambridge, UK, 1997.