An analytical solution for tidal fluctuations in unconfined aquifers with a vertical beach

(1)

An analytical solution for tidal fluctuations in unconfined

aquifers with a vertical beach

H. D. Yeh,

1

C. S. Huang,

1

Y. C. Chang,

1

and D. S. Jeng

2

Received 3 October 2009; revised 27 May 2010; accepted 3 August 2010; published 22 October 2010.

[1]

The perturbation technique has been commonly used to develop analytical solutions

for simulating the dynamic response of tidal fluctuations in unconfined aquifers.

However, the solutions obtained from the perturbation method might result in poor

accuracy for the case of the perturbation parameter being not small enough. In this

paper, we develop a new analytical model for describing the water table fluctuations in

unconfined aquifers, based on Laplace and Fourier transforms. In the new approach, the

mean sea level is used as the initial condition and a free surface equation, neglecting the

second

‐order slope terms, as the upper boundary condition. Numerical results show that the

present solution agrees well with the finite different model with the second

‐order surface

terms. Unlike Teo et al.

’s (2003) approximation which restricts on the case of shallow

aquifers, the present model can be applied to most of the tidal aquifers except for the very

shallow one. In addition, a large

‐time solution in terms of sine function is provided and

examined graphically with four different tidal periods.

Citation: Yeh, H. D., C. S. Huang, Y. C. Chang, and D. S. Jeng (2010), An analytical solution for tidal fluctuations in unconfined aquifers with a vertical beach, Water Resour. Res., 46, W10535, doi:10.1029/2009WR008746.

1. Introduction

[2] The water table of the coastal unconfined aquifer

usu-ally fluctuates with periodical tides, which are in fact com-posed of several dozens of tidal components. The anisotropy of an aquifer in the horizontal and vertical directions is a consequence of long‐time sediment migration and deposi-tion. The unconfined aquifer has a free surface on the top of the aquifer system and this may introduce a vertical flow component. Most previous analytical approximations for unconfined or leaky coastal unconfined aquifer systems used a linear groundwater flow equation to approximate the unconfined flow [e.g., van der Kamp, 1972; Li and Chen, 1991; Jiao and Tang, 1999; Li et al., 2001; Li and Jiao, 2001; Jeng et al., 2002; Li et al., 2007; Xia et al., 2007; Chuang and Yeh, 2007; Sun et al., 2008]. Song et al. [2007] adopted the Boussinesq equation to derive the solution for groundwater fluctuations due to tidal waves in coastal unconfined aquifers. However, they neglected the vertical flow in their mathematical model. Li et al. [2000] developed the solution of one‐dimensional (1‐D) linearized Boussinesq equation for unconfined aquifers with moving shorelines in a sloping beach. Teo et al. [2003] presented a solution of two‐dimensional (2‐D) Laplace equation for unconfined aquifers with sloping beaches under the free surface bound-ary condition using the perturbation technique. They found that the second‐order perturbation term and beach slope are important in the development of the solution for water

table fluctuations. On the basis of the perturbation tech-nique, some researches [e.g., Dentz et al., 2006; Bolster et al., 2007; Chang et al., 2010] decoupled the flow and transport equations and developed analytical models to describe the behavior of groundwater flow in coastal aquifers.

[3] In this paper, a mathematical model, which is

com-posed of 2‐D groundwater flow equation with the mean sea level (MSL) as the initial condition and a first‐order free surface equation as the upper boundary condition, is developed for describing the water level fluctuation in tidal unconfined aquifers. The time domain solution of the model is developed by applying Laplace and Fourier sine transforms. This new analytical solution is compared with the solution presented by Teo et al. [2003] and finite dif-ference solution with the second‐order slope terms in the free surface equation. The influences of the anisotropy and specific storage on the dynamic groundwater fluctuations are investigated. In addition, the advantages and possible applications of the new solution are also addressed.

2. Theoretical Formulations

2.1. Boundary Value Problem

[4] The conceptual model for a 2‐D homogeneous,

anisotropic tidal unconfined aquifer is depicted in Figure 1, which is infinitely extended landward in horizontal direction and of a vertical beach. The origin of the coordinates is located at the intersection of the vertical beach and MSL which is considered as the reference datum. The thickness of the aquifer is D measured from the bottom of the aquifer to the MSL, and A is the amplitude of the tide.

[5] The flow is considered to be compressible,

isother-mal, and inviscid. The governing equation for describing

1_{Institute of Environmental Engineering, National Chiao Tung}

University, Hsinchu, Taiwan.

2_{Division of Civil Engineering, University of Dundee, Dundee, UK.}

(2)

the 2‐D transient head distribution (x, z, t) is [Neuman and Witherspoon, 1970; Batu, 1998]

Kx @2 @x2þ Kz @2 @z2 ¼ Ss @ @t; ð1Þ

where Kxand Kz are horizontal and vertical hydraulic

con-ductivities, respectively, and Ssis specific storage.

[6] The boundary condition for the bottom of the aquifer

is

@

@zðx; D; tÞ ¼ 0; ð2Þ

and the remote boundary condition is

lim

x!1

@

@x¼ 0: ð3Þ

The seepage face is neglected if the coastal aquifer is highly permeable. The tidal boundary is then expressed as

0; z; tð Þ ¼ A sin wtð Þ; ð4Þ

where w is tide frequency.

[7] The 2‐D equation describing the free surface for the

unconfined aquifer without a surface recharge can be written as [Batu, 1998, p. 107] Sy@8_@t ¼ Kx @8_@x 2 þ Kz @8_@z 2 Kz@8_@z at z ¼ 8: ð5Þ

Consider the case that the slope of the water table is small such that the second‐order terms of the hydraulic gradient on the right‐hand side (RHS) of equation (5) are negligible. Thus, equation (5) may be simplified as

Sy

@8 @t ¼ Kz

@8

@z at z ¼ 8: ð6Þ

Compared to the aquifer thickness (D), the tidal amplitude may be small and the actual saturated thickness can there-fore be approximated by the initial aquifer thickness (D). The head at the free surface (z =) is hence replaced by the head at the MSL (z = 0), implying that the head distribution at the MSL is the water table. The domain in z direction ranges

from 0 to−D. Equation (6) is then expressed as [Neuman, 1972]

Sy@_@t ¼ Kz@_@z at z¼ 0: ð7Þ

It is worth noting that equation (7) gives good approxima-tion whena = A/D (amplitude parameter) is small.

[8] Consider that the initial condition is equal to the MSL,

thus

x; z; 0ð Þ ¼ 0: ð8Þ

It is common to express the governing equation and boundary conditions in dimensionless form. Define the following dimensionless variables: X ¼x L; Z ¼ z D; F ¼ D; and T ¼ wt; ð9Þ

where L is the decay length defined as [Nielsen, 1990]

L¼ ffiffiffiffiffiffiffiffiffiffiffi 2KxD new r : ð10Þ

With equation (9), the governing equation (1) leads to

"2@2F @X2þ KD @2_F @Z2 ¼ as @F @T; ð11Þ

where " = D/L (perturbation parameter), KD = Kz/Kx, and

as= SwD/Kxwith S = SsD. The boundary conditions (2)–(4)

and (7)–(8) then become

@F @ZðX; 1; TÞ ¼ 0; ð12Þ lim X!1 @F @X¼ 0; ð13Þ F 0; Z; Tð Þ ¼ sin Tð Þ; ð14Þ as @F @TðX; 0; TÞ ¼ KD @F @ZðX; 0; TÞ; ð15Þ F X ; Z; 0ð Þ ¼ 0; ð16Þ wheres = S/Sy.

(3)

2.2. Previous Perturbation Solution

[9] Teo et al. [2003] presented a perturbation solution for

2‐D Laplace equation with a tidal boundary at the sloping beach and free surface boundary defined by equation (5) for isotropic unconfined aquifers. Their head solution for aquifers with a vertical beach can be expanded in powers of " and a. Using the MSL as the reference datum, the zero‐ order ("0

) solution can be expressed as

F0¼ eXcos ð Þ þ 1 2 1 4 1 e 2X þ 1 2 e pffiffi2X cos ð Þ e2 2Xcos 2ð 1Þ h i ; ð17Þ where1= T− X and 2= 2T− ffiffiffi 2 p

X. Owing to the vertical beach, the first‐order ("1_{) solution is equal to the zero‐order}

solution, i.e.,

F1¼ F0: ð18Þ

Their second‐order ("2) solution can be expressed as

F2¼ F0 ffiffiffi 2 p 3 " 2_XeX cos 1 4 þ 1 3 2_"2 _{1 þ 1 þ}X 2 e2X_2Xepffiffi2X_cos 2 4 þ epffiffi2X_sin 2 ð Þ þpffiffiffi2Xe2X cos 21 4 e2X_{sin 2} 1 ð Þ : ð19Þ

Note thatF0,F1, andF2are expanded in powers ofa, and

thus, equation (19) gives accurate approximation only when both " and a are small.

2.3. New Analytical Solution

[10] Applying Laplace and Fourier sine transforms to

equations (11)–(16) and inverting the result leads to the following solution: F X ; Z; Tð Þ ¼ 4"2 Z 1 0 w₀ð Þ þ X 1 n¼1 wnð Þ " # sin X ð Þd; ð20Þ with

w₀ð Þ ¼ c01ð Þc 02ð Þ exp cf ½02ð Þ T þ cos Tð Þ þ c02ð Þ sin T ð Þg 1 þ c02ð Þ2

h i

ð21Þ wnð Þ ¼

cn1ð Þc n2ð Þ exp c f ½ n2ð Þ T cos Tð Þ þ cn2ð Þ sin T ð Þg

1 þ cn2ð Þ2 h i ; ð22Þ where cn1ð Þ ¼ K Dncos Zð nÞ þ ascn2ð Þ sin Z ð nÞ ½ nfKD2n2½ þ 1ð Þ þ n2 þ KD"22ð þ 2n2Þ þ "44g ; ð24Þ c₀₂ð Þ ¼ KD0 2_"22 as ; ð25Þ c_n2ð Þ ¼ KDn 2_{þ "}22 as : ð26Þ

For the detailed development of equation (20), readers are referred to section A1.

[11] Theb0in equations (23) and (25) andbnin equations

(24) and (26) are the roots of following two equations:

exp 2ð 0Þ ¼ KD02þ KD0þ "22 KD02þ KD0 "22 ð27Þ tan 2ð nÞ ¼ 2KD nascn2ð Þ KD2ð2n2þ n4Þ þ 2KD"22n2þ "44; ð28Þ

respectively. Note that the variables F, X, Z, and T have been defined as the dimensionless hydraulic head, inland distance, elevation, and time, respectively, in equation (9). In addition,z is a variable in Fourier sine transform domain defined in section A1.

[12] Equation (27) can be rearranged as

"22_¼KD0 e20 1 þ KD02 e20þ 1 e20þ 1 ð Þ : ð29Þ

Substituting equation (29) into equation (25) results in

c₀₂ð Þ ¼₀ KD0 e

20 1

asðe20þ 1Þ : ð30Þ

From equation (30), since the parameters KD,s, and as are

all positive, the term e2b0_{is less than one for a negative}b 0.

On the other hand, e2b0

is greater than one for a positiveb0.

This indicates that c02(b0) is always negative for any value

ofb0.

[13] A large‐time solution can be obtained by neglecting

the exponential term in equations (21) and (22) and repla-cing T by T′ which denotes the large‐time variable since c02(z) is always less than zero. The result is

F X ; Z; T_ð 0_{Þ ¼}0 sin Tð 0þ 8Þ; ð31Þ with 0_¼pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi_a₂_{þ b}₂_; _ð32Þ c₀₁ð Þ ¼ 1 þ exp 2 1 þ Zf ½ ð Þ0g 0exp Zð 0Þ Kf Dð 20Þ þ exp 2ð 0ÞfKD½ þ 20þ 20ð1 þ 0Þ 2"22gg ; ð23Þ

(4)

a¼ 4" 2 Z 1 0 c₀₁ð Þc ₀₂ð Þ 1 þ c02ð Þ 2 þX 1 n¼1 c_n1ð Þc _n2ð Þ 1 þ cn2ð Þ 2 " # sin X ð Þd; ð33Þ b¼ 4" 2 Z 1 0 c₀₁ð Þc ₀₂ð Þ2 1 þ c02ð Þ2 þX 1 n¼1 c_n1ð Þc _n2ð Þ2 1 þ cn2ð Þ2 " # sin X ð Þd; ð34Þ ¼ cos1 a 0; ð35Þ

where is the phase shift of the head fluctuation. Equation (31) indicates that the temporal head fluctuation along inland distance is simply a harmonic motion.

2.4. Finite Difference Solution

[14] To examine the effect of neglecting the second‐order

terms in equation (5), an implicit finite difference scheme is used to approximate equations (1)–(5) and (8). The finite difference equation for equation (1) is expressed as

Kx nþ1 i1; j 2nþ1i; j þ nþ1i1; j Dx ð Þ2 þ Kz nþ1 i; jþ1 2nþ1i; j þ nþ1i; j1 Dz ð Þ2 ¼ Ss nþ1 i; j ni; j Dt ; ð36Þ

where Dx and Dz are the nodal spacings in the x and z directions, respectively,Dt is the time increment, and i, j n

is the head at the nodal point (i, j) and time level n. The ordered integer pair (i, j) is a distance (i− 1)Dx in positive x direction and ( j− 1)Dz in negative z direction. In addition,

the superscript n+1 denotes the time level of the head which is one step later than the time level n.

[15] The boundaries of (2)–(4) are approximated as nþ1 i;nz1¼ nþ1i;nzþ1; ð37Þ nþ1 nx1; j¼ nþ1nxþ1; j; ð38Þ nþ1 i¼1; j¼ A sin w t½ ðnþ DtÞ; ð39Þ

where nx and nz represent the total nodal numbers in x and z axis of the problem domain, respectively.

[16] The boundary describing the free surface at the MSL

can be approximated as Sy nþ1 i;2 ni;2 Dt ¼ Kx nþ1 iþ1;2 nþ1i;2 Dx !₂ þ Kz nþ1 i;3 nþ1i;2 Dz !₂ Kz nþ1 i;3 nþ1i;1 2Dz : ð40Þ

Furthermore, equation (40) is rearranged as

nþ1 i;1 ¼ nþ1i;3 þ 2D z Kz Sy Dt nþ1 i;2 nþ1i;2 Dt 2D z Kz Kx Dx ð Þ2 nþ1iþ1;2 nþ1i;2 ₂ 1 2Dz nþ1 i;3 nþ1i;2 2 : ð41Þ

TheDx and Dz are chosen as 1 m, and Dt is set as 0.01 min in the numerical simulations. The dynamic fluctuation is assumed to propagate to a distance of 300 m in x direction, and this distance is chosen as the infinite boundary. The convergence of the finite difference solution is discussed in section B1.

3. Results and Discussion

[17] The following sections give the comparisons of the

present solution with the results of equation (19), Teo et al.’s [2003] perturbation solution, and finite difference solution, i.e., equation (41). The effect of the anisotropy on the head distributions is investigated fors = 5 × 10−3, Z = 0, as =

8.73 × 10−4, a = 0.05, and " = 0.3. The parameter as is

defined as Sw D/Kx, and its value can be obtained with

following parameters: S = 10−3, w = 2p(12 × 60) min−1, D = 10 m, and Kx= 0.1 m min−1.

3.1. Comparison With the Finite Difference Solution [18] Figure 2 shows the spatial tidal fluctuations at T =

p/2, p, 3p/2, and 2p graphically based on the present solution, equation (20), and the finite difference solution. As shown in Figure 2, the present solution has a good agreement with the finite difference solution. This implies that the neglect of the second‐order terms in equation (5) has no significant effect on the aquifer head distribution, and the motion of the water table is mainly vertical. Define the relative difference as

D Xð Þ ¼ FNumð Þ FX Anað ÞX FNumð ÞX

100%; ð42Þ

Figure 2. Spatial head fluctuations at times 0.5p, p, 1.5p, and 2p plotted on the basis of the present solution and finite difference solution.

(5)

where FNum(X) and FAna(X) are the finite difference

solu-tion and present solusolu-tion, respectively, at inland distance X. Then, the maximum relative difference between the present solution and finite difference solution is D(0.83)ﬃ 5% at the largest tide period, i.e., T = 0.5p. This discrepancy comes from the fact that the horizontal hydraulic gradient near the free surface is larger at large‐tide period than those at other tide periods, and therefore, the neglect of second‐order terms of the hydraulic gradient in equation (5) results in small error.

3.2. Application of the Large‐Time Solution

[19] Consider that the threshold of large‐time T′ is a

dimensionless time when the sum of the exponential terms in equation (20) reaches a very small value; say, e.g., 10−4. Then those exponential terms are negligible. Figure 3 shows the threshold of large‐time T′ versus dimensionless inland distance, indicating that the threshold time increases with dimensionless inland distance. Figure 4 exhibits the head fluctuations versus dimensionless time at different dimen-sionless inland distances. Figure 4 indicates that storage coefficient plays the role of dumping effect at early time, and the head fluctuation has a tendency toward a harmonic motion because of the decrease in exponential terms. When the dimensionless time exceeds the threshold of large time, the motion of the water table is almost a harmonic motion, and the effect of storage coefficient on the head fluctuation is negligible.

3.3. Comparison WithTeo et al.’s [2003] Perturbation Solution

[20] The large‐time solution, equation (31), is compared

with Teo et al.’s [2003] perturbation solution for a vertical beach. Note that the present solution gives good accuracy

whena is not very large while their solution is accurate only for small" and a. Figures 5–8 illustrate the present solution for head distribution versus horizontal distance from the sea shore, the first‐order and second‐order solutions of Teo et al. [2003] for two different tidal periods. Figures 5–8 dem-onstrate that the difference between Teo et al.’s [2003] first‐order and second‐order solutions increases with the perturbation parameter ", indicating that the high‐order terms have the significant influence on the groundwater flow in a coastal aquifer.

[21] The finite difference solutions are included in

Figures 5–8 to compare the present solution with Teo et al.’s [2003] solution fora = 0.1 and a = 0.4, respectively, with various". In Figures 5a and 5b, the difference between the first‐order and second‐order solutions of Teo et al. [2003] is very small when" = 0.1 and 0.2 for a = 0.1. This indicates that their solution gives accurate results when both" and a are small. In addition, the present solution is of a good agreement with Teo et al.’s [2003] solution when a = 0.1 as shown in Figures 5a and 5b for" = 0.1 and 0.2, respectively. The difference between first‐order and second‐order solu-tions by Teo et al. [2003] increases with" when a = 0.1. On the other hand, the present solution gives accurate results as shown in Figures 5c and 5d even" is very large, say " ≥ 0.5, whena = 0.1.

[22] Compared Figures 5c with 6c and 5d with 6d, since

Teo et al.’s [2003] solutions are expanded in power of a, the difference between the first‐order and second‐order solu-tions increases witha for a fixed ". Therefore, their results are in poor accuracy whena is large. However, compared Figures 5a with 6a and 5b with 6b, there is no obvious difference between the first‐order and second‐order solu-tions of Teo et al. [2003] for variousa when " is fixed. This is because " gives a second‐order effect in equation (19), and thus, the difference between first‐order and second‐ order solutions by Teo et al. [2003] is small for small" even Figure 3. The threshold of large‐time T′ versus

dimen-sionless inland distance.

Figure 4. Head fluctuations versus dimensionless time at various dimensionless inland distances.

(6)

whena is large. This reflects that " has a larger effect on the accuracy of their solution than that ofa.

[23] Define the maximum absolute difference between the

present solution and Teo et al.’s solution as a function of " and a

D_maxð"; Þ ¼ Max Fj _Anað"; Þ F_Teoð"; Þj; ð43Þ

where FAna(", a) and FTeo(", a) are the hydraulic heads

predicted by the present solution and Teo et al.’s [2003] solution, respectively. The maximum absolute difference increases with" and a as shown in Figure 9. Obviously, the present model can provide more accurate predictions than Teo et al. [2003], especially when both " and a are large.

3.4. Effect of Aquifer Anisotropy

[24] The effect of aquifer anisotropy on the head

fluctu-ation is examined in Figure 10. The long‐time sediment depositional process may result in an aquifer with the smaller vertical hydraulic conductivity in comparison with the horizontal hydraulic conductivity. The ratio of the ver-tical to horizontal hydraulic conductivity (conductivity ratio)

may range from 10% to 33% [Freeze and Cherry, 1979], and therefore, the tidal aquifer is very likely anisotropic. Figure 10 shows that the anisotropy of the hydraulic con-ductivity has obvious influences on the head fluctuation near the sea shore when the conductivity ratio is significantly smaller than one and the head increases slowly in response to the rise of the sea level for small KD. It reflects that the

head has a slow response to the tide when KDis small.

3.5. Potential Applications of the Proposed Model [25] The proposed solution provides better predictions of

tide‐induced groundwater fluctuations in coastal aquifers and wider range of applications, compared with the previous perturbation approximation proposed by Teo et al. [2003]. The present analytical solutions can be applied to several practical engineering problems. For example, storm surge induced groundwater fluctuations, which normally associates with large‐wave amplitude related to the thickness of aqui-fers. In this case, the previous perturbation solution cannot provide good predictions, but the present solution will do. Since the prediction of groundwater fluctuations in coastal aquifers has significant effects on the understanding of processes of saltwater intrusion and biologic activities in coastal regions, the present model will provide an effective Figure 5. Comparison of the present solution with the result of Teo et al.’s [2003] perturbation solution

for various thickness of the tidal aquifer at one‐quarter periods (T = p/2). (a) " = 0.1, (b) " = 0.2, (c) " = 0.5, and (d)" = 1 for a = 0.1.

(7)

tool for environmental scientists who involve in the assess-ment of environassess-mental impacts in coastal regions.

4. Concluding Remarks

[26] A two‐dimensional mathematical model for

describ-ing dynamic head response in a tidal unconfined aquifer with a vertical beach has been developed. The aquifer is anisotropic, homogenous and of an impermeable bottom. In addition, a first‐order free surface equation is chosen as the upper boundary condition and the MSL as the initial con-dition for the aquifer. The analytical solution of the model is developed by using Fourier sine transforms and Laplace transforms. In addition, a large‐time solution is obtained from neglecting the exponential terms in the solution. The large‐time solution is compared with the finite difference solution and Teo et al.’s perturbation solution, where both solutions include second‐order hydraulic gradient terms at the free surface boundary. The following major conclusions can be drawn: (1) The neglect of the second‐order terms of the hydraulic gradient in the free surface equation has no obvious effect on the predicted head distribution in tidal unconfined aquifers because the motion of the water table is mainly vertical. (2) Unlike the previous approximation [Teo et al., 2003], the present model does not have any restriction on the value of". (3) The influence of the storage coefficient

on dynamic response of the tidal aquifer is very small and thus negligible. (4) The anisotropic hydraulic conductivity has significant influence on the predicted head distribution if the ratio of the vertical to horizontal hydraulic conductivity is smaller than one.

Appendix A: Derivation of Equation (20)

[27] Applying Fourier sine transforms and Laplace

transforms to the equations (11)–(16) results in an ordinary differential equation and boundary conditions in terms of Z. Solving the differential equation, one can obtain

Fsð; Z; pÞ ¼ Fs1ð; Z; pÞ þ Fs2ð; pÞ; ðA1Þ Fs1ð; Z; pÞ ¼ ffiffiffi 2 r "2_{p e} ð2þZÞ _{þ e}Z KD 2ðp2þ 1Þ f pð Þ ; ðA2Þ Fs2ð; Z; pÞ ¼ ffiffiffi 2 r "2 KD 2ðp2þ 1Þ ; ðA3Þ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi "22_{þ a} sp KD s ; ðA4Þ

Figure 6. Comparison of the present solution with the result of Teo et al.’s [2003] perturbation solution for various thickness of the tidal aquifer at one‐quarter periods (T = p/2). (a) " = 0.1, (b) " = 0.2, (c) " = 0.5, and (d)" = 1 for a = 0.4.

(8)

f pð Þ ¼ KD as þ p e2 KD as þ p; ðA5Þ

where p andz are the variables of the Laplace and Fourier sine transforms, respectively.

[28] Taking inverse Laplace transform of equation (A1)

results in Fs1þ Fs2¼ L 1_½_{G p}_{ð Þ H p}_{ð Þ}_{þ L}1 _F s2 ; ðA6Þ where G pð Þ ¼ 1 p2þ 1 ð Þ ðA7Þ H pð Þ ¼ ffiffiffi 2 r "2_{p e} ð2þZÞ _{þ e}Z KD 2f pð Þ : ðA8Þ

[29] Applying Cauchy's residue theorem to G(p) andFs2

obtains g Tð Þ ¼ L1½G pð Þ ¼ sin Tð Þ; ðA9Þ L1F_s2 ¼ ffiffiffi 2 r "2_a s "44_{þ a}_s2 exp "22_T as cos Tð Þ þ "22 as sin T ð Þ ; ðA10Þ respectively.

[30] Let h(T) be the inverse Laplace transforms of H(p)

and the result after applying the Bromwich integral is

h Tð Þ ¼ 1 2 i

Z rþi1 ri1

H pð Þ epT _dp_; _ðA11Þ

where i is an imaginary unit and r is a real constant, which is so large that all of the real parts of the poles are smaller than it. The graph of the Bromwich integral contains a close contour with a straight line parallel to the imaginary axis and two semicircles. According to Jordan's Lemma, the values of the integration for the semicircles tend to zero when radius R approaches infinity. Consequently, one can obtain

h Tð Þ ¼X

1 n¼1

Resjpn; ðA12Þ

where pnare the singularities in complex plane.

Figure 7. Comparison of the present solution with the result of Teo et al.’s [2003] perturbation solution for various thickness of the tidal aquifer at half periods (T =p). (a) " = 0.1, (b) " = 0.2, (c) " = 0.5, and (d) " = 1 for a = 0.1.

(9)

[31] There are infinite singularities in H( p) and obviously

one pole occurs at p =−"2_z2_/a

s. In order to determine other

singularities, we have f ( p) = 0 for equation (A5). There are two ways to satisfy equation (A5): one is h = b0 or p =

c02(z) and the other is h = ibnor p =−cn2(z), where both b0

and bnare real variables. Note that both c02(z) and −cn2(z)

are the singularities in terms of b0 and bn, respectively.

Substitutingh = b0and p = c02(z) into equation (A5) results

in equation (27). Similarly, substituting h = ibn and p =

−cn2(z) into equation (A5) yields equation (28). Observe

that there are infinite rootsbnin equation (28). The residues

of the singularities can be determined from the following formula [Kreyszig, 1999]: Re sjp¼pn ¼ lim_p!p n H pð Þ epT _p_p n ð Þ: ðA13Þ

Substituting equation (A8) into equation (A13), one obtain

Re sjp¼pn¼ lim_p!p n ffiffiffi 2 r "2_{p e} ð2þZÞ _{þ e}Z KD 2f pð Þ epTðp pnÞ: ðA14Þ

Figure 8. Comparison of the present solution with the result of Teo et al.’s [2003] perturbation solution for various thickness of the tidal aquifer at half periods (T =p). (a) " = 0.1, (b) " = 0.2, (c) " = 0.5, and (d) " = 1 for a = 0.4.

Figure 9. The maximum absolute difference between the present solution and Teo et al.’s [2003] solution for various " and a.

(10)

By L’Hopital’s rule, equation (A14) can be reduced to

[32] To determine the residues of the pole at p =−"2z2/as,

we substitute pn = −"2z2/as into equation (A14). One can

then find Re sjp¼"22_=a_s¼ ffiffiffi 2 r "2 as exp "22 as T : ðA16Þ

Similarly, pnis replaced by c02(z) and −cn2(z) in (A15), and

the results are

Re sjp¼c02ð Þ ¼ ffiffiffi 2 r 2"2_c 01ð Þ c 02ð Þ exp c ½02ð Þ T ðA17Þ Re sjp¼cn2ð Þ ¼ ffiffiffi 2 r 2"2cn1ð Þ c n2ð Þ exp c ½ n2ð Þ T ; ðA18Þ

respectively. Note that there are infinite singularities in −cn2(z) due to the infinite roots of bn.Therefore,

h Tð Þ ¼ Resj_p¼"22_=a_sþ Resj_p¼c₀₂_{ð Þ} þ X1 n¼1Re sjp¼cn2ð Þ ! : ðA19Þ

By the convolution theorem,

Fs1 ¼

Z T 0

gð Þh T ð Þd; ðA20Þ

where the functions g(t) and h(T − t) are defined in equations (A9) and (A19), respectively, with the

argu-Re sjp¼pn¼ lim_p!p n ffiffiffi 2 r 2 "2_{p e} ð2þZÞ _{þ e}Z _epT 3 1 þ eð 2 ÞK_D þ 2"22½_{1 þ e}2 ð_{1 þ}Þ þ 2a_s_{p 2 þ e}½ 2 ð_{2 þ þ}Þ : ðA15Þ

Figure 10. The head fluctuation versus the inland distance with various vertical hydraulic conductivities at various times.

(11)

ment T replaced by t and T − t. Equation (A20) then becomes Fs1ð; Z; TÞ ¼ ffiffiffi 2 r 2"2 _w 0ð Þ þ X1 n¼1 wnð Þ " # ffiffiffi 2 r "2_a s "44_{þ a}2 s exp "22 as T cos Tð Þ þ "22 as sin T ð Þ ; ðA21Þ

where w0(z) and wn(z) are defined in equations (21) and

(22), respectively. Note that the second term on the RHS of equation (A21) is identical to equation (A10) except that one has a minus sign. Accordingly, the solution in time domain is Fsð; Z; TÞ ¼ ffiffiffi 2 r 2"2 w₀ð Þ þ X 1 n¼1 wnð Þ " # : ðA22Þ

The final solution given as equation (20) can be obtained after taking inverse Fourier sine transform of equation (A22).

Appendix B: Convergence of the Finite Difference

Solution

[33] The validity of the numerical solution can be checked

by comparing with an analytical solution. Unfortunately, the analytical solution is not available because the problem has a nonlinear free surface boundary (i.e., equation (5)). Therefore, the convergence of the numerical solution is checked via the use of successively smaller nodal spacing. Figure B1 presents the head fluctuations graphically at times 0.5p, p, 1.5p, and 2p with different spatial and temporal increments. It demonstrates that the result of the simulation

withDx = Dz = 2 m and Dt = 0.02 min is close to that with Dx = Dz = 1 m and Dt = 0.01 min. Furthermore, the result is almost unchanged even when the spatial and temporal increments are reduced to Dx = Dz = 0.5 m and Dt = 0.005 min, respectively, indicating that the good accuracy in the numerical solution has been reached.

[34] Acknowledgments. This study was partly supported by the Taiwan National Science Council under the grant NSC 96‐2221‐E‐009‐ 087‐MY3. The authors would like to thank the Associate Editor, Xavier Sanchez‐Vila, and three anonymous reviewers for their valuable and con-structive comments.

References

Batu, V. (1998), Aquifer Hydraulics: A Comprehensive Guide to Hydro-geologic Data Analysis, John Wiley, New York.

Bolster, D. T., D. M. Tartakovsky, and M. Dentz (2007), Analytical models of contaminant transport in coastal aquifers, Adv. Water Resour., 30(9), 1962–1972, doi:10.1016/j.advwatres.2007.03.007.

Chang, Y. C., D. S. Jeng, and H. D. Yeh (2010), Tidal propagation in an oceanic island with sloping beaches, Hydrol. Earth Syst. Sci., 14, 1341–1351, doi:10.5194/hess-14-1341-2010.

Chuang, M. H., and H. D. Yeh (2007), An analytical solution for the head distribution in a tidal leaky confined aquifer extending an infinite dis-tance under the sea, Adv. Water Resour., 30(3), 439–445, doi:10.1016/ j.advwatres.2006.05.011.

Dentz, M., D. M. Tartakovsky, E. Abarca, A. Guadagnini, X. Sanchez‐ Vila, and J. Carrera (2006), Variable density flow in porous media, J. Fluid Mech., 561, 209–235, doi:10.1017/S0022112006000668.

Freeze, R. A., and J. A. Cherry (1979), Groundwater, Prentice Hall, Upper Saddle River, N. J.

Jeng, D. S., L. Li, and D. A. Barry (2002), Analytical solution for tidal propagation in a coupled semi‐confined/phreatic coastal aquifer, Adv. Water Resour., 25(5), 577–584, doi:10.1016/S0309-1708(02)00016-7. Jiao, J. J., and Z. Tang (1999), An analytical solution of groundwater

response to tidal fluctuation in a leaky confined aquifer, Water Resour. Res., 35(3), 747–751, doi:10.1029/1998WR900075.

Kreyszig, E. (1999), Advanced Engineering Mathematics, John Wiley, New York.

Li, G., and C. Chen (1991), Determining the length of confined aquifer roof extending under the sea by the tidal method, J. Hydrol., 123(1–2), 97–104. Li, H., and J. J. Jiao (2001), Tide‐induced groundwater fluctuation in a coastal leaky confined aquifer system extending under the sea, Water Resour. Res., 37(5), 1165–1171, doi:10.1029/2000WR900296. Li, H., G. Li, J. Cheng, and M. C. Boufadel (2007), Tide‐induced head

fluc-tuations in a confined aquifer with sediment covering its outlet at the sea floor, Water Resour. Res., 43, W03404, doi:10.1029/2005WR004724. Li, L., D. A. Barry, F. Stagnitti, J. Y. Parlange, and D. S. Jeng (2000),

Beach water table fluctuations due to spring‐neap tides: Moving bound-ary effects, Adv. Water Resour., 23(8), 817–824, doi:10.1016/S0309-1708(00)00017-8.

Li, L., D. A. Barry, and D. S. Jeng (2001), Tidal fluctuations in a leaky con-fined aquifer: Dynamic effects of an overlying phreatic aquifer, Water Resour. Res., 37(4), 1095–1098, doi:10.1029/2000WR900402. Neuman, S. P. (1972), Theory of flow in unconfined aquifers considering

de-layed response of the water table, Water Resour. Res., 8(4), 1031–1045, doi:10.1029/WR008i004p01031.

Neuman, S. P., and P. A. Witherspoon (1970), Variational principles for confined and unconfined flow of groundwater, Water Resour. Res., 6(5), 1376–1382, doi:10.1029/WR006i005p01376.

Nielsen, P. (1990), Tidal dynamics of the water table in beaches, Water Resour. Res., 26(9), 2127–2134.

Song, Z., L. Li, J. Kong, and H. Zhang (2007), A new analytical solution of tidal water table fluctuations in a coastal unconfined aquifer, J. Hydrol., 340(3–4), 256–260, doi:10.1016/j.jhydrol.2007.04.015.

Sun, P., H. Li, M. C. Boufadel, X. Geng, and S. Chen (2008), An analytical solution and case study of groundwater head response to dual tide in an island leaky confined aquifer, Water Resour. Res., 44, W12501, doi:10.1029/2008WR006893.

Teo, H. T., D. S. Jeng, B. R. Seymour, D. A. Barry, and L. Li (2003), A new analytical solution for water table fluctuations in coastal aquifers with sloping beaches, Adv. Water Resour., 26(12), 1239–1247, doi:10.1016/j.advwatres.2003.08.004.

Figure B1. The convergence of the finite difference solution.

(12)

van der Kamp, G. (1972), Tidal fluctuations in a confined aquifer extend-ing under the sea, in Proceedextend-ings of the 24th International Geological Congress, edited by J. E. Gill, pp. 101–106, Montreal, Que., Canada. Xia, Y., H. Li, M. C. Boufadel, Q. Guo, and G. Li (2007), Tidal wave prop-agation in a coastal aquifer: Effects of leakages through its submarine

outlet‐capping and offshore roof, J. Hydrol., 337(3–4), 249–257, doi:10.1016/j.jhydrol.2007.01.036.

Y. C. Chang, C. S. Huang, and H. D. Yeh, Institute of Environmental Engineering, National Chiao Tung University, 1001 University Rd., Hsinchu 300, Taiwan. ([email protected])

D. S. Jeng, Division of Civil Engineering, University of Dundee, Dundee DD1 4HN, UK.