A drawdown solution for constant-flux pumping in a confined anticline aquifer

(1)

A drawdown solution for constant-ﬂux pumping in a conﬁned anticline aquifer

Yen-Ju Chen

a

, Hund-Der Yeh

a,⇑

, Chia-Chen Kuo

b

a

Institute of Environmental Engineering, National Chiao Tung University, 1001 University Road, Hsinchu 300, Taiwan

b

National Center for High-Performance Computing, National Applied Research Laboratories, Hsinchu 300, Taiwan

a r t i c l e

i n f o

Article history:

Received 26 October 2010

Received in revised form 15 April 2011 Accepted 29 May 2011

Available online 12 June 2011

This manuscript was handled by P. Baveye, Editor-in-Chief, with the assistance of Renduo Zhang, Associate Editor Keywords:

Anticline aquifer Pumping test

Partially penetrating well Anisotropy

Integral transform

s u m m a r y

An anticline, known as a convex-upward fold in layers of rock, commonly is formed during lateral com-pression, which may be elected as a potential site for carbon sequestration. A mathematical model is developed in this study for describing the steady-state drawdown distribution in an anticline aquifer in response to the constant-flux pumping. The topographical shape of the anticline is mimicked by three successive blocks. The solution is obtained by applying the infinite Fourier transform and the finite Fou-rier cosine transform in each blocks and acquiring the hydraulic continuities between the blocks. Simu-lated results reveal that the introduction of a thin-limbs or narrow-ridged anticline would produce a much greater head drop in the ridge zone. For a well of constant pumping rate, the dimensionless draw-down around the well increases with decreasing well screen length or/and aquifer anisotropy ratio. An examination of the effect of well location on the drawdown reveals that the partially penetrating well located at the top-middle of the ridge zone produces the largest drawdown. The simulation of the flow in an anticline aquifer based on MODFLOW results in slightly smaller drawdown values in most regions when compared with those predicted by the present solution. The present solution can also be used to simulate the flow in a slab-shaped aquifer or a hillslope aquifer. It can be applied to determine the aquifer parameters if coupled with an optimization scheme and to provide the basis for selecting a potential site for carbon sequestration in the future as well.

1. Introduction

An anticline, as a result of lateral compression in crustal defor-mation, is a convex-upward fold in layers of rock. A well-struc-tured anticline formation may be chosen as a candidate site for carbon sequestration. The movement of groundwater may carry the contaminants; therefore, explicit information such as geologi-cal structure and hydrogeologigeologi-cal data are necessary to assess the applicability of the potential storage sites or to predict the migration of the contaminant plume in the site.Ashjari and Raeisi (2006) indicated that the anticline structure of aquifers and the geometry of bedrocks primarily dominate the direction of regional groundwater ﬂow in an inquiry into the groundwater ﬂow in Zagros anticlines in Iran. In recent years, several articles (e.g., För-ster et al., 2006; Bergmann et al., 2010) have been devoted to the CO2SINK project at Ketzin site in Germany. The Ketzin site is

situ-ated in the eastern part of the Roskow–Ketzin double anticline and selected to inject CO2to investigate the in situ physical, chemical,

and biological process for geological carbon sequestration. Recently, numerical or analytical solutions were developed to investigate the head responses in anticline reservoirs due to the well injection or pumping.Al-Mohannadi et al. (2007)used a

ﬁ-nite-difference method to simulate the transient pressure re-sponses to horizontal wells in anticline reservoirs and curved wells in slab reservoirs.Yeh and Kuo (2010)proposed a steady-state analytical solution for a constant-head injection at a fully pe-netrating well into a heterogeneous, anisotropic, and dome-like anticline reservoir.

In engineering practices, a constant-flux pumping test with-draws water at a constant flow rate from the test well during the test period and measures the drawdown responses in one or more observation wells in the vicinity. Generally, a drawdown solution is either incorporated with an optimization technique or applied to generate the type curves for aquifer test interpretation to deter-mine the best-fit aquifer parameters from the observed drawdown data. TheTheis equation (1935)is the most popular tool used to estimate the drawdown distribution or to determine the aquifer parameters in an inverse problem for a constant-flux pumping in a confined aquifer. It would however not be appropriate to use the well and aquifer assumptions inherent in developing these equations to describe the flow in an anticline aquifer.

The integral transform method is used widely to deal with the groundwater flow problems edged with peculiar boundaries. For example,Chan et al. (1976)used the finite Fourier transform to develop the transient and steady-state drawdown solutions for pumping in a rectangular aquifer.Chan et al. (1978)andYeh and Chang (2006) applied the finite sine transform and Hankel

⇑Corresponding author. Tel.: +886 3 5731910; fax: +886 3 5725958. E-mail address:[email protected](H.-D. Yeh).

Contents lists available atScienceDirect

Journal of Hydrology

(2)

transform to develop the transient and steady-state analytical solutions for head distribution in a wedge-shaped aquifer. On the other hand, some drawdown solutions accounting for various topography boundaries in ﬂow systems are developed based on the image-well method. The method removes aquifer boundaries and place pumping or recharging image wells at judicious loca-tions. The drawdown in an observation well is calculated by sum-ming up the drawdown or buildup due to the real well and image wells (Ferris et al., 1962; Streltsova, 1988; Chen et al., 2009).

The domain decomposition method (DDM) can be applied to handle the problem with complex geometry or mix-typed bound-ary by splitting the problem domain into several subdomains. The solutions for each subdomain are developed to satisfy the corre-sponding boundary conditions as well as the continuities of head and flux at the interface between the connected elements. The con-cept of DDM was first presented inKirkham (1957)to calculate the electrostatic potential between two concentric coaxial capped cyl-inders. The procedure was further extended inKirkham (1959)to develop the hydraulic head solution for the flow toward a partially penetrating well in a confined aquifer. Later,Javandel and Zaghi (1975)used a similar procedure to obtain the potential distribution in a confined aquifer due to the pumping at a well vertically and fully penetrating the aquifer and of radially finite extension on the bottom of the aquifer. A similar decomposition concept was also applied by Connell et al. (1998)for solving the problem of

topographically driven ﬂow in hillslope aquifers by dividing the problem domain into several rectangular elements.

The objective of this study is to develop a mathematical model for describing the steady-state drawdown distribution in response to a constant-flux pumping in an anticline aquifer. A pumping well of infinitesimal diameter partially or fully penetrates the aquifer. The anticline aquifer is homogeneous, anisotropic and confined by a curved layer on the top and a horizontal impermeable layer at the bottom. Based on the DDM, the solution of the model is ob-tained by applying the integral transform techniques including Fourier transform (FT) and finite Fourier cosine transform (FFCT) within each block. The solution is useful in predicting the spatial drawdown distribution in a wide variety of anticline aquifer sys-tems and investigating the influences of well location, screen length, aquifer geometry and anisotropy on the groundwater flow system. Moreover, the present solution may be applied to simulate the flow in hillslope and slab-shaped aquifers by assuming some of the adjacent blocks with the same heights. In addition to the ana-lytical approach, the numerical model, MODFLOW, is used to per-form simulations and the results are compared with those predicted by the present solution. The present solution can serve as an invaluable tool for gaining physical insight of the behavior of groundwater flow affected by geologic and geometric settings and for determining the aquifer parameters in an inverse problem if integrated with an optimization algorithm.

Notation

ai distance from the origin to the outer boundary of zone i

in x-direction (L)

An function of xDdeﬁned by Eq.(23)

bi height of zone i (L)

Bk function of xDdeﬁned by Eq.(27)

C1m, C2m, C4m, C6m, C7m functions of xDdeﬁned by Eqs.(41), (42),

(44), (46), and (47), respectively

C3n, C5k, C8m, C9m constants deﬁned by Eqs.(43), (45), (A22), and

(A23), respectively

D00, D0i, Dn0, Dni constants deﬁned by Eqs.(A2), (A3), (A4), and

(A5), respectively

E00, E0j, En0, Enj constants deﬁned by Eqs. (A6), (A7), (A8), and

(A9), respectively

F00, F0i, Fk0, Fki constants deﬁned by Eqs.(A10), (A11), (A12), and

(A13), respectively

G00, G0j, Gk0, Gkj constants deﬁned by Eqs.(A14), (A15), (A16), and

(A17), respectively

kx, ky, kz hydraulic conductivities in the x, y and z directions,

respectively (L/T) l screen length (L)

q volumetric pumping rate per unit length of the pumping well (L2_/T)

qD dimensionless volumetric pumping rate per unit length

of the pumping well

QD dimensionless volumetric pumping rate of the pumping

well

rD dimensionless radial distance from the pumping well to

the observation well

RD relative difference calculated by Eq.(58)

si drawdown in zone i (L)

sDi dimensionless drawdown in zone i

sDi dimensionless drawdown for zone i in Fourier domain

^

sD1P; ^sD1N dimensionless drawdown for zones 1P and 1N in Fourier

and ﬁnite Fourier cosine domain, respectively S storativity of the aquifer

S0, Sn constants deﬁned by Eqs.(A18) and (A19)

tD dimensionless time deﬁned by Eq.(57)

T transmissivity of the aquifer (L2_/T)

T0, Tk constants deﬁned by Eqs.(A20) and (A21)

up, uim, uin dimensionless variables in well functions, deﬁned by

Eqs. (54), (55), and (56), respectively U unit step function

Uc, Uc0, Ucm constants deﬁned by Eq.(32)

V0, Vn coefﬁcients in Eqs.(22), (38), and (39)

W0, Wk coefﬁcients in Eqs.(26), (38), and (39)

W well function

x0, y0, z0 coordinate of the top point of the pumping well

x0D, y0D, z0D dimensionless coordinate of the top point of the

pumping well

xD, yD, zDdimensionless coordinate variables

xDai dimensionless x-direction distance from the origin to

the outer boundary of zone i zDbi dimensionless height of zone i

zDl dimensionless screen length of the pumping well

a

n constant deﬁned by Eq.(24)

bk constant deﬁned by Eq.(28)

v

yx,

v

zx anisotropy ratios

d Dirac delta function

e

Fourier transform variable /(m, n) constant deﬁned byEq. (48)

c

m constant deﬁned by Eq.(40)

km ﬁnite Fourier cosine transform variable used with

re-spect to Eq.(1)and deﬁned by Eq.(31)

h angle between the positive xD-axis and the line

connect-ing the pumpconnect-ing and observation well #(m, k) constant deﬁned byEq. (49)

x

n ﬁnite Fourier cosine transform variable used with

re-spect to Eq.(10)and deﬁned by Eq.(25)

fk ﬁnite Fourier cosine transform variable used with

(3)

2. Methodology

2.1. Mathematical modeling of the ﬂow problem

Fig. 1sketches the configuration for a well in an anticline aqui-fer. Assume that the line sink, i.e., the pumping well of an infinites-imal radius, is extended along the z direction with length l from the point (x0, y0, z0) = (0, 0, z0). The anticline aquifer is of a finite width

in the x direction, a finite thickness in the z direction, but infinite extent in both ±y directions. In addition, the aquifer is confined, homogeneous, and anisotropic with the hydraulic conductivities of kx, kyand kzin the x, y and z directions, respectively. To simplify

the ﬂow problem, three successive blocks with different heights and widths are used to mimic the topographical shape of anticline aquifer. The height of the middle block is determined by the acme of the anticline structure, while those of the adjacent blocks are designated by the corresponding margins of the limbs. The adopted widths of the blocks should make the simulated aquifer have the same volume as the original one as possible. Furthermore, the anti-cline aquifer is decomposed into four subdomains, i.e., zones 1P, 1N, 2 and 3, according to the shapes of blocks and the well location. The mathematical model is developed in a dimensionless form to produce the simulated results in the most general way. The height of the middle block, b1, is chosen as a reference length to

nondimensionalize other variables. The dimensionless variables and parameters are deﬁned as follows: xD= x/b1, yD= y/b1, and

zD= z/b1 denoting the dimensionless coordinate variables;

x0D= x0/b1, y0D= y0/b1, and z0D= z0/b1 representing the top point

of the pumping well in the dimensionless form; xDai= ai/b1

repre-senting the dimensionless distance in x-direction of the outer boundary from the origin in zone i; zDbi= bi/b1deﬁning the

dimen-sionless height of zone i, except that zDb1= 1 standing for those in

zones 1P and 1N; sDi= si/b1denoting the dimensionless drawdown

in zone i, where the notation siis the drawdown in zone i (L);

zDl= l/b1 representing the dimensionless screen length of the

pumping well; qD= q/kxb1expressing the dimensionless

volumet-ric pumping rate per unit length of the pumping well, where the notation q is the volumetric pumping rate per unit length (L2_T1_);

_v

yx= ky/kx and

v

zx= kz/kx representing the anisotropy

ratios.

2.1.1. Formulation for ﬂow in zone 1

In the development of the mathematical model, the middle block (shown inFig. 1) is regarded as zone 1, which includes zones

1P and 1N. The steady-state groundwater ﬂow to the pumping well in zone 1 is governed by @2sD1 @x2 D þ

v

yx @2sD1 @y2 D þ

v

zx @2sD1 @z2 D ¼ qDfU½zD ðz0D zDlÞ

UðzD z0DÞgdðxD x0DÞdðyD y0DÞ; xDa1N6xD6xDa1P;

1 6 yD61; 0 6 zD61 ð1Þ

where U and d are the unit step function and Dirac delta function, respectively. The sink term in Eq.(1)implies that the ﬂux through the screen is of uniform strength. The boundary conditions at inﬁn-ity from the sink in the y direction are assumed to be

sD1ðxD;1; zDÞ ¼ 0 ð2Þ

and

@sD1ðxD;1; zDÞ

@yD

¼ 0 ð3Þ

For a conﬁned aquifer, the conditions at the top and bottom imper-meable boundaries are respectively written as

@sD1ðxD;yD;1Þ @zD ¼ 0 ð4Þ and @sD1ðxD;yD;0Þ @zD ¼ 0 ð5Þ

The continuities of ﬂux and drawdown at the right-hand edge of zone 1 are respectively given as

@sD1ðxDa1P;yD;zDÞ @xD ¼ @sD2ðxDa1P;yD;zDÞ @xD ; 0 6 zD<zDb2 ð6aÞ @sD1ðxDa1P;yD;zDÞ @xD ¼ 0; zDb26zD61 ð6bÞ and sD1ðxDa1P;yD;zDÞ ¼ sD2ðxDa1P;yD;zDÞ; 0 6 zD<zDb2 ð7Þ

Similarly, for the left-hand edge of zone 1, the following conditions should be satisﬁed: @sD1ðxDa1N;yD;zDÞ @xD ¼ @sD3ðxDa1N;yD;zDÞ @xD ; 0 6 zD<zDb3 ð8aÞ @sD1ðxDa1N;yD;zDÞ @xD ¼ 0; zDb36zD61 ð8bÞ

Fig. 1. Schematic representation of a groundwater ﬂow problem in an anticline aquifer with a line sink located along the z axis. The anticline aquifer is approximately divided into three blocks.

(4)

and

sD1ðxDa1N;yD;zDÞ ¼ sD3ðxDa1N;yD;zDÞ; 0 6 zD<zDb3 ð9Þ

The steady-state groundwater ﬂow equation in zone 2 is ex-pressed as: @2sD2 @x2 D þ

v

yx @2sD2 @y2 D þ

v

zx @2sD2 @z2 D ¼ 0; xDa1P6xD6xDa2; 1 6 yD61; 0 6 zD6zDb2 ð10Þ

The boundary conditions at inﬁnity in the ±y directions require that

sD2ðxD;1; zDÞ ¼ 0 ð11Þ

and

@sD2ðxD;1; zDÞ

@yD

¼ 0 ð12Þ

The no-ﬂow conditions hold at the top and bottom boundaries respectively as @sD2ðxD;yD;zDb2Þ @zD ¼ 0 ð13Þ and @sD2ðxD;yD;0Þ @zD ¼ 0 ð14Þ

Assume a constant-head boundary located at the lateral distance of xDa2from the pumping well; that is,

sD2ðxDa2;yD;zDÞ ¼ 0 ð15Þ

Note that Eqs.(6a) and (7), representing the continuity conditions of ﬂux and drawdown at the interface between zones 1P and 2, are the left-hand boundary conditions of zone 2.

The steady-state groundwater ﬂow equation in zone 3 is given by @2sD3 @x2 D þ

v

yx @2sD3 @y2 D þ

v

zx @2sD3 @z2 D ¼ 0; xDa36xD6xDa1N; 1 6 yD61; 0 6 zD6zDb3 ð16Þ

The boundary conditions at inﬁnite distance in the ± y directions re-quire that sD3ðxD;1; zDÞ ¼ 0 ð17Þ and @sD3ðxD;1; zDÞ @yD ¼ 0 ð18Þ

The top and bottom conditions in zone 3 are respectively given as

@sD3ðxD;yD;zDb3Þ @zD ¼ 0 ð19Þ and @sD3ðxD;yD;0Þ @zD ¼ 0 ð20Þ

A constant-head condition is applied at the lateral distance of xDa3

from the pumping well, which is described as

sD3ðxDa3;yD;zDÞ ¼ 0 ð21Þ

Furthermore, Eqs.(8a) and (9)state the continuity requirements of ﬂux and drawdown at the interface between zones 1N and 3.

2.2. Analytical solutions

2.2.1. Dimensionless drawdown solutions in Fourier domain for zones 2 and 3

To solve the partial differential Eqs.(1), (10), and (16)with their corresponding boundary conditions, the techniques of FT and FFCT are applied to the variables yDand zD, respectively, to obtain the

or-dinary differential equations (ODEs) in terms of xD. Note that

yxþ

x

2n

v

zx q ; n ¼ 0; 1; 2; 3; . . . ð24Þ

In Eq.(24),

e

is the FT transform variable;

x

n, which is the transform

variable applied to Eq. (10)in the FFCT for the integral interval [0, zDb2], is deﬁned as

x

n¼

n

p

zDb2

; n ¼ 0; 1; 2; 3; . . . ð25Þ

Note that the coefﬁcients V0and Vnare the constants needed to be

determined by the remaining boundary conditions(6a) and (7). As for zone 3, the dimensionless drawdown solution of Eq.(16)

in Fourier domain satisfying conditions(17)–(21)is expressed in the series form as:

sD3ðxD;zDÞ ¼ W0B0ðxDÞ þ X1 k¼1 WkBkðxDÞ cosðfkzDÞ ð26Þ where BkðxDÞ ¼ sinh½bkðxDa3 xDÞ

sinh½bkðxDa3 xDa1NÞ

; k ¼ 0; 1; 2; 3; . . . ð27Þ with bk¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

e

2

_v

yxþ f 2 k

v

zx q ; k ¼ 0; 1; 2; 3; . . . ð28Þ

and fk, the transform variable used in the FFCT to Eq.(16)for the

integral interval [0, zDb3], deﬁned as

fk¼

k

p

zDb3

; k ¼ 0; 1; 2; 3; . . . ð29Þ

The coefﬁcients W0and Wkin Eq.(26)are the remaining

undeter-mined constants.

2.2.2. Dimensionless drawdown solution in Fourier domain for zone 1 The FT with respect to yDand FFCT with respect to zDare

ap-plied to Eq.(1)and the result is

d2^_s_D1 dx2D ð

e

2

_v

yxþ k 2 m

v

zxÞ^sD1¼ qDUcdðxD x0DÞ; xDa1N6xD 6x_Da1P _ð30Þ

(5)

where ^sD1is the dimensionless drawdown in Fourier and ﬁnite

Fou-rier cosine domain; kmis the transform variable used in the FFCT to

Eq.(1)for the integral interval [0, 1], which is deﬁned as

km¼ m

p

; m ¼ 0; 1; 2; 3; . . . ð31Þ In addition, Uc¼ Z z0D z0DzDl cosðkmzDÞdzD ð32Þ

which can be reduced to Uc= Uc0= zDlwhen m = 0; otherwise,

Uc¼ Ucm¼

sinðkmz0DÞ sin½kmðz0D zDlÞ

km

; m ¼ 1; 2; 3; . . . ð33Þ

To solve Eq.(30)with the term of Dirac delta function, we consider the following sets of ODEs by dividing zone 1 into zones 1P and 1N as: d2^_s_D1P dx2D ð

e

2

_v

yxþ k 2 m

v

zxÞ^sD1P¼ 0; 0 < xD6xDa1P ð34Þ and d2^_s_D1N dx2D ð

e

2

_v

yxþ k 2 m

v

zxÞ^sD1N¼ 0; xDa1N6xD<0 ð35Þ

The boundary condition at xD= 0 due to the continuity requirement,

which is expressed as

^

sD1Pð0þÞ ¼ ^sD1Nð0Þ ð36Þ

Integration of Eq.(30)with respect to xDalong 0to 0+yields the

second boundary condition as

d^sD1Pð0þÞ dxD d^sD1Nð0 _Þ dxD ¼ qDUc ð37Þ

The dimensionless drawdown solutions for zones 1P and 1N in Fou-rier domain can be obtained by taking the inversion of FFCT to the solutions of Eqs.(34) and (35)with conditions(36) and (37). Apply-ingEqs. (6) and (8)to the solutions of zones 1P and 1N, respectively, yields dimensionless drawdowns as

sD1PðxD;zDÞ ¼ qDUc0 2c0 ½C10ðxDÞ þ C20ðxDÞ zDb2C30C40ðxDÞV0þ zDb3C50C60ðxDÞW0 þX 1 m¼1 qDUcm cm ½C1mðxDÞ þ C2mðxDÞ 2 cm sinðkmzDb2Þ km h i a0C30V0þP 1 n¼1 /ðm;nÞanC3nVn C4mðxDÞ þ2 cm sinðkmzDb3Þ km h i b0C50W0þP 1 k¼1 #ðm;kÞbkC5kWk C6mðxDÞ 8 > > > > > < > > > > > : 9 > > > > > = > > > > > ; cosðkmzDÞ ð38Þ in zone 1P and sD1NðxD;zDÞ ¼ qDUc0 2c0 ½C10ðxDÞ þ C70ðxDÞ zDb2C30C40ðxDÞV0þ zDb3C50C60ðxDÞW0 þX 1 m¼1 qDUcm cm ½C1mðxDÞ þ C7mðxDÞ 2 cm sinðkmzDb2Þ km h i a0C30V0þP 1 n¼1 /ðm;nÞanC3nVn C4mðxDÞ þ2 cm sinðkmzDb3Þ km h i b0C50W0þP 1 k¼1 #ðm;kÞbkC5kWk C6mðxDÞ 8 > > > > > < > > > > > : 9 > > > > > = > > > > > ; cosðkmzDÞ ð39Þ in zone 1N, where

c

m¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

e

2

_v

yxþ k 2 m

v

zx q ; m ¼ 0; 1; 2; 3; . . . ð40Þ C1mðxDÞ ¼cosh½

c

mðxDa1Pþ xDa1N xDÞ

sinh½

c

mðxDa1P xDa1NÞ

; m ¼ 0; 1; 2; 3; . . . ð41Þ C2mðxDÞ ¼

cosh½

c

mðxDa1P xDa1N xDÞ

sinh½

c

mðxDa1P xDa1NÞ

; m ¼ 0; 1; 2; 3; . . . ð42Þ C3n¼ coth½

a

nðxDa2 xDa1PÞ; n ¼ 0; 1; 2; 3; . . . ð43Þ

C4mðxDÞ ¼

cosh½

c

mðxDa1N xDÞ

sinh½

c

mðxDa1P xDa1NÞ

; m ¼ 0; 1; 2; 3; . . . ð44Þ C5k¼ coth½bkðxDa3 xDa1NÞ; k ¼ 0; 1; 2; 3; . . . ð45Þ

C6mðxDÞ ¼

cosh½

c

mðxDa1P xDÞ

sinh½

c

mðxDa1P xDa1NÞ

; m ¼ 0; 1; 2; 3; . . . ð46Þ C7mðxDÞ ¼

cosh½

c

mðxDa1P xDa1Nþ xDÞ

sinh½

c

mðxDa1P xDa1NÞ

n ð48bÞ 8 > > > < > > > : and #ðm; kÞ ¼ sin½ðkmþ fkÞzDb3 2ðkmþ fkÞ þzDb3 2 ; for km¼ fk ð49aÞ sin½ðkmþ fkÞzDb3 2ðkmþ fkÞ þsin½ðkm fkÞzDb3 2ðkm fkÞ ; for km–fk ð49bÞ 8 > > > < > > > :

Substituting Eqs.(22) and (38)into Eq.(7), the coefﬁcients V0and Vn

can be expressed in terms of the functions of V0, Vn, W0and Wkvia

the determination of the coefﬁcients in the Fourier cosine series. Similarly, the coefﬁcients W0and Wkare related to V0, Vn, W0and

Wkby substituting Eqs.(26) and (39)into Eq.(9). The coefﬁcients

V0, Vn, W0and Wkcan then be solved in the matrix form as

pre-sented inAppendix A(i.e., Eq.(A1)). 2.2.3. Inverse Fourier transform

The FT of function f(yD) with respect to the variable yDis deﬁned

as (Jeffrey and Dai, 2008):

_{f ð}

_e

_{Þ ¼}Z 1

1

f ðyDÞeieyDdyD ð50Þ

where fð

e

Þ is the transformed function and its inversion is ex-pressed as

Fig. 2. The dimensionless drawdown distributions predicted by the present solution and the image-well method (Ferris et al., 1962) for pumping at the middle of a slab-shaped aquifer bounded by two parallel constant-head boundaries.

(6)

f ðyDÞ ¼ 1 2

Þ cosð

e

yDÞd

e

ð52Þ

The numerical calculation of Eq. (52)is achieved by the routine DQDAWF ofIMSL (2003), which has the ability to cope with inte-grals of semi-inﬁnite interval and of cosine or sine integrands. 3. Results and discussion

3.1. Special cases

3.1.1. Slab-shaped aquifer bounded by parallel recharge boundaries The present solution can be simpliﬁed to describe the pumping in an isotropic slab-shaped aquifer bounded by two parallel con-stant-head boundaries along y-direction if the three blocks are of the same height, i.e., zDb1= zDb2= zDb3= 1. In this section, we

as-sume a fully penetrating well pumped at a dimensionless ﬂow rate of QD= qDzDl= 1 and located at the middle of the slab-shaped

aqui-fer with xDa2= 1 and xDa3= 1. The same problem refers toFerris

et al. (1962, Fig. 42), who illustrates the application of image-well method for the pumping in an aquifer bounded by two parallel boundaries. If an observation well is located at the dimensionless distance of rD from the pumping well, the dimensionless

draw-down at the well can be formulated by superposition ofTheis solu-tion (1935)as sDðrD;tDÞ ¼ QD 4

p

WðupÞ þ X1 m¼1 ð1ÞmWðuimÞ þ X1 n¼1 ð1ÞnWðuinÞ " # ð53Þ

where W is the well function; up, uim, and uinare the dimensionless

variables respectively deﬁned as

up¼ r2 D 4tD ð54Þ uim¼ r2 Dþ 4m2 4mrDcos h 4tD ð55Þ and uin¼ r2 Dþ 4n2 4nrDcosð

p

hÞ 4tD ð56Þ

in which h is the angle between the positive xD-axis and the line

connecting the pumping and observation wells; tDis the

dimension-less time deﬁned as

Fig. 3. Plots for the pumping at a fully penetrating well in a hillslope aquifer. The simulations were carried out by the present solution and MODFLOW for (a) dimensionless drawdown contours and (b) relative difference map in a step-like aquifer. In case (c), MODFLOW with multiple steps is used to approximate the inclined boundary.

(7)

tD¼

Tt b21S

ð57Þ

with T and S representing the transmissivity and storativity of the aquifer, respectively. Note that when tD is large enough so that

the ﬂow system reaches the steady state, the dimensionless draw-downs calculated by the image-well method can then be compared to the simpliﬁed solution.Fig. 2compares the dimensionless draw-down calculated by the present solution and Eq.(53)for pumping at the middle of the slab-shaped aquifer bounded by two parallel con-stant-head boundaries. The dimensionless drawdown are calculated along the radial direction for h = 0,

p

/4, and

p

/2. The dimensionless drawdown distributions predicted by the present solution in this special case match very well with those predicted by Eq.(53). 3.1.2. Hillslope aquifer

Fig. 3a shows the dimensionless drawdown distribution pre-dicted by the present solution and MODFLOW for ﬂow in a hill-slope conﬁned aquifer due to the pumping. The hillhill-slope aquifer is mimicked by setting two of the adjacent blocks with the same height, which has the geometry of xDa1P= 0.5, xDa1N= 0.5,

xDa2= 2, xDa3= 1, zDb1= zDb3= 1, and zDb2= 0.5. In addition, a fully

penetrating well pumped at a dimensionless pumping rate of QD= 1. In the simulation achieved by MODFLOW, the aquifer is

as-sumed bounded by two parallel constant-head boundaries with the distance of 30 m in width. The aquifer thickness varies from 10 m to 5 m in the hillslope. Note that the inﬁnite boundaries in the ±y-directions are replaced by two assigned constant-head boundaries located at ±40 m from the pumping well. The hydraulic conductivities are 104_{m/s in the x-, y- and z-directions. The}

pumping rate at the well is 102_m3_{/s so that the dimensionless}

pumping rate will be QD= 1. The length of time is set to

9.46728 107s (1095.75 day) for the steady-state simulation. The model domain has been discretized using a uniform cubic grid with a step of 0.5 m. That is, the aquifer is discretized with 21 lay-ers, 61 columns, and 161 rows.

Fig. 3a shows that the dimensionless drawdown contours are influenced by the inclination in the hillslope aquifer. The vertical flow appears around the concave corner of the top boundary, while the flow seems horizontal elsewhere. The figure indicates that a slight difference in the drawdown distribution occurs near the con-cave corner of the top boundary between the present solution and MODFLOW.Fig. 3b is the relative difference map with the relative difference values calculated by

RDð%Þ ¼sD;MODFLOW sD;present solution

sD;present solution 100 ð58Þ

InFig. 3b, the relative difference contours with negative values rep-resent that the drawdown predicted by MODFLOW is smaller than that estimated by the present solution. The largest relative differ-ence is up to 8% (absolute value) happening near the concave corner of the top boundary.Fig. 3c presents the simulated results from MODFLOW in a hillslope aquifer with more elaborate representa-tion on the sloping boundary. The significant flexure contours be-tween 0.2 and 0.6 reflect the influence of inclined top boundary on the flow pattern.

Fig. 4. (a) Dimensionless drawdown contours and (b) relative difference map on xD–zDplane with yD= 0 for pumping at a fully penetrating well in an isotropic and step-like

(8)

3.2. Effect of anticline aquifer geometry 3.2.1. Base case of anticline aquifer

To investigate the inﬂuence of aquifer geometry on the ﬂow pattern, we assume a simple case of a fully penetrating well lo-cated at an isotropic anticline aquifer. The aquifer geometry is characterized by xDa1P= 0.5, xDa1N= 0.5, xDa2= 1, xDa3= 1, and

zDb2= zDb3= 0.5. The well pumps at a dimensionless ﬂow rate of

QD= 1. Fig. 4a depicts the dimensionless drawdown distribution

in the xD–zD plane for yD= 0 predicted by the present solution

and MODFLOW. In MODFLOW simulation, the anticline aquifer is bounded by two parallel constant-head boundaries with a distance of 20 m in width. The acme of the anticline structure is 10 m in height; the limbs intersect with the two parallel constant-head boundaries at 5 m in height. The settings of the y-direction bound-aries, hydraulic conductivities, pumping rate and time length are the same as those used in the simulation of the hillslope aquifer ﬂow. The aquifer is discretized into a mesh with 21 layers, 41 col-umns and 161 rows.

Fig. 4a reveals that most of water ﬂows horizontally around the well screen and in the limbs; however, conspicuous vertical ﬂow appears around the concave corner of the top boundary. The dimensionless drawdown contours intersect with the top and bot-tom impervious boundaries at right angles. The drawdown con-tours predicted by MODFLOW seem similar to those by the present solution except that slight differences can be observed near the concave corner of the top boundary and the pumping well. InFig. 4b, the relative difference map clearly indicates that MOD-FLOW results in smaller dimensionless drawdown around the

upper part of the pumping well and limb zones than the present solution. The drawdown predicted by MODFLOW is larger than that estimated by the present solution around the bottom part of the pumping well and near the concave corner of the top boundary. The largest relative difference is up to 5 % occurring near the con-cave corner of the top boundary.

Fig. 5a compares simulations of the dimensionless drawdown contours carried out by the present solution in a step-like aquifer (base case in Fig. 4a) and MODFLOW with multiple steps to approximate the inclined boundary. The MODFLOW settings of the x-direction and y-direction boundaries, hydraulic conductivi-ties, pumping rate, time length, and grid discretization are the same as those used in the simulation ofFig. 4a.Fig. 5a shows the smooth and curved drawdown contours in the yD–zDproﬁle for

xD= 0. The contours intersected with the top boundaries at bevel

angles due to the coarse discretization on the model grid.Fig. 5b shows the relative difference map of dimensionless drawdowns predicted by MODFLOW and the present solution on the over-lapped region. The ﬁgure shows that MODFLOW gives smaller dimensionless drawdown than the present solution for the most part in the anticline and larger dimensionless drawdown at the upper parts of the limbs, when the simulation is achieved by approximating the top boundary of aquifer with multiple steps. The difference in the approximations of aquifer geometry would lead to the relative difference of the predicted dimensionless draw-downs between MODFLOW and the present solution up to 30% (absolute value). Overall, one may overestimate the dimensionless drawdown in most regions when applying the simple one-step like top boundary to simulate the anticlinal geometry.

Fig. 5. (a) Dimensionless drawdown contours and (b) relative difference map on xD–zDplane with yD= 0 for pumping at a fully penetrating well in an isotropic anticline

(9)

3.2.2. Anticline aquifers of thin limbs and narrow ridge

Fig. 6examines results obtained from ﬂow in two other aqui-fers with different geometries by considering the aquifer

por-trayed in Fig. 4a as a base case. Fig. 6a illustrates the dimensionless drawdown contours for an anticline aquifer of thin limbs. The height of the limbs are reduced to half of that in the base case, i.e., zDb2= zDb3= 0.25. Fig. 6b depicts the case of

nar-row-ridged anticline with xDa1P= 0.25 and xDa1N= 0.25. Both

fig-ures show that most of water flows horizontally around the well and in the limbs. Significantly vertical flow appears around the concave corner of the top boundary in the ridge zone due to the geometric variation. In addition, the introduction of anticline aquifers with thin limbs or narrow ridge both produces a sharp head drop in the ridge zone in comparison with that of the base case at steady state.

3.3. Effect of well partial penetration

3.3.1. Effect of screen length and aquifer anisotropy

Fig. 7 demonstrates the largest dimensionless drawdown at (xD, yD) = (0.001, 0) along the z-direction for the cases of different

dimensionless screen length of zDl= 0.2, 0.4, 0.6, 0.8 and 1.0 and

various aquifer anisotropy ratios of

v

zx= 0.3, 1.0, and 3.0. The wells

are screened from the top of the aquifer with a constant dimen-sionless pumping rate of QD= 1. Among these cases, the largest

dimensionless drawdown appears at the case of the smallest

v

zx

and zDl(

v

zx= 0.3 and zDl= 0.2) near the top of the aquifer.

More-over, the inﬂuence of aquifer anisotropy on the drawdown in-creases with the decrease of screen length.

Figs. 8a and b display the dimensionless drawdown contours for the anisotropic cases of

v

zx= 0.3 and 3, respectively. InFig. 8a,

sig-nificant vertical flow can be observed in the ridge and limb zones. The contours for sD= 0.3 to 0.8 are nearly horizontal, which reflect

Fig. 6. Plots of dimensionless drawdown contours and ﬂow ﬁelds for pumping at a fully penetrating well in an isotropic aquifer of (a) thin limbs and (b) narrow ridge.

Fig. 7. A comparison of largest dimensionless drawdown at (xD, yD) = (0.001, 0) for

the cases of different screen lengths and aquifer anisotropy ratios. The wells are screened from the top-middle of the anticline aquifer. The geometry of the aquifer is the same as the base case shown inFig. 4a.

(10)

apparent vertical ﬂow cross this region. Because kxis greater than

kzin this case (

v

zx= 0.3), most of ﬂow goes through the horizontal

path surrounding the well, leading to the larger drawdown in the upper zone 1 (i.e., 0.8 6 zD< 1). The resultant hydraulic gradient

as well as the boundary restriction causes an obvious vertical ﬂow below this zone. On the other hand,Fig. 8b represents an uncom-mon case of kzgreater than kx (

v

zx= 3). In this case, the vertical

path for flow into the well is superior to the horizontal one, which results in the contours around the well look like half ellipses. The flow in the limbs is mainly horizontal. However, obviously vertical flow can be observed in the central zone 1 below the well (zD< 0.8

and |xD| < 0.2) and around the concave corner of the top boundary

in the anticline.

3.3.2. Effect of well location

Fig. 9illustrates the inﬂuence of well location on the ﬂow pat-tern. The dimensionless screen length of the partially penetrating well is considered to be zDl= 0.2; additionally, the isotropic

anti-cline aquifer has the same geometry as the base case. Figs. 9a and b display the dimensionless drawdown contours at yD= 0 for

the pumping at a partially penetrating well located at the top-mid-dle and bottom-midtop-mid-dle of the aquifer, respectively. The flow pat-terns on the profile are symmetric to the midline of the aquifer. Most of water flows horizontally in the limbs except in the zone near the concave corner of the top boundary. Obviously vertical flows occur upward and downward in the aquifer as shown in

Figs. 9a and b, respectively, especially in the zone toward the extremity of the well. The present solution can simulate the dimensionless drawdown for an arbitrarily located pumping well in the ridge zone. In the case ofFig. 9c, the partially penetrating

well with the dimensionless screen length of 0.2 is located at a dimensionless distance of 0.25 from the midline of the anticline aquifer. The figure shows an asymmetrical flow pattern affected by the well location and aquifer geometry. Considerable vertical flow appears in the ridge zone and in the right limb, where xD< 0.35. In addition, among these three cases, the well located

at the top-middle of the aquifer, as shown inFig. 9a, has the largest dimensionless drawdown around the well.

4. Conclusions

A mathematical model has been developed for describing the steady-state flow caused by the constant-flux pumping in an anti-cline aquifer. The proposed model accounts for the flow in re-sponse to partially or fully penetrating wells of infinitesimal diameter and with uniform inflow flux along the well screen. The anticline aquifer is homogeneous, anisotropic and confined with a shape mimicked by three consecutive blocks. The integral trans-form techniques FT and FFCT are applied to develop the steady-state solutions in transform space. The coefficients in the solutions require solving a system of linear equations represented in a ma-trix form. Finally, the Fourier inversion is applied to obtain the drawdown solution in real space.

The present solution is applicable to simulate the ﬂow in a slab-shaped aquifer or a hillslope aquifer by assuming that two or three successive blocks are of the same height. For a slab-shaped aquifer, the simulated drawdown responses based on the present solution are identical to those evaluated by the image-well method when the well is fully penetrating and the aquifer is homogeneous,

iso-Fig. 8. Plots of dimensionless drawdown contours and ﬂow ﬁelds for pumping at a partially penetrating well in the aquifers with the anisotropy ratios of (a)vzx= 0.3 and (b)

(11)

tropic, confined and bounded by two parallel constant-head boundaries. Both the present solution and the numerical model, MODFLOW, are applied to simulate the case of flow in a hillslope aquifer. The grid settings allow MODFLOW to simulate the sloping boundary in a more realistic manner. Commonly, the dimension-less drawdown predicted by MODFLOW is slightly smaller than that by the present solution. In addition, the solution is used to investigate the influences of the aquifer geometry and anisotropy as well as the well partial penetration and location on the stea-dy-state flow pattern. The results obtained from these cases exhibit significant vertical flow around the concave corner of the top boundary for a fully penetrating well or a partially penetrating well located at the hump zone of the anticline. The constant-flux pump-ing in a thin-limbs or narrow-ridged anticline would cause a much

sharper head drop in the ridge zone. The influence of aquifer anisotropy on the observed drawdown cannot be ignored when the pumping is carried out in a partially penetrating well, espe-cially for the well of short open screen. When the screen length or/and the anisotropy ratio decreases, the dimensionless draw-down around the pumping well increases under the same constant pumping rate. Finally, the present solution can simulate the flow field for an arbitrarily located pumping well. In inspecting the ef-fect of well location, we find that the well located at the top-middle of the aquifer would produce larger drawdown around the well due to the boundary restriction on the anticline shape. The model MODFLOW, which can provide a better approximation on the curved boundary, is employed to simulate the flow field of the anticline aquifer. The simulated results are compared with

Fig. 9. Plots of dimensionless drawdown contours and ﬂow ﬁelds for pumping at a partially penetrating well with the dimensionless screen length of 0.2. The wells are located at (a) z0D= 1.0 and (b) z0D= 0.2 on the midline of the anticline aquifer and (c) z0D= 0.8 at a dimensionless xDdistance of 0.25 from the midline of the anticline aquifer.

(12)

the predicted results from the present solution for ﬂow toward a fully penetrating well in an anticline aquifer. MODFLOW gives slightly smaller dimensionless drawdown than the present solu-tion in most regions, while the simulasolu-tion is achieved by approxi-mating the top boundary of aquifer with multiple steps. The drawdown solution developed in this study can be further applied to identify the aquifer parameters if integrated with an optimiza-tion algorithm and to perform preliminary assessment for the selection of a potential carbon sequestration site.

Acknowledgements

The authors are grateful for support from Taiwan National Sci-ence Council under the Projects NSC99-2221-E-009-062-MY3, NSC98-3114-E-007-015, and NSC99-NU-E-009-001.

Appendix A

The coefﬁcients V0, Vn, W0and Wkin Eqs.(22), (26), (38), and

mzDb3 C1mðxDa1NÞ ðA21Þ where

C8m¼ coth½

c

mðxDa1P xDa1NÞ; m ¼ 0; 1; 2; 3; . . . ðA22Þ

and

C9m¼ csch½

c

mðxDa1P xDa1NÞ; m ¼ 0; 1; 2; 3; . . . ðA23Þ

The subroutine DLSLRG ofIMSL (2003)is used to solve Eq.(A1)by setting i = j = k = n up to 100; accordingly, 202 linear equations should be solved simultaneously.

References

Al-Mohannadi, N., Ozkan, E., Kazemi, H., 2007. Pressure-transient responses of horizontal and curved wells in anticlines and domes. SPE Reserv. Eval. Eng. 10 (1), 66–76.

Ashjari, J., Raeisi, E., 2006. Inﬂuences of anticlinal structure on regional ﬂow, Zagros, Iran. J. Cave Karst Stud. 68 (3), 118–129.

1 þ D00 D01 D02 . . . D0i E00 E01 E02 . . . E0j D10 1 þ D11 D12 . . . D1i E10 E11 E12 . . . E1j D20 D21 1 þ D22 . . . D2i E20 E21 E22 . . . E2j .. . .. . .. . .. . .. . .. . .. . .. . . . . ... Dn0 Dn1 Dn2 . . . 1 þ Dni En0 En1 En2 . . . Enj F00 F01 F02 . . . F0i 1 þ G00 G01 G02 . . . G0j F10 F11 F12 . . . F1i G10 1 þ G11 G12 . . . G1j F20 F21 F22 . . . F2i G20 G21 1 þ G22 . . . G2j .. . .. . .. . .. . .. . .. . .. . .. . . . . ... Fk0 Fk1 Fk2 . . . Fki Gk0 Gk1 Gk2 . . . 1 þ Gkj 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 V0 V1 V2 .. . Vi W0 W1 W2 .. . Wj 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ¼ S0 S1 S2 .. . Sn T0 T1 T2 .. . Tk 2 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 4 3 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 5 ðA1Þ

(13)

Bergmann, P., Lengler, U., Schmidt-Hattenberger, C., Giese, R., Norden, B., 2010. Modelling the geoelectric and seismic reservoir response caused by carbon dioxide injection based on multiphase ﬂow simulation: results from the CO2

SINK project. Chem. Erde 70 (S3), 173–183.

Chan, Y.K., Mullineux, N., Reed, J.R., 1976. Analytical solutions for drawdowns in rectangular artesian aquifers. J. Hydrol. 31, 151–160.

Chan, Y.K., Mullineux, N., Reed, J.R., Wells, G.G., 1978. Analytic solutions for drawdowns in wedge-shaped artesian aquifers. J. Hydrol. 36, 233–246. Chen, Y.J., Yeh, H.D., Yang, S.Y., 2009. Analytical solutions for constant-ﬂux and

constant-head tests at a ﬁnite-diameter well in a wedge-shaped aquifer. J. Hydraul. Eng. ASCE 135 (4), 333–337.

Connell, L.D., Jayatilaka, C., Bailey, M., 1998. A quasi-analytical solution for groundwater movement in hillslopes. J. Hydrol. 204 (1–4), 108–123. Ferris, J.G., Knowles, D.B., Brown, R.H., Stallman, R.W., 1962. Theory of Aquifer Tests,

Water-Supply Paper 1536-E, 104 p.

Förster, A., Norden, B., Zinck-Jørgensen, K., Frykman, P., Kulenkampff, J., Spangenberg, E., Erzinger, J., Zimmer, M., Kopp, J., Borm, G., Juhlin, C., Cosma, C., Hurter, S., 2006. Baseline characterization of the CO2SINK geological storage

site at Ketzin, Germany. Environ. Geosci. 13 (3), 145–161.

IMSL, 2003. IMSL Fortran Library User’s Guide Math/Library, vol. 2 of 2, Version 5.0, Visual Numerics, Houston, Texas.

Javandel, I., Zaghi, N., 1975. Analysis of ﬂow to an extended fully penetrating well. Water Resour. Res. 11 (1), 159–164.

Jeffrey, A., Dai, H.H., 2008. Handbook of Mathematical Formulas and Integrals, fourth ed. Elsevier, 541 pp.

Kirkham, D., 1957. Potential and capacity of concentric coaxial capped cylinders. J. Appl. Phys. 28 (6), 724–731.

Kirkham, D., 1959. Exact theory of ﬂow into a partially penetrating well. J. Geophys. Res. 64 (9), 1317–1327.

Streltsova, T.D., 1988. Well Testing in Heterogeneous Formations. John Wiley & Sons, New York, 413 pp.

Theis, C.V., 1935. The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using groundwater storage. Trans. Am. Geophys. Union 16 (1), 519–524.

Yeh, H.D., Chang, Y.C., 2006. New analytical solutions for groundwater ﬂow in wedge-shaped aquifers with various topographic boundary conditions. Adv. Water Resour. 29 (3), 471–480.

Yeh, H.D., Kuo, C.C., 2010. An analytical solution for heterogeneous and anisotropic anticline reservoirs under well injection. Adv. Water Resour. 33 (4), 419–429.