Analytical Study on the Tidal Wave Propagation Coupled with Vertical Leakage from a Vertically Heterogeneous Aquitard in a Coastal Leaky Aquifer System

Chao Zhuang; Xiangdong Xu; Zhi Dou; Chun Zhu; Yun Yang; Jinguo Wang; Zhifang Zhou; Xiaoding Xu

doi:10.2113/2022/5856296

Abstract

Analytical studies on the tidal wave propagation in a coastal leaky aquifer commonly assume that the low-permeability aquitard is homogeneous. The aquitard is, however, vertically heterogeneous in nature due to varying soil types along the burial depth as can be frequently seen from borehole logs. In this study, an innovative analytical model is developed to explore the vertical heterogeneity in aquitard hydraulic conductivity (⁠ $K$ ⁠) and specific storage (⁠ $S_{s}$ ⁠) on the tidal wave propagation coupling with vertical leakage from the aquitard. The novelty behind the newly derived analytical solution is attributed to the fact that the aquitard along the burial depth can be divided into any number of homogeneous zones, each of which is associated with distinct $K$ and $S_{s}$ values, so that arbitrary vertical heterogeneity pattern of aquitard $K$ and $S_{s}$ can be captured. Theoretical analysis results reveal that an interlayer within the aquitard with a smaller $K$ can significantly enhance the amplitude and phase shift of periodic groundwater head fluctuations in the leaky aquifer, while a larger $S_{s}$ of the interlayer will weaken the amplitude and phase shift provided a relatively small aquitard $K$ ⁠. Buried locations of the interlayer also implement nonnegligible effects on the tidal wave propagation. For the scenario of exponentially decaying aquitard $K$ and $S_{s}$ ⁠, which is commonly encountered for a thick aquitard, a larger decay exponent results in smaller amplitude and phase shift. This analytical study highlights the importance of vertical aquitard heterogeneity on tidal wave propagation in a coastal leaky aquifer system.

1. Introduction

Sea tide can cause periodic groundwater level fluctuations in coastal leaky aquifers. This phenomenon can be termed tidal wave propagation and has attracted wide attention in the research area of coastal hydrogeology (e.g., [1, 2–4]). Earlier analytical models for the tidal wave propagation mainly considered the aquifer system to be homogeneous (e.g., [4–6]). In fact, the heterogeneity is a common place in lithosphere including the coastal leaky aquifer system.

Alluvial plains are often associated with spatial variability of sediments in terms of soil type and grain size. Particularly, 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 [7, 8]. To deal with such a horizontal heterogeneity in a coastal leaky aquifer system, Trefry [9] developed an analytical solution for a finite composite coastal aquifer, providing prescribed head, prescribed flux, and mixed inland boundary conditions. Trefry and Bekele [10] then applied this analytical solution to handle tidal wave propagation signals monitored at multiple wells placed in a finite-length island aquifer system of Garden Island, Australia. Monachesi and Guarracino [11] and Wang et al. [12] analytically investigated the effects of horizontally linear heterogeneity in aquifer hydraulic conductivity on the tidal wave propagation.

Groundwater leakage from the aquitard has been found to have critical effects on the tidal wave propagation [13]. In nature, the aquitard usually consists of relatively fine sediments associated with lower permeability while larger storativity comparing to sandy aquifers [14]. Consequently, the volume of groundwater released from aquitard storage is a significant portion of groundwater leakage that recharges the aquifer [15]. From this perspective, vertical heterogeneity in aquitard hydraulic properties (i.e., the hydraulic conductivity $K$ and specific storage $S_{s}$ ⁠) can implement nonnegligible effects on tidal wave propagation in the adjacent leaky aquifer. However, currently available analytical models for coastal leaky aquifer system have not dealt with the vertical heterogeneity in aquitard $K$ and $S_{s}$ ⁠.

Vertical heterogeneity in aquitard layers is rather common in the lithosphere and can be directly identified according to lithology analysis data. Generally speaking, the aquitard can be divided into multiple zones along the burial depth profile, and each subunit has a distinct soil type associated with independent $K$ and $S_{s}$ ⁠. Such a type of vertical aquitard heterogeneity can also be identified via interpreting in situ pressure and compaction data [16–19].

In addition to vertical heterogeneity of multiple zones, the aquitard can also contain a sand layer, gravel interlayer, and cracks. In other words, there could exist local anomaly of the soil type in the aquitard, leading to local heterogeneity in aquitard $K$ and $S_{s}$ ⁠. Even the small-scale heterogeneity could significantly control groundwater flow in the subsurface. For instance, Meriano and Eyles [20] found that the existence of aquitard interlayers played an important role in supplying the aquifer below. Gerber et al. [21] found that the existence of crack and sand gravels in the aquitard produced abnormal vertical groundwater flow velocity in the aquitard. Harrington et al. [22] studied the effect of sandy interlayer in the aquitard on solute transport and developed a model to explain the radial component of solute diffusion when the cross section of sandy interlayer is circular or elliptical, and the model was successfully applied to explain the source of high chloride ion concentration in the clay-rich area in southern Saskatchewan, Canada.

At the regional scale, the porosity of geologic media usually shows a decreasing trend with the burial depth, which can also be indirectly manifested as the depth-decreasing hydraulic properties [23–26]. Many studies have reported that the aquifer $K$ exponentially decays with burial depth [15, 27–30]. Comparatively speaking, a thick aquitard could also demonstrate exponentially decaying hydraulic properties with the burial depth. Ferris et al. [31] measured the aquitard $K$ values at 15 sites in south-central Saskatchewan, Canada, and further compiled the $K$ data through 21 studies to characterize the hydrogeology of aquitard on a regional scale. Their results indicated that the aquitard $K$ followed an exponentially decaying pattern with depth index, and the variation in aquitard $K$ reached 5 orders of magnitude in the range of 80 m. Previous studies also found that the aquifer $S_{s}$ also decreases with depth [32–34]. Recently, Kuang et al. [35] proposed an empirical model of aquifer $S_{s}$ exponentially decaying with depth and also concluded that $S_{s}$ has an exponential decay with depth in the shallow crust. Without losing generality, we believe that the aquitard $S_{s}$ could also exponentially decay with depth at a regional scale and from the statistical perspective, even though few literatures provided relevant evidence.

The objective of this study is to explore the influences of vertically heterogeneous aquitard $K$ and $S_{s}$ on the tidal wave propagation in a coastal leaky aquifer system. Firstly, we construct a quasi-2D model for a horizontally extended coastal leaky aquifer system with a finite length and with the consideration of depth-dependent aquitard $K$ and $S_{s}$ ⁠. Secondly, the analytical solution associated with tidal-head inland boundary condition will be derived. Lastly, tidal wave propagation bias of the groundwater level fluctuations caused by the vertical heterogeneity of aquitard $K$ and $S_{s}$ will be discussed.

2. Mathematical Model and Analytical Solutions

2.1. Mathematical Model

Here, we consider a horizontally extended coastal leaky aquifer system as shown in Figure 1. To consider the unconfined aquifer overlain, the aquitard of interest has a large specific yield to effectively damp the tidal effect so that the water table in the unconfined aquifer almost remains constant (e.g., [3, 13]). Besides, the difference in density between seawater and fresh groundwater is neglected. It has been verified that ignoring the density difference when one derives the analytical solution results in an error not greater than 2.5% of the tidal amplitude and the error decreases landward [36].

Let the $x$ -axis be positive in the landward direction and coincide with the contact interface between the aquitard and leaky confined aquifer. Let the $z$ -axis be vertical and positive upward. The origin of the Cartesian coordinate system is located at the coastline (Figure 1). The tidal boundary is located at $x = 0$ and the inland-side boundary is located at $x = l$ ⁠. Vertical Cartesian coordinate of the zone interfaces inside the aquitard are indicated by $z_{1}, z_{2}, \dots, z_{n}, \dots, z_{N}$ ⁠. The aquitard $K$ is at least two orders of magnitude smaller than the underlying confined aquifer, so that the groundwater flow is essentially vertical within the aquitard (e.g., [2, 37]). In comparison, groundwater flow is horizontal in the confined aquifer.

Under the above assumptions, the governing equations for groundwater flow within the confined aquifer and the aquitard are

\begin{matrix} (1) & S \frac{\partial h (x, t)}{\partial t} = T \frac{\partial^{2} h (x, t)}{\partial x^{2}} + K_{1}^{'} {\frac{\partial h_{1}^{'} (x, z, t)}{\partial z}|}_{z = 0}, \\ (2) & S_{s, n}^{'} \frac{\partial h_{n}^{'} (x, z, t)}{\partial t} = K_{n}^{'} \frac{\partial^{2} h_{n}^{'} (x, z, t)}{\partial z^{2}}, z_{n - 1} \leq z \leq z_{n}, \end{matrix}

where

h = h (x, t)

is the aquifer hydraulic head [

L

] at location

x [L]

and time

t [T]

⁠,

T [L^{2} T^{- 1}]

and

S [‐]

are the transmissivity and storage coefficient of the leaky confined aquifer, respectively;

h_{n}^{'} = h_{n}^{'} (x, z, t)

is the aquitard hydraulic head at location

(x, z)

and time

t

corresponding to zone

n

⁠;

K_{n}^{'} [{L T}^{- 1}]

and

S_{s, n}^{'} [L^{- 1}]

are the vertical hydraulic conductivity and specific storage of the aquitard corresponding to zone

n

⁠, respectively.

Based on the Euler formula, the tidal boundary for the confined aquifer at

x = 0

as shown in Figure 1 is expressed as

\begin{matrix} (3) & {h (x, t)|}_{x = 0} = A^{0} \exp [i (w t + c)], \end{matrix}

where

i = \sqrt{- 1}

⁠,

A^{0}

⁠,

w

⁠, and

c

are the amplitude [

L

], angular frequency [

T^{- 1}

], and initial radian phase [-] of the tidal fluctuations, respectively.

We also consider a tidal-head boundary at the inland side (i.e.,

x = l

⁠). Here, we consider different amplitudes and initial phase shift at the inland-side boundary comparing with that at the coastline. Therefore, one has

\begin{matrix} (4) & {h (x, t)|}_{x = l} = A^{l} \exp (i θ) \exp [i (w t + c)], \end{matrix}

where

i = \sqrt{- 1}

⁠,

h_{msl}^{l}

⁠,

A^{l}

⁠, and

θ

are the mean sea level [

L

], amplitude [

L

], and initial phase difference [-] at the inland-side boundary, respectively. Note that Equation () is typically representative for an island responding to dual tides [10, 38, 39], where

A^{l} = A^{0}

seems to coincide with the real-world situation. If the hydraulic head does not tidally fluctuate (i.e.,

A^{l} = 0

⁠), Equation () degenerates to a constant-head inland boundary condition, which could be suitable for a horizontally extensive coastal aquifer system.

In terms of the aquitard, at the interface between zones

n

and

n + 1

(i.e.,

z = x_{n}

⁠), the groundwater flow within the aquitard has to meet the continuity requirement of head and flux. As a consequence, one has

\begin{matrix} (5) & {h_{n}^{'} (x, z, t)|}_{z = z_{n}} = {h_{n + 1}^{'} (x, z_{n}, t)|}_{z = z_{n}}, \\ {K_{n}^{'} \frac{\partial h_{n}^{'} (x, z, t)}{\partial z}|}_{z = z_{n}} = {K_{n + 1}^{'} \frac{\partial h_{n + 1}^{'} (x, z, t)}{\partial z}|}_{z = z_{n}} . \end{matrix}

In addition, at the interface between the aquitard and the underlying confined aquifer of interest, as well as the constant-head overlying unconfined aquifer, hydraulic head is also continuous; thus, one has

\begin{matrix} (6) & {h_{1}^{'} (x, z, t)|}_{z = 0} = h (x, t), \\ (7) & {h_{N}^{'} (x, z, t)|}_{z = x_{N}} = 0 . \end{matrix}

2.2. Analytical Solution

Based on the separation-of-variable method, which has been widely employed to handle analytical models of tide-induced aquifer head fluctuations (e.g., [4, 39]), the periodic hydraulic head of the confined aquifer, as well as that of the aquitard, can be, respectively, expressed in the following complex form as

\begin{matrix} (8) & h (x, t) = h^{p} (x) \exp [i (w t + c)], \\ (9) & h_{n}^{'} (x, z, t) = {h^{'}}_{n}^{p} (x, z) \exp [i (w t + c)], \end{matrix}

where

h^{p} (x)

contains the information of amplitude and phase shift of periodic aquifer hydraulic head fluctuations at the horizontal location

x

⁠, and similarly,

{h^{'}}_{n}^{p} (x, z)

contains the information of amplitude and phase shift of periodic aquitard hydraulic head fluctuations at the horizontal location

x

and vertical location

z

within the zone

n

⁠. It is easy to find that the derivation of the analytical solution for

h (x, t)

is equivalent to the derivation of

h^{p} (x)

⁠.

Substituting Equation () into Equation () yields

\begin{matrix} (10) & \frac{\partial^{2} h^{p} (x)}{\partial x^{2}} + \frac{K_{1}^{'}}{T} {\frac{\partial {h^{'}}_{1}^{p} (x, z)}{\partial z}|}_{z = 0} - \frac{i w S}{T} h^{p} (x) = 0 . \end{matrix}

Substituting Equation () into Equation () yields

\begin{matrix} (11) & \frac{\partial^{2} {h^{'}}_{n}^{p} (x, z)}{\partial z^{2}} - ξ_{n}^{2} {h^{'}}_{n}^{p} (x, z) = 0, \end{matrix}

where

ξ_{n} = \sqrt{i w S_{s, n}^{'} / K_{n}^{'}}

⁠. The general solution of Equation () is

\begin{matrix} (12) & {h^{'}}_{n}^{p} (x, z) = E_{n} (x) \cosh (ξ_{n} z) + F_{n} (x) \sinh (ξ_{n} z) . \end{matrix}

The matching conditions at the interface between zone

n

and zone

n + 1

are

\begin{matrix} (13) & {{h^{'}}_{n}^{p} (x, z)|}_{z = z_{n}} = {{h^{'}}_{n + 1}^{p} (x, z)|}_{z = z_{n}}, \\ (14) & K_{n}^{'} {\frac{\partial {h^{'}}_{n}^{p} (x, z)}{\partial z}|}_{z = z_{n}} = K_{n + 1}^{'} {\frac{\partial {h^{'}}_{n + 1}^{p} (x, z)}{\partial z}|}_{x = x_{n}} . \end{matrix}

Substituting Equations () and () into Equation () yields

\begin{matrix} (15) & [\begin{matrix} E_{n} (x) \\ F_{n} (x) \end{matrix}] = Ψ_{n} [\begin{matrix} E_{n + 1} (x) \\ F_{n + 1} (x) \end{matrix}], \end{matrix}

where

\begin{matrix} (16) & Ψ_{n} = {[\begin{matrix} \cosh (ξ_{n} z_{n}) & \sinh (ξ_{n} z_{n}) \\ K_{n}^{'} ξ_{n} \sinh (ξ_{n} z_{n}) & K_{n}^{'} ξ_{n} \cosh (ξ_{n} z_{n}) \end{matrix}]}^{- 1} \cdot [\begin{matrix} \cosh (ξ_{n + 1} z_{n}) & \sinh (ξ_{n + 1} z_{n}) \\ K_{n + 1}^{'} ξ_{n + 1} \sinh (ξ_{n + 1} z_{n}) & K_{n + 1}^{'} ξ_{n + 1} \cosh (ξ_{n + 1} z_{n}) \end{matrix}] . \end{matrix}

Thus, one has

\begin{matrix} (17) & [\begin{matrix} E_{1} (x) \\ F_{1} (x) \end{matrix}] = Θ_{n - 1} [\begin{matrix} E_{n} (x) \\ F_{n} (x) \end{matrix}], \end{matrix}

where

\begin{matrix} (18) & Θ_{n} = Π_{k = 1}^{n} Ψ_{k} = Ψ_{1} Ψ_{2} \dots Ψ_{n} . \end{matrix}

It is noteworthy that if the aquitard is vertically homogeneous,

Ψ_{n}

and

Θ_{n}

will become the identity matrix and one has

E_{1} (x) = E_{n} (x)

and

F_{1} (x) = F_{n} (x)

⁠. In order to handle the zonational vertical heterogeneity of the aquitard, i.e., to determine

E_{n} (x)

and

F_{n} (x)

in Equation (), the boundary conditions for zones 1 and

N

(i.e., Equations () and ()) must be incorporated. In terms of zone 1, substituting Equations () and () into Equation () yields

\begin{matrix} (19) & {{h^{'}}_{1}^{p} (x, z)|}_{z = 0} = E_{1} (x) = h^{p} (x) . \end{matrix}

In terms of zone

N

⁠, substituting Equations () and () into Equation () yields

\begin{matrix} (20) & {h^{'}}_{N}^{p} (x, z_{N}) = E_{N} (x) \cosh (ξ_{n} z_{N}) + F_{N} (x) \sinh (ξ_{N} z_{N}) = 0 . \end{matrix}

Substituting Equation () into Equations () and (), one can further derive

\begin{matrix} (21a) & E_{N} (x) = - \frac{\tanh (ξ_{N} z_{N})}{Θ_{N - 1} (1, 2) - \tanh (ξ_{N} z_{N}) Θ_{N - 1} (1, 1)} h^{p} (x), \\ (21b) & F_{N} (x) = \frac{1}{Θ_{N - 1} (1, 2) - \tanh (ξ_{N} z_{N}) Θ_{N - 1} (1, 1)} h^{p} (x), \\ (21c) & F_{1} (x) = \frac{Θ_{N - 1} (2, 2) - \tanh (ξ_{N} z_{N}) Θ_{N - 1} (2, 1)}{Θ_{N - 1} (1, 2) - \tanh (ξ_{N} z_{N}) Θ_{N - 1} (1, 1)} h^{p} (x) . \end{matrix}

After substituting Equations () and () into Equation (), one can further obtain

\begin{matrix} (22) & {\frac{\partial {h^{'}}_{1}^{p} (x, z)}{\partial z}|}_{z = 0} = ξ_{1} \frac{Θ_{N - 1} (2, 2) - \tanh (ξ_{N} z_{N}) Θ_{N - 1} (2, 1)}{Θ_{N - 1} (1, 2) - Θ_{N - 1} (1, 1) \tanh (ξ_{N} z_{N})} h^{p} (x) . \end{matrix}

Substituting Equation () into () yields

\begin{matrix} (23) & \frac{\partial^{2} h^{p} (x)}{\partial x^{2}} - η^{2} h^{p} (x) = 0, \end{matrix}

where

\begin{matrix} (24) & η^{2} = \frac{i w S}{T} - \frac{K_{1}^{'} ξ_{1}}{T} \frac{Θ_{N - 1} (2, 2) - \tanh (ξ_{N} z_{N}) Θ_{N - 1} (2, 1)}{Θ_{N - 1} (1, 2) - \tanh (ξ_{N} z_{N}) Θ_{N - 1} (1, 1)} . \end{matrix}

The general solution of Equation () is

\begin{matrix} (25) & h^{p} (x) = C \cosh (η x) + D \sinh (η x), \end{matrix}

where

C

and

D

are constants that need to be determined by incorporating the boundary conditions at the coastline and inland side. Substituting Equations () and () into Equation () yields

\begin{matrix} (26) & C = A^{0}, \\ D = A^{l} \exp (i θ) csch (η l) - A^{0} \coth (η l) . \end{matrix}

As a result, the analytical solution for

h^{p} (x)

is

\begin{matrix} (27) & h^{p} (x) = A^{0} \cosh (η x) + [A^{l} \exp (i θ) csch (η l) - A^{0} \coth (η l)] \sinh (η x) . \end{matrix}

The normalized amplitude

α

(i.e., the ratio of the amplitude of the periodic head fluctuation to that of the tide) and radian phase shift (⁠

τ

⁠) of aquifer head fluctuations are computed as

\begin{matrix} (28) & α (x) = \frac{|h^{p} (x)|}{A^{0}}, \\ τ (x) = \arctan [\frac{Im (h^{p} (x))}{Re (h^{p} (x))}], \end{matrix}

where

|h^{p} (x)|

is the modulus of

h^{p} (x)

that represents the amplitude at the horizontal location

x

⁠, and Im and Re denote the imaginary and real parts of the expression, respectively.

3. Effects of Vertical Aquitard Heterogeneity on Tidal Wave Propagation

Previous studies have figured out the effects of horizontal aquifer heterogeneity on the tide-induced aquifer head fluctuations. This study mainly focuses on analyzing the effects of vertical heterogeneity in the aquitard $K$ and $S_{s}$ ⁠. Here, we consider a coastal leaky aquifer system of a long but narrow island with a horizontal extent of 500 m. To the best of our knowledge, the analytical model with the consideration of vertically heterogeneous aquitard $K$ and $S_{s}$ has not been yet developed for an island leaky aquifer system. By referring to the ranges of hydraulic parameters in leaky confined systems summarized by Li and Jiao [13], the transmissivity and storage coefficient of the confined aquifer are set to be $T = 1.5 \times 1 0^{2}$ m²/d and $S = 1 \times 1 0^{- 4}$ ⁠, respectively. The aquitard thickness is set to be 10 m. For the tidal waves, the amplitude is $A_{0} = A_{N} = 1$ m, the initial phase shift is $c = 0$ ⁠, the phase difference is $θ = 0$ ⁠, and the angular frequency is $w = 6.283$ d^-1 corresponding to the period of 1 day. Two scenarios of vertical aquitard heterogeneity are considered in the following analysis. The first scenario is related to the aquitard interlayer, which is associated with $K$ and $S_{s}$ anomaly, indicating the local heterogeneity. Another scenario is related to the exponentially decaying aquitard $K$ and $S_{s}$ ⁠, indicating the generally varying heterogeneity at the regional scale.

3.1. Interlayer with $K$ and $S_{s}$ Anomaly

To explore the effects of an aquitard interlayer on the tidal wave propagation, we first assume that the interlayer is 0.5 m thick and located at the middle of the aquitard. Consequently, the aquitard of interest can be divided into three zones (i.e., $N = 3$ ⁠) with $z_{1} = 4.75$ m, $z_{2} = 5.25$ m, and $z_{3} = 10$ m. Zones 1 and 3 are prescribed with the same hydraulic parameters that are $K_{1}^{'} = K_{3}^{'} = 1.0 \times 1 0^{- 2}$ m/d and $S_{s, 1}^{'} = S_{s, 3}^{'} = 1.0 \times 1 0^{- 4}$ m^-1. The interlayer (i.e., zone 2) is prescribed with dramatically different $K$ (case 1) or $S_{s}$ (case 2) values to illustrate local anomaly of aquitard hydraulic parameters.

For case 1, we consider aquitard $K$ anomaly, while the aquitard keeps vertical homogeneity in $S_{s}$ ⁠, i.e., $S_{s, 2}^{'} = 1.0 \times 1 0^{- 4}$ m^-1. Two $K$ values are, respectively, prescribed to the interlayer, which are $K_{2}^{'} = 1.0 \times 1 0^{- 0}$ m/d (a high- $K$ interlayer) and $K_{2}^{'} = 1.0 \times 1 0^{- 4}$ m/d (a low- $K$ interlayer). In addition, the situation associated with the homogeneity in aquitard $K$ ⁠, i.e., $K_{2}^{'} = 1.0 \times 1 0^{- 2}$ m/d, is also considered for the aim of comparison. A high- $K$ interlayer is generally corresponding to a thin subunit consisting of higher-content sands or lower-content clay in the aquitard. Comparatively, a low- $K$ interlayer is generally corresponding to a thin subunit that consists of higher-content clay. Both of the two situations are frequently seen in practical hydrogeological investigations. The situation of inland constant head boundary was also considered to reveal the generality of the newly developed analytical solution and to demonstrate the difference in the tidal wave propagation between the common coastal leaky aquifer system and the island leaky aquifer system.

Figure 2 shows that a high- $K$ interlayer seems to produce negligible effects on tidal wave propagation from both the perspectives of amplitude and phase shift. However, a low- $K$ interlayer can produce a remarkable bias of tidal wave propagation. To be specific, as shown in Figure 2, the low- $K$ interlayer will significantly increase $α$ ⁠, as well as $τ$ (refer to corresponding absolute value hereinafter). This is mainly due to the fact that the low- $K$ interlayer blocks the groundwater leakage into the confined aquifer, and the tidal wave propagation tends to be increasingly dominated by the confined aquifer.

As for the coastal aquifer system, the pattern of head fluctuations in the confined aquifer is similar to that of the bilaterally fluctuating island aquifer system, but the amplitude variation is small. Since the inland is a constant head boundary, the fluctuation amplitude gradually decreases inland until it is 0.

For case 2, we consider the $S_{s}$ anomaly, while the aquitard keeps vertical homogeneity in $K$ (i.e., $K_{2}^{'} = 1.0 \times 1 0^{- 2}$ m/d). Two $S_{s}$ values are, respectively, prescribed to the interlayer, which are $S_{s, 2}^{'} = 1.0 \times 1 0^{- 3}$ m^-1 (a high- $S_{s}$ interlayer) and $S_{s, 2}^{'} = 1.0 \times 1 0^{- 5}$ m^-1 (a low- $S_{s}$ interlayer). In addition, the situation associated with the homogeneity in $S_{s}$ (i.e., $S_{s, 2}^{'} = 1.0 \times 1 0^{- 4}$ m^-1) is also considered for the aim of comparison. Generally speaking, the specific storage values should be also highly correlated to the clay content for an aquitard. In this study, a high- $S_{s}$ interlayer may correspond to the deposits with more soft clay particles showing higher compressibility, and vice versa.

Figure 3 shows that a low- $S_{s}$ interlayer seems to produce negligible effects on tidal wave propagation, especially the phase shift. However, a high- $S_{s}$ interlayer will cause a decrease in $α$ and $τ$ ⁠. The main reason for this phenomenon is that the higher the $S_{s}$ of the interlayer, the more water it can release or store. As a consequence, relatively more leakage into the confined aquifer is expected. In other words, the tidal wave propagation tends to be decreasingly dominated by the confined aquifer.

Now, we analyze the effects of the locations of the low- $K$ interlayer and the high- $S_{s}$ interlayer on the tidal wave propagation. Here, we consider the $K$ values of zones 1 and 3 to be $K_{1}^{'} = K_{3}^{'} = 1.0 \times 1 0^{- 1}$ m/d, while zone 2 is associated with $K_{2}^{'} = 1.0 \times 1 0^{- 3}$ m/d (a low- $K$ interlayer). For comparison analysis, we consider the interlayer locating in the upper aquitard associated with $z_{1} = 7.25$ m and $z_{2} = 7.75$ m and the interlayer locating in the lower aquitard associated with $z_{1} = 2.25$ m and $z_{2} = 2.75$ m.

Figure 4 demonstrates the effects of the low- $K$ interlayer location on the tidal wave propagation. Specifically, when the low- $K$ interlayer locates closer to the confined aquifer, the amplitude $α$ will be significantly increased, while the radian phase shift $τ$ will be also increased. This is mainly due to the fact that the closer of the low- $K$ interlayer to the confined aquifer, the weaker the hydraulic connection between the confined aquifer and the overlying unconfined aquifer and the smaller the amount of aquitard depletion, so that the vertical leakage from the aquitard is far more decreased, making the tidal effect tend be dominated by confined aquifer rather than be influenced by the leakage. On the contrary, when the low- $K$ interlayer locates in the upper aquitard, the amplitude $α$ will decrease and the radian phase shift $τ$ will increase. This is mainly because the low- $K$ interlayer that locates further to the confined aquifer can make the aquitard release a certain amount of groundwater from storage and leake to the confined aquifer to mitigate the tidal wave propagation.

Figure 5 demonstrates the effects of high- $S_{s}$ interlayer location on the tidal wave propagation. Specifically, when the high- $S_{s}$ interlayer locates closer to confined aquifer, the decrease in the amplitude $α$ is almost negligible and the radian phase shift $τ$ will be visibly increased. The cause for this phenomenon is that the high- $S_{s}$ interlayer locates closer to the confined aquifer, and the vertical leakage will be enhanced by more amount of groundwater released from aquitard storage, inhibiting the groundwater level fluctuations in the confined aquifer. However, for a relatively larger aquitard $K$ (i.e., $1.0 \times 10^{- 1}$ m/d considered here), suggesting stronger hydraulic connection between the confined aquifer and the overlying unconfined aquifer, the amount of groundwater released from aquitard storage is insignificant compared to the vertical leakage from the unconfined aquifer, so the amplitude $α$ is almost unchanged with the location of high- $S_{s}$ interlayer. However, we note that when the aquitard $K$ is relatively small (i.e., $1.0 \times 10^{- 2}$ m/d considered before), the presence of a high- $S_{s}$ interlayer decreases the radian phase shift $τ$ (i.e., Figure 3(b)), while when the aquitard $K$ is relatively large (i.e., $1.0 \times 10^{- 1}$ m/d considered here), a high- $S_{s}$ interlayer increases the radian phase shift $τ$ (i.e., Figure 5(b)).

Overall, either a low- $K$ or a high- $S_{s}$ interlayer in the aquitard can result in significant bias of tidal wave propagation in the confined aquifer, implying that the two aquitard heterogeneity scenarios cannot be ignored and this is particularly true when the low- $K$ or high- $S_{s}$ interlayer locates closer to the confined aquifer.

3.2. Exponentially Decaying $K$ and $S_{s}$

To consider the scenario of vertical aquitard heterogeneity in $K$ and $S_{s}$ that can be empirically expressed as an exponential function of depth (e.g., [26, 40, 41]). In this study, the function $f (z) = f_{0} \exp (λ z)$ is applied to express the exponential variation pattern, where $f_{0}$ is the corresponding value at the contact interface between the aquitard and confined aquifer and $λ$ is the decay exponent expressing the exponential relationship. Note that a positive $λ$ produces a depth-decaying exponential relationship of $K$ and $S_{s}$ versus depth. Particularly, $λ = 0$ refers to the homogeneous situation.

It should be noted that the exact analytical solution with the incorporation of exponential relationship mentioned above is beyond the scope of this study. To serve our analysis, we seek to uniformly divide the entire soil depth profile of the aquitard into enough zones, each of which is assumed to be homogeneous in $K$ and $S_{s}$ ⁠. After that, specific values are assigned to corresponding zones to approximate the exponential relationship. For the 10 m thick aquitard of interest, one can see from Figure 6 that a uniform discretization of 20 zones is accurate enough to capture the exponential relationship of $K$ and $S_{s}$ versus depth. To consider the situations of slightly, moderately, and highly heterogeneous aquitard $K$ and $S_{s}$ ⁠, we set $λ = 0.2$ m^-1, 0.3 m^-1, and 0.4 m^-1, respectively. Meanwhile, the homogeneous situation, i.e., $λ = 0.0$ m^-1, is also included for the aim of comparison.

Figure 7 shows that the tidal wave propagation in the confined aquifer is greatly affected by the exponentially decaying aquitard $K$ ⁠. In particular, a larger $λ$ results in a smaller $α$ ⁠, as well as $τ$ ⁠. Such a phenomenon is easy to understand. A larger $λ$ means the groundwater fluctuations in the confined aquifer can penetrate to a larger distance from the aquitard bottom. Consequently, the aquitard can supply more leakage for the confined aquifer to mitigate periodic head variations.

Similarly, in terms of the exponentially decaying $S_{s}$ ⁠, a larger $λ$ means the aquitard can release or store more water provided the same groundwater level fluctuations inside. Thus, more leakage from the aquitard can be supplied for the confined aquifer and smaller $α$ and $τ$ will be produced (Figure 8).

Overall, we can see from Figures 2–8 that the vertical heterogeneity in aquitard $K$ and $S_{s}$ has significant effects on the tidal wave propagation in the coastal leaky aquifer system, implying the importance of practical hydrogeological investigations of vertical aquitard heterogeneity. In addition, the newly proposed analytical model can be applied to determine hydraulic parameters of a leaky aquifer system via interpreting in situ hydraulic response data of surface water-groundwater interaction.

4. Conclusions

In this study, we derive an analytical solution for the tidal wave propagation in a coastal leaky aquifer system with the consideration of vertical heterogeneity in aquitard $K$ and $S_{s}$ ⁠. Two scenarios of vertical aquitard heterogeneity in $K$ and $S_{s}$ (interlayer with $K$ and $S_{s}$ anomaly and exponentially decaying $K$ and $S_{s}$ ⁠) are considered, and their effects on the tidal wave propagation are analyzed. Main conclusions are summarized as follows:

(1)
A low- $K$ aquitard interlayer makes the amplitude and the radian phase shift of the confined aquifer increase obviously, while the high- $K$ aquitard interlayer has negligible effects on the tidal wave propagation. A low- $S_{s}$ aquitard interlayer results in negligible effects on the tidal wave propagation, especially the phase shift. However, a high- $S_{s}$ aquitard interlayer will cause a decrease in the amplitude and phase shift
(2)
If the low- $K$ aquitard interlayer locates closer to the confined aquifer, it could significantly increase the amplitude while decrease the phase shift. On the contrary, if the low- $K$ aquitard interlayer locates further to confined aquifer, the amplitude decreases while the phase shift increases. When the high- $S_{s}$ aquitard interlayer locates closer to the confined aquifer, the amplitude decreases and the phase shift increases. On the contrary, if the high- $S_{s}$ aquitard interlayer locates further to the confined aquifer, the amplitude increases while the phase shift decreases
(3)
Exponentially decaying aquitard $K$ and $S_{s}$ is expected and has significant effects on the tidal wave propagation. A large decay exponent, corresponding to the highly heterogeneous aquitard $K$ and $S_{s}$ ⁠, yields a remarkable decrease in the amplitude and the phase shift. Practical hydrogeological investigations of vertical aquitard heterogeneity are highly recommended

Data Availability

The source code data of the analytical solutions used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declared that they have no conflicts of interest to this work. We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.

Acknowledgments

This research was partially supported by the Key Research and Development Program of National Natural Science Foundation of China (2019YFC1804301), the National Natural Science Foundation of China (41902244), the Basic Research Program of Jiangsu Province (Natural Science Foundation) (BK20190497), and the China Postdoctoral Science Foundation (Grant No. 2021M690866).

Exclusive Licensee GeoScienceWorld. Distributed under a Creative Commons Attribution License (CC BY 4.0).