Interactive comment on “ Groundwater fluctuations in heterogeneous coastal leaky aquifer systems ”

The article presents a mathematical derivation of a simple equation to model tidal wave propagation in different subsurface geometries when leakage may take place between an aquifer and its confining unit. After deriving this equation, the authors present a number of different tests where they compare tidal wave propagation in different hydraulic settings, to help understand qualitatively how the model behaves. Essentially, they carry out a sensitivity analysis to the different parameters in their equation. Personally I like this approach, as it helps to understand the equations qualitatively. I think therefore the most important merit of this paper is this sensitivity analysis.


Introduction
The groundwater near the coast may fluctuate with periodical tides.Such a phenemenon is an interesting and practical issue for hydrogeologists.The coastal leaky aquifer is referred to an aquifer system consisting of an upper unconfined aquifer bounded from below by one or more lower aquifers (Chen and Jiao, 1999).Many researches revealed that the head fluctuation in the unconfined aquifer due to tides is significantly damped by its storage (e.g., Millham and Howes, 1995;Chen and Jiao, 1999).White and Roberts (1994) indicated that the farthest distance that tidal propagation can reach, defined as the intrusion distance, is usually not over Correspondence to: H.-D. Yeh (hdyeh@mail.nctu.edu.tw)30 m if the head fluctuation is heavily damped.The water table is therefore considered to maintain constant.
Several previous researches indicated that the effect of the leakage on the groundwater in coastal leaky aquifer systems is significant.Jiao and Tang (1999) developed an analytical solution for describing the head fluctuation in coastal confined aquifers, where the aquitard infinitely extends from the coastal line to the inland.They indicated that the leakage diminishes the amplitude of the head fluctuation and the intrusion distance.Li and Jiao (2001) presented an analytical solution for describing leaky aquifer systems with considering the leakage and the effect of aquitard storage.They found that both the leakage and the effect of aquitard storage on the head fluctuation in the lower aquifer are negligible when the aquitard storage is small and the storage ratio of the aquitard to the underlying aquifer is less than 0.5.Guo et al. (2010) developed an analytical solution for describing the groundwater head response to tidal fluctuation in a coastal aquifer consisting of two zones with different hydraulic properties.The solution was used to investigate the behaviors of amplitude and phase shift of the head fluctuation in the lower aquifer.In addition, the solution was also used to estimate the aquifer parameters in a real coastal two-zone aquifer.
The variations in depositional and post depositional processes over a very long period of time may result in heterogeneous hydraulic conductivity and nonuniform thickness of the aquitard (Cherry et al., 2006) and thus produce various amounts of leakage to the lower aquifer.The objective of this paper is to develop a mathematical model for describing the head fluctuation in the lower aquifer due to the tidal effect in a coastal leaky aquifer system.This system is divided into several horizontal regions for the heterogeneous aquitard or lower heterogeneous aquifer.The head of the upper unconfined aquifer is assumed constant.The analytical solution of the model is developed by the direct Fourier method.The effects of the length and location for a discontinuous aquitard on the amplitude and phase shift of the head fluctuation in 14 Figure 1.Schematic diagram of a heterogeneous coastal leaky aquifer system the lower aquifer are investigated.The present solution is applied to simulate the head distribution for the coastal aquifer system in the Chek Lap Kok (CLK) airport, Hong Kong (Jiao and Tang, 1999) and the results are compared with the observed data taken at this airport.In addition, the influence of formation heterogeneity on the spatial head distribution in the lower aquifer is also examined.

Conceptual model
Figure 1 shows a coastal leaky aquifer system in which the aquitard and underlying aquifer are heterogeneous in the horizontal direction.The system is divided to N different horizontal regions and the formation materials of the aquitard and underlying aquifer in each region are homogeneous.The origin of the coordinates in the x direction is located at the coastal line.The distance measured from the origin to the interface between regions 1 and 2 is d 1 .Similarly, the distance measured from the origin to the interface between regions n − 1 and n is d n−1 .The mean sea level (MSL) is chosen as reference datum.Consider that the unconfined aquifer may store a large amount of water which can damp the tidal effect so that the water table fluctuation in the unconfined aquifer is negligible in comparison with that in the lower one.Therefore, the water table is assumed to be the MSL.
The governing equation describing the head distributions in these regions is expressed as where n is an integer from 1, 2, 3 . . .N; h n , T n and S n are the hydraulic head, transmissivity, and storativity for the lower aquifer in n-th region, respectively; and L n is the leakage for the aquitard in n-th region.The tidal boundary is expressed as where A and ω are the amplitude and frequency of the tide, respectively.The remote boundary condition is where h N is the hydraulic head for the farthest lower aquifer from the coastal line.The continuities of the head and flux required at the interfaces d n are respectively and (5)

Solution for heterogeneous aquifer systems
Based on the direct Fourier method, the solutions for describing the head distributions in the regions demonstrated in Fig. 1 can be expressed as where Re means the real part of the complex expression and i is equal to √ −1.Substituting Eq. ( 6) into Eq.( 1) and eliminating the exponential terms can obtain ordinary differential equations in terms of X n (x).Substituting Eq. ( 6) into Eq.( 2) with n = 1 and eliminating the exponential terms acquire the boundary condition expressed as X 1 (0) = 1.Similarly, the equation lim x→∞ ∂X N /∂x = 0 can be obtained from Eqs. (3) and (6).The continuity requirements (i.e., X n = X n+1 and T n ∂X n /∂x = T n+1 ∂X n+1 /∂x) at x = d n can also be acquired from Eqs. (4)-(6) in a similar manner.The general solutions for X n (x) can then be solved as where the parameter a n = √ ωS n /2T n is a reciprocal of the decay length for the n-th lower aquifer and u n = L n /ωS n is the dimensionless leakage.The coefficients, c1 n and c2 n , for n-th lower aquifer can then be determined by the boundary conditions X 1 (0) = 1 and lim x→∞ ∂X N /∂x = 0 as well as the continuity requirements X n = X n+1 and T n ∂X n /∂x = T n+1 ∂X n+1 /∂x.
The equations for solving c1 n and c2 n are expressed in matrix form as and where both C and F are (N −1)×N matrixes and identical to matrixes B and E, respectively, except having a minus sign before the exponent in the exponential terms.
Based on Cramer's rule, the coefficients c1 n and c2 n can be determined as and where det is an abbreviation for determinant of the matrix; D n and D N+n can be obtained by replacing the n-th and (N + n)-th columns of matrix D, respectively, with a column matrix that the top element is one and the others are zero.Define the variable a X n and b X n representing the real and imaginary parts of the result calculated from Eq. ( 7), respectively, and thus one can have = cos(ωt) − i sin(ωt) into Eq.( 6) and taking the real part of Eq. ( 6) lead to the result in terms of cosine function as and where φ n is the phase shift and c n is the amplitude coefficient.

Special cases
Jiao and Tang (1999) presented an analytical solution for describing the head fluctuation in a coastal leaky aquifer system.Both the aquifer and aquitard in their system are homogeneous and of infinite extent from the coastal line.Their solution (1999, Eq. 4) can be obtained from the present solution for N = 1 and written, in our notation, as where and In addition, Guo et al. (2010) presented an analytical solution to describe the head fluctuations in a coastal aquifer system consisting of two zones which have different hydraulic properties in the horizontal direction.When N = 2 and u 1 = u 2 = 0, the present solution reduces to Guo et al.'s (2010) solution (2010, Eq. 7a for the first aquifer and Eq.7b for the second aquifer) which is expressed, in our notation, as where e a 1 (2d 1 −x) + e a 1 (2d 1 +x) + e 4d 1 a 1 − e 2a 1 x sin(a 1 x) + e 4d 1 a 1 − e 2a 1 x sin(a 1 x) Accordingly, the solutions derived by Jiao and Tang (1999) and Guo et al. (2010) can be considered as our special cases.

The effect of aquitard length
Consider the case that the leaky aquifer system with a homogeneous lower aquifer, overlain by a heterogeneous aquitard, can be divided into two regions (i.e., regions 1 and 2).The aquitard in region 1 has semi-permeable with leakage u 1 = 5 and a finite length d 1 while the aquitard in region 2 is impermeable (u 2 = 0) and of infinite extent in the x direction.The purpose of this case is to investigate the effect of aquitard length on the amplitude and phase shift of the head fluctuation in the lower aquifer.Figure 2 shows the curves of normalized amplitude of the head fluctuation versus dimensionless inland distance for a 1 = a 2 = 1 × 10 −2 m −1 and aquitard length d 1 of 0, 30, 50, 300, 850 m.The normalized amplitude is defined as the ratio of the amplitude of the head fluctuation to that of the tide.The data represented by the open circle is obtained by the present solution with no leakage (i.e., u 1 = u 2 = 0).The curves shown in Fig. 2 indicate that a larger d 1 has a smaller normalized amplitude and a shorter tidal intrusion distance.In addition, the present solution approaches Jiao and Tang's(1999) solution for u 1 = 5 and a 1 = 1×10 −2 m −1 when d 1 goes large (say d 1 = 850 m).This is because the aquitard length is larger than the tidal intrusion distance.

Effect of aquitard location
In contrast to the previous case, the aquitard in region 1 with a length d 1 is now treated as impermeable, i.e., u 1 = 0, while the aquitard in region 2 has u 2 = 5.The aquitard and lower aquifer in region 2 are considered to be of infinite extent.The purpose of this case is to study the effect of the value of d 1 on the amplitude and phase shift of the head fluctuation in the lower aquifer.show that those dashed lines approach the line with large d 1 , indicating that the effects of the aquitard on the amplitude and phase shift decrease with increasing d 1 .Once the d 1 is larger than the tidal intrusion distance, the effects of the leakage on both amplitude and phase shift become negligible.

Effect of heterogeneity in aquitard or aquifer
Consider the case that a leaky aquifer system has a heterogeneous aquitard.The aquitard is divided into three regions and the hydraulic conductivity of the aquitard in each region is homogeneous.The lower aquifer is considered to be homogeneous in order to investigate the effect of aquitard heterogeneity on the aquifer head fluctuation.Figure 6 shows the spatial head distribution of the lower aquifer with a 1 = a 2 = a 3 = 5×10 −3 m −1 for the aquitard having four different scenarios at the high-tide period (i.e., ωt = 2π).The aquitard is heterogeneous in the first two scenarios and homogenous but with different values of hydraulic conductivity in the other two scenarios.The normalized hydraulic head is defined as the ratio of the hydraulic head of the lower aquifer to the amplitude of the tide.For scenario 1, the leakages in regions 1-3 are chosen as u 1 = 10, u 2 = 5 and u 3 = 1, respectively.In contrast, the leakages in these three regions are selected as u 1 = 1, u 2 = 5 and u 3 = 10 for scenario 2. In addition, the leakages are chosen as u 1 = u 2 = u 3 = 1 for scenario 3 and u 1 = u 2 = u 3 = 10 for scenario 4. The head distribution curve of scenario 1 is below that of scenario 3 and close to that of scenario 4, indicating that a large leakage occurring near the coastal line (region 1) significantly diminishes the hydraulic head.On the other hand, the curve of scenario 2 in region 1 is below that of scenario 3 even when the leakage in region 1 of scenario 2 is the same as that of scenario 3.This is because large leakages occurring in regions 2 and 3 will influence the hydraulic head distributions in regions 1 and 2. Accordingly, the aquitard heterogeneity has significant effect on the hydraulic head distribution.
Alluvial fans formed at the base of mountains usually have coarser sediment at the upstream part of the fan and finer sediment at the downstream part.The formation of the coastal leaky aquifer in the alluvial fan may exhibit the phenomenon of trending heterogeneity (Freeze and Cherry, 1979).The lower aquifer is divided into three regions for simulating trending heterogeneity (i.e., T 1 = 10 m 2 /day, T 2 = 50 m 2 /day and T 3 = 100 m 2 /day).The leakages in these three regions are chosen the same for assessing the effect of aquifer heterogeneity.Figure 7 shows the spatial head distributions at ωt = 2π and ωt = 0.5π for the homogeneous aquifer and heterogeneous trending aquifer.Following parameter values are used in the analyses: u 1 = u 2 = u 3 = 5, S 1 = S 2 = S 3 = 10 −4 , d 1 = 100 m, d 2 = 200 m and ω = 2π/0.5day −1 .The hydraulic conductivities are considered as T 1 = T 2 = T 3 = 50 m 2 /day for the homogeneous aquifer and T 1 = 10 m 2 /day, T 2 = 50 m 2 /day and T 3 = 100 m 2 /day for the trending aquifer.The figure indicates that the trending aquifer has a smaller tidal intrusion distance than the homogeneous one.This is because the region 1 has a smaller hydraulic conductivity.In addition, the slopes of the normalized head distribution are markedly different near x = 100 m and 200 m because the hydraulic conductivities in these regions are different.Obviously, the aquifer heterogeneity also has an impact on the hydraulic head distribution.

Comparison of present solution with observed data
Figure 8 shows a geological cross section of the leaky aquifer system in the CLK airport, Hong Kong (Jiao and Tang, 1999).The lower aquifer is homogeneous and lies on the The thickness of the aquitard in region 1 is about a half of that in region 2; therefore, the leakage in region 1 is about two times that in region 2. Jiao and Tang (1999) applied their solution to simulate the head fluctuation for the aquifer system in the CLK airport.The observed well was located at 271 m distance from the coastal line.The lower aquifer is considered to be homogeneous with a 1 = 7.65 × 10 −3 m −1 .The thickness of the aquitard is assumed constant in the horizontal direction and the leakage is chosen as u 1 = 9.38 × 10 −3 .The amplitude and frequency of the tide are chosen as A = 0.8 m and ω = 2π/0.5day −1 , respectively.
Figure 9 shows the observed data taken from the aquifer system in the CLK airport, Hong Kong and the present solution computed with the parameters d 1 = 300 m, a 1 = a 2 = 7.65×10 −3 m −1 , ω = 2π/0.5day −1 , A = 0.8 m, u 1 = 9.38× 10 −3 and u 2 = 0.5u 1 .This figure shows that the present solution is identical to Jiao and Tang's (1999) solution and has a good agreement with the observed data.Such a result can be attributed to the fact that both leakages in regions 1 and 2 are very small and the length of the thin aquitard exceeds the tidal intrusion distance (i.e., 150 m), which is estimated based on Eq. ( 18).

Concluding remarks
Based on the direct Fourier method, a one-dimensional analytical model is developed for describing the head fluctuation due to the tidal effect in a coastal leaky aquifer system.The aquitard and underlying aquifer can be considered to be heterogeneous and divided into several regions and each region has a homogenous hydraulic conductivity.The water table of the upper unconfined aquifer is assumed constant while the hydraulic head of the lower aquifer is subject to the effect of tidal fluctuation.The solution developed by Jiao and Tang (1999) can be considered as a special case of the present solution when N = 1 or the length of the first aquifer is larger than the tidal intrusion distance.The present  (Jiao and Tang, 1999).Based on the present solution, the effects of formation heterogeneity as well as the length and location of the aquitard on the head fluctuation in the lower aquifer can be concluded as follows: 1.The aquitard with a large leakage and/or long length near the coastal line will diminish the amplitude of the head fluctuation in the lower aquifer.
2. The phase shift of the head fluctuation in the lower aquifer increases slowly with the inland distance for the case of large leakage and dramatically for the case of no leakage.
3. A lower aquifer with trending heterogeneity in a coastal leaky aquifer system usually has a small tidal intrusion distance and low head distribution.

Figure 6 .Fig. 6 .
Figure 6.The spatial head distributions in the lower aquifer with heterogeneous trending aquitards for scenarios 1 and 2 and with homogeneous aquitards but having different conductivity values for scenarios 3 and 4 at high-tide period ( π ω 2 = t ) Fig. 6.The spatial head distributions in the lower aquifer with heterogeneous trending aquitards for scenarios 1 and 2 and with homogeneous aquitards but having different conductivity values for scenarios 3 and 4 at high-tide period (ωt = 2π).

Fig. 8 .
Figure 8.A geological cross section in the CLK Airport, Hong Kong (modified from Jiao and Tang, 1999) Fig. 8.A geological cross section in the CLK Airport, Hong Kong (modified from Jiao and Tang, 1999).

Figure 9 .Fig. 9 .
Figure 9. Predicted hydraulic heads by the present solution and Jiao and Tang's solution 2 (1999) and the observed data measured at the CLK airport, Hong Kong 3 4 5 6