# Hydro-Mechanical Characterization of a Fractured Aquifer Using Groundwater Level Tidal Analysis: Effect of Pore Pressure and Seismic Dynamic Shear Stresses on Permeability Variations

Augustin Thomas, Jerome Fortin, Benoit Vittecoq, Hideo Aochi, Sophie Violette, Julie Maury, Frédéric Lacquement, Adnand Bitri

► To cite this version:

Augustin Thomas, Jerome Fortin, Benoit Vittecoq, Hideo Aochi, Sophie Violette, et al. Hydro-Mechanical Characterization of a Fractured Aquifer Using Groundwater Level Tidal Analysis: Effect of Pore Pressure and Seismic Dynamic Shear Stresses on Permeability Variations. Journal of Geophysical Research: Solid Earth, 2024, 129 (8), 10.1029/2024jb028847. http://dx.doi.org/10.1029/2024jb028847

HAL Id: hal-04669643

https://brgm.hal.science/hal-04669643v2

Submitted on 2 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# JGR Solid Earth

## RESEARCH ARTICLE

10.1029/2024JB028847

**Key Points:**

* The 15-years' time series of permeability and shear modulus of the aquifer is deduced from the analysis of solid earth and barometric tides

* Co-seismic irreversible changes of fractured aquifer permeability are induced by earthquakes dynamic shear stresses

* Interseismic permeability variations of the aquifer are controlled by pore pressure, that is, the hydraulic head

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

A. Thomas,
augustin.thomas.at@gmail.com

**Citation:**

Thomas, A., Fortin, J., Vittecoq, B., Aochi, H., Violette, S., Maury, J., et al. (2024). Hydro-mechanical characterization of a fractured aquifer using groundwater level tidal analysis: Effect of pore pressure and seismic dynamic shear stresses on permeability variations. *Journal of Geophysical Research: Solid Earth*, 129, e2024JB028847. https://doi.org/10.1029/2024JB028847

Received 29 JAN 2024

Accepted 20 JUL 2024

**Author Contributions:**

**Conceptualization:** A. Thomas, J. Fortin, B. Vittecoq, H. Aochi, S. Violette, J. Maury

**Data curation:** A. Thomas, H. Aochi, F. Lacquement, A. Bitri

**Formal analysis:** A. Thomas, H. Aochi, J. Maury, A. Bitri

**Funding acquisition:** J. Fortin, B. Vittecoq, S. Violette

**Investigation:** S. Violette, J. Maury, F. Lacquement

**Methodology:** A. Thomas, J. Fortin, B. Vittecoq, H. Aochi, S. Violette, J. Maury

© 2024. The Author(s). This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

## Hydro-Mechanical Characterization of a Fractured Aquifer Using Groundwater Level Tidal Analysis: Effect of Pore Pressure and Seismic Dynamic Shear Stresses on Permeability Variations

A. Thomas¹,² ID, J. Fortin¹, B. Vittecoq³ ID, H. Aochi¹,⁴ ID, S. Violette¹,⁵ ID, J. Maury⁴ ID, F. Lacquement⁴ ID, and A. Bitri⁴

¹Laboratoire de Géologie, Ecole Normale Supérieure - PSL Research University, CNRS, UMR, Paris, France, ²Now at BRGM, Orléans, France, ³BRGM, Fort-de-France, Martinique, ⁴BRGM, Orléans, France, ⁵Sorbonne University, UFR 918, Paris, France

**Abstract** Groundwater level tidal analysis is a powerful technique to monitor aquifer's permeability and hence its change over time. Earthquakes are known to affect aquifer's properties, in their vicinity through static stress changes but also further away through dynamic stresses. Most often changes are in the form of permeability increases, but sometimes decreases; the changes can be either permanent or transient. These observations are relatively well documented but the physical processes behind these changes are not well understood. By combining solid-earth and barometric tidal groundwater level responses in a borehole in a coherent poroelastic theoretical framework, and a bi-layer hydrogeological model, we recover a 15 years-long time series evolution of aquifer transmissivity and shear modulus. This study showcases the full potential of the tidal analysis method, coupling pore pressure diffusion and rock deformation, at the frontier of hydrogeology and rock physics. This unprecedented measurement of permeability and shear modulus evolution by tidal analysis reveals, during interseismic period, high sensitivity of this shallow aquifer to effective stress, and thus to pore pressure. Thanks to additional finite element simulation of seismic wave propagation, we explore the different mechanisms affecting permeability and shear modulus in the studied fractured andesite aquifer. This study confirms the predominant role of seismic dynamic stresses, and more precisely of dynamic shear stresses, in the change of permeability following an earthquake.

**Plain Language Summary** Tidal oscillations of groundwater level, observable in boreholes all around the world, are information-rich signals for hydrogeologist to infer the properties of surrounding aquifers, and how they evolve over time. Most importantly, it allows to monitor aquifer permeability and how it evolves under the influence of earthquakes, especially through static stresses changes or seismic wave propagation. Here we push the limits of what tidal analysis can reveal by also retrieving the shear modulus of the studied aquifer. This new observable, combined with the computation of regional earthquakes stresses, allows us to understand better how aquifer properties are modified by earthquakes. The analysis reveals first that dynamic shear stresses are the most probable cause of permeability changes, and second that the sensitivity of our fractured aquifer permeability to pore pressure is high.

## 1. Introduction

The tidal response of groundwater level in boreholes has long been identified as a useful tool to characterize hydrogeological properties of aquifers, starting with the pioneer work of Hsieh et al. (1987) that models the response to tides in a confined aquifer, due to borehole storage effect. It enables to recover transmissivity from the phase lag computation between earth tidal strain and the piezometric level in wells. This model arrived at a time when the poroelasticity theory was still building up, and theoretical considerations about the formalism were still discussed (Bredehoeft, 1967; Hsieh et al., 1988; Narasimhan et al., 1984; Van Der Kamp & Gale, 1983). Questions about the drained or undrained nature of the tidal response, or the definition of specific storage used were already raised. The technique of tidal analysis has spread and new models were conceived to represent a broader range of observations (Brodsky & Prejean, 2005; Odling et al., 2015; Roeloffs, 1996; Rojstaczer & Riley, 1990; Thomas et al., 2023; C.-Y. Wang et al., 2018). Yet the potential of this technique is still widely
**Project administration:** J. Fortin,

B. Vittecoq, S. Violette

**Resources:** H. Aochi, F. Lacquement,
A. Bitri

**Software:** A. Thomas, B. Vittecoq,

H. Aochi, J. Maury, A. Bitri

**Supervision:** J. Fortin, B. Vittecoq,
S. Violette

**Validation:** A. Thomas, J. Fortin,

B. Vittecoq, H. Aochi, S. Violette,

J. Maury, F. Lacquement, A. Bitri

**Visualization:** A. Thomas, B. Vittecoq,
F. Lacquement, A. Bitri

**Writing – original draft:** A. Thomas,
H. Aochi

**Writing – review & editing:** A. Thomas,
J. Fortin, B. Vittecoq, H. Aochi,
S. Violette, J. Maury, F. Lacquement,
A. Bitri

underutilized given the amount of available borehole data (McMillan et al., 2019), and this despite the relative availability of useful software and tools (e.g., Doan et al., 2008).

More specifically, this technique was used repeatedly to infer aquifer permeability changes following earthquakes (Elkhoury et al., 2006; Vittecoq et al., 2020; Xue et al., 2013) thanks to the continuity of the measurement over time. It resulted in the production of permeability time series which finally enabled accurate quantification of the permeability changes. These permeability changes were known to occur thanks to indirect observations like well level changes, or spring discharge changes. Manga et al. (2012) gathered such observations, and backed by laboratory studies (e.g., Elkhoury et al., 2011), proposed several mechanisms for permeability changes, including particle mobilization or more generally flow pathways clogging/unclogging (also found in Barbosa, Hunziker, et al., 2019; Boeut et al., 2020; Brodsky et al., 2003; Candela, 2014; Elkhoury et al., 2011; Thomas et al., 2023; Vittecoq et al., 2020; H. Wang et al., 2020), gas bubble mobilization and finally elastic processes like change in fracture density or aperture with fluid pressure (Faoro et al., 2012) or rock fracturation and healing (Xue et al., 2013). In particular, particle mobilization is today acknowledged to be the main process at play, particularly in the case of fractured aquifers, where fracture clogging/unclogging is almost universally cited as the responsible mechanism.

In a previous study performed on another aquifer from La Martinique (Thomas et al., 2023), we were able to invert the diffusivity and the loading efficiency of an aquitard thanks to tidal analysis of earth and barometric signals. Like similar studies with the same methodology (Valois et al., 2022; Xue et al., 2013), we quantified the impacts of earthquakes on permeability. Yet the physical mechanisms of the observed change in permeability remained at the stage of hypotheses. In this study, we wish to find quantitative evidence in the field to decipher between physical processes, especially between particle mobilization and poroelastic processes due to change in the fluid pressure. To that end, we aim to not only infer aquifer permeability but also shear modulus from the tidal signal inversion. Indeed shear modulus, which can be related to fracture density (David & Zimmerman, 2011; Fortin et al., 2007, 2011; Kachanov, 1993; O'Connell & Budiansky, 1974), will provide a new window to the micromechanics at play in the aquifer.

To reach this objective, we will adapt the idea of Van Der Kamp and Gale (1983) to formulate the earth and barometric source terms such that shear modulus is invertible from the model, which has never been done to the best of our knowledge. We will then apply this formalism to a bilayer aquifer geometry, unifying in a more general way the models of Odling et al. (2015) for barometric pressure and C.-Y. Wang et al. (2018) for earth tide. Thanks to this formalism, combined to a tidal analysis of the response of a well in Martinique Island already studied in Vittecoq et al. (2020), we will be able to understand better the processes controlling hydraulic conductivity changes linked to earthquakes. To make the most of the shear modulus and transmissivity time series, we also track the evolution of effective stress in the aquifer, and numerically simulate major regional earthquake static and dynamic stresses, in order to understand why the aquifer responses to each earthquake were different.

## 2. The Martinique Galion Aquifer and Its Regional Context

Martinique is an island located in the Lesser Antilles archipelago and is formed through the subduction of the Atlantic lithosphere beneath the Caribbean plate. It is the largest volcanic island in the archipelago, covering an area of 1,080 square kilometers. The island has experienced volcanic activity, primarily andesitic, for at least 25 million years (Westercamp et al., 1990). The landscape of Martinique consists of mountainous terrain in the northern region, with the highest and youngest volcano named Montagne Pelée reaching an elevation of 1,397 m (Vittecoq et al., 2023). The climate of Martinique is characterized as humid tropical, with a rainy season occurring from July to November. The trade winds bring precipitation, which can range from 5,000 to 6,500 mm per year in the northern summits and between 1,200 and 1,500 mm per year in the south. Temperature ranges from 18 to 32 degrees Celsius in Fort-de-France. Between 2007 and 2019, the Volcanological and Seismological Observatory of Martinique recorded 128 earthquakes felt on the island, with magnitudes ranging from 2.1 to 7.4 (Figure 1). Five of them did impact the piezometric level in the studied aquifer: M7.4—**2007/11/29** (hypocentral distance—hd: 163 km), M6.5—**2014/02/18** (hd: 203 km), M5.6—**2017/02/03** (hd: 74 km), M7.3—**2018/08/21** (hd: 485 km), and M5.4—**2018/09/28** (hd: 90 km). They will be referred to by their year of occurrence in the following, with a precision of the month to distinguish the 2018 earthquakes. The Galion borehole is located on the eastern coast of the island, 1.3 km inland and situated in a small alluvial plain, at an elevation of 6 m amsl (above mean sea level).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** From Vittecoq et al. (2020): (a) Localization of the 128 felt earthquakes in Martinique reported by the OVSM-IPGP observatory during the period 2007–2019. Moment tensors of the five strongest earthquakes: M7.4–2007/11/29 (hypocentral distance—hd: 163 km), M6.5–2014/02/18 (hd: 203 km), M5.6–2017/02/03 (hd: 74 km) M7.3 2018/08/21 (hd: 485 km) and M5.4 2018/09/28 (hd: 90 km). Digital elevation model from the General Bathymetric Chart of the Ocean (gebco.net), vector of convergence, main faults, ridges and subduction trench. (b) Localization of the Galion Borehole, the ocean tide gauge (Fort-de-France harbor), the barometer station (barometric pressure data), the rain gauge station, and the seismic station of the volcanic and seismologic observatory of Martinique (OVSM-IPGP). Figure performed using ArcGis 10.5.1 (https://www.esri.com/).</figcaption></figure>

The geological formations in the area of the Galion borehole (Figure 1) consist of, from bottom to top: 15.5 Ma old hyaloclastites which form the bedrock, fractured lava flows of “Vert-pré” basalt, dated from 14 Ma ago (Westercamp et al., 1990), themselves covered by clays due to weathering, all that topped by alluvium in the form of mixed sand and clay. Borehole description is consistent with surface data (Lacquement et al., 2022). What will be designated as the “aquifer” in what follows corresponds to the fractured basalt, which can be found between 21 and 50 m below the surface around our studied borehole (Figure 2). This aquifer is semi-confined, leaking through the clay aquitard (–10 to –21 m) to an unconfined aquifer, in the alluvium (0 to –10 m), as identified in resistivity measurements (Arnaud et al., 2011) as well as with recently bored piezometers at 40 and 10 m depth (PZ40 and
<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Schematic Galion aquifer cross section. (b) Geological borehole log and description of periodic flow response to tides.</figcaption></figure>

PZ10 in the following). The main borehole reaches a depth of 50 m and is open to the aquifer from -21 to -50 m, its drilling radius ($r_w$) is 165 mm, and its casing radius ($r_c$) is 125 mm. Heliborne geophysics studies in Martinique have revealed that aquifer properties are compartmentalized (Vittecoq et al., 2015, 2019), which explains why no ocean tide load is observed, due to lower transmissivities in the aquifer between the borehole and the sea (Figure 2a). A thorough analysis of the piezometric and tidal sources signals in the frequency domain (amplitude spectrums are presented in Supporting Information S1) confirmed the negligible impact of oceanic tide given the low diurnal amplitudes. Recorded piezometric level are shown on Figure 3, exhibiting both the dynamic response of the aquifer to rainfall (presence of many rainfall peaks) and its tidal response at the diurnal and bi-diurnal frequencies (see spectrums on Figure S1 in Supporting Information S1).
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Piezometric level, recorded hourly at the Galion Site, expressed in meters above mean sea level. Surface elevation at the borehole site is 6 m amsl, so the aquifer is temporarily artesian after some heavy rainfall events. (a) General overview, vertical brown bars signal the 5 main earthquakes (b) Zoom on a month of data to visualize the tidal signal, (c) Zoom on the piezometric data in piezometers PZ40 (orange—same aquifer as the main borehole) and PZ10 (yellow—water table aquifer), as well on the rainfall recorded at La Caravelle, 10 km away from the borehole (same position as the barometer on Figure 1). No tidal signal is recorded in PZ10, but the time lag of rainfall peaks validates a leakage between the aquifer and the aquitard (Figure 2a).</figcaption></figure>

# 3. Methodology

## 3.1. Tidal Analysis

We collect in situ measurements of piezometric level and barometric pressure at a 1 hr time rate, and 1 h synthetic earth areal tidal strain data computed with the SPOTL software (Agnew, 2012). Each time series was filtered with a band pass zero-shift fourth order Butterworth filter to extract specifically the bi-diurnal frequencies. Phase shift and amplitude ratio were computed with Harmonic least squares analysis. Contributions from earth tide and barometric pressure (which occur both at the M2 and S2 frequencies) were distinguished from one another thanks to the assumption of a similar response at these two near frequencies, as in Valois et al. (2022) or Thomas et al. (2023). These calculations result in timeseries of amplitude ratio (Figure 4a) and phase lags (Figure 4b), between earth tide and piezometric level (red) and atmospheric pressure and piezometric level (green).

Frequency spectrums of piezometric levels, earth areal tidal strain, barometric pressure and ocean tide are presented in the supplementary material (Figure S1 in Supporting Information S1). Following sections will present the theoretical framework and the hydrogeological flow model needed to interpret these phases and amplitude ratios as hydrodynamic or poroelastic parameters. The earth tide response of this aquifer was already studied in Vittecoq et al. (2020). They demonstrated the impact of several earthquakes' dynamic stresses on aquifer permeability above a threshold on the magnitude/distance relationship as well as the influence of rainfall above a certain amount per days. However, the addition of (a) barometric data analysis, as well as (b) new pumping test data in boreholes in the semi-confined aquifer and in the water table aquifer enabled us to understand further the elastic properties of the aquifers and their temporal changes under stress.
<figure><img src="image_6.png" /><figcaption>Figure 4. (a) Phase lag (°) and (b) Amplitude ratio evolution over time. Amplitude ratios are given in (Pa/nstrain) for earth tide (red) and are adimensional for barometric tide (green).</figcaption></figure>

## 3.2. Poroelastic Response of a Uniform Isotropic Medium to Areal Strain and Uniaxial Load

To fully characterize an aquifer layer hydromechanically, the first step is to choose the correct set of poroelastic parameters that match available data sets. In the framework of the linear theory of poroelasticity, most of the recent literature about earth tide response (Doan et al., 2008; Roeloffs, 1996; C.-Y. Wang et al., 2018) explicitly or implicitly use the formulation

$$ \zeta = S_{\varepsilon}p + S_{\varepsilon} \cdot BK_{u}\varepsilon, \qquad (1) $$

as poroelastic constitutive equation (H. F. Wang, 2000), with $\zeta$ the increment of fluid content (adimensional), $p$ pore pressure (Pa), $\varepsilon$ earth volumetric tidal strain, $S_{\varepsilon}$ unconstrained specific storage (Pa$^{-1}$), $B$ the Skempton coefficient (–) and $K_{u}$ the undrained bulk modulus (Pa). Consequently, it raises several issues. It first raises the question of the difference between $S_{\varepsilon}$ and $S_{h}$ the uniaxial specific storage (Pa$^{-1}$), which is often not discussed, a question also raised by Bastias Espejo et al. (2022). Second the amplitude of the signal is seldom utilized in model inversion because values of the Skempton coefficient is uncertain which implies increased uncertainty on the value of $K_{u}$. In addition, $\varepsilon$ is most often not measured but deduced from the areal strain in the horizontal direction $\varepsilon_{T}$ using the free surface boundary condition (Hsieh et al., 1987; C.-Y. Wang et al., 2018), which requires to assume a value of the Poisson ratio $\nu$. All in all, most of the earth tidal response literature focuses on the phase response, which is not affected by elastic coefficient values. In cases where a significant barometric response is available in the same aquifer, we suggest to rather use the formulation depending on the horizontal areal strain $\varepsilon_{T}$, as introduced by Van Der Kamp and Gale (1983) for earth tide:

$$ \zeta = S_{h}p + S_{h} \cdot 2\gamma G \varepsilon_{T}, \qquad (2) $$

with $G$ the shear modulus (Pa) and $\gamma$ the loading efficiency (or uniaxial Skempton coefficient). This allows (a) to be free of the Poisson ratio value hypothesis to determine $\varepsilon$, as well as (b) to replace $B$ by $\gamma$, which can be determined by the barometric response, and (c) ultimately leads to the recovery of the shear modulus. The barometric response can be formulated as follows:

$$ \zeta = S_{h}p - S_{h} \cdot \gamma P_{atm}, \qquad (3) $$

where $P_{atm}$ is atmospheric pressure (Pa). Relationships 1–3 are derived in Appendix A. The only required hypothesis are $\sigma_{zz} = 0$ at the surface (Free surface boundary condition) in the case of earth tide, and $\sigma_{zz} = -P_{atm}$ at the surface; $\varepsilon_{xx} = \varepsilon_{yy} = 0$ (no horizontal strain) in the case of barometric pressure. This last hypothesis is valid if the
aquifer is sufficiently extensive laterally. We precise that these local relationships, under the previous hypothesis, are always valid, and are not bound to any “drained” or “undrained” aquifer behavior, and that the adequate storage coefficient is, as classically in hydrogeology, the uniaxial one, (and not the unconstrained one), as we can then retrieve S once given an aquifer thickness b using $S = \rho g b \cdot S_h$.

Finally, to solve for pore pressure evolution and obtain the adequate diffusion equation, we substitute Equations 2 and 3 in:

$$ \frac{\partial \zeta}{\partial t} - \frac{K_h}{\rho g} \nabla^2 p = 0, \qquad (4) $$

(H. F. Wang, 2000), with $K_h$ the aquifer hydraulic conductivity (m s⁻¹), ρ water density (kg m⁻³) and g gravity acceleration (m s⁻²).

### 3.3. Hydrogeologic Response of a Borehole-Aquifer System to Tides

Following geological observations, the aquifer was modeled by a highly transmissive layer of fractured andesite, where flow is assumed horizontal (Dupuit assumption), topped by a clay aquitard layer where flow is assumed vertical (leakage), itself topped with an unconfined aquifer where the hydraulic head is assumed constant (Figure 2b). The previous study by Vittecoq et al. (2020), chose to use the leaky aquifer model taken from C.-Y. Wang et al. (2018), which was solved in the case of earth tide. We extended it to the case of barometric loading (which changes the top boundary condition), yet we remarked that we could not explain the observed difference between earth tide and barometric phase lag with this model. We thus had to reconsider the hypothesis of a negligible aquitard storativity (Figure 2b), which lead us to a model similar to that of Odling et al. (2015), developed in the case of barometric loading, which we extended to earth tide. The main equations governing flow in the aquifer for atmospheric loading (Equation 5) and earth tide (Equation 6) are:

$$ T \nabla^2 h + q = S \left( \frac{\partial h}{\partial t} - \frac{\gamma}{\rho g} \frac{\partial P_{atm}}{\partial t} \right) \quad (\text{Atm. Pres.}) \qquad (5) $$

$$ T \nabla^2 h + q = S \left( \frac{\partial h}{\partial t} - \frac{2\gamma G}{\rho g} \frac{\partial \epsilon_T}{\partial t} \right) \quad (\text{Earth Tide}) \qquad (6) $$

With T the transmissivity (m²/s), S the hydrogeological storativity (–), h the hydraulic head (m) in the aquifer layer, and q is the leakage flow (specific discharge—m/s) from the aquitard. q is computed from the solve of the diffusion equation in the top layer, which is provided in Appendix B (Equation B1). If we define

$$ Q = \frac{b'^2 \omega}{2D'} \quad \text{and} \quad W = \frac{r_w^2 \omega}{2T}, \qquad (7) $$

where, ω is the angular frequency, b' is the thickness of the aquitard, D' its diffusivity. W quantifies the borehole storage effect: when W ≫ 1, the effect is strong with important phase lag and attenuation, and conversely when W ≪ 1, the aquifer and the borehole are “coupled”: there is no phase lag nor attenuation and the hydraulic head in the aquifer is uniform. Q quantifies leakage and more precisely how “coupled” the aquifer and aquitard are: Q ≫ 1 corresponds to decoupling and the absence of leakage, while Q ≪ 1 corresponds to high leakage and a simultaneous pore pressure response in the aquifer and the aquitard (coupling). The adimensional number ξ is defined as:

$$ \xi = 1 + i \cdot c_r \frac{W K_0(\beta r_w)}{\beta r_w K_1(\beta r_w)}, \qquad (8) $$

where $K_0$ and $K_1$ are the Bessel function of the modified Bessel functions of the second kind, respectively of the zeroth and first order and $c_r = (r_c/r_w)^2$, with $r_c$ the casing inner radius and $r_w$ the radius of the screened portion of the well (Figure 2b). The adimensional number ξ was defined by C.-Y. Wang et al. (2018) following the model from Hsieh et al. (1987) model and characterize the storage effect due to the borehole geometry. The quantity β is defined as:
$$ \beta r_w = \sqrt{W} \left( \frac{(1+i) \tanh((1+i)\sqrt{Q})}{\sqrt{Q}} S' + 2iS \right)^{\frac{1}{2}}, \quad (9) $$

where S and S' are their hydrogeological storage coefficients (–) of the aquifer and aquitard respectively. Solving for both pore pressure in the aquifer and aquitard yields, in complex notation, for the hydraulic head in the well $h_w$:

$$ \frac{\rho g h_w}{-P_{atm}} = \frac{1}{\xi} \cdot \frac{1 - \gamma + (1 - \gamma') \frac{S'}{2S\sqrt{Q}} \left( \frac{(1-i)}{\tanh((1+i)\sqrt{Q})} - \frac{(1-i)}{\sinh((1+i)\sqrt{Q})} \right)}{1 + \frac{S'}{2S\sqrt{Q}} \frac{(1-i)}{\tanh((1+i)\sqrt{Q})}} \quad (10) $$

$$ \frac{\rho g h_w}{\varepsilon} = \frac{1}{\xi} \cdot \frac{2\gamma G + 2\gamma' G' \frac{S'}{2S\sqrt{Q}} \left( \frac{(1-i)}{\tanh((1+i)\sqrt{Q})} - \frac{(1-i)}{\sinh((1+i)\sqrt{Q})} \right)}{1 + \frac{S'}{2S\sqrt{Q}} \frac{(1-i)}{\tanh((1+i)\sqrt{Q})}} \quad (11) $$

Where G, G' are respectively the shear moduli of the aquifer and aquitard (Pa), γ, γ' their respective loading efficiencies (–). Equation 10 gives the change of the hydraulic in the case of atmospheric loading, whereas Equation 11 is the response of the aquifer to earth tide. In order to invert the model analytically, several simplifications can be done and are detailed in the next section.

### 3.4. Model Calibration and Inversion: Pumping Test Data and Coseismic Piezometric Level Drop

Two recently bored piezometers in the aquifer (PZ40 at 40 m depth) and in the water table aquifer (PZ10 at 10 m depth) enabled, with a pumping test performed in 2022, to calibrate values of storativities S and S' which are poorly constrained by tidal effects (Hsieh et al., 1987; C.-Y. Wang et al., 2018). The test was interpreted with the MLU software (Carlson & Randall, 2012), a free software based on an analytical solution (Maas, 1987) adapted to layered aquifer. MLU parameter optimizer was used to fit analytically modeled curves to measured data in PZ10 and PZ40 during a 3-day long pumping test (Figure S2 in Supporting Information S1). The pumping test also gave access to a punctual measurement of transmissivity T and aquitard hydraulic conductivity K' which will not be used prior to model inversion, and thus serve as model validation.

For model inversion, obtained storativities S and S' are supposed constant. This is a fair approximation as already pointed by Hsieh et al. (1987). Leakage quantification (in the form of adimensional variable Q or vertical diffusivity D'), is based on the observation of phase lag difference between earth tide and barometric pressure $\Delta\phi = \phi_{MT} - \phi_{PA}$.

$$ \phi_{MT} - \phi_{PA} = \arg\left(1 + \frac{1 - \gamma'}{1 - \gamma} \frac{S'}{2S\sqrt{Q}} (1 - i)\right) \quad (12) $$

The variable we actually invert is $\frac{1-\gamma'}{1-\gamma} \frac{1}{\sqrt{Q}}$, which must not vary much to match the observation of almost constant $\Delta\phi$ over the 15 years of observations (within the noise range of barometric phase lag, ±5°). Thus, we hypothesize a constant ratio $\frac{1-\gamma'}{1-\gamma}$, which we calibrate at the moment of the pumping test thanks to the measurement of D' combined with phase lag measurement. All in all, we only let this diffusivity D' vary over time to match the remaining variations of $\Delta\phi$, but because $\phi_{PA}$ is quite noisy (Figure 4a), due to the low amplitude of the barometric signal in the piezometric level, we choose to average $\Delta\phi$ on each interseismic period so that we catch the possible impacts of earthquakes on the aquitard. To check the impact of this assumption, an inversion of the model with a constant D' over the whole period will also be presented. To validate further the value of D', we then compared it to
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Coseismic water level drops following the five strongest earthquakes compared to a fixed “0 m” level set 48 hr after.</figcaption></figure>

the one obtained via coseismic piezometric level drop fitting. Indeed, following each major earthquake (2007, 2014, 2017, 2018-08, and 2018-09), the piezometric level decreased exponentially of a few centimeters, for the 2014 and the two 2018 events, up to around 0.5 m for the 2007 and 2017 events. The fact that these drops are irreversible (designated as type 3 in C.-Y. Wang et al., 2004 typology) indicates that a new permanent regime, different from before the earthquake, is reached, validating the observation of variable hydraulic conductivities. We recover characteristic times $\tau$ of these exponential decays, like already done in C.-Y. Wang et al. (2004), which remarked that for a sufficiently long time after the earthquake, this can be related to hydraulic properties independently of the geometry, such that $\tau \approx \frac{4L^2}{\pi D}$, with D the diffusivity (m²/s) and L the characteristic length scale. The choice of this length scale is debatable and was chosen to be the length scale of the aquifer (1.3 km) in the previous study (Vittecoq et al., 2020). Yet inferred diffusivities (1–6 m²/s) differed by more than one order of magnitude to pumping test derived aquifer diffusivity (0.05 m²/s). What's more, changes in aquifer transmissivities, which were both increasing and decreasing, were incapable of predicting all coseismic water level variations (which consistently dropped, reflecting more complex water flow). Thus, we set the length scale to be the

thickness of the aquitard (11 m), and recover values of vertical diffusivity in the aquitard by least squares fitting of the long-term decay curves (Figure 5). Fitting starts 5 hr after earthquakes to prevent geometrical effect, and stop 35 hr after earthquakes, before tidal fluctuations become significant again. Inferred diffusivities match better with the pumping test derived diffusivity ($7 \cdot 10^{-4}$ m²/s) with only a factor 2 difference, and variations over time match the phase lag variations.

The only other required hypothesis is that the clay aquitard shear modulus $G'$ is negligible as compared to the andesitic aquifer shear modulus $G$, which is backed by previously reported values for these lithologies (H. F. Wang, 2000) and in situ S wave velocity measurements (Figure S4 in Supporting Information S1). S wave speed was measured with the Multichannel Analysis of Surface Waves (MASW) technique (Park et al., 1998), giving a S wave speed profile and thus enabling to validate the hypothesis $G' \ll G$, since $G = \rho V_s^2$ with $\rho$ the density of the rock and $V_s$ the medium S wave speed. It yields an average value of $G = 0.97$ GPa in the aquifer and $G = 0.15$ GPa in the aquitard. Given all these calibrated parameters the model (details in Appendix B) was inverted analytically into three subsequent steps: first the transmissivity was inferred from earth tide phase lag (given calibrated $S$, $S'$, and $K'$), second the loading efficiencies were inferred from the atmospheric pressure amplitude ratio (given the ratio $1-\gamma'$), and finally the shear modulus $G$ was inverted using earth tide amplitude ratio.

## 3.5. Dynamic and Static Stress Computations

To investigate the effect of earthquakes on the properties of the aquifer, static stress were computed in Vittecoq et al. (2020), using the Okada model for a point source (Okada, 1992). However, the changes of static stress were found too small to explain the change of aquifer's transmissivity, and no correlation was found between change of static-stress and change of the aquifer's transmissivity. In this study, we investigate the effects of the dynamic stress induced by the major earthquakes on the Caribbean plate. Dynamic stress was calculated synthetically, since no direct observation exists. The equation of motion in a 3D continuum medium is written between velocity component $v_i(i = 1,2,3)$ and stress component $\sigma_{ij}(i, j = 1,2,3)$ in Cartesian coordinate $(O - x_1x_2x_3)$ as

$$ \rho \frac{\partial v_i}{\partial t} = \frac{\partial}{\partial x_j} \sigma_{ij} + F_i \quad (13) $$

where $\rho$ is medium density and $F_i$ is a body force. Here, the summation convention about the same index is assumed for compactly presenting the equations. The isotropic, elastic medium is assumed, as widely accepted to analyze the seismic wave propagation (e.g., Aki & Richards, 2002), and Hooke's law governs the constitutive relation between stress and strain $\epsilon_{ij} = (\partial u_i / \partial x_j + \partial u_j / \partial x_i)/2$ with two elastic coefficients (Lamé constants) $\lambda$ and $\mu$:

$$ \sigma_{ij} = c_{ijkl} \epsilon_{kl} \text{ where } c_{ijkl} = \lambda \delta_{ij} \delta_{kl} + \mu (\delta_{ik} \delta_{jl} + \delta_{il} \delta_{jk}). $$
Note that $\delta_{ij} = 1$ if $i = j$, otherwise $\delta_{ij} = 0$. In this study, we adopt a finite difference method (Aochi & Madariaga, 2003a) to calculate the seismic wave propagation from a given earthquake source to the site location of interest. As the above equations are differential equations in time and space, it is easy to apply the framework of finite difference method, especially when taking a staggered grid scheme (Virieux & Madariaga, 1982). Grid size $\Delta s$ and time step $\Delta t$ are required to satisfy the stability condition $\Delta t \le c \cdot \frac{\Delta s}{V_{max}}$, where $V_{max}$ is the maximum wave velocity in the medium. As the earthquakes of interest (Figure 1) are largely different in space, we prepare a model dimension of 100 × 200 × 150 km with a grid size of $\Delta s = 500$ m and $\Delta t = 0.01$ s. We use a known 1D structure of the region (Dorel, 1978), consisting of 10 layers supposing the flat ground surface without sea layer. The magnitudes of the earthquakes vary from 5.4 to 7.4. For simplicity, we suppose a point source, as the distance between source to site is sufficiently larger than the dimension of the finite source area. We then give a source time function STF(t), equivalent to the moment magnitude of each earthquake, whose bell-shape form is characterized by a duration T (Aochi et al., 2013).

$$T = \frac{L}{V_R} \qquad (14)$$

where L is the fault length, which is scaled with magnitude (Wells & Coppersmith, 1994) and we assume a typical rupture velocity $V_R = 2.2$ km/s. For example, for the 2007 M7.4 earthquake, we assume $L = 69$ km and then $T = 31.4$ s. This source time function is expressed as the seismic moment tensor $M_{ij}(t) = f(focalmechanism) \cdot STF(t)$, and this is decomposed as a combination of body force $F_i(t)$ for Equation 13 (Graves, 1996). The function $f(focalmechanism)$ is a coefficient for each component of tensor in function of the fault strike, dip and rake. The simulations of the seismic wave propagation though the finite difference scheme in the 3D medium for the given sources and stations are practical as we get the direct output of three components of velocity and six components of stress field at any point. We extract the values of stress components at the site position below the ground surface by one or a half grid, at a depth of 500 or 250 m. The boundary condition of the free surface leads to the zero-stress component concerning the vertical axis, namely $\sigma_{ij} = 0$ ($i = 3$) at the ground surface. Six components of the stress tensor allow to compute the volumetric stress and the second deviatoric stress invariant. Note that we observe no significant difference on the Martinique Island, as there is no complexity in model parameters of earthquake source and medium. We thus consider that the estimation is a regionally typical incident wave to the site of interest.

# 4. Results

## 4.1. Model Inversion Results

Transmissivity and shear modulus evolution were reconstructed over the 2007–2022 period (Figure 6). Transmissivity evolves between a minimum of $2 \cdot 10^{-5}$ m$^2$ s$^{-1}$ just after the November 2007 earthquake, and a maximum around $3 \cdot 10^{-4}$ m$^2$ s$^{-1}$ in the 2014–2017 interval. Shear modulus evolves between 1.1 GPa before 2007 earthquake, to a maximum around 1.8 GPa in 2013. Between these extremes, most changes occur at the date of major earthquakes some appearing reversible over a 6–7 months period, like in 2007 and 2017, some irreversible like in 2014, when the stabilized value differs from the previously measured values. All these results are validated by numerous methods which results also appear on Figure 6: yellow circles indicate values of transmissivity and vertical diffusivity inferred from pumping tests, and the purple circle is the value of shear modulus computed with the MASW technique. The value given by MASW is twice lower than what we measure, yet the order of magnitude of 1–2 GPa, very low for basalts, even highly fractured, is coherent.

Two features appear strikingly in the results, and will be discussed thoroughly in Sections 5.1 and 5.2. First the anti-correlation of variations of transmissivity and shear modulus, in terms of coseismic variations (2007 and 2014) but also in interseismic periods. When transmissivity increases, the shear modulus decreases and conversely, which is expected in a fractured medium. Second, the differences in the response of transmissivity to earthquakes: in 2007 & 2017 on the one hand, it appears the impacts on transmissivity were reversible with a recovery over a 6–12 months period. On the other hand, in 2014, the change in transmissivity appears to be irreversible with little variations after the earthquake. These differences cannot be explained by simple considerations of magnitude or peak ground velocity (PGV) measurements. Regarding vertical diffusivity $D'$ which was let “free” in the inversion, its variations are far less significant, its values being concentrated between $3 \cdot 10^{-4}$ and
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Results of model inversion (blue dots) in terms of aquifer transmissivity T and aquifer shear modulus G. Both correspond to the fractured andesite layer. Yellow circles indicate results from pumping test analysis, and the purple circle MASW-derived shear modulus.</figcaption></figure>

$6 \cdot 10^{-4} \text{ m}^2\text{/s}$. Consequently, assuming a constant vertical diffusivity over the whole interval only causes minor changes to the results, as can be seen on Figure S3 in Supporting Information S1.

*Uncertainty* increases for transmissivity after 2014 while it decreases for shear modulus. It is not linked to a difference in the input data but in the range of values and model inversion. Indeed, after the 2014 earthquake, transmissivity is near the threshold where it implies no borehole storage effect, thus the model inversion uncertainty rises. On the contrary, concerning shear modulus, the value is more uncertain before 2014, where the attenuation of the borehole storage effect adds up to the amplitude uncertainty. After 2014, the borehole storage effect attenuation fades and there only remain uncertainty on the input amplitude, which explains why the shear modulus is better constrained after the 2014 earthquake. One should also consider the uncertainty introduced by the calibration procedure and the hypotheses required by the inversion process (mainly the assumption of constant storativities *S* and *S*'). Given the relative sensitivity of the modeled phase lag to transmissivity and storativity, the calibration may introduce a systematic error on transmissivity no more than half the error on storativity (in logarithmic scale). This systematic error does not challenge the implications discussed later which concentrate on the relative changes.

## 4.2. Dynamic Stresses Model Results

From the stress tensor calculations described in Section 3.5, noting its principal components $\sigma_{n=1,2,3}$, we first extract two stress metrics: volumetric stress, or hydrostatic stress $\sigma_{kk} = \frac{\sigma_1 + \sigma_2 + \sigma_3}{3}$ (Figure 7a) and the second deviatoric stress invariant $J_2 = \sqrt{\frac{(\sigma_1 - \sigma_2)^2 + (\sigma_2 - \sigma_3)^2 + (\sigma_1 - \sigma_3)^2}{6}}$ (Figure 7b), which quantifies mean shear stress (without preferential directions). Volumetric stress is relevant because it could pressurize water in the aquifer and possibly cause rapid water flow and unclog fractures. Shear stress is also important because of its potential to create microfractures in the aquifer, or, in the presence of rough fractures, modify their aperture and therefore their transmissivity. To make the most of the simulation results, a survey of preferential fracture planes appearing in nearby outcrops was performed (see namely Arnaud et al., 2011; Lacquement et al., 2022). Two principal direction of planes stand out with respective N20° and N60° strikes (Figure 8), dominated by high dips (average >80°).

Given this information we computed the shear stresses amplitude relative to the identified fault planes N20° and N60° and plotted them on Figures 7c and 7d. Normal stresses to the fault planes were not plotted since they did not show significant differences with the volumetric stresses (Figure 7a). The earthquakes main characteristics and peak stress values are listed in Table 2, along with their impact on transmissivity and shear modulus.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Results of earthquake dynamic stresses simulation. Noting the stress tensor components $\sigma_{ij}$ and its principal components $\sigma_{n=1,2,3}$, plots show: (a) Volumetric stress $\sigma_{kk} = \frac{\sigma_1 + \sigma_2 + \sigma_3}{3}$; (b) Second deviatoric stress invariant</figcaption></figure>

$$J_2 = \sqrt{\frac{(\sigma_1 - \sigma_2)^2 + (\sigma_2 - \sigma_3)^2 + (\sigma_1 - \sigma_3)^2}{6}}; \quad (c) \text{Horizontal shear stress amplitude in the N20° fault plane } \|\vec{\tau}_{\text{N20}}\|; \quad (d) \text{Horizontal shear stress amplitude in the N60° fault plane } \|\vec{\tau}_{\text{N60}}\|.$$

As expected, earthquake dynamic stresses differ in frequency, with lowest dominant frequency for highest magnitudes, but the highest values of stresses are not reached for the strongest earthquake (2007). The 2014 earthquake, while weaker and located further away than the 2007, generated more stress because of its shallower hypocenter which allowed better generation of surface waves. It impacted the simulated stress at the surface more, especially regarding shear stresses. The same effect remains true for the 2017 and 2018/09 earthquakes, which generated significantly lower peak volumetric stresses but comparable peak shear stresses than the 2007 one, especially in the identified principal fault planes.

**5. Discussion**

In this section, we interpret the changes in transmissivity and shear modulus reported previously on a mechanical point of view. These interpretations imply the fractured nature of the aquifer. However, the hydrogeological model we used for tidal analysis assumed a porous media like flow. This difference is justified as long as the R.E.V. (representative elementary volume) scale of our aquifer is small as compared to the characteristic length of pore pressure diffusion around the borehole. Dimensional analysis reveals this characteristic length can be estimated to be $\frac{1}{|\vec{\rho}|}$ where $\beta$ is defined in Equation 9. Given the values given in Table 1, $\frac{1}{|\vec{\rho}|} \approx 12$ m for the Galion borehole. On the other hand, the scale of the fracture network is supposed smaller than the meter scale, meaning a R.E.V at a scale less than $1 \text{ m}^3$, which is coherent with the high fracture density required to explain the low shear modulus measured. Thus, flow at the scale of the aquifer can be approximated as a porous media flow, even if the porosity is the result of a fracture network. This porous equivalent medium representation is expected to be valid only for small deviations around a particular state of stress, strain and pore pressure (a hypothesis similar as in
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Polar diagram of the azimuthal and dip angles distributions of fractures in outcrops in the vicinity of the Galion borehole.</figcaption></figure>

linear poroelasticity). Typically, in the case of tidal fluctuations, we expect them to affect fracture parameters such as aperture negligibly, justifying the porous model for tidal response. On the contrary, under large stress or pore pressure changes (typically with earthquakes or large recharge events due to rainy season rainfall), the fractures parameters change significantly and affect the hydraulic behavior of the aquifer.

## 5.1. Reversible Changes: Transmissivity and Shear Modulus Dependence to Effective Stress

The correlation between transmissivity and pore pressure during inter-seismic periods is striking at the Galion borehole. Periods of low pore pressure after the 2007 and 2017 earthquakes are concurrent with a low transmissivity, and the recovery of transmissivity to the pre-earthquake value is concurrent with the increase of pore pressure with the rainy season (Figure 9).

Models exist to describe permeability evolution with effective pressure in a porous media, like namely one proposed by Walsh (1981), adapted to fractured media. It links fracture permeability change to effective pressure in a fractured media through the relationship

$$ \left(\frac{k}{k_0}\right)^{\frac{1}{n}} = 1 - \frac{\sqrt{2h}}{a_0} \ln\left(\frac{\sigma^*}{\sigma_0^*}\right) \quad (15) $$

Where *k* is the permeability of the fractured medium at the effective stress σ*, the subscript 0 designating a reference state, *a*₀ mean fracture aperture in the reference state and *h* the root mean square value of the fracture surface elevation, quantifying fracture rugosity. Because transmissivity is proportional to permeability, we can apply this model, given σ* = σ - αp, with σ mean stress and p aquifer mean pore pressure, where we assume Biot Willis coefficient α ≅ 1, and that the static stress variations are negligible in comparison with pore pressure variations (which is justified by Okada model calculations), Δσ ≪ Δp. If we work on an interseismic period, like the 2007–2014 interval, to avoid any influence of dynamic stresses or irreversible effects, we infer by least square fit a value of h/a₀ of approximately 2.8, describing a rough fracture surface (Figure 10a).

This finding is interesting when we expand the found relationship between transmissivity and piezometric level to our complete time series. We obtain an artificial transmissivity de-influenced from pore pressure (Figure 10b). If the absolute value of this variable has no physical meaning, its variation reveals the remaining factors that control the transmissivity. Most coseismic variations of transmissivity seem explained by piezometric level change, particularly post 2007- and 2017-earthquakes recoveries. Piezometric level drops, most marked in 2007 and 2017 are thus probably the first order cause of the decrease in aquifer transmissivity through effective stress increase. This observation is unprecedented, possibly because the study site specificities enable this strong sensitivity of

<table><caption>**Table 1**
List of Parameters Used in the Model Calibration or Inversion</caption><thead><tr><th>Parameter (unit)—(source)</th><th>Aquitard</th><th>Aquifer</th><th>Use</th></tr></thead><tbody><tr><td>Thickness (m)—(drilling data)</td><td>11</td><td>29</td><td></td></tr><tr><td>Hydraulic conductivity (m/s)—(pumping)</td><td>1.3 · 10<sup>-6</sup></td><td>/</td><td>Calibration</td></tr><tr><td>Transmissivity (m<sup>2</sup>/s)—(pumping)</td><td>/</td><td>1.6 · 10<sup>-4</sup></td><td>Validation</td></tr><tr><td>Diffusivity (m<sup>2</sup>/s)—(calibrated Figure 5)</td><td>1.4–3.6 · 10<sup>-4</sup></td><td>/</td><td>Calibration</td></tr><tr><td>Storativity S (–)—(pumping)</td><td>0.0354</td><td>0.0035</td><td>Calibration</td></tr><tr><td>Loading efficiency ratio <sup>1-γ'</n>/<sub>1-γ</sub> (–)—(calibrated)</td><td>/</td><td>0.76</td><td>Calibration</td></tr><tr><td>Shear modulus G(GPa)—(masw)</td><td>0.15</td><td>0.97</td><td>Calibration (aquitard) and validation (aquifer)</td></tr></tbody></table>

*Note.* Hydraulic conductivity, storativities and transmissivity are results of parameters optimization by the MLU software on the 2022 pumping test. Shear Modulus was measured with the MASW technique.
<table><caption>**Table 2**

*Compilation of Earthquakes Main Features*</caption><thead><tr><th>Year</th><th>Mag.</th><th>Depth (km)</th><th>PGV (cm/s)</th><th>Dominant frequency (Hz)</th><th>Peak σ<sub>kk</sub> (kPa)</th><th>Peak J<sub>2</sub> (kPa)</th><th>ΔT/T (%)</th><th>ΔT<sub>d</sub>/T<sub>d</sub> (%)</th><th>de-influenced</th><th>ΔG/G (%)</th></tr></thead><tbody><tr><td>2007</td><td>7.4</td><td>160.0</td><td>0.71</td><td>0.02</td><td>5.1</td><td>8.1</td><td>-66%</td><td>+20%</td><td></td><td>+17%</td></tr><tr><td>2014</td><td>6.5</td><td>15.0</td><td>0.38</td><td>0.1</td><td>8.1</td><td>15.9</td><td>+110%</td><td>+233%</td><td></td><td>-25%</td></tr><tr><td>2017</td><td>5.6</td><td>37.8</td><td>0.11</td><td>0.4-0.5</td><td>1.7</td><td>6.0</td><td>-83%</td><td>-55%</td><td></td><td>NSC</td></tr><tr><td>2018/09</td><td>5.4</td><td>48.8</td><td>0.059</td><td>0.5-0.7</td><td>1.1</td><td>2.9</td><td>NSC</td><td>+17%</td><td></td><td>-8%</td></tr></tbody></table>

*Note.* Magnitude and localization and PGV from of the volcanic and seismologic observatory of Martinique (OVSM-IPGP) and its seismic station (Figure 1). Main frequencies and peak stresses from finite element model simulation. Transmissivity and shear modulus relative variations from results of model inversion (Figure 6). For “De-influenced transmissivity,” see Section 5.1 and Figure 10. “NSC”: no significant change.

permeability (highly fractured medium with high compliance, and shallow depth which ensure effective stress relative variations are highly sensitive to pore pressure). The notable exceptions are the 2014 and 2017 earthquakes, which caused irreversible changes to aquifer transmissivity, which cannot be explained by pore pressure variations.

In terms of physical mechanism, this observation leads us to revisit the previous hypothesis of the influence of heavy rainfall events on aquifer properties. The previously envisioned process implied “piston stroke effects” (Vittecoq et al., 2020) following rainfall events, because the main hypothesis in the literature considered clogging/unclogging effects. If they are not ruled out by our observation, such clogging and piston stroke effects are not required because variations could be explained totally by mechanical effects. In other words, rainfall does have an impact on transmissivity but probably only through mean piezometric level, which does not require extreme events but is rather sensitive to the dry/wet changing seasons in the case of the Galion aquifer.

Adopting another point of view, we are able to estimate mean fracture compliance in our aquifer. In a isotropic, fractured medium with randomly oriented cracks, Gueguen and Dienes (1989) modeled the dependence of permeability *k* to crack density *ρ*, mean crack aperture *w*, mean crack length *L*, and connectivity *f*:

$$k = \frac{2}{15} \frac{f}{L} w^3 \rho \propto w^3 \rho \qquad (16)$$

Assuming no crack density or connectivity change during interseismic intervals, we can recover the evolution of mean crack aperture and plot it against pore pressure similarly as in Figure 10a. The slope of the inferred relationship between pore pressure (which results as an extensive normal stress on fractures) and hydraulic fracture aperture constitutes a measurement of mean fracture normal compliance $S_n = \frac{\delta w}{\delta \sigma_n} = 1 \cdot 10^{-8}$ m/Pa. This value is among the highest reported in the literature, as namely in Barbosa, Caspari, et al. (2019), and far higher than those measured in laboratory experiments (Gallagher et al., 2022). Yet it is coherent with what is expected for our site, since fracture compliance scales linearly with the size of fractures (Worthington & Lubbe, 2007) and decrease with effective stress (Brown & Scholz, 1986). Furthermore, because we compute it through effective pressure on long time scales, it corresponds to a drained fracture compliance, while in the literature the technique used is

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Transmissivity and piezometric level variations following the 2007 and 2017 earthquakes (vertical dark gray lines).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** (a) Walsh (1981) Model applied to our data on the interseismic interval 2007–2014. *T*₀ is fitted such that *T*/*T*₀ lie below 1 and *σ*⁰* corresponds to effective stress when T = T₀. (b) De-influenced transmissivity time series compared to the transmissivity.</figcaption></figure>

seismic attenuation, effectively giving undrained fracture compliances, which are expected to be one order of magnitude lower (Worthington & Lubbe, 2007).

Shear modulus interseismic variations are less obvious to analyze, as they are less closely correlated to piezometric level, yet it is notable that post-earthquake variations are always in the direction of a recovery toward pre-earthquake values, without reaching complete recovery. The long-term variations match with variations of transmissivity, being anti correlated, which can also be explained by piezometric level fluctuations (opening fractures tends to decrease the shear modulus -and increases transmissivity-while closing fractures increases the shear modulus and decreases permeability). If shear modulus variations are more complex than transmissivity, it is probably because it is not only sensitive to connected fractures and effective stress but to the state of damage in the whole rock.

**5.2. Earthquake Dynamic Stresses Impact on Transmissivity and Shear Modulus**

An unexpected result is the stronger impact on aquifer properties of the 2014 and 2017 earthquakes compared to the 2007 one, while it showed lower PGVs (peak ground velocity) and seismic energy densities more than 1 order of magnitude smaller. Seismic energy density, as proposed by C.-Y. Wang and Manga (2010), is well suited to predict piezometric level fluctuations, especially those linked to static stress, as it decays with the cube of the
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Relative intensities of peak stresses (3 first left bar groups) and property change (2 last bar groups) between the main earthquakes (2007, 2014, 2017, and 2018/09) (Normalized by the strongest, i.e. 2014 earthquake). "Peak shear stress on fracture plane" is the mean of peak shear stresses on fracture planes of azimuth N20° and N60°.</figcaption></figure>

distance to the earthquake. Yet permeability changes, away from the near field, are recognized to be linked to dynamic stresses, and thus PGV is widely used in the literature as a predictor to permeability change. Comparing results of our stress simulation and de-influenced transmissivity changes (Table 2), it appears that the best predictor of transmissivity variation is peak shear stress, since it matches best the relative amplitude of transmissivity changes (Figure 11). It is coherent with previous literature which uses PGV as a proxy for stress, yet our simulations show that one should be cautious when using PGV, as it may not always be reliable depending on local lithologies or influenced by site effects. Indeed, in our case PGV was measured from a seismic station 20 km away, on the other side of Martinique. Shear stress applied on rough fractures are expected to change its transmissive properties even more than normal stresses. Numerical simulations of fracture networks even show that at low stress, it is possible to link linearly fracture permeability and shear stresses (Yang et al., 2022), which is consistent with our observations where transmissivity changes scale linearly with peak dynamic shear stresses.

In terms of magnitude, it is interesting to note that peak dynamic stresses are comparable to effective stresses variations ($1 m_{H2O} = 10$ kPa), while changes in transmissivity where also comparable in magnitude between static effective stress effects and dynamic stress effects. It seems to imply a similar magnitudes of stresses impact transmissivity similarly, yet this consideration is possibly misleading as processes are very different, namely in terms of stress rates (which peak at $10 \text{ kPa s}^{-1}$ for earthquake dynamic stresses vs. $8 \cdot 10^{-4} \text{ kPa s}^{-1}$ for piezometric fluctuations).

A particularly interesting behavior is seen following the 2014 earthquake, an earthquake which did trigger a very small coseismic piezometric level drop yet the biggest change in transmissivity and shear modulus. Transmissivity does not vary or recover after the earthquake contrary to others, while shear modulus partially recovers over 10–12 months. We can explain this by adopting the same point of view as Illien et al. (2022) who analyzed time series of S-wave velocity variations $\delta v/v$, which is equivalent to our time series, considering that $G = \rho v^2$, thus $\frac{\delta G}{G} = 2 \frac{\delta v}{v}$. He decomposed $\delta v$ as two components, one hydrological contribution $\delta v_h$, linked to effects we mentioned above, and a mechanical contribution which he attributes to nonlinear mesoscopic elasticity $\delta v_{NLME}$, that is, a coseismic drop in velocity/elastic moduli, followed by an exponential recovery, a phenomenon which is observed on time scales of seconds to years, also called relaxation or healing. For the sake of simplicity, it is possible to write $\delta v = \delta v_h + \delta v_{NLME}$, yet reality is probably more complex and effects entangled. In our case, if we assume the hydrologic part $\delta G_h$ is linked to variations of transmissivity, we can isolate $\delta G_{NLME}$, which is the varying part. The best fit on $\delta G_{NLME}$ recovery yields a 20% coseismic drop (corresponding to a 10% drop in velocity) and a characteristic time of recovery of 116 days, values coherent with literature observations (Brenguier et al., 2008; Illien et al., 2022, 2023; Richter et al., 2014; Viens et al., 2018), which is satisfying given the difference in measurement technique and the specificity of our study site.
### 5.3. Estimation of Crack Density Variations

To interpret the value of the shear modulus variation in time in terms of fracture-scale processes occurring during and between earthquakes, we used the differential effective medium theory developed by David and Zimmerman (2011) which models the elastic properties of a rock sample containing spheroidal pores. It establishes relationships between pore properties (crack aspect ratio $\xi$ and crack density $\rho$) and the effective elastic properties of a rock. It is reasonable to assume the average aspect ratio in our aquifer $\xi \le 10^{-3}$, especially given the high compliance we computed in Section 5.1. In this case, the shear modulus only depends on crack density and not on the aspect ratio anymore. The semi-analytical relationship provided is cumbersome, but it is possible from it to approximate a relationship between the relative variations of crack density $\rho$ and shear modulus $G$, valid in our estimated range of parameters $\xi \le 10^{-3}$, $0.2 < \rho < 0.6$:

$$ \frac{\delta\rho}{\rho} \approx (-3 \pm 0.5) \frac{\delta G}{G} \qquad (17) $$

Because effective pressure changed the aperture of fracture, this analysis is only applicable for the 2014 earthquake, when we could separate hydrodynamic effects from presumable microcracking effects. As we obtained in previous section an estimated relative change in shear modulus by 20%, this last equation implies the crack density in our medium could have increased temporarily by as much as 60%(±10%), before retrieving partly its previous level, following shear modulus evolution. If we assume a crack density of $\rho = 0.2$, it means an increase - following the earthquake-to $\rho = 0.32$. As extreme as it may appear, it is reasonable when compared to the recorded increase in aquifer transmissivity from $9 \cdot 10^{-5} \text{ m}^2/\text{s}$ to $2 \cdot 10^{-4} \text{ m}^2/\text{s}$.

### 6. Conclusion

Thanks to an adequate theoretical poroelastic formulation we were able to bring to light the complete potential of groundwater level tidal analysis. Such analysis requires a relatively good knowledge of the study site, but most importantly a relative equilibrium in the amplitude of the response to different sources, here barometric pressure and earth tidal strain. Thus, it may be only applicable in intermediary sites, in between deep, confined aquifers which are expected to respond mostly to earth tide, and shallow, unconfined aquifer which generally tend to show a stronger response to barometric loading (e.g., in Thomas et al. (2023)). Yet, when such a favorable site is available, and if its geometry is constrained enough to easily model groundwater pore pressure diffusion, this technique has a tremendous potential in characterizing aquifers, both in terms of hydrodynamic properties (permeability) and poroelastic ones (shear modulus). Its passive character (it only requires hourly groundwater level monitoring in borehole and atmospheric pressure recordings) makes it cheap and promising as it produces results otherwise difficult to obtain (with invasive, costly techniques such as pumping tests or seismic surveys). One last parameter remains difficult to constrain accurately without additional hypothesis, the storage coefficient. It is a key parameter for resource managers, but has low impact on tidal response and thus may still require the use of other techniques (classical pumping tests or borehole teleseismic signal analysis).

Our results, validated with classical techniques, are comparable to those of tidal analyses carried out in previous studies with regard to the variable permeability of the aquifer, and to those of seismic velocity studies with regard to their variations. We reveal the high sensitivity of parameters to effective pressure in shallow aquifers, as well as the importance of shear stresses to disrupt aquifer properties. In comparable shallow aquifers, earthquake depths seem to be important to favor stronger responses of aquifer to earthquakes, through the amount of shear stresses generated. In terms of physical mechanism, if it doesn't rule out fracture clogging or unclogging, the main hypothesis in the literature in terms of occurrence, this study shows that such processes are not necessary. Purely elastic processes are sufficient to explain most of the strong variations observed. Processes evoked here, that is, the sensitivity to effective stress of field scale media, the potential of shear stress to disrupt existing fractures, as well as provoke microcracking are probably undermined, yet more surveys and comparison with other locations are necessary to assess the applicability of these results in other aquifer lithologies. Our results also question the choice of PGV to relate earthquake impact to shallow aquifer, as our results tend to favor peak stress as a better predictor of permeability changes. Indeed, the only earthquake that provoked irreversible changes to our aquifer properties is not the one with highest PGV but with highest peak stress.
# Appendix A: Poroelasticity Theory Framework

## A1. Notations and Constitutive Relations

$\zeta$ increment of fluid content (-)

$p$ pore pressure (Pa)

$\epsilon = \epsilon_{11} + \epsilon_{22} + \epsilon_{33}$ volumetric strain (positive for extension) (-)

$\sigma = \frac{\sigma_{11} + \sigma_{22} + \sigma_{33}}{3}$ mean stress (Pa, positive for tensile stress)

$S_{\epsilon} = \frac{\delta\epsilon}{\partial p|_{\epsilon=0}}$ constrained specific storage (Pa$^{-1}$)

$\alpha = \frac{\delta\zeta}{\partial\epsilon|_{p=0}}$ Biot-Willis coefficient (-)

$B = \frac{\partial p}{\partial\sigma|_{\zeta=0}}$ Skempton coefficient (-)

$K = \frac{\delta\sigma}{\delta\epsilon|_{p=0}}$ drained bulk modulus (Pa)

$G = \frac{\delta\sigma_{ij}}{2\delta\epsilon_{ij}|_{i\neq j}}$ shear modulus (Pa)

$K_u = \frac{\delta\sigma}{\delta\epsilon|_{\zeta=0}}$ undrained Bulk modulus (Pa)

$\delta_{ij}$ Kronecker delta

$K_v = K + \frac{4}{3}G$ drained uniaxial bulk modulus (Pa)

$K_v^{(u)} = K_u + \frac{4}{3}G$ undrained uniaxial bulk modulus/P-wave modulus (Pa)

$\gamma = \frac{\partial p}{\partial\sigma_{33}|_{\zeta=\epsilon_{11}=\epsilon_{22}=0}}$ loading efficiency/Uniaxial Skempton coefficient (-)

$S_h = \frac{\delta\epsilon}{\partial p|_{\epsilon_{11}=\epsilon_{22}=\sigma_{33}=0}}$ uniaxial specific storage (Pa$^{-1}$)

$k$ permeability (m$^2$)

$K_h = \frac{\rho gk}{\mu}$ hydraulic conductivity (m s$^{-1}$)

Poroelasticity constitutive relations:

(Mixed Stiffness formulation—H. F. Wang (2000) eq 2.25 or 2.50)

$$ \zeta = S_{\epsilon}p + \alpha\epsilon \quad (A1) $$

(H. F. Wang, 2000) eq 2.53

$$ \sigma_{ij} = \left( K_u - \frac{2}{3}G \right) \delta_{ij} \epsilon + 2G\epsilon_{ij} - BK_u \zeta \delta_{ij} \quad (A2) $$

Yields with Equation A1

$$ \sigma_{ij} = \left( K_u - \frac{2}{3}G - BK_u \alpha \right) \delta_{ij} \epsilon + 2G\epsilon_{ij} - \alpha p \delta_{ij} \quad (A3) $$

## A2. Earth Tide

We want to express our constitutive relation as a function of the areal tidal strain which is known from theoretical tide potential $\epsilon_T = \epsilon_{11} + \epsilon_{22}$. $\epsilon_{33}$ is unknown but is constrained by the free surface boundary condition as $\sigma_{33} = 0$:

Equation A3 with $i = j = 3$ and $\sigma_{33} = 0$ yields:

$$ \left( K_u - \frac{2}{3}G - \alpha BK_u \right) \epsilon + 2G(\epsilon - \epsilon_T) - \alpha p = 0 \quad (A4) $$
Given that $BK_u = \gamma K_v^{(u)}$ (H. F. Wang, 2000), eq 3.79, we obtain:

$$ \varepsilon \left( K_u + \frac{4}{3} G \right) (1 - \alpha\gamma) = 2G\varepsilon_T + \alpha p \quad (A5) $$

We let $\tau = 1 - \alpha\gamma$, and we get:

$$ \varepsilon = \frac{2G}{\tau K_v^{(u)} \varepsilon_T} + \frac{\alpha}{\tau K_v^{(u)} p} \quad (A6) $$

Injecting Equation A6 in Equation A1 yields:

$$ \zeta = S_{\varepsilon} p + \alpha \left( \frac{2G}{\tau K_{v}^{(u)}} \varepsilon_{T} + \frac{\alpha}{\tau K_{v}^{(u)}} p \right) \quad (A7) $$

that is:

$$ \zeta = \left( S_{\varepsilon} + \frac{\alpha^2}{\tau K_{v}^{(u)}} \right) p + \frac{2G\alpha}{\tau K_{v}^{(u)} \varepsilon_{T}} \quad (A8) $$

We notice in that

$$ \frac{\alpha}{K_u + \frac{4}{3}G} = \frac{\alpha}{K_v^{(u)}} = \gamma S_{\varepsilon} \quad (A9) $$

H. F. Wang (2000) eq 3.56 and 3.85

yielding

$$ \zeta = S_{\varepsilon} \left(1 + \frac{\alpha\gamma}{\tau}\right) p + \frac{2G\gamma S_{\varepsilon}}{\tau} \varepsilon_{T} \quad (A10) $$

Note that the “classical” uniaxial hydrogeologic storage $S$ is related to $S_{\varepsilon}$ by

$$ \frac{S_{\varepsilon}}{S} = \frac{K_v}{K_v^{(u)}} = 1 - \alpha\gamma = \tau \quad (A11) $$

(H. F. Wang, 2000) eq 3.54 and 3.81).

And that $(1 + \frac{\alpha\gamma}{\tau}) = \frac{1 - \alpha\gamma + \alpha\gamma}{\tau} = \frac{1}{\tau}$ finally giving:

$$ \zeta = Sp + 2G\gamma S\varepsilon_T \quad (A12) $$

Eventually, to obtain the diffusion equation verified by pore pressure, we use (H. F. Wang (2000) 4.62):

$$ \frac{\partial \zeta}{\partial t} - \frac{k}{\mu} \nabla^2 p = Q \quad (A13) $$

Which finally gives:

$$ S \frac{\partial p}{\partial t} + 2G\gamma S \frac{\partial \varepsilon_T}{\partial t} - \frac{k}{\mu} \nabla^2 p = Q \quad (A14) $$

Thus, generalizing the result of Van Der Kamp and Gale (1983). Note that with the convention $\varepsilon_T > 0$ for compression sometimes used in the literature, we need to add a minus sign before $\varepsilon_T$.
### A3. Atmospheric Pressure

In the case of atmospheric pressure, the boundary conditions are the following:

$\sigma_{33} = -P_{atm}; \varepsilon_{11} = \varepsilon_{22} = 0$ (aquifer laterally extensive).

Thus $\varepsilon = \varepsilon_{33}$ and Equation A3 with $i = j = 3$ yields:

$$ -P_{atm} = \left( K_u - \frac{2}{3}G - \alpha BK_u + 2G \right) \varepsilon - \alpha p \quad (A15) $$

that is:

$$ -P_{atm} = \tau K_v^{(u)} \varepsilon - \alpha p \quad (A16) $$

We use Equation A9 again to get:

$$ \varepsilon = \frac{\gamma S_{\varepsilon}}{\tau \alpha} (\alpha p - P_{atm}) \quad (A17) $$

Which we inject in A1 to get

$$ \zeta = S_{\varepsilon} \left( 1 + \frac{\alpha \gamma}{\tau} \right) p - \frac{S_{\varepsilon} \gamma}{\tau} P_{atm} \quad (A18) $$

Which yields as between A10 and A12:

$$ \zeta = Sp - S\gamma P_{atm} \quad (A19) $$

Which gives in terms of diffusion equation, using A13:

$$ S\left(\frac{\partial p}{\partial t} - \gamma \frac{\partial P_{atm}}{\partial t}\right) - \frac{k}{\mu} \nabla^2 p = Q \quad (A20) $$

Which is similar to A14. This result was already established, namely in Roeloffs (1996).

## Appendix B: Pore Pressure Diffusion Model

### B1. Model Derivation

Model derivation is inspired from Odling et al. (2015); Thomas et al. (2023); and C.-Y. Wang et al. (2018). The present model is intermediate between the three as it matches the boundary conditions in Wang et al. while it considers a coefficient $S' > 0$, like in Thomas et al. and Odling et al.

The diffusion equations are solved separately in the aquifer and the aquitard, coupled by the boundary conditions. In the aquitard, pore pressure $p$ follows:

$$ \frac{\partial^2 p}{\partial z^2} = \frac{S'}{K' b'} \left( \frac{\partial p}{\partial t} - \frac{\partial}{\partial t'} \right) \quad (B1) $$

Where $f'$ is the tidal forcing in the aquitard ($\gamma' P_{atm}$ for atmospheric pressure or $2\gamma'G'\varepsilon_T$ for earth tide). Boundary conditions are:

- At $z = 0$ (at the aquifer/aquitard boundary), $p = \rho gh$, with $h$ the hydraulic head in the aquifer

- At $z = b'$ (at the water table aquifer boundary), $p = p_{top}$, with $p_{top} = P_{atm}$ (Atmospheric pressure) or $p_{top} = 0$ (Earth tide).

In the aquifer, in cylindrical coordinates (cylindrical symmetry is assumed):
$$T \nabla^2 h(r) + q(r) = S \left( \frac{\partial h}{\partial t} - \frac{1}{\rho g} \frac{\partial f}{\partial t} \right) \qquad (B2)$$

Where f is the tidal forcing in the aquifer ($\gamma P_{atm}$ for atmospheric pressure or $2\gamma G\epsilon_T$ for earth tide) and $q(r)$ is the leakage flow (m/s) from the aquitard. Boundary conditions are:

- Away from the well ($r \to +\infty$), $\frac{\partial h}{\partial r} = 0$, such that $h(r = +\infty) = h_\infty$

- At the well, ($r = r_w$), $h(r_w) = h_w + \frac{p_{op}}{\rho g}$, with $h_w$ the measured water level in the well.

Solving Equation B1 yields, in complex notations

$$q = \frac{K' \partial p}{\rho g \partial z_{|z=0}} = \frac{K'}{b'} \left( \frac{p_{top} - f'}{\rho g} \frac{(1+i)\sqrt{Q}}{\sinh((1+i)\sqrt{Q})} - \left(h - \frac{f'}{\rho g}\right) \frac{(1+i)\sqrt{Q}}{\tanh((1+i)\sqrt{Q})} \right) \qquad (B3)$$

Injecting Equation B3 in Equation B2 and solving it yields:

$$\frac{\rho g h_w}{-p_a} = \frac{1}{\xi} \cdot \frac{1 - \gamma + (1 - \gamma') \frac{S'}{2S\sqrt{Q}} \left( \frac{(1-i)}{\tanh((1+i)\sqrt{Q})} - \frac{(1-i)}{\sinh((1+i)\sqrt{Q})} \right)}{1 + \frac{S'}{2S\sqrt{Q}} \tanh((1+i)\sqrt{Q})}$$

$$\frac{\rho g h_w}{\epsilon} = \frac{1}{\xi} \cdot \frac{2\gamma G + 2\gamma' G' \frac{S'}{2S\sqrt{Q}} \left( \frac{(1-i)}{\tanh((1+i)\sqrt{Q})} - \frac{(1-i)}{\sinh((1+i)\sqrt{Q})} \right)}{1 + \frac{S'}{2S\sqrt{Q}} \tanh((1+i)\sqrt{Q})}$$

With

$$Q = \frac{b^2 \omega}{2D'}; W = \frac{r_w^2 \omega}{2T}; c_r = \left(\frac{r_c}{r_w}\right)^2$$

$$\xi = 1 + i \cdot c_r \frac{W K_0(\beta r_w)}{\beta r_w K_1(\beta r_w)}$$

$$\beta r_w = \sqrt{W} \left( \frac{(1+i) \tanh((1+i)\sqrt{Q})}{\sqrt{Q}} S' + 2iS \right)^{\frac{1}{2}}$$

## B2. Model Inversion

Because $G' \ll G$ (Section 3.4),

$$\phi_{MT} - \phi_{PA} = \arg\left(\frac{\rho g h_w}{\epsilon}\right) - \arg\left(\frac{\rho g h_w}{-p_a}\right) = -\arg\left(\frac{1-\gamma}{1-\gamma'} + \frac{S'}{2S\sqrt{Q}}\left(\frac{(1-i)}{\tanh((1+i)\sqrt{Q})} - \frac{(1-i)}{\sinh((1+i)\sqrt{Q})}\right)\right)$$

Thanks to the measurement of $\phi_{MT} - \phi_{PA}$, and considering the calibration of S, S' & Q values thanks to 2022 pumping test, we recover the value of the ratio $\frac{1-\gamma}{1-\gamma'}$ which we assume constant. It enables to recover the total evolution of Q over the whole period. We then verify that $\tanh((1+i)\sqrt{Q}) \approx 1 \& \sinh((1+i)\sqrt{Q}) \gg 1$ in our ranges of values for Q, such that, as presented in Section 3.4:

$$\phi_{MT} - \phi_{PA} = \arg\left(1 + \frac{1 - \gamma'}{1 - \gamma} \frac{S'}{2S\sqrt{Q}}(1 - i)\right)$$
## Data Availability Statement

Input data sets, analysis code, intermediary results and pumping tests data useful for this research are available in Thomas et al. (2024). The finite difference code used for computing dynamic stress change is available in Aochi and Madariaga (2003b). Piezometric data comes from the BRGM piezometric monitoring network of Martinique, set up since 2005 as part of the European Water Framework Directive and is funded by the French Ministry of Environment and BRGM. The borehole national reference is BSS002NNZL (1175ZZ0154/NF4). The barometric data have been recovered from the French meteorological agency (Météo-France) soil station network. Barometric data were recovered between 1 January 2007 and 31 December 2022, from meteorological station Caravelle, which number is 97230001. These data are available in the French national platform (Météo France, 2024). Detailed instructions on how to recover barometric data on this platform, as well as station metadata are given in Thomas et al. (2024). Theoretical tidal strain was computed thanks to the SPOTL program (Agnew, 2012).

## Acknowledgments

This research project, HydroSeisme (2021–2023), was funded by ERC Persismo (Grant 865411), led by Harsha S. Bhat (ENS), by the initiative Ecosphere Continentale et Cotière (EC2CO) of the Institut National des Sciences de l’Univers/ Centre National de la Recherche Scientifique (INSU/CNRS), and from BRGM Sism’Eau Project (An internal M4 research program). We wish to thank Emile Gros from BRGM Martinique who participated to the pumping test campaign.

## References

Agnew, D. C. (2012). SPOTL: Some programs for ocean-tide loading [Computer software]. Scripps Institution of Oceanography. http://escholarship.org/uc/item/954322pg

Aki, K., & Richards, P. G. (2002). *Quantitative seismology* (2nd ed.). University Science Books. Retrieved from www.uscibooks.com

Aochi, H., Ducellier, A., Dupros, F., Delatre, M., Ulrich, T., de Martin, F., & Yoshimi, M. (2013). Finite difference simulations of seismic wave propagation for the 2007 Mw 6.6 Niigata-ken Chuitsu-Oki earthquake: Validity of models and reliable input ground motion in the near-field. *Pure and Applied Geophysics*, *170*(1), 43–64. https://doi.org/10.1007/s00024-011-0429-5

Aochi, H., & Madariaga, R. (2003a). The 1999 Izmit, Turkey, earthquake: Nonplanar fault structure, dynamic rupture process, and strong ground motion. *Bulletin of the Seismological Society of America*, *93*(3), 1249–1266. https://doi.org/10.1785/0120020167

Aochi, H., & Madariaga, R. (2003b). The 1999 Izmit, Turkey, earthquake: Nonplanar fault structure, dynamic rupture process, and strong ground motion [Software]. *Bulletin of the Seismological Society of America*, *93*(3), 1249–1266. https://doi.org/10.5281/zenodo.10225172

Arnaud, L., Coppo, N., Mathieu, F., & Lhotelin, M. (2011). *Prospections hydrogéologiques sur le site du Galion, commune de la Trinité (Martinique)*. Rapport final (RP-60419-FR). BRGM. Retrieved from http://ficheinfoterre.brgm.fr/document/RP-60419-FR

Barbosa, N. D., Caspari, E., Rubino, J. G., Greenwood, A., Baron, L., & Holliger, K. (2019). Estimation of fracture compliance from attenuation and velocity analysis of full-waveform sonic log data. *Journal of Geophysical Research: Solid Earth*, *124*(3), 2738–2761. https://doi.org/10.1029/2018JB016507

Barbosa, N. D., Hunziker, J., Lissa, S., Saenger, E. H., & Lupi, M. (2019). Fracture unclogging: A numerical study of seismically induced viscous shear stresses in fluid-saturated fractured rocks. *Journal of Geophysical Research*, 72. https://doi.org/10.1029/2019JB017984

Bastias Espejo, J. M., Rau, G. C., & Blum, P. (2022). Groundwater responses to Earth tides: Evaluation of analytic solutions using numerical simulation. *Journal of Geophysical Research: Solid Earth*, 127(10). https://doi.org/10.1029/2022JB024771

Boeut, S., Fujii, Y., Kodama, J.-I., Fukuda, D., Dassanayake, A., & Alam, B. A. K. M. (2020). Laboratory investigation on the permeability variation of fractured Inada Granite by multiple transient axial stress disturbances. *Pure and Applied Geophysics*, *177*(11), 5385–5396. https://doi.org/10.1007/s00024-020-02565-2

Bredehoeft, J. D. (1967). Response of well-aquifer systems to Earth tides. *Journal of Geophysical Research*, 72(12), 3075–3087. https://doi.org/10.1029/JZ072i012p03075

Brenguier, F., Campillo, M., Hadzioioannou, C., Shapiro, N. M., Nadeau, R. M., & Larose, E. (2008). Postseismic relaxation along the San Andreas Fault at Parkfield from continuous seismological observations. *Science*, *321*(5895), 1478–1481. https://doi.org/10.1126/science.1160943

Brodsky, E. E., & Prejean, S. G. (2005). New constraints on mechanisms of remotely triggered seismicity at Long Valley Caldera: Triggering constraints. *Journal of Geophysical Research*, 110(B4). https://doi.org/10.1029/2004JB03211

Brodsky, E. E., Roeloffs, E., Woodcock, D., Gall, I., & Manga, M. (2003). A mechanism for sustained groundwater pressure changes induced by distant earthquakes. *Journal of Geophysical Research*, 108(B8), 10. https://doi.org/10.1029/2002JB002321

Brown, S. R., & Scholz, C. H. (1986). Closure of rock joints. *Journal of Geophysical Research*, 91(B5), 4939–4948. https://doi.org/10.1029/JB091iB05p04939

Candela, T. (2014). Laboratory evidence for particle mobilization as a mechanism for permeability enhancement via dynamic stressing. *Earth and Planetary Science Letters*, 13. https://doi.org/10.1016/j.epsl.2014.02.025

Carlson, F., & Randall, J. (2012). MLU: A windows application for the analysis of aquifer tests and the design of well fields in layered systems. *Groundwater*, 50(4), 504–510. https://doi.org/10.1111/j.1745-6584.2012.00945.x

David, E. C., & Zimmerman, R. W. (2011). Elastic moduli of solids containing spheroidal pores. *International Journal of Engineering Science*, 49(7), 544–560. https://doi.org/10.1016/j.ijengsci.2011.02.001

Doan, M.-L., Brodsky, E. E., Prioul, R., & Signer, C. (2008). *Tidal analysis of borehole pressure: A tutorial* (p. 60). University of California, Santa Cruz; Schlumberger-Doll Research.

Dorel, J. (1978). *Sismicité et structure de l'arc des Petites Antilles et du bassin atlantique* (Vol. VI). University of Paris.

Elkhoury, J. E., Brodsky, E. E., & Agnew, D. C. (2006). Seismic waves increase permeability. *Nature*, 441(7097), 1135–1138. https://doi.org/10.1038/nature04798

Elkhoury, J. E., Niemeijer, A., Brodsky, E. E., & Marone, C. (2011). Laboratory observations of permeability enhancement by fluid pressure oscillation of in situ fractured rock. *Journal of Geophysical Research*, 116(B2), B02311. https://doi.org/10.1029/2011JB007759

Faoro, I., Elsworth, D., & Marone, C. (2012). Permeability evolution during dynamic stressing of dual permeability media: Permeability evolution. *Journal of Geophysical Research*, 117(B1), n/a-n-a. https://doi.org/10.1029/2011JB008635

Fortin, J., Guéguen, Y., & Schubnel, A. (2007). Effects of pore collapse and grain crushing on ultrasonic velocities and $V_p/V_s$. *Journal of Geophysical Research*, 112(B8), B08207. https://doi.org/10.1029/2005JB004005

Fortin, J., Stanchits, S., Vinciguerra, S., & Guéguen, Y. (2011). Influence of thermal and mechanical cracks on permeability and elastic wave velocities in a basalt from Mt. Etna volcano subjected to elevated pressure. *Tectonophysics*, 503(1–2), 60–74. https://doi.org/10.1016/j.tecto.2010.09.028
Gallagher, A., Fortin, J., & Borgomano, J. (2022). Seismic dispersion and attenuation in fractured fluid-saturated porous rocks: An experimental study with an analytic and computational comparison. *Rock Mechanics and Rock Engineering*, *55*(7), 4423-4440. https://doi.org/10.1007/s00603-022-02875-y

Graves, R. W. (1996). Simulating seismic wave propagation in 3D elastic media using staggered-grid finite differences. *Bulletin of the Seismological Society of America*, *86*(4), 1091-1106. https://doi.org/10.1785/BSSA0860041091

Gueguen, Y., & Dienes, J. (1989). Transport properties of rocks from statistics and percolation. *Mathematical Geology*, *21*, 1-13. https://doi.org/10.1007/bf00897237

Hsieh, P. A., Bredehoeft, J. D., & Farr, J. M. (1987). Determination of aquifer transmissivity from Earth tide analysis. *Water Resources Research*, *23*(10), 10-1832. https://doi.org/10.1029/WR023i010p01824

Hsieh, P. A., Bredehoeft, J. D., & Rojstaczer, S. A. (1988). Response of well aquifer systems to earth tides: Problem revisited. *5*. https://doi.org/10.1029/WR024i003p00468

Illien, L., Sens-Schönfelder, C., Andermann, C., Marc, O., Cook, K. L., Adhikari, L. B., & Hovius, N. (2022). Seismic velocity recovery in the subsurface: Transient damage and groundwater drainage following the 2015 Gorkha earthquake, Nepal. *Journal of Geophysical Research: Solid Earth*, *127*(2). https://doi.org/10.1029/2021JB023402

Illien, L., Sens-Schönfelder, C., & Ke, K.-Y. (2023). Resolving minute temporal seismic velocity changes induced by earthquake damage: The more stations, the merrier? *Geophysical Journal International*, *234*(1), 124-135. https://doi.org/10.1093/gji/ggad038

Kachanov, M. (1993). In J. W. Hutchinson & T. Y. Wu, (Eds.), *Elastic solids with many cracks and related problems* (Vol. 30, pp. 259-445). Elsevier. https://doi.org/10.1016/S0065-2156(08)70176-5

Lacquement, F., Baptiste, J., Bellier, V., Bernachot, I., Bernard, J., Bialkowski, A., et al. (2022). Carte des formations du régolithe de la Martinique comme appui pour la cartographie de l'aléa glissement de terrain. Rapport final V1. (BRGM/RP-71574-FR) (p. 8).

Maas, C. (1987). Groundwater flow to a well in a layered porous medium: 2. Nonsteady multiple-aquifer flow. *Water Resources Research*, *23*(8), 1683-1688. https://doi.org/10.1029/WR023i008p01683

Manga, M., Beresnev, I., Brodsky, E. E., Elkhoury, J. E., Elsworth, D., Ingebritsen, S. E., et al. (2012). Changes in permeability caused by transient stresses: Field observations, experiments, and mechanisms: Dynamic permeability. *Reviews of Geophysics*, *50*(2). https://doi.org/10.1029/2011RG000382

McMillan, T. C., Rau, G. C., Timms, W. A., & Andersen, M. S. (2019). Utilizing the impact of Earth and atmospheric tides on groundwater systems: A review reveals the future potential. *Reviews of Geophysics*, *57*(2), 281-315. https://doi.org/10.1029/2018RG000630

Météo France. (2024). Données climatologiques de base—Horaires. [dataset]. https://www.data.gouv.fr/fr/datasets/donnees-climatologiques-de-base-horaires/

Narasimhan, T. N., Kanehiro, B. Y., & Witherspoon, P. A. (1984). Interpretation of Earth tide response of three deep, confined aquifers. *Journal of Geophysical Research*, *89*(B3), 1913-1924. https://doi.org/10.1029/JB089iB03p01913

O'Connell, R. J., & Budiansky, B. (1974). Seismic velocities in dry and saturated cracked solids. *Journal of Geophysical Research*, *79*(35), 5412-5426. https://doi.org/10.1029/JB079i035p05412

Odling, N. E., Perulero Serrano, R., Hussein, M. E. A., Riva, M., & Guadagnini, A. (2015). Detecting the vulnerability of groundwater in semi-confined aquifers using barometric response functions. *Journal of Hydrology*, *520*, 143-156. https://doi.org/10.1016/j.jhydrol.2014.11.016

Okada, Y. (1992). Internal deformation due to shear and tensile faults in a half-space. *Bulletin of the Seismological Society of America*, *82*(2), 1018-1040. https://doi.org/10.1785/BSSA0820021018

Park, C. B., Miller, R. D., & Xia, J. (1998). *Ground roll as a tool to image near-surface anomaly* (pp. 874-877). SEG Technical Program Expanded Abstracts. https://doi.org/10.1190/1.1820627

Richter, T., Sens-Schönfelder, C., Kind, R., & Asch, G. (2014). Comprehensive observation and modeling of earthquake and temperature-related seismic velocity changes in northern Chile with passive image interferometry: Richter et al. *Journal of Geophysical Research: Solid Earth*, *119*(6), 4747-4765. https://doi.org/10.1002/2013JB010695

Roeloffs, E. (1996). Poroelastic techniques in the study of earthquake-related hydrologic phenomena. In *Advances in geophysics* (Vol. 37, pp. 135-195). Elsevier. https://doi.org/10.1016/S0065-2687(08)60270-8

Rojstaczer, S., & Riley, F. S. (1990). Response of the water level in a well to Earth tides and atmospheric loading under unconfined conditions. *Water Resources Research*, *26*(8), 1803-1817. https://doi.org/10.1029/WR026i008p01803

Thomas, A., Fortin, J., Vittecoq, B., Aochi, H., Violette, S., Maury, J., et al. (2024). Hydro-mechanical characterization of a fractured aquifer using groundwater level tidal analysis: Effect of pore pressure and seismic dynamic shear stresses on permeability variations (Version 2) [dataset]. *Zenodo*. https://doi.org/10.5281/zenodo.12073253

Thomas, A., Fortin, J., Vittecoq, B., & Violette, S. (2023). Earthquakes and heavy rainfall influence on aquifer properties: A new coupled earth and barometric tidal response model in a confined bi-layer aquifer. *Water Resources Research*, *59*(4), e2022WR033367. https://doi.org/10.1029/2022WR033367

Valois, R., Rau, G. C., Vouillamoz, J.-M., & Derode, B. (2022). Estimating hydraulic properties of the shallow subsurface using the groundwater response to Earth and atmospheric tides: A comparison with pumping tests. *Water Resources Research*, *36*(5). https://doi.org/10.1029/2021WR031666

Van Der Kamp, G., & Gale, J. E. (1983). Theory of earth tide and barometric effects in porous formations with compressible grains. *Water Resources Research*, *19*(2), 538-544. https://doi.org/10.1029/WR019i002p00538

Viens, L., Denolle, M. A., Hirata, N., & Nakagawa, S. (2018). Complex near-surface rheology inferred from the response of greater Tokyo to strong ground motions. *Journal of Geophysical Research: Solid Earth*, *123*(7), 5710-5729. https://doi.org/10.1029/2018JB015697

Virieux, J., & Madariaga, R. (1982). Dynamic faulting studied by a finite difference method. *Bulletin of the Seismological Society of America*, *72*(2), 345-369. https://doi.org/10.1785/BSSA0720020345

Vittecoq, B., Fortin, J., Maury, J., & Violette, S. (2020). Earthquakes and extreme rainfall induce long term permeability enhancement of volcanic island hydrogeological systems. *Scientific Reports*, *10*(1), 20231. https://doi.org/10.1038/s41598-020-76954-x

Vittecoq, B., Reninger, P.-A., Bellier, V., Taïlamé, A.-L., Nascimento, L., Gros, E., & Violette, S. (2023). Montagne Pelée volcano (Martinique, in the French Lesser Antilles) hydrogeological system revealed by high-resolution helicopter-borne electromagnetic imagery. *Hydrogeology Journal*, *31*(5), 1331-1352. https://doi.org/10.1007/s10040-023-02642-5

Vittecoq, B., Reninger, P.-A., Lacquement, F., Martelet, G., & Violette, S. (2019). Hydrogeological conceptual model of andesitic watersheds revealed by high-resolution heliborne geophysics. *Hydrology and Earth System Sciences*, *23*(5), 2321-2338. https://doi.org/10.5194/hess-23-2321-2019

Vittecoq, B., Reninger, P. A., Violette, S., Martelet, G., Dewandel, B., & Audru, J. C. (2015). Heterogeneity of hydrodynamic properties and groundwater circulation of a coastal andesitic volcanic aquifer controlled by tectonic induced faults and rock fracturing - Martinique island (Lesser Antilles - FWI). *Journal of Hydrology*, *529*, 1041-1059. https://doi.org/10.1016/j.jhydrol.2015.09.022
Walsh, J. B. (1981). Effect of pore pressure and confining pressure on fracture permeability. *International Journal of Rock Mechanics and Mining Sciences*, 18(5), 429-435. https://doi.org/10.1016/0148-9062(81)90006-1

Wang, C.-Y., Doan, M.-L., Xue, L., & Barbour, A. J. (2018). Tidal response of groundwater in a leaky aquifer—Application to Oklahoma. *Water Resources Research*, 15. https://doi.org/10.1029/2018WR022793

Wang, C.-Y., & Manga, M. (2010). Hydrologic responses to earthquakes and a general metric. *Geofluids*, 10(1-2), 206-216. https://doi.org/10.1111/j.1468-8123.2009.00270.x

Wang, C.-Y., Wang, C.-H., & Kuo, C.-H. (2004). Temporal change in groundwater level following the 1999 (Mw = 7.5) Chi-Chi earthquake, Taiwan. *Geofluids*, 4(3), 210-220. https://doi.org/10.1111/j.1468-8123.2004.00082.x

Wang, H., Xin, J., Zheng, X., Li, M., Fang, Y., & Zheng, T. (2020). Clogging evolution in porous media under the coexistence of suspended particles and bacteria: Insights into the mechanisms and implications for groundwater recharge. *Journal of Hydrology*, 582, 124554. https://doi.org/10.1016/j.jhydrol.2020.124554

Wang, H. F. (2000). *Theory of linear poroelasticity with applications to geomechanics and hydrogeology*. Princeton University Press.

Wells, D. L., & Coppersmith, K. J. (1994). New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. *Bulletin of the Seismological Society of America*, 84(4), 974-1002. https://doi.org/10.1785/bssa0840040974

Westercamp, D., Pelletier, B., Thibault, P. M., & Traineau, H. (1990). *Carte géologique de la Martinique au 1/50000ème*. Editions BRGM. Retrieved from http://infoterre.brgm.fr/formulaire/telechargement-cartes-geologiques-drom

Worthington, M. H., & Lubbe, R. (2007). The scaling of fracture compliance. *Geological Society, London, Special Publications*, 270(1), 73-82. https://doi.org/10.1144/GSL.SP.2007.270.01.05

Xue, L., Li, H.-B., Brodsky, E. E., Xu, Z.-Q., Kano, Y., Wang, H., et al. (2013). Continuous permeability measurements record healing inside the Wenchuan earthquake fault zone. *Science*, 340(6140), 6-1559. https://doi.org/10.1126/science.1237237

Yang, T., Wang, S., Wang, P., & Zhang, Z. (2022). Hydraulic and mechanical coupling analysis of rough fracture network under normal stress and shear stress. *KSCE Journal of Civil Engineering*, 26(2), 650-660. https://doi.org/10.1007/s12205-021-0660-2