Impact of Depth‐Dependent Heterogeneity in Aquifers on the Response of Unsaturated‐Saturated Flow to Earth Tides

Heterogeneity is a natural characteristic of aquifers, which has a profound influence on groundwater flow and solute transport. However, the heterogeneity in an aquifer is difficult to characterize and its impact on groundwater response to Earth tides remains insufficiently explored. A common heterogeneity in thick aquifers is the decrease in porosity and hydraulic conductivity with depth. In the present study, we present an analytical model that includes this particular heterogeneity in both the unsaturated and saturated zones to explore its impact on the tidal response. Our results reveal that, as the hydraulic conductivity decays with depth, the amplitude ratio of the tidal response increases and its phase shift decreases, and the impact of the capillary zone on tidal responses is primarily restricted to the shallow region of the aquifer. Neglecting the decay of conductivity in the previous models would result in underestimation of the amplitude ratio and overestimation of the phase shift. The solutions are applied to the field data to effectively estimate the decay in conductivity with depth. This study highlights the importance of considering aquifer heterogeneity in the analysis of the response of saturated and unsaturated flow to Earth tides.


Introduction
Gravitational interactions between the Earth, the Moon, and the Sun cause periodic deformation of the Earth's interior (Earth tides) and cause the groundwater level to fluctuate, making it a crucial element in the study of the dynamics of groundwater.Understanding the amplitude and phase of these groundwater fluctuations is of utmost importance, as they are intricately linked to aquifer parameters, with hydraulic conductivity playing a key role.As such, numerous studies have focused on estimating aquifer parameters by examining the correlation between groundwater fluctuations and Earth tides (Bredehoeft, 1967;Cutillo & Bredehoeft, 2011;Lai et al., 2013;Lu et al., 2022;Sun & Xiang, 2019;Vanderkamp & Gale, 1983;Xue et al., 2016).Over the years, various analytical models have been developed to investigate Earth tidal phenomena in aquifers (Detournay & Cheng, 1993;Hsieh et al., 1987;Roeloffs, 1996;H. Wang, 2000; C.-Y. Wang et al., 2018;G. L. Wang et al., 2023).While these models have provided valuable insights, the effect of the unsaturated zone has often been overlooked.The unsaturated zone serves as a critical conduit between the surface and subsurface environments, playing a crucial role in groundwater dynamics and substantially influencing aquifer responses to Earth tides (C.-Y.Wang et al., 2019).
Recent studies conducted by C.-Y. Wang et al. (2019), Liang et al. (2022), and Zhang et al. (2023) have highlighted the significant influence of the unsaturated zone on the Earth tide response of groundwater.C.-Y. Wang et al. (2019) employed a numerical approach that integrated groundwater flow with poroelastic strain in both the saturated and unsaturated zones.Their findings revealed that the presence of the unsaturated zone, when the water table drops below the surface, leads to deviations in predicted amplitude ratios and phase shifts compared to the traditional models.To better understand the influence of the unsaturated zone on tidal response and parameter estimation of aquifers, Liang et al. (2022) developed an novel analytical model that specifically accounts for the effect of the capillary zone in a half-space unconfined aquifer.This model effectively explains the observed seasonal variations in phase shift and amplitude ratio.Additionally, Liang et al. (2022) demonstrated that traditional models used for unconfined aquifers, which neglect the presence of the unsaturated zone, tend to overestimate the phase shift and underestimate the amplitude ratio.Extending this line of research, Zhang et al. (2023) expanded the analytical model by incorporating additional factors like finite aquifer thickness, anisotropy, and wellbore storage, further underscoring the importance of the unsaturated zone and its interplay with aquifer characteristics.
Despite these significant advances, a critical aspect remains unexplored: the impact of aquifer heterogeneity.Heterogeneity is an inherent property of aquifers, posing considerable challenges in accurately predicting groundwater flow and contaminant transport in subsurface environments (Freeze & Cherry, 1979;Yeh et al., 2015).Unfortunately, characterizing aquifer heterogeneity remains a daunting task, and its influence on groundwater responses to Earth tides remains insufficiently explored.
In light of this knowledge gap, our study seeks to extend the analytical model proposed by Liang et al. (2022) and examine the effects of heterogeneity on tidal analyses.In this paper, we present a mathematical model with a depth-dependent aquifer heterogeneity for unsaturated-saturated flow driven by Earth tides.By integrating this heterogeneity, our research extends our understanding on the intricate dynamics between Earth tides and water flow in aquifers.We first validate our model using numerical simulations and assess its performance in heterogeneous aquifers.We then apply the model to both synthetic groundwater data and to data from a field site.Through this effort, we aim to unravel the influence of the depth-dependent heterogeneity on the tidal response of groundwater and the utility of the tidal method toward a more comprehensive understanding of groundwater response to Earth tides.

Heterogeneity of Hydraulic Conductivity With Depth
The porosity of sedimentary rocks is known to decrease with depth.Based on the analysis of numerous sedimentary cores, Athy (1930) proposed an exponential relation to describe the decrease of porosity with depth, often referred to as the Athy's law: where φ denotes porosity [ ], φ o is the porosity at the surface [ ], z is depth [L], and c is a material constant [L 1 ].For example, Figure 1a shows the changes in porosity with depth for sand (c = 4.20 × 10 4 m 1 ) and shale (c = 5.53 × 10 4 m 1 ) (Proshlyakov, 1960) and their fits to the Athy's relation.(Proshlyakov, 1960).(b) Estimating Depth-related hydraulic conductivity for sand and shale using Athy's relation and the Kozeny-Carman equation.

Water Resources Research
10.1029/2023WR036310 The interconnection between porosity and hydraulic conductivity is crucial, as the decay of porosity with depth inherently leads to a corresponding decrease in hydraulic conductivity.However, the extent of this effect can be influenced by various hydrogeological conditions.In practice, the Kozeny-Carman equation, a semi-empirical formulation, is often utilized to establish the relationship between permeability and porosity (Carman, 1937(Carman, , 1956;;Kozeny, 1927).Furthermore, Bourbié et al. (1987) proposed a modified version of the Kozeny-Carman equation, incorporating a variable power on porosity, thereby removing the limitations of the equation and enhancing the accuracy of the permeability estimation: where k p [L 2 ] is permeability; C is an empirical coefficient [ ]; d is the average diameter of solid grains [L]; φis the porosity; n is an empirical exponent [ ], with a value between 3 and 8 (Xu & Yu, 2008).Assuming that the average diameter of solid grains does not change significantly with depth, we may combine the above two relations to obtain where Equation 3 demonstrates that permeability (and equivalently, hydraulic conductivity K s (φ)) decays exponentially with depth with a decay coefficient c s = nc.This exponential decay of hydraulic conductivity with depth is illustrated in Figure 1b, where the normalized hydraulic conductivity (k p /k 0 ) displays distinct rates of decrease in different lithologies, with n set to 3.
Based on Equation 3, we describe the porosity-dependent component of the hydraulic conductivity K φ (φ) for the unsaturated and the saturated zones using the following exponential model, where K 0 is the hydraulic conductivity near the interface between the unsaturated and saturated zone (z = 0) [LT 1 ]; c u is the decay coefficient of the K φ (it usually refer to the saturated hydraulic conductivity) in the unsaturated zone (0 ≤ z ≤ b) [L 1 ]; c s is the decay coefficient of the conductivity due to porosity change in the saturated zone ( L ≤ z ≤ 0) [L 1 ], which could be differ from that of the unsaturated zone.In the unsaturated zone, the hydraulic conductivity is also a function of the degree of saturation, as discussed in the next section.
The depth-dependence of aquifer conductivity in the form of Equation 4 is commonly used in previous studies, as compiled in Table 1.The table shows that the c s values range from 0.001 to 0.03 m 1 , which is consider as a reference in the later evaluation of our mathematical model.

Mathematical Model
Figure 2 illustrates the schematic diagram of groundwater flow in an unsaturated-saturated zone induced by Earth tides.The vertical z-axis represents upward direction, simplifying the model to a laterally infinite extent with only one-dimensional vertical flow.The saturated zone represents an unconfined aquifer with a finite thickness L and infinite lateral extent.The vertical hydraulic conductivity of the saturated zone varies with depth, and the saturated conductivity of the unsaturated zone is constant.The unsaturated zone is characterized by an average thickness of b, and the average water table is located at the origin (z = 0).
The undrained response of hydraulic heads in the saturated zone to volume strain induced by the diurnal and semidiurnal Earth tides can be described by the following governing equation (Liang et al., 2022;H. Wang, 2000;C.-Y. Wang et al., 2019) where h is the hydraulic head in the saturated zone [L]; K φ (z) is the depth-dependent hydraulic conductivity [LT 1 ], which is described by Equation 4; S s is the specific storage [L 1 ]; B is the Skempton's coefficient of the aquifer [ ]; K u is the undrained bulk modulus of the aquifer [Pa]; ρ is the density of water [ML 3 ]; g is the acceleration of gravity [LT 2 ]; ε is the poroelastic volumetric strain produced by Earth tides [ ]; ξ is the instantaneous depth of the moving water table.
The undrained response of hydraulic heads in the unsaturated zone to strain forced by the Earth tides can be described by the following Richards equation (Liang et al., 2022;C.-Y. Wang et al., 2019) a M represents the measured data, and S represents the synthetic data.b The c s is calculated using the Athy's Law and the Kozeny-Carman Equation with the depth-related porosity (Proshlyakov, 1960).It should be noted that the periodic tidal deformation occurring in the Earth's solids causes a change of water content in the pore volume of the aquifer, which is accounted for in our tidal model by the source terms in Equations 5a and 6a.This change of porosity is extremely small (on the order of 10 8 ) and hence is elastic and oscillatory (H.Wang, 2000); it does not contribute to a permanent change in porosity or its associated hydraulic conductivity.

Water Resources Research
Equation 6b is the Robin boundary condition, which has a capacity to adjust to the dynamic changes of boundary conditions.It will reduce to a no-flow boundary when the water table depth is large and to a constant head boundary (u = 0) when the water table rises to the ground surface.Detailed explanations of this boundary condition were given in earlier studies (Chui & Freyberg, 2009;Liang et al., 2022;Shoushtari et al., 2015;C.-Y. Wang et al., 2019) and is also discussed in Supporting Information S1.
In this study, the unsaturated medium properties are described by the Gardner-Kozeny model (Gardner, 1958;Kroszynski & Dagan, 1975;Mathias & Butler, 2006).Liang et al. (2022) compared this model with the van Genuchten-Mualem model (van Genuchten, 1980) and showed that both models are equivalent and can be interconverted.where θ r is the residual water content [ ]; θ s is the saturated water content [ ]; a c and a k are the constitutive exponents of the models of water retention and relative hydraulic conductivity, respectively [L 1 ]; S y = θ s θ r , is the drainable porosity or the specific yield [ ]; ψ a is the pressure head threshold at which desaturation of an aquifer begins and denotes the minimum pressure head required for air infiltration in a saturated medium [L]; ψ k is the pressure head above which relative hydraulic conductivity is effectively equal to unity, which differ to the air entry pressure head ψ a .The Gardner-Kozeny model is physically sound and amenable to analytical treatment.Therefore, it has been widely adopted in the analytical solutions of the Richards equations (Chang et al., 2018;Liang et al., 2017aLiang et al., , 2017bLiang et al., , 2018aLiang et al., , 2018b;;Lin et al., 2017;Malama, 2014;Mathias & Butler, 2006;Mishra & Neuman, 2010;Qi et al., 2020;Tartakovsky & Neuman, 2007;Zhang et al., 2023;Zhou et al., 2023).
The unsaturated and saturated zone flow regimes are coupled by interface conditions representing continuity of pressure and normal flux across the water table, which take the following forms: The coupled Equations 5-8 present challenges in analytical solutions due to the nonlinear nature of the Richards equation and the unknown position of the moving water table ξ.However, it is worth noting that the changes in hydraulic heads caused by Earth tides are relatively small compared to most aquifer thickness.Therefore, we may employ the perturbation method to linearize the Richards equation and fix the water table at a known position (Kroszynski & Dagan, 1975), as done in previous studies in solving the coupled unsaturatedsaturated flow equations (Chang et al., 2018;Liang, Zhan, Zhang & Liu, 2017;Liang, Zhan, Zhang, & Schilling, 2017;Liang, Zhan, Zhang, et al., 2018;Liang, Zhan, & Zhang, 2018;Lin et al., 2017;Mathias & Butler, 2006;Mishra & Neuman, 2010;Tartakovsky & Neuman, 2007;Zhang et al., 2023).Moreover, Liang et al. (2022) conducted a thorough comparison between the analytical solution obtained through linearizing the Richards equation and the numerical solution of the nonlinear Richards equation.The study showed a well agreement between the linearized solution and the nonlinear solution for groundwater flow induced by Earth tides.Following the process of linearization in Liang et al. (2022), we express the governing Equation 6a to the following simplified form: where where k 0 (z) and C 0 (z) are the zero-order approximations of the relative conductivity k (θ) and the specific moisture capacity C(θ), respectively, and can be written in the following forms, It should be noted that if we neglect the porosity-dependence of the hydraulic conductivity in the unsaturated zones, the parameters α * k and γ* will reduce to a k and γ, respectively.
With the linearized treatment, the water table is fixed at the average water table depth, that is, ξ = 0. Therefore, the interface conditions of the unsaturated and saturated flow regimes reduce to the following forms

Analytical Solutions
We assumed that the oscillatory strain induced by earth tides has the form ε = ε 0 e iωt , where ε 0 is the complex amplitude, ω = 2π/τ is the angular frequency [T 1 ], τ is the period of the tidal oscillation [T], and i = ̅̅̅̅̅̅ 1 √ .Therefore, the stationary periodic solutions of Equations 5a and 9 have the following forms where h 0 (z) is the complex amplitude of h, u 0 (z) is the complex amplitude of u.Substituting the complex forms of h(z, t) and u(z, t) into Equations 5 and 9-11, respectively, yields the following solutions and where and J -1 and Y -1 are Bessel function of first and second kind of order 1, respectively; J n and Y n are Bessel function of first and second kind of order n = α * k /λ, respectively.The definitions of the other parameters C 1 ,C 2 ,C 3 , andC 4 in above equations and the details of derivation are given in Supporting Information S1 (Texts S2 and S3).
Based on Equation 13, the amplitude ratio of the tidal response of hydraulic head in the saturated zone can be written as (Rojstaczer & Riley, 1990; C.-Y. Wang et al., 2018;G. L. Wang et al., 2023) and the phase shift of the tidal response of hydraulic head in saturated zone can be written as where arg[x] is the argument of the complex number x.
It should be noted that Equations 14a and 14b are the amplitude ratio and phase shift of the hydraulic head at an observed point within the aquifer.However, the hydraulic head measured in the observation well reflects the average hydraulic head across the entire well screen.Therefore, in practical terms, the amplitude of the hydraulic head in an observation well should be expressed as: where z b z a is the length of the observation well screen [L].It should be noted that Equation 15 contains an integration that is challenging to evaluate analytically.Therefore we will evaluate it using the numerical integration methods.

Water Resources Research
10.1029/2023WR036310 WANG ET AL.
The amplitude ratio and phase shift of the tidal response of hydraulic head in an observed well can be written as For the homogeneous aquifer with a constant hydraulic conductivity K 0 , that is, c s = c u = 0, the solutions of h 0 (z) and u 0 (z), respectively, reduce to where and the definitions of the other parameters C′ 1 ,C′ 2 ,C′ 3 ,and C′ 4 in Equation 17and the details of their derivation are given in Supporting Information S1 (Text S4).The parameter δ usually denotes as the "skin depth" of the tidal responses of groundwater flow (Detournay & Cheng, 1993), which regulates the impacts of the aquifer thickness (Zhang et al., 2023).The corresponding amplitude ratio and phase shift of the tidal response of hydraulic head in a homogeneous aquifer are, respectively, and It should be noted that the solution Equation 18 was derived for a finite homogeneous aquifer, whereas the solution in Liang et al. (2022) was derived for a homogeneous half-space aquifer.We will discuss the differences between these solutions in the subsequent sections.

Comparison With Numerical Solution
In this section, we validate our analytical solution by comparing it with both linearized and nonlinear numerical solutions obtained using COMSOL Multiphysics, a Galerkin finite-element software package.The purpose of comparing it with the linearized numerical solution is to verify the accuracy of the analytical solution since both are solutions to the same governing Equations 5 and 9.The comparison with the nonlinear numerical solution of governing Equations 5 and 6, on the other hand, aims to assess the errors introduced by the linearization assumption.First, we compare our solution Equations 13a and 13b with the linearized numerical solution.For illustration, the response to the semi-diurnal tide (M 2 ) is used in the following diagrams.The model parameter values used in these tests are summarized in Table 2. Figure 3 illustrates the variations in the complex amplitude of hydraulic head with depth for different values of c s , c u , and α k .The solid curves represent the predictions from our analytical solutions, while the circle symbols represent the predictions obtained from the linear numerical solution of COMSOL.The comparison shows a good agreement between our analytical solutions and the numerical results.
Second, we compared the solution of the hydraulic heads (h = h 0 exp(iωt)) with the nonlinear numerical solution of governing Equations 5 and 6. Figure 4 shows a good agreement between the present solution (solid curves) and the numerical solution (circle symbols) for the fluctuations of the hydraulic heads at z = 600 m with time, computed with the GK model ( 7) and different values of c s , c u , and α k .The model parameter values used in these tests are the same as those in Figure 3 and summarized in Table 2.These comparisons confirm that the analytical solution proposed in this study can serve as a tool for tidal analysis in aquifers with depth-dependent heterogeneity in conductivity.

Effects of Depth-Dependent Saturated Conductivity in the Unsaturated Zone
The decay coefficient c u represents the decline in saturated conductivity (K φ ) in the unsaturated zone due to depthdependent changes in porosity, which may be different from the decay coefficient of the conductivity in the saturated zone (c s ).To discuss the effect of c u on the tidal response of the aquifer, we assume that the decay  b(m) 1 0.5, 1 and 2 1 1 0, 0.3, 0.6, 0.9, 1.2, and 1.5   coefficient c s remains constant at 0.005 m 1 , which allows us to isolate the impact of c u on the tidal behavior of the aquifer.Furthermore, we assume that the K φ at the ground surface (z = b) remains constant but decreases with depth from the same initial value under different decay coefficient c u .
Figure 5 depicts the changes of the amplitude ratio and phase shift with respect to normalized depth (z/L for the saturated zone and z/b for the unsaturated zone) for different decay coefficients c u (0, 0.03, and 0.3 m 1 ) in the unsaturated zone, where the parameters of model are summarized in Table 2.It is worth noting that the half-space solution presented by Liang et al. (2022) disregards the aquifer thickness L. In this context, employing z/L represents a dimensionless normalization of the z-coordinate, normalized by a specified value of L, to facilitate comparison with the present solution.It shows that when the thickness of the unsaturated zone is thin (b = 0.5 m), the curves both for the amplitude ratios and the phase shifts corresponding to different values of c u converge closely, indicating that the impact of c u on the tidal response in the unsaturated zone and saturated zone may be neglected.When b = 1 m, there is a discernible deviation in the amplitude ratio and the phase shift in the saturated zone for different c u suggesting that c u may exert a discernible impact on the tidal response of the aquifer.Specifically, a larger c u leads to an increased amplitude ratio and a decreased phase shift.However, as the thickness of the unsaturated zone continues to increase (b ≥ 2 m), the effect of c u will diminish (Figures 5e and 5f).This attenuation is attributed to the possibility that in thicker unsaturated zones, the capillary fringe is not truncated by the surface, a phenomenon documented by Zheng et al. (2022).As a result, the impact of the unsaturated zone is reduced, leading to the hydraulic heads fluctuate synchronously with Earth tides, characterized by an amplitude ratio of 1 and a phase shift of 0, as shown in Figures 5e and 5f.This is consistent with the previous research findings (Liang et al., 2022).

Effects of Depth-Dependent Hydraulic Conductivity in the Saturated Zone
In this section, we discuss the impact of conductivity change with depth in the saturated zone due to porosity change (different values of c s ) on the response of groundwater to Earth tides.To streamline the discussion, we assume that K φ in the unsaturated zone is constant (i.e., c u = 0 m 1 ) and focus our analysis on the impacts of conductivity decay in the saturated zone.As a result, the parameters α * k and γ* are reduced to a k and γ, respectively, throughout the following sections of this study.This simplifying assumption allows us to concentrate on the specific effects of depth-dependent conductivity in the saturated zone while acknowledging that the unsaturated zone's influence remains consistent with respect to these parameters.2, and with the half-space solution of Liang et al. (2022) (dark dashed curves).It shows that most of amplitude ratios generally increase and most of phase shifts decrease with depth, indicating a gradual weakening of the influence of the unsaturated zone as depth increases, consistent with the result in Liang et al. (2022).Additionally, some curves both for the amplitude ratios and the phase shifts in Figure 6 have protrusions and return to stability as the depth increases.The initial increase and then decrease with decreasing z/L of the amplitude ratio in the tidal response of groundwater in the saturated zone are due to the poroelastic coupling between pore pressure and deformation.The thickness of this anomalous zone scales with the parameter δ (Equation 17), which is known as the "skin depth" of the tidal response of groundwater in the saturated zone (Detournay & Cheng, 1993;H. Wang, 2000).At greater depths, the tidal oscillations of groundwater tend to synchronize with the Earth tide, with the amplitude ratio becoming 1 and the phase shift becoming 0. It is worth noting that the curves both for the amplitude ratios and the phase shifts of the unsaturated zone corresponding to different values of c s converge closely (z/b ≥ 0), indicating that the impact of c s on the tidal response of the unsaturated zone may be neglected.However, the impacts of the unsaturated zone are regulated by the vertical heterogeneity of aquifers.Figure 6 highlights the significant influence of the decay coefficient c s on both the amplitude ratio and phase shift.In a homogeneous aquifer (c s = 0 m 1 ), the amplitude ratio is relatively small and the phase shift is relatively large.As c s increases, both the amplitude ratio and phase shift progressively deviate from those of the homogeneous aquifer.When c s is very large (0.02 m 1 ), the amplitude ratio reaches 1 and the phase shift approaches 0 in the deep region of the aquifer (z/L ≤ 0.55), indicating that the hydraulic head fluctuates synchronously with the Earth's tides.This can be attributed to the faster decay rate of hydraulic conductivity with depth in the saturated zone as c s increases, leading to lower hydraulic conductivity (≤1.5 × 10 7 m/s) and correspondingly slower flow velocity at the deep region of the aquifer.Consequently, the influence from the unsaturated zone becomes less significant in this region.These observations suggest that the impacts of the unsaturated zone on the tidal responses of groundwater flow is limited in a shallow region of the aquifer due to the decrease in hydraulic conductivity with increasing depth.This results may offer a criterion to estimate the desirable depth for installing a pressure gauge in a well given the decay rate constant c s , or the impact of the unsaturated zone on the measured tides, given the depth of the pressure gauge.Furthermore, it is worth noting that the present solution for the homogeneous aquifer (c s = 0 m 1 ) deviates from the half-space solution proposed by Liang et al. (2022), particularly in the phase shift prediction.This discrepancy can be attributed to the consideration of a finite aquifer thickness in the present solution, while Liang et al.'s solution assumes an semi-infinite aquifer.In Figure 6, the "skin depth" (δ = 1,193 m) exceeds the aquifer thickness (L = 800 m), resulting in a noticeable deviation between the two solutions.This implies that the aquifer thickness plays a role in the tidal response.Zhang et al. (2023) interpreted this phenomenon by suggesting that when the "skin depth" exceeds the aquifer thickness, the discrepancy between the half-space model and the finitelayer model becomes more significant.However, as the aquifer thickness exceeds the three times of the value of "skin depth," this difference gradually diminishes.
Figure 7 depicts the variations of the amplitude ratio and phase shift with respect to normalized depth (z/L) for different decay coefficients c s in the saturated zone and different α k in the unsaturated zone, where the parameters of the model are summarized in Table 2. Figures 7a and 7b show the predictions based on the solution of Liang et al. (2022).The results highlight the significant influence of the constitutive exponent α k on both the amplitude ratio and phase shift.Smaller values of α k lead to lower amplitude ratios and larger phase shifts, consistent with previous findings (Liang et al., 2022;Zhang et al., 2023).However, the influence of α k on the amplitude ratios and phase shifts is regulated by the decay coefficient c s of the hydraulic conductivity.In the case of a homogeneous aquifer (c s = 0 m 1 ), the amplitude ratio generally increases with depth, and the phase shift decreases and approaches 0 at the bottom of the aquifer (z/L = 1).With an increase in c s , all amplitude ratios progressively tend toward 1 at depths greater than half the aquifer's thickness (z/L < 0.5), as illustrated in Figure 7g.Similarly, all phase shifts exhibit a tendency to converge progressively toward 0, as depicted in Figure 7h.This suggests that a faster decay of conductivity with depth leads to a more rapid increase in the amplitude ratio and a faster decrease in the phase shift.It implies that the influence of the capillary zone on tidal responses is primarily limited to the shallow region of the aquifer due to the decreasing hydraulic conductivity with depth.It is worth noting that the solution proposed by Liang et al. (2022), which neglects the impacts of conductivity decay with depth, may result in underestimation of the amplitude ratio and overestimation of the phase shift, particularly in the deep region of the aquifer.In other words, the influence of the capillary zone on the tidal responses of groundwater in the deep aquifer, when its conductivity decays with depth, might be overestimated in Liang et al.'s solution.
The water table depth or the vadose zone thickness (b) is another important parameter, distinct from the constitutive exponent α k .It is typically driven by external factors such as rainfall infiltration and undergoes temporal fluctuations.Therefore, understanding the influence of aquifer heterogeneity on the tidal response of groundwater at different b values is crucial.Figure 8 illustrates the variations in amplitude ratio and phase shift with respect to normalized depth (z/L) for different b values and different decay coefficients (c s ) in the saturated zone, with model parameters summarized in Table 2.The findings reveal that vadose zone thickness (b) significantly affects both the amplitude ratio and phase shift of the tidal response of groundwater.Smaller b values result in lower amplitude ratios and larger phase shifts, consistent with previous studies (Liang et al., 2022;Zhang et al., 2023).Additionally, b influences the impact of the decay coefficient (c s ) of hydraulic conductivity on the tidal responses of groundwater.For relatively small b values (b < 0.9 m), larger c s values lead to a more pronounced increase in the amplitude ratio and a faster decrease in the phase shift with depth.Conversely, for relatively large b values (b ≥ 1.5 m), the amplitude ratio approaches 1 and the phase shift approaches 0, indicating that the influence of aquifer heterogeneity becomes negligible.This highlights that the significance of aquifer heterogeneity is influenced by the vadose zone thickness, with thinner vadose zones having greater effects of aquifer heterogeneity on the tidal responses of groundwater.

Applications
To demonstrate the application capacity of the present solution, we first apply it to a set of synthetic data and then to a field case as discussed below.

Application to Synthetic Data
In this section, we utilize synthetic data created through a numerical model to demonstrate the practicality of our analytical solution and to evaluate the potential errors in the homogeneous aquifer solutions (Liang et al., 2022) in parameter estimation.The use of synthetic data in this context is advantageous because all the parameters in the synthetic example are known, enabling a more accurate assessment of the potential errors in the parameter estimation from our analytical solutions and the previous solution.We use COMSOL Multiphysics to numerically solve the nonlinear governing Equations 5 and 6 with the Gardner-Kozeny model (7).The model parameters used in the simulation are: b = 1 m, S y = 0.426, c u = 0 m 1 , γ = 0.036 m, S s = 1 × 10 5 m 1 , L = 800 m, G = 1.5 × 10 6 m, ε 0 = 1 × 10 8 , and τ = 0.518 days.The remaining parameters K 0 , c s , α k , and α c used in the simulation are provided in Table 3.In the numerical simulation, heterogeneity is represented by the depthdependent conductivity described in Equation 4, with given decay rates of c s .These parameters will be estimated using the analytical solutions.It should be mentioned that the impact of the decay of the hydraulic conductivity in the unsaturated zone has already been discussed previously.To simplify the analysis of the results, in this synthetic application, the decay of the hydraulic conductivity in the unsaturated zone is neglected by assuming c u = 0 m 1 .The focus here is primarily on the decay of the hydraulic conductivity in the saturated zone.

10.1029/2023WR036310
We first use the nonlinear numerical model to generate the synthetic time series of the hydraulic head induced by the M 2 tide at nine different depths (as depicted in Figure S1 in Supporting Information S1).Subsequently, we utilize the Baytap08 software (Tamura et al., 1991) to extract the tidal signals from these hydraulic heads and compute the corresponding amplitude ratio and phase shift at each depth as the synthetic data for the subsequent analysis.
Figure 9 displays the amplitude ratio and phase shift of the hydraulic head responses to M 2 tides at different depths so synthesized (circle symbols).The synthetic data are respectively fitted by the present solution (red solid curves), the present solution for a homogeneous aquifer (c s = 0 m 1 , blue dashed curves), and the solution of Liang et al. (2022) (black dashed curves).The parameters (K 0 , c s , α k , and α c ) are estimated by minimizing the sum of squared differences between the model-calculated and synthetic data, while maintaining the other parameter values the same as those used in the numerical model for both the present solution and Liang et al.'s solution.
Figure 9 demonstrates that the present solution aligns closely with the synthetic data, while both the homogeneous aquifer solutions much less so.Specifically, while the present solution for a homogeneous aquifer agrees reasonably with the synthetic data at c s ≤ 0.0025 m 1 (RMSE of amplitude ratio = 0.03 and RMSE of phase shift = 3.6°, Table 3), Liang et al.'s solution fails to match the synthetic phase shift (RMSE of phase shift = 10.7°).At c s = 0.01 m 1 , however, both the homogeneous aquifer solutions deviate significantly from the synthetic (RMSE of amplitude ratio = 0.07 and 0.07, RMSE of phase shift = 15.2°and14.9°, respectively).This result demonstrates that it may be important to include the vertical heterogeneity of aquifers in the model in order to realistically capture the tidal response of groundwater, especially for aquifers with more rapid decrease of conductivity with depth.
Table 3 summarizes the estimated parameters obtained using the present solution and the solution of Liang et al. (2022).It indicates that, compared to the homogeneous aquifer solutions, the present solution provides more accurate parameter estimates including K 0 , c s , α k , and α c .The slight deviation of the estimated parameters using the present solution from the known values in the numerical simulation may be attributed to the model's nonlinearity, which was linearized in the present governing Equations 9-11.Thus our result demonstrates the effectiveness of the present solution in handling aquifer heterogeneity and in obtaining reliable parameter estimations.
To illustrate the capacities of the models to estimate the average conductivity of a heterogeneous aquifer, we compare the estimated K of the homogeneous aquifer solutions to the mean value, K = 1 L ∫ 0 L exp (c s z) dz, of the heterogeneous conductivity.For the synthetic data with a small c s (0.0025 m 1 ), the estimated K (5.2 × 10 4 m/s) using homogeneous conductivity in the present solution (c s = 0 m 1 ) is slightly larger than the mean conductivity (K = 4.3 × 10 4 m/s) of the heterogeneous aquifer, while the estimated homogeneous K (4.1 × 10 4 m/s) using Liang et al.'s solution is slightly smaller than the mean conductivity of the heterogeneous aquifer.However, for the synthetic data with a large c s (0.01 m 1 ), the homogeneous aquifer solutions significantly underestimate K and α c , and overestimate α k .These findings suggest that when the decay rate of hydraulic conductivity is relatively small (c s = 0.0025 m 1 ), the solution for homogeneous aquifers can provide an approximate estimation of the mean conductivity.However, it cannot accurately estimate the parameters of the unsaturated zone (α k and α c ).On the other hand, when the decay rate of hydraulic conductivity is relatively large (c s = 0.01 m 1 ), the solution for homogeneous aquifers significantly underestimates the mean hydraulic conductivity of a heterogeneous aquifer and cannot accurately estimate the parameters of the unsaturated zone.In short, the solution for homogeneous aquifers may be useful for estimating the mean values of the parameters of heterogeneous aquifers with smaller decay rates of hydraulic conductivity, they fall short in accurately capturing the complexities of heterogeneous aquifers with larger decay rates.Our findings emphasize the importance of considering aquifer heterogeneity in Earth tidal analysis and highlight the effectiveness of the present model in handling such complexities.

Application to Field Data
Here we use the field data documented at Beishan, an arid area in the Gansu Province in western China (Figure 10a), to further demonstrate the applicability of the present solution.Beishan has been identified as the first priority area for establishing a high-level radioactive waste (HLW) repository in China.The safe disposal of HLW has garnered significant attention around the world, and geological disposal has emerged as a feasible and secure option.A series of comprehensive investigations has been conducted in the Beishan region, including hydraulic testing and long-term groundwater level monitoring (J.Wang et al., 2018).
In our application, we utilized the water level documented in borehole BSQ08 (Figure 11a).This borehole was drilled to a depth of 86.1 m and, as displayed in Figure 10b, it traversed through primarily granodiorite (74.56% of all the rock strata encountered) and minor amount of breccia (12.31%), lamprophyre (6.86%) and monzonitic granite (6.27%).The borehole reveals strata lacking significant aquitards, with groundwater across various lithologies showing connectivity.The hydrogeological conditions at this site indicate the presence of an unconfined aquifer (Guo et al., 2000).The following Earth tide analysis further supports the classification of this aquifer as unconfined and indicates that we should use the unconfined aquifer model to analyze the field data, evidenced by observed water levels displaying positive phase shifts in response to Earth tides.Additionally, a slug test employing a plugging packer system conducted in this borehole determined that the hydraulic conductivity of the strata at the borehole's bottom is extremely low (K ∼ 10 10 m/s, as elaborated below).This finding categorizes the bottom of the borehole as the aquifer's base.Water level in the borehole was documented using a pressure transducer at an hourly sampling frequency from 1 July 2018, to 30 July 2019 (Figure 11a), except from 30 October 2018 to 30 Mach 2019 (gray interval in Figure 11a), during which the water level was sampled daily and cannot be used for the extraction of tidal signal.The mean depth of the water table is 26.8 m, resulting in an average aquifer thickness of 59.3 m.The borehole BSQ08 screen fully penetrated the aquifer, which means that the water levels observed in the borehole represent the average hydraulic head across the entire aquifer.The hydraulic conductivity at different depths of the borehole was determined using the slug tests.A plugging packer system was employed at specific depth intervals of the aquifer in the slug tests, and the hydraulic conductivity was estimated for each interval.The conductivities within each interval from the slug tests show a general decline in conductivity with depth and are plotted in Figure 13.The estimated conductivities obtained through the slug tests will serve as a benchmark for evaluating the performance of the presented solution in terms of conductivity estimation.
The water level fluctuations in the borehole BSQ08, when examined in a magnified view in Figure 11a and the periodic components (Figure 11b) extracted from the raw water levels data, show a clear response to Earth tides.Spectral analysis of the raw water level data also shows a clear correspondence between the tidal components in the water level data and the Earth tidal components at the field location (Figure 11c) calculated with the Mapsis software (Lu et al., 2002).The distinct peak at the M 2 frequency (Figure 11c) is used for constraining our inversion of the hydraulic parameters discussed below.
The software Baytap08 (Tamura et al., 1991) was used to calculate the phase shift and the amplitude ratio of the water-level response to the M 2 tide in the borehole BSQ08 (Figure 11b).We have selected a window size of 30 days, with an overlap of 15 days between successive sampling windows, enabling the separation of the semidiurnal M 2 and S 2 tidal frequencies (Allègre et al., 2016;Xue et al., 2016).The calculated amplitude ratio and phase shift of the water-level response are plotted as functions of time as circles with error bar in Figure 12.The Water Resources Research 10.1029/2023WR036310 significant decreases in the phase shift and the amplitude ratio near October 2018 may be due to the gap between October 2018 and March 2019 where the tidal data are not available.Fitting the remaining data with the present solutions (solid curves) by least square method leads to a RMSE of 0.078 in amplitude ratio and a RMSE of 5.79°in phase shift.The unsaturated parameters (c u , a c , a k , γ, and S y ) and the saturated parameters (c s , K 0 and S s ) are estimated by minimizing the sum of the squared differences between modelcalculated and observed values and are summarized in Table 4.The comparison between the present model estimates (solid red lines in Figure 12) and the slug test results shows good agreement both in the unsaturated (c u = 0.5 m 1 ) and in the saturated zone (c s = 0.08 m 1 ).
Fitting of the data with the homogeneous solution of Liang et al. (2022) (dashed black line in Figure 12) leads to a RMSE of 0.081 in amplitude ratio and a RMSE of 5.81°in phase shift, nearly identical to that with the present model, as noted above.On the other hand, the homogeneous solution naturally cannot simulate the measured decreasing hydraulic conductivity with depth.The fact that both the homogeneous solution (Liang et al., 2022) and the present solution fit the tidal data equally well (Figure 12), despite the large differences in their conductivity profiles (as shown in Figure 13), suggests that the effectiveness of the tidal method in discriminating between different conductivity profiles may be limited.This limitation is primarily due to the deficiency in the observation data.Specifically, the hydraulic head measured in the borehole represents the average head across the entire aquifer.Consequently, the observed phase shift and amplitude ratios (Figure 12) are averaged values over the aquifer.The homogeneous solution can fit the tidal data but, by its nature, cannot estimate the depth-dependent heterogeneous hydraulic conductivity.Therefore, to accurately characterize the heterogeneous hydraulic conductivity using the tidal method, detailed observational constraints at different depths within the aquifer are required.Water Resources Research 10.1029/2023WR036310 Finally, the results of the present inversion of the hydraulic parameters may be useful for future applications such as for the simulation of the hydrogeology of the Beishan area in the Gansu Province.

Limitation of This Study
This study highlights the significant impact of the depth-dependence of hydraulic conductivity, an important form of aquifer heterogeneity, on the groundwater response to Earth tides.It should be noted, however, that real aquifers exhibit much greater heterogeneous complexities, such as layered arrangements of aquifers and aquitards.Freeze and Cherry (1979) also noted that the hydraulic conductivity may not always decrease with depth.Additionally, geological processes like structural disruptions and intrusive igneous rocks can introduce even more intricate lithological variations.Previous studies employ stochastic fields to describe aquifer heterogeneity (Yeh et al., 2015;Zhang, 2002).Therefore, the applicability of our analytical approach to such heterogeneous aquifers is subject to limitations.Nevertheless, our study holds value as a preliminary attempt to incorporate heterogeneity in an analytical model, shedding light on its influence on groundwater response to Earth tides.

Concluding Remarks
This study highlights the significant impact of the depth-dependence of hydraulic conductivity, an important form of aquifer heterogeneity, on the groundwater response to Earth tides.A novel analytical model is presented to explore the groundwater flow induced by Earth tides in both the unsaturated and the saturated zones with specific consideration of the impacts of the depth-decay hydraulic conductivity.The solutions are applied to field data to estimate the decay in conductivity with depth.The results indicate the importance of considering the depthdependence of hydraulic conductivity in tidal analyses.Our finding reveals that the impact of the capillary zone on the tidal response of groundwater is regulated by the heterogeneity in hydraulic conductivity and is primarily restricted to the shallow region of the aquifer.Neglecting the depth-dependence of conductivity in the existing models may lead to underestimation of the amplitude ratio and overestimation of the phase shift.The significant differences between the tidal response of a homogeneous aquifer and that of a heterogeneous aquifer, as shown in this study, suggest the potential of applying the present solution to the tidal response of an aquifer to characterize its depth-dependent heterogeneity.

•
Decay in conductivity with depth causes an increase in amplitude ratio and a decrease in phase shift • Influence of capillary zones on tidal responses is restricted to shallow aquifers due to decrease in conductivity with depth • Homogeneous solutions estimate mean conductivity well with smaller decay rates of conductivity, but underestimate with larger decay rates Supporting Information: Supporting Information may be found in the online version of this article.

Figure 2 .
Figure 2. A schematic diagrams showing unsaturated-saturated flow induced by Earth tides and the change of hydraulic conductivity with depth.

a
The K φ = 1E 3 m/s at the ground surface (z = b) remains constant and K 0 = K φ exp( c u b).

Figure 3 .
Figure 3.Comparison of u 0 or h 0 profiles of the numerical solution (circle symbols) with the present solution (solid curves) for (a) different c u = c s (0.004, 0.005, and 0.006 m 1 ) and a given α k = 3 m 1 , and (b) different α k (2, 3, and 4 m 1 ) and given c u (0.005 m 1 ) and c s (0.08 m 1 ).

Figures
Figures 6a and 6b respectively present the changes of the amplitude ratio and phase shift with the normalized depth (z/L for the saturated zone and z/b for the unsaturated zone) for different values of c s , where the parameters

Figure 4 .
Figure 4. Comparison of the fluctuating hydraulic head (h) at z = 500 m of the numerical solution using GK model with the present solution (solid curves) for (a) different c u = c s (0.004, 0.005, and 0.006 m 1 ) with a given α k = 3 m 1 and (b) different α k (2, 3, and 4 m 1 ) with given c u (0.005 m 1 ) and c s (0.08 m 1 ).

Figure 5 .
Figure 5.The profiles of the amplitude ratio (left column) and phase shift (right column) of the tide response of hydraulic head in an aquifer with different c u and b in the unsaturated zone, where c s = 0.005 m 1 .

Figure 6 .
Figure 6.The profiles of the amplitude ratio (a) and phase shift (b) of the tide response of the hydraulic head in the saturated zone for different c s .Also shown is the model for a half-space aquifer (black dashed curves, Liang et al., 2022) with a constant conductivity of 1 × 10 3 m/s.

Figure 7 .
Figure 7.The profiles of the amplitude ratio (left column) and phase shift (right column) of the tide response of hydraulic head in an aquifer with different α k in the unsaturated zone, where panels (a) and (b) display the predictions by the solution of Liang et al. (2022), while the remaining graphs demonstrate the predictions by the present solution using different values of c s .The other paramters are as follows: b = 1 m, α c = 1 m 1 , L = 800 m, K 0 = 1E 4 m/s, S s = 1E 5 m 1 , and S y = 0.426.

Figure 8 .
Figure 8.The profiles of the amplitude ratio (left column) and phase shift (right column) of the response of hydraulic head in an aquifer with different thicknesses of the vadose zone, where panels (a) and (b) display the predictions by the solution of Liang et al. (2022), while the remaining graphs demonstrate the predictions by the present solution using different values of c s .The other paramters are as follows: α k = 3 m 1 , α c = 1 m 1 , L = 800 m, K 0 = 1E 4 m/s, S s = 1E 5 m 1 , andS y = 0.426.

Figure 9 .
Figure 9.The amplitude ratio (a, c) and phase shift (b, d) of hydraulic head responses to M 2 tides at different normalized depths, as predicted by the numerical model (circle symbols).The synthetic data predicted by the numerical model are respectively fitted by the present solution (red solid curves), the present solution for a homogeneous aquifer (c s = 0 m 1 , blue dashed curves), and the solution of Liang et al. (2022) (black dashed curves).

Figure 10 .
Figure 10.The basic information of the studied boreholes BSQ08: (a) Location and elevation map; (b) Lithological log.

Figure 11 .
Figure 11.Time series of (a) raw water level data from borehole BSQ08 and (b) residual tidal signals in water levels; and (c) comparisons of the spectra between water levels in borehole BSQ08 and Earth tides (ET, shown in black) where the x-axis represents frequency in cycles per day (cpd).

Figure 12 .
Figure 12.Comparisons of (a) phase shift and (b) amplitude ratio (circle symbols with error bars) for the tidal response to the M2 tide in BSQ08 with calculated results using the present solution (solid curves) and Liang et al.'s solution (dashed curves).

Figure 13 .
Figure 13.Estimated hydraulic conductivity at varying depths obtained through slug tests conducted in borehole BSQ08 and the corresponding estimation results for hydraulic conductivity using the present solution (solid curves) and Liang et al.'s solution (dashed curves).

Table 2
Model Parameters Used in Figures3-8

Table 3
Liang et al. (2022)meters Used in the Numerical Model for the Synthetic Case and the Corresponding Estimated Parameters Obtained Using the Present Solutions and the Solution ofLiang et al. (2022) WANG ET AL.

Table 4
Liang et al. (2022)mated Parameters Obtained Using the Present Solutions and the Solution ofLiang et al. (2022)in Figure12