Groundwater Responses to Earth Tides: Evaluation of Analytical Solutions Using Numerical Simulation

Harmonic Earth tide components in well water levels have been used to estimate hydraulic and geomechanical subsurface properties. However, the robustness of various methods based on analytical solutions has not been established. First, we review the theory and examine the latest analytical solution used to relate well water levels to Earth tides. Second, we develop and verify a novel numerical model coupling hydraulics and geomechanics to Earth tide strains. Third, we assess subsurface conditions over depth for a range of realistic properties. Fourth, we simulate the well water level response to Earth tide strains within a 2D poroelastic layered aquifer system confined by a 100 m thick aquitard. We find that the non‐linear inversion of analytical solutions to match two observations (amplitudes and phases) to multiple unknown parameters is sensible to the initial guess. We reveal that undrained, confined conditions are necessary for the analytical solution to be valid. This occurs for the dominant M2 frequency at depths >50 m and requires specific storage at constant strain of Sϵ ≥ 10−6 m−1, hydraulic conductivity of the aquitard of kl ≤ 5 ⋅ 10−5 ms−1 and aquifer of ka ≥ 10−4 ms−1. We further illustrate that the analytical solution is valid in unconsolidated systems, whereas consolidated systems require additional consideration of the Biot modulus. Overall, a priori knowledge of the subsurface system supports interpretation of the groundwater response. Our results improve understanding of the effect of Earth tides on groundwater systems and its interpretation for subsurface properties.

10.1029/2022JB024771 2 of 22 (Rojstaczer, 1988;Rojstaczer & Riley, 1990).For a given frequency, the ratio between the measured head at the observation well and the confined pore pressure produced by the change in strain is known as amplitude ratio (i.e., h 1 /P 1 in Figure 1), whereas the time lag between these parameters is known as phase shift (i.e., t 2 − t 1 in Figure 1).If field data is available, harmonic least squares (HALS) with the main tidal components can be directly applied to the time series to obtain the amplitude ratio and phase shift between both signals (Schweizer et al., 2021;Turnadge et al., 2019).
Field measurements of the groundwater response to Earth tides has resulted in negative and positive phase shifts between strain and well water levels.Negative phase shifts are interpreted as fully confined conditions and horizontal flow only (Hsieh et al., 1987).For example, hydraulic conductivity and specific storage were estimated from negative phase shifts (Roeloffs et al., 1989).However, positive phase shifts in the field were also observed and attributed to vertical flow through leaky aquitards (Allègre et al., 2016;Xue et al., 2016).This was interpreted using the analytical solution of vertical flow in an homogeneous aquifer caused by a harmonic load or stress ignoring the effect of the observation well (Wang, 2017).
Early research developed and tested an analytical solution to estimate subsurface properties from the relationship between strain and water levels in observation wells in a fully confined aquifer (Cooper Jr. et al., 1967;Hsieh et al., 1987).This was extended to include leaky conditions and allow concurrent use of multiple frequencies (Rojstaczer, 1988;Rojstaczer & Riley, 1990).Rojstaczer and Agnew (1989) studied the dependency of porosity and elastic parameters to a real deformation of a poroelastic medium.High sensitivity was reported when the applied strains occurs in low porosity and the increase with decreasing compressibility (inverse of the bulk modulus) of the solid matrix.
Recent research included modifications to the original analytical solution by Hsieh et al. (1987) to account for more realistic conditions.Most notable is the work by Wang et al. (2018) who developed an extended analytical solution which includes vertical leakage to model a two-layered aquifer system.Gao et al. (2020) investigated the well skin effect which originates from the fact that the formation around a well is disturbed and well water storage on the well water level response to Earth tides.The authors found that the skin effect may significantly delay the well water level phase response to Earth tides.In addition, Guo et al. (2021) developed a model for tidal response with a fault passing through the aquifer based on a fault-guided fracture network to estimate fracture properties.They found that the hydraulic diffusivity in the fault damage zone higher than previously established Figure 1.Representation of the pore pressure change and well water level in a semi-confined aquifer due to Earth tide strains.Earth tides induce subsurface stress which results in strains that generate changes in the pore water pressure.This leads to pressure gradients between the subsurface and observation wells that cause water movement in or out of the well.The ratio between the Earth tide strain and the well water level response can be expressed as amplitude ratio whereas the time difference between both signals as phase shift.values, but also that it remains below estimates based on induced seismicity migration.Sun et al. (2020) reviewed four of the most common analytical models to estimate hydraulic properties with Earth tidal analysis.They estimated hydraulic properties from a real data set and provide a range of applicability of the different models based on the transmissivity of the aquifer.
While analytical models offer a convenient approach to estimate hydraulic properties, their applicability is limited through simplifying assumptions arising from fundamental physics, conceptual model or boundary conditions.These include for example, only radial flow, negligible horizontal displacement of the aquifer, confined and undrained conditions, unconsolidated materials and no gravity.Moreover, it has been reported that some of these assumptions may significantly influence the estimated subsurface properties (Wang et al., 2019;Zhu & Wang, 2020).
Numerical modeling of tidal effects is common in coastal (Abarca et al., 2013;M. Zhang et al., 2021) and adjacent settings (Alcaraz et al., 2021;Jardani et al., 2012;Pendiuk et al., 2020), likely because the loading effect of ocean tides is much larger than that caused solely by Earth tide forces.So far, modeling the subsurface response to ocean tides has considered poroelastic conditions and harmonic loading by the weight of the water, therefore, solved as a consolidation problem.In contrast, the pressurization of an inland aquifer is produced by changes in the pore space volume of the porous material due to strains (i.e., eigenstrains) caused by the gravitational influence from the movement of celestial bodies.Moreover, the change of the confined pore pressure is generally measured using an observation well which causes fluid movement.
To the best of our knowledge, only hydraulic modeling approaches neglecting any geomechanical effects, that is, groundwater flow without coupled poroelasticity, have so far been used to investigate the groundwater response to Earth tides.For example, Wang et al. (2019) simulated the effect of capillarity of the unsaturated zone in one dimension.They found out that the assumption of fixed water table can lead to erroneous estimation of subsurface properties with analytical solutions.Zhu and Wang (2020) simulated a multi-layered system to study the effect of leakage through an aquitard and concluded that simplifications in the analytical model lead to overly conservative estimates of vertical flux between layers.Wang and Manga (2021) provide a summary of these works.
The confined pore pressure generated as result of Earth tide strains is a mechanical phenomenon caused by the elastic deformation of the porous matrix.Furthermore, unlike for traditional hydraulic testing approaches, there is a general lack of work investigating the effect of realistic conditions and assumptions on interpretations using analytical solutions.Investigating the influence of limiting assumptions and realistic subsurface conditions to better understand the applicability and robustness of analytical solutions requires development of more advance numerical models that also consider coupling with geomechanics.
The objective of this study is therefore to (a) critically examine assumptions upon which analytical solutions are based, (b) develop a numerical model for the groundwater responses to Earth tides, which couples hydraulic and geomechanical processes, (c) investigate and compare response conditions as well as the influence of geomechanical properties.Thus, our work significantly improves our understanding of the coupled physics, which controls the well response in a poroelastic medium.These results can act as a practical guide for improved estimation of aquifer properties due to the groundwater response to Earth tides.

Fundamental Theory of the Groundwater Response to Earth Tides
Earth tides are displacements of the solid Earth's crust caused by the gravitational forces of celestial bodies that move in relation to the Earth.Such displacements are typically expressed as harmonic signals that can be predicted from well-known astronomical relationships.Earth tide forces are dominant at distinct frequencies within the semi-diurnal and diurnal range, for example, M 2 at 1.97322 cpd or S 1 at 1.0 cpd.Under tidal forcing, the poroelastic space and the porous material elastically deforms depending on the mechanical properties of the system resulting in a small change of volume.If the subsurface is saturated, the filling fluid has to adapt to the new available pore space which raises or lowers the pore pressure.The processes can be mathematically represented by the Biot consolidation theory.Biot (1941) developed the constitutive laws which relate the applied forces (stresses) with deformation (strains) and motion within a elastically compressible porous medium.For the purpose of modeling, these laws are BASTIAS ESPEJO ET AL.

10.1029/2022JB024771
4 of 22 formulated in the form of mathematical equations which consist of four basic variables; total stress (σ ij ), strain (ϵ ij ), pore pressure (p f ) and increment of fluid content (ξ).The mechanical variables (stress or strain) can be related with one of the fluid quantity (pore pressure or fluid content) to form independent variable and therefore mathematical equations.For the particular case of Earth tides, is convenient to express the poroelastic equations in terms of total stress and pore pressure as independent variables, also termed pure stiffness formulation.The basic relation between total stress and pore pressure is, where   ′  is the effective stress, δ ij is the Kronecker delta and α is the Biot poroelastic coefficient, K s is the solid material bulk modulus, K is the porous medium bulk modulus.The latter is related to the undrained bulk modulus, K u , as here, M, is the Biot modulus which is defined as, where n is the porosity, K f is the bulk modulus of the fluid and B is the Skempton coefficient.
The effective stress,   ′  , is related to elastic strain by the generalized Hooke's law through the compliance tensor C ijlk (Berryman, 1999;Dropek et al., 1978), where G is the solid material shear modulus, ϵ is the volumetric strain (ϵ = ϵ xx + ϵ yy + ϵ zz ).Combining Equations 1 and 5 results and the dynamic property p f can be related to the volume change as Note that Equation 6is the effective stress equation, but for convenience is expressed in terms of total stress and the effective stress is reduce to the first and second terms of the right hand side.
Fluid transport is modeled with the fluid balance equation as where k p , ij is the porous medium permeability tensor, μ is the fluid viscosity, Q represents sinks or sources, ϵ ij = ∇u i relates strain with displacement, often preferred in simulation codes (Flemisch et al., 2011;Keilegavlen et al., 2021;Kolditz et al., 2012;Verruijt, 2013), and the effective stress can be related to the strain though the Hooke's law (i.e., Equation 5).S ϵ is the specific storage at constant strain and is related to the Biot modulus as Equations 6-8 can be solved in a coupled manner with appropriate boundary conditions and represent the elastic deformation and fluid movement in a porous medium.

Analytical Solution
When uniaxial-vertical strain and zero incremental vertical stress are assumed (this occurs only when one of the principal stresses is non-zero and the stress does not change with depth), the subsurface is mechanically restricted to move only in the vertical direction, for example, land surface subsidence due to consolidation occurs primarily in the vertical direction (Herrera-García et al., 2021).Under such conditions ϵ xx = ϵ yy = 0 and σ zz = 0 which leads to a simplified version of Equations 6 and 7 as and Combining Equations 10 and 11 to eliminate ϵ zz gives where This is the definition of storage coefficient in hydrology (Cheng, 2016;Verruijt, 2013;Wang, 2017).The specific storage, S s , is obtained when the specific weight of the fluid is considered as where ρ f is the density of the filling fluid and g is the Earth's gravitational acceleration constant.
With this derivation, we stress the conceptual difference between the specific storage at constant strain (Equation 9) and the storage coefficient with uniaxial strain (Equation 13).Please note that S approaches S ϵ when K ≫ K f , hence the second term of Equations 4 and 13 go toward zero.Physically, in such cases the amount of fluid coming out of storage will only depend on the fluid compressibility and the porosity of the porous medium, as the porous material is rigid (Freeze & Cherry, 1979;Lambe & Whitman, 1991;Verruijt & Van Baars, 2007).
Hydraulic head,   , can be used as a proxy for pore pressure in Equation 8,  =  −1   −1 .Moreover, in a confined aquifer with a constant thickness H a , hydraulic conductivity, k a = k p,ij ρ −1 g −1 , can be express in terms of transmissivity (T = k a H a ) and specific storage at constant strain in terms of storativity at constant strain S ϵ,t = S ϵ H a as (Cheng, 2016;Verruijt, 2013;Wang, 2017), For practical reasons, Equation 15 can be reformulated into cylindrical coordinates assuming only radial flow (C.E. Jacob, 1946;C. Jacob & Lohman, 1952), also the effect of a leaky layer on top of the aquifer can be included as vertical leakage in the sink/source in terms of hydraulic conductivity of the layer on top (k l ) and thickness of such layer (H l ) expecting In this work we use the terms aquifer and aquitard to reflect layers of higher and lower hydraulic conductivity, respectively, as is consistent with the terminology used in previous works.Equation 16was used to derive analytical solutions that describe the groundwater response to Earth tides in a fully vertical and horizontal confined aquifers (Hsieh et al., 1987) and in aquifers with vertical leakage (Wang et al., 2018).A detailed derivation of the more versatile leaky solution is presented in Wang et al. (2018), but is also included in Appendix A of this work.The solution describes the water level in a well (h w ) in terms of the hydraulic conductivity of the aquifer (k a ), hydraulic conductivity of the aquitard (k l ), the specific storage at constant strain of the aquifer (S ϵ ) and the BASTIAS ESPEJO ET AL. 10.1029/2022JB024771 6 of 22 geometry of the well (Equations A2, A3, and A4).In this formulation, any effects arising from well skin or storage are neglected, but such effects have been investigated in the literature (Gao et al., 2020).Fluid level changes in the well are caused by forces generated at the far field (far away from the radius of influence of the observation well).Assuming undrained conditions (ξ = 0, in Equation 7) and α = 1 representing unconsolidated systems, (e.g., sands, gravels and clays) such changes can be quantified as where ϵ G is a external volumetric strain.The change of water level in a well due to an areal strain is graphically shown in Figures 2a and 2b for two given times from t 0 = 0 to t = t with gravitational strains from ϵ G = 0 to ϵ G = ϵ G (t).
The amplitude ratio and phase shift between the piezometric head at a distance from the well beyond its radius of influence and the water level in the well are expressed as (Hsieh et al., 1987), and Here, ϵ 0 is the amplitude of the volumetric strain signal, which can be obtained from software based on tidal catalogs, for example, PyGtide (Rau, 2018), ETERNA PREDICT (Wenzel, 1996), TSoft (Van Camp & Vauterin, 2005) and h w is the fluid level in the well obtained by the analytical solution (Equation A5).The analytical solution is subject to the assumptions under which it was derived: (a) undrained conditions, (b) the confining layer has negligible specific storage, (c) the flow is horizontal in the aquifer, (d) the well is represented by a line with length matching the aquifer extent, (e) the deformation is only vertical, (f) no external forces such as gravity.
Moreover, field measured amplitude ratio, A, and phase shift, Δϕ, can be obtained from applying HALS to the Earth tide strain and hydraulic head time series (Schweizer et al., 2021).The obtained phase shift varies between −π ≥ ϕ ≥ π, thus, only fluid capable to flow from the aquifer to the observation well within that time frame (half a day for M2) is going to contribute in a change of the amplitude and phase shift.The latter bounds the scale of the method as high conductivity aquifers can move fluid from higher distances than low conductivity aquifers.
Using the results obtained from HALS, Equations 18 and 19 can be inverted to estimate constant values of Equation 16using any approach suitable for non-linear algorithm estimation (gradient methods), for example, least-squares.In fact, this approach has been used to estimate aquifer hydraulic conductivity, specific storage and aquitard leakage from the amplitude and phase response of groundwater to Earth tides (Rau, et al., 2020a;Y. Zhang et al., 2021).However, the task of solving these non-linear equations is an ill-posed problem, because the solution might not be unique.The computation of meaningful approximate solutions of inversion of Equations 18 and 19, therefore, can be quite challenging and strongly dependent on the a suitable initial guess feed to the non-linear algorithm.Moreover, it is not always a simple matter to decide which initial guess to choose.We note that the implications of initial guess nor the performance of iterative methods have been investigated for practical use of this solution.
Gradient methods such as the Levenburg-Marquardt often use in least-squares to numerically search for the nearest (local) minimum to the given initial condition and are readily implemented in SciPy (Virtanen et al., 2020).The function fitting model finds the local or global minimum, which depends on the feed initial guess.To investigate the performance of least-squares by the initial guess on the parameter estimation of Wang et al. (2018) (Equations 18 and 19), we systematically explore the solution space of a fitting function.The fitting function was formulated as and where A(k a , S ϵ , k l ) and Δϕ(k a , S ϵ , k l ) are Equations 18 and 19.The aquifer thickness and aquitard depth were arbitrarily defined to be 1 and 100 m, respectively; the radius of the well to 0.2 m, A obs and ϕ obs are objective amplitude ratio and phase shift values given by the user.The Earth tide frequency was set to 1.93 cpd, previous studies have suggested that by considering an extra tide frequency, for instance, 1 cpd might constrain better the minimization problem (Y.Zhang et al., 2021).But this is not explored in this study.
The least-squares solver minimizes the difference between A obs − A(k a , S ϵ , k l ) and ϕ obs − Δϕ(k a , S ϵ , k l ) of Equations 20 and 21 by iterating through a combination of values of k a , S ϵ , and k l .An array consisting of discrete values within realistic ranges for amplitude ratio (0.001 ≤ A ≤ 1) and phase (−90° ≤ Δϕ ≤ 90°) where generated.
For each pair of amplitude ratio and phase shift, the solution space was solved using least-squares of SciPy.The tolerance for termination by the change of the cost function was set to be 3 ⋅ 10 −6 and units for 20 and 21 where set to days so as to increase the numeric values and avoid errors.The sensitivity of the method was tested by generating 1,000 samples of k a , S ϵ and k l generated by a random log uniform distribution ranging from 1 ⋅ 10 −7 ≤ k a ≤ 1 ⋅ 10 −3 ms −1 , 1 ⋅ 10 −7 ≤ S ϵ ≤ 1 ⋅ 10 −5 m −1 and 1 ⋅ 10 −8 ≤ k l ≤ 1 ⋅ 10 −4 ms −1 .Each trio of samples was set as initial condition in Equations 20 and 21.Starting from different initial values allows the solver to find potentially different local minima.The outputs were stored and the maximal difference between each estimated property was used as a proxy for performance of the non-linear search.This approach allows identification of values within the solution space where several local minimum might be encounter.

Numerical Model
The generic equations presented in Section 2.2 can be solved analytically for specific boundary conditions (Equation 16), but analytical solutions do not accounts for the elastic deformation of the porous medium (Section 2.1, BASTIAS ESPEJO ET AL. 10.1029/2022JB024771 8 of 22 Equations 6-8).The coupled physics between fluid movement and mechanical deformation can be solved numerically.Here, we develop a novel numerical approach for simulating the groundwater response to Earth tides.This allows a more realistic physical representation compared to the limiting assumptions of the analytical solution presented in Section 2.2 and advances the previous study by Wang & Manga (2021).The aim is to investigate and establish robustness of the analytical solution when interpreting the groundwater response to Earth tides.
When modeling Earth tides, an external gravitational strain ϵ G (x, y, z, t) is applied to deform the Earth's crust and the resulting fluid pressure p f (x, y, z, t), and the displacement vector u ii (x, y, z, t) is calculated.Under the free surface condition, the normal stress along the radius is zero.Hence, gravitational strain can be decomposed into its vertical ϵ h (x, y, t) and horizontal component ϵ v (z, t) (Agnew, 2005), Earth tides induce an eigenstrain, that is, a strain that does not result directly from an applied force.Qu and Cherkaoui (2006) describes the differences and relationships between total, elastic and eigenstrains.To simulate the effect in a realistic well-aquifer system, in which the hydraulic and geomechanical properties of the material may vary in space, we applied vertical and horizontal strain as displacement boundary conditions.This fixes the internal strain throughout the model as a function of the filling material elastic tensor as (Wang, 2017), and where R and L are the horizontal and vertical lengths of the model, respectively.A constant atmospheric pressure (i.e., drained condition) is assumed at the top of the modeling domain and at the top of the observation well.Note that we exclude consideration of barometric variations, such as could be caused by atmospheric tides, which is a valid assumption for the M2 frequency discussed here (Rau, Cuthbert, Acworth, & Blum, 2020b), The governing equations follow the traditional Biot (1941) theory of a linear elastic, saturated and deformable porous medium with water as the fluid (Cheng, 2016;Wang, 2017).The strong form of the general equations in Section 2.1 can be converted into the respective weak form and discretized before solving with the finite element (FE) method.In this study, we adopt the continuum representation of an elementary volume (REV) in a porous medium.To solve the numerical system the Real Heterogeneity App (RHEA), a FE application based on the Multiphysics Object Oriented Simulation Environment (MOOSE) was used (Permann et al., 2020).A detailed description of the system of equations to be solved as well as further information of the tight coupling and numerical description of the FE implementation utilized in this study can be found in (Bastías Espejo et al., 2021;Wilkins et al., 2020Wilkins et al., , 2021)).For consistency, we keep the original notation used in Bastías Espejo et al. ( 2021), where the field variables are the fluid pressure p f and the displacement vector u ii .

Assessing the Subsurface Response Conditions
Within the theory of linear poroelasticity (Equations 6-8), one can distinguish between two end-members that describe the type of pore pressure response to stresses and strains: undrained and drained.When a deformation in a porous medium occurs, under drained conditions, the rate of applied distortion is slow in relation to the ability of the porous medium to allow dissipation of the pressure gradient.This results in the flow of fluid caused by the buildup of pressure differences.Under undrained conditions, the rate of applied distortion is fast enough for an instantaneous pore pressure response to external deformations and fluid cannot flow in response generating a constant confined pore pressure throughout the porous medium.The type of response is of importance as analytical solution often assume undrained systems, however most aquifers in nature are draining in at least one boundary, further its relevance has not yet been studied.
BASTIAS ESPEJO ET AL. 10.1029/2022JB024771 9 of 22 The type of subsurface response is represented by Equation 7, in which ξ describes an increment of change in fluid content.Under undrained conditions, ξ = 0 and Equation 7 reduces to Equation 17.While under drained conditions ξ ≠ 0. Hence, the confined pore pressure as a response to Earth tides can be obtained with Equation 7if the Earth tide strain (ϵ G ) is known.
For Earth tides, the applied strain depends on the frequency of the harmonic and the type of response is a function of the hydro-geomechanical subsurface properties as well as depth.To assess the type of response under realistic conditions, we numerically modeled an infinitely long 1D (5,000 m) section of the subsurface (Figure 2c) using a harmonic function as displacement boundary condition as follows Here, ϵ 0 is the amplitude of the Earth strain,   2 is the frequency of the M 2 component and t is the time in days.We computed the increment of fluid content, ξ, over depth 0 ≤ z ≤ 1,000 m and repeated the calculations by setting assumed, but realistic discrete values of specific storage at constant strain S ϵ , (1 ⋅ 10 −7 m −1 , 1 ⋅ 10 −6 m −1 , and 1 ⋅ 10 −5 m −1 ), and bulk modulus K, (1 ⋅ 10 9 Pa, 1 ⋅ 10 10 Pa, and 1 ⋅ 10 11 Pa).These values represent realistic conditions as reported in the geoscience literature (Cheng, 2016;Das & Das, 2008;Lade, 2001;Wang, 2017).We note that we use the term aquitard even though the assigned values of its hydraulic conductivity are relatively high (1 ⋅ 10 −7 ≤ k l ≤ 1 ⋅ 10 −4 ms −1 ).In our work the term aquitard reflects the fact that it is used as a layer with lower hydraulic conductivity values compared to the aquifer and its use is consistent with the terminology and values found in previous literature, for example, (Wang et al., 2018).
When analyzing the groundwater response to Earth tides, undrained conditions have to be given for the analytical solution to be valid (Appendix A).However, the type of response to Earth tides has not been assessed before and is therefore unknown.To assess whether the subsurface response is sufficiently undrained, we can use Equation 7assuming α = 1.Under undrained conditions ξ = 0, such condition can be assumed when ϵ ≫ ξ.Hence, the effect of ξ can be neglected in Equation 7.

Numerical Model of a Coupled Well-Aquifer System
To investigate the limitations of the analytical model presented in Section 2.2, a 2D axial-symmetric cylindrical model was developed.The model accounts for poro elastic aquifer (Section 2.1) bounded by a low permeable aquitard on top and by a rigid aquiclude on the bottom.Gravity, vertical and horizontal flows are allowed (Figure 2d).Boundary conditions are set as described in Section 2.5.The applied Earth strain (ϵ G ) was calculate theoretically with PyGtide (Rau, 2018), a Python wrapper for ETERNA PREDICT 3.4.We chose the city of Berlin (Germany) and a signal frequency of 2 cpd for simplicity as it closely resembles M 2 with a duration of 30 days.While the location is arbitrary, it does not change the conclusions because the context of this study is generic.
The borehole-subsurface system is modeled as a 1D element outside the 2D system.To relate the porous medium and the borehole-subsurface, the boundary at r = 0 is modeled as a free drainage boundary, that is, a sink boundary where the flux is computed in function of the pressure at r = 0 (p r = 0 ) and at the bottom of the borehole-subsurface (p w ) where    is the mass flux, C the conductance (i.e., how efficiently fluid is transported though a boundary) between the borehole-subsurface and the model.For this study C = 10 −3 m 2 s −1 which is high enough to ensure that p f = p r = 0 at the end of each non-linear iteration.
As a result, the mass flux through the boundary between the model and the well is computed in every non-linear iteration, which fixes the pore pressure at the boundary for the next iteration.This way, the fluid level in the well is tightly coupled to the pressure at the well boundary of the model as the fluid level in the well is computed in the same Jacobian matrix with the numerical model as, BASTIAS ESPEJO ET AL.
10.1029/2022JB024771 10 of 22 Here, h f represents the water level in the well , A c is the external area of a cylinder and A w is the inner area of a cylinder.Since the model is linear elastic, typically, only two non-linear iterations are needed.
The model domain is R = 5,000 m long in the r direction and L = 101 m in depth in the z direction.The aquitard is 100 m thick, whereas the aquifer is 1 m thick, and we assume that the well is screened throughout this unit.The 101 m long well is located at the left boundary of the modeling domain and the well has r w = 0.2 m of radius (Figure 2d).The geometry complies with previous studies (Hsieh et al., 1987;Wang et al., 2018) and therefore enables a comparison.The finite elements were discretized using the built-in mesh generator of MOOSE and the element size increases logarithmically along the r-axis away from the well.The mesh is vertically and logarithmically discretized 5 times across the aquifer, 20 times across the top layer and 100 times in the horizontal direction.The material properties of the model are summarized in Table 1.The values of Table 1 were assumed in previous studies (Wang et al., 2018) and extracted from literature (Cheng, 2016;Das & Das, 2008;Lade, 2001;Wang, 2017).
The initial pore pressure condition is set as   0  = −  and the effective initial stress as   ′0  = ( −  )  , where   ′0  is the vertical component of the effective stress at time zero.Again, we apply a harmonic displacement function with the M 2 frequency computed with a tidal catalog, the amplitude of the strain is ϵ 0 = 1.2 ⋅ 10 −8 .The model runs until it reaches quasi steady-state, at which point the well physics as well as tidal forcing as boundary conditions are activated.This approach minimizes potential numerical overshooting produced by the free drainage boundary between the porous medium and the well.We verify this numerical implementation using the analytical solution of Wang et al. (2018) (Subsection 2.2) with the aquitard permeability set to zero, that is, the model represents only one layer.

Analytical Solutions
Previous research has used the analytical solutions by Hsieh et al. (1987) and H. F. Wang (2000), Wang et al. (2018) to estimate hydraulic properties for negative and positive phase shifts between groundwater and Earth tides, respectively.We note that the Hsieh et al. (1987) and Wang et al. (2018) require undrained conditions, whereas H. F. Wang (2000) drained conditions.The robustness of these assumptions have not been investigated for Earth tide frequencies.For drained conditions the relationship between stress and strain is no longer linear, as pore pressure also plays a role bearing loads (see Equation 7).Furthermore, while H. F. Wang (2000)  one dimensional poroelastic aquifer, it neglects the influence of an observation well.As shown in Figure 3, a well generates phase shifts between the confined far distance pore pressure and the water level in the well as the fluid requires time to move in and out of the well.Strictly speaking, this solution was derived for surface loads, such as exerted from atmospheric pressure, but not for Earth tide strains.These aspects illustrate that H. F. Wang (2000) has limited use when estimating hydraulic properties from the groundwater response to Earth tides.Wang et al. (2018) provides an extended formulation to Hsieh et al. (1987) considering vertical aquitard leakage accounting for both negative and positive phase shifts.It is useful to illustrate the solution space of Wang et al. (2018) by providing an overview of amplitude ratios and phase shifts (Equations 18 and 19) as a function of realistic ranges of the aquifer hydraulic conductivity (k a ) and specific storage at constant strain (S ϵ ) as well as discrete values of leakage, see Figure 3.Note that this is based on the dominant harmonic signal frequency of 2 cpd, a well and screen radius of 0.2 m, a screen length of 1 m and an aquitard thickness of 100 m.The first row, Figures 3a and 3b, shows the case where there is no vertical leakage leading to negative phase shifts only.
We confirm the reports by Wang et al. (2018) that the analytical solution matches the previous solution by Hsieh et al. (1987) when the aquitard hydraulic conductivity is set to zero.
The solution space shows that vertical leakage causes positive phase shifts at relatively high aquifer hydraulic conductivity, that is, k a > 1 ⋅ 10 −5 ms −1 in Figure 3d.This threshold is even more clear for vertical leakage larger than k l > 1 ⋅ 10 −6 ms −1 where the transition to positive phase shift is almost linear.Moreover, in Figures 3d and 3f, the phase shift behavior is very similar for the lower part of the specific storage at constant strain under study (S ϵ < 1 ⋅ 10 −5 m −1 ).A similar case is observed in Figures 3c and 3e where the amplitude response of the analytical solution shows very similar results.
The solution space illustrated in Figure 3 shows that the functions are non-linear and, therefore, to estimate subsurface parameters a gradient root finding method is necessary.This is based on an iterative method resulting in an approximate solution only.The search is based on an initial guess and the method may find different potential solutions depending on the function gradients.We investigate the effect of the initial guesses on the solution space by visualizing the difference in the found local minimum.Figure 4 shows the variability in estimated parameters due to providing different initial guesses for least-squares solving.Blue color means less variability whereas red color shows a larger difference in the solution space and therefore higher number of local minimum.
In general, higher sensitivity to the initial values is observed in the phase shifts compared to the amplitude ratios.
The sensitivity for inverting the hydraulic conductivity of the aquifer k a is relatively low (Δk a < 5 ⋅ 10 −4 ms −1 ) at negative phases and increases as the phase shift approaches 0°.For values higher than zero degrees the sensitivity decreases again until approximately 40° where it starts to increase again (Figure 4a).The highest sensitivity is found around amplitude ratio of one and positive phase shift.Situations where the amplitude ratio is one and the phase shift much higher than zero are not realistic and should be disregarded.
Specific storage at constant strain shows a high contrast in solution variability with values that are very close to ΔS ϵ = 1 ⋅ 10 −5 m −1 for most of the solution space (Figure 4b).At low phase shifts (Δϕ < −70°) the variability significantly reduces to ΔS ϵ = 1 ⋅ 10 −7 m −1 .In practice, most of the realistic cases will fall within the high variability zone.Vertical leakage shows relatively low sensitivity where the phase shift is negative (Δk l < 5 ⋅ 10 −5 ms −1 , Figure 4c).However, at positive phase sifts the variability increases up to two orders of magnitude, demonstrating the effect of the phase shift on vertical leakage.
In the illustrated case, the sensitivity of the specific storage at constant strain is constant throughout the solution space and therefore its initial value does not play a significant role on finding different solutions.Therefore, is likely that for each specific storage at constant strain study a local minimum was found during the minimization.We, therefore, advise special attention when selecting the initial value of the specific storage in order to obtain a meaningful result.This value can be bound if knowledge of porosity is available (Section 3.2).Negative phase shifts show a low sensitivity to the initial condition and will likely result in an accurate inversion of the hydraulic properties without a priori knowledge of the subsurface properties as mentioned in Section 2.2.For positive phase shifts a handle on at least one of the properties is necessary as the vertical leakage significantly increases its variability.Further, the hydraulic conductivity increases its variability toward high amplitudes which is the range where Earth tide methods work best.Overall, this complies with the previous finding that positive phase shifts can robustly be interpreted as vertical leakage (Wang et al., 2018).

Notes on the Specific Storage
An interesting implication of Section 2.1 arises when α = 1.The latter refers to systems where the compressibility of porous medium is small compared to the compressibility of grains, such as is the case for unconsolidated materials.Here, the specific storage at constant strain (the inverse of the Biot modulus, Equation 4) reduces to Since the bulk modulus of water is known (K f = 2.2 ⋅ 10 9 Pa), the porosity of the material can also be estimated from the groundwater response to Earth tides.However, for consolidated materials the Biot coefficient may be smaller than one.This can help to constrain the expected values of the specific storage at constant strain.For instance, if the subsurface material is unconsolidated and has realistic porosity values, that is, 0.01 ≤ n ≤ 0.3, then the specific storage at constant strain is constrained to We note that previous studies which estimated the specific storage from the groundwater response to Earth tides have not considered the appropriate context for this property.The result is referred to as "specific storage at constant strain" (S ϵ ) and it can vary significantly from the specific storage generally used in hydrogeology (S s , see Equation 13) (Hantush, 1960).The difference between both coefficients originates from the underlying assumptions.The specific storage at constant strain is defined in conditions in which the volume of the porous frame is maintained constant but the fluid volume is not, which induces changes in the pore volume because fluid has to be accommodated.In contrast, for the specific storage used in hydrogeology the porous frame is allowed to deform in the vertical direction.This is mathematically represented by the second term of Equation 13.Thus, when the subsurface material is much less compressible than the filling fluid and the pore space the second term of Equation 13 tends to zero because no deformation of the frame takes occurs, hence S ≈ S ϵ .Moreover, note that S ϵ ≤ S. Thus, as demonstrated in this study, attention must be paid to the conceptual difference between these two parameters.
Analytical models typically assume that the leaky layers have zero specific storage.Zhu and Wang (2020) numerically investigated the effect of specific storage on Earth tide analysis in leaky layers.The authors concluded that the assumption of leaky layers with zero specific storage may lead to wrong estimations of subsurface properties as the specific storage changes the phase shift.As shown in this work, in unconsolidated systems the specific storage at constant strain depends on porosity only.Therefore, the porosity of the aquifer has to exceed that of the leaky layer based on the results of Zhu and Wang (2020).

Numerical Modeling of the Groundwater Response to Earth Tides
The fluid continuity equation (Equation 8) has been solved in previous studies assuming that the strain term ϵ(t) is known and solely time-dependent with adequate boundary conditions.This equation is an inhomogenous diffusion equation for which the change of volumetric strain is mathematically equivalent to a sink/source in the aquifer storage term.Therefore, changes of strain result in changes of pore pressure in the entire model domain.When mechanical coupling is included, the continuity equation needs to be coupled to the state of stress, hence the strain is tightly coupled to pore pressure.Thus, the strain term in Equation 8 is no longer uniform over the entire model and may vary depending on the amount of change of fluid.For instance, changes in pore pressure (for earth tides within the radius of influence of the well) induces changes in the volumetric strain, which generate drained conditions (Section 3.4).Therefore, assuming that the applied strain is constant within the model domain is inaccurate for our purposes since, as explained before, the strain is function of the pore pressure.Another common way to express the coupling between pore pressure and strain is by rearranging Equations 6 and 7 as BASTIAS ESPEJO ET AL. 10.1029/2022JB024771 14 of 22 The relative movement of celestial bodies in relation to Earth induces variations in the gravitational force which results in small deformation of the Earth's crust.Such deformations are not caused by an applied stress.In continuum mechanics, this problem is known as eigenstrain, and it is very common in heat transport, for example, dilatation caused by heating of materials.The relationship between deformation and eigenstrain can be obtained experimentally leading to a constitutive model (Qu & Cherkaoui, 2006).In this work we apply a simpler approach.As FE implementation typically requires displacement or load as boundary conditions, we set displacement as boundary condition and directly applied the strain obtained from an Earth tide catalog multiplied with the length of the model, that is, While this approach is convenient it has limitations.If the mechanical properties change over the modeling domain (composite material), the displacement will not be uniformly distributed across the domain and therefore the resulting strain will also be non-uniform.This would produce larger displacements in soft layers resulting in higher pore pressure.One way to solve this problem is to assume vertical heterogeneity and to apply the total volumetric strain only at the horizontal boundaries.This would result in a uniform displacement distribution in the horizontal axis and therefore result in an appropriate pore pressure response.We note that the effects of distributed mechanical heterogeneity are not further explored in this work.
Initialization of the numerical model is not trivial since the initial hydrostatic and mechanical states (initial pore pressure and stresses) has to be in equilibrium (Settari & Walters, 2001).This challenge applies in particular for heterogeneous distributions of material properties and transient boundary conditions.Achieving mechanical equilibrium at time t = 0 is difficult and may in most cases require a separate initialization step during the simulation (Chen et al., 2009).We recommend to first simulate steady-state conditions which generates the stress and pore pressure distribution within the modeling domain.
From a numerical point of view, the simulator is setting the force balance as an approximation, that is, ∇σ = 0.
In practice, a non-linear step is finished when the force balance falls below a threshold close to zero but residual errors always remain.Earth tides generate only small changes in pore pressure which are close to the residual error.For example, if the acceptable error is e = 1 Pam −2 , then in our case the area is 50,50,000 m 2 leading to total residuals up to R ≈ 5,050 kPa at the bottom of the model.Since Earth tides generate pore pressure change in smaller magnitude, minimizing the error is an important consideration when modeling.Numerical modeling of Earth tides therefore requires attention to decreasing the tolerance of the numerical solver (e.g., by increasing the number of linear steps), increasing space discretization (e.g., by increasing the size of the Jacobian matrix) or decreasing time discretization (e.g., by increasing the number of time steps).

Are Conditions for the M 2 Earth Tide Drained or Undrained?
When a stress is applied to an undrained subsurface system, the load is shared by the bulk material, the grains and the pore fluid.The balance between these three responses results in instantaneous deformation of the pore space and a change in fluid pressure.If the rate of the applied deformation is slow enough then fluid can flow out of the system which result in a change of the pore pressure.The balance between the rate of Earth tide stress and realistic hydro-geomechanical subsurface properties is rarely known.Moreover, fluid movement (i.e., drained conditions) may be given leading to ξ ≠ 0. Under such conditions the assumptions of the analytical solutions are violated potentially leading to errors when interpreting the groundwater level response to Earth tides.
To assess the conditions under which an undrained response occurs for the M 2 frequency, we numerically simulate a 1D vertical column with depth 0 ≤ z ≤ 5,000 m (Figure 2b) and with a range of realistic hydraulic and geomechanical properties.Equation 7 can be solved for fluid quantity (ξ) assuming the worst scenario (ϵ G = 1 ⋅ 10 −8 corresponding to a low tide amplitude) and α = 1, Figure 5 shows the results of our numerical model which calculates ξ up to 1,000 m depth and for a range of realistic hydraulic properties as well as discrete values of the bulk modulus.As typical Earth tide amplitudes vary between 1 ⋅ 10 −7 ≤ ϵ G ≤ 1 ⋅ 10 as a condition for an undrained response for which the analytical solution is valid, that is, no pore pressure changes occur under this value.This is highlighted in Figure 5 and allows an assessment of the conditions for which the analytical solution should be valid.
Figure 5 shows that undrained conditions are more likely the deeper a system.Further, when the hydraulic conductivity of the leaky layer (k l ) increases, the system behaves more drained.This is expected as the system becomes more permeable and therefore allows flow in response to pressure gradients.This results in fluid movement which causes increased drainage.Similarly, as the specific storage at constant strain increases (rows of plots in Figure 5) the level of drainage decreases.This is because as S ϵ increases the volume of fluid that the system contains due to deformation increases leading to less fluid moving out of the system.This can also be explained using Equation 8 when dividing by S ϵ Equation 36 illustrates that the hydraulic diffusivity of a system  ( decreases with the increase of S ϵ .As the bulk modulus (K) increases, see columns in Figure 5, the system becomes more drained.An explanation for this is that as the filling material becomes stiffer, the mechanical coupling becomes less relevant and the system approaches an incompressible porous skeleton.Under such conditions, only a drained response is allowed and an instantaneous pneumatic response of the system is no longer possible.This can also be explained by revisiting the definition of the Skempton coefficient (B).Assuming α = 1, then which illustrates that when the bulk modulus increases B decreases.This results in a reduction of the overall storativity of the system and consequently also drainage.Another way to understand this result is by considering the coupling of equations.For this simulation we assumed mechanical stress balance as follows Here, Equation 6 must remain constant when the total stress increases (first two terms of Equation 6) thus the amount of fluid leaving the system must increase (third term of Equation 6).
In general, a larger porosity will increase the value of the specific storage at constant strain, which will decrease the level of drainage.Our assessment shows that when the hydraulic conductivity of the leaky layer exceeds k l > 10 −5 ms −1 , this leads to drained conditions and could result in errors when the analytical solution is used to estimate the properties of the aquifer.However, it is worth noting that the level of drainage depends on the geomechanical properties of the system, as well as depth and frequency of the signal.The amplitude of the signal, ϵ 0 , for field measurements, as higher amplitudes will generate higher confined pore pressure and facilitate detection of fluid level changes inside the observation well.
As shown in Figure 5, the level of confinement depends on the hydraulic and geomechanical properties of the subsurface under consideration.Consequently, defining conditions under which an undrained response exists depends on the particular field conditions, for example, depth of the borehole and some knowledge of the subsurface properties.Figure 5 can be used as a preliminary guide for assessing whether or not it is appropriate to apply the analytical solution for interpreting the groundwater level response to Earth tides.
We note that Figure 5 represents the subsurface response to the M 2 frequency.Whether or not a porous medium is drained or undrained depends, among other things, on the frequency of the applied strain.In general, the slower the frequency the deeper the transition between drained and undrained.Consequently, if two Earth tide components were used to estimate properties (Equations 20 and 21) then the observation must be deep enough to ensure undrained conditions for both components.

Robustness and Limits of Analytical Earth Tide Interpretations
Determining subsurface properties from Earth tide responses requires system confinement as a basic condition.To study the effects of a realistic well-aquifer system and the effect of (un)drained conditions, the level of confinement is gradually relaxed in a layered 2D model in this section.The red dots in Figure 6 shows the results from our numerical model (Section 2.5) compared to the analytical solution without vertical leakage (k l = 0 ms −1 ) for a tidal signal with 2 cpd frequency (Section 2.2).The good agreement of amplitude ratios and phase differences verifies our coupled numerical modeling approach.This allows a rigorous hydraulic and geomechanical assessment of how realistic conditions (e.g., subsurface layered heterogeneity) affects the groundwater response to Earth tides.
The effect of the Biot modulus is also shown in Figure 5. Discrete values of α = 0.75, α = 0.5, and α = 0.25 were considered and show the effect over the confined pore pressure generated by a Earth tide deformation.In an undrained system, the Biot modulus represents the ratio of deformation between the porous space and the porous material.Thus it becomes relevant in situations where porous material rearrangement is not possible, such as in Similar results are observed when the specific storage at constant strain of the aquifer is S ϵ = 10 −5 m −1 .When the level of leakage is k l ≤ 5 ⋅ 10 −5 ms −1 the system is undrained or in the transition zone and the numerical results comply with the analytical solution.For lower levels of confinement (i.e., k l > 5 ⋅ 10 −5 ms −1 ), the system becomes drained and the simulation results, once again, differ compared to the analytical solution.
In the particular case when the hydraulic conductivity of the aquifer is k a = 10 −5 ms −1 , the numerical result do not comply with the ones obtained with the analytical solution even under undrained conditions (Figures 7c and 7f).
As the hydraulic conductivity of the aquifer decreases, the finite time to move fluid in or out of the well increases.Hence, it is likely that the fluid velocity is much more influenced by high gradients generated by the drained top boundary rather than by the gradients produce inside the open well under such conditions.
In all three columns of Figure 5, as the level of confinement provided by the aquitard decreases, the simulated results of phase shift tend toward the same value (90° for the simulated system).This means that, as the drainage from the aquitard increases, the effect of the top drained boundary over the pore pressure in the aquifer increases.
The same effect is observable on the amplitude ratio, where the final value of the amplitude ratio is the same in every specific storage at constant strain under study.This effect indicates that the hydraulic conductivity of the aquifer loses relevance as the drainage from the aquitard increases.And, therefore, if the system is draining, at low levels of confinements (k l > 5 ⋅ 10 −5 ms −1 ), the groundwater level measured in the field can potentially result in very similar values regardless of the aquifer hydraulic conductivity.
The effect of drained conditions on the amplitude ratio and phase shift can be better understood when streamlines (i.e., the Darcy velocity field) of the system are plotted.Figure 8 shows streamlines in an area close to the open well when the amplitude of the Earth tide strain is at maximum.This provides understanding of how flow paths change as the level of confinement decreases at fixed aquifer specific storage at constant strain (S ϵ = 10 −6 m −1 ) and hydraulic conductivity (k a = 10 −4 ms −1 ).At the high confinement (k l = 10 −7 ms −1 ) the flow within the aquifer is horizontal and pore pressure gradients are directed toward the well.This complies with the assumption of horizontal flow inherent to the analytical solution.As confinement decreases (i.e., increasing leakage of the aquitard), the velocity field shows increasing flow in the vertical direction through layers which reduces the radius of influence of the well.With the smallest confinement investigated (Figure 6c), vertical flow dominates in the aquifer and almost no horizontal flow is observable.This shows that the pressure wave produced by the open top boundary has strong effects on the amplitude ratio and phase shift at low confinement and dampens the pore pressure signal generated by Earth tides.
Figure 3 can be used in conjunction with the results shown in Figure 7 to assess the potential error due to requiring undrained conditions when the analytical solution is utilized.For example, assuming a typical specific storage at 10.1029/2022JB024771 19 of 22 constant strain of S ϵ = 10 −6 m −1 , leakage of k l = 10 −6 ms −1 and hydraulic conductivity of the aquifer k a = 10 −3 ms −1 , the simulated amplitude ratio is 0.61 and the phase shift is 55.0° (Figures 6a and 6d).Since the phase shift is positive there is leakage from and to the aquifer to the aquifer is occurring (i.e., k l ≠ 0 ms −1 Figures 2c-2h).When the hydraulic conductivity of the aquitard is between 0 ≤ k l ≤ 10 −8 ms −1 (Figures 2c and 2d) the positive phase shifts occur only at low amplitude 0 ≤ A ≤ 0.1.When k l is 10 −4 ms −1 (Figures 2g and 2h) the amplitude is smaller than 0.1 in the studied range.Thus, k l should range between 10 −8 < k l < 10 −4 ms −1 , 10 −4 ≤ k a ≤ 10 −2 ms −1 and the specific storage at constant strain between 10 −5 ≤ S ϵ ≤ 10 −4 m −1 (Figures 2e and 2f).
Our results show that a high specific storage at constant strain (S ϵ ≥ 10 −6 m −1 ) in combination with a high confinement (k l ≤ 5 ⋅ 10 −5 ms −1 ) and hydraulic conductivity of the aquifer (k a ≥ 10 −4 ms −1 ) allow application of the analytical solution.Application to real world system further requires a high contrast in hydraulic conductivity between the layers  (  −1  ≥ 10 3 ) with specific storage values that are typical (≈10 −6 m −1 ).In reality, the confined pore pressure is damped by the movement of fluid and fully undrained conditions may be rare.Any a priori knowledge of the formation (e.g., thicknesses and hydraulic properties) is key in the assessment of Earth tidal analysis, not only to have a good approximation when inverting equations, but also to approximate the level of drainage and therefore assess potential errors when the analytical solutions are utilized.
For unconsolidated systems the soil matrix is more compressible than the grains (i.e., an unconsolidated subsurface) which leads to α = 1.Moreover, hydraulic conductivity and porosity for what can be considered an aquifer varies between 10 −2 ms −1 ≤ k a ≤ 10 −4 ms −1 and 0.2 ≤ n ≤ 0.3, respectively.This means that the specific storage at constant strain varies between 9.1 ⋅ 10 −7 m −1 ≤ S ϵ ≤ 1.4 ⋅ 10 −6 m −1 (with the bulk modulus of water, K f = 2.2 ⋅ 10 −9 Pa).Considering these ranges and given sufficiently high confinement between layers, application of the analytical solution to well fluid levels is valid.For example, this would be the case for hydraulic properties of sands and gravels overlain by clays or silts.
For consolidated systems, Earth tidal analysis poses a challenge as the Biot coefficient generally is α < 1.In order to use the groundwater response to Earth tides, the Biot coefficient has to be known as it directly attenuates the pore pressure response to strain (Equation 7).Although some values of the Biot coefficient have been reported for different rock types varying from 0.1 to 1 (Cheng, 2016), real world measurements are difficult to find in the literature (Cosenza et al., 2002;C. Wang & Zeng, 2011).Our work shows that the Biot coefficient requires estimation when the groundwater response to Earth tides is quantitatively evaluated for wells screened in consolidated systems.Hence, for real systems, this leads to the following trade-off: As deeper wells are more likely to contain Earth tide influences because undrained conditions exists, but they are also more likely consolidated, in which case an estimate of the Biot modulus is required.Overall, our results show that a presence of Earth tide components in wells that are screened in deep and unconsolidated systems are likely to have undrained conditions and are therefore suitable for interpretation.

Conclusions
The amplitude and phase of the groundwater response to harmonic Earth tide components can be used to estimate hydraulic conductivity and specific storage values of aquifer systems.However, this approach is based on simplified analytical solutions to the groundwater flow equation, which has various assumptions that have not been tested yet.To assess the effect of such assumptions, we present a numerical method to simulate the groundwater response to Earth tides by coupling compressible flow to geomechanics.We demonstrated that this can be solved numerically using the Multi Object Oriented Simulation Environment (MOOSE) and verify this using a simplified analytical solution of the groundwater flow equation.We further use simulations to assess the conditions of validity for simplified analytical solutions when estimating hydraulic properties from the groundwater response to Earth tides.
By first focusing on the aquitard layer, we assess the subsurface response type, that is, drained or undrained conditions, to the dominant harmonic Earth tide component at M 2 with frequency of 1.93227 cpd for depths up to 5 km and a range of hydraulic conductivities.Based on typical Earth tide strains, we define that undrained conditions exist when the incremental of fluid content is smaller than 5 ⋅ 10 −11 for which the groundwater equation and associated analytical solution should be valid.Our results show that this is the case for specific storage at constant strain larger than 1 ⋅ 10 −6 m −1 and depths higher than 50 m for low conductivity systems (k a < 10 −7 ms −1 ) and depths up to 1 km for high conductivity systems (k a ≥ 10 −3 ms −1 ).
We revisited previously interpretations based on analytical solutions and showed that the specific storage has been often misinterpreted.Moreover, only an approximation of the exact solution of the non-linear analytical Figure1.Representation of the pore pressure change and well water level in a semi-confined aquifer due to Earth tide strains.Earth tides induce subsurface stress which results in strains that generate changes in the pore water pressure.This leads to pressure gradients between the subsurface and observation wells that cause water movement in or out of the well.The ratio between the Earth tide strain and the well water level response can be expressed as amplitude ratio whereas the time difference between both signals as phase shift.(a) Aquifer and well pressure are in hydro-static equilibrium.(b) Earth tidal strain changes confined pore pressure.(c) Pressure gradients move fluid to the observation well.(d) Confined pore pressure changes due to Earth tidal forcing and well water level changes over time.

Figure 2 .
Figure 2. Overview of the conceptual models used in this work.(a) Analytical model ofWang et al. (2018) when no external force is applied (b) Analytical model ofWang et al. (2018) when the confined pore pressure generated at the far field is generated and fluid flow toward the well.(c) A 1D column of the subsurface representative of the aquitard to assess the type of response.At t = 0 the column is equilibrium, at t > 0 a harmonic strain is applied which results in fluid movement in and out of the column.(d) A 2D model of the aquifer bounded by a aquitard and connected to a well.Earth tide strain is applied which moves fluid toward a well that is numerically modeled.

Figure 3 .
Figure 3. Amplitude and phase shift response of the analytical solution presented in Wang et al. (2018) for a realistic range of hydraulic conductivity and specific storage at constant strain values.(a and b) are representative of zero leakage through the aquitard corresponding to Hsieh et al. (1987).(c-f) consider distinct and increasing aquitard hydraulic conductivity values.The harmonic signal frequency is 2 cpd, the well and screen radius are 0.2 m, the screen length is 1 m and the aquitard has 100 m thickness.

Figure 4 .
Figure 4. Color map exploring the solution space, that is, the variability of parameters as a function of the initial guess, of the under-determined problem by(Wang et al., 2018).Estimating three hydraulic properties out of two measured parameters: (a) aquifer hydraulic conductivity, (b) specific storage at constant strain, (c) aquitard hydraulic conductivity.Note that each color scale has a different range.Blue indicates less variability, whereas red means more variability of the results.

Figure 5 .
Figure5.Change of fluid content ξ over depth and aquitard hydraulic conductivity for a 1D column (Figure2b) and Earth tide forcing with M 2 frequency.Rows correspond to different values of specific storage whereas columns are representative for different bulk moduli, for example, clay (a, d, g), sand (b, e, h) and hard rock (c, f, i).Values of ξ can be used to infer the depth at which the system response is undrained, that is, where application of the analytical solution (Equations 18 and 19) is valid.Value ranges of validity are delineated by the dashed line.

Figure 6 .
Figure6.Verification of the 2D numerical model againstHsieh et al. (1987) for a harmonic forcing signal with 2 cpd frequency.Here, the simulation of the hydraulic conductivity of the leaky layer was set to zero.The figure also shows the effect of the Biot modulus on the confined pore pressure.

Figure 8 .
Figure 8. Streamlines show the velocity field toward the observation well during maximum Earth tide strain for three different aquitards with varying hydraulic conductivities: (a) low, (b) medium, (c) high.This illustrates that the flow direction changes from horizontal to vertical as the leakage increases.

Table 1
considers vertical flow in a Overview of the Hydraulic and Mechanical Properties of the Numerical Model Note.The subscript l and a refer to aquitard and aquifer, respectively.