Articles | Volume 28, issue 4
Technical note
29 Feb 2024
Technical note |  | 29 Feb 2024

Technical note: Analytical solution for well water response to Earth tides in leaky aquifers with storage and compressibility in the aquitard

Rémi Valois, Agnès Rivière, Jean-Michel Vouillamoz, and Gabriel C. Rau

In recent years, there has been a growing interest in utilizing the groundwater response to Earth tides as a means of estimating subsurface properties. However, existing analytical models have been insufficient in accurately capturing realistic physical conditions. This study presents a new analytical solution to calculate the groundwater response to Earth tide strains, including storage and compressibility of the aquitard, borehole storage, and skin effects. We investigate the effects of aquifer and aquitard parameters on the well water response to Earth tides at two dominant frequencies (O1 and M2) and compare our results with hydraulic parameters obtained from a pumping test. Inversion of the six hydro-geomechanical parameters from amplitude response and phase shift in both semi-diurnal and diurnal tides provides relevant information about aquifer transmissivity, storativity, well skin effect, aquitard hydraulic conductivity, and diffusivity. The new model is able to reproduce previously unexplained observations of the amplitude and frequency responses. We emphasize the usefulness in developing a relevant methodology to use the groundwater response to natural drivers in order to characterize hydrogeological systems.

1 Introduction

Aquifer properties play a vital role in managing groundwater resources, particularly amid increasing anthropogenic groundwater use and the impact of climate change. While pump testing can be costly, there exists a cost-effective alternative for assessing aquifer hydraulic properties – analysing the groundwater response to Earth tides or atmospheric tides (McMillan et al., 2019; Valois et al., 2023). Observations of variations in groundwater level due to tidal fluctuations date back to the works of Klönne (1880), Meinzer (1939), and Young (1913). However, it was only later that hydro-geomechanical models were employed to elucidate these variations (Bredehoeft, 1967; Hsieh et al., 1987; Roeloffs, 1996; Wang, 2000; Cutillo and Bredehoeft, 2011; Kitagawa et al., 2011; Lai et al., 2013; Wang et al., 2018). This progression in understanding offers a valuable opportunity to evaluate aquifer hydraulic properties through the response of groundwater to Earth tide strain fluctuations within the aquifer system.

Hsieh and Bredehoeft (1987) introduced the horizontal-flow model, focusing on confined conditions influenced by tidal strain within the aquifer. Conversely, Roeloffs (1996) and Wang (2000) explored interactions within vertical flow under tidal fluctuations. Wang et al. (2018) expanded on these studies by incorporating a flow from an upper aquitard, albeit assuming negligible storage within it. Later, Gao et al. (2020) extended these models to include borehole skin effects. A definition of the skin effect can be found in Wen et al. (2011). Thomas et al. (2023) developed an Earth-tide–groundwater (ET-GW) model incorporating storage and strain response in the aquitard. They applied their model to a specific site to evaluate transmissivity variations and validated it using pumping tests.

Numerous studies have investigated aquifer hydromechanical properties by analysing groundwater (GW) level variations induced by Earth tides, employing the models mentioned in the previous literature (Narasimhan et al., 1984; Merritt, 2004; Fuentes-Arreazola et al., 2018; Zhang et al., 2019a; Shen et al., 2020). Some studies focused also on the tidal response in fractured rock aquifers (Carr and Van Der Kamp, 1969; Bower, 1983; Burbey et al., 2012; Rahi and Halihan, 2013; Sedghi and Zhan, 2016). However, only a limited number of validations have been conducted, which involve comparing the results with robust hydraulic assessments, such as hydraulic conductivity derived from slug testing (Zhang et al., 2019b) or specific storage and transmissivity characterizations obtained through long-term pumping tests (Allègre et al., 2016; Valois et al., 2022). The current evaluations predominantly focus on purely confined conditions, leaving a gap in knowledge regarding tidally induced GW responses in aquifers under semi-confined conditions.

As far as the authors are aware, the publications by Sun et al. (2020), Valois et al. (2022), and Thomas et al. (2023) are the sole references addressing a comparison for a leaky aquifer. Sun et al. (2020) found significant discrepancies between transmissivities obtained from Earth tide fluctuations and those derived from slug or pump tests. However, it is worth noting that the comparison may be subject to discussion, as the authors employed a leaky aquifer model for analysing tidally induced fluctuations, whereas they used a confined aquifer model for conducting slug and pump tests. In the study conducted by Valois et al. (2022), the existing ET-GW models, as described earlier, were unable to reproduce a low ratio of semi-diurnal to diurnal amplitude with positive phase shifts in conjunction with pumping test transmissivity data. This discrepancy highlights the complexity and challenges in modelling tidally induced groundwater responses in leaky aquifers and the need for further investigation in this area.

We note that our previous attempts to model the observed substantial amplitude decrease from the diurnal tide (O1) to the semi-diurnal tide (M2) frequency, combined with phase shifts close to zero, proved unsuccessful when using Earth tide models found in the literature. None of the existing models could provide satisfactory results. The pursuit of an explanation led to the realization that analytical solutions with more realistic assumptions are required. For example, aquifers are widely recognized to be influenced by aquitards, which often consist of highly porous and compressible clay materials, contributing significant amounts of stored water to the aquifer (Moench, 1985). Moreover, these aquitards are also impacted by Earth tide strains (Bastias et al., 2022).

Our first objective is to develop an analytical solution considering storage and tidal strain in the aquitard. Unlike Thomas et al. (2023), our model incorporates borehole skin effects and allows for a fixed hydraulic head at the top of the aquitard, broadening its applicability to a broader range of hydrogeological conditions. The second motivation of our study is to develop a model that better fits the observed results by considering aquitard storage, as evident in the pumping tests. Third, we compare the results obtained from our new ET-GW model with those derived from a pumping test in leaky aquifers with storage in the aquitard. Fourth, since most publications have predominantly focused on assessing hydraulic properties using the semi-diurnal tide (M2), our third motivation is to demonstrate the potential of using diurnal tides (O1) in combination with the semi-diurnal ones (M2 or N2) to provide a more comprehensive characterization of aquifer and aquitard hydromechanical properties. Our new development offers the potential to enhance hydraulic and geomechanical subsurface characterization by employing a more realistic model for the groundwater response to natural forces.

2 Groundwater response to Earth tides in a leaky aquifer with aquitard storage and strain

Hantush (1960) pioneered the modelling of aquitard storage by modifying the leaky aquifer theory to account for storage in the aquitard. In our study, we consider a semi-confined configuration (Fig. 1) where the target aquifer is overlain by an aquitard that allows for storage, strain, and vertical flux. Both layers are assumed to be slightly compressible, spatially homogenous, and infinite laterally and to have constant thickness. Building upon the work of Wang et al. (2018), our research incorporates Earth tide (ET) fluctuations into the leaky aquifer equations proposed by Moench (1985). Additionally, we incorporate the skin effect, as described by Gao et al. (2020). Two-dimensional cylindrical coordinates are used because of the radial symmetry caused by the well surrounded by hydrogeological material affected by tidal strain.

Groundwater flow and storage in an aquifer overlain by an aquitard can be described as (De Marsily, 1986):


Here, h (m) and h (m) are the hydraulic heads in the aquifer and the aquitard respectively; hj (m) is the fixed hydraulic head at the top of the aquitard; r (m) is the radial distance from the studied well; T (m2 s−1) and S are the aquifer transmissivity and storativity; B, B, Ku (Pa), and Ku (Pa) are Skempton's coefficient and the undrained bulk modulus of the aquifer and aquitard; ρ (kg m−3) and g (m s−1) are the water density and gravity constant; ε is the volumetric Earth tide strain; K (m s−1) is the aquitard hydraulic conductivity; Ss (m−1) is the specific storage of the aquitard; and D (m2 s−1) is aquitard hydraulic diffusivity (K/Ss ratio). Any natural regional groundwater flow is considered negligible.

Borehole drilling causes a zone of damage with a radius rs (see Fig. 1) that is responsible for the skin effect (Van Everdingen, 1953). A negative skin can be caused by a greater hydraulic conductivity around the well because of the material damaged by the drilling, while a positive skin can be associated by porosity clogging caused by the drilling mud. This is reflected in the well's pressure head Δhs. The skin factor (sk) can be defined as

(3) sk = Δ h s r h r r = r w .

Following above assumptions, the boundary conditions are


Here, t is the time (s); rw and rc are the radius of the well-screened portion and the radius of well casing in which the water level fluctuates; zi is the aquifer–aquitard interface elevation (see Fig. 1) and hw is the hydraulic head at rw.

Figure 1Semi-confined system with a compressible aquitard with storage.


Following Hsieh et al. (1987) and Wang et al. (2018), complex numbers were used to facilitate harmonic model development, and the solution is obtained by first solving Eq. (2) in the aquitard and then deriving the head response in the aquifer far away from the well (h), which is independent of the radial distance. h is the aquifer response to the tidal harmonic sources far from the well. Thus, h is the aquifer hydraulic head response without any disturbance from a well–aquifer system. Then, the well effect on the aquifer response is considered by using a flux condition at the well that accounts for wellbore storage. Since h, h, ε,hw, and h are all periodic functions, they can be expressed as


Here, i=-1; ε0 () is the ET strain amplitude and ω (s−1) is the angular frequency. In this case, Eq. (2) becomes

(14) 2 h 0 z 2 = 1 D i ω h 0 - i ω B K u ρ g ε 0 .

According to Wang (2000) and Roeloffs (1996) and as detailed in Appendix A, the solution of Eq. (2) is

(15) h 0 = A 1 e 1 + i δ ( z - z i ) + A 2 e - 1 + i δ ( z - z i ) + B K u ρ g ε 0 ,

where δ=2Dω1/2. Thus, at the interface between the aquifer and the aquitard (z=zi), we have pressure continuity as hzi,t=h0eiωt=h(t) which leads to

(16) h z z = z i = 1 + i δ h 0 - B K u ρ g ε 0 tan h ( 1 + i δ z i ) - h j - B K u ρ g ε 0 sin h ( 1 + i δ z i ) e i ω t .

Equation (16) is in agreement with Butler and Tsou (2003), where leakage is shown to be a scale-invariant phenomenon.

Equation (1) can be solved far away from the well using h, which is independent of the radial distance from the well, and by using the source term from h as follows:


The disturbance in water level due to the well can be expressed as

(19) s r , t = h r , t - h ( t ) .

By expressing Eq. (1) with the sum of s and h (Eq. 19) and using Eq. (16) to express the leaky component and using Eq. (17) to remove h, it follows that

(20) T 2 s r 2 + 1 r s r = S s t - K 1 + i δ s 1 tan h ( 1 + i δ z i ) ,

with the following boundary conditions:


The solution of this differential equation is sr=CII0βr+CKK0βr (Wang et al., 2018), where I0 and K0 are the modified Bessel functions of the first and second kind and the zeroth order respectively. Further,

(24) β = i ω S T - K T 1 + i δ 1 tan h ( 1 + i δ z i ) 1 / 2 .

The boundary conditions lead to CI=0 and CK=-iωrc2hw,02TβrwK1βrw because dK0(r)dr=-K1(r). Therefore, the final solution for the well water level is

(25) h w , 0 = B K u ρ g ε 0 S i ω + K 1 + i δ B K u B K u - 1 tan h ( 1 + i δ z i ) - h j ρ g B K u ε 0 - 1 sin h ( 1 + i δ z i ) σ S i ω - K 1 + i δ 1 tan h ( 1 + i δ z i ) ,


(26) σ = 1 + i ω r c 2 K 0 β r w 2 T β r w K 1 β r w + i ω r c 2 2 T sk .

By assuming hj=0 (i.e. the hydraulic head at the top of the aquitard corresponds to the unsaturated–saturated interface at z= 0), Eq. (25) can be reorganized to

(27) h w , 0 = B K u ρ g ε 0 S i ω + K 1 + i δ R K u B 1 - cos h ( 1 + i δ z i ) sin h ( 1 + i δ z i ) σ S i ω - K 1 + i δ 1 tan h ( 1 + i δ z i ) ,


(28) R K u B = K u B K u B .

By disregarding the BKuρg product which only controls the amplitude, the solution has six independent parameters which are T, S, K, D, sk, and the RKuB ratio.

Let us now define the amplitude response (or amplitude ratio), A, and phase shift,α, of the GW response to ET fluctuations:


Figure 2 shows the amplitude response and phase shift as a function of frequency using our new solution in comparison to key models reported in the literature. Aquitard parameters were set according to Batlle-Aguilar et al. (2016), while aquifer parameters were chosen according to the field application below. We validate the solution using a very low aquitard conductivity (10−14 m s−1), so that we can compare it to the horizontal flux with wellbore storage model (Hsieh et al., 1987). It shows a perfect match.

Figure 2Frequency variation in amplitude response and phase shift in the groundwater response to Earth tides. The transmissivity (T) is 10−5 m2 s−1, storativity (S) is 10−4, hydraulic conductivity of the aquitard (K) is 10−8 m s−1, aquitard hydraulic diffusivity (D) is 10−4 m2 s−1, the skin factor (sk) is 0, RKuB is 1.4, the well casing radius (rc) and screen radius (rw) are 6.03 cm. b was set to 5 m. Screen depth (z) was set to 23 m for the vertical-flow model of Wang (2000).


Because both horizontal and vertical flux models are associated with opposite phase shift signs (Fig. 2b), the latter can offer valuable insights for model selection (Allègre et al., 2014). Positive phase shifts in the vertical flux model are related to an increasing amplitude ratio with frequency, whereas the wellbore storage model exhibits the opposite behaviour. Wang et al. (2018) developed a leaky model capable of demonstrating both positive and negative phase shifts, where positive phase shifts correspond to an increasing amplitude ratio with frequencies and negative phase shifts are linked to a decreasing amplitude ratio.

Our new model showcases positive or negative phase shifts with an increasing or decreasing amplitude ratio over frequency, even allowing for amplitude ratios greater than 1. Notably, Wang (2000) observed a similar characteristic in the vertical flux model, with an amplitude just above 1 (1.06) for very specific conditions. At high frequencies, our model displays amplitude ratios and phase shifts similar to those of the leaky and wellbore storage models, reflecting the attenuation of high-frequency pore pressure fluctuations in the aquifer by well water.

At low frequencies (Fig. 2), the purely confined model exhibits a constant phase shift (0°) and amplitude ratio (Eq. 1). This constant behaviour is the signature of the absence of well impact on groundwater level fluctuations and the absence of phase shift between the Earth tide strain and the aquifer level fluctuations. It means that the groundwater fluctuations in the aquifer are the same as the groundwater fluctuations in the well (absence of amplification/attenuation and phase shifts) and that there is no phase shift between the strain and the water pressure variations inside the aquifer. For the “Leaky & Storage (present study)” model (Fig. 2), the leaky conditions do provoke a phase shift and an amplitude modification as compared to being purely confined. Such values of phase shift and amplitude modification do vary with frequency.

We explored the parameter space by focussing on the frequency-dependent amplitude response and phase shift responses for different sets of parameter values. The reference parameter set is the one described above. T, RKuB, and zi have a large impact on model shapes. As also observed by Hsieh et al. (1987), S does not have a major impact on the results (Fig. 3c and d). The skin effect does not play a large role in the useful frequency band (about 1 to 2 cpd, cycles per day) for amplitudes, but its influence is larger for the phase shifts when compared to the reference parameter set. K does not significantly influence the results with respect to the reference parameter set used in the study (Fig. 3e to f). K has the opposite role of S, and they appear to compensate each others' effects because of their respective role in Eq. (27).

Figure 3Illustration of the amplitude response (a, c, e, g, i, k, m) and phase shift (b, d, f, h, j, l, n) as a function of frequency for various parameters compared to the reference parameter set.


Figure 4 shows the impacts of considering leaky and skin effects for a realistic parameter set. Purely confined conditions (no leaky aquitard) do not create a phase shift between the volumetric strain in the aquifer (Fig. 4a) and the associated hydraulic head variation (Fig. 4b), while a compressible and leaky aquitard (Fig. 4b) or skin effects around the well (Fig. 4c) could involve positive phase shifts.

Figure 4Example of the volumetric strain time series generated by the M2 Earth tide in (a), which creates aquifer hydraulic head variations in (b), resulting in well water level variations in (c). The transmissivity (T) is 10−6 m2 s−1, storativity (S) is 7 × 10−4, hydraulic conductivity of the aquitard (K) is 10−6 m s−1, aquitard hydraulic diffusivity (D) is 10−4 m2 s−1, the skin factor (sk) is 5 m, zi is 10 m, RKuB is 0.3, and the well casing radius (rc) and screen radius (rw) are 6.03 cm. B is set to 0.8 and Ku to 10 GPa.


3 Application of the new model to a groundwater monitoring dataset from Cambodia

3.1 Field site and previous results

The field site in northwest Cambodia comprises three boreholes drilled into the subsurface, consisting of mudstone, claystone, siltstone, and sandstone. Time series and pumping test results have been reported in Valois et al. (2022), while details of the lithology can be found in Vouillamoz et al. (2012, 2016) and Valois et al. (2017, 2018, 2022). Pumping from the aquifer is limited by a very low specific yield, attributed to the presence of fine deposits such as clay and mudstone (Vouillamoz et al., 2012; Valois et al., 2018).

The boreholes were drilled to a depth of 31 m with a radius of 0.1524 m and equipped with 0.1016 m. There is PVC casing from top to bottom, featuring a 9 m long screen at the hole's base. The aquifer is situated within a hard rock medium, comprising either claystone or sandstone, located beneath a 10 m thick clay layer.

For the pumping tests, the wells were pumped for 3 d, and water levels were allowed to recover for 4 d in two observation wells. The interpretation of the pumping tests utilized AQTESOLV/Pro v4.5 software, employing the leaky aquifer with aquitard storage model (Moench, 1985) or a 3D flow using the generalized radial flow model (Barker, 1988). The selected solutions, compared to other models (Theiss, Hantush, without aquitard storage), demonstrated the best fit with a root mean square error (RMSE) of 0.02 m for Cambodia.

3.2 Well sensitivities and phase shifts to Earth tides

Between 2010 and 2015, well water levels were measured at 20, 40, or 60 min intervals using absolute pressure sensors (Diver data loggers, Eijkelkamp Soil & Water, the Netherlands). To compensate for barometric pressure (BP) effects, data from a barometer located a few kilometres away from the field site were utilized (Eijkelkamp Soil and Water, the Netherlands). A zero-phase Butterworth filter was employed to eliminate low-frequency content (periods longer than 10 d) from both groundwater (GW) and BP data.

For each site's geolocation (latitude, longitude, and height), ET strain time series were computed at 20 min intervals using the SPOTL software (Agnew, 2012). The time series were then modelled using harmonic least squares (HALS; Schweizer et al., 2021) with eight frequencies corresponding to the major tides (Table 1) following Merritt's description (2004). HALS provides amplitude and phase estimations for each tidal component and record.

Table 1Dominant tidal components that are generally found in groundwater measurements (adapted from MacMillan et al., 2019).

Download Print Version | Download XLSX

The results obtained from HALS were utilized to calculate the amplitude response and phase shift between GW and Earth tide (ET) for each tidal component. These amplitude responses are commonly known as “well sensitivities” to Earth tide strains (Rojstaczer and Agnew, 1989) and are summarized in Fig. 5 alongside the corresponding phase shifts.

Figure 5Amplitude responses and phase shifts as a function of frequency for the Cambodian site. S1 and S2 were excluded as they are not only generated by Earth tides. The circle size is proportional to the amplitude in the well water levels.


The well sensitivities to tides exhibit a frequency-dependent behaviour, resulting in similar values for neighbouring frequencies and a generally decreasing magnitude. The amplitudes of M2 and N2 are relatively straightforward to assess due to their significant magnitudes (11.2 and 2.1 mm respectively), and their amplitude responses and phase shifts are highly similar. The phase shifts for the tides of interest (O1, N2, and M2) are positive. However, the signs of the amplitudes for the other tides can be attributed to their low-amplitude responses, which are challenging to characterize using HALS.

3.3 Fitting the M2/O1 amplitude response ratio and phase shifts

The analysis is restricted to two types of tides: the semi-diurnal and diurnal tides. This limitation arises because the magnitude of Earth-tide-induced well water levels is significantly damped for higher frequencies, making it difficult to discern and analyse tides beyond these two types. Here, we use amplitude responses (AM2, AO1) and phase shifts (αM2, αO1) to estimate hydraulic subsurface properties. The N2 tide was not used since its response may be too similar to M2 and does not help with constraining the model. Amplitudes are influenced by geomechanical parameters (BKu) which are generally not considered in classical hydrogeology. Valois et al. (2022) previously illustrated that the M2-to-O1 amplitude response ratio can be computed because it is not directly multiplied by BKu and because it provides useful information about model choice. This leads to a system of three equations and six parameters (T, S, K, D, skin, and RKuB) by using the simplified model in Eq. (27) when the geometry of the well and the aquitard–aquifer system is known:


Table 2 displays the data to be fitted using the three equations above.

Table 2Data to be fitted using the ET-GW model.

Download Print Version | Download XLSX

A systematic exploration of the entire parameter space without any constraints other than the well and aquifer geometry was carried out. Hydraulic and geomechanical property ranges are chosen according to the literature, i.e. De Marsily (1986) and Domenico and Schwartz (1998). In order to assess the goodness of fit with the three observed parameters (Eqs. 31 to 33), the objective chi-square function is defined below:

(34) χ 2 = 1 / N i = 1 N ( Obs i - Mod i ) Error i 2 ,

where N is the number of observed parameters (3 here), Obsi, Modi, and Errori are the observed parameter, modelled parameter, and their errors respectively. Thus, this objective function takes into account errors in the observed parameters (De Pasquale et al., 2022). They were set to 0.1, 0.5° and 0.2 for αM2, αO1, and AM2AO1 respectively, according to Valois et al. (2022).

The model allows us to fit both O1 and M2 positive phases and the low M2-to-O1 amplitude ratio (misfit close to 0 in Fig. 6), whereas the model of Gao et al. (2020) cannot do this (misfit above 1; see Appendix B and Valois et al., 2020). The T value is in good agreement with the pumping test range (Fig. 6a). S is 0.5 orders of magnitude below the pumping test range (Fig 6b), whereas the storativity best fit for Gao et al. (2020) is 2 orders of magnitude above (Appendix B). The skin effect also shows acceptable values as compared to the pumping test (Fig. 6f). The parameter exploration shows best fits for K and D, whereas it is difficult to identify a clear best fit for the RKuB parameter. The values are within the expected range for the hydrogeological configuration: the mudstone aquitard has a lower hydraulic conductivity (10−8 m s−1) than the underlying claystone aquifer (coarser grain size than the aquitard, with a K value of about 10−7 m s−1 for an aquifer thickness of 22 m) and a diffusivity of about 10−4 m2 s−1. This is in agreement with the aquitard classification of Pacheco (2013).

Figure 6Results of full parameter exploration using the leaky-with-storage aquifer model for the Cambodian case study. Red circles represent the best fit.


4 Discussion

4.1 Uncertainties and discrepancies

There are several sources of uncertainty which originate from measurement and their propagation as well as uncertainties introduced by model assumptions. We believe that uncertainties linked to pressure sensor resolution (0.2 mm) and time resolution (20 min) as well as the HALS decomposition were too low to be worth considering, at least for the semi-diurnal tides. This can be deduced from the nearly identical amplitude responses at M2 and N2 at our field site. Because those responses are indeed identical, it means that errors in the raw dataset did not influence the response characterization. We note that amplitude responses and phase shifts show larger discrepancies for the diurnal tides. This could be linked to overall lower amplitudes which are generally more difficult to characterize. We therefore conclude that errors arising from uncertainties are negligible compared to the uncertainty introduced by model assumptions, in agreement with Sun et al. (2020).

Discrepancies between hydro-geomechanical properties derived from the groundwater response to Earth tides (termed “passive” and assuming a compressible matrix) and hydraulic testing (e.g. slug, pump and lab testing, termed “active” and generally assuming an incompressible matrix) have been reported in the literature and have not been appropriately reconciled. By fitting the amplitude response ratio and phase shifts (Sect. 3.3), a T value discrepancy of 1 order of magnitude can be observed between both approaches. We hypothesize that this is caused by parameter anisotropy.

Zhang et al. (2019b) pointed out differences in hydraulic conductivities of more than 1 order of magnitude between ET analysis and slug tests and attributed this to differences in the investigated scale. Allègre et al. (2016) reported much higher values of storativity derived from pumping test as compared to ET when using the vertical-flow model. Sun et al. (2020) showed that T values are frequency-dependent with differences of several orders of magnitude when comparing co-seismic, ET, and slug or pump test methods. The discrepancies can be explained by the different conceptual models used in the active (based on perfectly confined conditions) and passive methods (based on leaky conditions) or by the frequency dependency of hydraulic parameters. The literature illustrates that transmissivity, hydraulic conductivity, or specific storage can indeed vary depending on the frequency of the forcing (e.g. Cartwright et al., 2005; Renner and Messar, 2006; Guiltinan and Becker, 2015; Rabinovich et al., 2015). This demonstrates the need for attention when assessing hydraulic parameters using passive methods for semi-confined conditions. We specifically emphasize the need to use the same conceptual model (i.e. confined, leaky with or without storage, vertical flow) when comparing active and passive methods, as well as the need for preliminary hydrogeological knowledge of both the aquifer system (i.e. presence of an aquitard with or without storage) (Bastias et al., 2022) and the borehole skin effect.

4.2 The use of the leaky model with aquitard storage

Our new analytical solution describing the well water level response to harmonic Earth tide strains contains at least six hydro-geomechanical parameters that could be derived from only three features, e.g. the M2-to-O1 amplitude response ratio and M2 and O1 phase shifts. Applying this model to real-world cases to derive properties from amplitude responses and phase shifts provides relevant information on T, S, D, K, and the skin effect, but it is prone to non-uniqueness. Thus, a priori information may be needed depending on the capacity of the inverse problem to fit observed data (phases shifts and amplitude ratio). In our case study, parameter assessment would benefit from prior information on S (or K) and RKuB.

The model presented in this study can be useful when the hydrogeological configuration involves storage in the aquitard with fixed head (i.e. Dirichlet) boundary conditions and for cases where phase shifts and amplitude ratio do exemplify a specific pattern. For example, when compared to Gao et al. (2020) and using the parameterization of the present study (Fig. 7), our solution is able to model lower a M2-to-O1 amplitude ratio, lower phases, and a higher amplitude ratio for phases close to 0.

Figure 7Outputs of the models using the parameterization of the study (rc = rw = 6.08 cm, zi=-10 m).


While the amplitudes are controlled by the product of the Skempton coefficient and the undrained bulk modulus, these mechanical parameters also affect phase shifts. Therefore, further investigations are needed to assess these influences using other methods or to link them empirically with the hydraulic parameters. This is crucial to enhance confidence in utilizing the groundwater response to Earth tides as a valuable tool for better understanding and characterizing groundwater resources.

5 Conclusions

We have developed a new analytical solution for the well water level response to Earth tide strains. This solution considers a previously unprecedented physical reality, specifically, a leaky aquifer with aquitard storage, subject to Dirichlet boundary conditions under tidal strain. Additionally, our model considers the influence of borehole storage and skin effects, further improving the accuracy and comprehensiveness of the analysis. This model extends upon previous models and allows advanced characterization of the subsurface using the groundwater response to natural forces. The new model overcomes previous limitations, for example it explains very low M2-to-O1 amplitude ratios as well as a large phase shift difference between M2 and O1 tides. The model relies on six combinations of hydro-geomechanical parameters. In this study, we assess the most sensitive parameters to be the transmissivity, the well skin effect, the aquitard-to-aquifer mechanical parameter ratio (BKu/BKu), and the aquitard diffusivity and aquitard-conductivity-to-aquifer-storativity ratio.

We apply our new model to a groundwater monitoring dataset from Cambodia and compare the results with pumping tests undertaken in the same formation. We used the diurnal (O1) and semi-diurnal (M2) tides to better constrain the model. Results illustrate significant insight into subsurface properties. For example, we derive relevant information about T, S, D,K, and the skin effect, when compared to the pumping test results. Overall, our new model can be used to shed light on previously inexplicable well water level behaviour and can be paired with other investigation methods to enhance the understanding of subsurface processes.

Appendix A: The analytical solution in the aquitard

To solve Eq. (14) with the boundary conditions in Eqs. (7) and (8), we define

(A1) h ^ 0 = h 0 - B K u ρ g ε 0 .

Thus, the equation system becomes


The solution h^0 is of the following form:

(A5) h ^ 0 = A 1 e 1 + i δ ( z - z i ) + A 2 e - 1 + i δ ( z - z i ) .

It yields



(A8) h 0 = A 1 e 1 + i δ ( z - z i ) + A 2 e - 1 + i δ ( z - z i ) + B K u ρ g ε 0 .

Appendix B: Additional information on parameter exploration

Figure B1 shows the misfit functions using the model developed in this study and the model of Gao et al. (2020). Misfits are clearly higher for the older model, which does not consider storage and tidal response in the aquitard. The storativity best fit using the model of Gao et al. (2020) failed to reproduce pumping test values. Nevertheless, transmissivity and skin estimates fall within the pumping test range.

Figure B1Comparison of misfit function for the present model and the one of Gao et al. (2020) for the Cambodian case study. Only the first hundred best fits were plotted.


Code and data availability

The field data and code are hosted on the GitHub platform at (Valois, 2024).

Author contributions

RV did the field survey and analysis and coordinated the paper writing. AR verified the model theory and participated in the writing. JMV did the field work and participated in the writing. GCR participated in the data analysis and paper writing.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


We thank Sang Sokheng and Phan Sophoeun for their efficient assistance in field work.

Financial support

This work has been carried out in the framework of the Institut de Recherche pour le Développement and the French Red Cross collaborative project 39842A1 – 1R012-RHYD, with financial support from the European Commission (grant nos. DIPECHO SEA ECHO/DIP/BUD/2010/01017 and DCI-FOOD/2011/278-175).

Review statement

This paper was edited by Nadia Ursino and reviewed by Chi-Yuen Wang and one anonymous referee.


Agnew, D. C.: SPOTL: Some programs for ocean-tide loading, SIO technical report, Scripps Institution of Oceanography, UC San Diego, California, (last access: 22 February 2024), 2012. 

Allègre, V., Brodsky, E. E., Xue, L., Nale, S. M., Parker, B. L., and Cherry, J. A.: Using earth-tide induced water pressure changes to measure in situ permeability: A comparison with long-term pumping tests, Water Resour. Res., 52, 3113–3126, 2016. 

Barker, J. A.: A generalized radial flow model for hydraulic tests in fractured rock, Water Resour. Res., 24, 1796–1804, 1988. 

Bastias Espejo, J. M., Rau, G. C., and Blum, P.: Groundwater responses to Earth tides: Evaluation of analytical solutions using numerical simulation, J. Geophys. Res.-Sol. Ea., 127, e2022JB024771,, 2022. 

Batlle-Aguilar, J., Cook, P. G., and Harrington, G. A.: Comparison of hydraulic and chemical methods for determining hydraulic conductivity and leakage rates in argillaceous aquitards, J. Hydrol., 532, 102–121,, 2016. 

Bower, D. R.: Bedrock fracture parameters from the interpretation of well tides, J. Geophys. Res.-Sol. Ea., 88, 5025–5035, 1983. 

Bredehoeft, J. D.: Response of well-aquifer systems to earth tides, J. Geophys. Res., 72, 3075–3087, 1967. 

Burbey, T. J., Hisz, D., Murdoch, L. C., and Zhang, M.: Quantifying fractured crystalline-rock properties using well tests, earth tides and barometric effects, J. Hydrol., 414, 317–328, 2012. 

Butler Jr., J. J. and Tsou, M. S.: Pumping-induced leakage in a bounded aquifer: An example of a scale-invariant phenomenon, Water Resour. Res., 39, 1344,, 2003. 

Carr, P. A. and Van Der Kamp, G. S.: Determining aquifer characteristics by the tidal method, Water Resour. Res., 5, 1023–1031, 1969. 

Cartwright, N., Nielsen, P., and Perrochet, P.: Influence of capillarity on a simple harmonic oscillating water table: Sand column experiments and modelling, Water Resour. Res., 41, W08416,, 2005. 

Cutillo, P. A. and Bredehoeft, J. D.: Estimating aquifer properties from the water level response to Earth tides, Groundwater, 49, 600–610, 2011. 

De Marsily G.: Quantitative hydrogeology, Academic, San Diego, CA, ISBN 0-12-208915-4, 1986. 

de Pasquale, G., Valois, R., Schaffer, N., and MacDonell, S.: Contrasting geophysical signatures of a relict and an intact Andean rock glacier, The Cryosphere, 16, 1579–1596,, 2022. 

Domenico, P. A. and Schwartz, F. W.: Physical and chemical hydrogeology, 2nd edn., Wiley, Chichester, UK, ISBN 978-0-471-59762-9, 1998. 

Fuentes-Arreazola, M. A., Ramírez-Hernández, J., and Vázquez-González, R.: Hydrogeological properties estimation from groundwater level natural fluctuations analysis as a low-cost tool for the Mexicali Valley Aquifer, Water, 10, 586,, 2018. 

Gao, X., Sato, K., and Horne, R. N.: General Solution for Tidal Behavior in Confined and Semiconfined Aquifers Considering Skin and Wellbore Storage Effects, Water Resour. Res., 56, e2020WR027195,, 2020. 

Guiltinan, E. and Becker, M. W.: Measuring well hydraulic connectivity in fractured bedrock using periodic slug tests, J. Hydrol., 521, 100–107, 2015. 

Hantush, M. S.: Modification of the theory of leaky aquifers, J. Geophys. Res., 65, 3713–3725, 1960. 

Hsieh, P. A., Bredehoeft, J. D., and Farr, J. M.: Determination of aquifer transmissivity from Earth tide analysis, Water Resour. Res., 23, 1824–1832, 1987. 

Kitagawa, Y., Itaba, S., Matsumoto, N., and Koizumi, N.: Frequency characteristics of the response of water pressure in a closed well to volumetric strain in the high-frequency domain, J. Geophys. Res., 116, B08301,, 2011. 

Klönne, F.: Die periodischen Schwankungen des Wasserspiegels in den inundierten Kohlenschachten von Dux in der Periode, Sitzungsberichte der mathematisch-naturwissenschaftlichen Classe der Kaiserlichen Akademie der Wissenschaften, Wien, 8, 101–106, 1880. 

Lai, G., Ge, H., and Wang, W.: Transfer functions of the well-aquifer systems response to atmospheric loading and Earth tide from low to high-frequency band, J. Geophys. Res.-Sol. Ea., 118, 1904–1924, 2013. 

McMillan, T. C., Rau, G. C., Timms, W. A., and Andersen, M. S.: Utilizing the impact of Earth and atmospheric tides on groundwater systems: A review reveals the future potential, Rev. Geophys., 57, 281–315, 2019. 

Meinzer, O. E.: Ground water in the United States, a summary of ground-water conditions and resources, utilization of water from wells and springs, methods of scientific investigation, and literature relating to the subject, (Tech. Rep. 1938–39), U.S. Department of the Interior, Geological Survey,, 1939. 

Merritt, M. L.: Estimating hydraulic properties of the Floridan aquifer system by analysis of earth-tide, ocean-tide, and barometric effects, Collier and Hendry Counties, Florida (No. 3), US Department of the Interior, US Geological Survey,, 2004. 

Moench, A. F.: Transient flow to a large-diameter well in an aquifer with storative semiconfining layers, Water Resour. Res., 21, 1121–1131, 1985. 

Narasimhan, T. N., Kanehiro, B. Y., and Witherspoon, P. A.: Interpretation of earth tide response of three deep, confined aquifers, J. Geophys. Res.-Sol. Ea., 89, 1913–1924, 1984. 

Pacheco, F. A. L.: Hydraulic diffusivity and macrodispersivity calculations embedded in a geographic information system, Hydrol. Sci. J., 58, 930–944, 2013. 

Rabinovich, A., Barrash, W., Cardiff, M., Hochstetler, D. Ls., Bakhos, T., Dagan, G., and Kitanidis, P. K.: Frequency dependent hydraulic properties estimated from oscillatory pumping tests in an unconfined aquifer, J. Hydrol., 531, 2–16, 2015. 

Rahi, K. A. and Halihan, T.: Identifying aquifer type in fractured rock aquifers using harmonic analysis, Groundwater, 51, 76–82, 2013. 

Renner, J. and Messar, M.: Periodic pumping tests, Geophys. J. Int., 167, 479–493, 2006. 

Roeloffs, E.: Poroelastic techniques in the study of earthquake-related hydrologic phenomena, Adv. Geophys., 37, 135–195,, 1996. 

Rojstaczer, S. and Agnew, D. C.: The influence of formation material properties on the response of water levels in wells to Earth tides and atmospheric loading, J. Geophys. Res.-Sol. Ea., 94, 12403–12411, 1989. 

Schweizer, D., Ried, V., Rau, G. C., Tuck, J. E., and Stoica, P.: Comparing Methods and Defining Practical Requirements for Extracting Harmonic Tidal Components from Groundwater Level Measurements, Math. Geosci., 53, 1147–1169,, 2021. 

Sedghi, M. M. and Zhan, H.: Hydraulic response of an unconfined-fractured two-aquifer system driven by dual tidal or stream fluctuations, Adv. Water Resour., 97, 266–278, 2016. 

Shen, Q., Zheming, S., Guangcai, W., Qingyu, X., Zejun, Z., and Jiaqian, H.: Using water-level fluctuations in response to Earth-tide and barometric-pressure changes to measure the in-situ hydrogeological properties of an overburden aquifer in a coalfield, Hydrogeol. J., 28, 1465–1479,, 2020. 

Sun, X., Shi, Z., and Xiang, Y.: Frequency dependence of in-situ transmissivity estimation of well-aquifer systems from periodic loadings, Water Resour. Res., 56, e2020WR027536,, 2020. 

Thomas, A., Fortin, J., Vittecoq, B., and Violette, S.: Earthquakes and Heavy Rainfall Influence on Aquifer Properties: A New Coupled Earth and Barometric Tidal Response Model in a Confined Bi‐Layer Aquifer, Water Resour. Res., 59, e2022WR033367,, 2023. 

Valois, R., Vouillamoz, J. M., Lun, S., and Arnout, L., Assessment of water resources to support the development of irrigation in northwest Cambodia: a water budget approach, Hydrol. Sci. J., 62, 1840–1855, 2017.  

Valois, R., Vouillamoz, J. M., Lun, S., and Arnout, L.: Mapping groundwater reserves in northwestern Cambodia with the combined use of data from lithologs and time-domain-electromagnetic and magnetic-resonance soundings, Hydrogeol. J., 26, 1187–1200, 2018. 

Valois, R., Rau, G. C., Vouillamoz, J. M., and Derode, B.: Estimating hydraulic properties of the shallow subsurface using the groundwater response to Earth and atmospheric tides: a comparison with pumping tests, Water Resour. Res., 58, e2021WR031666,, 2022. 

Valois, R., Derode, B., Vouillamoz, J. M., Valerie Kotchoni, D. O., Lawson, M. A., and Rau, G. C.: Use of atmospheric tides to estimate the hydraulic conductivity of confined and semi-confined aquifers, Hydrogeol. J., 31, 2115–2128, 2023. 

Valois, R.: Remival/CambodiaData-for-leaky-and-compressible-aquiatrd: Aquifer with compressible and storing aquitard (Earth-tide), Zenodo [data set and code],, 2024. 

Van Everdingen, A. F.: The skin effect and its influence on the productive capacity of a well, J. Petrol. Geol., 5, 171–176, 1953. 

Vouillamoz, J. M., Sokheng, S., Bruyere, O., Caron, D., and Arnout, L.: Towards a better estimate of storage properties of aquifer with magnetic resonance sounding, J. Hydrol., 458, 51–58, 2012. 

Vouillamoz, J. M., Valois, R., Lun, S., Caron, D., and Arnout, L.: Can groundwater secure drinking-water supply and supplementary irrigation in new settlements of North-West Cambodia?, Hydrogeol. J., 24, 195–209, 2016. 

Wang, C. Y., Doan, M. L., Xue, L., and Barbour, A. J.: Tidal response of groundwater in a leaky aquifer – Application to Oklahoma, Water Resour. Res., 54, 8019–8033,, 2018. 

Wang, H. F.: Theory of linear poroelasticity with applications to geomechanics and hydrogeology, (Vol. 2), Princeton University Press, ISBN 9780691037462, 2000. 

Wen, Z., Zhan, H., Huang, G., and Jin, M.: Constant-head test in a leaky aquifer with a finite-thickness skin, J. Hydrol., 399, 326–334, 2011. 

Young, A.: Tidal phenomena at inland boreholes near Cradock, T. Roy. Soc. S. Afr., 3, 61–106, 1913. 

Zhang, H., Shi, Z., Wang, G., Sun, X., Yan, R., and Liu, C.: Large earthquake reshapes the groundwater flow system: Insight from the water-level response to earth tides and atmospheric pressure in a deep well, Water Resour. Res., 55, 4207–4219, 2019a. 

Zhang, S., Shi, Z., and Wang, G.: Comparison of aquifer parameters inferred from water level changes induced by slug test, earth tide and earthquake – A case study in the three Gorges area, J. Hydrol., 579, 124169,, 2019b. 

Short summary
Characterizing aquifer systems is challenging because it is difficult to obtain in situ information. They can, however, be characterized using natural forces such as Earth tides. Models that account for more complex situations are still necessary to extend the use of Earth tides to assess hydromechanical properties of aquifer systems. Such a model is developed in this study and applied to a case study in Cambodia, where a combination of tides was used in order to better constrain the model.