DOI 10.1007/s11242-006-9040-0 O R I G I NA L PA P E R
Upscaling, relaxation and reversibility of dispersive flow
in stratified porous media
C. W. J. Berentsen · C. P. J. W. van Kruijsdijk · M. L. Verlaan
Received: 19 August 2005 / Accepted: 16 July 2006 / Published online: 13 January 2007 © Springer Science+Business Media B.V. 2007
Abstract Dispersive tracer released in a unidirectional velocity field belonging to a stratified porous of finite height describes a transition, called relaxation, from a convective dominated behaviour for short times to Fickian behaviour for asymptotic long times. The temporal relaxation state of the tracer is controlled by the transverse mixing term. In most practical applications, the orders of the time and length scales of the relaxation mechanism are such that in an upscaled model of a stratified medium the dispersive flux is in a pre-asymptotic state. Explicit modelling of the relaxation of the dispersive flux in the pre-asymptotic region is required to improve the accuracy. This paper derives a pre-asymptotic one-dimensional upscaled model for the trans-verse averaged tracer concentration. The model generalises Taylor dispersion (Proc. R. Soc. London 219, 186–203 (1953)) and extends the method of Camacho (Phys. Rev. E 47(2), 1049–1053 (1993a); Phys. Rev. E 48 (1993b)) to dispersion tensors that may vary as function of the transverse direction. In the averaging step, the governing two-dimensional equation is first spectrally decomposed in terms of the eigenfunctions of the transverse mixing term. Next, the resulting modal relaxation equations are com-bined into an effective relaxation equation for the extended dispersive Taylor flux. Contrary to the one-dimensional Fickian approach, the upscaled model approximates the multi-scale relaxation behaviour as a single scale relaxation process and accounts
C. W. J. Berentsen· C. P. J. W. van Kruijsdijk · M. L. Verlaan
Department of Geotechnology, Delft University of Technology, Mijnbouwstraat 120, 2628 RX Delft, The Netherlands
C. W. J. Berentsen (
B
)Department of Earth Sciences, Hydrogeology group, Utrecht University, Budapestlaan 4, 3508 TA Utrecht, The Netherlands
e-mail: berentsen@geo.uu.nl C. P. J. W. van Kruijsdijk
Shell Canada Ltd., Calgary Research Centre, 3655 36 Street NW, Calgary, Alberta, Canada M. L. Verlaan
for the partial reversibility of convective dispersion upon reversal of the flow direc-tion. The upscaled model is evaluated against the original two-dimensional model by means of moment analysis. The longitudinal tracer variance predicted by our model is quantitatively correct in the short and long time limits and is qualitatively correct for intermediate times.
Keywords Relaxation· Reversibility · Taylor dispersion · Non-Fickian dispersion · Upscaling· Tracer flow · Stratified flow
List of symbols
B/ ˜B/BN [m/s] Modal velocity interaction matrix w.r.t.
{eigen-function base of −∂y(DT(y)∂y)/cosine Fourier
base/eigenfunction base of ˜DT,N}.
c [kg/m3] 2D concentration.
c0 [kg/m3] Transverse or height averaged concentration. cn/˜cn/cN,n [kg/m3] n-th modal concentration w.r.t. {eigenfunction
base of−∂y(DT(y)∂y)/cosine Fourier
base/eigen-function base of ˜DT,N}.
c/˜c [kg/m3] Infinite column vector with modal concentrations {cn/˜cn} for n= 1, . . ., ∞.
˜cN/cN [kg/m3] Finite column vector with modal concentrations
{˜cN,n/cN,n} for n = 1, . . . , N.
d [m] Total layer height or transverse width.
Deff [m2/s] Effective dispersion coefficient. DL [m2/s] Longitudinal dispersion coefficient.
DL,0 [m2/s] Height or transverse averaged of longitudinal
dispersion coefficient.
DL,n/ ˜DL,n [m2/s] n-th spectral mode of longitudinal dispersion
coefficient w.r.t {eigenfunction base of −∂y(DT(y)∂y)/cosine Fourier base}.
DL/ ˜DL/DL,N [m2/s] Modal interaction matrix resulting from spectral
decomposition of DL(y)c(x, y, t) w.r.t.
{eigenfunc-tion base of−∂y(DT(y)∂y)/ cosine Fourier
base/ei-genfunction base of ˜DT,N}.
dL/˜dL [m2/s] Infinite column vector with spectral modes of
lon-gitudinal dispersion coefficients{DL,n/ ˜DL,n} for
n= 1, . . . , ∞.
Dmol [m2/s] Molecular diffusion coefficient. DT(y) [m2/s] Transverse dispersion coefficient.
˜DT [m2/s] Matrix resulting from cosine Fourier
decomposi-tion of−∂y(DT(y)∂y).
˜DT,N [m2/s] ˜DTtruncated at N modes(= ˜DT(1 . . . N, 1 . . . N)).
Eξ,x,k k-th non-centred non-normalised spatial moment of quantityξ.
F [-] Tortuosity factor.
I [-] Unity matrix.
Mξ,x,k k-th non-centred normalised spatial moment of quantityξ.
M0 [kg] Tracer mass.
t/tr [s] {Time/Reversal time}.
Tcos,φ [-] Transformation matrix from cosine Fourier base to eigenfunction base of−∂y(DT(y)∂y).
TN [-] Transformation matrix from cosine Fourier base
truncated at N-modes to eigenfunction base of ˜DT,N.
v(y) [m/s] Longitudinal velocity.
v0 [m/s] Height or transverse averaged velocity.
vn/˜vn/vN,n [m/s] n-th modal velocity belonging to {eigenfunction
base of−∂y(DT(y)∂y)/cosine Fourier
base/eigen-function base of ˜DT,N}.
v/˜v [m/s] Finite column vector with modal velocities {vn/˜vn} for n = 1, . . ., ∞.
˜vN/vN [m/s] Infinite column vector with model velocities
{˜vn/vN,n} for n = 1, . . . , N.
x/y [m] {Longitudinal/transverse} co-ordinate. x0 [m] Initial longitudinal position.
Greek
αL/αT [m] {Longitudinal/Transverse} dispersivity.
φn(y) [-] n-th eigenfunction of operator−∂y(DT(y) ∂y).
λn/λN,n [s−1] n-th eigenvalue of {operator −∂y(DT(y)∂y)/
matrix ˜DT,N}.
/N [s−1] (Diagonal) matrix with eigenvalues of
{−∂y(DT(y)∂y) / ˜DT,N}.
µ Mean.
σ Co-variance or variance.
σ2
c,x [m2] Spatial variance belonging to concentration.
σuni [m2] Asymptotic deviation of variance from Fickian behaviour for particle distributions that are ini-tially distributed uniform over the height. τn/τN,n [s]= (λ−1n /λ−1N,n) n-th modal relaxation time.
τ [s] Relaxation time.
n, m With respect to the interaction of the n-th and m-th base functions (n, m> 0).
rev After flow reversal.
T Transverse.
TEL Generalised Telegraph equation.
v Velocity.
x Spatial.
∞ For asymptotic long times.
overbars
∼ [-] W.r.t cosine Fourier base.
1 Introduction
Convective dispersion is a physical mechanism that results from the interaction of a (transverse) mixing process with a spatially non-uniform velocity-field. The classi-cal convection–dispersion equation (CCDE) that is traditionally used to describe this mechanism from a one-dimensional point of view does not always lead to an adequate representation. In particular, the following situations are not well described by the 1D CDDE:
• In natural rocks often several scales of heterogeneities (e.g. layers) are encoun-tered which exhibit a different dispersion behaviour than those that can be ob-tained from the CCDE (Lake and Hirasaki 1981; Gelhar 1993; Gelhar et al. 1979). • The interaction between convection and molecular diffusion can be significant
(Taylor 1953, 1954; Aris 1956; Saffman 1959, 1960; De Josselin de Jong 1958). • When the flow direction is reversed, convective spreading will also reverse. The
CCDE describes the dispersion as an irreversible mixing mechanism (Hiby 1963; Rigord et al. 1990).
an upscaled model for dispersive tracer flow in unidirectional velocity fields in which the mixing may be an arbitrary function of the transverse co-ordinate. Contrary to Verlaan et al. (2000), Verlaan (2001) no assumptions are made with respect to the multi-scale relaxation behaviour prior to the upscaling step. Moreover, the param-eters within the height averaged upscaled model can again be directly expressed as function of the small scale properties. For a commonly used form of the mixing term, the model behaviour is analysed in more detail.
2 Model description and physical behaviour
2.1 Dispersive flow in unidirectional velocity field
Similar to Berentsen et al. (2005), we consider a macroscopically stationary unidirec-tional velocity field v(y), which corresponds to flow through an infinite two-dimensional porous media bounded in the transverse or y-direction, (Fig. 1). At t= 0 we inject a finite mass (M0) of ideal1 tracer uniformly or non-uniformly distributed over the height (y) at longitudinal position x0. In addition to convective transport the tracer particles are subject to a micro scale mixing mechanism that is assumed to exhibit un-correlated Fickian behaviour at the macro scale (i.e. scale of observation). Contrary to Berentsen et al. (2005) this mixing mechanism may now vary as function of the transverse (or y) co-ordinate. We further assume that macroscopically the convective flow is strictly parallel to the layer orientation.
2.2 Physical behaviour—relaxation from convective to Fickian behaviour and reversibility of dispersion upon reversal of the flow direction
In the absence of a mixing mechanism, each tracer particle moves with a fixed veloc-ity along its streamline. Consequently, the megascopic longitudinal particle variance increases quadratically with time. Microscopic mixing has two effects. The longitudinal
d=5m
x y
x0
Magnitude unidirectional velocity field v(y) Microscopic mixing
Initial tracer distribution at x0
0.9m/day
1.0m/day
0.8m/day 0.1m/day
0.3m/day
Fig. 1 Consider a layered porous medium of which the longitudinal co-ordinate (x) stretches from
-∞ to +∞ and the transverse co-ordinate (or height) (y) from [0, d]. A longitudinal velocity field v(y) is present in this medium. Tracer particles are initially released in this velocity field at longitudinal position x= x0uniformly or non-uniformly distributed over the transverse direction. The tracer particles are subject to microscopic mixing
component of the mixing process smoothes the sharp front propagated by convection. In addition, and more importantly, transverse dispersion moves particles away from their initial streamlines thereby changing the particle velocity. Over time, each tracer particle will sample each vertical position equally frequently, and the relation between a particle and its initial streamline is entirely lost. Consequently, transverse disper-sion smoothes concentration differences in the transverse direction, and as a result reduces the longitudinal spreading caused by the velocity profile. In the long time limit, a dynamic equilibrium between the convective and diffusive spreading mecha-nisms that exhibits Fickian behaviour is established, referred to as Taylor dispersion (see Taylor 1953, 1954; Aris 1956). The relaxation timeτ, which is the characteristic time for this transition process2, is directly related to the transverse mixing time. It is proportional to the (finite) field height squared and inversely proportional to the transverse mixing.
The (pre-) asymptotic behaviour of the spatial moments has been extensively stud-ied as function of the statistics of the underlying permeability field for flow in layered porous media of which the height is finite or infinite and for flow parallel and not parallel to the layer orientation (Marle et al. 1967; Gelhar et al. 1979; Matheron and de Marsily 1980).
Now consider the case where after a certain time we reverse the flow direction. Assuming convective transport only, the process is fully reversible and the particles are transported back to their original positions. If we reverse the flow and include micro scale mixing, the relaxation process starts all over again. However, as initially the particles return longitudinally along the same path as they arrived from, the transport process demonstrates a (partially) reversible behaviour. The variance of the height-averaged particle positions decreases as the particles turn back along their original paths. Due to transverse mixing, the velocity of a particle becomes less and less correlated to the velocity history of its forward movement. As a result convec-tive dispersion takes over and the variance once again increases monotonously. The process is identical to the original relaxation process. For a detailed description of the relaxation concept caused by isotropic molecular diffusion and the partial reversibility of dispersion for flow reversal, see Berentsen et al. (2005).
3 Mathematical model—the 2D uCCDE
For the situation and assumptions described in sub-section 2.1, the tracer concen-tration is governed by the macroscopic two-dimensional unidirectional conventional convection dispersion equation (uCCDE):
∂c ∂t + v(y) ∂c ∂x= ∂ ∂x DL(y)∂c ∂x + ∂ ∂y DT(y)∂c ∂y (1)
where the longitudinal DL(y) and transverse dispersion DT(y) are arbitrary piece
commonly defined at the macro scale as (see Bear 1972; Perkins and Johnston 1963): DL(y) = Dmol/F + αL|v (y)|
DT(y) = Dmol/F + αT|v (y)|
(2)
in which Dmolis the molecular diffusion, F is the tortuosity factor andαLandαTare
the longitudinal dispersivity and transverse dispersivity, respectively. In the examples that follow we use this specific definition for the dispersion coefficients. However, in the more general derivations these coefficients may be an arbitrary function of y. We assume zero mass-flux over the top and bottom boundary and an exponential convergence of c in the longitudinal (x) direction towards infinity:
xk∂xmc(±∞, y, t) = 0, ∀k, m ∈N (3)
which implies that the concentration behaviour sufficiently far away from the front remains unaffected. Throughout this paper, we restrict the examples to a tracer that is initially injected uniformly over the height:
c(x, y, t = 0) = M0δ (x − x0) . (4)
The objective is to derive an expression for the evolution of the megascopic transverse averaged tracer concentration c0(x, t) =d1
d
0 c(x, y, t) dy.
4 Spectral decompositions and truncation
The relaxation state of the 2D uCCDE in time is controlled by the transverse mixing term. The relaxation times or, equivalently, the characteristic transverse mixing times are the inverse values of the eigenvalues of the differential operator Ly = −∂y(DT∂y).
Section4.1gives the spectral expansion of the 2D uCCDE in terms of these eigen-values and the corresponding eigenfunctions. These eigenfunctions are the actual spatial scales within the 2D uCCDE. Since the temporal behaviour of each eigenfunc-tion is controlled by its own relaxaeigenfunc-tion time only, such a decomposieigenfunc-tion disentangles the multi-scale relaxation behaviour of the 2D uCCDE in time. Moreover, it serves as the theoretical starting point for the moment analysis and the upscaling part. Unfor-tunately, we can generally not explicitly determine these eigenfunctions. To overcome this, we first make a spectral cosine expansion of the 2D uCCDE (Sect.4.2). Section
4.3shows that “diagonalising” the infinite cosine expansion equals the desired eigen-function expansion. Section4.4shows how to obtain the diagonalised representation for a finite Fourier truncation.
4.1 Theoretical decomposition in terms of eigenvalues and eigenfunction
The relaxation times of the 2D uCCDE are inverses of the eigenvalues of the self-adjoint differential operator Ly = −∂y
DT∂y
. To obtain these eigenvalues, we look for finite square integrable functions f :[0, d] → ∞ that are solutions to the so-called Sturm-Liouville eigenvalue problem (EVP):
with DT(y) > 0 and subject to boundary conditions
∂f ∂y(0) =
∂f
∂y(d) = 0. (6)
In mathematical literature, this Sturm-Liouville problem is thoroughly analysed (Zauderer 1989). EVP (5), (6) has an infinite number of non-negative real-val-ued countable eigenvaluesλi, i = {1, 2, . . .} with λi ≤ λi+1 and the eigenfunctions φi(y), i = {0, 1, 2, . . .} form a orthogonal base, where the orthogonality is expressed
with respect to the following inproduct: φi,φj = 2 d d 0 φiφj dy= 2 i= j = 0 δij otherwise. (7)
Completeness guarantees that we may decompose each square (y-)integrable con-centration function c(x, y, t) that satisfies (6) as a convergent sum of eigenfunctions φi(y): c(x, y, t) = ∞ k=0 ck(x, t) φk(y) (8)
where the concentration modes ck(x, t) are defined as:
ck(x, t) = c (x, y, t) , φφ k(y) k,φk = 1 2c (x, y, t) , φk(y) , k = 0 c (x, y, t) , φk(y) , k > 0. (9) In particular,λ = 0 is the (smallest) eigenvalue of (5) and (6) with eigenfunction φ0(y) = 1 and the height averaged concentration c0(x, t) is its coefficient. After expanding c(x, y, t), DL(y) and v(y) in the uCCDE in terms of eigenfunctions3, we
multiply the result with eigenfunctionφnand average it over y∈ [0, d]. In doing so,
we transform the 2D uCCDE in an infinite set of 1D evolution equations for the spectral coefficients cn(x, t). For (n = 0), we obtain the evolution for the averaged
concentration: ∂c0 ∂t + v0∂c0 ∂x − DL,0∂ 2c 0 ∂x2 = − ∂JE,T ∂x (10)
where the extended Taylor flux (JE,T) and modal dispersive flux (JE,n) are defined as:
JE,T≡ ∞
n=1
JE,n with JE,n=
1 2 vncn(x, t) − DL,n∂cn(x, t) ∂x . (11)
For (n> 0), we find the evolution equation of concentration mode cn(x, t) belonging
to eigenfunctionφn: ∂cn ∂t + cn τn + v0∂cn ∂x − DL,0∂ 2c n ∂x2 +∞ m=1 Bn,m∂cm ∂x − DLn,m∂ 2c m ∂x2 = − vn∂c0 ∂x − DL,n∂ 2c 0 ∂x2 . (12)
Here,τn = λ−1n is the modal relaxation time. The Fourier coefficients ξk(withξ =
{v, DL} are similar to (9) given by:
ξk= ξ(y), φk(y)
φk(y), φk(y)
(13) and express the modal velocities (vk) and modal longitudinal dispersion coefficients
(DL,k) in terms of the eigenfunction base. Finally the coefficientsξn,mwithξ = {B, DL}
defined as,
ξn,m= (ξ(y) − ξ
0), φn(y)φm(y)
φn(y), φn(y) =
∞l=1ξlφl(y), φn(y)φm(y)
φn(y), φn(y)
(14) account for the modal interactions of the velocity (B) and longitudinal dispersion (DL) between two non-constant eigenfunctions. In (12), the following terms can be
distinguished. The first term between square brackets is the so-called relaxation term and accounts for the evolution of the concentration mode cn(x, t) in time. The second
term accounts for average convection and longitudinal dispersion, the third term for the modal interactions between different concentration modes. The right hand side of (12) is the driving force and expresses the interaction with the averaged concentration. This expansion decouples the 2D uCCDE in terms of relaxation times, while ei-genfunctions remain coupled via the spatial operators present in the 3rd term of (12). The advantage of this decomposition is that it disentangles the temporal relaxation behaviour. Moreover, it serves as starting point for an upscaled model and allows us to determine the spatial moment belonging to the 2D uCCDE sequentially. Besides, Sect.5 will show that the spatial modal coupling does not affect the longitudinal mean and variance within the tracer distribution for tracer uniformly injected in the transverse direction.
4.2 Decomposition in terms of cosine Fourier series
Unfortunately, it is in general not possible to derive explicit expressions for the ei-genfunctions. To overcome this problem, we decompose the 2D uCCDE first in terms of cosine series after which the system may be diagonalised (see next sub-section). Similar to Camacho (1993b) we replace the concentration, velocity and longitudinal dispersion (1) by cosine Fourier series:
ξ (x, y, t) = ξ0(x, t) + ∞ n=1 ˜ξn(x, t) cosn(y), (15) with cosn= 1, n= 0 cos(nπy/d) , n > 0 (16)
forξ = {v, DL, c}. Next, we multiply the result by cosnand integrate it over y. For
The evolution of the higher order concentration modes (n> 0) reads in matrix vector notation: ∂ ˜c ∂t+ ˜DT˜c + ˜v0I∂ ˜c ∂x− ˜DL,0I∂ 2˜c ∂x2 + ˜B ∂˜c ∂x− ˜DL∂ 2˜c ∂x2 = −˜v∂c0 ∂x + ˜dL∂ 2c 0 ∂x2 . (18)
Here, the vectors˜c =˜c1,˜c2,. . . ,˜v =˜v1,˜v2,. . . T and ˜dL= ˜DL,1, ˜DL,2,. . . T con-tain the cosine Fourier modes of c(x, y, t), v(y) and DL(y), respectively. Moreover, for
m, n> 0 the matrix entries ξn,m = (ξ(y) − ξ0), cosmcosn for ξ = {∂yDT(y)∂y, v(y),
DL(y)} can be determined explicitly and take the form:
m= n m= n ˜Bn,m ˜v2n/2 ˜vm+n+ ˜v|m−n| /2 ˜DLn,m ˜DL,2n/2 ˜DL,m+n+ ˜DL,|m−n| /2 ˜DTn,m (nπ/d) 2˜D T,0− ˜DT,2n/2 (π/d)2− (mn) ˜D T,m+n+ (mn) ˜DT,m−n /2. (19)
The structure of Eqs. (17) and (18) is identical to Eqs. (10) and (12) following from the eigenfunction decomposition. The main difference is that generally matrix ˜DTis
full, such that the relaxation behaviour of the different cosine Fourier modes are fully coupled in time. Only for a constant transverse dispersion coefficient, is matrix ˜DT
already diagonal and are the cosine and the eigenfunction base equivalent. 4.3 Diagonalisation of the infinite dimensional matrix ˜DT
The coupling of the temporal relaxation mechanisms within the cosine Fourier expan-sion is expressed by the fullness of matrix ˜DT. In theory it is possible to “diagonalise” the infinite dimensional matrix ˜DT, such that we are able to retrieve the eigenfunction
expansion. Completeness guarantees that any square integrable function f(y) satisfy-ing boundary conditions (6) can be fully expressed in terms of the eigenfunction base with coefficient vector cT= (c1, c2,. . .) or in terms of the {cosn} base with coefficient
vector˜cT =˜c1,˜c2,. . .
. Hence, each cosncan be expressed in terms of a convergent
eigenfunction series as well, and vice versa. As a consequence we can construct base transformation matrices Tφ,cosand Tcos,φwith entries defined as:
(Tφ,cos)k,n= cosn,φk = φk, cosn = (Tcos,φ)n,k (20)
with which we can transform the cosine-base vector˜c into the corresponding eigen-function-vector c and vice versa:
c= Tφ,cos˜c and ˜c = Tcos,φc (21) Applying (21) twice we also find:
c= Tφ,cos˜c= Tφ,cosTcos,φc=Tφ,cosTcos,φc= Ic. (22) Hence from (20) and (22) follows that Tφ,cos= T−1cos,φ = TTcos,φ. Equation (21) explicitly shows that each cosndescribes a linear combination of eigenfunctions each of which
has its own relaxation time. Now consider square integrable functions c(y) satisfying (6) for whichc, Lyc
is finite as well. If we expand the termc, Lyc
With the use of the base transformation matrices (21) we can rewrite this expression as an eigenfunction expansion: ˜cT ˜D T˜c = Tcos,φc T ˜D T Tcos,φc = cTT φ,cos˜DTTcos,φ c (24)
of which the structure of (Tφ,cos˜DTTcos,φ) is not yet known. Alternatively by using the
eigenvalue definition (5) we can expandc, Lyc
directly in terms of eigenfunctions: c, Lyc = 2 d d 0 cLyc dy= ∞ k,l=1 ckcl 2 d d 0 φkλlφldy= ∞ l=1 c2lλl= cTc < ∞. (25)
where is diagonal with the eigenvalues as entries (i.e. ii = λi= τn−1). Since (23),
(24) and (25) holds for all functions c(y) for whichc, Lyc
< ∞ we obtain the equality:
= Tφ,cos˜DTTcos,φ. (26)
This means that the infinite dimensional matrix ˜DT may be “diagonalised” with the
use of base transformation matrix Tφ,cos. Hence, by multiplying the cosine represen-tation (17), (18) with the infinite dimensional base transformation matrix Tφ,cos(Eq. 20), we obtain the eigenfunction representation (10), (12). Unfortunately the infinite dimensional base-transformation matrices are in practice not available.
4.4 Diagonalisation of a truncated cosine series
In practice, we have to work with a finite cosine Fourier series truncation. Diagonali-sation of such finite truncation results in an eigenspace representation of which the relaxation mechanism in time is decoupled between the specific eigenfunctions. How-ever, the eigenfunctions and eigenvalues belonging to the truncated case will generally differ from the infinite eigenfunction expansion. As the infinite dimensional matrix ˜DTis positive definite (see Appendix 9.1) and symmetric, each truncated sub-matrix
˜DT,N= ˜DT(1 : N, 1 : N) has these properties. Consequently, ˜DT,Nonly has real and
strictly positive eigenvalues. Since any finite symmetric matrix is diagonalisable, we can write ˜DT,Nas:
N = TTN˜DTTN with TTN= T−1N (27)
where the diagonal matrixN only has non-zero values on the main diagonal being
the eigenvalues of the truncated matrix ˜DT,N, sorted by increasing magnitude (ii<
jj, i< j). TN is a truncation mode (N) dependent transformation matrix of which
the columns contain the corresponding (normalised) eigenvectors. Exploring this, we change from the modal concentration vector˜cNto vector cNusing:
cN = TTN˜cN, ˜cN = TNcN. (28)
where the coefficients (without tilde∼) are defined as:
BN= TTN˜BNTN, DL,N = TTN˜DL,NTN, dL,N= TTN˜dL,N, vN = TTN˜vN (31)
in which the coefficients ˜ξN (with ξ = {DL, B, v}) are defined as in Sect.4.2, but
truncated at N modes. The relaxation matrixNis diagonal with entries:
(N)n,n= λN,n= τN,n−1. (32)
For N= ∞, the previous sub-section showed the equivalence between the diagona-lised set and the eigenfunction base. We were not able to give a formal proof that the eigen-functions and values belonging to the truncated cosine space converge in the limit for N→ ∞ to the infinite dimensional eigen-functions and values (for N = ∞). Therefore the “convergence of the series” is investigated in Sect.6.1by analysing the behaviour of the relaxation times and modal velocity for large N. Comparison with the infinite case is done by numerical comparison of the behaviour of concentration profiles and spatial moments between Fourier truncations and the full 2D solution in Sect. 6.2.
5 Spatial moments for unidirectional flow
The infinite eigenfunction representation is the starting point for the derivation of the exact spatial moments belonging to the 2D uCCDE. Moment analysis allows us to investigate the behaviour of any truncated Fourier set with the spatial moments belonging to the full 2D behaviour. Appendix 9.2 or Berentsen (2003) gives the der-ivation of the spatial moments in detail. In the remainder of this paper we limit the discussion, for simplicity, to tracer distributions that are initially uniform in the trans-verse direction4. For such uniform initial conditions, the transverse tracer distribution over the streamlines is not affected by redistribution of the tracer by transverse dis-persion. Hence, the mean tracer position (µc,x) moves with the height averaged fluid
velocity (µc,x = x0+ v0t). The (megascopic) spatial variance (σc,x2 ) belonging to the
height averaged tracer concentration reads (see Appendix 9.2 or Berentsen 2003): σ2 c,x= 1 M0 ∞ −∞c0(x, t) x− µc,x 2 dx= 2DL,0t+ ∞ n=1 v2nτn t+ τn e−τnt − 1. (33) The 1st contribution is the transverse averaged macroscopic longitudinal dispersion. The 2nd contribution is due to the interaction of the longitudinal velocity field and transverse dispersion. For short times (t<< τ1), this contribution is proportional to the time squared and the variance within the velocity field:
σ2 c,x,conv(t) = 2DL,0t+ 1 2 ∞ n=1 v2nt2≡ 2DL,0t+ σv2t2 (34)
where the variance of the velocity field is given (in equivalent expressions) by : σ2 v = 1 2 ∞ n=1 v2n= 1 2v Tv=1 2˜v T˜v (35)
where v=v1, v2, v3,. . . . T
the modal velocity vector. Over time (t>> τ1) the entire spatial spreading behaves Fickian and increases proportionally with time:
σ2
c,x= 2Deff,∞t+ σuni (36)
where the dispersion coefficient for asymptotic long time reads: Deff,∞= DL,0+ 1 2 ∞ n=1 v2nτn = DL,0+ 1 2v T−1v= D L,0+ 1 2˜v T ˜D−1 T ˜v (37)
and the (negative) constant contribution resulting from the relaxation process is given by: σuni= − ∞ n=1 v2nτn2. (38)
Figure 2 gives a typical example of the transition of the variance in time from an initial∼ t2 behaviour towards a∼ t behaviour for long times. The spatial moments belonging to the truncated and diagonalised Fourier set (29) and (30) are obtained by replacing infinity (∞) in the upper bounds by N, and the coefficients ξnbyξN,n(for
ξ = {τ, v, c, DL}). It is interesting to note that, when the tracer is distributed uniformly
over the transverse direction, the y-dependency of the small scale dispersion tensor results in exactly the same functional form of the spatial mean and variance as for the case of a constant dispersion coefficient (see Berentsen et al. 2005).
6 Robustness of the Fourier series truncation
6.1 Behaviour of the diagonalised set as function of the truncation mode
The eigen-values and eigen-functions belonging to a truncated cosine series depend on the truncation mode (N). This section investigates the sensitivity of the relaxation Fig. 2 Relaxation of the
variance and from
correlated-convective (∼ t2) for short times to Fickian behaviour (∼ t) for large times
times, the modal velocities and the effective dispersion coefficient to the truncation mode (N). We restrict the analysis to macroscopic dispersion coefficients of the form (2). Figure 3 is a typical example for the development of the relaxation timesτN,nas
function of the truncation mode N of the cosine Fourier series. It shows that the modal relaxation timeτN,nis nearly inversely proportional to the modal number (n) squared.
Moreover, Fig. 3 shows that the Fourier series truncation does not significantly affect the large relaxation times (provided that the truncation mode N is sufficiently large). This also holds for the corresponding modal velocities vN,n(Fig. 4). Finally, we
inves-tigate the behaviour of the dispersion coefficient (Deff,∞,N) of the truncated Fourier set for asymptotic long times as function of the mode N at which we truncate the Fourier series. The finite equivalent of (37) for finite N reads:
Deff,∞,N= DL,0+ N
n=1
v2nτn= DL,0+ vTN−1N vN = DL,0+ ˜vTN˜D−1T,N˜vN. (39)
Figure 5 shows a typical example of the convergence of the effective dispersion coeffi-cient towards some asymptotic value as function of the mode at which the Fourier series is truncated.
6.2 Comparison between the truncated Fourier set and the full 2D behaviour It remains to be investigated whether a finite Fourier truncation converges to the infinite Fourier case, or equivalently, to the full 2D behaviour. We compare the diag-onalised cosine Fourier set truncated at N= {10, 50, 100, 400} modes with the full 2D behaviour. The set of one-dimensional evolution Eqs. (29) and (30) is solved using a 3rd order accurate 1D finite difference solver. The full 2D model is evaluated using a 2nd order accurate finite difference scheme. We illustrate the comparison by means of the layered field shown in Fig. 1. The velocity field consists of five layers of equal width (0.1 m) with layer velocities v= [0.9, 0.1, 1.0, 0.3, 0.8] m/day. We inject the tracer uniformly over the transverse direction. We specify DLand DT as in (2) and set the
diffusion to Dmol= 10−6m2day−1and the dispersivities toαL= 1 mm, αT= 0.1 mm.
Figure 6 clearly shows that the concentration profile of the truncated Fourier set slowly
Fig. 3 Comparison of the
modal relaxation times (τN,n (32)) between the Fourier series truncated at five different modes
N= (100, 125, 150, 175, 400).
The underlying velocity field is shown in Fig. 1. The dispersion is defined according to (2) with
Fig. 4 Comparison of the
modal velocities (vN,n(31)) for five truncation modes
N= (100, 125, 150, 175, 400) at
which the cosine Fourier series is truncated. The underlying velocity field is shown in Fig. 1. The dispersion is defined according to (2) with αT= 10−4,αL= 10−3and Dmol= 10−6 0 10 20 30 40 50 60 10–3 10–2 10–1 mode n
modal velocity v [m/day]
N,n 400 100 125 150 175
Fig. 5 Total effective
dispersion coefficient (39) as function of truncation mode of the original Fourier series
0 25 50 75 100 125 150 175 0 1 2 3 4 5 modes [-] D [m /day] eff 2
converges towards the full 2D behaviour. After 400 modes, the agreement becomes acceptable.
6.3 Cosine Fourier series truncation versus eigenfunction truncation
Assuming that the infinite Fourier series is approximated accurately by a cosine Fourier series with 400 modes, we compare two additional truncation methods with the full 2D. Method (A) truncates the cosine-Fourier series while method (B) trun-cates the “eigen-function” series:
(A) Truncation of the Cosine Fourier series
Fig. 6 The concentration
profile at time t= 500 days of the full 2D model (2D) compared to the concentration profiles for a cosine Fourier series truncated at N= 10, 50, 100 and 400 modes 100 200 300 400 500 0 0.5 1 1.5 2 2.5 3 3.5 distance [m] concentration [kg/m ] 3 x 10-2 2D 10 50 100 400 preserved: ˆvA N,n= ⎧ ⎪ ⎨ ⎪ ⎩ vN,n, n≤ N vT400v400− vTNvN = 400 n=N+1˜v 2 n, n= N + 1. (40)
The relaxation time of the last mode is defined such that the (N+ 1) modes describe the same asymptotic dispersion coefficient for asymptotic long times as 400 Fourier modes do: ˆτA N,n= ⎧ ⎪ ⎨ ⎪ ⎩ τ400,n, n≤ N Deff,∞,400− Deff,∞,N 1 2(ˆv A N,N+1)2 = v T 400−1400v400− vTN−1N vN (ˆvA N,N+1)2 , n= N + 1. (41)
(B) Truncation of the “eigen-function” series
We first diagonalise the entire cosine Fourier set5. This approximately gives us the eigen-function representation of the EVP, (5) and (6). In method (B), the relaxation of each of these N “eigen-modes” is described explicitly. In addition, we use one single mode that accounts for the relaxation of the neglected eigen-modes (i.e. mode N+ 1, . . . , 400). The velocity assigned to this last mode is again defined in such a way that the total velocity variance described by 400 Fourier modes is preserved:
ˆvB N,n= ⎧ ⎪ ⎨ ⎪ ⎩ v400,n, n≤ N vT400v400− vT400,Nv400,N = 400 n=N+1 v2400,n, n= N + 1 (42) where v400,N = v400,1,. . . , v400,N
times: ˆτB N,n= ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ τ400,n, n≤ N Deff,∞,400− Deff,∞,400,N 1 2(ˆv C N,N+1) 2 = 400 n=N+1 v2400,nτn2 400 n=N+1 v2400,n , n= N + 1. (43)
In both methods, all interactions of the first N-modes with the additional mode via the matrices B, DT DL are set to zero. The comparison between both methods
and the full 2D behaviour is illustrated using the velocity field described in the pre-vious sub-section. Figure 7b shows that the “eigen-function” truncation method B already gives an excellent fit for N= 5, while the cosine-Fourier truncation method A (although improving for increasing N) is still significantly different from the 2D profile for N= 50 (Fig. 7a).
As mentioned in Sect. 3, each cosndescribes a fraction of each eigenfunctionφn
with relaxation timeτn. Hence, a low order cosine Fourier truncation will neglect a
significant part of eigenfunctions with large relaxation times. As a consequence, the additional mode (N+ 1) in method A still needs to describe a significant part of the full 2D relaxation behaviour and will be characterised by a “large” modal relaxation time. In contrast, method B truncates the diagonalised approximately-infinite Fourier set sorted on decreasing relaxation times. This causes the method to improve rapidly with increasing N, and the relaxation time to decrease rapidly. Figure 8a indeed shows that the modal relaxation time of mode N+ 1 rapidly decreases for method B while it remains “large” for method A. Similarly, Fig. 8b shows that fraction of the asymptotic dispersion coefficient that need to be described is negligible at N= 5 with method B, while with method A it still exceeds the 10%. Note that the scale of y-axis in Fig. 8a is logarithmic. 100 200 300 400 500 distance [m] 100 200 300 400 500 distance [m] x 10-2 x 10-2 2D 5 10 20 50 b a 2D 1 3 5 10 concentration [kg/m ] 3 concentration [kg/m ] 3 0 0.5 1 1.5 2 2.5 3 3.5 0 0.5 1 1.5 2 2.5 3 4 3.5
Fig. 7 The concentration profile at time t= 500 days of the full 2D model (2D) compared to the
0 10 20 30 40 100 101 102 103 truncation mode N 0 10 20 30 40 truncation mode N
relaxation time of mode N+1 [-]
Original Diagonalised 0 0.8 1 0.4 0.2 0.6 Fraction of D by mode N+1 eff, φ Original Diagonalised a b
Fig. 8 Comparison of the original and diagonalised Fourier set by means of the (a) the relaxation
time of the additional mode and (b) the fraction of the asymptotic dispersion coefficient described by the additional mode. Note the logarithmic axis of the relaxation time in plot (a)
Fig. 9 Typical shape of the
effective relaxation timeτeff(t) normalised byτeff,∞as function of the observation time (t/τ1) for the velocity field shown in Fig. 1 and with the transverse dispersion defined as in (2). Parameters that are used in this figure are
αT= 0.05 m, Dmol= 10−6m2/day (τ1= 214 days, τeff,init= 0.12 days, τeff,∞= 10.1 days) 10–6 10–4 10–2 100 102 10–1 100 observation time t/τ [–]
normalised eff. relaxation time [–]
eff eff,init eff, 1 τ τ τ ∞
7 The upscaled c–J models
The non-truncated eigenfunction representation is an infinite set of one-dimensional equations that describes the exact 2D behaviour. For computational convenience a single evolution equation for the dispersive Taylor flux (11) is desirable. To achieve this we follow an approach similar to Camacho’s (1993a, b) upscaling method. We write the total extended Taylor flux (11) as:
JE,T= ∞ n=1 LJ,ncn, with LJ,n= 1 2 vn− DL,n ∂ ∂x . (44)
all modes (see Appendix 9.3): ∂JE,T ∂t + JE,T τeff + (v0+ Beff) ∂JE,T ∂x − DL,0+ DL,eff ∂ 2J E,T ∂x2 = −σ2 v ∂c0 ∂x + σv,DL ∂2c 0 ∂x2 − σ 2 DL ∂3c 0 ∂x3 . (45)
Here,τeffis the effective relaxation time, Beff the coefficient that accounts for the modal interactions of the velocity and DL,eff the coefficient that accounts for the
modal interaction by longitudinal dispersion. Theoretically, these space and time dependent coefficients are defined analogously as in Camacho (1993b):
JE,T τeff ≈ ∞ n=1 JE,n τn , Beff∂JE,T ∂x ≈ ∞ n,m=1 ˆBn,mLJ,n∂cm ∂x , DL,eff∂ 2J E,T ∂x2 ≈ ∞ n,m=1 ˆDLn,mL˜J,n∂ 2c m ∂x2 . (46)
For practical purposes, these coefficients are assumed to be constant. The remaining deterministic coefficients in (45) are the velocity variance (σv2), the covariance of the velocity and the longitudinal dispersion (σv,DL) and the variance of the longitudinal dispersion (σD2
L), which are given by: σ2 v = 1 2 ∞ n=1 v2n= 1 2 ∞ n=1 ˜v2 n, σv,DL= ∞ n=1 vndL,n= ∞ n=1 ˜vn˜dL,n, σ2 DL= 1 2 ∞ n=1 d2L,n= 1 2 ∞ n=1 ˜d2 L,n. (47)
In case the dispersion is defined according to (2),σv,DLandσ
2
DLcan both be expressed in terms of the velocity variance:
σv,DL = 2αLσ
2
v, σD2L = α
2
Lσv2. (48)
The contribution of the extended Taylor flux to the evolution equation of c0is already given by (10). In Appendix 9.3, we eliminate the Taylor flux from the c–J model, i.e. expressions (10) and (45), and we obtain a 4th order expression that is the extension of the upscaled model of Camacho to dispersive mixing. Expressed in a moving frame of reference this equation reads:
∂2c 0 ∂t2 + 1 τeff ∂c0 ∂t + Beff∂ 2c 0 ∂x∂t = σ2 v + DL,0 ˜τeff ∂2c 0 ∂x2 + DL,0+ DL,eff ∂ 3c 0 ∂x2∂t +BeffDL,0− σv,DL ∂3c 0 ∂x3 +−DL,0DL,eff+ σD2L ∂4c 0 ∂x4. (49)
by DL,0<< σv2τeff). In this case, the contribution by small scale longitudinal disper-sion can be neglected (DL,0 = DL,eff = σv,DL = σ
2
DL = 0) in (see Berentsen et al. 2005), and we obtain a generalised Telegraph type of equation:
τeff∂ 2c 0 ∂t2 + ∂c0 ∂t + Beffτeff∂ 2c 0 ∂x∂t = σv2τeff∂ 2c 0 ∂x2 . (50)
7.1 Initial and boundary conditions for the c–J model
In practice, we solve for the first order Eqs. (10) and (45) rather than for the 4th order equation (49). The initial condition for JE,Tfollows immediately from the definition
of the Taylor flux (11), the Fourier base-transformation (27) and the (initial) value of the cosine Fourier coefficients{˜cn(x, 0), ˜vn, ˜DL,n}. For the initial condition (4) we get:
c0(x, 0) = M0δ (x − x0) , JE,T(x, 0) = 0. (51)
The boundary conditions for c0(x, t) and JE,T, being the equivalents of (3), read:
xk∂xmJE,T(±∞, y, t) = xk∂xmc0(±∞, y, t) = 0, ∀k, m ∈N. (52) For the 2nd order model (50), the initial condition is a result of its hyperbolic reversible nature and reads:
c0(x, 0) = M0δ (x − x0) , ∂tc(x, 0) = 0. (53)
7.2 Spatial variance of the upscaled models
For unidirectional flow and an initial distribution that is uniform in the transverse direction, the variance computed by the c–J model (Eqs. (10) and (45)) or equiva-lently the upscaled 4th order model (49) is of similar form as the variance of the full 2D problem (33) but restricted to a single relaxation scale (see Appendix 9.4):
σ2 c,x,Ap(t) = 2DL,0t+ 2σv2τeff t+ τeff e−t/τeff− 1. (54)
In the short time or convective limit (t << τeff), the relaxation part (2nd term) becomes proportional to t2and independent of the relaxation time:
lim t/τeff↓0σ
2
c,x,Ap(t) = 2DL,0t+ σv2t2. (55)
In the long time or Fickian limit(t/τeff >> 1), the whole variance becomes propor-tional to t: lim t/τeff→∞σ 2 c,x,Ap(t) = 2 DL,0+ σv2τeff t. (56)
The variance computed by the 2nd order Telegraph model is equal to (54) for DL,0= 0.
7.3 Comparison of the upscaled models with the full 2D model
et al. (2005), we define the effective relaxation time in the upscaled model as the averaged over the time domain [0,t] (see Appendix 9.5):
τeff(t) = 1 t t 0 ∞ n=1 v2nτn 1− exp−˜t/τn ∞ n=1 v2 n 1− exp−˜t/τn d˜t. (57)
For short (t<< τ1), respectively, long time (t>> τ1) this definition reduces to:
lim t/τ1↓0τeff≡ τeff,init= ∞ n=1 v2n ∞ n=1 v2 n/τn , lim t/τ1→∞τeff≡ τeff,∞= ∞ n=1 v2nτn ∞ n=1 v2 n . (58)
This time invariant relaxation timeτeff,∞equals the definition of Camacho (1993b), hence the approximation describes the proper variance for asymptotic long times. Figure 9 shows the typical form of the effective relaxation time as function of the observation time (t/τ1). With a small value for the effective relaxation time for small observation times, the model attempts to account for the relaxation of the smaller scales by means of a single scale relaxation process. Note however that the contribu-tion of the smaller scales is restricted outside the short time convective limit as their relaxation time approximately scales with 1/n2as shown in Fig. 3.
In the upscaled model, the modal interaction coefficient (Beff) is defined such that the 3rd moment from the approximation matches the exact 3rd moment in the short time limit (t<< τ1) (see Berentsen 2003):
Beff= Mv3/σv2 (59) where M3v= 1d0d(v(y) − v0)3dy= 14 ∞ m,n=1 vnvm[vm+n+v|m−n|,m=n] is the 3rd centred moment of the velocity. We a priori assume that the modal interactions of the higher components within DL(y) have a negligible effect on the development c0(x, t) and set (DL,eff) to zero. Remark that neither DL,effnor Beffaffect the spatial variance (54) of the upscaled concentration.
Figure 10a–c shows the variance of the models in the time domain t∈0, 1000days for flow including reversal of the flow direction at time tr = 500 days (see Appendix
0 500 1000 0 0.5 1 1.5 2 2.5 3 x 104 va ri ance [m 2] a 0 500 1000 0 0.5 1 1.5 2 2.5 x 104 time [days] time [days] time [days] b c 0 500 0 500 1000 1500 2000 2500 exact 4th order ( eff) 2nd order ( eff) 4th order (eff,) τ τ τ ∞
Fig. 10 Comparison of the variance of the 2nd-order (50) and 4th-order (49) upscaled models for τeff(t), the variance of 4th order upscaled model for τeff,∞with the exact variance for flow including reversal of direction (33) for Dmol= 10−6m2s−1and (a)αL= αT= 10−6m (τeff= 1.1 × 104days,
τeff,∞ = 8.8 × 104days), (b)αL = αT = 5 × 10−4m (τeff = 620 days, τeff,∞ = 1, 001 days) (c)
αL= αT= 0.05 m (τeff= 10 days, τeff,∞= 10.2 days).The underlying velocity field is shown in Fig. 1
100 200 300 400 500 0 0.2 0.4 0.6 0.8 1 distance [m] c o n c ent ra ti on [ k g /m 3] 100 200 300 400 500 0 0.2 0.4 0.6 0.8 1 distance [m] 200 250 300 350 400 0 0.2 0.4 0.6 0.8 1 distance [m] full 2D 4thorder, τeff 2ndorder,τeff a b c
Fig. 11 Comparison of the cumulative concentrations∫x+∞˜c(˜x, t)d˜x of the 2nd-order (50) and
4th-order (49) models and the full 2D model (1) at the time of reversal t= tr= 500. Dmol= 10−6m2s−1 and (a)αL = αT = 10−6m (τeff = 1.1 × 104days), (b)αL = αT = 5 × 10−4m (τeff = 620 days), (c)αL= αT= 0.05 m (τeff= 10 days).The underlying velocity field is shown in Fig. 1. Note that the (Gibbs) overshoots in the models are caused by the higher order upwind scheme that is used for the spatial discretisation
For the same parameter values as in Figs. 10, 11 shows that the concentration profiles of both upscaled models are also nearly identical. For short times the hyper-bolic part ˜τeff∂t2c0 in the Telegraph model dominates over the parabolic part (∂tc0) and its concentration profile shows two shocks, typical for a 2nd order hyperbolic equation. Even though the profiles of the upscaled models seem to give a poor rep-resentation of the exact profile (Fig. 11a), they match the exact solution up to the 2nd moment (Fig. 10a). For intermediate time the resemblance between the upscaled profiles and the concentration profile by the 2D model improves (Fig. 11b). In the long term, both upscaled models demonstrate Fickian behaviour, characterised by the S-shaped concentration profile and give a good match with the exact solution (Fig. 11c).
7.4 Relaxation mechanisms of the upscaled model
For DT(y) defined as (2), the transverse dispersion matrix ˜DT(19) following from the
˜DT= (DmolN/F + αTv0M) /d2, (60) where the form of the symmetric dimensionless M depends on the Fourier expansion of the velocity, and the dimensionless N is diagonal with entries Nnn = n2π2.
Con-sider the two limiting cases. If diffusion dominates the transverse mixing, ˜DTreduces
to ˜DT = DmolN/Fd2. The effective relaxation time (57) in the upscaled model is inversely proportional to Dmol:
τeff(t) ∼ (d/π)2 Dmol/F
(61) and the relaxation is purely a time dependent process independent of the applied velocity. For asymptotic long times (t>> τ1), the effective relaxation time reads
τeff,∞= d πneff,1 2 1 Dmol/F (62) where neff,1 is the effective scale in the velocity field in the diffusive limit. In case transverse convective dispersion dominates the transverse mixing, ˜DT reduces to
˜DT= αTv0M/d2. In this limit the effective relaxation time (57) in the upscaled model is inversely proportional to v0.
τeff(t) ∼ (d/π) 2 αTv0
(63) and the relaxation is a space dependent process. At mean velocity v0(= L/t) an averaged distance L is covered in time t and the relaxation state (t/τeff) is space dependent: t τeff ∼ t π2α Tv0 d2 = t π2α T(L/t) d2 = αTL (π/d)2. (64)
For asymptotic times (t>> τ1), the effective relaxation time reads τeff,∞= d πneff,2 2 1 αTv0 (65) where neff,2is the effective scale in the velocity field in the dispersive limit. For inter-mediate situations the relaxation time is a combination of diffusion and convective dispersion directly related to ˜D−1T . In the asymptotic limit (t >> τ1), the effective relaxation time (τeff,∞) takes the form:
τeff,∞= k1 2 d π 2 1
n2eff,1Dmol/F + n2eff,2αTv0
. (66)
Fig. 12 k2as function of the Péclet number
(NPe= αTv0/Dmol) for a parabolic velocity field and the five layer field shown in Fig. 1
10-2 100 104 106 1 1.05 1.1 1.15 1.2 1.25 1.3 k2 N Pe
five layer field parabolic
scale. The behaviour of the five layer field is controlled by a broad range of Fourier scales.
7.5 Practical relevance of the non-Fickian transition behaviour
The time needed to reach Fickian behaviour may easily be of the same order of mag-nitude as the time scale of modelling of the practical problems in the oil/gas industry and hydrology. Consequently, when upscaled models are applied one should be care-ful to properly account for the pre-asymptotic state. In case diffusion dominates the spreading in the transverse direction, the largest possible effective relaxation time is limited by the relaxation time of the largest Fourier scale (Berentsen et al. 2005) (neff,1= 1):
τeff,∞= D1 mol
d2
π2. (67)
Fluids in the gas phase: Imagine a layered porous medium with a total height of (d=)1 m for which the velocity field is dominated by its first order Fourier modes. A typical value for the molecular diffusion of two gasses at T= 293 K and P = 100 atm is Dmol= 10−2m2Day−1. The corresponding relaxation time (67) is 10 Days.
Fluids in the liquid phase: A typical diffusion coefficient for two component in the liquid phase is of the order Dmol= 1e − 10 m2s−1or equivalently Dmol = 1e − 5 m2Day−1. For a velocity field of height (d=)1 m that is dominated by its first Fourier component, the corresponding relaxation time is even 10,000 Days.
Consider again a velocity field that is dominated by its first order Fourier component and let us assume that the macroscopic dispersion coefficient obeys (2). At sufficiently large velocities the relaxation mechanism is dominated by transverse dispersion and the corresponding effective relaxation time reads (for neff,2= 1):
The time for the flow needed to travel over a distance L with mean velocity v0is
t= L/v0. (69)
From (68) and (69), it follows that the spatial relaxation length for whicht/τeff,∞= 1 is given by: Lt/τeff,∞=1= d 2 π2 1 αT . (70)
For a total layer height of (d=)1 m and a transverse dispersivity (αT) of 0.1 mm, the
relaxation length is Lt/τeff,∞=1≈ 103m. Typical groundwater velocities are in the range [0.1–10] mDay−1. Hence the relaxation time for a layer of 1 m is forαT = 0.1 mm in
the range [102–104] Days.
8 Conclusions
• We derived an upscaled model that explicitly accounts for (a) the dynamics causing the relaxation of the tracer dispersion from short time correlated convec-tive behaviour to asymptotic long time Fickian behaviour and (b) the dynamics causing the partial reversibility (e.g. reduction) of the spreading that occurs upon reversal of the flow direction.
• The transverse dispersion term entirely controls the temporal relaxation state. The relaxation times are inverses of the eigenvalues of the transverse dispersion term. In theory the multi-scale relaxation behaviour present in the two-dimensional problem can be disentangled by a spectral expansion of the 2D uCCDE in terms of the corresponding eigenfunctions. In practice these eigenfunctions cannot be determined explicitly and are obtained in an approximate sense by diagonalisation of a finite truncated cosine-Fourier series expansion.
• For a uniform initial tracer distribution the upscaled model gives a quantitatively accurate prediction of the longitudinal variance in the short and long term limit. For intermediate times the agreement is qualitatively correct.
• The upscaled model is able to account for the reduction of the dispersion under reversal of flow. This is a major advantage with respect to the 1D classical convec-tion dispersion equaconvec-tion.
• For particle distributions that are initially uniform in the transverse direction, the spectral eigenfunction representation that follows from a transverse varying dis-persion behaves the same with respect to the mean and variance as in the case of a constant transverse dispersion.
• The effective relaxation time in the upscaled model is a functional combination of transverse diffusion (Dmol) and transverse dispersion (αT v0). For increasing “Péclet” number (αT v0/Dmol, withαTthe transverse dispersivity) the
relaxation time shifts from a time dependent to a scale dependent process.
9 Appendices
9.1 Proof DTis positive definite
Consider the functionγ (x, t): γ (x, t) = d 0 c(x, y, t) ∂ ∂y −DT(y) ∂ ∂yc(x, y, t) dy . (71)
where c(x, y, t) = c (x, y, t) − c0(x, t) is the fluctuation on the height averaged con-centration. Partial integration of (71) with respect to y and using that flow over the y-boundaries is not possible, it turns out thatγ (x, t) is positive for non-zero concen-tration fluctuations (c= 0): γ (x, t) = − c(x, y, t) DT(y) ∂ ∂yc(x, y, t) d 0 + d 0 c(x, y, t) ∂ ∂y −DT(y) ∂ ∂yc(x, y, t) dy = 0+ d 0 DT(y) ∂c(x, y, t) ∂y 2 dy > 0. (72)
Moreover, by replacing the concentration(s) and the dispersion coefficient DT(y) by
cosine Fourier series and integration over y we find (see Berentsen 2003) as well that: γ (x, t) = d 0 c(x, y, t) − c0(x, t) ∂ ∂y DT(y) ∂ ∂y c(x, y, t) − c0(x, t) dy = ⎡ ⎣ ∞ k,m,n=1 cn(x, t) DT,kcm(x, t) d 0 cosnπy d ∂ ∂y cos kπy d ∂ ∂ycos mπy d dy = ˜c (x, t)T ˜D T˜c (x, t) (73) with ˜c = ˜c1,˜c2,. . . t
the vector containing the higher order concentration modes. Combining (71) and (73) shows that ˜DThas to be positive definite. With ˜DTpositive
definite, all its eigenvalues are strictly positive as well, since for each eigenvector x belonging to eigenvalueλ we find that:
0< xT
˜DTx
= xT(λx) = λ |x|2. (74)
9.2 Spatial moments for flow including reversal
The moments are derived in detail in Berentsen (2003). Here we provide a brief overview. Consider the release of a tracer with finite mass M0 in the spatial domain x∈ (−∞, +∞). The kth spatial moment of the nth concentration mode is defined as:
Ecn,x,k= +∞
−∞ x kc
The 0th moments or “mass” in the nth concentration mode is denoted by Mcn. Mcn,x,k is the tracer mass normalised kth moment of the nth concentration mode, (Ecn,x,k/M0). Mcc
n,x,k, for k> 1, is the centred normalised kth moment (Mcn,(x−Mc,x,1),k). For sim-plicity the mean particle position (Mc0,x,1) is written asµc,xand the particle variance
(Mcc
0,x,2) asσ 2 c,x.
9.2.1 The evolution of the spatial moments of the spectral equivalent
Assuming that each concentration mode and its spatial derivatives converge expo-nentially to zero for x→ ±∞, the moment equations belonging to the 2D uCCDE are obtained by multiplying Eqs. (10) and (12) of the spectral equivalent with xkand
averaging the result over x. The equation for the kth normalised moment belonging to c0(x, t) reads: ∂Mc0,x,k ∂t = v0kMc0,x,k−1+ k (k − 1) DL,0Mc0,x,k−2 +1 2 ∞ n=1 kvn Ezn,x,k−1 Mc0 + k(k − 1) dL,n Ezn,x,k−2 Mc0 . (76)
The evolution of the kth non-normalised moment of a higher concentration mode (n> 0) reads: ∂Ezn,x,k ∂t + Ezn,x,k ˜τn = kv0Ezn,x,k−1+ ∞ m=1 kβm,nEzm,x,k−1+ vnkMc0M c c0,x,k−1 + ∞ m=1 αLm,n k2− k Ezm,x,k−2+ dL,nk(k − 1) Mc0M c c0,x,k−2 . (77)
The terms between curly brackets are additional to the moment expression for isotro-pic diffusion (DL,DT= constant). The 0th and 1st moments are not affected by these
terms (since Ecn,x,k−2= 0, k ≤ 1). See Berentsen (2003, 2005) for more detail. 9.2.2 Mean particle position
For a tracer initially released at x= x0, the expression for the mean position of the averaged concentration (µc,x ≡ Mc0,x,1), is identical to the case of isotropic diffu-sion (DL, DT= constant), see Berentsen (2003). For unidirectional flow and for flow
including reversal of direction (i.e. v(y)t>tr = −v (y)t≤tr) at time t = trit reads
In (78) and (79) the 1st bracketed term expresses the particle movement for particle that are distributed uniformly in the transverse direction. The 2nd term accounts for the relaxation of the temporal mean particle velocity towards the mean fluid velocity for particles that are initially non-uniformly distributed over the height.
9.2.3 The spatial variance for non-uniform particle distributions The general expression for the variance (σ2
c,x = Mc0c,x,2) for tracers that are initially released at position x = x0 reads (see Berentsen 2003, 2005 for a derivation) for unidirectional flow: σ2 c,x(t) = 2DL,0t+ ∞ n=1 v2nτn t+ τn e−τnt − 1 + ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ ∞ n=1 βn,n Minitzn M0 vn τ2 n− (t + τn) τne− t τn + ∞ n,m=1,m=n Mzinit m M0 vnβn,m τnτm τm− τn τm 1− e−t/τm+ τ n 1− e−t/τn ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ − ⎡ ⎣ ! 1 2 ∞ n=1 (τnvn) Mzinit n M0 e−τnt −1" 2⎤ ⎦+ ∞ n=1 DL,nτn Minitz n M0 1−e−τnt. (80) The 1st term expresses the relaxation of the variance for tracer distributions that are initially uniform over the height. The 2nd, 3rd and 4th terms account for the (differ-ence in) relaxation for tracer distributions which are initially non-uniform over the height. The 2nd term accounts for the modal interaction of the higher order Fourier modes. The 3rd term expresses the effect of the relaxation of the mean particle veloc-ity on the variance and equals the 2nd term in (78) squared. The 4th term of (80) is the only term that is additional with respect to the case of isotropic diffusion (constant DL, DT) and expresses the effect of the higher order Fourier modes of DL(y) on the
variance. For tracer distributions that are initially uniform over the height the variance takes the form (33).
Variance for flow including reversal of direction (limited to uniform tracer distribu-tions)
For a tracer distribution that is initially uniform in the transverse direction, the spatial variance by the full 2D model for flow including reversal of the flow direction at time t= trreads (see Berentsen 2003 or Berentsen et al. 2005):
σ2 c,x,rev(t) = 2DL,0t + ∞ n=1 v2nτn tr+ τn e−tr/τn−1+(t − t r) + τn e−(t−tr)/τn− 1 +τn etr/τn− 1e−t/τn+ τ n e−tr/τn− 1 , t> tr. (81)
the time correlation of the velocity of a particle before flow reversal with the velocity of that particle before flow reversal. This term dominates the spreading behaviour directly after flow reversal and causes the variance to show a temporal decrease as the velocities alter sign. With time this term relaxes to a negative constant and the 2nd term within the relaxation part takes over.
9.3 Derivation of the upscaled c–J model and 4th order model
Starting point of the derivation forms the eigenfunction representation being the evolution equation for the averaged concentration (10) together with the evolution equation of the modal concentration (12) where the extended Taylor flux (JE,T) and
modal dispersive flux (JE,n) are defined as in (11). Multiplication of Eq. (12) with the
differential operator Ln= 12[vn− DL,n∂x] gives the evolution of modal dispersive flux
JE,n: ∂JE,n ∂t + JE,n τn + v0∂JE,n ∂x − DL,0∂ 2J E,n ∂x2 + ∞ m=1 ˆvn,m ∂ ∂x− ˆDL,n,m ∂ 2 ∂x2 Lncm = − v2n 2 ∂ ∂xc0− vnDL,n ∂ 2 ∂x2c0+ D 2 L,n ∂3 ∂x3c0 . (82)
The evolution of the extended Taylor flux is now obtained by summing (82) over n= 1, . . . , ∞ ∂JE,T ∂t + JE,T τeff + (v0+ Beff) ∂JE,T ∂x − DL,0+ DL,eff ∂ 2J E,T ∂x2 = −σ2 v ∂c0 ∂x + σv,DL ∂2c 0 ∂x2 − σ 2 DL ∂3c 0 ∂x3 (83) where the approximations (46) have been made and where the deterministic coeffi-cients σ2
v,σv,DL,σ
2
DL are defined as in (47). The c–J model consists of the evolu-tion equaevolu-tion for the averaged concentraevolu-tion (10) and the evoluevolu-tion equaevolu-tion of the extended Taylor flux (83), respectively. Mulitiplication of (83) with the operator (−∂/∂x) gives: ∂ ∂t −∂JE,T ∂x + 1 τeff −∂JE,T ∂x + (v0+ Beff) ∂ ∂x −∂JE,T ∂x + DL,eff ∂ 2 ∂x2 −∂JE,T ∂x = σ2 v ∂2c 0 ∂x2 − σv,DL ∂3c 0 ∂x3 + σ 2 DL ∂4c 0 ∂x4 . (84)
Replacing (−∂xJE,T) in (84) with the left hand side of (10) gives the 4th order upscaled
Eq. (49) after combining terms of equal order.
9.4 Spatial moment of the 4th order upscaled and the c–J models
