Tide-induced groundwater level fluctuation in a U-shaped coastal aquifer

Journal of Hydrology 530 (2015) 291–305

Contents lists available at ScienceDirect

Journal of Hydrology journal homepage: www.elsevier.com/locate/jhydrol

Tide-induced groundwater level fluctuation in a U-shaped coastal aquifer Fu-Kuo Huang a, Mo-Hsiung Chuang b, Grace S. Wang c, Hund-Der Yeh d,⇑ a

Department of Water Resources and Environmental Engineering, Tamkang University, 151 Yingzhuan Rd., Tamsui Dist., New Taipei City 251, Taiwan Department of Urban Design and Disaster Management, Ming-Chuan University, 5 Deming Rd., Guishan Dist., Taoyuan City 333, Taiwan c Department of Construction Engineering, Chaoyang University of Technology, 168 Jifeng E. Rd., Wufeng Dist., Taichung City 413, Taiwan d Institute of Environmental Engineering, National Chiao Tung University, 1001 University Rd., Hsinchu 300, Taiwan b

a r t i c l e

i n f o

Article history: Received 26 May 2015 Received in revised form 19 August 2015 Accepted 14 September 2015 Available online 1 October 2015 This manuscript was handled by Peter K. Kitanidis, Editor-in-Chief, with the assistance of Jian Luo, Associate Editor Keywords: Coastal aquifer U-shaped coastline Shoreline length Tidal wave interaction Sensitivity analysis

s u m m a r y The prediction of groundwater level fluctuation due to tidal waves propagation in coastal aquifers is important for the planning and management of water resources in coastal areas. A two-dimensional (2-D) analytical solution is derived to describe the tidal groundwater fluctuation in an aquifer bounded by three water–land boundaries that form a U-shaped coastline. Two opposite sides represent estuary– land boundaries on which the amplitude attenuation and phase shift of the tidal movement in the estuary are considered while the third side is an ocean–land boundary. The effects of wave interaction due to the propagation of oscillating oceanic tides in the cross-shore direction inland and the transmission of the two opposite estuarine tides in the along-shore direction are investigated. Three existing headfluctuation solutions can be considered as special cases of the present solution; one is for onedimensional flow and the other two are for 2-D flow. A transition distance ranging from 10 to 15 times of tidal propagation length along the shoreline can be estimated based on the solution. This distance can be used to judge whether the interaction among tides is significant. The influences of hydraulic properties on the tidal fluctuations within the aquifer can therefore be assessed quantitatively. Based on sensitivity analyses, one can conclude that the tidal head is most sensitive to the transmissivity and storativity of the aquifer, and least to the damping coefficient of tidal amplitude and wave number along the estuary. The sensitivities of head fluctuation to the changes of transmissivity and storativity depend on the shoreline length and whether the interaction among waves is significant. On the other hand, the sensitivities of head fluctuation to the changes of damping coefficient and wave number increase with diagonal distance from the entry of estuary and reach the largest magnitude near the estuary far away seashore. Ó 2015 Elsevier B.V. All rights reserved.

1. Introduction Tide-induced groundwater fluctuations in coastal aquifers have received much attention in past few decades. Dynamic behavior of the groundwater flow plays an important role in many coastal hydrological, environmental, and engineering problems. They are, for example, coastal aquifer characterization, saltwater intrusion, marine contamination, and stability of slopes and retaining structures along the coast. Since the 1950s, many researchers have carried out a variety of studies on those topics to explore associated phenomena and solve related problems. These research works

⇑ Corresponding author. E-mail addresses: [email protected] (F.-K. Huang), [email protected]. edu.tw (M.-H. Chuang), [email protected] (G.S. Wang), [email protected] (H.-D. Yeh). http://dx.doi.org/10.1016/j.jhydrol.2015.09.032 0022-1694/Ó 2015 Elsevier B.V. All rights reserved.

may be classified according to the aquifer configuration or type as: single layer aquifers (e.g., Jacob, 1950; Nielsen, 1990; Farrell, 1994; Sun, 1997; Li et al., 1999); aquifers with L-shaped boundary (Li et al., 2000, 2002; Li and Jiao, 2002) or with sinusoidal coastline (Jeng et al., 2005); submarine aquifers (e.g., Chuang and Yeh, 2007; Guo et al., 2007; Li et al., 2007, 2008; Xia et al., 2007; Geng et al., 2009; Dong et al., 2012; Guarracino et al., 2012; Wang et al., 2014); aquifers with sloping beaches (Teo et al., 2003) or with sloping bottom (Asadi-Aghbolaghi et al., 2012); island aquifers (Trefry and Bekele, 2004; Li et al., 2006; Sun et al., 2008; Huang et al., 2012) and inhomogeneous tidal aquifers (Chuang et al., 2010; Asadi-Aghbolaghi et al., 2014). Among those, Jacob (1950) developed a one-dimensional (1-D) analytical equation to describe the transmission of groundwater head fluctuation, in response to tidal oscillations in a coastal confined aquifer. His equation has been widely employed to evaluate field data in many case studies

292

F.-K. Huang et al. / Journal of Hydrology 530 (2015) 291–305

Nomenclature A Am Bm D E F Fs F
1 s H H1 H2 H2 Im K1() KL KR L Li Ln Lp NS O Pi Re S Si T U Um V Vm X Y

amplitude of tidal wave along ocean shore [L] coefficient in Eq. (A15) [–] coefficient in Eq. (A15) [–] =T/S, hydraulic diffusivity of the aquifer [L2 T
1] coefficient in Eq. (A29) [–] coefficient in Eq. (A29) [–] operator of Fourier sine transform [–] operator of inverse Fourier sine transform [–] dimensionless tidal head [–] component of dimensionless tidal head H in Eq. (A1) [–] component of dimensionless tidal head H in Eq. (A1) [–] value of Fourier sine transform for H2 [–] imaginary part of complex function modified second-kind Bessel function of first order =k1L in Eq. (A30) [–] =k2L in Eq. (A31) [–] length of the ocean shore between the opposite estuaries [L] =6Lp, intrusion distance [L] =L/Lp, normalized shoreline length [–] pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ 2T=xS, tidal propagation length [L] normalized sensitivity output function of the groundwater model ith input parameter of the groundwater model real part of complex function storativity of the aquifer [–] sensitivity of the ith input parameter Pi transmissivity of the aquifer [L2 T
1] separation variable of X in Eq. (A7) [–] Eq. (A13) [–] separation variable of Y in Eq. (A7) [–] Eq. (A14) [–] =x/L, dimensionless distance of x [–] =y/L, dimensionless distance of y [–]

(e.g., Ferris, 1951; Carr and van der Kamp, 1969; Fetter, 1994; Millham and Howes, 1995). Sun (1997) extended Jacob’s (1950) solution further into two-dimensional (2-D) groundwater flow and considered the head fluctuation within a confined aquifer adjacent to an estuary, where the tidal amplitude and phase vary along the estuarine bank. His solution, however, is 1-D, or called pseudo 2-D, with non-interacting tidal waves along the seashore because of not including the propagating effect of oceanic tides on the aquifer. Li et al. (2000) took this situation into account and developed a 2-D analytical solution based on Green’s function method for an aquifer with an L-shaped shoreline. Their solution demonstrated that the propagation of cross-shore oceanic waves may interfere with the estuarine tidal signals significantly in the coastal aquifer. Such an interaction behavior may be too complicated to be predicted by the solution of Sun (1997). To explore the influence of the coastline shape on tidal head fluctuation, Li et al. (2002) also developed a solution to the same problem as did in Li et al. (2000) but did not require the use of initial condition for the groundwater head distribution in the L-shaped coastal aquifer. At natural coastal areas, coastlines are often full of inlets, bays, headlands, etc. Some literatures indicated that the configuration of coastlines has a significant effect on the tidal head fluctuations in the coastal aquifer. However, it is not rare to see that the aquifer may be bounded by three water–land boundaries, which forms a U-shaped aquifer, where the two opposite sides are estuary–land boundaries and the other one is an ocean–land boundary. One

bm c cm dm h i k1 k1i k1r k2 k2i k2r ki kr m t u x y

DPi X

c cm

u k km

x

n

coefficient in Eq. (A13) [–] phase shift of tidal wave along ocean shore (radian) coefficient in Eq. (A14) [–] coefficient in Eq. (A14) [–] tidal head fluctuation [L] pffiffiffiffiffiffiffi ¼
1 =k1r + ik1i [L
1] wave number of tidal wave for the left estuary [L
1] damping coefficient of tidal amplitude for the left estuary [L
1] k2r + ik2i [L
1] wave number of tidal wave for the right estuary [L
1] damping coefficient of tidal amplitude for the right estuary [L
1] wave number of the tidal loading in the estuary [L
1] amplitude damping coefficient of the tidal loading in the estuary [L
1] index of the eigenvalue [–] time [T] qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ n2
X in Eq. (21) [–] along-shore distance from the left aquifer–estuary interface [L] cross-shore distance landward along the estuary from the ocean shore [L] =10
3Pi, increment of the ith parameter Pi ¼
ixL2 ðS=TÞ ¼
2iL2n in Eq. (13) [–] ¼
ðk þ XÞ in Eq. (A10) [–] ¼
ðkm þ XÞ in Eq. (A14) [–] phase shift of tidal head within the aquifer (radian) separation constant in Eq. (A8) [–] =
(mp)2, the mth eigenvalue in Eq. (A12) [–] angular frequency of the tidal wave [T
1] transform parameter of the Fourier sine transform [–]

example as shown in Fig. 1 has a nearly U-shaped aquifer system located at Linyuan District, Kaohsiung City. This site is a part of Pingtung Plain, which is an important groundwater source for aquaculture and coastal fishery in the southwestern Taiwan. It is important to note that the tidal propagation along the estuary exhibits the phenomena of amplitude attenuation and phase shift. Under those circumstances, the groundwater flow in the tidal aquifer system is naturally subject to the effects of those three water– land boundaries simultaneously. The characteristics of interaction behavior of tidal waves may be rather complicated and need to be investigated thoroughly. This study presents an analytical model considering three tidal boundary effects for U-shaped aquifers. The solution of the model is developed based on the methods of separation of variables and Fourier sine transform. The content of this study is organized as follows: a conceptual model and associated boundary conditions are first introduced, and the corresponding analytical solution is then derived. According to the solution, two special cases with different lengths of ocean shore are demonstrated for evaluating the head fluctuation within the aquifer. After that, the present solution is compared with the Jacob’s (1950) 1-D solution, Sun’s (1997) solution for pseudo 2-D flow, and Li et al.’s (2002) solution for 2-D flow. The influences of the normalized shoreline length (Ln) on the head fluctuation, and the variations of tidal head along the normalized diagonal line from the entry of estuary for different Ln are then examined. The former is to identify the transition distance of shoreline that judges whether the interaction among tides

F.-K. Huang et al. / Journal of Hydrology 530 (2015) 291–305

293

Fig. 1. The aquifer located at the Linyuan District, Kaohsiung City, within the Pingtung Plain, forms a U-shaped aquifer with shoreline length of about 1.4 km. The aquifer is bounded by three water–land boundaries (modified from Google Earth (2015)).

is significant, based on the properties of the aquifer and the tidal waves; whereas the latter is to investigate the influence of hydraulic properties upon the tidal head fluctuation within the aquifer. Finally, the sensitivity analysis is performed to assess the response of groundwater head fluctuation to the change in each of the hydraulic parameters and tidal parameters. The characteristics of spatial distribution of the sensitivity curves for tidal head fluctuation are then clearly explored.

2. Conceptual model Consider a horizontal non-leaky aquifer with a U-shaped coastline as displayed in Fig. 2, where three water–land boundaries are considered. The side contacted with the ocean, denoted as A1 A2 , is an ocean–land boundary, whereas the other two opposite sides, denoted as A1 A3 and A2 A4 , adjacent to estuaries represent estuary–land boundaries. The estuarine tidal waves propagating along

294

F.-K. Huang et al. / Journal of Hydrology 530 (2015) 291–305

hðL; y; tÞ ¼ Ae
k2r y cosðxt
k2i y þ cÞ ¼ A Re½e
k2 yþiðxtþcÞ 

estuarine tides

U-shaped aquifer

oceanic tides A1

ð4Þ

where k1r P 0 and k2r P 0 are the damping coefficients of tidal amplitude [L
1] for the left and right estuaries, respectively; k1i P 0 and k2i P 0 are the wave numbers of tidal wave [L
1] for the left and right estuaries, respectively; k1 = k1r + ik1i and k2 = k2r + ik2i; L is the length [L] of the ocean shore between the opposite estuaries.

A4

estuarine tides

estuary-land boundary

A3

estuary-land boundary

y or Y

On the inland side ðA3 A4 Þ far from the seashore, the no-flow boundary condition is assumed, i.e.,

lim

@h

ðx; y; tÞ ¼ 0

y!1 @y

ð5Þ

Physically, this condition implies that the tidal effect on the aquifer diminishes as y approaches infinity.

A2 x or X

3. Analytical solution

ocean-land boundary To solve Eq. (1) subject to the boundary conditions of Eqs. (2)–(5), two dimensionless variables are introduced as Fig. 2. Conceptual model of 2-D groundwater flow in a coastal aquifer bounded by U-shaped coastlines with four zones classified based on interaction behavior of the oceanic and estuarine tidal waves. X–Y is the dimensionless coordinate of x–y normalized with respect to L, which represents the length of ocean shore between the two opposite estuaries.

the ocean shore direction (i.e., parallel to A1 A2 ) will interact with the oceanic tidal waves transmitting in the cross-shore direction (i.e., parallel to A1 A3 or A2 A4 ) within the aquifer. Assume that the aquifer is homogeneous and isotropic with a uniform thickness, and directly connected to the tidal water. The water–land boundary is assumed vertical with neglecting the effects of beach slope and seepage face for simplicity. Under these conditions and Dupuit assumption (Bear, 1972), the governing equation describing the 2-D transient groundwater flow in the aquifer can be written as 2

2

@h T @ h @ h ¼ þ @t S @x2 @y2

!

as

ð2Þ

where A and x are amplitude [L] and angular frequency [T
1] of the tidal wave, respectively; c is the phase shift (radian) of the tidal pffiffiffiffiffiffiffi wave at the boundary; i ¼
1; and ‘Re’ represents the real part of the complex function. Note that only one wave component is considered here for simplicity, and the initial hydraulic head in the aquifer is assumed to be zero everywhere. Other tidal constituents can easily be included based on the principle of superposition because Eq. (1) is linear. The estuary–land boundaries, accounting for tidal amplitude damping and phase shift along the estuaries, are given as (Sun, 1997; Li et al., 2000, 2002)

ð3Þ

y L

ð6Þ

Hence, the governing equation and associated boundary conditions becomes:

@h T 1 @ 2 h @ 2 h þ ¼ @t S L2 @X 2 @Y 2

!

ð7Þ

hðX; 0; tÞ ¼ A Re½eiðxtþcÞ 

ð8Þ

hð0; Y; tÞ ¼ A Re½e
k1 LYþiðxtþcÞ 

ð9Þ

hð1; Y; tÞ ¼ A Re½e
k2 LYþiðxtþcÞ 

ð10Þ

lim

ð1Þ

Along the coastline ðA1 A2 Þ, the tidal boundary can be expressed

hð0; y; tÞ ¼ Ae
k1r y cosðxt
k1i y þ cÞ ¼ A Re½e
k1 yþiðxtþcÞ 

Y¼

@h

Y!1 @Y

where h(x, y, t) is the spatial and temporal distribution of groundwater head [L] for a confined aquifer, or the head of water table for an unconfined aquifer if the tidal amplitude is small compared to the aquifer thickness; x and y are the coordinates [L] (Fig. 2) with x being the distance along the ocean shore from the left aquifer– estuary interface, and y being the distance landward along the left estuary from the entry; t is the time [T]; and T and S are the transmissivity [L2 T
1] and storativity (dimensionless) of the aquifer, respectively.

hðx; 0; tÞ ¼ A cosðxt þ cÞ ¼ A Re½eiðxtþcÞ 

x X¼ ; L

ðX; Y; tÞ ¼ 0

ð11Þ

Now assume that

hðX; Y; tÞ ¼ A Re½HðX; YÞeiðxtþcÞ 

ð12Þ

where H(X, Y) is a complex function (dimensionless). And let

X ¼
ixL2

  S T

ð13Þ

Substituting Eq. (12) into Eqs. (7)–(11) and using Eq. (13) results in the following dimensionless equations:

@2H @X

2

þ

@2H @Y 2

þ XH ¼ 0

ð14Þ

HðX; 0Þ ¼ 1

ð15Þ

Hð0; YÞ ¼ e
k1 LY

ð16Þ

Hð1; YÞ ¼ e
k2 LY

ð17Þ

lim

Y!1

@HðX; YÞ ¼0 @Y

ð18Þ

Eqs. (15)–(17) are nonhomogeneous; the problem is then decomposed into two parts, each of which has homogeneous boundary conditions on parallel boundaries, as demonstrated in Fig. 3. The methods of separation of variables and Fourier sine transform are then employed to solve the boundary value problems according to the homogeneous state of boundaries with a finite domain in X direction and a semi-infinite domain in Y direction. Detail

295

F.-K. Huang et al. / Journal of Hydrology 530 (2015) 291–305

∂H1 =0 ∂Y

lim

Y →∞

Y

∂2H ∂2H + +ΩH =0 ∂X 2 ∂Y 2

ľ

B1'

(1)

B2'

∂ 2 H1 ∂ 2 H1 + + Ω H1 = 0 ∂X 2 ∂Y 2

Ĭ

B3'

B3

H =1

B1''

(2)

B2''

∂2H2 ∂2H2 + + Ω H2 = 0 ∂X 2 ∂Y 2 B3''

X 1

H 2 = e − k1LY

B2

H1 = 0

H = e − k2 LY

B1

∂H 2 =0 ∂Y B4''

B4'

B4

H = e − k1LY

lim

Y →∞

Y

H 2 = e − k2 LY

∂H =0 ∂Y

H1 = 0

lim

Y →∞

Y

X 1

H1 = 1

X

H2 = 0

1

Fig. 3. Decomposition of the problem into two parts, each of which has homogeneous boundary conditions on parallel boundaries.

derivations of the solution are presented in Appendix A, and the head fluctuation h(x, y, t) is given by (" 1 pffiffiffiffiffiffiffiffiffi y mpx X 2 hðx;y; tÞ ¼ ARe ð1
cosðmpÞÞe
km þXiðLÞ sin mp L m¼1     Z ux ux 2 1 ny dn eiðxtþcÞ Ecosh ð19Þ þ þ F sinh sin L L L p 0

amplitude and the wave numbers of tidal wave for the two opposite estuaries are assumed to be equal for simplicity. Now define a tidal propagation length [L]

where

where D (=T/S) is the hydraulic diffusivity of the aquifer [L2 T
1] and a normalized shoreline length Ln defined as L/Lp (dimensionless). For the present study, Lp is calculated to be 87.4 m based on the parameters given in Table 1. According to the solution of Jacob (1950) for 1-D groundwater flow in a coastal aquifer, the tidal waves decay inland exponentially and diminish over a certain distance. The distance inland at which the amplitude of head fluctuation reduced to 1% and 0.25% of that at the shore is about 4.6Lp (Li et al., 1999, 2000) and 6Lp, respectively. In this study, different lengths of ocean shore, L, in terms of various multiples of the tidal propagation length are used to study the interaction behavior of the tidal waves. To evaluate the tidal head variation, the spatial distributions of the amplitude and phase shift of the head fluctuation h(x, y, t) or h (X, Y, t) within the aquifer are calculated based on the solutions of Eqs. (19)–(24) or Eqs. (A36) and (A37). Note that the amplitude of qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi head is equal to ðReðHðX; YÞÞÞ2 þ ðImðHðX; YÞÞÞ2 (dimensionless)

km ¼
ðmpÞ2

ð20Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffi n2
X

ð21Þ

u¼ E¼

F¼

n

ð22Þ

K 2L þ n2 !

n K 2R þ n2

K L ¼ k1 L;

csch u


n K 2L þ n2

! coth u

K R ¼ k2 L

ð23Þ ð24Þ

4. Results and discussion 4.1. Amplitude and phase shift of the tidal head fluctuation Eq. (19) indicates that the head fluctuation in the aquifer can be divided into two parts. The first part is in terms of an infinite series, reflecting the influence of the cross-shore oceanic tide, whereas the second is an infinite integral, representing the effect of the alongshore estuarine tides. The interaction behavior of tidal waves driven from the oceanic tides and the estuarine tides within the aquifer is therefore rather complicated. To demonstrate the distribution of head fluctuation in the aquifer, a hypothetical example is considered and the parameter values listed in Table 1 are used. The damping coefficients of tidal

Table 1 Values of hydraulic and tidal parameters used in this study. Parameter

Value

Transmissivity, T Storativity, S Tidal amplitude along ocean shore, A Tidal frequency, x

100 m2/h 0.1 1m 0.2618 h
1 (i.e., a diurnal tide) 0 8  10
5 m
1

Phase shift of tidal wave along ocean shore, c Damping coefficient of tidal amplitude, k1r = k2r Wave number of tidal wave, k1i = k2i

6.67  10
7 m
1

Lp ¼

rffiffiffiffiffiffiffi rffiffiffiffiffiffiffi 2T 2D ¼ xS x

ð25Þ

and the phase shift of head fluctuation within the aquifer denoted as u is tan
1(Im(H(X, Y))/Re(H(X, Y))) (radian), in which ‘Re’ and ‘Im’ represent the real part and imaginary part of the complex function H(X, Y), respectively. Now consider that the U-shaped aquifer with L = 5Lp, or Ln = 5, the contour lines of the amplitude and phase shift for this case are shown in Fig. 4. As expected, the amplitude of the tidal head decays landward both in the crossand along-shore directions because of the effect of damping. On the other hand, the associated phase shift increases with the distance from ocean shore or estuarine shore when the waves propagate inland. Since the tidal waves in the estuary attenuate with the distance from the entry, the contour lines of amplitude tilt toward the estuarine bank as they extend inland from the seashore. But most importantly in this case, significant interactions between the estuarine tidal waves and the oceanic tidal waves near the two corners of the aquifer can be observed from the curved contour lines. Even though in the middle area around Y < 0.6 in the aquifer, the contour line for the amplitude of 0.2 also exhibits U shape indicating that the wave interferences still prevails over the central region. Not only the oceanic tides propagate into the inland aquifer, the estuarine waves from both left and right sides also transmit landward simultaneously. If the length of ocean shore, L, is not long enough, the waves inland from the three different directions will interfere with each other intensely. That is, the interaction effect among those waves on the head fluctuation in a tidal aquifer can be easily observed.

296

F.-K. Huang et al. / Journal of Hydrology 530 (2015) 291–305

(b) Phase shift of head fluctuation (L/Lp = 5)

1.4

1.4

1.2

1.2

1

1

Y (= y/L)

Y (= y/L)

(a) Amplitude of head fluctuation (L/Lp = 5)

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

0.8

0

1

Unit: degree

0

0.2

X (= x/L)

0.4

0.6

0.8

1

X (= x/L)

Fig. 4. Contour lines of (a) amplitude and (b) phase shift distributions of the head fluctuation for L/Lp = 5. Significant interactions between the estuarine tidal waves and the oceanic tidal waves near the U-shaped corner of the aquifer can be observed.

(b) Phase shift of head fluctuation (L/Lp=15)

1.4

1.4

1.2

1.2

1

1

Y (= y/L)

Y (= y/L)

(a) Amplitude of head fluctuation (L/Lp=15)

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0

0

0.2

0.4

0.6

X (= x/L)

0.8

1

0

Unit: degree

0

0.2

0.4

0.6

0.8

X (= x/L)

Fig. 5. Contour lines of (a) amplitude and (b) phase shift distributions of the tidal head fluctuation for L/Lp = 15.

1

F.-K. Huang et al. / Journal of Hydrology 530 (2015) 291–305

As the length of ocean shore is longer, for example L = 15Lp, the corresponding contour lines of amplitude and phase shift of the head fluctuation in the aquifer are demonstrated in Fig. 5. In contrast to the case of L = 5Lp, strong interactions of tidal waves still exist in the regions near those two U-shaped corners. The wave interaction is however weakened as the distance from the estuary or ocean increases. In areas far from the estuary, the head fluctuation of the aquifer is dominated by the oceanic tide while the head fluctuation in areas far from the ocean shore (i.e., large y or Y) is mainly dominated by the estuarine tide. In order to further study the tidal-wave behavior, the distance, 6Lp, at which the amplitude of head fluctuation reduced to 0.25% of that at the shore for 1-D groundwater flow is regarded as an intrusion distance [L], denoted as:

Li ¼ 6Lp

ð26Þ

Such a distance is similar to the critical distance defined in Li et al. (2000) for the head fluctuation reduces to 1%. Based on the contour lines plotted in Fig. 5 and the intrusion distance, the aquifer with different head fluctuations can be divided into four zones with qualitatively different interaction behavior. Two dividing lines in x direction, i.e., X (=x/L), in the dimensionless domain are located at 0.4 (X = Li/L = 6Lp/15Lp) and 0.6 because of symmetry of the aquifer when L = 15Lp. Similarly, the dividing line in y direction is located at Y = 0.4. According to the amplitude variations of the tidal head shown in Fig. 5(a), there are conspicuous interactions in Zone I where the oceanic and estuarine tidal waves interact strongly. The closer the curved contour lines, the stronger the interaction between waves near the regions of U-shaped corners. In Zone II, the head fluctuation is dominated by the oceanic tidal waves. Contour lines are parallel to the shoreline indicating that the damping effect for amplitude is uniform and the heads are unaffected by the estuarine tidal waves from both sides of estuarine. The response of aquifer to the tides in this zone is typical of 1-D mode and thus the head fluctuation can be predicted by the solution of Jacob (1950). In contrast to that, Zone III is predominantly influenced by the along-shore propagation of estuarine tidal waves. Along the estuary from the entry, the tidal loading varies and attenuates such that the aquifer’s response to tidal oscillation in this zone behaves 2-D. Hence the head fluctuation, with non-interacting wave, within the aquifer far away from the seashore can be predicted by the solution of Sun (1997). In Zone IV, the upper middle region of the aquifer, the tidal fluctuation is small and negligible. As to the characteristics of phase variation of the tidal head fluctuation, a similar behavior of aquifer response can also be observed in Fig. 5(b). In Zone I, significant interactions are evident between the oceanic and estuarine tidal waves, whose effects are remained and extended to Zone IV. These phenomena, however, become insignificant in Zones II and III except for small areas near the dividing lines. If the zoning criteria for the case of L = 15Lp is also applied to the case of L = 5Lp as shown in Fig. 4, one immediately finds that none of zones can be differentiated. Since the length of ocean shore is not large enough, the tidal interactions among the oceanic wave and two estuarine waves from both sides are significant over the whole aquifer. Under such a circumstance, only the present solution is applicable to predict the tidal head fluctuation in the coastal aquifer. However, the interactions become weakened when Y is greater than 1.2 (i.e., Y = Li/L = 6Lp/5Lp). 4.2. Comparison with the existing analytical solutions The predicted results from the present solution are compared with those predicted by the analytical solutions of Jacob (1950) for 1-D groundwater flow, Sun (1997) for 2-D flow due to tidal

297

loading in an estuary, and Li et al. (2000) and Li et al. (2002) for 2-D flow induced by both oceanic waves and estuarine tidal signals in an L-shaped aquifer. Table 2 summarizes those four solutions for groundwater head fluctuation in coastal aquifers. For the purpose of comparison, these four solutions are rewritten using the same notations as those in the present study. In the table, kr and ki represent the amplitude damping coefficient and the wave number of the tidal loading in the estuary, respectively, and K1() is the modified second-kind Bessel function of first order. In addition, the datum of the piezometric head in the aquifer is set to be the averaged mean sea level and assumed to be zero. 4.2.1. Comparison with the solution of Jacob (1950) The predicted result by the present solution for the case of L/Lp = 15 shown in Fig. 5 is used for demonstration. As illustrated in Zone II of Fig. 5, the aquifer’s response to tidal oscillation is dominated by the oceanic wave, which demonstrates the 1-D condition. Along the line of X = 0.5, i.e., the centerlines of Zones II and IV, the amplitude of the tidal head fluctuation within the aquifer is shown in Fig. 6. On the other hand, the result predicted from Jacob’s solution (1950) is also marked in the same figure. The figure shows close match between these two solutions in Zone II, indicating that the present solution can reduce to the 1-D solution of Jacob (1950). In addition, the amplitude of the tidal head almost diminishes when Y = 0.4 or y/Lp = 6 (bottom axis in Fig. 6). The area with a length larger than the intrusion distance belongs to Zone IV, where the tidal signals are very little and negligible. 4.2.2. Comparison with the solution of Sun (1997) In contrast to Zone II, the along-shore tidal waves from estuary dominate the response of aquifer in Zone III. Since the tide level varies along the estuary, the tidal head within the aquifer should be evaluated by a 2-D model. Notice that the 2-D solution of Sun (1997) is appropriate only for non-interacting tidal waves in the aquifer, the tidal amplitude at the location of intrusion distance, Y = 0.4 or y/Lp = 6, is therefore selected for comparison. The location of Y = 0.4 is also the margin to divide Zones I and III. The amplitudes of the tidal head, at Y = 0.4 along the ocean shore direction from left side up to X = 0.5, predicted by both the present solution and Sun’s solution (1997) are also shown in Fig. 6. This figure shows close match of the results in Zone III predicted by these two solutions, implying that the present solution can also reduce to Sun’s solution (1997). In addition, the figure indicates that the amplitude of the tidal head is almost vanished when at X = 0.4 or x/Lp = 6 (top axis in Fig. 6). The estuarine tidal signals at the location with a distance beyond the intrusion distance are also very little and negligible. 4.2.3. Comparison with the solution of Li et al. (2002) Apart from 1-D situation in Zone II and 2-D condition with noninteracting waves in Zone III, the groundwater fluctuation in Zone I is complicated due to the interacting tidal waves from both ocean and estuary. The solutions developed by Li et al. (2000) and Li et al. (2002) for an L-shaped aquifer can handle this scenario. The solution of Li et al. (2002) does not require an initial condition in deriving the solution; therefore, it is adopted to make comparison with the present solution. Now a reference line, mainly in Zone I, from the origin, with an angle of 45° with respect to the x-axis, is selected for calculating the head. The distance from the origin along the line can be reprepffiffiffi sented by the normalized diagonal distance, 2X. The curves of the head fluctuation amplitude, at locations along this line, predicted by the present solution and Li et al.’s solution (2002) are drawn in Fig. 7. The cases for shoreline lengths of L/Lp = 5 and 15 are considered for comparison. Note that four zones in the aquifer can be

298

F.-K. Huang et al. / Journal of Hydrology 530 (2015) 291–305

Table 2 Some existing analytical solutions for groundwater head in coastal aquifer. Source

Assumption

Groundwater head, h pffiffiffi ffi qffiffiffiffiffi   xS S hðy; tÞ ¼ Ae
2T y cos xt
x 2T y þ c

Jacob (1950)

1-D, straight coastline

Sun (1997)

Pseudo 2-D, straight coastline in an estuary

  r ki T hðx; y; tÞ ¼ Ae
px
kr y cos xt
ki y
xS
2k xþc 2pT where p ¼ p1ffiffi2

Li et al. (2000)

( 

2

2

ki
kr

2



2 þ xTS
2kr ki

)1=2

1=2

2

2

þ ki
kr

pffiffiffi ffi qffiffiffiffiffi   nR o xS t S hðx; y; tÞ ¼ Ae
2T y cos xt
x 0 ½f ðka ; yÞ
f ðka ;
yÞ
f ðkb ; yÞ þ f ðkb ;
yÞdt 0 2T y þ c þ A Re

2-D, L-shaped coastlines

where

     x2 0 Þþf pffiffiffiffiffiffiffiffiffiffiffiffi ;  exp n2 Dðt
t0 Þ þ ixt 0 þ nf
4Dðt
t  1 þ erf 2nDðt
t 0Þ 2 Dðt
t 0 Þ qffiffiffiffiffi qffiffiffiffiffi T x S x S D ¼ S ; ka ¼
ðkr þ iki Þ; kb ¼
2T þ i 2T f ðn; fÞ ¼

Li et al. (2002)

pffiffiffi x 4 p½Dðt
t 0 Þ3=2



hðx; y; tÞ ¼ A Re ½Iðax; ay; 1 þ i; 1 þ iÞ þ Iðay; ax; m þ in; 1 þ iÞ þ e
ð1þiÞay þ e
ke y
ðmþinÞax eiðxtþcÞ where h i 2 R1 s;kÞÞ K 1 ðqðn;gþs;kÞÞ
qðn;gþs;kÞ ds; Iðn; g; l; kÞ ¼
kpn 0 e
ls K 1qðqðn;ðn;gg

s;kÞ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qðn; g; kÞ ¼ k n2 þ g2 ; qffiffiffiffiffi S a¼ x 2T ;

2-D, L-shaped coastlines

ke ¼ kr þ iki "rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi #  2  2 2 2  2 2 2 1=2 k
k k
k m¼ 1
kar 2ki þ r2a2 i ;
r2a2 i n¼ This study

2-D, U-shaped coastlines

"rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi #  2  2 2 2  2 2 2 1=2 k
k k
k 1
kar 2ki þ r2a2 i þ r2a2 i

Eqs. (19)–(24)

U-shaped aquifer can only be evaluated by the present solution. Therefore, it is important to investigate how large the shoreline length is when the interaction behavior among tides becomes significant.

Normalized shoreline distance along seashore direction, x/L p (at Y=0.4) 0

1

2

3

4

5

6

7

Amplitude of head fluctuation

1.0 1 D (X=0.5; Jacob, 1950) 2 D (X=0.5; this study, L/L p=15) 2 D (Y=0.4; Sun, 1997) 2 D (Y=0.4; this study, L/L p=15)

0.8

4.3. Effect of shoreline length on the tidal interactions 4.3.1. Determination of transition distance for the shoreline Now consider a U-shaped aquifer with several different lengths of shoreline in terms of the normalized shoreline length Ln. The curves of amplitude and phase shift of the tidal head fluctuation versus Ln are calculated for various y/Lp distances at X = 0.5 demonstrated in Fig. 8. Fig. 8(a) shows that the amplitude decreases with Ln sharply for small Ln and gently for large Ln until Ln > 10, where the amplitude becomes unchanged with the length of ocean shore. In other words, there requires about 10Lp for the length of ocean

0.6

0.4

0.2

0.0 0

5

10

15

Normalized shoreline distance along estuary direction, y/L p (at X=0.5)

depicted clearly as exhibited in Fig. 5 if L = 15Lp; in contrast, the tidal wave interaction is too strong to distinguish from different zones as indicated in Fig. 4 if L = 5Lp. Fig. 7 shows that the predicted results from these two solutions yield close match for the case of L/Lp = 15, indicating that the solution of Li et al. (2002) can also be considered as our special case when the length of ocean shore L is large. However, Fig. 7 shows the discrepancies between the two solutions for the case of L/Lp = 5 due to the joint effect of oceanic and estuarine tidal waves. For the case at the location of X = 0.5 and Y = 0.5, the solution of Li et al. (2002) under-estimate the tidal amplitude up to 32.7%. In such situation, the tidal head in the

L/Lp

Amplitude of head fluctuation

Fig. 6. Spatial distribution curves of the head amplitude predicted by the present solution and Jacob’s (1950) solution at X = 0.5, and by the present solution and Sun’s (1997) solution at Y = 0.4 for L/Lp = 15.

1.0 5 ( this study ) 15 ( this study )

0.8

5 ( L -shaped, Li et al., 2002 ) 15 ( L -shaped, Li et al., 2002 ) Y

0.6 Y = 0.5

0.4 X = 0.5

1

X

0.2

0.0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Normalized diagonal distance Fig. 7. Curves of the head amplitude versus normalized diagonal distance predicted by the present solution and Li et al.’s (2002) solution for L/Lp = 5 and 15.

299

F.-K. Huang et al. / Journal of Hydrology 530 (2015) 291–305

(a) 1.0 Amplitude

0.8

0.6 X = 0.5

1

X

y/L p 1 2 3 5 10 15 20

0.4

0.2

0.0 0

5

10

15

20

Normalized shoreline length, L n (=L/L p) 630

Y

Phase shift y/L p 1 2 3 5 10 15 20

540 450 360

X = 0.5

1

4.4. Effects of hydraulic properties on the head fluctuation

X

Note that the normalized shoreline length Ln can also be considered as a dimensionless length reflects the influence of hydraulic properties on the tidal fluctuation. The response of groundwater flow near the U-shaped corner, mainly within Zone I where the tidal interactions are most intensely, is thus evaluated to explore

270 180 90 0 0

5

10

15

20

Normalized shoreline length, L n (=L/L p) Fig. 8. The curves of (a) amplitude and (b) phase shift versus normalized shoreline length at X = 0.5 for different y/Lp.

shore of the U-shaped aquifer to avoid the influence of tidal interaction when y/Lp 6 3. From this figure one can also observe that when the location is closer to the ocean shore with a low y/Lp, where oceanic tides dominate, the amplitude curve tends to be stabilized earlier with a higher magnitude for a lower Ln. On the other hand, if the locations are at a distance of y/Lp > 5, there are no significant difference between the amplitude curves except that Ln 6 3. These amplitude curves won’t be stabilized until Ln > 12–15 at which the amplitude is almost vanished. Fig. 8(b) displays that the phase shift increases with Ln and tends to be stabilized quickly, especially for low y/Lp, where the oceanic tides dominate. The length of ocean shore is still about Ln = 10 to keep away from the effect of tidal interaction for y/Lp 6 3. If the location is far away from the ocean shore with y/Lp > 5, the phase shift is mainly dominated by the estuarine waves. Under this circumstance, the phase shift is increased with Ln and has no apparent difference between various distances of y/Lp, especially when y/Lp P 15. Based on the results demonstrated above, the shoreline length of Ln = 10–15 can be regarded as a transition length that one can judge whether the interaction among tides is significant, depending on the distance of y/Lp from the ocean shore considered. When referring to Eq. (25), one can easily observe that if the hydraulic diffusivity D of the aquifer is large, the transition length of the ocean shore and the distance of y from the shore to get rid of the influence of tidal interaction are large, too. In addition, the transition length is longer for diurnal tides than for semidiurnal tides because the former has a smaller oscillatory frequency.

(a) Amplitude of head fluctuation

Phase shift of head fluctuation, ϕ(degree)

(b)

4.3.2. The case of Linyuan District, Kaohsiung City The aquifer system, mentioned in the Introduction section, located at the Linyuan District exhibits near a U shape with shoreline length of about 1.4 km. Due to heavy exploitation of groundwater for the fishery and industrial consumption, problems of seawater intrusion and land subsidence occurred within past few decades in this coastal region (Kuo and Wang, 2001; Hsu et al., 2007). Groundwater resources planning and management is therefore imperative there. The evaluation of tide-induced groundwater head fluctuations in the aquifer is thus of practical interest. The Pingtung Plain is underlain by unconsolidated sediments of Quaternary age, which constitute the main aquifer in the plain (Ting, 1997). Due to spatial variations, the transmissivity and specific yield of the aquifer are estimated to be 144 m2/h and 0.071, respectively, in the Linyuan site. The transition length thus determined corresponding to Ln = 10–15 is in the range of 1.25– 1.87 km when diurnal tide is considered. Accordingly, the tidal head fluctuations in the aquifer at Linyuan site subject to the effects of three water–land boundaries simultaneously can be evaluated by the present solution (i.e., Eqs. (19)–(24)).

1.0 Amplitude L n 1 2 3 5 10 15 20

0.8

0.6

0.4 Y

0.2

Y = 0.5

X = 0.5

0.0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

1

0.9

X

1

Normalized diagonal distance

(b) Phase shift of head fluctuation, ϕ (degree)

Amplitude of head fluctuation

Y

630 Phase shift L n 1 2 3 5 10 15 20

540 450 360 270 Y

180 Y = 0.5

90 X = 0.5

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

X

1

Normalized diagonal distance Fig. 9. The curves of (a) amplitude and (b) phase shift of the tidal head fluctuation versus normalized diagonal distance for various Ln.

300

F.-K. Huang et al. / Journal of Hydrology 530 (2015) 291–305

the influence of aquifer hydraulic properties. The curves of amplitude and phase shift of the tidal head fluctuation versus normalized diagonal distance are drawn in Fig. 9 from X = 0 to 0.5 for different Ln. Since Ln is equal to L/Lp, there are two ways to interpret its effect upon head for different values of Ln. One is to fix Lp, the influence of the length of ocean shore L on the groundwater head fluctuation can be observed. The other is to keep the length of ocean shore L as a constant, the effect of the hydraulic properties upon the response of aquifer can then be examined. In the first case of fixed Lp, Fig. 9(a) indicates that the head amplitude decreases markedly as the ocean shore L increases. The amplitude is just attenuated slightly for Ln = 1, but it is nearly reduced to zero for Ln = 20 while normalized diagonal distance ranges from 0 to 0.707. Obviously, the effect of the length of ocean shore on the head amplitude is significant. In the second situation of given L, the tidal amplitude drops dramatically for large Ln (i.e., small Lp) along the diagonal direction, implying that the damping effect is significant when the D is rather low. In other words, the tidal wave signals cannot propagate far away if the aquifer has a low T and/or a high S. Fig. 9(b) shows the plot for the phase shift of head fluctuation versus normalized diagonal distance. For the first case, the phase shift increases with normalized diagonal distance. The figure indicates that the longer the L is, the larger the phase shift will be. Whereas in the second situation, the smaller the Lp is, the larger the phase shift at the same location along the diagonal direction will be. That is, the phase lags more if the tidal amplitude decays larger. Substituting Lp in Eq. (25) into the parameter X in Eq. (13) yields X =
2iLn2. Because the head fluctuation h(x, y, t) in Eq. (19) is influenced by X, the effect of Ln in Fig. 9 on the amplitude and phase shift of the head can be examined explicitly. That is, the normalized shoreline length Ln which reflects the influence of hydraulic properties and the length of ocean shore can be used to evaluate the response of tidal head fluctuations. Following, the values of hydraulic and tidal parameters listed in Table 1 are used to explore the behavior of head fluctuations further. Based on Eq. (25), the tidal propagation length Lp is equal to 87.4 m. Two normalized shoreline length Ln, 5 and 15, are considered and thus the shoreline length are L = 5Lp = 437 m and L = 15Lp = 1311 m, respectively. To appreciate the influence of the magnitude of hydraulic parameters on the head fluctuation, other two transmissivity values (say, 500 and 1000 m2/h) and storativity values (0.01 and 0.001) are considered. Fig. 10 demonstrates the amplitude and phase shift of the tidal head fluctuation versus normalized diagonal distance with different parameter combinations for the cases of L = 5Lp and L = 15Lp, in which (A.) stands for amplitude denoted as solid lines and (P.) for phase shift with dash line. In the following, two ordinates are used for ease of comparison, the left ordinate corresponds to the amplitude while the right ordinate is for the phase shift. Fig. 10(a) shows that the amplitude attenuates and phase shift increases dramatically along the diagonal distance for small T and large S. Moreover, such a trend is more profound for large shoreline lengths for the case L = 15Lp shown in Fig. 10(b). It is interesting to note that the spatial distribution curves of the amplitude and phase shift for an aquifer with T = 1000 m2/h and S = 0.1 completely match with those with T = 100 m2/h and S = 0.01 as demonstrated in Fig. 10 because of the same hydraulic diffusivity value, i.e., D = 10,000 m2/h. 4.5. Sensitivity analysis 4.5.1. Definition of normalized sensitivity To further explore the effects of uncertainties in hydraulic parameters and tidal parameters on the head fluctuation, sensitivity analyses are performed based on the present solution. The

sensitivity that measures the effect of parametric variations on the model output can be expressed as (Yeh et al., 2008)

Si ¼

@O OðP i þ DPi ; Pjjj–i Þ
OðP 1 ; P 2 ; . . . ; Pn Þ ¼ @P i DP i

ð27Þ

where Si is the sensitivity of the ith input parameter Pi under consideration (i.e., T, S, k1r, k1i, k2r, or k2i); O is the output function of the model (i.e., the amplitude and phase shift of H(X, Y) shown in Eq. (A37)). To avoid the differences of unit and/or the magnitude of the parameters, the normalized sensitivity (NS) associated with different parameters is used and defined as (Liou and Yeh, 1997)

NSi ¼

@O @O ¼ Pi @Pi =Pi @Pi

ð28Þ

in which NSi is the NS of the ith input parameter Pi. The partial derivative in Eq. (28) is approximated by the forward finitedifference method as

NSi ¼ Pi

OðPi þ DPi Þ
OðPi Þ DP i

ð29Þ

where the increment DPi is approximated by DPi = 10
3Pi (Yeh et al., 2008). The NS will be employed to investigate the parameter sensitivity here. 4.5.2. Normalized sensitivity for amplitude and phase shift of head fluctuation The NS curves for the amplitude (left ordinate) and phase shift (right ordinate) of H(X, Y) versus normalized diagonal distance from 0 to 1.414 (i.e., X from 0 to 1) for various hydraulic parameters (T and S) and tidal parameters (kr and ki) are shown in Fig. 11(a) for L = 5Lp and in Fig. 11(b) for L = 15Lp. Note that the tidal parameters for the left and right estuaries are assumed equal, thus the subscripts 1 and 2 of k in the legend of the figure are dropped. The amplitude damping coefficients are abbreviated as kr and the wave numbers of tidal loading as ki. Fig. 11 indicates that the spatial distribution of each NS reflects the spatial change of the head fluctuation in response to the relative variation of each parameter. Among the four parameters, one can find that the amplitude and phase shift of the tidal head fluctuation are most sensitive to the changes in T and S and least sensitive to kr and ki in general. Interestingly, the sensitivity curves of T and S are symmetric in shape along the horizontal axis with same magnitude both for amplitude and phase shift, implying that the parameters T and S are both highly correlated. On the hydraulic parameters, Fig. 11(a) shows that the NS curve for the amplitude in response to the change of T has a highest value at the normalized diagonal distance near 0.56. A positive value of the NS indicates that the T has positive influence on the amplitude. Moreover, the larger the value of the NS, the higher the impact on the amplitude will be. As to the phase shift, its sensitivity curve associated with the change of T has a lowest value near the middle of the normalized diagonal distance (i.e., 0.707). Such a curve indicates that a positive change in T poses negative influence on the phase shift and the smallest value of the sensitivity has the highest negative influence. Interestingly, Fig. 11 indicates that T has a positive effect on the amplitude, but a negative effect on the phase shift. In contrast, S exerts a negative effect on the amplitude, but a positive effect on the phase shift. In view of Eq. (13), the parameter X in the solution of Eqs. (19) and (A37) is proportional to S, and inversely proportional to T. Obviously, the effect of T on the head H(X, Y) is opposite to that of S, the trend of NS for T is therefore reverse to that for S with the same shape. On the tidal parameters, the NS of amplitude for kr keeps zero initially and starts to turn to more negative value with increasing distance at the normalized diagonal distance of 0.8 for L = 5Lp shown in Fig. 11(a) and 1.0 for L = 15Lp in Fig. 11(b). Such a curve

301

F.-K. Huang et al. / Journal of Hydrology 530 (2015) 291–305

(a) 1.0 Amplitude of head fluctuation

T= 100, S=0.1 (A.) T= 100, S=0.1 (P.) T= 500, S=0.1 (A.) T= 500, S=0.1 (P.) T=1000, S=0.1 (A.) T=1000, S=0.1 (P.) S=0.01, T=100 (A.) S=0.01, T=100 (P.) S=0.001,T=100 (A.) S=0.001,T=100 (P.)

0.8

0.6

0.4

160

120

80

0.2

40

0.0

Phase shift of head fluctuation, ϕ (degree)

200 L = 5Lp (T: m2/hr)

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

Normalized diagonal distance

(b) 1.0 Amplitude of head fluctuation

T= 100, S=0.1 (A.) T= 100, S=0.1 (P.) T= 500, S=0.1 (A.) T= 500, S=0.1 (P.) T=1000, S=0.1 (A.) T=1000, S=0.1 (P.) S=0.01, T=100 (A.) S=0.01, T=100 (P.) S=0.001,T=100 (A.) S=0.001,T=100 (P.)

0.8

0.6

0.4

480

360

240

0.2

120

0.0

Phase shift of head fluctuation, ϕ (degree)

600 L = 15Lp (T: m2/hr)

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.1

1.2

Normalized diagonal distance Fig. 10. The amplitude and phase shift of the tidal head fluctuation versus normalized diagonal distance for various parameter combinations for the cases (a) L = 5Lp and (b) L = 15Lp.

indicates that a positive change in kr exerts negative influence on the amplitude. In addition, Fig. 11 shows that the NS of the phase shift for ki is too small to read its value. The effect of ki on the phase shift is therefore demonstrated in Fig. 12. The figure indicates that the sensitivity of the phase shift for ki increases with the normalized diagonal distance for both L = 5Lp and L = 15Lp except there is a little disturbance between 0.7 and 0.9 in the case of L = 15Lp. Such a disturbance is due to the effect of dynamic interaction among tidal waves as demonstrated in Fig. 5(b) where the target area is within or near the ‘nose’ of the contour lines. The phenomenon of different influence extent for kr and ki can be explained by Eqs. (3) and (4). It exhibits that the damping coefficient of tidal amplitude, kr, influences the amplitude only and no effect on phase shift. Inversely, the wave number of tidal wave, ki, influence the phase shift only and no effect on amplitude. As shown in Figs. 11 and 12, both NS of amplitude for kr and NS of phase shift for ki have the largest magnitude near the other estuary and far away the seashore. It is apparently that the shoreline length has a large impact on the normalized sensitivity. Taking the hydraulic parameters for example, the NS of head fluctuation for different parameters changes with shoreline length. For a shorter length like L = 5Lp shown in Fig. 11(a), the amplitude of head does not attenuate too low such that the magnitude of NS remains high in the middle

of the aquifer owing to the interactions of tidal waves. On the other hand, for L = 15Lp in Fig. 11(b), the magnitude of NS drops dramatically to near zero in the middle from its peak value at the normalized diagonal distance of 0.13, and then increases gradually to another peak at the normalized diagonal distance of 1.32 and descends again to zero near the other estuary. That is, for a larger shoreline length like L = 15Lp, the magnitude of NS of amplitude is high only at the areas near the estuary or seashore where the amplitude of head is large. The phenomena can be closely related to the contour lines distribution of amplitude of head shown in Figs. 4(a) and 5(a). When the NS of phase shift is concerned, one can observe from Fig. 11(a) that the sensitivity curve of the phase shift for L = 5Lp turns its direction at the middle of the diagonal distance because of the interaction behavior of tidal waves. As to L = 15Lp in Fig. 11(b), the curve exhibits slightly irregular for the normalized diagonal distance ranging from 0.7 to 0.9. This can also be attributed to the problem of the dynamic interaction among tidal waves as mentioned above. 5. Concluding remarks This paper derives a 2-D analytical solution to delineate the tidal groundwater fluctuation within a U-shaped aquifer bounded by three water–land boundaries. The interaction behaviors among

302

F.-K. Huang et al. / Journal of Hydrology 530 (2015) 291–305

(a) 0.4

80

T(A.) T(P.) S(A.) S(P.) k r (A.)

0.2

40

k i (P.)

0.0

0

-0.2

-40

-0.4

NS of phase shift of head (degree)

NS of amplitude of head

L = 5Lp (S=0.1, T=100m2/hr)

-80 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Normalized diagonal distance

(b) 0.4

240

T(A.) T(P.) S(A.) S(P.) k r (A.)

0.2

120

k i (P.)

0.0

0

-0.2

-120

-0.4

NS of phase shift of head (degree)

NS of amplitude of head

L = 15Lp (S=0.1, T=100m2/hr)

-240 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Normalized diagonal distance Fig. 11. The normalized sensitivities of amplitude and phase shift of the tidal head fluctuation versus normalized diagonal distance for the cases of (a) L = 5Lp and (b) L = 15Lp.

NS of phase shift of head (degree)

0.15 NS for ki (ki =6.67 ×10 m ) (S=0.1, T=100m2/hr) L = 5 Lp (P.) -7

0.10

-1

L =15 Lp (P.)

0.05

0.00

-0.05 0

0.2

0.4

0.6

0.8

1

1.2

1.4

1.6

1.8

2

Normalized diagonal distance Fig. 12. The NS of phase shift of the tidal head fluctuation to the wave number of tidal loading, ki, versus normalized diagonal distance for the cases of L = 5Lp and L = 15Lp.

the oceanic tides inland and the transmission of the two opposite estuarine tides are clearly explored. Four zones, I–IV, are identified to differentiate the response of the aquifer. It is demonstrated that

the present solution can reduce to Jacob’s (1950) solution for 1-D groundwater flow in Zone II, and Sun’s (1997) solution for 2-D situation with non-interacting wave in Zone III. Moreover, the solution of Li et al. (2002) for L-shaped coastlines, portraying the 2-D interaction of tidal waves near a cross-shore estuary in Zone I, is our special case if the length of ocean shore is large enough. Based on the present study, several important features and findings can be drawn below: Firstly, the present solution was employed to predict the tidal head distribution in a U-shaped costal aquifer with a finite shoreline. It has been found that 10–15 times of tidal propagation length Lp for the shoreline can be regarded as a transition distance to judge whether the interaction among tides over the cross-shore distance is significant. The transition length in the ocean shore and the intrusion distance from the shore will be large if the hydraulic diffusivity is large. Moreover, diurnal tides has a longer transition length than semidiurnal tides. Secondly, the influences of hydraulic properties and shoreline length, upon the joint effect of those three water–land boundaries for the aquifer water level fluctuation, were evaluated for various values of normalized shoreline length Ln quantitatively. The tidal amplitude drops dramatically for large Ln along the diagonal direction, demonstrating that the damping effect is significant for low transmissivity and/or high storativity. In addition, an aquifer with

303

F.-K. Huang et al. / Journal of Hydrology 530 (2015) 291–305

a higher Ln has a larger phase shift at the same location for the head fluctuation. In other words, a highly attenuated tidal amplitude has a large phase shift in the head fluctuation. Thirdly, sensitivity analysis was performed to investigate the tidal head fluctuation in response to the changes in parameters T, S, kr and ki. The spatial distribution of normalized sensitivities for amplitude and phase shift of head fluctuation indicates that the tidal head is most sensitive to T and S with opposite influence, and least to kr and ki in general. The magnitude of the sensitivities for T and S depends on the shoreline length and spatial location in the aquifer, and the interaction behavior among tidal waves. Both the sensitivities of amplitude for kr and phase shift for ki increase with diagonal distance and have the largest magnitude near the estuary and far away from the seashore.

where c ¼
ðk þ XÞ. The boundary conditions for U(X) can be obtained by substituting Eq. (A7) into Eqs. (A3) and (A4) as

Uð0Þ ¼ 0;

Uð1Þ ¼ 0

For U(X) of Eq. (A9) together with Eq. (A11), it constitutes a regular Sturm–Liouville problem. The eigenvalue and corresponding eigenfunction, i.e., the solution of U(X), are, respectively,

km ¼
ðmpÞ2 ;

This study was partly supported by the Taiwan Ministry of Science and Technology under the grants NSC 102-2221-E-009-072-MY2, MOST 103-2221-E-009-156, and MOST 104-2221-E-009-148-MY2. The authors would like to thank the associate editor and three anonymous reviewers for their valuable and constructive comments.

U m ðXÞ ¼ bm sinðmpXÞ;

ðA12Þ

m ¼ 1; 2; . . .

ðA13Þ

where bm is constant. As to V(Y) of Eq. (A10), the general solution, with k (and thus c) determined, is pffiffiffiffi cm Y

þ dm e


pffiffiffiffi

cm Y

ðA14Þ

where cm ¼
ðkm þ XÞ, cm and dm are constants. Then the solution of H1(X, Y), based on Eq. (A7) and superposition principle, is

H1 ðX; YÞ ¼

1 X U m ðXÞV m ðYÞ m¼1

¼

1 X 

Am e

pffiffiffiffi

cm Y

pffiffiffiffi

þ Bm e


cm Y



sinðmpXÞ

ðA15Þ

m¼1

Appendix A. Derivation of the solution Assume that the complex function H(X, Y) can be split into two parts:

HðX; YÞ ¼ H1 ðX; YÞ þ H2 ðX; YÞ

m ¼ 1; 2; . . .

V m ðYÞ ¼ cm e

Acknowledgements

ðA11Þ

ðA1Þ

The governing equation of Eq. (14) and boundary conditions of Eqs. (15)–(18) can also be divided into two parts such that each with homogeneous boundary conditions at parallel boundaries as exhibited in Fig. 3. The final solution is the sum of the solution obtained from those two individual parts. The derivations of the solution for H1(X, Y) and H2(X, Y) are as follows: (1) Derivation of the solution for H1(X, Y)

where Am = bm cm and Bm = bm dm. H1(X, Y) satisfies Eq. (A6), which yields Am = 0. Then Eq. (A5) gives

H1 ðX; 0Þ ¼

1 X Bm sinðmpXÞ ¼ 1

ðA16Þ

m¼1

The above infinite series is a Fourier sine series, thus the coefficients Bm for m = 1, 2, . . . can be given by

Bm ¼

2 1

Z

1

1  sinðmpXÞdX ¼

0

2 ð1
cosðmpÞÞ mp

ðA17Þ

From which, the solution of H1(X, Y) results in

H1 ðX; YÞ ¼

1 X pffiffiffiffi Bm e
cm Y sinðmpXÞ m¼1

The governing equation for H1(X, Y) is 2

@ H1 @X 2

2

þ

@ H1 @Y 2

þ XH1 ¼ 0

¼ ðA2Þ

1 pffiffiffiffiffiffiffiffiffi X 2 ð1
cosðmpÞÞe
km þXiY sinðmpXÞ mp m¼1

ðA18Þ

(2) Derivation of the solution for H2(X, Y)

with following boundary conditions

H1 ð0; YÞ ¼ 0

ðA3Þ

H1 ð1; YÞ ¼ 0

ðA4Þ

H1 ðX; 0Þ ¼ 1

ðA5Þ

@H1 ðX; YÞ ¼0 Y!1 @Y

ðA6Þ

lim

ðA7Þ

where U is a function of X and V is a function of Y. Substituting Eq. (A7) into Eq. (A2) and letting k as the separation constant obtains 00

00

U V ¼

X¼k U V

ðA8Þ

Accordingly, Eq. (A8) is equivalent to two linear ordinary differential equations (ODEs) expressed as

U 00
kU ¼ 0

ðA9Þ

V
cV ¼ 0

@X

2

þ

@ 2 H2 @Y 2

þ XH2 ¼ 0

ðA19Þ

H2 ð0; YÞ ¼ e
k1 LY

ðA20Þ

H2 ð1; YÞ ¼ e
k2 LY

ðA21Þ

H2 ðX; 0Þ ¼ 0

ðA22Þ

lim

Y!1

@H2 ðX; YÞ ¼0 @Y

ðA23Þ

On the interval [0, 1) in Y domain, when at Y = 0 and Y ? 1, the boundary conditions are homogeneous so that the Fourier sine transform, with respect to Y, is adopted to solve for H2(X, Y). We define

Z

F s ½H2 ðX; YÞ ¼

and 00

@ H2

with following associated boundary conditions

In view of the boundary conditions, H1(X, Y) can be determined by the method of separation of variables. Let

H1 ðX; YÞ ¼ UðXÞVðYÞ

The governing equation for H2(X, Y) is 2

1

H2 ðX; YÞ sinðnYÞdY ¼ H2 ðX; nÞ

0

ðA10Þ

where n is a transform parameter.

ðA24Þ

304

F.-K. Huang et al. / Journal of Hydrology 530 (2015) 291–305

Taking the Fourier sine transform of Eq. (A19) results in

Fs

" # @ 2 H2 @X 2

þ Fs

" # @ 2 H2 @Y 2

þ F s ½XH2  ¼ F s ½0

ðA25Þ

Suppose that H2(X, Y) satisfies the operational properties of the Fourier sine transform, it yields 2

d H2 ðX; nÞ dX

2


n2 H2 ðX; nÞ þ nH2 ðX; 0Þ þ XH2 ðX; nÞ ¼ 0

ðA26Þ

  hðX;Y;tÞ ¼ ARe HðX;YÞeiðxtþcÞ (" 1 pffiffiffiffiffiffiffiffiffi X 2 ð1
cosðmpÞÞe
km þXiY sinðmpXÞ ¼ ARe mp m¼1   Z 2 1 ðEcoshðuXÞ þ F sinhðuXÞÞ sinðnYÞdn eiðxtþcÞ þ ðA37Þ

p

0

Based on Eq. (6) and replacing X by x/L and Y by y/L, the solution h(x, y, t) reads as Eq. (19).

Substituting Eq. (A22) into (A26) gives 2

d H2 ðX; nÞ dX

2


ðn2
XÞH2 ðX; nÞ ¼ 0

ðA27Þ

Let u2 = n2
X, Eq. (A27) becomes 2

d H2 ðX; nÞ dX

2


u2 H2 ðX; nÞ ¼ 0

ðA28Þ

Since the domain of variable X is finite, the solution of the second order ODE can be expressed, with two undetermined coefficients E and F, as

H2 ðX; nÞ ¼ E coshðuXÞ þ F sinhðuXÞ

ðA29Þ

The transforms of Eqs. (A20) and (A21) are, respectively,

H2 ð0; nÞ ¼ F s ½H2 ð0; YÞ ¼ F s ½e
k1 LY  ¼

Z

1

e
k1 LY sinðnYÞdY

0

¼

n 2

ðk1 LÞ þ n2

¼

n

ðA30Þ

K 2L þ n2

H2 ð1; nÞ ¼ F s ½H2 ð1; YÞ ¼ F s ½e
k2 LY  ¼

Z

1

e
k2 LY sinðnYÞdY

0

¼

n 2

ðk2 LÞ þ n2

¼

n

ðA31Þ

K 2R þ n2

where KL = k1L and KR = k2L. Applying these conditions to the solution H2 ðX; nÞ in Eq. (A29) leads to

H2 ð0; nÞ ¼ E ¼

n

ðA32Þ

K 2L þ n2

H2 ð1; nÞ ¼ E cosh u þ F sinh u ¼

n K 2R þ n2

ðA33Þ

Substituting E in Eq. (A32) into Eq. (A33), the coefficient F is, in turn, determined as

F¼

!

n K 2R þ n2

cschu


n K 2L þ n2

!

coth u

ðA34Þ

Therefore, H2 ðX; nÞ in Eq. (A29) is obtained. Based on the inverse Fourier sine transform, one has

Z   2 1 H2 ðX; YÞ ¼ F
1 H ðX; nÞ ¼ H ðX; nÞ sinðnYÞdn 2 s p 0 2 Z 1 2 ðE coshðuXÞ þ F sinhðuXÞÞ sinðnYÞdn ¼

p

ðA35Þ

0

Finally, the solution of H(X, Y) in Eq. (A1) is the sum of H1(X, Y) in Eq. (A18) and H2(X, Y) in Eq. (A35), that is

HðX; YÞ ¼ H1 ðX; YÞ þ H2 ðX; YÞ 1 pffiffiffiffiffiffiffiffiffi X 2 ð1
cosðmpÞÞe
km þXiY sinðmpXÞ ¼ m p m¼1 Z 2 1 ðE coshðuXÞ þ F sinhðuXÞÞ sinðnYÞdn þ

p

ðA36Þ

0

According to Eq. (12), the head fluctuation in dimensionless X–Y domain is

References Asadi-Aghbolaghi, M., Chuang, M.H., Yeh, H.D., 2012. Groundwater response to tidal fluctuation in a sloping leaky aquifer system. Appl. Math. Model. 36 (10), 4750– 4759. Asadi-Aghbolaghi, M., Chuang, M.H., Yeh, H.D., 2014. Groundwater response to tidal fluctuation in an inhomogeneous coastal aquifer–aquitard system. Water Resour. Manage. 28, 3591–3617. Bear, J., 1972. Dynamics of Fluids in Porous Media. Elsevier, New York. Carr, P.A., van der Kamp, G., 1969. Determining aquifer characteristics by the tidal methods. Water Resour. Res. 5 (5), 1023–1031. Chuang, M.H., Yeh, H.D., 2007. An analytical solution for the head distribution in a tidal leaky confined aquifer extending an infinite distance under the sea. Adv. Water Resour. 30 (3), 439–445. Chuang, M.H., Huang, C.S., Li, G.H., Yeh, H.D., 2010. Groundwater fluctuations in coastal heterogeneous aquifer systems. Hydrol. Earth Syst. Sci. 14, 1819–1826. Dong, L.Y., Chen, J.Y., Fu, C.S., Jiang, H.B., 2012. Analysis of groundwater-level fluctuation in a coastal confined aquifer induced by sea-level variation. Hydrogeol. J. 20 (4), 719–726. Farrell, E.R., 1994. Analysis of groundwater flow through leaky marine retaining structures. Geotechnique 44 (2), 255–263. Ferris, J.G., 1951. Cyclic fluctuations of water level as a basis for determining aquifer transmissibility. Int. Assoc. Sci. Hydrol. Publ. 33, 148–155. Fetter, C.W., 1994. Applied Hydrogeology. Prentice-Hall, Englewood Cliffs, New Jersey. Geng, X.L., Li, H., Boufadel, M.C., Liu, S., 2009. Tide-induced head fluctuations in a coastal aquifer: effects of the elastic storage and leakage of the submarine outlet-capping. Hydrogeol. J. 17 (5), 1289–1296. Google Earth, 2015. Linyuan District, Kaohsiung City. 22.500° N and 120.377° E. Feb. 25, 2015. Guarracino, L., Carrera, J., Vazquez-Sune, E., 2012. Analytical study of hydraulic and mechanical effects on tide-induced head fluctuation in a coastal aquifer system that extends under the sea. J. Hydrol. 450–451, 150–158. Guo, Q.N., Li, H., Boufadel, M.C., Xia, Y.Q., Li, G.H., 2007. Tide-induced groundwater head fluctuation in coastal multi-layered aquifer systems with a submarine outlet-capping. Adv. Water Resour. 30, 1746–1755. Hsu, K.C., Wang, C.H., Chen, K.C., Chen, C.T., Ma, K.W., 2007. Climate-induced hydrological impacts on the groundwater system of the Pingtung Plain, Taiwan. Hydrogeol. J. 15 (5), 903–913. Huang, C.S., Yeh, H.D., Chang, C.H., 2012. A general analytical solution for groundwater fluctuations due to dual tide in long but narrow islands. Water Resour. Res. 48, W05508. http://dx.doi.org/10.1029/2011WR011211. Jacob, C.E., 1950. Flow of groundwater. In: Rouse, H. (Ed.), Engineering Hydraulics. John Wiley, New York, pp. 321–386. Jeng, D.-S., Teo, H.T., Barry, D.A., Li, L., 2005. Two-dimensional approximation for tidal dynamics in coastal aquifers: capillary correction. J. Eng. Mech. 131 (5), 534–541. Kuo, C.H., Wang, C.H., 2001. Implication of seawater intrusion on the delineation of the hydrogeologic framework of the shallow Pingtung Plain aquifer. West. Pac. Earth Sci. 1 (4), 459–472. Li, G.H., Li, H., Boufadel, M.C., 2008. The enhancing effect of the elastic storage of the seabed aquitard on the tide-induced groundwater head fluctuation in confined submarine aquifer systems. J. Hydrol. 350 (1–2), 83–92. Li, H., Jiao, J.J., 2002. Tidal groundwater level fluctuations in L-shaped leaky coastal aquifer system. J. Hydrol. 268, 34–243. Li, H., Jiao, J.J., Luk, M., Cheung, K., 2002. Tide-induced groundwater level fluctuation in coastal aquifers bounded by L-shaped coastlines. Water Resour. Res. 38 (3). http://dx.doi.org/10.1029/2001WR000556. Li, H., Jiao, J.J., Tang, Z.H., 2006. Semi-numerical simulation of groundwater flow induced by periodic forcing with a case-study at an island aquifer. J. Hydrol. 327, 438–446. Li, H., Li, G.Y., Cheng, J.M., Boufadel, M.C., 2007. Tide-induced head fluctuations in a confined aquifer with sediment covering its outlet at the sea floor. Water Resour. Res. 43, W03404. http://dx.doi.org/10.1029/2005WR004724. Li, L., Barry, D.A., Cunningham, C., Stagnitti, F., Parlange, J.-Y., 2000. A two dimensional analytical solution of groundwater response to tidal loading in an estuary and ocean. Adv. Water Resour. 23 (8), 825–833. Li, L., Barry, D.A., Stagnitti, F., Parlange, J.-Y., 1999. Tidal along-shore groundwater flow in a coastal aquifer. Environ. Modell. Assess. 4, 179–188. Liou, T.S., Yeh, H.D., 1997. Conditional expectation for evaluation of risk groundwater flow and solute transport: one-dimensional analysis. J. Hydrol. 199, 378–402.

F.-K. Huang et al. / Journal of Hydrology 530 (2015) 291–305 Millham, N.P., Howes, B.L., 1995. A comparison of methods to determine K in a shallow coastal aquifer. Ground Water 33 (1), 49–57. Nielsen, P., 1990. Tidal dynamics of the water table in beaches. Water Resour. Res. 26 (9), 2127–2134. Sun, H., 1997. A two-dimensional analytical solution of groundwater response to tidal loading in an estuary. Water Resour. Res. 33 (6), 1429–1435. Sun, P.P., Li, H., Boufadel, M.C., Geng, X.L., Chen, S., 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. http://dx.doi.org/10.1029/ 2008WR006893. Teo, H.T., Jeng, D.-S., Seymour, B.R., Barry, D.A., Li, L., 2003. A new analytical solution for water table fluctuations in coastal aquifers with sloping beaches. Adv. Water Resour. 26, 1239–1247.

305

Ting, C.S., 1997. Groundwater Resources Evaluation and Management for the Pingtung Plain, Taiwan. PhD Thesis, Free University, Amsterdam. Trefry, M.G., Bekele, E., 2004. Structural characterization of an island aquifer via tidal methods. Water Resour. Res. 40, W01505. http://dx.doi.org/10.1029/ 2003WR002003. Wang, C., Li, H., Wan, L., Wang, X., Jiang, X., 2014. Closed-form analytical solutions incorporating pumping and tidal effects in various coastal aquifer systems. Adv. Water Resour. 69, 1–12. Xia, Y.Q., Li, H., Boufadel, M.C., Guo, Q.N., Li, G.H., 2007. Tidal wave propagation in a coastal aquifer: effects of leakage through its submarine outlet and offshore roof. J. Hydrol. 337, 249–257. Yeh, H.D., Chang, Y.C., Zlotnik, V.A., 2008. Stream depletion rate and volume from groundwater pumping in wedge-shaped aquifers. J. Hydrol. 349 (3–4), 501–511.