Upscaling and reversibility of Taylor dispersion in heterogeneous porous media

Upscaling and reversibility of Taylor dispersion in heterogeneous porous media

C. W. J. Berentsen,*M. L. Verlaan,†and C. P. J. W. van Kruijsdijk‡

Department of Geotechnology, Delft University of Technology, Mijnbouwstraat 120, Delft, The Netherlands

共Received 9 July 2004; revised manuscript received 17 November 2004; published 29 April 2005兲

Tracer flow in stratified porous media is dominated by the interaction between convective transport and transverse diffusive mixing. By averaging the tracer concentration in the transverse direction, a one-dimensional non-Fickian dispersion model is derived. The model accounts for the relaxation process that reduces the convective transport to dispersive mixing. This process is 共short-兲 time correlated and partially reversible upon reversal of flow direction. For multiscale velocity fields, the relaxation is a multiscale process. To date only single scale processes have been successfully upscaled. Our procedure extends this to multiscale processes, using scale separation. The model parameters can be calculated a priori based on the velocity profile. For periodic flow reversal, the results are essentially the same. Despite the non-Fickian behavior during a cycle, the net contribution of each cycle to the spreading relaxes to a Fickian process in a similar way as for unidirectional flow. The cycle time averaged dispersion coefficient is a monotonically increasing function of the reversal time. It asymptotically converges towards the effective dispersion coefficient in the absence of flow reversal.

DOI: 10.1103/PhysRevE.71.046308 PACS number共s兲: 47.55.Mh, 47.55.Hd, 05.60.⫺k, 05.70.Ln

I. INTRODUCTION

Considering convective dispersion, there are a number of situations where the physical mechanisms are not well de-scribed by the classical convection-dispersion equation 共CCDE兲.

共i兲 In natural rocks, often several scales of heterogeneities 共e.g., layers兲 are encountered which exhibit a different dis-persion behavior from those that can be obtained from the CCDE关1,2兴.

共ii兲 The interaction between convection and molecular dif-fusion can be significant关3–8兴.

共iii兲 When the flow direction is reversed, convective spreading will also reverse. The CCDE describes the disper-sion as an irreversible mixing mechanism关9,10兴.

Numerous researchers modified the CCDE to fit their problem. The modifications range from extending the model for the macroscopic dispersion coefficient关1,11,12兴 to add-ing extra terms to the CCDE关2,13,14兴. However, none of the extended models is capable of handling all three situations outlined above. This paper describes a method for upscaling of dispersion in arbitrary layered porous media to obtain a macroscopically averaged model which incorporates the above-mentioned physical mechanisms. It extends the ap-proach taken by Camacho关15,16兴. The emphasis of this pa-per is on the validity of the upscaled model and the relation between the parameters of the model and the physics of the process.

II. PHYSICS—DISPERSION IN A UNIDIRECTIONAL FLOW FIELD

Consider a macroscopically stationary unidirectional ve-locity field v共y兲, which corresponds to flow through an infi-nite two-dimensional porous media bounded in the trans-verse or y direction, see Fig. 1. At t = 0, we inject a finite mass 共M0兲 of ideal 关24兴 tracer 共non兲uniformly distributed over the height共y兲 at longitudinal position x0. In addition to convective transport, the tracer particles are subject to a mi-croscale mixing mechanism, which we will limit in this pa-per to isotropic molecular diffusion 关25兴. The evolution of the 共tracer兲 concentration can be described, at the macros-cale, by the two-dimensional unidirectional convection-diffusion equation共2D uCDifE兲,

⳵c共x,y,t兲

⳵t +v共y兲

⳵c共x,y,t兲

⳵x = Dmol⌬c共x,y,t兲. 共1兲 Our objective is to describe the evolution of the transverse averaged concentration c0共x,t兲=共1/d兲兰0

d_{c共x,y,t兲dy.}

A. Convective displacement

First consider convective transport only. Each tracer par-ticle moves with a constant velocity in the x direction along its initial streamline. Hence, the mean position of all tracer particles moves linearly in time with a constant averaged particle velocity共u0兲. This velocity may differ from the mean fluid velocity共v0兲 if the particles are nonuniformly distrib-uted in the transverse direction. Similarly, the spreading of the particles in the longitudinal direction, expressed in terms of the variance, increases with time squared and is propor-tional to the variance of the tracer velocity determined by those streamlines occupied by tracer particles.

*Present address: Department of Earth Sciences, Utrecht Univer-sity. Electronic address: berentsen@geo.uu.nl

†

Present address: Horizon Energy Partners. Electronic address: marco.verlaan@horizon-ep.com

‡

Electronic address: c.p.j.w.vankruijsdijk@citg.tudelft.nl

B. Relaxation from convective towards uncorrelated-Fickian behavior

In addition to convection, consider mixing by molecular diffusion. The longitudinal component of diffusion smoothes the sharp front propagated by convection. More importantly, in the transverse direction diffusion moves particles away from their initial streamlines and the particle velocity will change. Over time, each tracer particle will sample each ver-tical position equally frequently, and the relation between a particle and its initial streamline is lost. Consequently, diffu-sion smoothes concentration differences in the transverse di-rection, and as a result reduces the longitudinal spreading caused by the velocity profile. In the long-time limit, a dy-namic equilibrium between the convective and diffusive spreading mechanisms is established, which exhibits Fickian behavior. Figure 2 visualizes this process for the velocity field shown in Fig. 1. The relaxation time␶is the character-istic time for the transition process 关26兴 also called relax-ation. This time is a function of the value of the molecular diffusion coefficient and the correlation of the longitudinal velocity in the transverse direction.

C. The effect of the initial transverse particle distribution If the particle distribution is uniform over the height, local redistribution of the particles by molecular diffusion does not affect the particle distribution over the streamlines. The mean particle velocity equals the mean fluid velocity at all times. If the initial particle distribution is not uniform, transverse dif-fusion will establish a uniform distribution over the height in time and the mean particle velocity evolves from its initial value u0to the mean fluid velocity v0.

III. SPECTRAL EQUIVALENT

Similar to Camacho关16兴, the starting point for the deri-vation of an evolution equation for the height averaged con-centration 关c0共x,t兲兴 is the spectral equivalent of the 2D uCDifE. We replace the concentration and velocity 共1兲 by cosine Fourier series, multiply the result by cos共n␲y / d兲, and integrate it over y. For n = 0, we obtain the evolution equation of c0共x,t兲, ⳵c₀ ⳵t +v0 ⳵c₀ ⳵x − Dmol ⳵2_c 0 ⳵x2 = − ⳵J_T ⳵x, 共2兲

where, following Camacho关16兴, the single valued dispersive Taylor flux JTis defined as

JT共x,t兲 = 1 2

兺

n=1

⬁

vncn共x,t兲. 共3兲 It originates from the convective term in Eq.共1兲 and ac-counts for the contribution of the higher-order modes to c0共x,t兲. The evolution of the higher-order modes 共n⬎0兲 is given by ⳵cn ⳵t + cn ␶n +v0 ⳵cn ⳵x +_m=1

兺

⬁ ␥mn vn vm ⳵cm ⳵x − Dmol ⳵2_c n ⳵x2 = −vn ⳵c0 ⳵x, 共4兲 where the modal relaxation time ␶n and modal interaction coefficient␥mn are defined as

FIG. 1. Layered porous media in the domain⍀=共−⬁, +⬁兲⫻关0,d兴 in which a longitudinal velocity field v共y兲 is present. Tracer particles are initially released in this field at longitudinal position x0uniformly or nonuniformly distributed in the transverse direction. Tracer particles are subject to isotropic molecular diffusion.

FIG. 2. Typical evolution of particle clouds in the共x,y兲 space when the particles are initially uniformly distributed over the transverse direction. Clouds at increasing position共x兲 represent clouds for increasing time. Times 共t兲 are in days. The molecular diffusion is Dmol = 1.25⫻10−4_m2_{/ d.}

␶n= d2 共n␲兲2_D mol 共5兲 and ␥mn=

冦

1 2 vn vm 关vm+n兴, m = n 1 2 vn vm 关vm+n+v兩m−n兩兴, m ⫽ n.

冧

共6兲

IV. MOMENT ANALYSIS FOR UNIDIRECTIONAL FLOW

The spectral representation also serves as a starting point for the derivation of the spatial moments belonging to the 2D uCDifE. As the 2D model is generally too involved to allow direct analytical solution, moment analysis allows us to quantify the accuracy of an upscaled model. We briefly dis-cuss the behavior of the zeroth, first, and second spatial mo-ments belonging to the 2D uCDifE. For a full analysis, see 关5,17,18兴.

A. General definition of the spatial moments

Consider the release of a tracer with finite mass M0in the spatial domain x苸共−⬁, +⬁兲. The kth spatial moment of the nth concentration mode is defined as

E_c n,x,k=

冕

−⬁ +⬁ xk_c n共x兲dx. 共7兲

The zeroth moments or “mass” in the nth concentration mode are denoted by M_c

n. Mcn,x,kis the tracer mass

normal-ized kth moment of the nth concentration mode,共Ec_n,x,k/ M0兲. M_c

n,x,k

c _{, for k⬎1, is the centerd normalized kth moment} 关Mc_n,_{共x−Mc,x,1兲,k}兴. For simplicity, the mean particle position 共Mc₀,x,1兲 is written as␮c,xand the particle variance共Mc₀,x,2

c _兲

as␴_c,x2 .

B. The evolution of the spatial moments of the spectral equivalent

To obtain the moment equations belonging to the 2D uCDifE, we multiply Eqs.共2兲 and 共4兲 of the spectral equiva-lent with xk and average the result over x. Assuming that each concentration mode and all its spatial derivatives con-verge exponentially to zero for x→ ±⬁, the evolution equa-tion for the kth normalized moment belonging to c₀共x,t兲 can be written as ⳵Mc₀,x,k ⳵t − kv0Mc0,x,k−1−共k 2_{− k兲D} molMc₀,x,k−2 =k 2

兺

_n=1 ⬁ vn Ec_n,x,k−1 M_c 0 . 共8兲

The evolution of the kth non-normalized moment of a higher concentration mode共n⬎0兲 reads

冋

␶n −1₊ ⳵ ⳵t

册

Ecn,x,k= kvnMc0Mc0,x,k−1+ k关2v0+v2n兴Ecn,x,k−1 +共k2− k兲DmolEc_n,x,k−2 + k 2_m=1

兺

_⫽n ⬁ 关vm+n+v兩m−n兩兴Ec_n,x,k−1. 共9兲

C. The transverse particle distribution

If the tracer is initially nonuniformly distributed over the height, the initial zeroth moment共“mass”兲 of a higher con-centration mode cn共n⬎0兲 may be nonzero,

M_c n init = 2

冕

−⬁ +⬁

冕

0 d c共x,y,0兲cos

冉

n␲y d

冊

dydx. 共10兲 Solving Eq. 共9兲 for k=0 with respect to initial condition 共10兲 gives the evolution of the zeroth moment of the higher-order modes共n⬎0兲,

Mc_n共t兲 = Mc_n init

e−t/␶n_. 共11兲

Equation共11兲 shows that the zeroth moment of c_nrelaxes to zero at a speed which only depends on its own relaxation time␶n. Its relaxation is not affected by other concentration modes.

D. The mean particle position and velocity

If the particles are distributed nonuniformly in the trans-verse direction, the mean particle velocity 关u共t兲兴 shows a relaxation from the initial mean particle velocity关u0= u共0兲兴, determined by the initially occupied streamlines, towards the mean fluid velocity共v0兲,

u共t兲 =

冕

−⬁ +⬁

冕

0 d vcdydx

冕

−⬁ +_⬁

冕

0 d cdydx =v0+ 1 2

兺

_n=1 ⬁ vn M_c n init M0 e−t/␶_{n. 共12兲}

Convergence from u0to v0is faster for increasing magni-tude of the molecular diffusion Dmol, as shown in Fig. 3.

The exact mean position of the particles is obtained from Eq.共8兲 with k=1 or by integration of the mean particle ve-locity共12兲 over time,

␮c,x=关x0+v0t兴 + 1 2

兺

_n=1 ⬁ ␶nvn Mc_n init M₀ 关1 − e −t/␶_{n兴. 共13兲} In addition to a contribution present for initial uniform particle distributions 共x0+v0t兲, a relaxation term expresses the contribution due to the redistribution of the particles over the streamlines. After relaxation, the mean moves with the mean fluid velocity as in the uniform case, but shows a de-viation⌬␮_relaxresulting from the relaxation of the mean par-ticle velocity,

␮c,x,disp= x0+v0t +

兺

n=1 ⬁ ␶nvn 2 Mc_n init M₀ = x0+v0t +⌬␮relax. 共14兲

E. The variance of the averaged concentration

The spatial variance is given by 共for derivation, see Ap-pendix A 2兲 ␴c,x 2 _{共t兲 = 共2D} molt兲 +

冉

兺

n=1 ⬁ vn 2_␶ n共t +␶n关e−t/␶n− 1兴兲

冊

−

冉

1 2

兺

_n=1 ⬁ _M c_n init M0 vn␶n关1 − e−t/␶n兴

冊

2 +

冉

1 2

兺

_n=1 ⬁ _M c_n init M₀ vnv2n共␶n 2 −关␶_n2+␶nt兴e−t/␶n兲 +1 2_n⫽m=1

兺

⬁ _M c_m init M0 vn关vm+n+v兩m−n兩兴 ␶n␶m ␶m−␶n ⫻ 共␶m关1 − e−t/␶m兴 −␶n关1 − e−t/␶n兴兲

冊

共15兲 The variance has four contributions. The first is by the longitudinal component of diffusion. The second expresses the interaction of the velocity field and transverse molecular diffusion. The third and fourth terms correct the second term when the particle distribution is initially nonuniform. They describe a contribution of the changing averaged particle ve-locity resulting from the transverse redistribution. The third term is a direct result of the relaxation of the mean共second term of Eq.共13兲兲. The fourth term accounts for interactions of the higher-order modes. The temporal particle velocity variance is given by ␴u 2_{共t兲 =} 1 M0

冕

−⬁ ⬁

冕

0 d

关v共y兲 − u共t兲兴2_{c共x,y,t兲dydx =␴}

v 2_+⌬␴ u 2_共t兲 =

兺

n=1 ⬁ vn 2 2 +

冋

−

冉

兺

_n=1 ⬁ vn 2 M_c n init M0 e−t/␶n

冊

2 +

兺

n=1 ⬁ vnv2n 4 Mc_n init M₀ e −t/␶_n +

兺

m⫽n=1 ⬁ _关v m+n+v兩m−n兩兴vm 4 M_c n init M0 e−t/␶n

册

_. 共16兲 ⌬␴u

2_{共t兲 is the deviation of the particle velocity variance} from the fluid velocity variance共␴_v2兲 for transverse nonuni-form particle distributions. For short times共tⰆ␶兲, the part of the spreading共15兲 induced by convection is proportional to time squared and to the variance of the initial particle veloci-ties,

␴c,x,conv

2 _{共t兲 = 2D} molt +␴u

2_共0兲t2_{. 共17兲} For uniform particle distributions, two velocity fields with the same velocity variance but with different modal compo-sition produce at early times the same amount of spreading. In the long-time limit共tⰇ␶兲, the variance behaves Fick-ian and is proportional to t,

␴c,x,disp 2 _{共t兲 = 2D} ef f,_⬁t +⌬␴uni 2 _+⌬␴ non 2 _{. 共18兲}

As ultimately the particle distribution becomes uniform over the height, only the first two terms of Eq.共15兲 contrib-ute to the effective asymptotic dispersion coefficient共Def f,⬁兲,

Def f,⬁= Dmol+ 1 2

兺

_n=1 ⬁ vn 2_␶ n. 共19兲

For fields with the same velocity variance, those domi-nated by lower-order modes 共vn兲 show ultimately more spreading than those dominated by higher-order modes. In Eq. 共18兲, ⌬␴_uni2 is the constant contribution to the variance resulting from the relaxation of the uniform part of the vari-ance关second term in Eq. 共15兲兴 defined as

⌬␴uni 2 = −

兺

n=1 ⬁ v_n2␶_n2. 共20兲 ⌬␴non 2

is a constant contribution to the variance that re-sults from the relaxation of the nonuniform part 关third and fourth terms of Eq.共16兲兴 given by

⌬␴non 2 _{= +}

兺

m⫽n=1 ⬁ _M c_m init M0 vn关vm+n+v兩m−n兩兴␶n␶m 2 −

冉

兺

n=1 ⬁ _M c_n init M0 vn␶n 2

冊

2 +

兺

n=1 ⬁ _M c_n init M0 vnv2n␶n␶2n 2 . 共21兲 Figure 4 shows the variance for a uniform initial particle distribution. It shows a relaxation from a t2 _{behavior for} short times towards Fickian behavior 共⬃t兲 for long times.

FIG. 3. Typical example of the relaxation of the mean particle velocity towards the mean fluid velocity v0for three values of the molecular diffusion.

The shift in the variance关⌬␴_uni2 , Eq.共20兲兴 resulting from the relaxation共Fig. 4兲 is negative as the variance time derivative monotonously increases towards its asymptotic value 共2Def f,⬁兲.

V. EVALUATION OF THE UPSCALED RELAXATION MODEL BY CAMACHO

The laminar flow problem studied by Camacho关16兴 ex-hibits a single dominant Fourier scale. He derived an evolu-tion equaevolu-tion for the Taylor flux by summing the relaxaevolu-tion of the modal Taylor fluxes共Jn=vncn/ 2兲 关27兴,

⳵JT ⳵t + JT ␶ +v0 ⳵JT ⳵x +␤ ⳵JT ⳵x − Dmol ⳵2_J T ⳵x2 = −␴v 2⳵c0 ⳵x. 共22兲 Subsequently, he replaced␶共x,t兲 and␤共x,t兲 by the con-stants ␶ef f and␤ef f yielding an effective single scale relax-ation. Combination of Eqs.共2兲 and 共22兲 yields a fourth-order PDE for c0共x,t兲. ⳵c0 ⳵t +v0 ⳵c0 ⳵x − Dmol ⳵2_c 0 ⳵x2 +␶ef f

冋

⳵2_c 0 ⳵t2 +关共v0+␤ef f兲v0−␴v 2_兴⳵ 2_c 0 ⳵x2 +共2v0+␤ef f兲 ⳵2_c 0 ⳵x⳵t

册

=␶ef f

冋

2Dmol ⳵3_c 0 ⳵x2⳵t+共2v0+␤ef f兲Dmol ⳵3_c 0 ⳵x3 − Dmol 2 ⳵ 4_c 0 ⳵x4

册

. 共23兲 Camacho关16兴 effectively reduced the general multiscale problem to a single scale problem. We investigate the devia-tions from the full problem by comparing the mean and vari-ance of the approximation with the exact results.

A. Comparison of the mean particle position

The expression for the mean particle position using the upscaled equation 共23兲 is similar to the full solution 关Eq. 共13兲兴 but restricted to a single relaxation scale 共see Appendix B 1兲,

␮c,x ap_{= x}

0+v0t +␶ef f共u0−v0兲共1 − e−t/␶ef f兲. 共24兲 Equation共24兲 only matches the mean of the full model for single scale problems or for initially uniform distributions. The convective limit is always correctly described as it fol-lows from the initial conditions. The constant contribution to the Fickian limit originating from the relaxation process can be made to match the exact value by defining the effective relaxation time as ␶ef f=

兺

n=1 ⬁ ␶nvn M_c n init M0 2共u0−v0兲 =

兺

n=1 ⬁ ␶nvn M_c n init M₀

兺

n=1 ⬁ vn M_c n init M0 . 共25兲

Note that this approach fails when the initial distribution is nonuniform and u0= v0.

B. Comparison of the variance

The variance calculated from the upscaled equations has again a structure similar to the full expression关Eq. 共15兲兴 but restricted to a single scale共see Appendix B 2兲,

␴c,x ap,2

= 2Dmolt + 2␴v

2

␶eff关␶eff共e−t/␶eff− 1兲 + t兴 + 2␶eff␤eff共u0−v0兲 ⫻关␶eff共e−t/␶eff− 1兲 + te−t/␶eff兴 − 兵␶eff共u0−v0兲关e−t/␶eff

− 1兴其2. 共26兲

For uniform initial distributions, the third and fourth terms on the right-hand side disappear. As only a single ef-fective relaxation process is present in the approximation, it is only exact for single scale problems. The short-time limit is again exact. In the long-time limit共tⰇ␶兲, the approxima-tion reduces to ␴c,x,disp ap,2 = 2共Dmol+␴v 2_␶ eff兲t. 共27兲

If we set the effective relaxation time共␶ef f兲 equal to the so-called standard relaxation time␶std共see Camacho 关16兴兲,

␶ef f=␶std⬅ 1 2␴_v2

兺

n=1 ⬁ v_n2␶n=

兺

n=1 ⬁ v_n2␶n

兺

n=1 ⬁ vn 2 , 共28兲

we also match the exact variance in the long-time limit共18兲. However, Sec. X will show that this definition does not give an accurate approximation for intermediate times. For an ini-tially nonuniform distribution, the situation is more complex. In the convective limit, the upscaled variance reads

␴c,x,conv

ap,2

= 2Dmolt + 2␴v

2

t2−␤eff共u0−v0兲t2−共u0−v0兲2t2. 共29兲 This expression does not have any information on the initial particle velocity variance which drives the variance evolution in the convective limit. However, we can define parameter ␤_eff to introduce this information in such a way that the short-time limit is again exact,

␤ef f= ␴v 2 −␴_u,init2 −共u0−v0兲2 u0−v0 . 共30兲

FIG. 4. Relaxation of the variance from convective共⬃t2_{兲 for} short times to Fickian behavior共⬃t兲 for large times.

Note that as observed before in the discussion of the mean, in the case u0= v0 the approach above fails. The up-scaled equations simply cannot handle this situation. In the long-time limit共tⰇ␶兲, the variance shows the correct Fick-ian behavior, ␴c,x,disp ap,2 = 2共Dmol+␴v 2 ␶eff兲t − 关2␴v 2 ␶eff 2

+ 2␶eff2␤eff共u0−v0兲 +␶eff2共u0−v0兲2兴. 共31兲 Ultimately, the growth of the variance is independent of the initial distribution. The relaxation constant, however, is different, since the relaxation process has evolved in a dif-ferent way. To match the approximation with the exact vari-ance in the Fickian limit, the uniform definition for the re-laxation time can be applied. The constant contribution due to the relaxation process, however, cannot be matched.

C. Relaxation time considerations

For nonuniform initial conditions, the effective relaxation time has to be tuned to the initial distribution共25兲 to obtain a match for the mean for long times. Turning to the effective dispersion coefficient in the Fickian limit, we can match the upscaled expression to the exact relationship by an appropri-ate definition of the effective relaxation time共28兲. However, as may be clear from the physics, this definition is based on the fluid velocity profile and is unrelated to the initial distri-bution. The upscaled model can only match one of the two conditions.

VI. GENERALIZED TELEGRAPH EQUATION

The fourth-order upscaled PDE Eq.共23兲 covers the whole spectrum from cases where 共longitudinal兲 molecular diffu-sion dominates the spreading to cases where the transverse variation in the velocity field dominates the spreading. If the contribution by longitudinal molecular diffusion to the be-havior of the fourth-order model can be neglected, the para-bolic fourth-order upscaled equation reduces to a second-order hyperbolic generalized Telegraph equation,

␶ef f ⳵2_c 0 ⳵t2 + ⳵c0 ⳵t +v0 ⳵c0 ⳵x +␶ef f共2v0+␤ef f兲 ⳵2_c 0 ⳵x⳵t −␶_{ef f}关␴_v2−共v₀+␤_{ef f}兲v₀兴⳵ 2_c 0 ⳵x2 = 0. 共32兲

This equation describes the same mean and variance as the fourth-order model if the contribution of the spreading caused by longitudinal molecular diffusion is ignored.

A. Concentration profiles

Unfortunately, the fourth-order upscaled PDE Eq. 共23兲 does not have an easily analyzable analytical solution. In contrast, for initial conditions c共x,y,0兲=H共x−x₀兲 and

⳵tc0共x,0兲=0, the solution to the Telegraph model Eq. 共32兲 reads for x苸共−⬁,⬁兲, see 关18兴,

c共x,t兲 =

冦

H共␹−⌫兲

冋

De−␹I0兵

冑

␹2−⌫2其 + e−⌫ 2 + ⌫ 2

冕

_⌫ ␹ _e−␰˜

冑

˜␰2_−⌫2 I1兵

冑

˜␰2−⌫2其d˜␰

册

, x⬎ v0t 1 + H共␹−⌫兲

冋

De−␹I0兵

冑

␹2−⌫2其 − e−⌫ 2 − ⌫ 2

冕

_⌫ ␹ _e−␰˜

冑

˜␰2_−⌫2 I1兵

冑

˜␰2−⌫2其d˜␰

册

, x⬍ v0t

冧

共33兲 with ␹= 2␺ ␧2_{+ 4␺␸}

冉

t

冋

1 − v0␧ 2␺

册

+ x ␧ 2␺

冊

, ⌫ = 兩x − v0t兩

冑

␧2_{+ 4␺␸}, D = −␧ 2

冑

␧2+ 4␺␸, 共34兲 and

␸=␶eff, ␧ = +␤eff␶eff, ␺=␴v

2

␶eff. 共35兲 For the velocity field shown in Fig. 1, we compare the concentration profiles by the fourth-order model 共23兲, the Telegraph equation共32兲, and the height averaged concentra-tion profile by the full 2D model共1兲. Both upscaled models

show nearly identical concentration profiles during the relax-ation from convective to Fickian behavior关Figs. 5, 6共a兲, and 6共b兲兴, even for relatively large values for the molecular dif-fusion 关Fig. 6共b兲兴. For short times, the hyperbolic part 共␶ef f⳵t

2

c0兲 in the Telegraph model dominates over the para-bolic part 共⳵c₀/⳵t兲 and its concentration profile shows two shocks, typical for a second-order hyperbolic equation. Even though the profiles of the upscaled models seem to give a poor representation of the exact profile 关Fig. 6共a兲兴, they match the exact solution up to the second moment. In the long-term limit, both upscaled models demonstrate Fickian behavior, characterized by the typical S-shaped concentration profile, and give an excellent match with the exact solution 关Fig. 6共b兲兴. Note that the propagation speed of the Telegraph model is constant and finite, contrary to the fourth-order model and classical convection diffusion equation关28兴.

VII. PHYSICS—REVERSED FLOW

Now consider the case where after a time t = t_rwe reverse the flow direction. Assuming convective transport only, the process is fully reversible and the particles are transported back to their original positions, see Fig. 7共a兲.

A. Reversal of flow and relaxation

If we reverse the flow and include molecular diffusion, the relaxation process starts all over again. However, the particle distribution at the time of reversal is shaped by the convection dispersion process. As discussed above, each par-ticle velocity in this process is correlated over a short time.

As a result, initially particles return longitudinally along the same path as they arrived from. Hence the transport process demonstrates a 共partially兲 reversible behavior. The variance of the height-averaged particle positions decreases as the par-ticles turn back along their original paths. Due to transverse diffusion, the velocity of a particle becomes in time less correlated to the velocity history of its forward movement. As a result, convective dispersion takes over and the variance once again increases monotonously. The process is identical to the original relaxation process. Consequently, the same relaxation time characterizes the interaction process of the reversed velocity field and molecular diffusion to form a dy-namic equilibrium exhibiting a Fickian behavior. Moreover, this共relaxation兲 time is independent of the moment in time at which we reverse the flow.

Consider a single scale relaxation process with relaxation time␶. We distinguish three cases based on the ratio of di-mensionless reversal time, trD= tr/␶. For trDⰆ1, the particles exhibit a fully correlated behavior similar to pure convection 关Fig. 7共a兲兴. In contrast, if trDⰇ1, the transport process is fully relaxed at the time of flow reversal 关Fig. 7共c兲兴. For intermediate times trD, the particles, at the time of flow re-versal, will be significantly correlated to their original streamlines and each particle velocity, at the time of reversal, will be correlated to the initial particle velocity in time. Con-sequently, the velocity profile may be clearly visible in the particle cloud关Fig. 7共b兲兴. Although the reversibility may be much more significant and visible for small trD⬍1, partial reversibility, however small, is observed for all values of trD. In general, the relaxation process may evolve over multiple scales. At times smaller than the relaxation time of the larg-est scale, the various modes will generally not be in the same state of relaxation.

B. Evolution of the variance

If we reverse the flow, two combined processes contribute to the variance. The first process expresses the correlation between forward and reversed velocities of a particle and causes a decrease in the variance. The second is a new re-laxation process identical to unidirectional rere-laxation, but with the particle distribution at reversal as initial distribution. Initially, the first process dominates and the reversed velocity of a particle is completely correlated to its last forward ve-locity. Consequently the time derivative of the variance changes sign关Fig. 8共b兲兴, which we experience, even in the

FIG. 5. Evolution of the concentration profiles of Camacho’s fourth-order model and the gener-alized Telegraph equation 共for D_mol= 1.25

⫻10−4_m2_d−1_{兲. For t=500 d, the profile obtained} from a full 2D random-walk simulation is shown for comparison.

FIG. 6. Comparison of the concentration profiles of the full 2D, Camacho’s fourth-order, and the Telegraph model at time t = 25 d.

共a兲 The concentration profile in the convective limit for Dmol= 2.5

⫻10−9_m2_d−1_{.共b兲 The concentration profile in the long-time} Fick-ian limit for D_mol= 6.25⫻10−4_m2_d−1_{. Note that the共Gibbs兲} over-shoots in the fourth-order model are caused by the third-order up-wind scheme that is used for the spatial discretization.

Fickian limit, as ”partial” reversibility. With time a particle velocity loses its correlation with the forward velocities, starting with the short-time correlations共or small scales兲. At a certain moment, the two processes counterbalance each other, marked by the point at which the variance reaches its minimum关Fig. 8共a兲兴 and its time derivative passes zero 关Fig. 8共b兲兴. From this point on, the second process dominates and relaxes to the same Fickian behavior as for unidirectional flow while the first process relaxes to a constant value. In the convective limit, the velocity of each particle is fully corre-lated共or constant兲 in time and the first process dominates up to the point that the time after reversal equals the reversal time关Fig. 8共a兲,␶= 25 000兴.

VIII. MOMENT ANALYSIS FOR FLOW INCLUDING REVERSAL OF DIRECTION

Here we analyze the spatial moments belonging to the 2D uCDifE for flow including reversal. The evolution expres-sions for the kth moments关29兴 M_c

0,x,k

rev _{共t兲 and E}

c_n,x,k rev _{共t兲 after}

flow reversal are similar to expressions共8兲 and 共9兲, respec-tively, with the signs of the modal velocities changed.

A. Transverse particle distribution and mean particle position Flow reversal does not affect the way transverse diffusion redistributes the particles in the vertical direction. Hence, flow reversal does not affect the zeroth moment of the

higher-order modes and is given by Eq.共11兲. The mean par-ticle velocity is directly related to the transverse distribution and only changes sign when the flow direction is reversed. The mean particle position for t⬎tr reads

␮c,x rev_{共t兲 = x} 0+v0共2tr− t兲 + 1 2

兺

_n=1 ⬁ ␶nvn M_c n init M0 ⫻关e−t/␶_{n− 2e}−t_r/␶_{n+ 1兴, t ⬎ t} r. 共36兲 In the convective limit, the mean is fully reversible,

␮c,x,conv

rev _{= x}

0+ u0共2tr− t兲. In the Fickian limit 共t⬎␶兲, the verti-cal particle distribution is degenerated to a uniform distribu-tion and the mean particle posidistribu-tion reads

␮c,x,disp rev _{共t兲 =␮} c,x,uni rev _{共t兲 + ⌬␮} relax rev =关x0+v0共2tr− t兲兴 +

冋

兺

n=1 ⬁ ␶nvn 2 M_c n init M0 共1 − 2e−t_r/␶_n兲

册

_{. 共37兲} The deviation of the particle mean from the fluid mean, expressed by the relaxation constant⌬␮_relaxrev , is a function of the reversal time tr and converges to⌬␮relax共14兲 for trⰇ␶.

B. The variance of the averaged concentration The exact variance for uniform initial distributions is given by共see appendix A 2兲

FIG. 7. Particle clouds belonging to the five-layer velocity field共Fig. 1兲 at the moment of flow reversal 共t=tr兲 and at time 共t=2tr兲 for three different relaxation regimes.共a兲 The convective limit 共Dmol= 0 m2d−1兲. 共b兲 Intermediate situation 共Dmol= 3⫻10−5m2d−1兲. 共c兲 The Fickian limit共D_mol= 6.3⫻10−4_m2_d−1_兲.

␴c,x,rev,uni 2 _{共t兲 = 2D} molt +

兺

n=1 ⬁ ␶nvn 2_⫻

关

_关t r+␶n共e−tr/␶n− 1兲兴 +关共t − tr兲 +␶n共e−共t−tr兲/␶n− 1兲兴

+关␶n共etr/␶n− 1兲e−t/␶n+␶n共e−tr/␶n− 1兲兴

兴

. 共38兲 The contribution by the longitudinal component of mo-lecular diffusion is independent of the flow direction. The contribution by convection and transverse diffusion has three subterms. The first expresses the value of the variance at the time of reversal. The second represents a new relaxation pro-cess, identical to the initial one. The third term is a ”demix-ing” or reversibility term expressing the effect of flow rever-sal on the variance induced by the forward convection-dispersion process. It is increasingly negative for t⬎tr and asymptotically reaches a constant once correlation with the flow before reversal is lost. For an arbitrary initial distribu-tion, the exact variance after flow reversal reads共Appendix A 2兲 ␴c,x,rev 2 _{共t兲 =␴} c,x,rev,uni 2 _{共t兲 −}

冉

1 2

兺

_n=1 ⬁ ␶nvn Mc_n init M₀ 关e −t/␶_{n− 2e}−t_r/␶_n + 1兴

冊

2 +1 2

兺

_n=1 ⬁ _M c_n init M0 vnv2n共2tr␶n关e −t/␶_{n− e}−t_r/␶_n兴 − t␶ne−t/␶n+␶n 2_{关1 − e}−t/␶_n兴兲 +1 2_m,n=1

兺

_共m⫽n兲 ⬁ _M c_m init M0 vn关vm+n+v_兩m−n兩兴 ␶n␶m ␶m−␶n ⫻ 共␶m关1 − e−t/␶m兴 −␶m+ 2␶n关e−tr/␶n− e−t␶/␶m兴 +␶_net␶/␶_n关2e−t_␶/␶_{m− e}−t_r/␶_n兴e−t/␶_{n兲 共39兲} The additional contribution to the variance only depends on the reversal time t_rif the distribution has not yet relaxed to a uniform distribution before reversal. Taking the time derivative of Eq.共39兲 shows that for any value of␶, imme-diately upon flow reversal the dispersive contribution to the variance is negative, lim t↓t_r

冉

⳵␴c,x,rev 2 ⳵t − 2Dmol

冊

= − limt↑t_r

冉

⳵␴c,x 2 ⳵t − 2Dmol

冊

. 共40兲 In the convective limit共t,t_rⰆ␶兲, no relaxation takes place and the variance by the convective-dispersive part displays fully reversible behavior, reducing to zero at t = 2tr,

␴c,x,rev,conv

2 _{共t兲 =␴} u

2_{共0兲共t − 2t}

r兲2, 共41兲

where␴_u2is given by Eq.共16兲. In the Fickian limit 共tⰇ␶兲, the variance takes the form

␴c,x,rev,disp

2 _{共t兲 = 2D}

ef f,⬁t +⌬␴uni,rev

2

+⌬␴_non,re2 _v. 共42兲 The asymptotic dispersion coefficient Def f,⬁ is indepen-dent of the flow direction and is given by Eq. 共19兲. The constant contribution to the variance resulting from the re-laxation process for initially uniform particle distribution, ⌬␴uni,rev

2 _{, varies with the reversal time,}

⌬␴uni,rev 2 _{= −}

兺

n=1 ⬁ ␶n 2 vn 2_{共3 − 2e}−t_r/␶_{n兲. 共43兲} Its magnitude is bounded between ⌬␴_uni2 共for trⰆ␶兲 and 3⌬␴_uni2 共for trⰇ␶兲. The constant contribution to the variance resulting from the relaxation of a transverse nonuniform to a transverse uniform particle distribution⌬␴_non,re2 _v reads

⌬␴non,rev 2 _{= −}

冉

1 2

兺

_n=1 ⬁ ␶nvn M_c n init M0 关1 − 2e−t_r/␶_n兴

冊

2 +1 2

兺

_n=1 ⬁ _M c_n init M0 vnv2n␶n关␶n− 2e−tr/␶ntr兴 +1 2_m⫽n=1

兺

⬁ _M c_m init M0 vn关vm+n+v兩m−n兩兴 ⫻

冉

␶n␶m+ 2␶n e−tr/␶n_{− e}−tr/␶m ␶n−␶m

冊

共44兲 If the transverse particle distribution has relaxed to uni-formity before flow reversal,⌬␴_non,re2 _v equals⌬␴_non2 共21兲.

IX. PERIODIC FLOW REVERSAL

If we keep reversing the flow each time t has increased by tr, the time derivative of the variance asymptotically turns

FIG. 8. 共a兲 Typical example of the variance for flow including reversal of direction at time t = 250 days and共b兲 the corresponding time derivative for four values of the relaxation time ␶ =共2.5,25,250, and 25 000兲 days.

periodic as well关Fig. 9共a兲兴. Ultimately the periodic growth of the variance becomes linear in time even if the variance itself does not relax in a period tr关Fig. 9共b兲兴. Each reversal a new共but identical兲 relaxation process starts while a second process expresses the correlations with the previous reversal cycles. After sufficient time has passed, the velocity of a particle in the nth reversal cycle is no longer correlated to the velocity in the first forward period. The time to lose correla-tion is again related to the relaxacorrela-tion time. After this time has passed, a particle experiences statistically the same velocity correlation in each subsequent reversal period.

A. Mathematical derivation

Moment analysis shows that the time derivative of the variance becomes periodic when the flow is reversed periodi-cally. Since nonuniform initial distributions relax ultimately to a uniform distribution, we limit the discussion to uniform initial particle distributions共M_c

n,0

init_{= 0兲. The time derivative of} the variance in the共NR兲th reversal cycle reads 共Appendix C兲

⳵␴c,x 2 ⳵t 共t˜;Nr,tr兲 = 2Dmol+_n=1

兺

⬁ vn 2_␶ n⫻

冋

1 +共− 1兲Nre−共t˜+Nrtr兲/␶n − 2e −共t˜−t_r兲/␶_{n+共− 1兲}N_re−共t˜+N_rt_r兲/␶_n etr/␶n_{+ 1}

册

, 共45兲

where t˜ is the time since the last flow reversal. For 共Nr

→⬁兲, this expression 共45兲 converges to ⳵␴c,x,lim 2 ⳵t 共t˜;tr兲 = limNr→⬁ ⳵␴c,x 2 ⳵t 共t˜;Nr,tr兲 = 2Dmol+

兺

n=1 ⬁ vn 2_␶ n

冋

1 − 2 e−共t˜−tr兲/␶n etr/␶n_{+ 1}

册

. 共46兲 Thus even if the variance does not relax in a period t_r 共tr⬍␶兲, the growth of the variance relaxes if we keep on reversing the flow each time t has increased with tr. The net cyclic dispersion coefficient Dcyclefor Nr→⬁ follows by in-tegration of Eq.共46兲 over a cycle 共t˜苸关0,tr兴兲 and division of the result by共2tr兲, Dcycle= 1 2tr

冕

0 t_r⳵␴ c,x,lim 2 ⳵t 共t˜;tr兲dt˜ = Dmol+ 1 2

兺

_n=1 ⬁ vn 2_␶ n

冋

1 − 2 ␶n共etr/␶n− 1兲 tr共etr/␶n+ 1兲

册

. 共47兲 Figure 10 shows the normalized modal cyclic dispersion coefficient 关term between square brackets in Eq. 共47兲兴. It increases monotonously from 0 for the limit of tr/␶nto zero towards 1 for the limit of t_r/␶_nto infinity. The corresponding expression of the variance for repetitive flow reversal reads

␴c,x 2 _共t;N r,tr兲 = 2Dmolt +

兺

n=1 ⬁ vn 2_␶ n 2_⫻

冉

t ␶n + 2

兺

k=0 N_r

⬘

共− 1兲N_r−k_共e−共t−kt_r兲/␶_{n− 1兲} −

兺

k=1 N_r 共− 1兲k_关4共N r− k兲 + 2兴共e−ktr/␶n− 1兲

冊

. 共48兲 The sums in Eq. 共48兲 increase linearly with the reversal cycle 共“time”兲 when the time exceeds the relaxation time 关see also Fig. 9共b兲兴. Figure 11 demonstrates the effect of a decreasing reversal time on the variance.

X. SEPARATION OF SCALES

A. Multiple scale treatment

For multiscale velocity fields, the fourth-order model is approximate. To improve its accuracy, there are two ways to account for the changing interaction of the scales present in the full problem over time. We can optimize the effective relaxation time or apply scale separation. We discuss their relative merit as a function of scale distribution and the ob-servation time. To demonstrate this, consider a bimodal ve-locity field v共y兲,

v共y兲 = vn₁cos

冉

n1␲

d y

冊

+vn2cos

冉

n2␲

d y

冊

. 共49兲

FIG. 9.共a兲 Convergence of the time derivative of the variance to periodic behavior for periodic reversal.共b兲 The corresponding evo-lution of the variance.

We define a dimensionless observation time as t_obs= t /␶₁, and the closely related modal observation time of mode n as tn,obs= t /␶n= n2t /␶1. In this section, we restrict ourselves to uniform initial particle distributions.

B. Effective relaxation time

Let both modes in Eq.共49兲 describe a significant relax-ation with respect to the observrelax-ation time. Take, for example, modes n1= 1 with␶1= 1500 d and n2= 5 with␶1/␶2= 25 and t = 50. Figure 12 shows that the variance by the fourth-order model matches the exact variance, using the relaxation time as fit parameter 共␶fit= 130兲. In contrast, the standard defini-tion 共28兲 overestimates the relaxation time 共␶std= 780兲 and erroneously predicts the variance. To provide an improved a priori estimate of the relaxation time, consider the evolution equation for the approximate variance共B5兲 for uniform con-ditions共u0= v0兲. Replacing the approximate variance 共␴app,x

2 _兲

by the exact variance for unidirectional flow共␴_c,x2 兲 共15兲 and rearranging terms gives the evolution of the effective relax-ation as a function of time,

␶ef f共t兲 = ⳵␴app,x 2 _共t兲 ⳵t − 2Dmol 2␶␴_v2−⳵ 2_␴ app,x 2 _共t兲 ⳵t2 =

兺

n=1 ⬁ vn 2_␶ n共1 − e−t/␶n兲

兺

n=1 ⬁ vn 2_{共1 − e}−t/␶_n兲 . 共50兲 This relaxation time matches␶std共28兲 in the Fickian limit and equals the harmonic, modal velocity weighted average of the modal relaxation times for short times,

lim t_↓0 ␶ef f共t兲 =

兺

n=1 ⬁ v_n2

兺

n=1 ⬁ vn 2 /␶n . 共51兲

This explains why␶std 共28兲 overestimates the relaxation time for short times. Let us now define␶avgas the average of

␶ef f共t兲 over the time domain 关tmin, tmax兴 of interest,

␶avg=

1

tmax− tmin

冕

t_min t_max

␶ef f共t兲dt. 共52兲 Using this relaxation time␶avgin the approximate model

yields a very good fit with the exact variance Fig. 12.

C. Separation of scales

Next we study a two-scale problem共e.g., n1= 1, n2= 100兲 and consider observation times in between the two modal relaxation times. Figure 13 shows that the standard relax-ation time, Eq. 共28兲, overestimates the variance. Equation 共28兲 describes an averaged relaxation of both scales, while only the larger scale still undergoes a relaxation for the time scale considered here. With scale separation, we put the small relaxed scale共n2兲 together with the molecular diffusion in a microscale dispersion coefficient关30兴,

Dmicro= Dmol+ 1 2vn2

2_␶

2, 共53兲

and only describe the relaxation of the larger scale. This reduces the velocity variance to the contribution by the scale undergoing relaxation, ␴v,rel 2 =1 2vn1 2_␶ 1, 共54兲

and changes the variance of the approximate model to

␴ap,x,rel 2 _{= 2D} microt +vn₁ 2_␶ 1共␶1关e −t/␶_{1− 1兴 + t兲. 共55兲} With this scale separation, we obtain an excellent fit with the exact variance共Fig. 13兲.

If we put a scale共say n兲 in a small-scale dispersion coef-ficient and neglect its relaxation, the共relative兲 error we make in the contribution of this scale to the exact variance共15兲 is bounded, Error =

冏

vn 2_␶ n共t + 关e−t/␶n− 1兴兲 − vn 2_␶ nt vn 2_␶ nt

冏

艋

冏

vn 2_␶ n 2 vn 2_␶ nt

冏

=␶n t = 1 t_n,obs. 共56兲

FIG. 10. The normalized modal cyclic dispersion coefficient

关the term between square brackets in Eq. 共47兲兴 as a function of the

dimensionless reversal time共tr/␶n兲.

FIG. 11. Effect of the dimensionless reversal time共t_r/␶₁兲 on the development of the variance共normalized at the variance for unidi-rectional flow兲.

This relative error decreases for increasing observation time. We now postulate a separation criterion based on this single scale analysis. Scale 共n兲 is added to the microscale dispersion coefficient in the approximation and its relaxation ignored if the relative error bound共56兲 drops below a certain error value␧, ␶n t = 1 n2 ␶1 t = 1 n2 1 tobs 艋 ␧⇔n艌

冑

␧/tobs. 共57兲 We define the observation time-dependent separation scale nsep共tobs;␧兲 as the real-valued 共n兲 for which the equality in Eq.共57兲 holds,

n_sep共t_obs;␧兲 =

冑

␧/t_obs. 共58兲

D. General approach

In an empirical approach, we combine the effective relax-ation time formulrelax-ation with scale separrelax-ation. Filter␻*共tobs; n兲 is the fraction of mode n that at observation time tobsis put into the relaxation part, and关1−␻*_共t

obs; n兲兴 is the fraction put into the microscale dispersion coefficient共Dmicro兲. We define

␻*_共t

obs; n兲 as

␻*_共t

obs,n兲 =

冦

1, n⬍ N+关nsep共tobs;␧兲兴 nsep−关N+共nsep兲 − 1兴, n = N+关nsep共tobs;␧兲兴 0, n⬎ N+关nsep共tobs;␧兲兴

冧

.

共59兲

FIG. 12. Comparison of the exact variance with the variance by the fourth-other upscaled model for a bimodal velocity field共49兲 with n1= 1 and n2= 5. In the fourth-order model, the relaxation time is used as fit parameter共␶_fit兲 or is defined as␶_std共28兲 and␶_avg共52兲, respectively.

FIG. 13. Comparison of the exact variance with the variance by our upscaled model for a bimodal velocity field共49兲, with n1= 1 and n₂= 100. For␶_std共28兲, the relaxation of both scales is described in an averaged sense. Scale separation puts the smaller scale共n2兲 in a microscale dispersion coefficient共53兲.

FIG. 14. Error⌬sq共66兲 for the velocity fields of Table I as a function of the observation time. 共a,d兲 No scale separation; optimization of

␶ only. 共b,e兲 Scale separation withnsep共58兲 and␶avg共64兲. 共c,f兲 Scale separation with nsep共58兲 and␶std共28兲. The errors 共a,b,c兲 are for unidirec-tional flow. The errors共d,e,f兲 for flow including reversal of direction. The␣ and ⑀values were found by trial and error. The results in 共b兲 and

Here the observation-time-dependent scale separator nsep共tobs;␧兲 is defined according to Eq. 共58兲 and the function N+ is the first integer larger than nsep 关N+共x兲=minn苸N兵n ⬎x其兴. Filter␻*_共t

obs; n兲 共59兲 is defined such that the contribu-tion of each mode to the relaxacontribu-tion part smoothly vanishes in time, in descending mode order. At a specific point in time, only a single scale关mode n=N+_共n

sep兲兴 is partitioned over the relaxation part and Dmicro. In the limit of the observation time to infinity, the separator scale nsep 共58兲 becomes zero and according to Eq.共59兲 all scales are described by an effective dispersion coefficient. However, for the accuracy of the ap-proximation it is always favorable to at least describe the relaxation of a part of the largest scales. We relate this part to the 共normalized兲 contribution of all modes to the effective dispersion coefficient in the Fickian limit. Let us define mode nmin共␣兲 as n_min共␣兲 = min n苸N

再

n

冏

兺

k=1 n vk 2_␶ k艌␣Def f,⬁

冎

. 共60兲 Given ␣, all modes up to mode nmin will be fully de-scribed by the relaxation part independent of the magnitude of the observation time. With this modification, the fraction of scale n关␻共tobs; n兲兴 that is put in the relaxation part finally reads ␻共tobs;n兲 =

再

1, n艋 nmin共␣兲 ␻*_共t obs;n兲, n ⬎ nmin共␣兲

冎

. 共61兲

In the upscaled approximation including scale separation, the microscale dispersion coefficient consists of the molecu-lar diffusion and the fractions of the scales that are relaxed,

Dmicro= Dmol+

兺

n=1 ⬁ 关1 −␻共tobs;n兲兴vn 2_␶ n. 共62兲

The relaxation part of the velocity variance 共␴_v,rel2 兲 de-scribes the fraction of those velocity modes which undergo relaxation, ␴v,rel 2 =1 2n=1

兺

⬁ ␻共tobs;n兲vn 2 . 共63兲

We compute the relaxation time of this velocity fraction analogously to Eq.共52兲, ␶avg= 1 ⌬t

冕

t_min t_max

兺

n=1 ⬁ ␻共tobs;n兲vn 2_␶ n共1 − e−t/␶n兲

兺

n=1 ⬁ ␻共tobs;n兲vn 2_{共1 − e}−t/␶_n兲 dt. 共64兲 With these definitions, the microscale dispersion coeffi-cient D_micro共62兲 replaces the molecular diffusion coefficient Dmolin Eq.共2兲. In Eq. 共22兲, the relaxation part of the vari-ance ␴_v,rel2 共63兲 is substituted for the total velocity variance

␴v

2

and␶avg共64兲 replaces the relaxation time ␶. The

corre-sponding time derivative of the variance reads

⳵t␴ap,x,rel 2 _{共t兲 = 2 ⫻}

再

Dmicro+␴v,rel 2 _␶ avg共1 − e−t/␶avg兲, t艋 tr Dmicro+␴v,rel 2 _␶ avg共1 − 关2e+tr/␶avg− 1兴e−t/␶avg兲, t ⬎ tr.

冎

共65兲

E. Optimization and evaluation

The upscaled model with scale separation may be inter-preted as a function of two independent variables, the sepa-ration scale 共N兲 and the relaxation time 共␶兲. We want to quantify the accuracy of the model with respect to the be-havior of the variance in time. The error measure⌬sq共nsep,␶兲 computes the relative error in the time derivative of the vari-ance that is made by the approximation共65兲 on a time do-main关tmin, tmax兴 as

⌬sq共N,␶兲 =

冕

t_min t_max 关⳵t␴c,x 2 _{共t兲 −⳵} t␴ap,x 2 _{共t;N,␶兲兴}2_dt

冕

t_min t_max 关⳵t␴c,x 2 _{共t兲 − 2D} mol兴2dt . 共66兲 Here ⳵t␴ap,x

2 _{is the time derivative of the variance by the} approximation共65兲 and⳵t␴c,x

2

the exact variance time deriva-tive共15兲.

F. Test fields

We test the upscaled model for a parabolic velocity field, two cosine velocity fields, and two layered velocity fields. The latter two consist of equally sized layers of different velocities, as shown in Table I. We consider flow with and without flow reversal. For each field, we vary the observation time t_obsin the range t_obs苸关10−3_{, 10}1_{兴. For each t}

obs, we com-pute共n_sep,␶_avg兲 关Eqs. 共58兲 and 共64兲兴 and 共n_sep,␶_std兲 关Eqs. 共58兲 and共28兲兴. We evaluate the approximations as function of tobs by comparing the error 共⌬sq兲 evaluated on the time range 关tmin, tmax兴=关0,2兴␶1tobs.

G. Results

Figures 14共b兲 and 14共e兲 show that the scale separator nsep 共58兲 combined with relaxation time ␶avg 共64兲 gives small

errors. Moreover, the results are almost indistinguishable from those obtained by independently optimizing n and ␶ with respect to ⌬_sq共N,␶兲. On the contrary, the same scale separator combined with the standard relaxation time ␶std 共28兲 gives errors that are much larger 关Figs. 14共c兲 and 14共f兲兴.

As mentioned, Eq.共28兲 overestimates the relaxation time by overemphasizing the relaxation times of large scales for short times. In case we do not apply scale separation, the errors are large 关Figs. 14共a兲 and 14共d兲兴. For unidirectional flow, the relaxation time has a bigger impact than scale separation. For flow including reversal, the largest errors are made if we do not separate scales关Fig. 14共d兲兴.

XI. DISCUSSION

Hassanizadeh关19兴 and Tompson and Gray 关20兴 each de-rived a macroscopic multidimensional non-Fickian disper-sion model. Hassanizadeh 关19兴 derived his equation by ex-ploring the mass and momentum balances for a solute and solvent at the macroscale. Tompson and Gray 关20兴 applied the method of volume averaging on the convection diffusion equation. Although the model of Camacho 关Eqs. 共8兲 and 共39兲兴 is derived from upscaling to the megascopic scale to a 1D representation, it has the same functional form. The equivalence of these models suggests that dispersion at the macroscale might be interpreted as multidimensional Taylor dispersion. In Verlaan 关21兴 and Berentsen 关18兴 it is shown that the fourth-order model is able to explain the scale de-pendency of the dispersivity observed in echo experiments measured by Rigord et al.关10兴. Moreover, they showed that in a qualitative sense the fourth-order model is able to repro-duce the dependency of the effective dispersion coefficient on the Péclet number共see Bear 关22兴兲.

XII. CONCLUSIONS

The fourth-order approximation as proposed by Camacho provides acceptable results only in the case of single scale problems and uniform initial conditions. We characterized the scale interaction by analyzing the moments derived through spectral analysis. The analysis allowed us to formu-late an upscaled model with effective parameters that can be calculated a priori from the velocity distribution. Separation of scales in combination with a new definition of the effec-tive relaxation time allowed us to extend the model to obtain good results also for multiscale problems. In addition, we studied the effects of periodic flow reversal. Periodic flow reversal results in relaxation of the moments to periodic be-havior in the same relaxation time as for the unidirectional flow, even if the reversal time is smaller than the relaxation time. The effective dispersion coefficient is a monotonously increasing function of the dimensionless cycle time共tr/␶兲. It asymptotically converges for increasing tr/␶towards the ef-fective dispersion coefficient in the absence of any flow re-versal.

APPENDIX A: EXACT SPATIAL MOMENTS OF THE 2D CONVECTION DIFFUSION EQUATION

1. Mean position of higher concentration modes Unidirectional flow

The evolution expression of the non-normalized first mo-ment of mode n, Eq.共9兲, for k=1, reads

冋

␶n −1 + ⳵ ⳵t

册

Ecn,x,1= 1 2关2v0+v2n兴Mcn, +vnM0+ 1 2_m=1

兺

_⫽n ⬁ 关vm+n+v兩m−n兩兴Mc_m. 共A1兲 We replace M_c

n共t兲 and Mcm共t兲 with Eq. 共11兲 on the RHS of

Eq. 共A1兲 and solve the result for the initial condition, Ec_n,x,1共0兲=x0Mc_n

init

, since all tracer is released at x = x0. For E_c n,x,1共t兲, one obtains Ec_n,x,1共t兲 = Mc_n init x0e−t/␶n+vnM0␶n关1 − e−t/␶n兴 +

兺

m=1⫽n ⬁ ␶n␶m ␶m−␶n vm+n+v兩m−n兩 2 Mcm init_关e−t/␶_{m− e}−t/␶_n兴 +

冋

v0+ v2n 2

册

Mcn init_te−t/␶_{n. 共A2兲} Flow reversal

The evolution expression of E_c

n,x,1

rev _{共t兲 for flow after}

rever-sal is equal to Eq.共A1兲 with the signs of the modal velocities changed,

冋

␶n −1 + ⳵ ⳵t

册

Ecn,x,1 rev _{= −}1 2关2v0+v2n兴Mcn共t兲 − vnM0 −1 2_m=1

兺

_⫽n ⬁ 关vm+n+v兩m−n兩兴Mc_m共t兲. 共A3兲 The solution to Eq.共A3兲 is found by replacing Mc_n共t兲 and M_c

m with Eq.共11兲 and using continuity of Ecn,x,1at the

mo-ment of reversal, E_c

n,x,1

rev _共t

r兲=Ec_n,x,1共tr兲, and yields E_c n,x,1 rev _{共t兲 = −}

冋

v0+ v2n 2

册

Mcn init_{共t − 2t} r兲e−t/␶n+ Mc_n init x0e−t/␶n −vnM0␶n关1 − 2e−共t−tr兲/␶n+ e−t/␶n兴 −

兺

m=1⫽n ⬁ ␶n␶m ␶m−␶n vm+n+v兩m−n兩 2 Mcm init_{⫻ 共e}−t/␶_m −关2e关共␶m−␶n兲/␶m␶n兴tr_{− 1}兴e−t/␶n兲. 共A4兲

TABLE I. Description of the velocity fields v共y兲 under consideration 共d is the height of the field兲 共in 10−1_{m / d兲.}

Parabolic cos 1 + cos 10 ⌺ cos共1¯10兲 5 layer 10 layer

15关1 −共2y

_Ⲑ

d − 1兲2兴 cos共␲y

Ⲑ

d 兲 +cos共10_␲y

_Ⲑ

d兲 兺 n=1 10 cos共 n␲

Ⲑ

d y兲 layer velocities v =关1,3,9,1,8兴 layer velocities v =关3,2,12,7,13,1,8,5,4,13兴

2. Variance of average concentration Unidirectional flow

The evolution equation for the second moment belonging to c0共x,t兲 reads 关Eq. 共8兲 for k=2兴

⳵Mc₀,x,2 ⳵t = 2Dmol+ 2v0␮c0,x+

兺

n=1 ⬁ vn M0 E_c n,x,1. 共A5兲 Using, ␴_c,x2 = Mc₀,x,2−␮c₀,x 2

, the evolution equation for the centered second共variance兲 yields

⳵␴c,x 2 ⳵t = 2Dmol+

冋

2v0− 2 ⳵␮c₀,x ⳵t

册

␮c0,x共t兲 +

兺

n=1 ⬁ vn M0 Ec_n,x,1共t兲. 共A6兲 We substitute the known expressions for ␮c₀,x 共13兲 and Ec_n,x,1共A2兲 on the right-hand side and find expression 共15兲 as the solution to the initial condition␴_c,x2 共0兲=0.

Flow reversal

The evolution expression for the variance for flow after reversal is equivalent to Eq.共A6兲 with the signs of the ve-locity changed, ⳵␴c,x,rev 2 ⳵t = 2Dmol+

冋

− 2v0− 2 ⳵␮c₀,x rev ⳵t

册

␮c0,x rev_共t兲 −

兺

n=1 ⬁ vn M0 E_c n,x,1 rev _{共t兲. 共A7兲}

The solution共39兲 to this ODE is obtained by substituting Eq.共36兲 for␮_c,xrev and Eq.共A4兲 for E_c

n,x,1

rev _{and using the}

con-tinuity of the variance at the moment of reversal,␴_c,x,re2 _v共tr兲 =␴_c,x2 共tr兲. For a more detailed derivation, see 关18兴.

APPENDIX B: SPATIAL MOMENTS OF THE FOURTH-ORDER UPSCALED EQUATION

Multiplication of the fourth-order equation共23兲 with xk_, integrating the result over the x domain, gives the following general evolution equation for the kth moment of the concen-tration: ␶⳵2Mc,x,k ap ⳵t2 + ⳵M_c,x,kap ⳵t = kv0Mc,x,k−1 ap + k共2v0+␤兲 ⳵M_c,x,k−1ap ⳵t +共k 2 − k兲关Dmol+␴v 2_␶−共v 0+␤兲v0␶兴Mc,x,k−1 ap + 2共k2_{− k兲D} mol␶ ⳵M_c,x,k−2ap ⳵t −

冉

兿

i=0 2 共k − i兲

冊

Dmol共␤+ 2v0兲␶Mc,x,k−3 ap −

冉

兿

i=0 3 共k − i兲

冊

Dmol 2 _␶ Mc,x,k−4 ap . 共B1兲

The moments of the Telegraph equation are obtained by taking Dmolequal to zero in Eq.共B1兲.

1. Mean particle position

The evolution expression for the mean position of the average concentration yields

␶⳵ 2_␮ c,x ap ⳵t2 + ⳵␮c,x ap ⳵t = −v0. 共B2兲

Initially all particles are released at x0. If the particles are initially distributed nonuniformly over the height, the initial particle velocity of the particles is u0, which may be different from v0, ␮c,x ap_{共0兲 = x} 0, ⳵␮c,x ap ⳵t 共0兲 = u0. 共B3兲 The solution of Eq.共B2兲 subject to result is given by Eq. 共24兲.

2. Spatial variance in the particle distribution The evolution expression for the second normalized mo-ment reads ␶⳵2Mc,x,2 ap ⳵t2 + ⳵Mc,x,2 ap ⳵t = 2共2v0+␤兲␶ ⳵␮c,x ap ⳵t + 2v0␮c,x ap_{+ 2关D} mol+␴v 2 ␶−共v0+␤兲v0兴. 共B4兲 Using ␴_c,xap,2= M_c,x,2ap −共␮_c,xap兲2 and the expressions for ␮_c,xap given in Eq. 共24兲, the following evolution equation for the variance is obtained: ␶⳵ 2_␴ c,x ap,2 ⳵t2 + ⳵␴c,x ap,2 ⳵t = 2共Dmol+␴v 2_{␶兲 + 2␶␤共u} 0−v0兲e−t/␶ − 2␶共u0−v0兲2e−2t/␶. 共B5兲 Initially, the variance is zero while the initial increase of the variance is by molecular diffusion only,

␴c,x

ap,2_{共0兲 = 0,} ⳵␴c,x ap,2

⳵t 共0兲 = 2Dmol. 共B6兲 The solution of Eq.共B5兲 with respect to Eq. 共B6兲 is given by Eq.共26兲.

APPENDIX C: PERIODIC FLOW REVERSAL

The differential equation共in a moving frame of reference兲 for the first moment of a higher-order concentration mode in the共Nr兲th reversal cycle reads

冋

⳵ ⳵t+

1

␶n

册

Ec_n,x,1共t兲 = 共− 1兲NrvnM0. 共C1兲 The first moment共Ec_n,x,1兲 is initially zero and continuous each time the flow is reversed. Consequently, the solution to Eq.共C1兲 is given by

Ec_n,x,1共t兲 = − 2M0vn␶n

兺

k=0 N_r