the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
In situ estimation of subsurface hydrogeomechanical properties using the groundwater response to semidiurnal Earth and atmospheric tides
Timothy C. McMillan
Martin S. Andersen
Wendy A. Timms
Subsurface hydrogeomechanical properties crucially underpin the management of Earth's resources, yet they are predominantly measured on core samples in the laboratory while little is known about the representativeness of in situ conditions. The impact of Earth and atmospheric tides on borehole water levels is ubiquitous and can be used to characterise the subsurface. We illustrate that disentangling the groundwater response to Earth (M_{2}) and atmospheric tidal (S_{2}) forces in conjunction with established hydraulic and linear poroelastic theories leads to a complete determination of the whole hydrogeomechanical parameter space for unconsolidated systems. Further, the characterisation of consolidated systems is possible when using literature estimates of the grain compressibility. While previous field investigations have assumed a Poisson's ratio from literature values, our new approach allows for its estimation under in situ field conditions. We apply this method to water level and barometric pressure records from four field sites with contrasting hydrogeology. Estimated hydrogeomechanical properties (e.g. specific storage; hydraulic conductivity; porosity; shear, Young's, and bulk moduli; Skempton's and Biot–Willis coefficients; and undrained or drained Poisson's ratios) are comparable to values reported in the literature, except for consistently negative drained Poisson's ratios, which is surprising. Our results reveal an anisotropic response to strain, which is expected for heterogeneous (layered) lithological profiles. Closer analysis reveals that negative Poisson's ratios can be explained by in situ conditions differing to those from typical laboratory core tests and the small strains generated by Earth and atmospheric tides. Our new approach can be used to passively, and therefore costeffectively, estimate subsurface hydrogeomechanical properties representative of in situ conditions and it improves our understanding of the relationship between geological heterogeneity and geomechanical behaviour.
A perpetual challenge for subsurface water, mineral resource, or geotechnical projects is a proper characterisation of the physical properties that may have bearings on the rate of resource extraction, operation, safety, and environmental impact of the project. The main reason for this challenge is the subsurface's heterogeneous nature and that the sampling density necessary to describe it may be prohibitively expensive (e.g. by drilling and testing of core). This issue is further exacerbated by the difficulty in approximating in situ environments in laboratory testing in regard to both scale and subsurface pressures (Hoek and Diederichs, 2006; Cundall et al., 2008; Bouzalakos et al., 2016). These difficulties may be overcome by in situ characterisation of hydrogeomechanical properties of the subsurface (Villeneuve et al., 2018). Here, the in situ pressure, stress conditions, and the scaling and inclusion of heterogeneities can achieve a more representative estimate than is possible from selective laboratory testing.
Detailed timeseries analysis of water table fluctuations in boreholes due to Earth and atmospheric tides (EAT) has been shown to be capable of deriving hydrogeomechanical properties (Hsieh et al., 1987; Rojstaczer and Agnew, 1989; Zhang et al., 2019). Further, with assumed values of some key variables, previous authors have also been able to extend the use of EAT to estimate additional subsurface hydrogeomechanical properties (Bredehoeft, 1967; Beavan et al., 1991; Cutillo and Bredehoeft, 2011). However, methods which use EAT, referred to as “tidal subsurface analysis” (TSA) techniques, remain underutilised.
EAT are natural phenomena, caused by celestial bodies (e.g. Sun or Moon), that occur throughout the Earth's crust, which have been measured and analysed in the subsurface since the mid20th century (McMillan et al., 2019). Traditionally these techniques have focused on either Earth tides (ET) (Bredehoeft, 1967; Hsieh et al., 1987; Cutillo and Bredehoeft, 2011; Zhang et al., 2019; Burbey, 2010), barometric pressure (Clark, 1967; Cutillo and Bredehoeft, 2011), or atmospheric tide (AT) loading (Acworth et al., 2016; McMillan et al., 2019; Rau et al., 2020) of the confined subsurface. Bredehoeft (1967) first proposed that once specific storage is calculated from the groundwater response to Earth tides, an aquifer porosity and compressibility can be determined from the formation pressure response to a uniformly distributed surface load such as caused by barometric pressure changes (Narasimhan et al., 1984; Rojstaczer, 1988; Rojstaczer and Riley, 1990; Ritzi et al., 1991; Burbey et al., 2012). This concept has been reiterated in the literature but, to the best of our knowledge, never solved without the use of either an assumed Poisson's ratio or bulk modulus (Cutillo and Bredehoeft, 2011) due to difficulties in attributing the superimposed EAT effects to their respective drivers (e.g. celestial body gravitational or atmospheric loading forces). Recent work by Rau et al. (2020) compares methods that estimate amplitudes and phases from monitoring datasets and concludes that harmonic least squares (HALS) is superior to fast Fourier transforms (FFT).
In this paper, the theory of the groundwater response to ET and AT is combined, thereby providing a new methodology for the estimation of the primary subsurface hydrogeomechanical properties (storage, hydraulic conductivity, and poroelastic properties). This new method improves upon the work of Cutillo and Bredehoeft (2011), as it quantitatively disentangles the groundwater response to EAT within the frequency domain, removing the influence from nonharmonic signals (e.g. precipitation and episodic recharge events) and enabling the separate and objective estimation of properties from each driver before combining the strain responses. Here, the hydraulic and linear poroelastic works of Hsieh et al. (1987), Rojstaczer and Agnew (1989), Beavan et al. (1991) and Rau et al. (2020) are integrated and combined, leading to a complete determination of the parameter space for unconsolidated systems. Further, the characterisation of consolidated systems is possible when using literature estimates of the grain compressibility (van der Kamp and Gale, 1983; Green and Wang, 1990). Finally, the new methodology is applied to groundwater and atmospheric pressure records in five boreholes from four sites to estimate hydrogeological and geomechanical properties of various consolidated and unconsolidated stratigraphies.
2.1 Extracting tidal components
Atmospheric heating and the gravitational pull of celestial bodies exert a loading of the Earth's crust (Agnew, 2010). The gravity variations and loading exerted by the movement of these celestial bodies (i.e. the Moon and Sun), as shown in Table 1, cause stress and strain responses in the Earth's crust. This causes a subsurface strain signal that is composed of numerous superimposed signals of various frequencies and amplitudes. For undrained conditions (pressurised) of either confined or semiconfined aquifers, this strain manifests as a groundwater pore pressure fluctuation (McMillan et al., 2019). A conceptual illustration of these processes is shown in Fig. 1.
Three variables are required to calculate subsurface properties using specific harmonic components (McMillan et al., 2019): (1) a computed dilatation or strain due to ETs (denoted by the superscript “ET”); (2) measured barometric pressure (denoted by the superscript AT in later equations); and (3) measured groundwater heads (denoted by the superscript “GW”). First, a moving average spanning across a period of 3 d is applied. The tidally induced frequency components are then extracted by using HALS to estimate the harmonic components caused by tides whose frequencies are well known (Hsieh et al., 1987; Xue et al., 2016; Rau et al., 2020; Schweizer et al., 2021). The moving average acts like a highpass filter and the extraction of the tidal components at specific frequencies means that lowerfrequency and episodic events (atmospheric pressure changes related to weather systems passing across a site, rainfall and recharge events, and recession from pumping or droughts) present in the groundwater level or pressure data are discarded and therefore of no consequence for the analysis. Further, since our analysis is exclusively valid for the semiconfined and confined subsurface, water movement in the overlying shallow vadose zone and unconfined saturated zone is irrelevant. The extracted frequency components are complex numbers at discrete frequencies (${\widehat{z}}_{\left(\mathrm{f}\right)}$, e.g. ${\widehat{z}}_{{M}_{\mathrm{2}}}$) for which amplitudes and phases can be calculated using the real and imaginary parts. The approach we use in our work (e.g. HALS) was comprehensively tested showing that it leads to accurate estimates of harmonics in the presence of noise levels not exceeding the measurement resolution (Schweizer et al., 2021). Highquality pressure transducers generally fulfil this criterion (Rau et al., 2019).
2.2 Earth tide influences on well water levels
2.2.1 Subsurface strain response to gravity changes
Rojstaczer and Agnew (1989) argued that for ET, horizontal areal strain is a sufficient approximation for depths of up to thousands of kilometres. This approximation is therefore sufficient for application to groundwater resources as they are generally much shallower. The strain is often referred to as “dilatation”, which is the total increase in volume of the material due to forcing by the ET (in this case the tidal pull). In porous media, assuming incompressible grains, this dilatation manifests as an opening of the total pore space, decreasing the water pressure within the material (Agnew, 2010).
In this paper, we focus on the ET component at the frequency of 1.932274 cycles per day (cpd; denoted by a subscript of its Darwin name M_{2}) and the combined ET and AT component at the frequency of 2 cpd (denoted by a subscript of its Darwin name S_{2}), described in Table 1. While other frequency components can also be used (Hsieh et al., 1988; Merritt, 2004; Cutillo and Bredehoeft, 2011), M_{2} and S_{2} generally have the strongest tidal impact and their forces remain constant over time (Acworth et al., 2016; McMillan et al., 2019).
In this paper the term “dilatation” is used broadly to represent both the dilation and compression due to the cyclical forcing of the tides, consistent with its previous use in the literature (Xue et al., 2016; Allègre et al., 2016). The distortions by dilatation can be estimated through the planar strain concept known as “tidal dilatation” (Schulze et al., 2000; FuentesArreazola et al., 2018). Tidal dilatation can be defined as
where V is the tidal potential (m^{2} s^{−2}), M_{2} is the tidal frequency, g is acceleration due to gravity (≈9.81 m s^{−2}), e^{v} is vertical displacement (–), e^{h} is horizontal displacement (–) (Agnew, 2010), R the average radius of the Earth (m) adjusted for any significant elevation, and V is the tidal potential (m^{2} s^{−2}) as defined in Table 1. The term (e^{v}−3e^{h}) may also be approximated by Love–Shida numbers where e^{v} can be replaced by ${}_{\mathrm{S}}^{\mathrm{L}}h$ with an assumed value of 0.6032 and e^{h} may be replaced with ${}_{\mathrm{S}}^{\mathrm{L}}l$ with an assumed value of 0.0839 (Agnew, 2010; Cutillo and Bredehoeft, 2011). Previous work has demonstrated the use of theoretical ET when analysing the groundwater response (Roeloffs, 1996; Xue et al., 2016; Allègre et al., 2016; McMillan et al., 2019). The terms e^{v} and e^{h} can be directly calculated from software that generates theoretical ET potential (ET^{pot}) or tidal dilatation (e) or tidal strains (ET^{ϵ}) based on geolocation, for example ETERNA (Wenzel, 1996), TSoft (Van Camp and Vauterin, 2005), or, as was done for this paper, using PyGTide (Rau, 2018). This is based on ETERNA and uses the Wahr–Dehant–Zschau model which assumes an elliptical, rotating, inelastic, and oceanless Earth (Wahr, 1981; Dehant and Zschau, 1989).
The first approach using ET to estimate specific storage used the potential for water movement from the tides to the corresponding water movement in a monitoring well in a confined aquifer for undrained conditions. Here, Bredehoeft (1967) defined specific storage (S_{s}) for a medium that is presumed to be incompressible as
where ν is the Poisson's ratio (generally assumed), ${h}_{\mathrm{S}}^{\mathrm{L}}$ and ${l}_{\mathrm{S}}^{\mathrm{L}}$ are dimensionless parameters describing Earth properties (Love–Shida numbers), and $\mathrm{\Delta}{A}_{{M}_{\mathrm{2}}}^{{\mathrm{ET}}^{\mathrm{pot}}}$ is the change in the tidal potential to the corresponding change in hydraulic head Δh. Here, the tidal dilatation (Eq. 1) has been incorporated into Eq. (2). Equation (2) was then used by Cutillo and Bredehoeft (2011) for is advantages over methods such as those provided by Hsieh et al. (1987), as it does not require the separation of individual tidal components or the knowledge of a well's dimensions. Progressive improvements in the precision and duration of gravity measurement methods have since allowed for more accurate decomposition and cataloguing of the various tidal components (Agnew, 2010). These established catalogues of precise frequencies provide the basis for component separation using harmonic filtering techniques. The full separation of ET and AT at one frequency enables their individual and combined use towards better in situ hydrogeomechanical characterisation (Rau et al., 2020).
We note that these codes do not account for ocean tide loading, i.e. the deformation of the subsurface due to the weight of the ocean tides. Ocean tide loading causes harmonic subsurface strain that is added to that imposed by ET. The actual subsurface strain amplitude variation depends on the phase of both signals and is, in the worst cases, either added to or subtracted from the ET. To understand the potential impact of this effect we used the ocean tide loading provided by Chalmers University of Technology (http://holt.oso.chalmers.se/loading/index.html, last access: 22 August 2022) to estimate the M_{2} aerial strain amplitude for our five locations with the stateoftheart finite element model FES2014b (Penna et al., 2008; Matviichuk et al., 2021). However, we further note that ocean tide loading is a complicated phenomenon (Jentzsch, 1997) and its detailed assessment is beyond the scope of this paper.
2.2.2 Well water level response to harmonically forced pore pressure
The relative amplitude response of the groundwater, as measured in a borehole in relation to the tidal dilatation or strain, can be expressed as (Hsieh et al., 1987; Xue et al., 2016; Allègre et al., 2016)
where ${\widehat{z}}_{{M}_{\mathrm{2}}}^{\mathrm{GW}}$ and ${\widehat{z}}_{{M}_{\mathrm{2}}}^{{\mathrm{ET}}^{\mathit{\u03f5}}}$ are the complex frequency component of the groundwater pressure head and tidal strain, respectively; ${A}_{{M}_{\mathrm{2}}}^{\mathrm{GW}}$ is the amplitude of the groundwater pressure head fluctuation and ${A}_{{M}_{\mathrm{2}}}^{\mathit{\u03f5}}$ is the amplitude of the tidal strain fluctuation, all at the frequency of the M_{2} tidal component. Note that ${A}_{{M}_{\mathrm{2}}}^{{\mathrm{ET}}^{\mathit{\u03f5}}}$ is also referred to as “areal strain sensitivity” (Hsieh et al., 1987).
It is important to note the difference presented in Eq. (3) from Xue et al. (2016) with the original dimensionless amplitude response calculated by Hsieh et al. (1987) as
where ${\widehat{z}}_{{M}_{\mathrm{2}}}^{\mathrm{p}}$ is the complex aquifer pore pressure response (superscript “p” reflects pore). Here, the denominator term has changed from the complex amplitude of the pressure fluctuation to the tidal dilatation, effectively incorporating Eq. (2). This key difference allows for the addition of the storage term S_{s} within the amplitude response equations due to the sensitivity of storage to the amplitude response for confined and leaky responses described in Sect. 2.2.3 and 2.2.4. ${A}_{{M}_{\mathrm{2}}}$ from Eq. (4) is dimensionless, with values $\mathrm{0}\le {A}_{{M}_{\mathrm{2}}}\le \mathrm{1}$.
The phase shift (or phase difference) is defined as the strain response observed as the complex groundwater pressure head (water level) fluctuation, minus the phase of the complex tidal dilatation (tidal forcing) stress, defined as
where ${\mathit{\varphi}}_{{M}_{\mathrm{2}}}^{\mathrm{GW}}$ is the phase angle expressed in the groundwater response and ${\mathit{\varphi}}_{{M}_{\mathrm{2}}}^{{\mathrm{ET}}^{\mathit{\u03f5}}}$ is the phase angle of the theoretical ET component, in this case at the frequency of the M_{2}. A negative phase shift is expressed where the groundwater response lags behind the induced strain (i.e. water level response lags behind the pressure head disturbance; Hsieh et al., 1987), whereas a positive phase shift indicates the groundwater response is leading the strain response.
It should be noted that in this method development, a homogeneous, isotropic aquifer of infinite lateral extent is assumed for all derivations (Hsieh et al., 1987). All derived hydrogeomechanical variables are treated as bulk properties averaged over a distinct but unknown volume, representative of the EAT area of influence around the monitoring well screened interval, including effects from geological heterogeneities and the well construction, such as the inclusion of a gravel pack. The exact nature and dimensions of the volume of influence, i.e. the volume of subsurface around the well being “sampled”, is currently unresolved. It is commonly assumed that the ET amplitude response is negligibly influenced by fluid flow when the screened aquifer is confined (Xue et al., 2016). Instead, it is predominantly controlled by a change in storage. This is used as a justification to modify the hydraulic diffusivity term in the amplitude response equations to $\mathrm{1}/Ss$ when including the ET strain estimation (Eqs. 6 and 13), i.e. the tidal dilatation (Hsieh et al., 1987; Wang, 2000; Xue et al., 2016).
2.2.3 Confined water level response
Positive and negative phase shifts are either leading (leaky) or lagging (confined), respectively, in relation to the strain response expressed by the water level in a well to formation tidal forcing. Hsieh et al. (1987) provided an analytical solution for the confined groundwater flow equation with harmonic forcing to describe the relationship between aquifer pore pressure and well water level. Their model assumes horizontal flow only and is formulated in terms of amplitude ratio and phase shift, thereby allowing for the solution of two properties, transmissivity and storativity from the amplitude and phase response. The confined (negative phase) model is defined by Hsieh et al. (1988) as
and
with the following parameters:
and
These include
and
where
Here, Ker, Kei and Ker_{1}, Kei_{1} are the real and imaginary parts of the Kelvin function of order 0 and 1, respectively. The storativity S and transmissivity T can be related to specific storage as S=S_{s}b and hydraulic conductivity as T=k^{f}b, respectively; b is the aquifer thickness, which when the aquifer thickness is unknown is approximated with the vertical screen length; r_{w} is the internal radius of the well screen (accounts for well storage); r_{c} is the outer radius of the casing. Ker and Kei are Kelvin functions of zero order, and Ker_{1} and Kei_{1} are Kelvin functions of the first order. Figure 2a and b show the amplitude and phase solution space when considering the strain response as well as separation of hydraulic properties.
2.2.4 Leaky water level response
The leaky water level model is based on the description of a periodic load on a halfspace, as described by Wang (2000), and is used for ET where a vertical pressure propagation and flow exist (Xue et al., 2016; Allègre et al., 2016). Equations (13) and (14) were derived from the force equilibrium equations (refer to Wang, 2000)
and
where z is depth of the midpoint open screen interval, ω is the angular frequency of the tidal component (M_{2}), and
Here, D_{h} is then the hydraulic diffusivity, defined as
where T is the transmissivity (T=k^{f}b), k is permeability, k^{f} is hydraulic conductivity, ρ_{w} is the density of water (0.9982 kg L^{−1} at 20 ^{∘}C), μ is the dynamic viscosity of water (8.90×10^{−4} Pa s), and S is storativity. Equations (13) and (14) require iterative solving for D_{h} and S_{s}.
Figure 3 shows the solution space for Eqs. (13) and (14). It is noteworthy that the amplitude can be attenuated significantly for high hydraulic diffusivities and shallow depths. Overall, confined conditions should prevail for low hydraulic diffusivities and larger depths where the amplitude ratio is ≈1. Semiconfined conditions are indicated by positive phase shifts. Note that hydraulic diffusivities in the range of ∼0.01–0.1 (m^{2} s) can increase the amplitudes to approx. 1.07 (Fig. 3a) for which phase shifts can assume negative values (e.g. −0.5^{∘} in Fig. 3b). Figure 2c and d show the solution space when considering the strain response as well as separation of hydraulic properties for leaky conditions.
2.2.5 Distinguishing between leaky and confined conditions
The sets of Eqs. (6) and (7) (Hsieh et al., 1987) describe lateral water movement between the well bore and surrounding subsurface, whereas Eqs. (13) and (14) (Wang, 2000) explain the positive phase shift by allowing vertical flow within the groundwater system. Both sets of equations have been used to estimate hydraulic conductivity and specific storage. This is achieved by decomposing the hydraulic diffusivity using the assumptions outlined at the end of Sect. 2.2.2.
The phase shift determines which of these sets of analytical solutions are appropriate. For a phase shift of $\mathrm{45}{}^{\circ}\le \mathrm{\Delta}{\mathrm{\Phi}}_{{M}_{\mathrm{2}}}<\mathrm{0}{}^{\circ}$, the confined response model is used, and for a phase shift of $\mathrm{0}\le \mathrm{\Delta}{\mathrm{\Phi}}_{{M}_{\mathrm{2}}}\le \mathrm{90}{}^{\circ}$, the leaky response model is applied (both are illustrated in Fig. 2). Note that the leaky model may result in a slight negative phase shift for certain parameter ranges. Consequently, there is a range of ambiguity in phase shift values between −1 and 0^{∘} in which both sets of solutions should be used, and the most physically plausible results should be selected (Xue et al., 2016). Note, the unit input as either pressure or hydraulic head will also be carried through the equations resulting in a unit difference where ${S}_{\mathrm{s}}^{\mathrm{p}}$ is specific storage expressed in 1 Pa^{−1} whereas S_{s} is specific storage expressed as 1 m^{−1}, as demonstrated in Eq. (16).
2.3 Atmospheric tide influences on well water levels
Methods that quantify the barometric efficiency of subsurface systems are based on quantifying the groundwater response magnitude to atmospheric pressure changes (Clark, 1967; Rasmussen and Crawford, 1997; Barr et al., 2000; Gonthier, 2003) or atmospheric tides (Acworth et al., 2016). Turnadge et al. (2019) reviewed these methods and concluded that the method by Acworth et al. (2016) was the most robust and reliable. However, their approach was limited by the assumption of an instantaneous and undamped response. Rau et al. (2020) developed a new method that completely disentangles the influences of EAT at the same frequency, e.g. S_{2}. This further considers the damping of the amplitude that can be caused by low hydraulic conductivity materials under confined conditions (Sect. 2.2.3) or attenuation of the amplitude under semiconfined conditions (Sect. 2.2.4). Their new approach is (Rau et al., 2020)
where
Here, ${A}_{{M}_{\mathrm{2}}}$ corrects for the damping of the subsurfacewell system and can be inferred from ET (Eq. 4) for confined conditions, e.g. for low hydraulic conductivity (Eq. 6), or damping under semiconfined conditions (Eq. 13). ${\widehat{z}}_{{S}_{\mathrm{2}}}^{\mathrm{GW}\cdot \mathrm{AT}}$ is the S_{2} component of the groundwater response to AT, and ${\widehat{z}}_{{S}_{\mathrm{2}}}^{\mathrm{AT}}$ is the S_{2} frequency component (AT) embedded in atmospheric pressure measurements. ${\widehat{z}}_{{S}_{\mathrm{2}}}^{\mathrm{GW}\cdot \mathrm{ET}}$ is the S_{2} component of the groundwater response to Earth tides (e.g., potential or strain). BE forms a stress balance, described as (Jacob, 1940)
where γ is the loading efficiency.
2.4 Combining ET and AT responses
2.4.1 General relationships
Within the following derivations it is assumed that ET only induce horizontal areal strain (${\mathit{\u03f5}}_{xy}={\mathit{\u03f5}}_{x}+{\mathit{\u03f5}}_{y}$) whereas AT only induce vertical strain (${\mathit{\u03f5}}_{z}={p}^{\mathrm{AT}}$) (Rojstaczer and Agnew, 1989; Cutillo and Bredehoeft, 2011), all of which are assumed to act instantaneously on the subsurface as is consistent with previous literature (Wang, 2000; Rau et al., 2018). Under such conditions, van der Kamp and Gale (1983) has shown that the rigidity modulus (also known as “the shear modulus”, G) can be estimated, with the previous outlined assumptions, from combined ET and AT influences as
where, ${A}_{{M}_{\mathrm{2}}}^{{\mathrm{ET}}^{\mathit{\u03f5}}}$ originates from ET (Eq. 3), whereas BE or γ is derived from AT (Eq. 19).
The disentanglement of ET and AT from the groundwater level response in a well, and the subsequent use of these separate frequency components to quantify hydrogeomechanical properties, allows further geomechanical derivations to be made. Two methods are presented below which solve for the assumption of either incompressible (suitable for unconsolidated material) or compressible grains (suitable for consolidated material). The choice between which method to use is established by examining an estimated Biot–Willis coefficient defined as (Wang, 2000)
where K is the bulk modulus, K_{s} is the bulk modulus of the solid grain, β is the compressibility, and β_{s} is the compressibility of the solid grain. For unconsolidated conditions, where the bulk modulus is much smaller than the bulk modulus of the grains (K≪K_{s}), it is possible to assume that the grains do not contribute to the overall compressibility, thus the grains are incompressible. The Biot–Willis coefficient α→1 shows that the contribution of the grains to the compressibility of the bulk material is insignificant (Rau et al., 2018). By contrast, in consolidated cases K becomes larger, leading to a coefficient that deviates appreciably from 1 (α<1). In such cases, the grain compressibility is a significant proportion of the total material compressibility and must be accounted for.
2.4.2 Unconsolidated systems
In unconsolidated systems, the compressibility of the grains is much smaller than that of the bulk leading to α=1. This enables simplification of the theory. In the following, undrained parameters are indicated by the superscript “u”. The uniaxial loading efficiency is related to the uniaxial bulk properties as (van der Kamp and Gale, 1983)
where β_{f} is the compressibility of the fluid ($\mathrm{4.59}\times {\mathrm{10}}^{\mathrm{10}}$ Pa^{−1} at 20 ^{∘}C), ${\mathit{\beta}}_{\mathrm{v}}^{\mathrm{u}}$ is the vertical undrained bulk compressibility, and θ is the total porosity of the formation. The uniaxial specific storage (assuming incompressible grains) is defined by Jacob (1940) as
Acworth et al. (2016) used Eq. (22) to simplify Eq. (23) as follows:
However, this requires a prior estimate of the porosity θ, which is often difficult to determine due to the lack of independent measurements.
Note that the above equations assume that barometric loading is uniaxial, and as such use vertical compressibility (${\mathit{\beta}}_{\mathrm{v}}^{\mathrm{u}}$) rather than the volumetric (bulk) compressibility (β^{u}). Here, we instead propose using the S_{s} derived from the response to ET (Sect. 2.2) to estimate the subsurface porosity by rearranging Eq. (24) (similar to Jacob, 1940) as
To achieve a similar outcome as Acworth et al. (2016) this porosity, in addition to the calculated S_{s}, can also be used to provide a uniaxial (vertical) bulk compressibility (inverse vertical undrained bulk modulus (${K}_{\mathrm{v}}^{\mathrm{u}}$)) of the subsurface defined as (Acworth et al., 2016)
This approach is similar to the one of Cutillo and Bredehoeft (2011) but uses the objective BE method developed by Rau et al. (2020) instead of the subjective correlation by Gonthier (2003). In this subsection it has been shown that it is possible to derive an estimate of porosity from a loading strain if the specific storage is known. This assumes incompressible grains (α≈1) and is therefore applicable for unconsolidated materials (Rau et al., 2019).
The assumption of incompressible grains allows for the removal of the grain compressibility and provides a simplification of the poroelastic parameter space. This step combined with the new derivation of the shear modulus enables a linear analytical solution of the remaining elastic variables in unconsolidated material (α≈1). The first step can be taken by deriving the undrained bulk modulus (K_{u}) with the ${K}_{\mathrm{v}}^{\mathrm{u}}$ from Acworth et al. (2016) as (Wang, 2000)
which makes it possible to solve the Skempton coefficient defined as (Rau et al., 2018)
Determination of the Skempton coefficient along with the loading efficiency unlocks the undrained Poisson's ratio using (Wang, 2000)
and drained Poisson's ratio as (Wang, 2000)
Determination of the drained Poisson ratio further unlocks all remaining poroelastic properties such as Young's modulus (E), defined as (Wang, 2000)
Equations (22)–(31) define the complete parameter space for unconsolidated materials.
2.4.3 Consolidated systems
To determine the poroelastic properties of consolidated materials, the grain compressibility must be considered leading to α<1. Further, the following two assumptions must be made: (1) although pore fluids technically respond to cubic strains, the areal strain can be used to estimate the subsurface strain from ET; (2) the system is homogeneous and laterally extensive, thus ignoring topographic effects and considering the barometric loading to be uniform. The equations that define the remaining elastic properties for such conditions are (Beavan et al., 1991) Skempton's constant
the total porosity
the Biot–Willis coefficient
and the specific storage
Equations (32)–(35) form a nonlinear system which must be solved by iteration.
If the petrology of the lithology is known, appropriate literature compressibility values of the dominant grain mineralogy (K_{s}) could be used. Quartz is the most common naturally occurring mineral and is also one of the least compressible (it is also applicable for most of our case sites), and it will therefore be used to define the upper bounds of K_{s} here. Richardson et al. (2002) summarised literature values of polycrystalline quartz for K_{s} to range between 36 and 40 GPa, and reported K_{s} values for the quartz Ottawa Sandstone to be in a range of 30–50 GPa. The average of these ranges has been summarised as 42 GPa (Rau et al., 2018) and will be used in this work.
With the established inputs of γ(BE), ${A}_{{M}_{\mathrm{2}}}$, G (Eq. 20), S_{s}, and an estimate of K_{s}, it is possible to simultaneously solve Eqs. (32)–(35) for Skempton's coefficient (B), porosity (θ), and Biot–Willis coefficient (α) (Beavan et al., 1991). This enables a complete calculation of all remaining poroelastic properties using the relationships that can be found in Wang (2000).
3.1 Field sites, datasets, and interpretation
To demonstrate the new method, groundwater and barometric pressure records from four field sites and five monitoring bores were used. These sites were selected based on three main criteria: (1) data availability and quality; (2) a significant groundwater response to ET (M_{2}); (3) providing a variety of hydrogeological settings with existing studies for parameter comparisons. Further details about each site are described in separate sections below. Specific bore geometries and measurements used in the analysis of these sites, such as depths and bore construction, casing and screen radius as well as screen depth and length, are summarised in Table 2.
Groundwater pressure head and barometric pressure time series were recorded at subhourly intervals at all sites (e.g. two sites shown in Fig. 4) for at least 3 months, which is longer than the ∼28 d being suggested as the minimum requirement of the HALS method used here (Schweizer et al., 2021). Geoposition of the boreholes, theoretical ET strain for the same duration, and sampling frequency of each site as well as aerial strain amplitudes from ocean tide loading for M_{2} are summarised in Table 2.
All time series were detrended by fitting a linear function through a moving 3 d window using the SciPy detrend function, and the main tidal components were extracted using HALS (Sect. 2.1). Tables 3 and 4 present the estimated hydrogeomechanical properties for all the field sites and are discussed in the sections below.
In this paper, all of the methodology and equations were implemented in the Python programming language, and joint iterative solving was completed using SciPy's leastsquares functionality. The following realistic parameter bounds were considered during root finding: $\mathrm{0}\le B\le \mathrm{1}$, $\mathrm{1}\le \mathit{\nu}\le \mathrm{0.5}$, $\mathrm{0.005}\le G\le \mathrm{40}$ GPa, $\mathrm{0}\le \mathit{\theta}\le \mathrm{0.5}$. We note that (1) the parameter units were scaled to avoid bias towards parameters with large values, (2) the solver was set to 64bit machine precision (epsilon $\sim \mathrm{1.11}\times {\mathrm{10}}^{\mathrm{16}}$), and (3) none of the estimated parameters reached or exceeded any of the prescribed solution constraints.
3.2 Cattle Lane (NSW, Australia)
Cattle Lane is located on the Liverpool Plains, NSW, eastern Australia. Erosion of the basaltic Liverpool Range to the south produced a succession of unconsolidated silts, clays, sands, gravel, and minor carbonate nodules within the Liverpool Plains. A thick sequence of claybound sediments overlie a gravel aquifer at 40 m. This aquifer has previously be shown to respond to loading by rainfall events (Timms and Acworth, 2005). The lithology of the 1 m screened interval was described by Acworth et al. (2015) as major basalt fragments mixed with coarse sand, shell, and carbonate nodules. The site has previously been cored to 31.5 m depth, lithologically logged and geophysically surveyed, confirming that it is horizontally extensive (Acworth et al., 2015). Crosshole seismics were also conducted by Rau et al. (2018) to the depth of 40 m (the screened interval of bore BH30061 is at 55 m depth, see Table 2), providing depth profiles of seismically inferred elastic variables that were used to constrain the pore pressure response to atmospheric tides analysis.
Further studies at this site include Acworth et al. (2016) and Acworth et al. (2017), which were precursors to Rau et al. (2018) in the investigation of pore pressure response to atmospheric tides, and Timms et al. (2018) on a core scale analysis of the site's laterally extensive and thick aquitard. Due to the sufficient M_{2} response, time series data of groundwater pressure heads measured between 21 January 2016 and 14 April 2018 were used from bore BH30061. The groundwater pressure heads were collected using vented InSitu Troll 700H series loggers at hourly intervals and submillimetre precision. Atmospheric pressure was measured by an InSitu Baro Troll absolute gauge transducer.
The borehole BH30061 from Cattle Lane produced positive M_{2} phase shifts (Table 2), with specific storage and hydraulic conductivity therefore being derived from the leaky model (Sect. 2.2.4). When applying the grain compressibility of quartz (K=42 GPa) a Biot–Willis coefficient of 0.99 is obtained and hence justifies the assumption of incompressible grains (α≈1) and the unconsolidated analytical model (Sect. 2.4.2).
The specific storage, hydraulic conductivity, porosity, shear modulus, and undrained Poisson's ratio from Cattle Lane are consistent with literature values for the sediment type (Bowles, 1996), and comply with previous estimates from higher in the stratigraphy at the same site obtained by crosshole seismics (Acworth et al., 2015, 2016; Rau et al., 2018). Young's modulus of 0.4 GPa deviates from the expected material range reported in the literature for an unconsolidated clay, sand, and gravel mixture of between 0.025 and 0.2 GPa, although the obtained value is reasonable when considering the degree of consolidation at 55 m depth and the in situ determination rather than those obtained from laboratory tests (Bouzalakos et al., 2016). The Poisson ratio of −0.29 is the only parameter that deviates significantly from the expected range of 0.2–0.5. This will be discussed later.
3.3 Thirlmere Lakes (NSW, Australia)
The Thirlmere Lakes are located in the southwest of the Sydney Basin, NSW. Both bores are located in the quartz arenite Hawkesbury Sandstone, which is about 100 m thick at the site. This sandstone is deposited by a braided river with the heterogeneous deposits showing overlapping and selfincised finingup sequences, with overbank deposited fines at paleochannel boundaries (Miall and Jones, 2003). There is evidence that the upper portion of bore Thirlmere 2 passes through a geological fault damage zone, with drilling fluid losses recorded above the screened interval due to fractures (Impax, 2019). Other studies in the same lithology include the work by Ross (2014), who investigated the potential for a bore field development within the Hawkesbury Sandstone; however, no publicly available studies exist for this lithology at the case site.
The time span and collection frequency of the time series data for the two bores differ. The time series for GW075409.1.2 (Russell, 2012) covers the period from 3 July to 14 December 2018 and was downloaded from the WaterNSW realtime data portal with 15 min intervals. For this bore, coinciding barometric time series data were obtained from a weather station approx. 500 m away. The bore Thirlmere 2 is located about 2 km from GW075409.1.2 and geocoordinates are shown in Table 2. The time series data for Thirlmere 2 were collected for this study using a vented InSitu Troll 700H series pressure transducer every 5 min between 23 July and 29 October 2019. The coinciding barometric time series was collected using a Solinst Barologger installed in the airspace of the borehole. The field data are summarised in Fig. 4 together with the theoretical ET data for the site.
The borehole GW075409.1.2 at Thirlmere Lakes produced positive M_{2} phase shifts (Table 2), with specific storage and hydraulic conductivity therefore being derived from the leaky model (Sect. 2.2.4). Both of the quartz sandstone bores returned Biot–Willis coefficients of 0.96 which requires a value for the grain compressibility (Sect. 2.4.3).
Estimates of hydrogeomechanical parameters – S_{s} of $\mathrm{3.2}\times {\mathrm{10}}^{\mathrm{6}}$ and $\mathrm{2.8}\times {\mathrm{10}}^{\mathrm{6}}$ (1 m^{−1}); k^{f} of $\mathrm{4.8}\times {\mathrm{10}}^{\mathrm{7}}$ and $\mathrm{1.9}\times {\mathrm{10}}^{\mathrm{5}}$ (m s^{−1}) – for the two sandstone bores are considered realistic for a quartz sandstone in this area. The higher k^{f} for Thirlmere 2 is believed to be indicative of enhanced hydraulic conductivity due to fractures. For this sandstone formation, SCA (2005, 2006) has previously reported S_{s} values of $\mathrm{2.49}\times {\mathrm{10}}^{\mathrm{6}}$ to $\mathrm{2.41}\times {\mathrm{10}}^{\mathrm{4}}$ (1 m^{−1}) and k^{f} of $\mathrm{1.15}\times {\mathrm{10}}^{\mathrm{6}}$ to $\mathrm{3.36}\times {\mathrm{10}}^{\mathrm{6}}$ (m s^{−1}) within this formation, including fracture networks (Ross, 2014). Our estimate of the shear modulus of 2.6 GPa marginally exceeds the expected range of 1–2 GPa (Bertuzzi, 2014; Zhang et al., 2016; Zhang and Lu, 2018). Conversely, the bulk modulus and Young's modulus both fall within the expected ranges of 2.6–5.3 and 3–8 GPa, respectively. The estimated Poisson's ratios of −0.64 and −0.03 are low compared to values between 0.2 and 0.3 typically measured in the laboratory (McMillan et al., 2019).
3.4 Dodowa (Ghana)
Dodowa is located in the Shai Osudoku District in the southeastern part of the Greater Accra Region, Ghana. The local geology consists of the Togo Structural and Dahomeyan Structural units. The Togo is composed of a series of metamorphic and folded quartzites, phyllites, and schists, while the Dahomeyan is composed of altered belts of acid and basic gneisses. BH11 used in this paper is located in a Dahomeyan gneiss (Attoh et al., 1997). All units within the region appear highly weathered, resulting in a 5 m unconsolidated regolith, confining the underlying fractured igneous and metamorphic units.
BH11 was installed and previously studied by Foppen et al. (2020), including AT analysis. The time series for the water levels of BH11 was collected at 20 min intervals between 18 October 2016 and 7 June 2017 using MiniDiver DI501 (Schlumberger pressure transducers), with atmospheric pressure being recorded with a MiniDiver DI500 (Schlumberger barometric diver), located above ground (Foppen et al., 2020).
The gneiss bore returned a Biot–Willis coefficient of 0.84, and therefore required a value for the grain compressibility (Sect. 2.4.3).
The hydrogeomechanical estimates of hydraulic conductivity of $\mathrm{2.6}\times {\mathrm{10}}^{\mathrm{6}}$ (m s^{−1}) and specific storage of $\mathrm{7.2}\times {\mathrm{10}}^{\mathrm{7}}$ (1 m^{−1}) are comparable with the values for the Togo Structural Unit from Foppen et al. (2020) derived from pumping and slug tests, which indicated ranges of 10^{−5}–10^{−6} (m s^{−1}) and $\mathrm{2.3}\times {\mathrm{10}}^{\mathrm{7}}$–$\mathrm{7.7}\times {\mathrm{10}}^{\mathrm{8}}$ (1 m^{−1}), respectively. The estimated porosity of 0.08 for BH11 slightly exceeds the range of 0.005–0.05 in Foppen et al. (2020). Comparison of elastic modulus is problematic for schists, as values are dependent on the original protolith and may vary significantly, and because schistose rock masses are known for high values of anisotropy (Hoek and Diederichs, 2006). For example, Young's modulus for a schist, as in the screened interval of BH11, can vary significantly between 21 and 117 GPa depending on mineralogy and foliation orientation (Condon et al., 2020). Our estimated value of 34.2 GPa falls within this range. However, detailed mineralogy does not exist for this bore to enable a closer comparison with literature values.
3.5 Death Valley (California, USA)
The Death Valley site is located in the western part of the United States on the border of Nevada and California. Bore BLM1 is located in Paleozoic carbonate rock and was left as an open well. The same time series record was also used in Rau et al. (2020) and it is the same bore for which data were analysed in Cutillo and Bredehoeft (2011). Data were recorded at 15 min intervals using an InSitu Troll with a vented cable and an InSitu Barotroll. The time series spans between 25 June and 16 December 2009.
As the Death Valley dataset has previously been analysed by Rau et al. (2020), who justified the use of confined ET, this method was used for the −1^{∘} phase shift of the Death Valley dataset which is located in the limit between leaky and confined models.
The estimated hydraulic conductivity of $\mathrm{4.3}\times {\mathrm{10}}^{\mathrm{6}}$ (m s^{−1}) is in agreement with the ET analysis derived value of $\mathrm{1.3}\times {\mathrm{10}}^{\mathrm{6}}$ (m s^{−1}) by Cutillo and Bredehoeft (2011). By contrast, the estimated specific storage value of $\mathrm{6.7}\times {\mathrm{10}}^{\mathrm{7}}$ (1 m^{−1}) is an order of magnitude smaller than the value of $\mathrm{7.3}\times {\mathrm{10}}^{\mathrm{6}}$ (1 m^{−1}) from Cutillo and Bredehoeft (2011). However, the specific storage and hydraulic conductivity values are both consistent with the values published by Rau et al. (2020) for the same dataset, using a method based on ET. The determined porosity (0.06) also aligns with the lower end of the range proposed by Cutillo and Bredehoeft (2011), and it is reasonable to assume the calculated Young's and shear modulus values of 28.28 and 24.10 GPa are similarly plausible as they compare with literature values (Parent et al., 2015). We note that the derived Poisson's ratio of −0.23 differs significantly from the value of 0.25 which was merely assumed in Cutillo and Bredehoeft (2011).
4.1 Influences on the quantification of properties
While quantifying hydrogeomechanical properties it is beneficial to consider factors that influence the accuracy of parameter estimation. We note that absolute errors (i.e. offsets) in the measurement of atmospheric or groundwater pressure are irrelevant for the presented methodology. This is because all parameters are based on the amplitude and phase of tidal components embedded within the measurements, i.e. the relative change characteristic of the harmonics. Here, the resolution of the measurement device directly determines detection and quantification of the responses to tidal forces. Schweizer et al. (2021) demonstrated that extraction of harmonics using HALS is accurate if the resolution of the pressure transducer is larger than twice the amplitude of the tidal component under consideration. All instruments deployed in the field examples of this work comfortably fulfilled this criterion. Schweizer et al. (2021) further noted that HALS outperforms the discrete Fourier transform, but also that devising an objective error estimation for HALS is difficult, as it depends on the nature of the residuals (difference between measurement and model), and this requires further investigation. We consider that the accuracy of HALS is at least as good as that resulting from fitting a conceptual model to measurements obtained during standard aquifer testing.
Previous works have illustrated that quantifying BE by disentangling the groundwater response to EAT based on theoretical ET does not lead to additional uncertainty in parameter estimation since it evaluates the relative responses between ET and GW (Acworth et al., 2016). However, this observation is valid only in a subsurface where the hydraulic conductivity is $<K>\mathit{\u2a86}\mathrm{1}\times {\mathrm{10}}^{\mathrm{5}}$ m s^{−1} (Rau et al., 2020). The borehole water level response in lower $<K>$ environments is damped and shifted compared to the pore pressure response outside the bore. A correction requires knowledge of both $<K>$ and S_{s}, which can be quantified using calculated ET strains. While this has been done before (Hsieh et al., 1987; Xue et al., 2016; Allègre et al., 2016; Rau et al., 2020), there is no literature investigating its associated uncertainties.
4.2 Harmonic disentanglement enables estimation of the poroelastic parameter space
In this work, we make use of recent advances that enable quantitative disentanglement of the groundwater response to both ET and AT forces. Since the two mechanisms act differently on the subsurface, the disentangled responses can be recombined to maximise the outcome. This allows us to solve the complete poroelastic space for unconsolidated systems entirely based on time series of measured groundwater pressure heads, atmospheric pressure, and theoretical ET. For consolidated systems, the complete poroelastic space can also be solved through a system of nonlinear equations by assuming a value for grain compressibility, an approach previously used in Rau et al. (2018). While the results in Tables 3 and 4 provide reasonable values, except for Poisson's ratio which is discussed separately below, they are difficult to validate because independent measurements of poroelastic properties are rare.
It is wellestablished in the literature that a negative phase shift between strain and borehole water levels is representative of confined conditions and only horizontal flow between the formation and the bore (Bower, 1983; Hsieh et al., 1987; Kümpel, 1997; Schulze et al., 2000). Conversely, the meaning of positive phase shifts is not well established in the literature. Although Sect. 2.2.4 is based on the assumption that positive phase shifts relate to a component of vertical flow between the borehole screen and the water table, i.e. semiconfined or leaky conditions (Roeloffs et al., 1989), other explanations for positive phase shifts exist in the literature. These include the influence of fracture transmissivity and length as well as the influence of ocean loading, heterogeneous material properties, and topographic effects (Roeloffs, 1996; Burbey, 2010). Here, positive phases from either vertical flow or fracture flow describe a process in which pressure is able to be propagated rapidly, either to the water table or along a highly transmissive fracture (Bower, 1983). Other mechanisms for phase shifts have also been explored in the broader literature, such as by Hanson and Owen (1982), who related fracture orientation (strike and dip) to either positive or negative phase shifts.
Our results from various field sites show both positive and negative phase shifts. A comprehensive understanding of the causes and interpretation of phase shifts is still lacking in the scientific literature. Shi and Wang (2016) observed that negative phase shifts were indicative of predominantly horizontal groundwater flow in a completely undrained system, while a positive phase shift was indicative of a vertical flow in a semiconfined or unconfined system. The method by Hsieh et al. (1987) outlined above as the confined model (negative phase shift), which was used by Shi and Wang (2016), is based on the assumption of radial horizontal flow into a well. If a positive rather than negative phase shift is used as an input into the system of equations provided by Hsieh et al. (1987), the results will not be sensible. As such, a positive phase shift model is required. For this project the method provided by Wang (2000), and adapted by Xue et al. (2016) and Allègre et al. (2016), was implemented to account for vertical flow. Both Eqs. (13) and (14) were developed for harmonic loading (i.e. ocean or barometric loading) where strain is produced at the Earth's surface and propagated down (Wang and Davis, 1996). ET do not act by surface loading; instead the mechanism is tidal dilatation, where gravitational forces act directly on all mass across the vertical profile. This points to possible issues with simplified conceptual models and the validity of their assumptions. Further research is required to independently validate results derived from passive methods that are based on simplified conceptual understanding and their analytical solutions and to test the influence of different and more complex realworld conditions, such as geological heterogeneity at different scales.
4.3 Strain responses reveal subsurface heterogeneity and anisotropy
Combining ET and AT responses in the subsurface analysis is based on the principle that EAT induce strains with a different directionality. ET are fundamentally cubic, but are approximated as planar (tidal dilatation or strain) (Schulze et al., 2000; FuentesArreazola et al., 2018). However, Rojstaczer and Agnew (1989) stated that the use of the horizontal areal strain from ET is a sufficient approximation for subsurface depths of up to thousands of kilometres, which should cover depths that are relevant to most groundwater systems. For ET, the well water levels must respond to strain in the vicinity of the well bore screen, although the sensitivity to strain that is more distant from the screen is unknown. The subsurface strain response to ETinduced stress depends on the elastic properties, which are highly heterogeneous on a small scale. However, the pore pressure response as measured by a well intersects a larger volume and should therefore be representative of the theoretical values derived from ET calculations.
Rojstaczer and Agnew (1989) predicted that the groundwater response to ET (areal strain) should be high for low porosity and low compressibility. Similarly, for such conditions, the barometric efficiency should approach 1 (BE→1, or equivalently γ→0). However, this is not reflected in our results for Death Valley and Dodowa where the groundwater response magnitude to ET is large but BE is significantly smaller than unity. This phenomenon can be explained by the fact that BE is estimated vertically across a typical geological profile as a surface load, uniaxially compressing the subsurface. Here, consolidation generally increases with depth and we hypothesise that the AT response vertically integrates the material properties above the well screen, i.e. the result is representative of the vertical heterogeneity in elastic properties encountered. The precise geometry of the representative volume from either ET or AT is currently unknown, but it is assumed to be equivalent. However, if this assumption is flawed and the representative volumes of ET or AT significantly differ, strain anisotropy may exist between these two forces and complicate their joint interpretation. Detailed field experimentation or coupled hydraulicgeomechanical modelling would be required to explore these processes.
4.4 Possible reasons for discrepancy in poroelastic properties
Our results in Tables 3 and 4 largely comply with previously established values (Wang, 2000), except for the observation of negative Poisson's ratios. It is important to note that previous studies typically assume a literature value for Poisson's ratio when calculating geomechanical properties (Cutillo and Bredehoeft, 2011). Our new approach based on tidal disentanglement removes the need for such an assumption. However, the negative Poisson's ratios are a surprising result and require explanation. We investigated whether the influence from ocean tide loading could be responsible. However, considering the maximum possible influence from ocean tide loading for each site did not lead to positive Poisson's ratios.
It is theoretically possible for Poisson's ratio to range at $\mathrm{1}\le \mathit{\nu}\le \mathrm{0.5}$ (Lakes, 1991; Lakes and Witt, 2002). Here, materials with a negative Poisson's ratio are described as auxetic, i.e. materials that become thicker parallel to the direction of the stress. The occurrence of auxetic behaviour in rocks was discussed by Gercek (2007), who summarised that as a Poisson's ratio becomes increasingly negative ($\mathit{\nu}\to \mathrm{1}$), the material becomes highly resistant to shear deformations but is easy to deform volumetrically. Ji et al. (2018) succinctly describe this relationship for conditions where the shear modulus is much greater than the bulk modulus, defined as $K<\mathrm{2}G/\mathrm{3}$, and geologically is most likely associated with highly anisotropic rocks. This ratio between the bulk and shear modulus is consistent with all results presented in this paper. As such, the negative Poisson's ratios are indicative of the subsurface laterally contracting while being vertically compressed, following the theory of linear poroelasticity.
Negative Poisson's ratios for standard uniaxial core sample testing in thermally induced microcracked granites have previously been reported (HomandEtienne and Houpert, 1989; Zhao et al., 2020). However, auxetic behaviour in rock is predominantly found in studies involving low strains and low confining pressures. For example, in the Berra Sandstone, Handin et al. (1963) observed that small compressive strains (i.e. less than 200 Bar≈2000 m H_{2}O or 20 MPa) for confining pressure conditions cause the dilatation of pore spaces. Similarly, observations of pore volumes remained constant for moderate strains (20–50 MPa) and reduced in volume for large strains (>50 MPa). Ji et al. (2018) have recently examined auxetic behaviour over a broad range of lithologies and pressures. They concluded that negative Poisson's ratios are possible in crystalline igneous and metamorphic rocks (nonfractured) for confining hydrostatic pressures less than 5 MPa, and less than 200–300 MPa for more quartzrich sedimentary rocks such as silt stones and sand stones. Further, Ji et al. (2018) observed that the porosity of sedimentary rocks plays an important role in controlling auxetic effects, similar to the nanoscale fabric in artificial auxetic materials (e.g. metallic foams).
The results in this paper are obtained in situ for fully saturated and leaky or confined conditions and caused by smallmagnitude strains. Such conditions differ considerably from those used in traditional laboratory techniques for determining elastic moduli (i.e. E, G, K, ν). Compared to the conditions experienced during a compressive laboratory test, or those described above by Ji et al. (2018), the strains caused by EAT are very small. For example, the loading variations caused by the AT component S_{2} is typically less than $\mathrm{9}\times {\mathrm{10}}^{\mathrm{5}}$ MPa (0.1 m H_{2}O), and the confining pressure caused by an artesian standing water level of 100 m H_{2}O equates to a confining pressure of only 0.98 MPa. Laboratory results are also well known for demonstrating bias in the sample strength, with the strength decreasing with the sample's increasing physical size. It has been found that this occurs due to the incorporation of more heterogeneity in the sample at larger scales, such as minor lithological changes or discontinuities due to fracturing or jointing (Cundall et al., 2008; Masoumi et al., 2016).
Alternative in situ methods, such as seismicbased methods (Rau et al., 2018), derived positive Poisson's ratios when passing through the same heterogeneous material at the same confining pressures. However, elastic moduli have previously been shown to be frequency dependent when saturated and under confining pressure (Wang, 1993; Tutuncu et al., 1998). Here, we hypothesise that the low frequency of the EAT induced stresses ($\sim \mathrm{2}\phantom{\rule{0.125em}{0ex}}\mathrm{cpd}\approx \mathrm{2.3}\times {\mathrm{10}}^{\mathrm{5}}$ Hz), compared to seismic propagated waves (1–100 Hz≈86 400–8 640 000 cpd), causes a highly relaxed response which allows sufficient time for pressure redistribution (Tutuncu et al., 1998). By contrast, the seismic frequency produces a localised unrelaxed or undrained response as the seismic waves pass through the subsurface, where this effect has been shown to change with the frequency (Pimienta et al., 2016). Both Tutuncu et al. (1998) and Pimienta et al. (2016) provide evidence of decreasing Poisson's ratios with decreasing frequency when below the typical undrained response domain ($<\mathrm{10}\phantom{\rule{0.125em}{0ex}}\mathrm{Hz}\approx \mathrm{864}\phantom{\rule{0.125em}{0ex}}\mathrm{000}$ cpd).
For small strains, as relevant for this study, Zaitsev et al. (2017) have shown that the occurrence of negative Poisson's ratios is not as exotic as previously thought. Considering the context of Cundall et al. (2008), Gercek (2007), and Ji et al. (2018), the negative Poisson's ratios derived by TSA in this paper seem plausible. We propose the following possible reasons which could lead to negative Poisson ratio hypotheses for conditions, such as the scale of the effective sample size, anisotropic strain responses from heterogeneities, low confining pressures, the low frequency and small strains caused by EATs, and boundary effects. Meeting the requirements of a negative Poisson's ratio at these small strains is defined by Lakes (1991)as nonaffine deformation (nonuniform between scales), noncentral forces, and in a state of preexisting strain (e.g. from overburden). The geomechanical derivations of this paper (Sect. 2.4) are based on linear poroelasticity. However, the auxetic responses observed by Ji et al. (2018) occur both linearly and nonlinearly within the negative Poisson's ratio space, depending on the confining pressure and the type of material (Zaitsev et al., 2017). Currently, no relationships between EAT and nonlinear poroelastic theory have been established in the literature. Future work in this area should therefore consider the integration of nonlinear geomechanics (Khan et al., 1991; Johnson and Rasolofosaon, 1996).
To the best knowledge of the authors, no explicit or robust relationships exist in the literature between elastic moduli results obtained in the field to those estimated from the laboratory testing of core (Leriche, 2017). Similarly, no in situ method currently exists that can derive elastic estimates representative of the large volume of material, such as surrounding a well bore screen, as has been proposed for ET (Zhang et al., 2019). It is likely that heterogeneity within almost any geological media will produce an anisotropic strain response to either ET or AT over such a large volume. Anisotropy may result in apparently atypical properties, e.g. negative Poisson's ratios, and should be investigated for the validity of generic assumptions common to most hydro or geomechanical investigations of a homogeneous, isotropic aquifer of infinite lateral extent.
4.5 Implications for passive quantification of subsurface hydrogeomechanical properties
There are several uncertainties associated with the findings of this paper, with implications for passive quantification of subsurface hydrogeomechanical properties. These uncertainties and limitations of the method are as follows:

Although subjective estimates have been attempted (Zhang et al., 2019), the size and scale of the volume of influence from either ET or AT are unknown. It is also possible that there is a difference between the size of influence for ET and AT. Further research is required to elucidate the zone of influence for which the derived properties are representative, such as by numeric modelling.

Currently the poroelastic response to EAT is considered to be linear. However, rocks have previously been shown to respond in a nonlinear manner for undrained, triaxially loaded laboratory settings, particularly at small strains (Johnson and Rasolofosaon, 1996; Zaitsev et al., 2017). As in situ derivations of rock mass (or sediment) poroelastic values without assuming one or more of the primary values (E, G, K, ν) is relatively novel, the implication of assuming linearity for the analysis of in situ properties remains unknown.

The mechanism behind leaky responses is believed to be due to a partial drained response in the subsurface. However, the exact causes of such responses are still unknown. In order for the validity of a positive phase shift model to be proven, a more comprehensive understanding of such mechanisms must be further developed.

Skin and well bore storage effects have been assumed to be negligible in this paper. However, these two effects will alter the phase responses to either ET or AT, as was shown in the recent work by Gao et al. (2020). It is important to note that any phase uncertainties mainly influence the hydraulic conductivity values. However, additional consideration of skin and well bore storage effects will increase the accuracy and confidence in results.

We note that there is very little literature reporting values let alone ranges of grain compressibility for mineralogy other than quartz, as has been discussed by Rau et al. (2018). Since this is the only real unknown, further work is required to elucidate the effect of grain compressibility uncertainty on the accuracy of the parameter estimation.
Passively characterising the subsurface using the groundwater response to natural signals may improve our understanding of the subsurface. For example, the possibility of auxetic behaviour of subsurface materials in undrained conditions (i.e. hydraulically coupled) will have implications for assessing compaction from groundwater estimates or the stability of slopes and cuttings. Here, the low strain elastic estimates from TSA may provide a lower bounding endmember for plausible ranges of properties. With further study, it may be possible to infer poroelastic properties at different confining pressures and frequencies or to provide a more accurate in situ determination of geomechanical rock properties (e.g. specific storage, strength, etc.) prior to excavation and construction of civil and mining projects. However, further research is required to elucidate the scope of validity (space and time) and transferability of hydrogeomechanical properties derived from different methods.
The method developed in this paper provides a comprehensive approach to estimate in situ hydrogeomechanical properties using tidal subsurface analysis (TSA), i.e. from the monitored groundwater response to Earth and atmospheric tides (EAT). Our new method first objectively disentangles the groundwater response to Earth tides (ET) and atmospheric tides (AT) for the dominant response frequencies (M_{2} and S_{2}). Secondly, the approach uses the amplitude and phase responses to ET and AT to determine the complete hydrogeomechanical parameter space: specific storage; hydraulic conductivity; porosity; shear, Young's, and bulk modulus; undrained and drained Poisson's ratio; and Skempton's and Biot–Willis coefficients. Unlike previous research, our new approach does not require an a priori estimate of the Poisson's ratio. However, the application to consolidated systems requires an estimate for the grain compressibility for which literaturebased values can be considered if data is unavailable.
Application of our new method to five groundwater and barometric pressure records from four different hydrogeological settings delivers physically realistic results that are consistent with previous estimates. However, we show that the in situ estimates of Poisson's ratio are consistently negative indicating auxetic behaviour. A closer look at the literature reveals that this is not unrealistic and can be attributed to an interplay between simultaneous in situ conditions that differ from those of established laboratory tests. These include a larger effective sample size with scaling effects, anisotropic strain responses due to heterogeneities (e.g. microcracking), significantly lower confining pressures, and the small strains at low frequencies caused by the EATs.
Our approach enables the estimation of the complete hydrogeomechanical parameter space in a passive way, i.e. from monitoring records of groundwater pressure head, measured atmospheric pressure, and calculated ET. The primary advantage is that all parameters are determined for the same in situ conditions and that the estimated values therefore should be internally consistent. The new method provides hydrogeomechanical properties of the larger rock mass. This is a clear advantage to methods that require taking samples to the laboratory where replicating field conditions such as in situ confining pressure and representative scale can be problematic. When combined with laboratory estimates on intact rock, it enables evaluation of scalespecific heterogeneity. Further, our method enables more monitoring bores to be tested for hydrogeomechanical properties at a lower cost compared to conventional aquifer pump testing. There is thus the possibility of better characterising the heterogeneity of aquifer properties. However, our method also raises the need for further research in key areas where significant uncertainties remain, e.g. the possibility for nonlinearity of the poroelastic response to surface loading and ET forces. Addressing the identified uncertainties could contribute towards improving subsurface monitoring and characterisation in both consolidated and unconsolidated systems.
The dataset is available on Figshare under https://doi.org/10.6084/M9.FIGSHARE.20353209.V1 (Rau et al., 2022).
GCR and TCM conceived the idea for this paper. GCR and TCM analysed the data and made the figures. MSA and WAT contributed with datasets and suggested improvements. GCR completed major revisions.
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 in published maps and institutional affiliations.
Thanks are due to Francis Andorful, George Lutterodt, and Jan Willem Foppen from the TGroUP project for drilling observation well BH11, installing a diver, and permitting the use of the Dodowa groundwater hydrograph data. We thank Paula Cutillo and Shannon Mazzei from the National Park Service (NPS) in California (USA) for providing the barometric and groundwater pressure dataset for BLM1. We thank HansGeorg Scherneck and Machiel Bos from Chalmers University of Technology for providing the ocean tide loading strains used in this manuscript. Some of the data used in this paper were collected with equipment provided by the Australian Federal Governmentfinanced National Collaborative Research Infrastructure Strategy(NCRIS).
This research has been supported by the H2020 Marie SkłodowskaCurie Actions (grant no. 835852). Timothy C. McMillan was partly supported by an Australian Government Research Training Program (RTP) scholarship. We acknowledge support by the KITPublication Fund of the Karlsruhe Institute of Technology.
The article processing charges for this openaccess publication were covered by the Karlsruhe Institute of Technology (KIT).
This paper was edited by Philippe Ackerer and reviewed by Severine Rosat and one anonymous referee.
Acworth, R. I., Timms, W. A., Kelly, B. F., Mcgeeney, D. E., Ralph, T. J., Larkin, Z. T., and Rau, G. C.: Late Cenozoic paleovalley fill sequence from the Southern Liverpool Plains, New South Wales – implications for groundwater resource evaluation, Aust. J. Earth Sci., 62, 657–680, https://www.tandfonline.com/doi/abs/10.1080/08120099.2015.1086815, 2015. a, b, c
Acworth, R. I., Halloran, L. J. S., Rau, G. C., Cuthbert, M. O., and Bernardi, T. L.: An objective frequency domain method for quantifying confined aquifer compressible storage using Earth and atmospheric tides, Geophys. Res. Lett., 43, 611–671, https://doi.org/10.1002/2016GL071328, 2016. a, b, c, d, e, f, g, h, i, j, k
Acworth, R. I., Rau, G. C., Halloran, L. J. S., and Timms, W. A.: Vertical groundwater storage properties and changes in confinement determined using hydraulic head response to atmospheric tides, Water Resour. Res., 53, 2983–2997, https://doi.org/10.1002/2016WR020311, 2017. a
Agnew, D. C.: Earth Tides, in: Geodesy: Treatise on Geophysics, Elsevier, p. 163, eBook ISBN 9780444535795, 2010. a, b, c, d, e, f
Allègre, V., Brodsky, E. E., Xue, L., Nale, S. M., Parker, B. L., and Cherry, J. A.: Using earthtide induced water pressure changes to measure in situ permeability: A comparison with longterm pumping tests, Water Resour. Res., 52, 3113–3126, https://doi.org/10.1002/2015WR017346, 2016. a, b, c, d, e, f
Attoh, K., Dallmeyer, R. D., and Affaton, P.: Chronology of nappe assembly in the PanAfrican Dahomeyide orogen, West Africa: evidence from ^{40}Ar^{39}Ar mineral ages, Precamb. Res., 82, 153–171, https://doi.org/10.1016/S03019268(96)000319, 1997. a
Barr, A. G., van der Kamp, G., Schmidt, R., and Black, T. A.: Monitoring the moisture balance of a boreal aspen forest using a deep groundwater piezometer, Agr. Forest Meteorol., 102, 13–24, 2000. a
Beavan, J., Evans, K., Mousa, S., and Simpson, D.: Estimating aquifer parameters from analysis of forced fluctuations in well level: An example from the Nubian Formation near Aswan, Egypt: 2. Poroelastic properties, J. Geophys. Res.Solid, 96, 12139–12160, https://doi.org/10.1029/91JB00956, 1991. a, b, c, d
Bertuzzi, R.: Sydney sandstone and shale parameters for tunnel design, Aust. Geomech. J., 49, 1–39, 2014. a
Bouzalakos, S., Crane, R. A., McGeeney, D., and Timms, W. A.: Stressdependent hydraulic properties of clayeysilt aquitards in eastern Australia, Acta Geotech., 11, 969–986, 2016. a, b
Bower, D. R.: Bedrock fracture parameters from the interpretation of well tides, J. Geophys. Res.Solid, 88, 5025–5035, https://doi.org/10.1029/JB088iB06p05025, 1983. a, b
Bowles, L. E.: Foundation analysis and design, McGrawHill, ISBN 13 9780071188449, 1996. a
Bredehoeft, J. D.: Response of wellaquifer systems to Earth tides, J. Geophys. Res., 72, 3075–3087, https://doi.org/10.1029/JZ072i012p03075, 1967. a, b, c, d
Burbey, T. J.: Fracture characterization using Earth tide analysis, J. Hydrol., 380, 237–246, https://doi.org/10.1016/j.jhydrol.2009.10.037, 2010. a, b
Burbey, T. J., Hisz, D., Murdoch, L. C., and Zhang, M.: Quantifying fractured crystallinerock properties using well tests, earth tides and barometric effects, J. Hydrology, 414, 317–328, https://doi.org/10.1016/j.jhydrol.2011.11.013, 2012. a
Clark, W. E.: Computing the barometric efficiency of a well, J. Hydraul. Div., 93, 93–98, 1967. a, b
Condon, K. J., Sone, H., Wang, H. F., AjoFranklin, J., Baumgartner, T., Beckers, K., Blankenship, D., Bonneville, A., Boyd, L., Brown, S., Burghardt, J. A., Chai, C., Chen, Y., Chi, B., Condon, K., Cook, P. J., Crandall, D., Dobson, P. F., Doe, T., Doughty, C. A., Elsworth, D., Feldman, J., Feng, Z., Foris, A., Frash, L. P., Frone, Z., Fu, P., Gao, K., Ghassemi, A., Guglielmi, Y., Haimson, B., Hawkins, A., Heise, J., Hopp, C., Horn, M., Horne, R. N., Horner, J., Hu, M., Huang, H., Huang, L., Im, K. J., Ingraham, M., Jafarov, E., Jayne, R. S., Johnson, S. E., Johnson, T. C., Johnston, B., Kim, K., King, D. K., Kneafsey, T., Knox, H., Knox, J., Kumar, D., Lee, M., Li, K., Li, Z., Maceira, M., Mackey, P., Makedonska, N., Mattson, E., McClure, M. W., McLennan, J., Medler, C., Mellors, R. J., Metcalfe, E., Moore, J., Morency, C. E., Morris, J. P., Myers, T., Nakagawa, S., Neupane, G., Newman, G., Nieto, A., Oldenburg, C. M., Paronish, T., Pawar, R., Petrov, P., Pietzyk, B., Podgorney, R., Polsky, Y., Pope, J., Porse, S., Primo, J. C., Reimers, C., Roberts, B. Q., Robertson, M., Roggenthen, W., Rutqvist, J., Rynders, D., Schoenball, M., Schwering, P., Sesetty, V., Sherman, C. S., Singh, A., Smith, M. M., Sone, H., Sonnenthal, E. L., Soom, F. A., Sprinkle, P., Strickland, C. E., Su, J., Templeton, D., Thomle, J. N., Tribaldos, V. R., Ulrich, C., Uzunlar, N., Vachaparampil, A., Valladao, C. A., Vandermeer, W., Vandine, G., Vardiman, D., Vermeul, V. R., Wagoner, J. L., Wang, H. F., Weers, J., Welch, N., White, J., White, M. D., Winterfeld, P., Wood, T., Workman, S., Wu, H., Wu, Y. S., Yildirim, E. C., Zhang, Y., Zhang, Y. Q., Zhou, Q., Zoback, M. D., and CollabTeam, E. G. S.: Low Static Shear Modulus Along Foliation and Its Influence on the Elastic and Strength Anisotropy of Poorman Schist Rocks, Homestake Mine, South Dakota, Rock Mech. Rock Eng., 53, 5257–5281, https://doi.org/10.1007/s00603020021824, 2020. a
Cundall, P. A., Pierce, M. E., and Mas Ivars, D.: Quantifying the Size Effect of Rock Mass Strength, in: SHIRMS 2008: Proceedings of the First Southern Hemisphere International Rock Mechanics Symposium, Australian Centre for Geomechanics, Perth, 3–15, https://doi.org/0.36487/ACG_repo/808_31, 2008. a, b, c
Cutillo, P. A. and Bredehoeft, J. D.: Estimating Aquifer Properties from the Water Level Response to Earth Tides, Ground Water, 49, 600–610, https://doi.org/10.1111/j.17456584.2010.00778.x, 2011. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p
Dehant, V. and Zschau, J.: The Effect of Mantle Inelasticity On Tidal Gravity: A Comparison Between the Spherical and the Elliptical Earth Model, Geophys. J. Int., 97, 549–555, https://doi.org/10.1111/j.1365246X.1989.tb00522.x, 1989. a
Foppen, J. W., Lutterodt, G., Rau, G. C., and Minkah, O.: Groundwater flow system analysis in the regolith of Dodowa on the Accra Plains, Ghana, J. Hydrol.: Reg. Stud., 28, 100663, https://doi.org/10.1016/j.ejrh.2020.100663, 2020. a, b, c, d
FuentesArreazola, M., RamírezHernández, J., and VázquezGonzález, R.: Hydrogeological Properties Estimation from Groundwater Level Natural Fluctuations Analysis as a LowCost Tool for the Mexicali Valley Aquifer, Water, 10, 586, https://doi.org/10.3390/w10050586, 2018. a, b
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, https://doi.org/10.1029/2020WR027195, 2020. a
Gercek, H.: Poisson's ratio values for rocks, Int. J. Rock Mech. Min. Sci., 44, 1–13, https://doi.org/10.1016/j.ijrmms.2006.04.011, 2007. a, b
Gonthier, G.: A Graphical Method for Estimation of Barometric Efficiency from Continuous Data – Concepts and Application to a Site in the Piedmont, Air Force Plant 6, Marietta, Georgia, Tech. rep., US Geological Survey, https://doi.org/10.3133/sir20075111, 2003. a, b
Green, D. H. and Wang, H. F.: Specific storage as a poroelastic coefficient, Water Resour. Res., 26, 1631–1637, https://doi.org/10.1029/WR026i007p01631, 1990. a
Handin, J., Hager Jr., R. V., Friedman, M., and Feather, J. N.: Experimental deformation of sedimentary rocks under confining pressure: pore pressure tests, AAPG Bull., 47, 717–755, 1963. a
Hanson, J. M. and Owen, L. B.: Fracture orientation analysis by the solid earth tidal strain method, in: vol. 1982, Proceedings – SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, New Orleans, Louisiana, p. 18, https://doi.org/10.2118/11070MS, 1982. a
Hoek, E. and Diederichs, M. S.: Empirical estimation of rock mass modulus, Int. J. Rock Mech. Min. Sci., 43, 203–215, https://doi.org/10.1016/j.ijrmms.2005.06.005, 2006. a, b
HomandEtienne, F. and Houpert, R.: Thermally induced microcracking in granites: characterization and analysis, Int. J. Rock Mech. Min. Sci. Geomech. Abstr., 26, 125–134, https://doi.org/10.1016/01489062(89)900016, 1989. a
Hsieh, P. A., Bredehoeft, J. D., and Farr, J. M.: Determination of aquifer transmissivity from Earth tide analysis, Water Resour. Res., 23, 1824–1832, https://doi.org/10.1029/WR023i010p01824, 1987. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q
Hsieh, P. A., Bredehoeft, J. D., and Rojstaczer, S. A.: Response of well aquifer systems to Earth tides: Problem revisited, Water Resour. Res., 24, 468–472, https://doi.org/10.1029/WR024i003p00468, 1988. a, b
Impax: Permeability testing in borehole Thirlmere 2, Tech. rep., Sigra Job Reference Number 489, Sigra Pty Ltd., 2019. a
Jacob, C. E.: On the flow of water in an elastic artesian aquifer, Eos Trans. Am. Geophys. Union, 21, 574–586, 1940. a, b, c
Jentzsch, G.: Earth tides and ocean tidal loading, Springer, Berlin, Heidelberg, 145–171, https://doi.org/10.1007/BFb0011461, 1997. a
Ji, S., Li, L., Motra, H. B., Wuttke, F., Sun, S., Michibayashi, K., and Salisbury, M. H.: Poisson's Ratio and Auxetic Properties of Natural Rocks, J. Geophys. Res.Solid, 123, 1161–1185, https://doi.org/10.1002/2017JB014606, 2018. a, b, c, d, e, f
Johnson, P. A. and Rasolofosaon, P. N. J.: Nonlinear elasticity and stressinduced anisotropy in rock, J. Geophys. Res.Solid, 101, 3113–3124, https://doi.org/10.1029/95JB02880, 1996. a, b
Khan, A. S., Xiang, Y., and Huang, S.: Behavior of Berea sandstone under confining pressure part I: Yield and failure surfaces, and nonlinear elastic response, Int. J. Plasticity, 7, 607–624, https://doi.org/10.1016/07496419(91)900462, 1991. a
Kümpel, H.J.: Tides in water saturated rock BT – Tidal Phenomena, Springer, Berlin, Heidelberg, 277–291, https://doi.org/10.1007/BFb0011467, 1997. a
Lakes, R.: Deformation mechanisms in negative Poisson's ratio materials: structural aspects, J. Mater. Sci., 26, 2287–2292, https://doi.org/10.1007/BF01130170, 1991. a, b
Lakes, R. S. and Witt, R.: Making and characterizing negative Poisson's ratio materials, Int. J. Mech. Eng. Educ., 30, –58, 2002. a
Leriche, A.: Stress estimation from borehole scans for prediction of excavation overbreak in brittle rock, PhD thesis, Queen's University, Kingston, Ontario, Canada, http://hdl.handle.net/1974/22825 (last access: 22 August 2022), 2017. a
Masoumi, H., Douglas, K. J., and Russell, A. R.: A Bounding Surface Plasticity Model for Intact Rock Exhibiting SizeDependent Behaviour, Rock Mech. Rock Eng., 49, 47–62, https://doi.org/10.1007/s0060301507448, 2016. a
Matviichuk, B., King, M. A., Watson, C. S., and Bos, M. S.: Limitations in OneDimensional (an)Elastic Earth Models for Explaining GPSObserved M_{2} Ocean Tide Loading Displacements in New Zealand, J. Geophys. Res.Solid, 126, e2021JB021992, https://doi.org/10.1029/2021jb021992, 2021. a
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, https://doi.org/10.1029/2018RG000630, 2019. a, b, c, d, e, f, g, h, i
Merritt, M. L.: Estimating hydraulic properties of the Floridan Aquifer System by analysis of earthtide, oceantide, and barometric effects, Collier and Hendry Counties, Florida, Tech. rep., USGS, https://doi.org/10.3133/wri034267, 2004. a
Miall, A. D. and Jones, B. G.: Fluvial architecture of the Hawkesbury sandstone (Triassic), near Sydney, Australia, J. Sediment. Res., 4, 531–545, 2003. a
Narasimhan, T. N., Kanehiro, B. Y., and Witherspoon, P. A.: Interpretation of Earth tide response of three deep, confined aquifers, J. Geophys. Res.Solid, 89, 1913–1924, https://doi.org/10.1029/JB089iB03p01913, 1984. a
Parent, T., Domede, N., Sellier, A., and Mouatt, L.: Mechanical characterization of limestone from sound velocity measurement, Int. J. Rock Mech. Min. Sci., 79, 149–156, https://doi.org/10.1016/j.ijrmms.2015.08.009, 2015. a
Penna, N. T., Bos, M. S., Baker, T. F., and Scherneck, H. G.: Assessing the accuracy of predicted ocean tide loading displacement values, J. Geod., 82, 893–907, https://doi.org/10.1007/s0019000802202, 2008. a
Pimienta, L., Fortin, J., and Guéguen, Y.: Effect of fluids and frequencies on Poisson’s ratio of sandstone samples, Geophysics, 81, D183–D195, https://doi.org/10.1190/geo20150310.1, 2016. a, b
Rasmussen, T. C. and Crawford, L. A.: Identifying and Removing Barometric Pressure Effects in Confined and Unconfined Aquifers, Ground Water, 35, 502–511, https://doi.org/10.1111/j.17456584.1997.tb00111.x, 1997. a
Rau, G. C.: PyGTide: A Python module and wrapper for ETERNA PREDICT to compute synthetic model tides on Earth, Zenodo [code], https://doi.org/10.5281/zenodo.1346260, 2018. a
Rau, G. C., Acworth, R. I., Halloran, L. J. S., Timms, W. A., and Cuthbert, M. O.: Quantifying Compressible Groundwater Storage by Combining CrossHole Seismic Surveys and Head Response to Atmospheric Tides, J. Geophys. Res.Earth, 123, 1910–1930, https://doi.org/10.1029/2018JF004660, 2018. a, b, c, d, e, f, g, h, i, j
Rau, G. C., Post, V. E. A., Shanafield, M., Krekeler, T., Banks, E. W., and Blum, P.: Error in hydraulic head and gradient timeseries measurements: a quantitative appraisal, Hydrol. Earth Syst. Sci., 23, 3603–3629, https://doi.org/10.5194/hess2336032019, 2019. a, b
Rau, G. C., Cuthbert, M. O., Acworth, R. I., and Blum, P.: Technical note: Disentangling the groundwater response to Earth and atmospheric tides to improve subsurface characterisation, Hydrol. Earth Syst. Sci., 24, 6033–6046, https://doi.org/10.5194/hess2460332020, 2020. a, b, c, d, e, f, g, h, i, j, k, l, m
Rau, G., McMillan, T., Andersen, M. S., and Timms, W. A.: Dataset: Insitu estimation of subsurface hydrogeomechanical properties using the groundwater response to semidiurnal Earth and atmospheric tides, figshare [data set], https://doi.org/10.6084/M9.FIGSHARE.20353209.V1, 2022. a
Richardson, N. D., Williams, K. L., Briggs, K. B., and Thorsos, E. I.: Dynamic measurement of sediment grain compressibility at atmospheric pressure: acoustic applications, IEEE J. Ocean. Eng., 27, 593–601, https://doi.org/10.1109/JOE.2002.1040941, 2002. a
Ritzi, R. W., Sorooshian, S., and Hsieh, P. A.: The estimation of fluid flow properties from the response of water levels in wells to the combined atmospheric and Earth tide forces, Water Resour. Res., 27, 883–893, https://doi.org/10.1029/91WR00070, 1991. a
Roeloffs, E.: Poroelastic Techniques in the Study of EarthquakeRelated Hydrologic Phenomena, Adv. Geophys., 38, 135–195, https://doi.org/10.1016/S00652687(08)602708, 1996. a, b
Roeloffs, E. A., Burford, S. S., Riley, F. S., and Records, A. W.: Hydrologic effects on water level changes associated with episodic fault creep near Parkfield, California, J. Geophys. Res., 94, 12387, https://doi.org/10.1029/jb094ib09p12387, 1989. a
Rojstaczer, S.: Determination of fluid flow properties from the response of water levels in wells to atmospheric loading, Water Resour. Res., 24, 1927–1938, https://doi.org/10.1029/WR024i011p01927, 1988. a
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., 94, 12403, https://doi.org/10.1029/JB094iB09p12403, 1989. a, b, c, d, e, f
Rojstaczer, S. and Riley, F. S.: Response of the water level in a well to Earth tides and atmospheric loading under unconfined conditions, Water Resour. Res., 26, 1803–1817, https://doi.org/10.1029/WR026i008p01803, 1990. a
Ross, J. B.: Groundwater resource potential of the Triassic Sandstones of the Southern Sydney Basin: an improved understanding, Aust. J. Earth Sci., 61, 463–474, https://doi.org/10.1080/08120099.2014.910548, 2014. a, b
Russell, G.: Thirlmere Lakes drilling report, Tech. rep., NSW Office of Water, Sydney, https://www.industry.nsw.gov.au/__data/assets/pdf_file/0004/153436/thirlmere_lakes_drilling_report.pdf (last access: 22 August 2022), 2012. a
SCA: Metropolitan Water Plan. Priority Groundwater Investigations for Emergency Drought Relief, Area 2: Upper Nepean Catchment, Tech. rep. 05GL31A/2, 2005. a
SCA: Upper Nepean Groundwater Pilot Studies – Pumping Test Interpretation and Data Logger Installation, Tech. rep., Project 5078 RP01, 2006. a
Schulze, K. C., Kümpel, H.J., and Huenges, E.: InSitu Petrohydraulic Parameters from Tidal and Barometric Analysis of Fluid Level Variations in Deep Wells: Some Results From KTB BT – Hydrogeology of Crystalline Rocks, Springer Netherlands, Dordrecht, 79–104, https://doi.org/10.1007/9789401718165_4, 2000. a, b, c
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, https://doi.org/10.1007/s11004020099159, 2021. a, b, c, d, e
Shi, Z. and Wang, G.: Aquifers switched from confined to semiconfined by earthquakes, Geophys. Res. Lett., 43, 111–166, https://doi.org/10.1002/2016GL070937, 2016. a, b
Timms, W. A. and Acworth, R. I.: Propagation of pressure change through thick clay sequences: an example from Liverpool Plains, NSW, Australia, Hydrogeol. J., 13, 858–870, https://doi.org/10.1007/s1004000504367, 2005. a
Timms, W. A., Acworth, R. I., Crane, R. A., Arns, C. H., Arns, J., McGeeney, D. E., Rau, G. C., and Cuthbert, M. O.: The influence of syndepositional macropores on the hydraulic integrity of thick alluvial clay aquitards, Water Resour. Res., 54, 3122–3138, 2018. a
Turnadge, C., Crosbie, R. S., Barron, O., and Rau, G. C.: Comparing Methods of Barometric Efficiency Characterization for Specific Storage Estimation, Groundwater, 57, 844–859, https://doi.org/10.1111/gwat.12923, 2019. a
Tutuncu, A. N., Podio, A. L., Gregory, A. R., and Sharma, M. M.: Nonlinear viscoelastic behavior of sedimentary rocks, Part I: Effect of frequency and strain amplitude, Geophysics, 63, 184–194, https://doi.org/10.1190/1.1444311, 1998. a, b, c
Van Camp, M. and Vauterin, P.: Tsoft: graphical and interactive software for the analysis of time series and Earth tides, Comput. Geosci., 31, 631–640, https://doi.org/10.1016/j.cageo.2004.11.015, 2005. a
van der Kamp, G. and Gale, J. E.: Theory of earth tide and barometric effects in porous formations with compressible grains, Water Resour. Res., 19, 538–544, https://doi.org/10.1029/WR019i002p00538, 1983. a, b, c
Villeneuve, M. C., Heap, M. J., Kushnir, A. R. L., Qin, T., Baud, P., Zhou, G., and Xu, T.: Estimating in situ rock mass strength and elastic modulus of granite from the SoultzsousForêts geothermal reservoir (France), Geoth. Energy, 6, 11, https://doi.org/10.1186/s4051701800961, 2018. a
Wahr, J. M.: Body tides on an elliptical, rotating, elastic and oceanless earth, Geophys. J. Roy. Astron. Soc., 64, 677–703, https://doi.org/10.1111/j.1365246X.1981.tb02690.x, 1981. a
Wang, H. F.: Quasistatic poroelastic parameters in rock and their geophysical applications, Pure Appl. Geophys., 141, 269–286, https://doi.org/10.1007/BF00998332, 1993. a
Wang, H. F.: Theory of Linear Poroelasticity with Applications to Geomechanics and Hydrogeology, Princeton University Press, ISBN 9780691037462, 2000. a, b, c, d, e, f, g, h, i, j, k, l, m, n
Wang, K. and Davis, E. E.: Theory for the propagation of tidally induced pore pressure variations in layered subseafloor formations, J. Geophys. Res.Solid, 101, 11483–11495, 1996. a
Wenzel, H.G.: Accuracy assessment for tidal potential catalogues, Bulletin d'Informations des Marées Terrestres, 124, 9394–9416, 1996. a
Xue, L., Brodsky, E. E., Erskine, J., Fulton, P. M., and Carter, R.: A permeability and compliance contrast measured hydrogeologically on the San Andreas Fault, Geochem. Geophy. Geosy., 17, 858–871, https://doi.org/10.1002/2015GC006167, 2016. a, b, c, d, e, f, g, h, i, j, k
Zaitsev, V. Y., Radostin, A. V., Pasternak, E., and Dyskin, A.: Extracting realcrack properties from nonlinear elastic behaviour of rocks: abundance of cracks with dominating normal compliance and rocks with negative Poisson ratios, Nonlin. Processes Geophys., 24, 543–551, https://doi.org/10.5194/npg245432017, 2017. a, b, c
Zhang, C. and Lu, N.: What Is the Range of Soil Water Density? Critical Reviews With a Unified Model, Rev. Geophys., 56, 532–562, https://doi.org/10.1029/2018RG000597, 2018. a
Zhang, C., Mitra, R., Oh, J., and Hebblewhite, B.: Analysis of Mininginduced Valley Closure Movements, Rock Mech. Rock Eng., 49, 1923–1941, https://doi.org/10.1007/s0060301508801, 2016. a
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, https://doi.org/10.1016/j.jhydrol.2019.124169, 2019. a, b, c, d
Zhao, Z., Xu, H., Wang, J., Zhao, X., Cai, M., and Yang, Q.: Auxetic behavior of Beishan granite after thermal treatment: A microcracking perspective, Eng. Fract. Mech., 231, 107017, https://doi.org/10.1016/j.engfracmech.2020.107017, 2020. a