RESEARCH ARTICLE
10.1002/2013WR014212
Groundwater flow to a pumping well in a sloping fault zone
unconfined aquifer
Ching-Sheng Huang1, Shaw-Yang Yang2, and Hund-Der Yeh1
1Institute of Environmental Engineering, National Chiao Tung University, Hsinchu, Taiwan,2Department of Civil
Engineering, Vanung University, Chungli, Taiwan
Abstract
This study develops a mathematical model for simulating the hydraulic head distribution in response to pumping in a sloping fault zone aquifer under a water table boundary condition. A two-dimensional equation with a sink term representing the pumping is used for describing the head distribu-tion in the aquifer. In addidistribu-tion, a first-order free surface equadistribu-tion is adopted to represent the change in water table at the outcrop. The analytical solution of the model, derived by the Laplace and finite Fourier cosine transforms, is expressed in terms of a double series. A finite difference solution within a deformable grid framework is developed to assess the solution obtained by specifying the free surface equation at the outcrop. Based on the analytical solution, we have found that the model’s prediction tends to overestimate drawdown in a late pumping period. The temporal head distribution is independent of the aquifer slope if the water table change is small, and exhibits a double-humped shape due to the effect of the free surface. The temporal drawdown predicted from the analytical solution is further compared with those measured from a pumping test conducted in northern Portugal.1. Introduction
There are two primary types of sloping unconfined aquifers: one is an aquifer with the water table as the upper boundary and a sloping impervious bed at the bottom, while the other is a sloping fault zone aquifer with the water table at the outcrop. The former involves water table change when subject to rainfall recharge. Chapuis [2011] reviewed analytical solutions for the problem involving steady state groundwater seepage due to recharge. He assessed the validation of those solutions using a finite element solution con-sidering both saturated and unsaturated flows. On the other hand, the second type of aquifer, also called a sloping fracture, represents an aquifer in a rock fracture or a fault that connects the atmosphere with the free surface at the outcrop. The permeability of the surrounding rock in that case is generally much lower than that of the aquifer.
Analytical solutions which may be applicable to the sloping fault zone aquifer are reviewed as follows. Han-tush [1962] considered a wedge-shaped aquifer connected to an external reservoir, where the water level varies temporally. He presented an analytical solution describing hydraulic head in the aquifer and indicated that the wedge-shaped aquifer cannot be approximated by a uniform one. Latinopoulos [1985] considered a rectangular aquifer where each side can be subject to either the Dirichlet, no-flow, or Robin boundary con-dition. He developed analytical solutions describing hydraulic head in the aquifer with different combina-tions of the boundary condicombina-tions. Pacheco [2002] considered a sloping fault zone aquifer extending semi-infinitely from the free surface at the outcrop. He used Cooper and Jacob [1946] solution to describe draw-down in the aquifer before the drawdraw-down cone reaches the free surface. Free surface movement is approxi-mated as a known function of time [Pacheco, 2002, Equation (7), p. 119] once the drawdown cone reaches the free surface. Nevertheless, such an approximation neglects the effect of gravity drainage from the free surface decline.
This current paper develops a mathematical model for describing two-dimensional (2-D) groundwater flow induced by pumping in a sloping fault zone aquifer with the water table at the outcrop. A first-order free surface equation is employed to describe a water table decline. A sink term in the 2-D flow equation repre-sents a constant pumping rate for a vertical well. The analytical solution of the model, expressed in terms of a double series, is developed by applying the Laplace transform and finite Fourier cosine transform. The finite difference solution of the model within a deformable grid framework is also developed for
Key Points:
An analytical solution for head in a sloping unconfined aquifer is developed
The effect of the aquifer slope on temporal and spatial head is addressed
The simplified free surface equation specified at the outcrop is assessed
Correspondence to:
H.-D. Yeh,
Citation:
Huang, C.-S., S.-Y. Yang, and H.-D. Yeh (2014), Groundwater flow to a pumping well in a sloping fault zone unconfined aquifer, Water Resour. Res., 50, 4079–4094, doi:10.1002/ 2013WR014212.
Received 3 JUNE 2013 Accepted 30 APR 2014
Accepted article online 5 MAY 2014 Published online 20 MAY 2014
Water Resources Research
comparison with the analytical solution. The effect of the aquifer slope on temporal and spatial head distri-butions is investigated. Additionally, the temporal drawdown predicted by the analytical solution is com-pared with the field drawdown data observed in a pumping test by Pacheco [2002].
2. Mathematical Development
2.1. Conceptual Model
Rocks surrounding a sloping fault zone unconfined aquifer are regarded as an impermeable stratum when the ratio of the rock’s hydraulic conductivity to that of the aquifer is less than 1029[Huang et al., 2012]. Fig-ure 1 illustrates the schematic diagram of such a sloping fault zone unconfined aquifer with a pumping well. This aquifer extends finitely from the ground, outcrops with a free surface, slants with an angle h, and has a thickness B as shown in Figure 1a. The x and y axes are horizontal, and the z axis is vertical in the Car-tesian coordinate system. The x0and z0axes are parallel and perpendicular to the aquifer, respectively, in the sloping coordinate system. The well position is (0, y0) as shown in Figure 1b. The aquifer has finite Wx
and Wyin width in the x0-direction and y-direction, respectively. The distance between the well and the free
surface is W1in the x0-direction and W1cos h in the x-direction.
Consider 2-D groundwater flow in the sloping aquifer. The flow equation describing transient hydraulic head hðx; y; tÞ in response to constant pumping can be expressed as
T@ 2h @x021T @2h @y25S @h @t1Qd x 0 d y2yð 0Þ (1)
where T is the aquifer transmissivity, S is the aquifer storage coefficient, Q is a pumping rate, and d() is the Dirac delta function. Equation (1) is applicable to 2-D flow in the sloping aquifer except in the case that the aquifer is horizontal or vertical. The case of the horizontal aquifer is trivial. For the vertical aquifer, one may refer to Anderson [2006].
The Dirac delta functions in equation (1) represent an infinitesimal well radius with negligible effects of a finite well radius and wellbore storage on the head. Papadopulos and Cooper [1967] mentioned that those effects greatly diminish when t > 2:53102r2
c=T , where rcis the inner radius of a well. In addition, Yeh and
Chang [2013] also stated that the effects are negligible for a small-diameter well with a finite radius varying between 0.05 and 0.25 m.
The flow is static prior to pumping, and the hydraulic head is uniform over the whole domain. With the ref-erence datum chosen at the free surface, the initial condition is therefore expressed as
h50 at t50 (2)
That is, the atmospheric pressure is set to zero. Based on equation (2), the hydraulic head h is negative for pumping, and its absolute value represents drawdown.
The first-order free surface equation describing vertical water table movement induced by pumping is writ-ten as
K@h @z52Sy
@h
@t at z5W1sin h1h and xo x xf (3)
where K5T/B is the aquifer hydraulic conductivity, Syis the specific yield, xo5W1cos h2B=ð2sin hÞ, and xf5W1cos h1B=ð2sin hÞ. Introduce the relation between the Cartesian and sloping coordinate systems as
x05xcos h1zsin h y05y z052xsin h1zcos h
(4)
Based on equation (4), equation (3) can be written as
K @h @x0 @x0 @z1 @h @z0 @z0 @z ! 52Sy @h @t (5)
where @x0=@z5sin h. The derivative term @h=@z0approaches zero since the groundwater moves along the x0-direction. Under this condition, equation (5) reduces to
Ksin h@h @x052Sy @h @t at x 0 5W11h=sin h (6)
Note that W11h=sin h represents a water table position in x0-direction and makes equation (6) nonlinear due to the presence of the unknown water table position. Equation (6) can be linearized by neglecting the last term h=sin h as
Ksin h@h @x052Sy @h @t at x 0 5W1 (7)
which implies that the boundary condition is fixed at x05W1. This treatment is similar to the one taken by Neuman [1972] in developing an analytical solution for a horizontal unconfined aquifer.
The edges of the aquifer are considered under the no-flow condition as
@h=@x050 at x052W2 (8)
@h=@y50 at y56Wy=2 (9)
where W2is the distance between the well and the bottom edge shown in Figure 1a. The beginning time
accounting for the boundary effect on the head can be estimated by R2S=ð16TÞ where R represents a short-est distance measured from the well to the boundaries of the aquifer [Wang and Yeh, 2008, Table II]. In
addition, each series of the model’s solution, equations (14) and (15), converges faster when setting smaller widths of Wxand Wy.
Define dimensionless variables and parameters below
hD5 Th Q; tD5 Tt SW12 ; x0D5 x0 W1 ; xD5 x W1 ; yD5 y W1 ; y0D5 y0 W1 ; r5 SyB SW1sin h ; s5 Q TW1sin h v5W2 W1 ; a5Wx W1 ; n5Wy W1 (10)
where the subscript D is used to denote the dimensionless variables. Based on equation (10), equations (1) and (2) and (6)–(9) can be written, respectively, as
@2h D @x02D1 @2h D @y2 D 5@hD @tD 1d x0D d yð D2y0DÞ (11) hD50 at tD50 (11a) @hD @x0 D 52r@hD @tD at x0D511shD (11b) @hD @x0 D 52r@hD @tD at x0D51 (11c) @hD=@x0D50 at x0D52v (11d)
@hD=@yD50 at yD56n=2 (11e)
2.2. Analytical Solution
Applying the finite Fourier cosine transform to yDand the Laplace transform to tDin equations (11) and
(11a) and (11c)-11e) results in an ordinary differential equation (ODE) in terms of x0D. Solving the ODE simul-taneously with the transformed boundary conditions yields the head solution in the Laplace and Fourier domain as hD1ðx 0 D;wm;pÞ5uðv; 12x 0 DÞ for 0 x 0 D 1 (12) hD2ðx 0 D;wm;pÞ5uðv1x 0 D;1Þ for 2v x 0 D 0 (13) with
uða; bÞ52coshðakÞ½kcosh ðbkÞ1rpsinh ðbkÞ
pk½rpcosh ðakÞ1ksinh ðakÞ cos½wmðy0D1 n 2Þ (13a) k5 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi w2 m1p=# q (13b) wm5pm=n (13c)
where p is the Laplace transform parameter, and m僆 1, 2, . . .1 is the finite Fourier cosine transform parameter. A detailed derivation of equations (12) and (13) is given in Appendix A. Application of complex analysis to the inverse Laplace transform leads to the head solution in the Fourier domain illustrated in Appendix B. Introducing the formula, equation (A3), for the inverse finite Fourier cosine transform to the Fourier-domain solution yields the head solution as
hD1ðx 0 D;yD;tDÞ5Uðv; 12x 0 DÞ for 0 x 0 D 1 (14)
hD2ðx 0 D;yD;tDÞ5Uðv1x 0 D;1Þ for 2v x 0 D 0 (15) with Uða; bÞ51 n utða; bÞ1X 1 n51 unða; b; 0Þ12X 1 m51
usða; b; wmÞ1u0ða; b; wmÞ
½ YðmÞ 12X 1 m51 X1 n51 unða; b; wmÞYðmÞ 8 > > > > < > > > > : 9 > > > > = > > > > ; (15a) utða; bÞ52½3l1ða21b2Þ16l1ðtD1rbÞ2l3=ð6l21Þ (15b)
usða; b; wÞ52cosh ðawÞcosh ðbwÞ=½#wsinh ðawÞ (15c)
u0ða; b; wÞ52cosh ðab0Þexp ½2k0ðwÞtD½b0coshðbb0Þ2rk0ðwÞsinh ðbb0Þ=g0ðwÞ (15d) unða; b; wÞ52cos ðabnÞexp ½2knðwÞtD½bncosðbbnÞ2rknðwÞsin ðbbnÞ=gnðwÞ (15e)
YðmÞ5cos ½pmð0:51yD=nÞcos ½pmð0:51y0D=nÞ (15f)
g0ðwÞ5k0ðwÞfl2b0coshðab0Þ1½12rak0ðwÞsinh ðab0Þg (15g) gnðwÞ5knðwÞfl2bncosðabnÞ1½12raknðwÞsin ðabnÞg (15h)
k0ðwÞ5w22b20;knðwÞ5w21b2n (15i)
l15a1r; l25a12r; l35a2ða13rÞ (15j)
where wmis defined in equation (13c) and b0and bnare the roots of equations denoted, respectively, as
expð2ab0Þ5 2rb201b01rw2 m rb201b02rw2 m (16) and tanðabnÞ52rðb2 n1w2mÞ=bn (17)
Notice that equation (16) has only one positive root, whereas equation (17) has infinite positive roots caused by the periodic function of tanðabnÞ. Estimates for b0and bnare given in section 2.3. The solution,
equations (14) and (15), is composed of four terms. The first term is the closed form expression for /t; the
second and third terms contain the simple series expanded in terms of bnand m, respectively; the last term
is the double series expanded in terms of m and bn. Equation (15b) is a first-order polynomial in time,
indi-cating that there is no steady state head distribution. The water table position can be approximated as
x0Dapprox:511shD1 (18)
where hD1, a function of yDand tD, is the head predicted from equation (14) with x
0
D51. The absolute value of hD1represents drawdown. The term hD1=sin h reflects water table movement along the negative direc-tion of x0D-axis. Thus, the problem domain falls in the range of x
0
D x
0
Dapprox:. The x
0
Dapprox:will be com-pared with 11shDwhere hDis predicted by the finite difference solution developed based on equation
(11b) in section 2.4. This comparison will be illustrated in section 3.3.
2.3. Calculation of b0and bn
The roots in equations (16) and (17), b0and bn, can be found by Newton’s method with appropriate initial
guess values. The roots are the intersection points of the left-hand side (LHS) and right-hand side (RHS) functions of equation (16) for b0and (17) for bn. The root b0is very close to the vertical asymptote of the
RHS function of equation (16). The initial guess for b0is chosen as the position of the asymptote derived by
setting the denominator of the RHS function to zero, and in turn expressed as
b0initial5d1ð211
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 114ðrwmÞ2 q
Þ=ð2rÞ (19)
where d is a small positive value, say d51028, to prevent the denominator from being zero in the iteration process. Similarly, the roots bnare also close to the vertical asymptotes of the LHS function, tanðabnÞ. The initial guesses for bnare therefore chosen as
bninitial5d1ð2n21Þp=ð2aÞ (20) where n僆 1, 2, . . .1.
2.4. Finite Difference Solution
An implicit finite difference method is applied to approximate equations (11)–(11b) and (11d) and (11e) for comparison with the present analytical solution developed based on equation (11c). A nonuniform finite difference grid is used to discretize the problem domain, and its nodal points are shown in Figure 2. Assume that the aquifer has a dimensionless width of 10 in the x0D-direction and yD-direction, and the
problem domain falls in the range of 29 x0D 1 and 25 yD 5. The pumping well is located at the origin (x0D50 and yD50). The region near the well has small grid sizes and away from the well has large grid sizes.
Equation (11) is approximated as
hDn11i11;j22hDn11i;j 1hDn11i21;j ðDx0
Dk;jÞ2 1
hDn11i;j1122hDn11i;j 1hDn11i;j21 ðDyDi;kÞ2 5hD
n11 i;j 2hDni;j
DtD 1
½ðDx0D1;jÞðDyDi;1Þ21at well position
0 else where
8 < :
(21)
where DtDis a dimensionless time step, tD5DtD3n is the present time, hDn11i;j is the hydraulic head at nodal point (i, j), and time tD, hDni;jis the head at one step earlier than hDn11i;j , and k僆 1, 2, and 3. The integer pair (i, j) is ordered from the LHS boundary in the x0D-direction and from the lower boundary in the yD-direction.
The grid sizes Dx0D1;j, Dx0D2;j, and Dx0D3;jin the x0D-direction are 0:02 for 20:2 x
0 D 0:2; 0:1 for 21 x0D 20:2 and 0:2 x 0 D 1, and 0:4 for 29 x 0
D 21, respectively, as shown in Figure 2. Similarly, the grid sizes DyDi;1, DyDi;2, and DyDi;3in yD-direction are 0:02 for 20:2 yD 0:2; 0:1 for 21 yD 2 0:2 and 0:2 yD 1, and 0:4 for 25 yD 21 and 1 yD 5, respectively.
The initial condition of equation (11a) can be denoted as
hD1i;j50 at each ði; jÞ (22)
Equation (11d) representing the no-flow boundary at the first i is approximated as
hDn11i21;j5hDn11i11;j at i51 (23) Similarly, equation (11e) representing the no-flow boundaries for the first and last j is approximated, respec-tively, as
hDn11i;j215hDn11i;j11at j51 (24) hDn11i;j215hDn11i;j11at j5ny (25) where nydenotes the total number of nodes in yD-direction, and (0, j), (i, 0), and (i, ny11) are the positions
of fictitious nodes outside the no-flow boundaries.
Consider a deformable grid framework in the region of 0:2 x0D 1 as shown in Figure 2, where water table movement is described by equation (11b). In this region, the grid size remains constant before the movement. The grid sizes DyDi;1, DyDi;2, and DyDi;3remain constant while the grid size Dx0D2;jdecreases to match a new location of the water table. The deformable grid size Dx0Dn11j is then adjusted according to the water table location at the previous time step as
Dx0Dn11j 5
Wd1shDni;j nd
at i5nx (26)
where nxis the total number of the grids in x
0
D-direction, Wd(the region width) is 0.8, and nd(the number
of the grids inside the region) is 8. Note that the grid size Dx0Dn11j will be changed at each time step, and the nxand ndmaintain constant during the entire simulation time. Equation (11b), which describes the
water table movement, can be approximated as
hDn11 i;j 2hDn11i21;j Dx0 Dn11j 52rhD n11 i;j 2hDni;j DtD at i5nx (27)
where Dx0Dn11j has been defined by equation (26). In the numerical simulation, the time increment DtDis set to 1, and the simulation period is 0 tD 1460 (365 day).
3. Results and Discussion
This section analyzes hydraulic head predicted by the analytical solution and the finite difference solution. In section 3.1, the influence of the aquifer slope on temporal and spatial head distributions is assessed. In section 3.2, the effects of the no-flow and water table boundaries on a temporal head distribution is investi-gated. In section 3.3, the use of equation (11c) to develop the analytical solution is evaluated and discussed. In these three sections, default values of variables and hydraulic parameters are given in the second column of Table 1. In section 3.4, the application of the present solution to a real-world case is presented.
3.1. Effect of Aquifer Slope on Head
The aquifer slope affects late parts of a temporal head distribution. Figure 3 illustrates the temporal distribu-tions of hydraulic head hD2(20.2, 0, tD) and hD2(21, 0, tD) predicted by equation (15) when h 5 20, 45, and
60. The hydraulic head hD1(1, 0, tD) predicted by equation (14) is also considered, and its absolute valuej
hD1j is used to approximate a water table decline as illustrated in Figure 3. When tD<6, there is no differ-ence in predicted head hD2for various h. When tD 6, the water table represented by jhD1j declines
signifi-cantly, and a larger h leads to a smaller head hD2at a specific time. It seems reasonable to conclude that the
change in the aquifer slope does not affect the temporal head distributions before the occurrence of a sig-nificant water table decline.
Figure 4 shows the spatial head distribution predicted by equation (14) for h 5 0, 20, 45, and 60. We con-sider tD50.04 for avoiding the boundary effect. The Theis [1935] solution for a horizontal confined aquifer
is taken for comparison. The figure indicates that the spatial head distributions plotted in the horizontal coordinate system for varied h are different. The reason for difference from them is because the varying value of h causes different x0Dby the relationship x
0
D5xD=cos h for a fixed xD. Accordingly, the Theis [1935]
solution for a horizontal confined aquifer can be used to account for the slope effect prior to the boundary effect when replacing the radial variable r by r=cos h. The Theis [1935] solution gives the same predicted result as the present analytical solution for each h.
3.2. Boundary Effect on Temporal Head
The free surface, described by equation (11c), causes temporally varied flow in a double-humped shape, which exhibits unconfined behavior. Figure 5 demonstrates the temporal head distributions of hD1(1, 0, tD)
predicted by equation (14) and hD2(21, 0, tD) predicted by equation (15) for r 5 40 (Sy50.1) and 80
(Sy50.2). The absolute value of the headjhD1j denotes a vertical water table decline, as illustrated in Figure
5. That figure shows that the temporal head distribution of hD2reaches its flat stage in the period of
4 tD 10, when r540 or 4 tD 20 when r580. The flat stage is due to gravity drainage at x0D51, and
Table 1. Default Values to Variables and Hydraulic Parameters
Notation Default Value (unit)
Field Data (unit)
Description
E3 E5
h none none none Hydraulic head
(x0, y, t) (250 m, 0, 0.01 day) (0.07m, 0, 280 min) (0.07m, 0, 110 min) Variables of sloping Cartesian coordinate,
and time variable
y0 0 0 0 Well position in y-direction
B 30 m 35 m 35 m Aquifer thickness in z0-direction
T 30 m2 /d 1.2 3 1023 m2 /min 3.5 3 1024 m2 /min Transmissivity S 3 3 1023 6.1 3 1024 6.1 3 1024 Storage coefficient Q 100 m3/d 0.0342 m3/min 0.0354 m3/min Pumping rate
K 1 m/d 3.4 3 1025
m/min 1.0 3 1025
m/min Hydraulic conductivity defined as K 5 T/B (Sy, h) (0.2, p/6) (0.01, p/4) (0.01, p/4) Aquifer specific yield and slope, respectively
(W1, W2) (50 m, 450 m) (178 m, 322 m) (184 m, 316 m) Distance measured from the well to the top and
bottom boundaries, respectively, in x0-direction
(Wx, Wy) (500 m, 500 m) (500 m, 500 m) (500 m, 500 m) Aquifer widths in x0-direction and y-direction,
respectively, and Wx5W11W2
(hD1, hD2) none none none Th/Q for 0 x
0 D 1 and 2v x 0 D 0, respectively ðx0D;yD;tDÞ (21, 0, 0.04) (3.9 3 1024, 0, 0.017) (3.8 3 1024, 0, 1.9 3 1023) ðx 0 =W1;y=W1;Tt=ðSW21ÞÞ y0D 0 0 0 y0/W1 (r, s) (80, 0.133) (4.56, 0.226) (4.41, 0.777) ðSyB=ðSW1sin hÞ; Q=ðTW1sin hÞÞ (v, a, n) (9, 10, 10) (1.81, 2.81, 2.81) (1.72, 2.72, 2.71) (W2/W1, Wx/W1, Wy/W1)
a larger r leads to a longer flat stage. The flat stage occurs at tD54, with a significant water table decline
represented byjhD1j. On the other hand, the no-flow boundary at the edges of the aquifer makes the tem-poral head distribution decrease dramatically. Figure 6 shows the temtem-poral distributions of hD2(21, 0, tD) for
the aquifer widths of a 5 n 5 10 (500 m) and 20 (1000 m). Both curves deviate at tD520, indicating the
boundary effect on the head distribution in the case of a 5 n 5 10. Note that the effect of the no-flow boundary is more apparent at the later period than that of the water table boundary because the pumping well is closer to the outcrop, as shown in Figure 6.
3.3. Free Surface Equation
The finite difference solution is used to assess a water table decline predicted based on equation (11c), which is commonly used to develop analytical solutions [e.g., Neuman, 1972; Zhan and Zlotnik, 2002; Huang et al., 2012]. Figure 7 illustrates the spatial head distributions predicted by the analytical solution and finite difference solution at times of 4 (1 day), 400 (100 day), and 1460 (365 day). The numerical solution considers the free surface boundary, described by equation (11b) specified at the dynamic water table. On the other hand, the analytical solution is developed based on the free surface equation fixed at the outcrop (i.e., x0D51). Consequently, the location of the water table can be approximated by equation (18). The water table declines predicted by both solutions are along the aquifer slope, as shown in the figure. When tD54,
the predicted head distributions from both solutions agree well and are near symmetrical to the well because of the insignificant water table decline at xD50:86. The figure however indicates that the analytical solution overestimates the drawdown at the late pumping times of tD5400 and 1460. Furthermore, the
temporal head distributions at the locations of (26.6, 0) and (23, 0) and at the outcrop of (1, 0) predicted by both solutions are drawn in Figure 8. The results predicted by the analytical solution are valid before
Figure 3. Temporal distributions of hydraulic head at the observation wells of (21, 0) and (20.2, 0) and at the outcrop of (1, 0) predicted by the analytical solution for various aquifer slopes.
tD510 that is 2.5 day calculated based on t5ðSW12ÞtD=T for the values of T, S, and W1(550 m) given in
Table 1. After tD510, the analytical solution underestimates the hydraulic head (or overestimates the
draw-down) because of the simplification of equation (11b) to (11c) by specifying the free surface equation at the outcrop (i.e., x0D51). However, the location to specify this equation should follow the water table, which moves toward the well as time increases. In engineering practice, a large distance between a pumping well and the outcrop should be seriously considered for much groundwater exploitation. Accordingly, the ana-lytical solution can give fairly good prediction before t 5 100 day for the given default values except that W15320 m.
3.4. Comparison With Field Data
Pacheco [2002] reported a pumping test conducted in a geologic formation consisting mostly of serpentine rocks splitting with several fractures in the study region. Crushed materials with clay minerals fill these frac-tures. The fractures distribute anomalously and crisscross with each other. Eleven pumping wells were installed and named as serial numbers from E1 to E11. One fracture with E3 and E5 wells corresponds with the present model considering a single fracture with a pumping well [Pacheco, 2002, FZ3 in Figure 5, p. 124]. The pumping rates at E3 and E5 wells are 0.0342 and 0.0354 m3/min, respectively. The radius for both wells is 0.07 m. The drawdown data observed from E3 and E5 wells are shown in Figure 9a. The fracture inclines approximately 45toward the southwest, and its average thickness is 35 m. A schematic diagram of the fracture with the wells is shown in Figure 9b.
Pacheco [2002] indicated that the values of T and S range from 8.4 3 1025to 1.2 3 1022m2/min and from 2 3 1025to 6.1 3 1022, respectively. Accounting for the aquifer heterogeneity, the values of T are chosen as 1.2 3 1023and 3.5 3 1024m2/min, and the values of K are 3.4 3 1025and 1.0 3 1025m/min for the local aquifer near E3 and E5 wells, respectively. The S is set as 6.1 3 1024and Syis set as 0.01 because of
the presence of the clay formation. The hydraulic parameters mentioned above are shown in Table 1.
The time-drawdown curves predicted by equation (14) are compared with those predicted by the Pacheco [2002] solution and the field data observed at E3 and E5 wells [Pacheco, 2002] shown in Figure 9a. The E3 and E5 wells operated for 280 and 110 min, respectively. The drawdown cones induced by E3 and E5 wells have not reached the outcrop at t 5 280 and 110 min, respectively, as illustrated in Figure 9b, implying that the flow is still under the confined condition. Figure 9a reveals that both solutions overestimate the draw-down over the whole pumping period. This discrepancy may be attributed to water recharge from adjacent small fractures and a nearby fault zone [Pacheco, 2002, FZ4 in Figure 5, p. 124] mentioned in Pacheco and Van der Weijden [2012]. In addition, the difference in the predicted drawdown curves is mainly because the Pacheco solution is developed based on the Cooper-Jacob formula.
4. Concluding Remarks
A mathematical model has been developed for simulating groundwater flow due to pumping in a sloping fault zone unconfined aquifer outcropping with a free surface. A 2-D equation with a sink term representing a pumping well is used for describing the flow in the sloping aquifer. A first-order free surface equation specified at the outcrop describes water table movement while the no-flow boundary is imposed for the other edges of the aquifer. The analytical solution of the model, derived by the Laplace transform and finite Fourier cosine transform, is expressed as a double series expanded in terms of integers as well as roots b0
and bn. The roots are determined by Newton’s method with the initial guesses represented by analytical
expressions. Furthermore, the finite difference solution of the model is developed within a deformable grid framework to simulate the position of the dynamic water table at the outcrop. Additionally, the drawdown distribution predicted by the analytical solution is compared with field data observed from a pumping test reported in Pacheco [2002]. Main conclusions providing new physical insights can be drawn below:
1. The aquifer slope does not affect a temporal head distribution before a significant water table decline tak-ing place.
Figure 6. Temporal distributions of hydraulic head at the observation well of (21, 0) predicted by the analytical solution for different aquifer widths.
2. The spatial head distribution in the horizontal coordinate system depends on the aquifer slope because of the relationship of x0D5xD=cos h.
3. The Theis [1935] solution, replacing r by r=cos h, is applicable to evaluating the effect of the slope on a spatial head distribution before a significant water table decline taking place.
Figure 8. Temporal distributions of hydraulic head at the observation wells of (26.6, 0) and (23, 0) and at the outcrop of (1, 0) predicted by the analytical solution and finite difference solution.
Figure 9. Comparison of the temporal drawdown distributions predicted by the analytical solution and Pacheco [2002] solution with the field drawdown data reported by Pacheco [2002].
4. The time-drawdown curve for a sloping fault zone unconfined aquifer exhibits a double-humped shape due to gravity drainage once the water table declines.
5. The predicted head distribution obtained based on the assumption that the free surface equation is fixed at the outcrop is valid for a short dimensionless pumping period and tends to overestimate drawdown for a long one.
Appendix A: Derivation of Equations (12) and (13)
The definition of the finite Fourier cosine transform is expressed as
hDðmÞ5 ðn=2 2n=2 hDcos wmðyD1 n 2Þ dyD (A1)
where wmis defined by equation (13c). The transform has the property of
ðn 0 @2h D @y2 D cosðwmyDÞdyD5ð21Þm@hD @yD y D5n=2 2@hD @yD y D52n=2 2wm2hDðmÞ (A2)
The first and second RHS terms equal zero due to equation (11e). The formula for the inverse finite Fourier cosine transform can be written as
hD5 1 n hDð0Þ1 2 n X1 m51 hDðmÞcos wmðyD1 n 2Þ (A3)
On the other hand, the definition of the Laplace transform is expressed as
hDðpÞ5 ð1
0
hDexpð2ptDÞdtD (A4)
It has the property of
ð1 0
@hD @tD
expð2ptDÞdtD5phD2hDjtD50 (A5)
where hDis the head transformed by the finite Fourier cosine transform. The second RHS term equals zero because of equation (11a).
After applying these two transforms to equations (11), (11a), and (11c)-11e), the resultant ODE and trans-formed boundary conditions in terms of x0Dare rewritten as
@2hD @x02 D 2ðp1w2 mÞhD5 1 pcos wmðy0D1 n 2Þ d x0D (A6) @2hD @x0 D 52rphDat x 0 D51 (A7) @hD=@x 0 D50 at x 0 D52v (A8)
where p is the Laplace transform parameter. Equation (A6) can be separated according to the Dirac delta function d x 0D into two homogeneous ODEs as
@2hD @x02 D 2ðp1w2 mÞhD150 for 0 x 0 D 1 (A9) @2hD @x02 D 2ðp1wm2ÞhD250 for 2v x0D 0 (A10)
Two requirements at x0D50 are needed to solve the ODEs. One is the continuity of the transformed hydrau-lic head denoted as
hD15hD2at x
0
D50 (A11)
The other can be obtained by integrating equation (A6) to x0Dfrom x
0 D502to x 0 D501as @2hD1 @x0 D 2@ 2h D2 @x0 D 51 pcos wmðy0D1 n 2Þ at x0D50 (A12)
which reflects the discontinuity of the flux at x0D50 due to d x
0
D
. Note that the integration of the second LHS term in equation (A6) equals zero because of equation (A11). Solving equations (A9) and (A10) simulta-neously with equations (A7), (A8), (A11), and (A12) leads to equations (12) and (13).
Appendix B: Inverse Laplace Transform of Equation (13a)
The function / defined by equation (13a) is a value function of p in a complex plane. The single-value function gives the only result to a specific p and has no branch cut that reflects spatial discontinuity between p1and p2. Therefore, uðp1Þ equals uðp2Þ for any complex number p. One can confirm uðp1Þ5u ðp2Þ in the following way. Let p1and p2be in terms of the polar coordinate with the origin at p52w2
m obtained from the root of kðpÞ50 as
p15rexpðidÞ2w2
m (B1)
p25rexp½iðd22pÞ2w2
m (B2)
where r is a radius from the origin, d is an argument between 0 and 2p, and i is the imaginary unit. Note that equations (B1) and (B(2)) have the same result in terms of a complex number. Substituting equations (B1) and (B2) into equation (13b), respectively, yields
k5pffiffirexpðid=2Þ5pffiffir½cos ðd=2Þ1isin ðd=2Þ (B3) k5pffiffirexp½iðd22pÞ=252pffiffir½cos ðd=2Þ1isin ðd=2Þ (B4) Notably, equations (B3) and (B4) differ in the negative sign. The results of substituting equations (B3) and (B4) separately into equation (13a) are the same, indicating uðp1Þ5uðp2Þ and no branch cut.
Since / is a single-value function, an integral contour for its inverse Laplace transform consists of a straight line and a semicircle. The straight line extends from p5r2i1 to p5r1i1 where r herein represents a posi-tive value to contain all of the poles of the function /. The semicircle has an infinite radius and connects the straight line to enclose those poles. The integral along the semicircle is zero because of the infinite radius. According to the residue theory, the time domain inversion equals the sum of the residues of the inside poles. Each residue can be estimated by the following formula
Res5 1 ðm21Þ!limp!} dm21 dpm21 uðpÞ3ðp2}Þ m ½ (B5)
where } represents the location of the pole of the function /, and m is the order of the pole. The locations of the poles are determined by setting the denominator of equation (13a) to be zero as
pk½rpcosh ðakÞ1ksinh ðakÞ50 (B6) where k5 ffiffiffiffiffiffiffiffiffiffiffiffiffiw2
m1p p
. The roots of the equation represent the locations of the poles in a complex plane and appear only at the negative part of the real axis. Note that p52w2
mderived from k50 is not a pole because the power of ffiffiffiffiffiffiffiffiffiffiffiffiffiw2
m1p p
is 1/2 rather than an integer. Equation (B6) will be analyzed according to the value of wm.
When wm>0, equation (B6) obviously has one simple pole at p 5 0. According to equation (B5), the residue equals the product of cos½wmðy0D1n=2Þ and /sdefined by equation (15c). In addition, other poles can be
obtained numerically by
rp coshðakÞ1k sinh ðakÞ50 (B7)
which is derived from setting the term in the bracket of equation (B6) to be zero. One simple pole is located at p 5 p0between p 5 0 and p52w2m. Let k5
ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi w2
m1p0 p
5b0for conciseness, and thus one obtains p05b202w2
m. Substituting k5b0and p5p0into equation (B7) and rearranging the result yields equation (16). With the value of b0, the location of the simple pole can be obtained by p05b202wm2. Its residue is derived by applying equation (B5) as the product of cos½wmðy0D1n=2Þ and /0defined by equation (15d). On the
other hand, behind p52w2
m, there are infinite simple poles at p5pnwhere n僆 1, 2, 3, . . .1. For preventing
the residues of these poles from the imagery unit i, we let k5 ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiw2 m1pn p
5ibnand have pn52b2n2w2 m. Substi-tuting k5ibnand p5pninto equation (B7) and rearranging the result leads to equation (17). With the values of bn, the locations of those simple poles can be estimated by pn52b2n2wm2. The residues are expressed by applying equation (B5) as the product of cos½wmðy0D1n=2Þ and /ndefined by equation (15e).
When wm50, equation (B6) reduces to
p2½rkcosh ðakÞ1sinh ðakÞ50 (B8)
where k5 ffiffiffipp. Apparently, one second-order pole is at p 5 0. Applying equation (B5) results in its residue as /tdefined by equation (15b). Furthermore, infinite simple poles are at p 5 pnbehind p 5 0. The residues
can be obtained directly by substituting wm50 into equation (15e).
References
Anderson, E. I. (2006), Analytical solutions for flow to a well through a fault, Adv. Water Resour., 29(12), 1790–1803, doi:10.1016/ j.advwatres.2005.12.010.
Chapuis, R. P. (2011), Steady state groundwater seepage in sloping unconfined aquifers, Bull. Eng. Geol. Environ., 70(1), 89–99, doi:10.1007/ s10064-010-0282-2.
Cooper, H. H., and C. E. Jacob (1946), A generalized graphical method for evaluating the formation constants and summarizing well field history, Trans. AGU, 27, 526–534.
Hantush, M. S. (1962), Flow of ground water in sands of nonuniform thickness Part 1. Flow in a wedge-shaped aquifer, J. Geophys. Res., 67(2), 703–709.
Huang, C.-S., H.-D. Yeh, and C.-H. Chang (2012), A general analytical solution for groundwater fluctuations due to dual tide in long but nar-row islands, Water Resour. Res., 48, W05508, doi:10.1029/2011WR011211.
Latinopoulos, P. (1985), Analytical solutions for periodic well recharge in rectangular aquifers with third-kind boundary conditions, J. Hydrol., 77, 293–306.
Neuman, S. P. (1972), Theory of flow in unconfined aquifers considering delayed response of the water table, Water Resour. Res., 8(4), 1031–1045.
Pacheco, F. A. L. (2002), Response to pumping of wells in sloping fault zone aquifers, J. Hydrol., 259(1-4), 116–135, doi:10.1016/S0022-1694(01)00584-4.
Pacheco, F. A. L., and C. H. Van der Weijden (2012), Weathering of plagioclase across variable flow and solute transport regimes, J. Hydrol., 420-421, 46–58, doi:10.1016/j.jhydrol.2011.11.044.
Papadopulos, I. S., and H. H. Cooper (1967), Drawdown in a well of large diameter, Water Resour. Res., 3(1), 241–244.
Theis, C. V. (1935), The relation between the lowering of the piezometric surface and the rate and duration of discharge of a well using ground-water storage, Trans. AGU, 16, 519–524.
Wang, C. T., and H. D. Yeh (2008), Obtaining the steady-state drawdown solutions of constant-head and constant-flux tests, Hydrol. Proc-esses, 22, 3456–3461, doi:10.1002/hyp.6950.
Yeh, H. D., and Y. C. Chang (2013), Recent advances in modeling of well hydraulics, Adv. Water Resour., 51, 27–51, doi:10.1016/ j.advwatres.2012.03.006.
Zhan, H., and V. A. Zlotnik (2002), Groundwater flow to a horizontal or slanted well in an unconfined aquifer, Water Resour. Res., 38(7), 1108, doi:10.1029/2001WR000401.
Acknowledgments
The authors are very much grateful to the esteemed reviewers for their valuable comments and suggestions. Research leading to this paper has been partially supported by the grants from Taiwan National Science Council under the contract NSC 101-2221-E-009-105-MY2 and NSC 102-2221-E-009-072-MY2.
