Two-dimensional stochastic analysis of flow in leaky confined aquifers subject to spatial and periodic leakage

Two-dimensional stochastic analysis of flow in leaky confined aquifers subject

to spatial and periodic leakage

Hund-Der Yeh

*

_{, Ching-Min Chang}

Institute of Environmental Engineering, National Chiao Tung University, No. 75, Po-Ai Street, Hsinchu 300, Taiwan

a r t i c l e

i n f o

Article history:

Received 12 February 2009

Received in revised form 16 August 2009 Accepted 17 August 2009

Available online 22 August 2009

Keywords:

Nonstationary spectral approach Nonstationary head fields Leaky confined aquifers Heterogeneous media Coefficient of leakage

a b s t r a c t

This paper addresses the issue of flow in heterogeneous leaky confined aquifers subject to leakage. The leakage into the confined aquifer is driven by spatial and periodic fluctuations of water table in an over-lying phreatic aquifer. The introduction of leakage leads to non-uniformity in the mean head gradient and results in nonstationarity in hydraulic head and velocity fields. Therefore, a nonstationary spectral approach based on Fourier–Stieltjes representations for the perturbed quantities is adopted to account for the spatial variability of nonstationary head fields. Closed-form expressions for the variances of hydraulic head and specific discharge are developed in terms of statistical properties of hydraulic param-eters. The results indicate that the spatiotemporal variations in leakage leads to enhanced variability of the hydraulic head and of the specific discharge, which increase with distance from any arbitrary refer-ence point. The coefficient of leakage and the spatial structure of log transmissivity field and of the ampli-tude of water table fluctuation are critical in quantifying the variability of the hydraulic head and of the specific discharge.

Ó 2009 Elsevier Ltd. All rights reserved.

1. Introduction

Most of the existing stochastic analyses of groundwater flow in heterogeneous formations rely on the assumption of a constant mean head gradient in the treatment of flow problem. It is recog-nized that groundwater recharge may cause non-uniformity in the mean hydraulic head gradient, which results in nonstationarity in the statistics of random head fields in heterogeneous media

[7,17]. Several studies[3,9,12,18]have revealed that the assump-tion of a uniform mean head gradient leads to the soluassump-tion that fails to capture the covariances and cross-covariances in random nonstationary velocity fields. It is well known that the log conduc-tivity-head cross-covariance is of significance in controlling the variability of flow-dependent variables such as hydraulic head. This implies that the applicability of solutions developed based on the uniform mean head gradient (or uniform mean flow) assumption to the quantification of the flow perturbation caused by the presence of recharge is excluded. Motivated by this, this study is devoted to the development of closed-form analytic expressions for quantifying the variability of hydraulic head and specific discharge in a leaky confined aquifer, where the stochastic nature of hydraulic head is unstationary.

The changes in groundwater table in response to a given tempo-ral forcing are a common occurrence in many ground water basins. This temporal effect can directly affect the velocity variation and, therefore, the migration potential of contaminant plumes

[6,10,16]. Therefore, to make correct groundwater management decisions, it is important to understand the temporal effect of re-charge (leakage) on the mean behavior of flow system and the var-iability of the output processes (such as the varvar-iability of head and specific discharge). Several stochastic investigations of problems of flow and solute transport in a leaky confined aquifer have been presented over the past years [8,13,14,22,23]. However, to our knowledge, the application of the nonstationary spectral approach

[11,12] to quantify the variability of hydraulic head and specific discharge in a leaky confined aquifer subject to spatial and periodic leakage has not been done so far. We hope that the approach in this study provides a basic framework for quantifying and understand-ing field-scale flow processes in heterogeneous leaky confined aquifers and the findings will be useful in stimulating further re-search in this area.

2. Mathematical formulation of the problem

Groundwater flow through a leaky confined aquifer overlain by a leaky phreatic aquifer is considered to be essentially horizontal, so that it can be modeled using the vertically integrated form of the continuity equation combined with Darcy’s law[2]

0309-1708/$ - see front matter Ó 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.advwatres.2009.08.004

* Corresponding author. Fax: +886 35 726050.

E-mail addresses: [email protected] (H.-D. Yeh), [email protected]. edu.tw(C.-M. Chang).

Contents lists available atScienceDirect

Advances in Water Resources

j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / a d v w a t r e s

@ @Xi TðXÞ@/ðX; tÞ @Xi   /ðX; tÞ  /pðX; tÞ f2 ¼ S @/ðX; tÞ @t ; i ¼ 1; 2 ð1Þ where the leakage into the confined aquifer takes place through a leaky confining layer whose hydraulic conductivity is much smaller than that of the main aquifer. In Eq.(1), / is the hydraulic head in the confined aquifer, /pis the hydraulic head in the phreatic aqui-fer, T is the transmissivity, S is the storativity and f2is the coeffi-cient of leakage[2], which is defined as the root of the ratio of the thickness of the confining layer to its hydraulic conductivity. The conductivity of the confining layer does not vary significantly in space compared to the spatial variation of hydraulic conductivity in the confined aquifer. Therefore, the conductivity of the confining layer is assumed to be uniform. Expanding these terms in Eq.(1)

and dividing the equation by the nonzero transmissivity yields @2/ @X2 i þ@ln T @Xi @/ @Xi 1 T / /p f2 ¼ S T @/ @t ð2Þ

The hydraulic head in the phreatic aquifer and log transmissiv-ity ln T are regarded as realizations of stationary random fields per-turbed into means and corresponding zero-mean perturbations /pðX; tÞ ¼ h/pðX; tÞi þ hpðX; tÞ ¼ Hpþ hpðX; tÞ

ln TðXÞ ¼ hln TðXÞi þ f ðXÞ ¼ F þ f ðXÞ ð3Þ where hi stands for the expected value operator. The random spatial fluctuations of lnT and the water table fluctuations in space and time result in random spatiotemporal variability of aquifer heads /ðX; tÞ ¼ h/ðX; tÞi þ hðX; tÞ ¼ HðX; tÞ þ hðX; tÞ ð4Þ where hðX; tÞ is the zero-mean perturbation.

Substituting perturbation expansions(3) and (4)into Eq. (2), dropping all products of perturbations and subsequently taking the expected value in the resulting equation, one obtains a first-or-der approximation of the flow equation governing the mean hydraulic head @2H @X2 i H  Hp eF_f2 ¼ S eF @H @t ð5Þ

where the variability of S is assumed negligible. Subtracting this mean equation from Eq.(2)leads to the following equation @2h @X2i h  hp eF_f2 þ @f @Xi @H @Xiþ f H  Hp eF_f2 ¼ S eF @h @t fS eF @H @t ð6Þ

With reference to Eq.(5), a first-order partial differential equation relating the perturbations f and h is found as

@2h @X2 i h  hp eF_f2 þ @f @Xi @H @Xi þ f@ 2_H @X2 i " # ¼ S eF @h @t ð7Þ

If the X1direction is selected to be in the direction of the mean flow, the equations for the mean head and perturbation can then be sim-plified to @2H @X2 1 H  Hp eF_f2 ¼ S eF @H @t ð8Þ @2h @X2 i h  hp eF_f2 þ @f @X1 @H @X1 þ f@ 2_H @X2 1 " # ¼ S eF @h @t ð9Þ

Eqs.(8) and (9)provide the framework required to develop the first two moments of the hydraulic head and specific discharge in terms of the statistics of the input hydraulic parameters.

The perturbation expansion (first-order analysis) is an efficient and powerful method for solving the stochastic equation. This method is formally limited to relatively small variance (

r

2

f  1, where

r

2

f is the variance of lnT). However, Zhang and Winter

[21]found it to be accurate for the head variance solutions for

r

2 f as high as 4.38. A similar finding was reported in Gelhar[4]. 3. Spectral solutions

The approach followed is to solve Eq.(9)to fully characterize the second moments of hydraulic head and specific discharge. However, the mean Eq.(8)must be solved first in order to develop expressions for the products of mean hydraulic head gradients in Eq.(9).

3.1. Head variance

Suppose that the water table starts to fluctuate in response to periodic forcing (e.g., seasonal recharge or tides) and it would con-sist of a steady component (mean water table) plus a periodic per-turbation. As such, the mean Eq.(8)for the flow would be time invariant and the solution to Eq.(8)can be shown to be

HðX1Þ ¼ Hpþ ðH0 HpÞ coshð

g

X1Þ 

J0

g

sinhð

g

X1Þ ð10Þ where

g

¼ 1=ðfeF=2_{Þ and H}

0and J0are the reference mean head and negative mean head gradient at arbitrary location X1¼ 0, respec-tively. From Eq.(10), we immediately have

@H @X1¼

g

ðH0 HpÞ sinhð

g

X1Þ  J0coshð

g

X1Þ ð11Þ @2H @X21 ¼

g

2_ðH 0 HpÞ coshð

g

X1Þ  J0

g

sinhð

g

X1Þ ð12Þ

Superposition is a useful tool in analyzing linear groundwater problems[19,20]. The principle of superposition states that a com-plex equation can be divided into sub-equations and the solution to the original equation is then obtained by summing the individ-ual solution to each of the sub-equations. Based on this principle, the head perturbation in Eq.(9)can be separated conveniently into steady-state and time-varying components[19,20]such that hðX; tÞ ¼ hsðXÞ þ hsðX; tÞ ð13Þ This leads to separate differential equations for the steady and peri-odic components of head

@2hs @X2i  hs eF_f2þ @f @X1 @H @X1 þ f@ 2_H @X21 " # ¼ 0 ð14Þ @2hs @X2i hs hp eF_f2 ¼ S eF @hs @t ð15Þ

The solutions to Eqs.(14) and (15) can be determined using Fourier–Stieltjes representations for the perturbed quantities

[1,11,12]. By using this approach, the lnT perturbation field f is as-sumed to be a second-order stationary random field and repre-sented by the following two-dimensional wave number integral: f ðXÞ ¼

Z 1

1

eiKX_dZ

fðKÞ ð16Þ

where dZfðKÞ is the complex Fourier amplitude of lnT process, K ¼ ðK1;K2Þ is the wave number vector, and K2¼ K21þ K

2 2. The mean hydraulic head gradient is dependent of X as indicated in Eq.(11). This space-dependent mean head gradient causes the head random perturbations in Eq.(14)to be nonstation-ary. However, the head perturbed quantities in Eq.(14) can be presented using the nonstationary spectral representation[11,12]

as hsðXÞ ¼

Z 1

1

whereUhfðX; KÞ is a transfer function to be given. Thus substituting Eqs.(16) and (17)into Eq.(14)and invoking the uniqueness of the spectral representation gives the following equation

@2

U

hf

@X2 i



g

2

_U

hf¼

g

½iK1ðH0 HpÞ þ J0 sinhð

g

X1Þ þ ½

g

2ðH0 HpÞ

þ iJ0K1 coshð

g

X1ÞgeiKX ð18Þ

The corresponding solution is

U

hfðX; KÞ ¼ J0

g

R1ð

g

X1Þ K2  2K21 K4þ 4

g

2_K2 1 þ iR2ð

g

X1ÞðK 2 þ 2

g

2_ÞK 1 K4þ 4

g

2_K2 1 ( ) eiKX ð19Þ where R1ð

g

X1Þ ¼

mg

coshð

g

X1Þ  sinhð

g

X1Þ, R2ð

g

X1Þ ¼

mg

sinhð

g

X1Þ  coshð

g

X1Þ,

m

¼ ðH0 HpÞ=J0. Note that taking the advantage of a closed-form expression, the boundary effects on the perturbation hsis assumed negligible[12,14]in obtaining the solution of(19). It is expected that the boundary effect is largely limited to a small zone next to the medium boundary.

Suppose that when t > 0, the water table starts to fluctuate (as in response to a periodic forcing) according to

hpðX; tÞ ¼ hp0ðXÞ sinð

x

tÞ ð20Þ

where hp0 is the amplitude of water-level fluctuation in the leaky

phreatic aquifer and

x

is the angular frequency. We consider hp0

and hs as spatially correlated, stationary, random processes. Stationarity of the hsand hp0 processes allows the Fourier–Stieltjes

representations hsðX; tÞ ¼ Z 1 1 eiKX_dZ hsðK; tÞ ð21Þ hp0ðXÞ ¼ Z1 1 eiKX_dZ hpðKÞ ð22Þ

where dZhsand dZhp are the complex Fourier amplitudes of hs and

hp0 processes, respectively.

The transient-state part of the spectral relation follows from Eq.

(15) through the application of Eqs. (20)–(22) and the use of uniqueness of the representations:

d dtdZhsðK; tÞ þ eF Sð

g

2_{þ K}2 ÞdZhsðK; tÞ ¼ eF

_g

2 S sinð

x

tÞdZhpðKÞ ð23Þ

It is assumed that at t ¼ 0 the transient-state part of the head dis-tribution is smooth, that is, hsðX; t ¼ 0Þ ¼ 0. Thus, dZhs¼ 0 at t ¼ 0. The solution for(23)is then

dZhs ¼ eF

_g

2 S 

x

cosð

x

tÞ þ

v

sinð

x

tÞ þ

x

evt

v

2_þ

_x

2 dZhp ð24Þ where

v

¼ eF_ðK2

þ

g

2_{Þ=S. With Eq. (24), the time-varying} com-ponent of head perturbation(21)is given by

hsðX; tÞ ¼

g

2 Z 1 1 eiKX ðK2þ

g

2_Þ2_þ

_-

2 sinð

x

tÞ 

-tanð

x

tÞþ K 2 þ

g

2    þ

-

evt  dZhpðKÞ ð25Þ

where

-

¼

x

S=eF_{. The last term of the integrand in Eq.(25)which} dies away as t becomes large.

Finally, the head perturbation in Eq.(13)is determined by inte-grals(17), (19) and (25)and the corresponding solutions for the steady and time-varying components of head variance can be ex-pressed, respectively, as

r

2 hsðXÞ ¼ J 2 0

g

2_R2 1 Z 1 1 ðK2 2K21Þ 2 ðK4þ 4

g

2_K2 1Þ 2SffðKÞdK ( þ R22 Z1 1 ðK2þ 2

g

2_Þ2 K21 ðK4þ 4

g

2_K2 1Þ 2SffðKÞdK ) ð26Þ

r

2 hsðtÞ ¼

g

4Z 1 1 sinð

x

tÞ 

-tanð

x

tÞþ K 2 þ

g

2   þ

-

evt  2  ShphpðKÞ ½ðK2þ

g

2_Þ2 þ

-

2_2dK ð27Þ

where SffðKÞ is the spectrum of lnT and ShphpðKÞ is the spectrum of hp. By noting that hsand hsare statistically independent, the com-plete solution for the head variance is then given by

r

2 hðX; tÞ ¼

r

2 hsðXÞ þ

r

2 hsðtÞ ð28Þ

Eq.(28)is a local head variance relationship illustrating that the head variability is determined by the statistical properties of input hydraulic parameters.

3.2. Variance of specific discharge

The first-order equation for the specific discharge perturbation, which can be determined using Darcy’s equation, takes the form

[5] q0 i¼ TG di1JðX1Þf  @h @Xi   ð29Þ where TG¼ exp½F and

JðX1Þ ¼ 

@H @X1¼ J0

coshð

g

X1Þ 

g

ðH0 HpÞ sinhð

g

X1Þ ð30Þ

defined previously by Eq.(11). From Eqs.(13), (17), (19) and (25), the last term on the right-hand side of Eq.(29)in the X1direction can be written as @h @X1 ¼ J0 Z 1 1 ð

X

1þ i

X

2ÞeiKXdZfðKÞ þ i

g

2 Z 1 1 eiKX_K 1 ½ðK2þ

g

2_Þ2_þ

_-

2_2  sinð

x

tÞ 

-tanð

x

tÞþ K 2 þ

g

2   þ

-

evt   dZhpðKÞ ð31Þ where

X

1¼ R2ðn1Þ K2K2 1þ

g

2_K2  4

g

2_K2 1 K4þ 4

g

2_K2 1 ð32Þ

X

2¼ 2

g

R1ðn1Þ K2_K 1 K31þ

g

2K1 K4þ 4

g

2_K2 1 ð33Þ Similarly, in the X2direction

@h @X2 ¼ Z 1 1 iK2

U

hfðX; KÞdZfðKÞ þ i

g

2 Z 1 1 eiKX_K 2 ½ðK2þ

g

2_Þ2_þ

_-

2_2  sinð

x

tÞ 

-tanð

x

tÞþ K 2 þ

g

2   þ

-

evt   dZhpðKÞ ð34Þ

whereUhfis defined previously by Eq.(19).

Substituting Eqs.(31) and (34)and the Fourier–Stieltjes repre-sentations of specific discharge perturbations, i.e.,

q0 i¼

Z1

1

exp½iK  XdZqiðKÞ ð35Þ

where dZqiðKÞ is the complex Fourier–Stieltjes amplitude of q

0 i, into the specific discharge perturbation Eq. (29) and invoking the uniqueness of the spectral representation gives the following spe-cific discharge spectra in the longitudinal and transverse directions, respectively,

Sq1q1ðKÞ ¼ Sq1sðKÞ þ Sq1sðKÞ ¼ T2G J 2 0½R 2 2ð1 þ

K

1Þ2þ R21

K

2 2SffðKÞ n þ sinð

x

tÞ 

-tanð

x

tÞþ K 2 þ

g

2   þ

-

evt  2 

g

4_K2 1ShphpðKÞ ½ðK2þ

g

2_Þ2_þ

_-

2_2 ) ð36Þ Sq2q2ðKÞ ¼ Sq2sðKÞ þ Sq2sðKÞ ¼ T2G J 2 0½R 2 2

K

2 3þ R 2 1ðn1Þ

K

24SffðKÞ n þ sinð

x

tÞ 

-tanð

x

tÞþ K 2 þ

g

2   þ

-

evt  2 

g

4_K2 2ShphpðKÞ ½ðK2þ

g

2_Þ2 þ

-

2_2 ) ð37Þ where

K

1¼ K2K21þ

g

2K 2  4

g

2_K2 1 K4þ 4

g

2_K2 1 ð38aÞ

K

2¼ 2

g

ðK2 K21þ

g

2ÞK1 K4 þ 4

g

2_K2 1 ð38bÞ

K

3¼ K1ðK2þ 2

g

2ÞK2 K4þ 4

g

2_K2 1 ð39aÞ

K

4¼

g

ðK2 2K21ÞK2 K4 þ 4

g

2_K2 1 ð39bÞ The expressions in Eqs.(36) and (37)contain two principal compo-nents depending on spatial fluctuations alone and both spatial and temporal fluctuations. Within the spectral framework, the specific discharge variance can be evaluated as

r

2 qi¼

r

2 qissþ

r

2 qis¼ Z1 1 ½SqisðKÞ þ SqisðKÞdK ð40Þ

where the longitudinal and transverse components of integrand are defined by Eqs.(36) and (37), respectively.

4. Closed-form solutions

The variances of hydraulic head(28)and specific discharge(40)

requires the knowledge of spectral density functions of f and hp, we need to select them before using Eqs.(28) and (40). We consider the case where the random lnT perturbation field can be repre-sented by the following spectral density function[12,13,15]

SffðKÞ ¼ 3

r

2 f

a

2K 4

p

ðK2þ

a

2_Þ4 ð41Þ where

a

¼ 3

p

=ð16kÞ and

r

2

f and k are the variance and integral scale of lnT, respectively. Note that the choice of(41)is to meet the math-ematical requirement of zero spectral amplitude at zero wave num-ber, and, therefore, in order to produce a finite-variance head process[4]. In addition, the spectral density function to characterize the amplitude of water table fluctuation in Eq.(20)is assumed to be described by Li and Graham[13]

ShphpðKÞ ¼ 3

r

2 hp

a

2 pK 4

p

ðK2þ

a

2 pÞ 4 ð42Þ

where

a

p¼ 3

p

=ð16kpÞ and

r

2hpand kpare the variance and integral scale of the amplitude of water table fluctuation, respectively.

To simplify the analysis we assume that sufficient time has elapsed for the exponential term in Eq. (25) to die away (i.e., evt_{! 0). Because of}

_-

_{ 1, corresponding to ranges of}

_x

_{;S and} eF_{that are likely to be of interest, and the form of the denominator}

in Eq.(27), the approximate solution for the time-varying compo-nent of head variance(27)may be developed by neglecting the

-

2 term in the denominator. As a result

r

2 hsðtÞ ¼

g

4_sin2 ð

x

tÞ Z 1 1 _tanð-xtÞþ K 2 þ

g

2 h i2 ðK2þ

g

2_Þ4 ShphpðKÞdK ð43Þ

Similarly, the specific discharge spectra in the longitudinal and transverse directions (i.e., Eqs.(36) and (37)) can be approximated, respectively, as Sq1q1ðKÞ ¼ T 2 G J 2 0½R 2 2ð1 þ

K

1Þ2þ R21

K

2 2SffðKÞ þ

g

4sin2ð

x

tÞ 8 > < > : K 2 1 tanð-xtÞþ K 2 þ

g

2 h i2 ðK2þ

g

2_Þ4 ShphpðKÞ 9 > = > ; ð44Þ Sq2q2ðKÞ ¼ T 2 G J 2 0 R 2 2

K

2 3þ R 2 1ðn1Þ

K

2 4 h i SffðKÞ þ

g

4sin2ð

x

tÞ 8 > < > :  K2 2 tanð-xtÞþ K 2 þ

g

2 h i2 ðK2þ

g

2_Þ4 ShphpðKÞ 9 > = > ; ð45Þ 4.1. Head variance

Substituting Eqs. (41) and (42) into Eqs. (26) and (43) and integrating leads to the following expressions for the steady and time-varying components of head variance, depending on spatial fluctuations alone and both spatial and temporal fluctuations, respectively:

r

2 hsw¼ 16 3

p

 2 J2 0k 2

_r

2 f R 2 1ðbn1Þ 1 ð

C

2_{ 4Þ}4 1 2

C

6 þ 13

C

4þ 42

C

2þ 12   " (  12 ð

C

2  4Þ9=2

C

5 þ 3

C

3þ 8

C

4

C

ln 1 2

C

þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C

2=4  1 q  # þ R2₂ðbn1Þ 1 ð

C

2 4Þ4 1 4

C

8  5

C

6 10

C

4 4

C

2þ 12   " þ 6 ð

C

2 4Þ9=2 5

C

5_{ 2}

_C

3_{ 4}

_C

_þ8

C

  ln 1 2

C

þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C

2_{=4  1} q  #) ð46Þ

r

2 hs¼

r

2 hpsin 2 ð

x

tÞ

e

2 ð

e

 1Þ7 ð

e

 1Þ½1 þ ð8 þ 20

l

þ

l

2_Þ

_e

 þ ð18 þ 18

l

þ 29

l

2_Þ

_e

2_{þ ð8  36}

_l

_{þ 29}

_l

2_Þ

_e

3 þ ð1 þ

l

Þ2

e

4_{  6}

_e

_ln½

_e

_½ð

_e

_{ 1Þ}2 ð1 þ

e

Þ þ ð1 þ 5

e

 3

e

2_{ 3}

_e

3_Þ

_l

_{þ 2}

_e

_{ð1 þ 3}

_e

_þ

_e

2_Þ

_l

2_{ ð47Þ}

where R1 and R2 are redefined, respectively, as R1ðbn1Þ ¼

m

bcosh ðbn1Þ=k  sinhðbn1Þ, R2ðbn1Þ ¼

m

bsinhðbn1Þ=k  coshðbn1Þ, b ¼ k

g

;n1 ¼ X1=k;C¼ 3

p

=ð16bÞ,

e

¼ ½ð16=3

p

Þ

q

2;

q

¼ kp

g

;

g

¼ 1=ðfeF=2Þ and

l

¼ S

x

f2=tanð

x

tÞ.

The behavior of the dimensionless steady component of head variance in Eq.(46)as a function of dimensionless position for var-ious b is presented graphically inFig. 1a. It indicates that the head variance increases with position, while it decreases with the coef-ficient of leakage f2_{, which is inversely related to bð¼ k=ðfe}F=2_{ÞÞ, for} fixed values of eF_{and k at a fixed location. Note that a larger} coef-ficient of leakage leads to less leakage into the confined aquifer. The increase in the head variance with decreasing f at a fixed loca-tion can be explained by the fact that a decrease in f produces

more persistence of head fluctuations, which leads to larger deviations of the head from the mean head surface.Fig. 1b depicts the behavior of the dimensionless time-varying component of head variance in Eq.(47)as a function of dimensionless time. It can be clearly seen the reduction in hydraulic head variability with

increasing f for fixed values of eF _{and k}

p at a specified time. As expected, the increase in the time-varying component of head variance with kp is shown inFig. 1c. The increase in the scaled time-varying component of head variance is generally with

q

, due to either decreasing coefficient of leakage, because of greater Fig. 1. (a) Dimensionless steady component of head variance as a function of dimensionless position. Dimensionless time-varying component of head variance as a function of (b) dimensionless time and (c) dimensionless integral scale of the amplitude of water table fluctuation.

communication with the overlying water table, or due to increas-ingly spatially correlated water table fluctuations. The latter implies that water table fluctuations are either consistently above or below zero, thus contributing to a greater head variance. Finally, we note that in the limit of b ! 0 (i.e., f ! 1), corresponding to the case of no leakage, Eq.(46)converges to ½8=ð3

p

Þ2J20

r

2fk

2 , which is identical to Eq.(21)of Mizell et al.[15]. On the other hand, for the no leakage case ðf ! 1Þ the two aquifer are hydraulically dis-connected, therefore, the time-varying component of head vari-ance in Eq.(47)reduces to zero.

4.2. Variance of specific discharge

With Eqs.(41) and (42), closed-form solutions for the steady and time-varying components of specific discharge variance in the longitudinal and transverse directions are obtained by substi-tuting Eqs.(36) and (37)into Eq.(40), respectively, and integrating them over the wave number domain

r

2 q1s¼ T 2 GJ 2 0

r

2f R 2 2ðbn1Þð1 þ

U

1þ

U

2Þ þ R21ðbn1Þ

U

3 n o ð48Þ

r

2 q2s¼ T 2 GJ 2 0

r

2 f R 2 2ðbn1Þ

U

4þ R21ðbn1Þ

U

5 n o ð49Þ

r

2 q1s¼

r

2 q2s ¼ T 2 G

r

2 hpsin 2 ð

x

tÞ

g

2

_e

4ð

e

 1Þ7  ð

e

 1Þ½1 þ 2ð5 þ 2

l

Þ

e

þ 2

l

ð36 þ 11

l

Þ

e

2 þ ð26  36

l

þ 76

l

2_Þ

_e

3_{þ ð17  40}

_l

_{þ 22}

_l

2_Þ

_e

4_ þ 6

e

2_ln½

_e

_½ð

_e

_{ 1Þ}2 ð3 þ

e

Þ  2

l

ð3  3

e

þ 5

e

2_þ

_e

3_Þ þ ð1 þ 9

e

þ 9

e

2_þ

_e

3_Þ

_l

2_{ ð50Þ} where

U

1¼ 2 þ 11 4

C

2 þ 3

C

2_ln½

_C

 þ 1 ð

C

2  4Þ3  11 4

C

8 þ 29

C

6 86

C

4 16

C

2 12   þ 1 ð

C

2  4Þ7=2  3 2

C

9 þ 21

C

7  105

C

5 þ 210

C

3     ln 1 þ 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C

2=4  1 q =

C

   ln 1  2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C

2=4  1 q =

C

     168

C

þ48

C

  ln

C

2þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C

2_{=4  1} q   þ 2 ð

C

2  4Þ4ð

C

8  17

C

6 þ 118

C

4  264

C

2 Þ ð51aÞ

Fig. 2. Dimensionless steady components of specific discharge variance (a) in the longitudinal direction and (b) in the transverse direction as a function of dimensionless position. Dimensionless time-varying component of specific discharge variance in the longitudinal direction as a function of (c) dimensionless time and (d) dimensionless integral scale of the amplitude of water table fluctuation.

U2¼  11 4C 2_þ13 16C 4_{þ 3C}2_ðC2_{=2  1Þ ln½C þ} 1 ðC2  4Þ4  13 16C 12 þ57 4C 10 775 8 C 8 þ1237 4 C 6  409C4 þ 42C2 þ 222 þ12 C2   þ 1 ðC2 4Þ9=2  3 4C 13_{þ 15C}11_243 2 C 9_{þ 504C}7_2205 2 C 5_{þ 1260C}3     ln 1 þ 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C2=4  1 q =C    ln 1  2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C2=4  1 q =C      300C3_{þ 552Cþ}144 C  48 C3   ln C 2þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C2_{=4  1} q  

ð51bÞ

U3¼ C2  13 8C 2 þ33 16þ 9 4 3C 2   ln½C   þ 1 ðC2 4Þ4 13 8 C 12 401 16C 10 þ573 4 C 8 1421 4 C 6 þ 303C4  þ 59C2 120 þ12 C2  þ 1 ðC2 4Þ9=2  3 2C 13_225 8 C 11_þ837 4 C 9_3087 4 C 7_þ2835 2 C 5_{ 1260C}3     ln 1 þ 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C2_{=4  1} q =C    ln 1  2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C2_{=4  1} q =C     þ 348C3þ 384C þ24 Cþ 48 C3   ln C 2þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi C2=4  1 q   ð51cÞ

U

4¼ 11 8

C

2 13 16

C

4 þ3 2

C

2 ð1 

C

2Þ ln½

C

 þ 1 ð

C

2  4Þ4  13 16

C

12 103 8

C

10 þ615 8

C

8 829 4

C

6 þ449 2

C

4  25

C

2  52   þ 1 ð

C

2 4Þ9=2 3 4

C

13 57 4

C

11 þ 108

C

9819 2

C

7   þ1575 2

C

5  630

C

3  ln 1 þ 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C

2=4  1 q =

C

    ln 1  2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C

2_{=4  1} q =

C

  þ 72

C

3_{þ 276}

_C

_þ48

C

   ln

C

2þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C

2=4  1 q   ð52aÞ

U

5¼  33 16

C

2 þ13 8

C

4 þ 3

C

2

C

23 4   ln½

C

 þ 1 ð

C

2  4Þ4  13 8

C

12_þ401 16

C

10_573 4

C

8_þ711 2

C

6_{ 307}

_C

4_{ 47}

_C

2_{þ 12}   þ 1 ð

C

2_{ 4Þ}9=2  3 2

C

13 þ225 8

C

11 837 4

C

9 þ3087 4

C

7   2835 2

C

5 þ 1260

C

3  ln 1 þ 2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C

2=4  1 q =

C

    ln 1  2 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C

2_{=4  1} q =

C

   342

C

3 þ 324

C

48

C

   ln

C

2þ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

C

2=4  1 q   ð52bÞ Fig. 2 (continued)

Fig. 2a and b shows how the dimensionless steady components of specific discharge variance in the longitudinal and transverse directions, respectively, vary with dimensionless position for vari-ous b. These figures suggest that decreasing the coefficient of leak-age or increasing geometric mean of aquifer transmissivity increase the variability of specific discharge in the longitudinal and transverse directions at a fixed location. The increase in the dimensionless time-varying component of specific discharge vari-ance in the longitudinal direction with

q

is displayed inFig. 2c.

Fig. 2d depicts the dependence of time-varying component of spe-cific discharge variance in the longitudinal direction on the integral scale of the amplitude of water table fluctuation and indicates that the variability of time-varying component of specific discharge in-creases with the integral scale of the amplitude of water table fluc-tuation or with decreasing coefficient of leakage.

It is clear that in the case of no leakage ðb ! 0Þ, the steady com-ponents of specific discharge variance, Eqs.(48) and (49), respec-tively, reduce to

r

2 q1s¼ 3 8T 2 GJ 2 0

r

2f ð53Þ

r

2 q2s¼ 1 8T 2 GJ 2 0

r

2 f ð54Þ

which are well-known expressions for two-dimensional flow re-ported in the literature. In addition, the time-varying component of specific discharge variance reduces to zero as f ! 1.

5. Conclusions

We have analyzed groundwater flow in heterogeneous leaky confined aquifers subject to leakage from a stochastic point of view. The leakage is driven by spatial and periodic fluctuations of the water table in an overlying phreatic aquifer. The presence of leakage affects the mean hydraulic head gradient and thereby causes nonstationarity in the statistics of hydraulic head fields. Closed-form expressions for the variances of hydraulic head and specific discharge are developed in terms of statistical properties of hydraulic parameters based on the perturbation approximation and spectral Fourier–Stieltjes nonstationary representations for the perturbed quantities.

The analyses in this study are limited to small perturbations in hydraulic properties, assuming that

r

2

f is smaller than unit so that second-order terms in the flow equation can be neglected. For a lar-ger

r

2

f, the adequacy of the first-order approximation is uncertain. The results indicate that the introduction of spatial and tempo-ral variations in leakage leads to enhanced and periodic variability of the hydraulic head and of the specific discharge, which increase with distance from any arbitrary reference point. The larger the coefficient of leakage, the less variability of the hydraulic head and of specific discharge. The variability of the time-varying com-ponent of hydraulic head and of specific discharge increases with the integral scale of the amplitude of water table fluctuation.

Acknowledgements

Research leading to this work has been supported by ‘‘Aim for the Top University Plan” of the National Chiao Tung University and Ministry of Education, Taiwan. We are grateful to the five anonymous reviewers for constructive comments that improved the quality of the work.

References

[1] Bakr AA, Gelhar LW, Gutjahr AL, MacMillan JR. Stochastic analysis of spatial variability in subsurface flows: 1. Comparison of one- and three-dimensional flows. Water Resour Res 1978;14(2):263–71.

[2] Bear J. Hydraulics of groundwater. New York, NY: McGraw-Hill; 1979. [3] Chang C-M, Yeh H-D. Large-time behavior of macrodispersion in

heterogeneous trending aquifers. Water Resour Res 2007;43:W11501. doi:10.1029/2007WR006017.

[4] Gelhar LW. Stochastic subsurface hydrology. Englewood Cliffs, New Jersey: Prentice Hall; 1993.

[5] Gelhar LW, Axness CL. Three-dimensional stochastic analysis of macrodispersion in aquifers. Water Resour Res 1983;19(1):161–80. [6] Dagan G, Bellin A, Rubin Y. Lagrangian analysis of transport in heterogeneous

formations under transient conditions. Water Resour Res 1996;32(4):891–9. [7] Destouni G, Graham W. Solute transport through an integrated heterogeneous

soil–groundwater system. Water Resour Res 1995;31(8):1935–44.

[8] Hantush MM, Marino MA. One-dimensional stochastic analysis in leaky aquifers subject to random leakage. Water Resour Res 1994;30(2):549–58. [9] Indelman P, Rubin Y. Flow in heterogeneous media displaying a linear trend in

the log conductivity. Water Resour Res 1995;31(5):1257–65.

[10] Kim K, Anderson MP, Bowser CJ. Enhanced dispersion in groundwater caused by temporal changes in recharge rate and lake levels. Adv Water Resour 2000;23(6):625–35.

[11] Li S-G, McLaughlin D. A nonstationary spectral method for solving stochastic groundwater problems: unconditional analysis. Water Resour Res 1991;27(7):1589–605.

[12] Li S-G, McLaughlin D. Using the nonstationary spectral method to analyze flow through heterogeneous trending media. Water Resour Res 1995;31(3):541–51.

[13] Li L, Graham WD. Stochastic analysis of solute transport in heterogeneous aquifers subject to spatiotemporal random recharge. Water Resour Res 1999;35(4):953–71.

[14] Lu G, Zhang D. Nonstationary stochastic analysis of flow in a heterogeneous semiconfined aquifer. Water Resour Res 2002;38(8):1155. doi:10.1029/ 2001WR000546.

[15] Mizell SA, Gutjahr AL, Gelhar LW. Stochastic analysis of spatial variability in two-dimensional steady groundwater flow assuming stationary and non-stationary heads. Water Resour Res 1982;18(4):1053–67.

[16] Rehfeldt KR, Gelhar LW. Stochastic analysis of dispersion in unsteady flow in heterogeneous aquifers. Water Resour Res 1992;28(8):2085–99.

[17] Rubin Y, Bellin A. The effects of recharge on flow nonuniformity and macrodispersion. Water Resour Res 1994;30(4):939–48.

[18] Rubin Y, Seong K. Investigation of flow and transport in certain cases of nonstationary conductivity fields. Water Resour Res 1994;30(11):2901–11. [19] Townley LR. The response of aquifers to periodic forcing. Adv Water Resour

1995;18(3):125–46.

[20] Trefry MG. Periodic forcing in composite aquifers. Adv Water Resour 1999;22(6):645–56.

[21] Zhang D, Winter CL. Moment equation approach to single phase fluid flow in heterogeneous reservoirs. Soc Petrol Eng J 1999;4(2):118–27.

[22] Zhu J. Specific discharge covariance in a heterogeneous semiconfined aquifer. Water Resour Res 1998;34(8):1951–7.

[23] Zhu J, Sykes JF. Head variance and macrodispersivity tensor in a semiconfined fractal porous medium. Water Resour Res 2000;36(1):203–12.