Unsaturated Flow Influences the Response of Leaky Aquifer to Earth Tides

It is important for us to continuously monitor the confinement of aquifers as the confinement may be transient and can be destroyed by various disturbances, such as earthquakes, which will cause the wastewater leakage and pollute the environment [1, 2]. The response of water level in wells to Earth tides has recently become popular to evaluate the permeability of leaky aquifers [3–5]. Yet, these studies are flawed as they did not consider the effect of unsaturated zone, which might change the water movement in the saturated zone and thereby induce watertable fluctuations [6–8]. The unsaturated zone can store or release water, leading to the watertable fluctuations. When the unsaturated zone is excluded, the estimated permeability from the tidal response of groundwater could be essentially wrong [9, 10].

There have been studies revolving around how the unsaturated flow influences the response of coastal aquifers to the ocean tides [11–15]. According to experiments of Cartwright, in a shallow aquifer, the unsaturated zone may be potentially truncated by the ground surface and hence change the overall water storage in the unsaturated zone [16, 17]. Currently, considerable attentions have been paid to how the unsaturated flow influences the response of the ocean tides. Nevertheless, research regarding how the unsaturated flow influences the response of groundwater to the Earth tides is quite rare. By studying the water level in a deep well in China, Liao and Wang [18] identified that the amplitude of the tidal response decreases sharply, while its phase shift increases saliently when the watertable nears the ground surface. Wang et al. [9] made some preliminary explanations for the above phenomena by considering the capillary effect on groundwater response, based on one-dimensional Richards equation for describing watertable fluctuations in an unconfined aquifer. Actually, deep wells cannot be deemed as purely unconfined aquifers, so a multilayer model with the unsaturated topmost unconfined aquifer, which is much closer to the actual situation, needs to be established. Considering that the fundamentals of the ocean tides and the earth tides differ vastly, accordingly, their mathematical formulation differs markedly. To interpret the response to Earth tides in a multilayer model, a two-dimensional Richards equation coupling poroelastic deformation of the aquifer to groundwater flow is essential.

Here, we establish a two-dimensional multilayered finite element model to investigate the influence of unsaturated flow on the leaky aquifer to the Earth tides. The numerical model is verified by an analytical solution without considering the unsaturated flow [2] and a one-dimensional numerical model considering the unsaturated flow [9]. We then explore how the unsaturated flow affects the tidal response of a leaky aquifer. Finally, we illuminate the tidal response of a well installed in Lijiang, Yunnan province, China, by employing the numerical model.

We employ COMSOL Multiphysics to solve the differential equation (1), constrained by the boundary conditions mentioned above and by the material properties to be set later. Small element size (0.1 (m)) and time steps (500 (s)) are required to guarantee accuracy and numerical stability.

We make comparisons between the simulated results for a leaky aquifer and the analytical solution in Wang et al. [2]. It is derived from the classical leaky aquifer model of Hantush and Jacob [20], which assumes a relatively thin aquifer confined above by a semiconfining aquitard with no storativity. Leakage is supposed to occur by vertical flow across the semiconfining aquitard and is included as a volumetric source $Kz2hp/b2$⁠.

To verify this 2D simulated model with the 1D analytical leaky model of Wang et al. [2], first, we set the watertable is the ground surface, so the whole model is saturated. Then, we choose a thickness of 10 (m) for the aquifer and 100 (m) for the aquitard, and we select a range of horizontal transmissivity and storativity for the aquifer and zero storativity for the aquitard as in the analytical model. Figure 2 shows that the simulated results (circles) agree with the analytical solution (curves) for the studied ranges of transmissivity and storativity.

We next compute the tidal response of the one-dimensional unconfined silt aquifer and the sand aquifer (Table 1). Here, $Tzi$ and $Ssi$ are assigned the same value in all layers as $T$ and $S$⁠, respectively, to simulate a half-space unconfined aquifer, and the van Genuchten parameters and hydrology parameters of sand and silt are shown in Table 1. The simulated phase shift and amplitude ratio are plotted against $z/δ$ (Figures 3(a) and 3(b)), where $δ≡2D/ω$⁠, z is the depth of the screening interval of a cased well, $D=T/S$ is the hydraulic diffusivity, and $ω$ is the angular frequency of the $M2$ tide. By setting the average watertable at different depths, we obtain the set of curves in Figure 3. The figure shows that, with decreasing screened interval depth $z$⁠, the amplitude of the oscillation decreases and the phase increases.

Compared to the results by Wang et al. [9], which is a one-dimensional unconfined aquifer with saturated and unsaturated zones. When the watertable is at the ground surface (i.e., the whole aquifer is saturated), both simulations can be simplified to Doan’s solution. However, when the watertable is below the ground surface (i.e., the topmost of the aquifer is unsaturated), there are some obvious differences between these two simulated results, especially the sand aquifer. The main reason is that the mixed boundary condition at the ground surface is different. In Wang et al.’s solution [9], the average watertable is set at the ground surface. Yet here, the average watertable is fixed to the initial given watertable. The ground surface boundary in this paper is much closer to the reality, so the results in this paper are much credible.

We may thus use the code to explore the effect of the unsaturated zone on the response of the leaky aquifer to the Earth tides.

Here, we set a thickness of 10 (m) for the aquifer and 100 (m) for the aquitard. To compare the result with the analytical solution of Wang et al. [2], the vertical permeability of the aquifer and the aquitard is the same which ranges from $10−16m2$ to $10−8m2$⁠. The specific storage of the aquitard is 0. The topmost unconfined aquifer is assumed to have high hydraulic conductivity 0.1 (m/s) [3], and we set the permeability of the vertical and horizontal of the saturated zone $10−8m2$ and the specific storage $10−51/m$ [3].

To investigate the effect of the thickness of the unsaturated zone on aquifer’s response to the Earth tides, firstly, a fine-grained topmost unconfined aquifer with thick unsaturated zone was calculated. Experimental data for $θ$⁠, $hp$⁠, and $k3θ$ have been accumulated and documented, and several empirical relations have been developed to accurately describe $C3θ$ and $k3θ$ as functions of $θ$ from the measured water retention data [25, 27]. For the present study, we use the empirical van Genuchten model. The residual water content $θr$ and the saturated water content $θs$ are 0.034 and 0.46, respectively. The fitting parameters $α$⁠, $l$⁠, and $n$ are for matching the model to experimental data which are 1.6, 0.5, and 1.37, respectively. For comparison, we also compute a coarse-grained unconfined aquifer with a thin unsaturated zone. The fitting parameters $α$⁠, $l$⁠, and $n$ are for matching the model to experimental data which are 14.5, 0.5 and 2.68, respectively. Here, the horizontal saturated permeability and specific storage for the aquifer are fixed as $10−12m2$ and $10−51/m$⁠. According to equation (1), the hydraulic properties of the topmost unconfined aquifer in the saturated and unsaturated zone are calculated, and the figures of both conditions described above are shown in Figures 4(a) and 4(b), respectively. It shows that the unsaturated zone of a fine-grained topmost unconfined aquifer, in this case, is largely truncated by the surface, while that of a coarse-grained topmost unconfined aquifer is relatively unaffected, the difference directly related to the different thicknesses of the unsaturated zone in the two aquifers.

From Figure 5, we see that when the watertable is at the surface, the unsaturated zone is removed and the simulated results match the classical model perfectly; but when the watertable is below a threshold depth beneath the surface, which depends on the thickness of topmost unsaturated zone, the presence of the unsaturated flow substantially alters the tidal response of the aquifer; the response and the tidal forcing become more closely coupled, and the simulated amplitude ratio approaches 1, and the simulated phase shift approaches 0. The unsaturated zone in the topmost aquifer acts as a seal for the leaky aquifer. The thicker the unsaturated zone, the better the seal effect.

Comparing the tidal response of the fine-grained and coarse-grained upmost unsaturated zone (shown in Figures 5(a) and 5(b)), we see that the fine-grained topmost aquifer has a thicker unsaturated zone (5 m), and the coarse-grained topmost aquifer has a thinner unsaturated zone (0.3 m). These phenomena are resulted from the different hydraulic properties of the topmost aquifer (shown in Figure 4), and the capillary effect in fine-grained unsaturated zone is stronger than coarse-grained unsaturated zone.

In the following discussion, we explore the sensitivity of the tidal response of a semiconfined aquifer to the VM parameters for the unsaturated zone. Only a preliminary exploration is provided; exhaustive exploration is left to future studies. We first fix $n$ at 1.6 and $Kz3$ at 0.1 (m/s) but vary the parameter $α$ with watertable at a depth of 2.5 (m) below the ground surface. Figure 6(a) shows that, when$α≤0.5$⁠, the simulated phase shift and amplitude ratio are close to those of the analytical solutions. But as $α$ increases to 0.9 and beyond, the simulated phase shift and amplitude ratio begin to deviate drastically from the analytical solutions, suggesting that the influence of the unsaturated flow is evident.

We next fix $α$ at 0.3 and $Kz3$ at 0.1 (m/s) but change the values of $n$ with a watertable depth of 2.5 (m) below the ground surface. Figure 6(b) shows that, when$n≥1.1$⁠, the simulated phase shift and amplitude ratio are close to the analytical solutions; but when$n≤1.1$⁠, the simulated phase shift and amplitude ratio deviate noticeably from the analytical solutions.

Finally, we fix the VM parameters $α$ at 0.3 and $n$ at 1.6 but vary the saturated vertical conductivity of the topmost unconfined aquifer (⁠$Kz3$⁠) from 0.0005 to 0.1 (m/s). Figure 6(c) shows that, when $Kz3≥0.02$ (m/s), the phase shift and amplitude ratio are similar to the analytical solution and the effect of the unsaturated flow may be neglected; but when $Kz3$ decreases to $Kz3≤0.01$ (m/s), the difference between the simulated results and the analytical solution becomes obvious. This result clearly shows that the unsaturated flow causes a decrease in the vertical conductivity, which may be compensated by increasing the saturated vertical conductivity $Kz3$⁠.

Barry et al. [10] introduced a dimensionless value $ω∗=ωB/K$ to correct the capillary effect in an unconfined aquifer, where $ω$ is the frequency of the earth tide, $B$ is the capillary parameter, and $K$ is the hydraulic conductivity. It induced that when $ω∗≪1$⁠, the capillary effect can be neglected. As for the leaky aquifer, from Sections 3.1 and 3.2, we can see that many factors can influence the tidal response of the water level in a well, for example, the watertable, the VG parameters, and topmost vertical conductivity. We also do a series of computations of different aquitard thicknesses, shown in supplementary material (available here), and find that the effect of the unsaturated flow is not influenced by the thickness of the aquitard. Hence, like Barry et al. [10], we also introduce a dimensionless quality $ω∗=ωhc/Kz3k3θ$ as a criterion to evaluate the effect of the unsaturated flow in a leaky aquifer, where $hc$ is the unsaturated zone height, $Kz3$ is the topmost vertical hydraulic conductivity, and the $k3θ$ is the relative hydraulic conductivity.

Considering Figure 6, we find that when$α≤0.5$⁠, $n≥1.1$⁠, or $Kz3≥0.01$ (m/s), the phase shift and amplitude ratio are close to the analytical solution. In other words, in these situations, the influence of unsaturated flow on the tidal response is none. We do not have to consider the unsaturated flow. However, in some situations just like$α>0.5$⁠,$n<1.1$⁠, or$Kz3<0.01$ (m/s), the tidal response of the groundwater differs from the analytical solution [2]. We interpret that the phenomenon above depends on the dimensionless quality $ω∗$ at the surface. Figures 7(a)–7(c) show the soil water retention curves of the above three situations, respectively. They show that, at a fixed watertable, the $ω∗$ decreases as $α$ and $Kz3$ increase and increases as $n$ increases. The watertable of the upmost unconfined aquifer is 2.5 m below the surface, and $ω∗$ of the surface is shown as circles and asterisks in Figure 7. The common character of three situations is that, when $ω∗<10−0.5$⁠, the unsaturated flow has little effect on the tidal response. However, when $ω∗>10−0.5$⁠, the effect of the unsaturated flow act on the tidal response becomes evident. The bigger the $ω∗$ is, the larger the effect would be.

The Lijiang well locates in Yunnan province, southwest of China. Liao and Wang [18] analyzed several years of the tidal response of an aquifer in SW China. The hydrogeology of this well consists of an aquifer of Triassic carbonate rocks overlain by an aquitard of Tertiary siltstone that, in turn, is overlain by a thin layer of Quaternary sediments. They found that, when the watertable is more than 2 (m) below the surface (dry seasons), the amplitude ratio is ~1 and the phase shift to the $M2$ tide is ~0° (Figure 8). The small difference between the phase shifts from zero is likely to be within observational errors because, due to a lack of strain measurement near the Lijiang well, the theoretical tidal strain was used as the reference [18], which may differ by a few degrees from the actual strain [28]. When the watertable rises close to the surface, on the other hand, the phase shift increases from ~0° to ~20° and the amplitude ratio decreases from ~1 to ~0.7.

Previous interpretations of this change in tidal response were based on a model for unconfined aquifer [18]. But the lithological log of this well (Figure 9(a)) clearly shows that the aquifer is not unconfined but is overlain by a thick (167 (m)) aquitard. Hence, the previous interpretations may be inappropriate to explain the phenomenon above. Here, we establish a two-dimensional multilayer model according to the Lijiang well completion diagram, and the model has three layers (Figure 9(b)). The depth of the whole model is 347.3 (m), and the thickness of each layer from bottom to top are, respectively, 179.8 (m), 103.4 (m), and 64.1 (m) (Figure 9(a)). The topmost aquifer is siltstone with a certain height of unsaturated zone. The permeability and storativity of siltstone are $10−10$ (m²) and $2×10−6$ (1/m) [29]. After a series of simulations, we infer that the permeability of the bottom and middle aquifers, such as the horizontal and vertical permeability of the bottom aquifer are, respectively, $5×10−14$ (m²) and $10−13$ (m²). The storativity of the bottom aquifer is $4.5×10−7$ (1/m). The permeability of the middle aquifer is $10−13$ (m²). The residual water content $θr$ and the saturated water content $θs$ are 0.034 and 0.46, respectively. The fitting parameters $α$⁠, $l$⁠, and $n$ are for matching the model to experimental data which are 0.15, 0.5 and 1.37, respectively.

The simulated results are shown in Figure 8, in which the simulated phase shift and amplitude ratio agree well with the observed results, which illustrates that the unsaturated zone above the watertable has an obvious effect on the tidal response of the well water level. It caused the phase shift and amplitude ratio, respectively, decreasing and increasing sharply when the watertable is between the ground surface and 2 (m), and the phase shift and amplitude ratio remain almost the same when the watertable is below 2 (m). According to the above simulation, the unsaturated zone above the watertable maybe is a reason which leads to the phase shift increasing and amplitude ratio decreasing when the watertable approaches 2 m below the ground surface.

In this article, we prove that the unsaturated flow can largely influence the tidal response of leaky aquifers; the assessment of aquifer leakage from the observed tidal response cannot be ignored. Particularly, we have expounded on the effects of unsaturated flow on the groundwater table fluctuations in these aquifers vulnerable to periodical tidal forcing. The results revealed that the thickness of the unsaturated zone basically determines the effects of tidal response. Compared to that estimated by the traditional leaky model, the phase shift turned smaller, and the amplitude ratio is bigger. Additionally, a dimensionless quality $ω∗$ is introduced to estimate the effect of the unsaturated flow; when $ω∗<10−0.5$⁠, the effect of the unsaturated flow can be neglected, otherwise, the capillary fringe has great effect on the tidal response of the water level; in other words, the estimated result of aquifer leakage differs vastly when the unsaturated flow is not considered when $ω∗>10−0.5$⁠. Using this two-dimensional Richards model, the tidal response of the Lijiang well can be explained perfectly, which reveals that the unsaturated flow may be an important factor which influences the water level of the well. Besides, for the sake of obtaining much more precise permeability of the aquifer from the tidal response of the well, the multilayer two-dimensional Richards model is much better than the one-dimensional Richards model. Our study suggests that the shallow and deep hydrological processes may interact across thick layers of intervening aquitards and that the analysis of the tidal response may be an effective means to evaluate such interaction.

The data used in this paper are all from references.

The author declares that there are no conflicts of interest.

This work was supported by the National Key R&D Program of China (No. 2018YFC1503200 and No. 2018YFC0603502) and the Special Fund of the Institute of Geophysics, China Earthquake Administration (No. DQJB19A0124). We thank Prof. Chi-Yuen Wang’s and Yan Zhang’s helpful advice and suggestion.