# Analysing 50 yr of the Lacq induced seismicity (Southwestern, France) highlights the role of fluid injection

Laeticia Jacquemond, J. Letort, F. Cotton, M. Causse, J. R. Grasso, G.

Senechal, J. B. Ammirati, B. Derode, F. Grimaud, H. Pauchet, et al.

► To cite this version:

Laeticia Jacquemond, J. Letort, F. Cotton, M. Causse, J. R. Grasso, et al. Analysing 50 yr of the Lacq induced seismicity (Southwestern, France) highlights the role of fluid injection. Geophysical Journal International, 2024, 238, pp.214-234. 10.1093/gji/ggae119 . insu-04833967

HAL Id: insu-04833967

https://insu.hal.science/insu-04833967v1

Submitted on 12 Dec 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
# Analysing 50 yr of the Lacq induced seismicity (Southwestern, France) highlights the role of fluid injection

L. Jacquemond,<fnref n="1" />,<fnref n="2" /> J. Letort,<fnref n="1" /> F. Cotton,<fnref n="3" /> M. Causse,<fnref n="4" /> J.R. Grasso,<fnref n="4" />,* G. Senechal,<fnref n="5" /> J.B. Ammirati,<fnref n="1" /> B. Derode,<fnref n="1" /> F. Grimaud,<fnref n="1" /> H. Pauchet,<fnref n="1" /> S. Benhamed<fnref n="1" /> and M. Sylvander<fnref n="1" />

Accepted 2024 March 20. Received 2024 January 19; in original form 2023 May 31

## SUMMARY

The Lacq area in southwest France has been associated with continuous moderate induced seismic activity since 1969. However, the mechanisms driving this induced seismicity are not fully understood: reservoir depletion has been proposed as the main factor, and more recently wastewater injection has been suggested to play a more important role. The interpretation of these mechanisms relies heavily on the quality of earthquake locations, which we prove to be weak due to a lack of local instrumentation for several years. In order to provide the most complete and reliable induced event catalogue for the studies of the Lacq induced seismicity mechanisms and seismic hazard, we made an exhaustive compilation, analysis and improvement of all available catalogues. We also provided new earthquake detections and relocations in a 3-D velocity model from past and present temporary deployments never used for studying the Lacq area. Important remaining location uncertainties lead us to also carefully sort the events according to their location confidence, defining 3 classes of events (unconstrained location, location constrained within 2–3 km and 1–2 km, respectively). This new harmonized catalogue and the identification of well-constrained events, covering 50 yr of induced seismicity, allow us to propose that wastewater injection is almost certainly the main mechanism driving the seismicity, with (i) most of the constrained events located within the reservoir boundaries and (ii) the released seismic energy variations following variations in injection operations at different scales. In particular, we have also highlighted a change in the injection–seismicity relationship around 2010–2013. From 2013, despite lower injection volumes, seismicity remained persistent and some clusters of earthquakes were detected predominantly in spring, summer and early autumn, except in winter periods. From 2016, we observed a strong temporal relationship between days with higher rate/volume injections (approximately above 400 m<fnref n="3" /> d⁻<fnref n="1" />) and both clustered events and higher magnitude earthquakes (greater than 2.4).

**Key words:** Europe; Earthquake hazards; Earthquake source observations; Induced seismicity.

## 1 INTRODUCTION

Understanding induced seismicity is a pivotal concern for various geo-resource applications, including mining, fluid withdrawal, gas and oil extraction, underground gas storage and geothermal energy (e.g. Keranen & Weingarten 2018; Foulger et al. 2018), as

well as for risk mitigation. The study of anthropogenic seismicity has become crucial for European countries that have recently witnessed a surge in deep geothermal energy and carbon sequestration projects in densely populated areas (Grigoli et al. 2017). The felt induced ground shaking by local populations has recently forced the suspension or significant reduction of certain activities, such

<footnote n="1">¹IRAP, UMR 5277, Observatoire Midi-Pyrénées, Université Paul-Sabatier, CNRS, 31062, Toulouse, IRAP, France. E-mail: jean.letort@irap.omp.eu</footnote>

<footnote n="2">²Université Côte d'Azur, CNRS, Observatoire de la Côte d'Azur, IRD, Laboratoire Geoazur, 06560, Valbonne, France</footnote>

<footnote n="3">³Helmholtz-Centre Potsdam, German Research Center for Geosciences GFZ, Telegrafenberg, 14473 Potsdam, Germany</footnote>

<footnote n="4">⁴ISTerre, Univ. Grenoble Alpes, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, ISTERRE, 38000, Grenoble, France</footnote>

<footnote n="5">⁵Université de Pau et des Pays de l'Adour, E2S UPPA, CNRS, TOTAL, LFCR, UMR 5150, 64000, Pau, France</footnote>

<footnote>*Retired</footnote>
as the Groningen gas field in the Netherlands (Van Thienen-Visser & Breunese 2015; Grigoli et al. 2017), the Strasbourg geothermal project in France (Schmittbuhl et al. 2021), and the project in Basel, Switzerland, which incurred economic losses amounting to $9 million (Giardini 2009). Among these three well-known recent cases in Europe, the maximum-recorded magnitude has remained moderate ($M_w$ 3.6). Substantial efforts are being made and need to be continued to better characterize these areas, establish the link between seismicity and ground activities, and assess the associated seismic risk.

Another major European industrial site, the Lacq area in southwestern France, has received relatively less attention recently, despite experiencing significant induced seismic activity, with a local magnitude reaching approximately 4 in 2016 (Aochi & Burnol 2018). Lacq field has been exploited since 1950, initially for oil extraction from a shallow reservoir (0.7 km deep) and subsequently, starting in 1957, for gas extraction from a deep reservoir (4–5 km deep, Fig. 1b). This deep reservoir is a 20-km-long anticlinal reservoir, which culminates at a depth of 3200 m and is constituted by lower Cretaceous and upper Jurassic limestones, sealed under duc tile Cretaceous marls (Fig. 1). Gas field production has resulted in a cumulative 64.5 MPa depletion within the gas reservoir rocks (Grasso et al. 2021). The initial gas pressure at a depth of 3.7 km was 66.1 MPa (Segall et al. 1994). The intensive extraction period concluded in 2013, when the reservoir pressure had decreased to approximately 1.5 MPa (CLSIC4000 2019). In parallel with oil and gas extraction operations, wastewater disposal activities began in the shallow reservoir in 1957 and continued just below the deep gas field from 1974 at 4.5 km deep (Wittlinger 1980; Grasso & Feignier 1990; Grasso & Wittlinger 1990; Maury et al. 1992; Segall et al. 1994; Bardainne 2005; Grasso et al. 2021). These wastewater injections persist to this day, with millions of cubic metres injected into the shallow reservoir and tens of millions into the deeper one (Grasso et al. 2021). Grasso et al. (2021) noted four main phases for the operations (Table 1).

Around Lacq site, earthquakes are detected regularly, even nowadays, with tens of magnitude above 1 yr⁻¹, and usually few events above magnitude 2–3 yr⁻¹. This Lacq seismicity is well separated from the seismicity occurring along the Pyrenean range (Fig. 1a). The faults located a few kilometres south of the Lacq reservoir (e.g. Lagor fault, Fig. 1b), part of the North Pyrenean Frontal Thrust, are considered inactives (Choukroune 1992). Deformation rate through the North Pyrenean fault is not well resolved (less than 1 mm yr⁻¹) from geodetic surveys (Rigo et al. 2015). No known historical natural earthquakes have been documented in this Lacq area (historians and medieval writings make no mention of such events) (Rothé 1970, 1977; Grasso & Wittlinger 1990). Two seismic events occurred in 1969 and 1970, and are considered to be the onset of the Lacq seismicity (Grasso et al. 2021). These occurred 14 yr after the beginning of the Lacq exploitation. This was followed by a relatively calm period until 1976. Starting from 1976, moderate-magnitude events ($M_L$ 2.5+) were observed almost annually, peaking at eight events per year in the 1990s (Grasso et al. 2021), with six events registering a local magnitude above 4 (in the years 1972, 1978, 1978, 1981, 2013 and 2016). Several authors have proved that this seismicity is related to the Lacq operations (Wittlinger 1980; Grasso & Feignier 1990; Grasso & Wittlinger 1990; Maury et al. 1992; Segall et al. 1994; Bardainne 2005; Grasso et al. 2021). The time lag (14 yr) between the first operations and the beginning of the seismicity is not surprising in regard to other examples of induced seismicity in intraplate context where strain rates and background seismic activity are low (e.g. wastewater injection in Oklahoma:

induced seismicity started after a 13 yr lag, Zhai et al. 2019; extraction Groningen gas field: seismicity started 25 yr after, Smith et al. 2019).

## 1.1. State of the art and mechanisms driving the Lacq induced seismicity

Until recently, the Lacq seismicity has been largely interpreted as a consequence of the important depletion process from reduction of pore pressure from gas extraction. A possible influence of wastewater injections was considered to be negligible (Rothé 1977; Grasso & Wittlinger 1990; Segall et al. 1994; Bardainne et al. 2008). To support this hypothesis, a possible migration for seismicity from the centre towards the edge of the reservoir was suggested around the Lacq field, observed from 1974 until 1996 (from the analysis of a local network), as expected in the case of a depletion process (e.g. Bardainne et al. 2008). Moreover, the location of earthquakes on faults above the reservoir, and not inside, were also considered as in favour of this depletion more than a fluid related process. More recently, however, a more important role of injections has been proposed by Grasso et al. (2021). Taking into account the more recent seismicity (1996–2016), they did not confirm any lateral migration, and noted that 75 per cent of the 2013–2016 events cluster within 4–8 km depths, that is, close to and below the 4.5 km deep injection well. Grasso et al. (2021) therefore re-analysed entirely the Lacq seismicity and found, among other things, that the cumulative injected volume was high enough to trigger the 1969 $M_w$ 3 event at the onset of the Lacq seismicity, that the 1974 switch of injection depth from shallow to deep reservoir occurred before a sharp increase in seismicity and that an increase in injection rate in the year 2010 was followed by the two largest events that did occur in 2013 and 2016. They proposed that the interrelation between injection and extraction could be the most probable cause of the Lacq seismicity from 1969 until the present day.

## 1.2. Data availability

These previous studies and their interpretations depend strongly on the earthquake localization quality. We will show in this study that the uncertainties are important among the different existing catalogues that have pushed the interpretations toward depletion and/or injection processes. For instance, from the analysis of an exact same network (the IPGS-LGIT network, Fig. 2), different catalogues co-exist with different event locations: according to Bardainne et al. (2008), 70 per cent of the seismicity is above the deep gas reservoir where Guyoton et al. (1992) and Volant (1993) report a catalogue with events inside the deep reservoir or at its close boundary. Grasso et al. (2021) rely partly on the industry catalogue (CLSIC4000 2019), but this catalogue is built from only one instrumented well (Fig. 2) and focuses only on the seismicity close to the injection wells. The location quality of the catalogue is reduced, especially for events at distance from the borehole seismometers. Finally, for a long time period (2000–2020), no dense enough networks were deployed in the vicinity of the Lacq site (Fig. 2), and the national/regional catalogues (OMP 2021; RENASS 2021; CEA 2021; BCSF-Rénass 2022) are unable to constrain the earthquake locations for this area more precisely than approximately 10 km. To advance on understanding the Lacq seismic properties (earthquake locations, inside or outside the reservoir, presence of clusters, existence of a migration pattern), and the mechanisms driving this seismicity, an homogeneous database of locations is needed. It should
<figure><img src="image_2.png" /><figcaption>**Figure 1.** (a) Recent seismicity in southern France and northern Spain from 2019 until 2023 (French national RENASS catalogue). Circles indicate quarry blasts (dark red), natural seismicity (blue) and induced or triggered events (green). The Lacq seismicity (delimited by the black rectangle) is well separated from the seismicity related to the Pyrenean range. (b) Schematic cross section of the Lacq anticline structure, modified from Grasso *et al.* (2021), showing the position of the shallow oil reservoir (dark green zone) and the deep gas reservoir (red zone). The thick vertical (purple) bars are the injection well locations (injection occurs just below the gas reservoir at 4.5 km). Black dots are well located events from this study (1974–2021 Class-2 events—see main text, Section 2).</figcaption></figure>

<table><caption>**Table 1.** Lacq fluid manipulations (1955–present), from Grasso *et al.* 2021.</caption><thead><tr><th>Fluid manipulation phase</th><th>Extraction operation (depth, type)</th><th>Injection operation (depth, type)</th></tr></thead><tbody><tr><td>Phase I (1955–1974)</td><td>3–5 km, Gas reservoir</td><td>0.7 km, Wastewater</td></tr><tr><td>Phase II (1974–2006)</td><td>3–5 km, Gas reservoir</td><td>4.5 km, Wastewater</td></tr><tr><td>Phase III (2006–2012)</td><td>3–5 km, Gas reservoir</td><td>0.7 and 4.5 km, Wastewater</td></tr><tr><td>Phase IV (2013–present)</td><td>Negligible</td><td>4.5 km, Wastewater</td></tr></tbody></table>

be as robust as possible. Especially well-constrained (<1 km) depth estimates on the order of the deep reservoir width (500 m, Guyoton *et al.* 1992; CLSIC4000 2019) and their associated uncertainties are crucial.

To build this catalogue, we first compiled and compared all existing catalogues. We also relied on nearby temporary networks that have not been used yet in the Lacq context, especially the two dense networks from the MAUPASSAC (2017) (Lehujeur *et al.* 2021) and PYROPE (2013) (Chevrot *et al.* 2017) projects, deployed nearby for seismic imagery objectives. They provide additional detected events and constrained locations. We also deployed a temporary

network, in 2021 (Letort *et al.* 2023) allowing a precise view of the recent seismicity location 10 yr after the end of the extraction period. We relocated some events and analysed uncertainties. All new provided locations are made on a 3-D velocity model (Section 1). Despite our efforts, most of the events remain poorly constrained because of a lack of local stations. We then define different classes of events, according to their uncertainties (Section 2). Finally, for the 2016–2020 period, we had access to the daily wastewater injected volume allowing a detailed analysis of the link between seismicity and injection processes (Section 3).

## 2 RELOCATION USING A 3-D VELOCITY MODEL

The use of a 3-D velocity model during the location process for the Lacq area is necessary to provide improved earthquake locations (Volant 1993; Bardainne 2005; Bardainne *et al.* 2008). Without the use of a 3-D velocity model, depth estimation is particularly challenging (Volant 1993). Throughout this study, we rely on the 3-D velocity model for *P* waves originally developed by Guyoton
<figure><img src="image_3.png" /></figure>

**Figure 2.** (a) Networks and stations deployed nearby Lacq, coloured according to time. (b) Station locations. Yellow: Geopetrol sensors (3-components), all located in the same observation well LA405. Violet: stations (1-components) deployed by LGIT and IPGS, used by the Bardainne *et al.* (2008) and the EPOS (2018) catalogues (see main text). Orange: SISLACQ-UPPA stations (2002–2006). Black: SISLACQ2021 stations (3-components), still operating. Blue: stations of the PYROPE array (3-components), located at the western part of the reservoir. Green: 27 stations (3-components) from the MAUPASSAC experiment deployment were located at the south and west, less than 10 km of the Lacq reservoir, and included in the Lacq 3-D velocity model (dashed black box). Magenta: the only OMP-RENASS permanent station URDF (3-components) located in the Lacq area. Other permanent stations (PYPC, ATE, ORDF and REYF) are further than 20 km.

*et al.* (1992) from active seismic studies and modified by Bardainne (2005). This velocity model (called 3-D model, Figs 3a and b) covers most of the paths of the events/stations used in this study (for 3549/4114 events, all paths are included in the 3-D model). This velocity model is assumed to deviate less than 10 per cent from the true *P* velocities (Bardainne 2005).

However, 85 events detected in 2013 by the PYROPE experi-
ment (Fig. 3a) and 89 from the regional/national catalogues that we
aimed to relocate in the 3-D model during this study (see Section
2) are detected from stations outside the 3-D model (Fig. 3a), to
the south and to the west. To allow the use of stations deployed
further off the borders of the Lacq reservoir, we therefore extended
this model. This new 3-D velocity model (3-D Extended, Fig. 3a)
is estimated by first relying on the Bardainne (2005) extrapolation
toward the east and the south (Model 1, Fig. 3a), keeping the same
constant *P*-velocity value of the border of the original grid. Then, as
the geologic/tectonic structures are mainly oriented east–west, we
used a lateral extrapolation by stretching the velocities toward the
west to extend the Bardainne (2005) model and cover the PYROPE
stations. For the stations of the permanent network, further south,
with latitude lower than 43.2, we rely on the 1-D velocity model
used routinely from the regional/national localisation procedures
(Pauchet *et al.* 1999) and an interpolation between this 1-D model
in the south and ‘Model 1’ allows us to finally build the model ‘3-D
extended’. This extension of the velocity model was made to im-
prove the location of 174 events and only 17 of them, located from
the PYROPE experiment (red triangles, Fig 3a) will be considered
as well-located events (see Section 2).

This 3-D velocity model is valid for *P* waves traveltimes only.
We therefore linked the *S*-wave velocity to the *P*-wave velocity

by the use of a constant $V_p/V_s$ ratio. We defined this ratio from a
Wadati plot (Fig. 4), finding high ratios of around $1.9 \pm 0.1$ for local
stations in the Lacq area, and $1.8 \pm 0.1$ for stations further south
(MAUPASSAC data, Fig. 4). The large dispersion of observed S-P
delays is indicating a complex propagation model, with probably
local variations of $V_p/V_s$ ratios that we cannot take into account.
The observed high $V_p/V_s$ ratios are consistent with a fluid saturated
fault zone usually observed around injection wells (e.g. Tselentis
et al. 2011; Improta et al. 2017; Nur 1972).

3 BUILDING THE MOST COMPLETE
CATALOGUE FOR THE LACQ AREA
(1974–2021)

3.1. 1974–1997: the IPGP-LGIT networks

From 1974 to 1997, 5 yr after the first earthquake felt, the Univer-
sity of Strasbourg (IPGS) installed a network of seven, 1-component
stations and one, 3-component station (Wittlinger 1980). The net-
work was extended by the University of Grenoble (LGIT) in 1991
(Bardainne 2005). The arrival times from this IPGP-LGIT network
have been processed to build up two earthquake catalogues; (i) the
Volant catalogue, 1974–1995 (Volant 1993; Volant & Grasso 1994,
extended to 1995 by Boyer 1996; Lahaie & Grasso 1999) and (ii)
the Bardainne *et al.* (2008) catalogue, 1974–1997.

The arrival times from this IPGP-LGIT network have been pro-
cessed to build up two earthquake catalogues; (i) the Volant cata-
logues, 1974–1995 (Volant 1993; Volant & Grasso 1994), extended
to 1995 by Boyer (1996) & Lahaie & Grasso (1999), compiled and
provided on the EPOS platform (EPOS 2018). (ii) The Bardainne
<figure><img src="image_4.png" /><figcaption>**Figure 3.** (a) Map of networks and area covered by the 3-D velocity model. The central red box (called 3-D model) shows the limitations of the original 3-D *P*-wave velocity model derived from active seismic (Guyoton *et al.* 1992). The black box called Model 1 represents the limit of the extrapolation made by Bardainne (2005), and the large blue box (3-D Extended) shows the limits of the extrapolation of the 3-D velocity model made for this study. The MAUPASSAC experiment (green unfilled triangles) is located southwest of Lacq, only 27 stations (green filled triangles) are used for this study, all within the boundaries of 3-D model. The blue filled triangles indicate the PYROPE stations, the violet filled triangles the LGIT-EOST stations and in dark filled triangles the SISLACQ2021 stations. Pink squares are stations from the OMP permanent network. (b) *P*-velocities for the 3-D model derived from Guyoton *et al.* (1992) (3-D model), using a Lambert III projection, showing the anticlinal fold structure of the reservoir.</figcaption></figure>

*et al.* (2008) catalogue, 1974–1997. The Bardainne *et al.* (2008) catalogue is made from a compilation and a review of the IPGP-LGIT database. It is composed of more events (1759), covering a longer time period (1974–1997), including 594 common events with the EPOS (2018) catalogue (807 events, 1974–1995). For both catalogues, earthquakes have been localized in a similar 3-D earth model. However, we noted a systematic difference for the depth estimation of common events, with a quasi-systematic shift toward surface of +1.5 km + −2.7 km for Bardainne *et al.* (2008) compared to EPOS (2018) (Fig. 5). To investigate these differences, we relocated the Bardainne *et al.* (2008) catalogue, providing 1652 revised locations (some events located in the original catalogue are found not locatable by our new procedure). We used the NonLinLoc location algorithm (Lomax *et al.* 2000), introducing a $V_p/V_s$ ratio of 1.9, instead of 1.73 for Bardainne *et al.* (2008) (an example of a NonLin-Loc control file can be found Suppl. S1). From the Wadati diagram, we also identify incorrect *S* wave picks (Fig. 4b) (that we removed before relocation) when the *S* arrival time is above $t_P*0.95 + 0.2$ s, or below $t_P*0.85–0.2$ s, *tP* being the *P*-arrival time. *S*-wave picks on 1 component stations can be biased by converted phases of *P* coda. The 3-component stations of the subsequent networks (SISLACQ-UPPA and SISLACQ2021) clearly show that the *S*-wave arrival is often visible only on the horizontal components (Fig. S2). The relocated catalogue (Bardainne *et al.* 2008; Fig. 5) yields earlier origin times (−0.3 + −0.4 s) and deeper events compared to the original (−1.4 km + −2.7 km). The depths of the relocated catalogue are consistent with the EPOS (2018) depth estimations (+0.05 km + −1.4 km). The epicentral locations have small differences between these 3 catalogues (each time below 0.3 + −1.5 km), but once

again the relocated and the EPOS (2018) catalogues are closer (0.1 + −1 km). This coherency led us to choose to rely preferentially on the relocated and the EPOS (2018) catalogue for the construction of the final complete catalogue.

## 3.2. 1997–2020: lack of dense local networks

From the end of the 1974–1997 period, no long duration networks were deployed in the Lacq area. Only one permanent station from the RESIF national network has been available from 2012 (URDF, Fig. 2). Hence, for a time period of more than 20 yr, the Lacq area was not instrumented enough to allow high quality locations (Fig. 2). However, for short time periods, Lacq was covered by close or nearby temporary deployments.

### 3.2.1. SISLACQ, UPPA (2002–2006)

This temporary network of five stations (3-components) was deployed by the University of Pau (UPPA) between 2002 and 2006 (Figs 2 and 3). The bulletin, with phase arrival picks, compiled by Dubos-Sallée *et al.* (2006), was made available by UPPA. We relocated 166 events, using the NonLinLoc inversion procedure described previously (3-D velocity model, $V_p/V_s$ ratio of 1.9). In July–August 2005, one cluster was activated in the southern part of the reservoir, and recorded by this network (Figs 6a and c). This cluster is found shallow (<2–3 km) in the original catalogue (Dubos-Sallée *et al.* 2006), very close to the North Pyrenean front faulting system. A possible activation of this system would have an important impact on the seismic hazard. However, our relocation
<figure><img src="image_5.png" /><figcaption>**Figure 4.** (a) Wadati diagram for the four catalogues: Bardainne *et al.* (2008) (red), Maupassac (green), Pyrope (blue) and SISLACQ, 2021 (black). (b) Wadati diagram for the two catalogues Bardainne *et al.* (2008) (red) and SISLACQ, 2021 (Black). The dashed ellipse shows an area with very small S-P delays for the Bardainne *et al.* (2008) catalogue. These small delays are probably due to incorrect S wave identification using vertical 1-component sensors (see main text).</figcaption></figure>

procedure shows important location differences, especially for the longitude (depth differences 0.1 + -2.9 km, latitude differences -0.1 + -2.9 km, longitude differences -1 + -4.4 km, T0 differences 0.3 + -0.4 s). The events are found badly located from our NLL relocation, with an average vertical uncertainty of 2.5 km and an average horizontal uncertainty of 3.8 km (events are considered as Class-0 events, see next paragraph 2E). Some events of the July–August 2005 swarm become deeper, around 6 km deep (Fig. 6c). Therefore, these uncertainties don't allow us to confirm the existence of a shallow seismicity on the North Pyrenean front, nor if the seismicity, even possibly deeper, extends further south of the reservoir boundaries (similarly to the deep cluster, approximately -0.58, 43.7, well located during the IPGP-LGIT network period, Fig. 5a)

### 3.2.2. B2-PYROPE (2013)

During the PYROPE deployment (PYRenean Observational Portable Experiment), aiming to realize a seismic imaging of the lithospheric structure of the Pyrénées, some 3-component stations (velocimeters, Güralp CMG40-T) were installed nearby Lacq at its

western part (Fig. 6). Combined with the use of the Lacq permanent station, URDF, it allowed us to detect and locate 86 events in 2013. One part of these events were located in a deep seismic cluster (10–13 km) at the northwestern part of the reservoir. This cluster depth has to be interpreted precariously as an extrapolation of the 3-D velocity model had to be made to include the PYROPE stations. Depth estimation is usually strongly dependent on the velocity model accuracy (e.g. Laporte 2022). We can therefore assume that the low depth estimation uncertainties provided by NonLinLoc, less than 2 km, are underestimated.

### 3.2.3. The Geopetrol catalogue (2010–2017)

A seismic catalogue also exists for the Lacq area derived by the company in charge of the reservoir exploitation. Four 3-component borehole geophone arrays, all installed in the same observation well between 300 and 570 m deep (LA405, Fig. 6a), have been used to track the microseismicity near the injection wells during wastewater operations. The waveforms and arrival times are not freely available, and only an event catalogue covering the period 2010–2017 is provided (CLSIC4000 2012; CLSIC4000 2019). This catalogue focuses on events at a short distance from the injection wells
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Epicentral localisations (a) and S-to-N cross section (b) for common events located from the 3 following catalogues: Bardainne *et al.* (2008) (red), EPOS (2018) (dark green) and the relocation from this study (violet). Larger dots are the events which are considered well located in the EPOS (2018) catalogue (RMS < 0.08 s & Quality = 1). The black contour lines represent the contour of the deep gas reservoir. Events located above the reservoir from the Bardainne *et al.* (2008) catalogue, are located deeper, inside the reservoir or below, in the 2 others catalogues. The vertical violet thick line is the position of the injection well LA109. The black rectangle is the position of a cluster, assumed well located, south of the reservoir boundaries.</figcaption></figure>

(Fig. 6a), and is not exhaustive in space or time. The robustness of the reported localisations is also uncertain as locations are derived from sensors located in one single observational well (a large azimuthal gap) using a 1-D velocity model (CLSIC4000 2012) and a $V_p/V_s$ ratio of 1.82. Despite these important uncertainties and limitations, we added the Geopetrol catalogue to our final compilation for completeness, but we record and acknowledge the higher location uncertainties of these events.

using the same localisation procedure (3-D model, $V_p/V_s$ ratio of 1.9). In addition to the 10 short-period sensors and the permanent station URDF, 5 RaspberryShakes were deployed by OMP and UPPA during the year, and contribute to the arrival times. Almost all events located during this year are inside or below the reservoir (Fig. 7). The catalogue provides constrained localisations (uncertainties estimated through NonLinLoc are lower than 1–2 km), and thus allows a reliable comparison between the seismic properties of the intensive extraction period (<2000) and 10 yr after its end.

### 3.2.4. The MAUPASSAC catalogue (2017)

Finally, the MAUPASSAC very dense seismic experiment (Lehujeur *et al.* 2021), aimed to image the crust in the Mauleon Basin, south from Lacq, and provides an earthquake catalogue encompassing the Lacq area (Fig. 6). Twenty-seven MAUPASSAC stations close to the Lacq area (Fig. 3a) have been used, as well as the permanent station URDES, to relocate this catalogue of 117 events found to occur in the Lacq area. We use the 3-D lacq velocity model, with a constant $V_p/V_s$ ratio of 1.8, estimated from the Wadati diagram, because stations are located primarily to the south of the study area (see discussion of $V_p/V_s$ in section 1, Fig. 4). During this time period, one cluster was activated in the deep northwestern part of the reservoir, close to the one identified from the PYROPE network.

### 3.3. 2021: Deployment of a new network (SISLACQ2021, Letort *et al.* 2023)

From January 2021, a new deployment of 10 short period 3-components stations, provided by GFZ and deployed by the university of Toulouse (OMP), allowed us to detect and locate 295 events

### 3.4. Regional and national catalogue

We also added 614 events, which have been detected only by the regional (OMP 2021) and/or the national (RENASS 2021; CEA 2021; BCSF-Rénass 2022) catalogues, for the period 1969–2020. These locations are made on a 1-D model aiming to locate natural Pyrenean earthquakes, and which is therefore not well suited for Lacq events. These locations have larger uncertainties (>5 km). When enough phase pickings for close stations, including URDF, were available, we relocated the event using the 3-D velocity model. 214 events were thus relocated, but the 3-D velocity model is unadapted for stations in the Pyrénées range, and these new locations are not well constrained with average epicentral location uncertainty provided by NLL of 9 km.

Aiming to improve these unconstrained locations further, a difficult task in regard to the network configuration and the lack of close stations, a relative localization procedure was applied (Suppl. S3), fixing as master events those located during the SISLACQ2021 period, when we have reliable reference localizations. Consequently,
<figure><img src="image_7.png" /><figcaption>**Figure 6.** (a) Map of located earthquakes, included in the final catalogue from temporary networks: for the PYROPE (2013) experiment (blue), MAUPASSAC (2017) (green), Géopetrol 2010–2017 (yellow) and from the relocation of SISLACQ1 (2002–2006) (brown). The isolines are the contours of the reservoir. The size of the dots represent the class of events (big dots = class-2 = best). The red square represents the observation well LA405, and the black square the injection wells. (b) S-to-N cross section. The vertical purple line is the deep injection well. (c) Comparison of the Dubos-Sallée *et al.* (2006) locations (orange dots) and the relocation of SISLACQ1 (brown dots). The events are plotted on the schematic geological cross section extracted from Grasso *et al.* (2021). The main faults are represented by black lines. The green and red contours represent the oil and gas reservoir, respectively.</figcaption></figure>

134 events were relocated relatively to these master events using the HypoDD code (Waldhauser & Ellsworth 2000).

## 3.5. Final catalogue: compilation and homogenization

We compiled all previous described catalogues in a single one (4114 events from 1969 to 2021). When an event is located from different catalogues, we choose the location from the most trusted catalogue. The order of catalogue trustworthiness is defined empirically (1 = best catalogue to 15 = less well constrained catalogue), taking into account the network geometry, the presence of close stations, the use or not of the 3-D velocity model, etc. (Table 2). Hence, we defined the best catalogues as being the local catalogues, with at least 5 stations covering the Lacq area. Then, we use the temporary networks. Then, comes the 9- Geopetrol catalogue, with only one instrumented well. Finally, the less constrained catalogues are the regional/national catalogues.

### 3.5.1. Magnitude estimation

It is crucial to have a homogeneous magnitude estimation protocol for the combined catalogue in order to be able to study earthquake evolutions. However, the diversity of catalogue and data makes this task challenging. For older data, waveforms are no longer available (Bardainne *et al.* 2008), and we can only rely on the published catalogue magnitudes. We choose as a reference magnitude the local magnitude used at the French Atomic Agency (called MLCEA), as it is the reference magnitude chosen for the Sihex, compiled French national catalogue (Cara *et al.* 2015). This MLCEA local magnitude

is also the one used at the Toulouse University (Sylvander *et al.* 2021).

Then, all RENASS local magnitudes are converted to MLCEA local magnitude following Cara *et al.* (2015) (Suppl. S4). For the Geopetrol catalogue, no waveforms are available. We then rely on the local magnitude provided by CLSIC4000 (2019), which have been calibrated using the higher magnitudes from the national RENASS catalogue (CLSIC4000 2019). We then assumed that they are RENASS-type local magnitudes, and we converted them into MLCEA local magnitudes.

For the SISLACQ2021 and MAUPASSAC catalogues, built during this study from new detections/localisations, we directly computed MLCEA from the waveforms following Sylvander *et al.* (2021). For PYROPE events, no magnitude has been estimated, due to issues in the response of the temporary stations. Therefore, if a PYROPE event is included in another catalogue, its magnitude will be taken from the other catalogue; otherwise, no magnitude estimation is provided.

Finally, for the Bardainne *et al.* (2008), events have been calibrated from the RENASS catalogue (Bardainne *et al.* 2008). We then apply the MLRENASS/MLCEA conversion (Suppl. S3).

### 3.5.2. Uncertainty estimation—defining 3 classes of constrained events

The final combined catalogue still contains large uncertainties. In order to be able to correctly interpret the locations in terms of process, we put efforts on estimating well located events, defining 3 classes (Class-0: not well located, Class-1: events located at 2–3 km,
<figure><img src="image_8.png" /><figcaption>**Figure 7.** (a) Map of the SISLACQ 2021 stations (triangles) and located events during the year 2021 (dots). The white square delimits the location of a cluster (SWE—see later) activated in August 2021 at the northwestern part of the reservoir. (b) Cross section (south-north) of the event location compared to the position of the reservoir (red) and the well LA109 (purple tick line). The 2021 seismicity is plotted on the schematic geological cross section, extracted from Grasso *et al.* (2021). The white rectangle shows the position of the August 2021 cluster.</figcaption></figure>

Class-2: events located at 1–2 km). The 15 catalogues/relocation procedures do not constrain the event localisations equally, as can be seen by the observed important localisation differences between the 15 catalogues (Fig. S5). Locations are relatively coherent among catalogues before the year 2000, with differences usually less than 2–3 km (cases 2–4, Table 2). Then, from the year 2000, the comparison between the different catalogues highlights larger differences, especially for the regional/national catalogues (cases 13–15, Table 2, Fig. S5) for which the localisations can differ more than 5 km to those estimated from the temporary networks (Fig. S5).

A robust estimation of uncertainties is needed. The network geometries were not sufficient to provide low uncertainty localisations (Geopetrol, RENASS, OMP catalogues, cases 9, 10, 13, 14, 15, Table 1), and hence those catalogues were down-weighted in the ranking (classification) scheme. For other catalogues, uncertainty criteria such as RMS values, number of used phases, and azimuthal gap are no longer available (cases 4, 9). For a robust estimation of uncertainties, we will mostly rely on the 3-D relocations made for this study, and also rely on the proposed EPOS (2018) location uncertainty criteria (EPOS 2018).

The procedure to identify well located events is the following :

3.5.2.1 *Class-2 events (events well located at 1–2 km)* For all relocations relying on the NonLinLoc relocalisation process (cases 1, 3, 5, 6, 7, 11, 12, Table 2), we rely on the horizontal and vertical uncertainties provided by the location algorithm. We assume an event is well located (Class-2 events) if the following criteria are met:

(i) NLL horizontal uncertainty lower than 1.5 km.

(ii) NLL vertical uncertainty lower than 2 km.

(iii) At least 6 phases used.

(iv) RMS value lower than 0.2 s.

(v) Azimuthal gap below 200°.

We also add as Class-2 events those from the EPOS (2018) catalogue, for which a constrained events list is provided as well as RMS values. Two categories are present in the EPOS (2018) catalogue (category 1: constrained events—category 0: non-constrained events), defined following Guyoton *et al.* (1992): well-constrained category 1 events are events (1) recorded at least by 6 stations, (2) a rms lower than 0.2 s, (3) the location algorithm used (REL3D, Roecker 1982) must converge with less than 0.5 km change in location on the final iteration and (4) the condition number (i.e. the ratio of largest to smallest eigenvalues) of the matrix of partial derivatives must be less than 50.
<table><caption>**Table 2.** Catalogues used during this study, sorted according to their Quality Metrics.</caption><thead><tr><th>Catalogue rank order</th><th>Catalogue</th><th>Used velocity model</th><th>New locations from this study</th><th>Catalogue time range</th><th>Magnitude MLCEA estimation method</th><th>Number of stations above the reservoir / Nearby stations (less than 15 km)</th></tr></thead><tbody><tr><td>1</td><td>SISLACQ2021</td><td>3-D</td><td>yes</td><td>2021</td><td>On waveforms, following (Cara et al. 2015; Sylvander et al. 2021)</td><td>10/15</td></tr><tr><td>2</td><td>EPOS (2018)</td><td>3-D</td><td>no</td><td>1975–1995</td><td>No conversion made (original EPOS magnitude)</td><td>8/12</td></tr><tr><td>3</td><td>Relocation Bardainne et al. (2008)</td><td>3-D</td><td>yes</td><td>1974–1997</td><td>Converted from RENASS magnitude to CEA magnitude (Cara et al. 2015)</td><td>8/12</td></tr><tr><td>4</td><td>Original Bardainne et al. (2008)</td><td>3-D</td><td>no</td><td>1974–1997</td><td>Converted from RENASS magnitude to CEA magnitude (Cara et al. 2015)</td><td>8/12</td></tr><tr><td>5</td><td>MAUPASSAC (2017)</td><td>3-D</td><td>yes</td><td>2017</td><td>On waveforms, following (Cara et al., 2015; Sylvander et al. 2021)</td><td>1/27</td></tr><tr><td>6</td><td>PYROPE (2013)</td><td>3-D Extended</td><td>yes</td><td>2013</td><td>No Magnitude estimation</td><td>1/6</td></tr><tr><td>7</td><td>Relocation Dubos-Sallée et al. (2006)</td><td>3-D</td><td>yes</td><td>2000–2006</td><td>Converted from RENASS magnitude to CEA magnitude (Cara et al. 2015)</td><td>3/5</td></tr><tr><td>8</td><td>Dubos-Sallée et al. (2006)</td><td>3-D</td><td>no</td><td>2000–2006</td><td>Converted from RENASS magnitude to CEA magnitude (Cara et al. 2015)</td><td>3/5</td></tr><tr><td>9</td><td>Geopetrol catalogue</td><td>1-D</td><td>no</td><td>2010–2018</td><td>Converted from RENASS magnitude to CEA magnitude (Cara et al. 2015)</td><td>1/1 (4 in a same well)</td></tr><tr><td>10</td><td>Relocation hypoDD RENASS</td><td>1-D</td><td>yes</td><td>1980–2021</td><td>Converted from RENASS magnitude to CEA magnitude (Cara et al. 2015)</td><td>1/1</td></tr><tr><td>11</td><td>Relocation OMP</td><td>3-D Extended</td><td>yes</td><td>1998–2014</td><td>On waveforms, following (Cara et al. 2015; Sylvander et al. 2021)</td><td>1/1</td></tr><tr><td>12</td><td>Relocation RENASS</td><td>3-D Extended</td><td>yes</td><td>1998–2021</td><td>Converted from RENASS magnitude to CEA magnitude (Cara et al. 2015)</td><td>1/1</td></tr><tr><td>13</td><td>Original OMP</td><td>1-D</td><td>no</td><td>1997–2014</td><td>On waveforms, following (Cara et al. 2015; Sylvander et al., 2021)</td><td>1/1</td></tr><tr><td>14</td><td>BCFS-RENASS 2022 (Sihex extended)</td><td>1-D</td><td>no</td><td>1969–1996 and one event in 2012</td><td>(Cara et al., 2015)</td><td>1/1</td></tr><tr><td>15</td><td>Original RENASS</td><td>1-D</td><td>no</td><td>1980–2021</td><td>Converted from RENASS magnitude to CEA magnitude (Cara et al. 2015)</td><td>1/1</td></tr></tbody></table>
For Class-2 events selection, we were more restrictive here by selecting EPOS constrained events (category 1 events), but only with RMS values lower than 0.03 s.

364 events only (9 per cent) belong to the Class-2 constrained events, including 156 events in 2021 (SISLACQ2021), 83 from the EPOS (2018) catalogue, 108 from the Bardainne et al. (2008) relocation, and only one from MAUPASSAC (2017) and 17 from PYROPE (2013).

### 3.5.2.1 Class-1 events (events well located at 2–3 km)
As for the Class-2 events, we rely on the NonLinLoc horizontal and vertical uncertainties (cases 1, 3, 5, 6, 7, 11, 12, Table 1). We assume an event is satisfyingly located (Class-1 events) if the following criteria are met:

(i) NLL horizontal uncertainty lower than 2 km.

(ii) NLL vertical uncertainty lower than 3 km.

(iii) At least 6 phases used.

(iv) RMS value lower than 0.3 s.

(v) Azimuth gap below 270°.

Category 1 events from the EPOS (2018) catalogue, with RMS values between 0.03 and 0.05 s are kept as Class-1 events. 289 events belong to these Class-1 constrained events. Including the Class-2 events, it is then only 16 per cent of the whole final catalogue that can be considered as well located at 2–3 km. All other events are considered as Class-0 events (3461 events: 84 per cent).

Despite the lack of continuous dense instrumentation, temporary networks have been installed at various time periods. Hence well constrained events span a time period broad enough (Figs 8 and 9) to study the temporal evolution of the seismicity in regards to extraction and injection operations (Section 4). After removing the Class-0 events, almost all are located around the reservoir boundaries (Fig 8b). It also becomes clear that the seismicity is divided in several clusters/swarms (Fig. 8b). This coherency between Class-1 and Class-2 events locations on tight clusters are a good clue of the success of this events classification process. Remember that we could have underestimated the uncertainties for the deep northwestern cluster in 2013 (PYROPE), which falls in the Class-2 category but could be better ranked in the Class-1. The second deep northwestern cluster in 2017 reached the Class-1 category but with the maximum allowed uncertainties (3 km vertical—2 km horizontal). Therefore, the exact position of these two clusters should be interpreted with care. Their epicentral location at the northwestern part of the reservoir seems reliable at ±2 km. They could be outside or just inside the lateral extent of the reservoir boundaries. They seem to be related to deep faults far below the reservoir, at 10–12 km, but we cannot discard the hypothesis that they are closer to the reservoir boundaries, at 7–8 km.

## 4 IDENTIFICATION AND INTERPRETATION OF CLUSTERS

The improved earthquake locations highlight the presence of clusters and swarms (Fig. 8) as already noted by Wittligner (1980), Grasso & Wittlinger (1990), Volant & Grasso (1994) and Bardainne et al. (2008). Swarms/clusters are observed in several tectonic contexts, such as tectonic plate boundaries (e.g. Nadeau & McEvilly 1999), glaciers (Helmstetter et al. 2015) and volcanoes (Thelen et al. 2011). They are usually explained by aseismic processes or/and fluid/gas diffusion, the two mechanisms not excluding one another (Vidale & Shearer 2006; Roland & McGuire 2009; Chen & Shearer 2011; Chen et al. 2012; De Barros et al. 2020). Induced seismicity generally has a higher rate of repeating events and a larger proportion of small clusters compared to natural events (Zaliapin & Ben-Zion 2016). In the context of the Lacq gas field, with deep injected wastewater, it is therefore important to better characterize these clusters and to study the link between them and fluid/aseismic processes.

Clusters can be defined based on their waveform similarity (e.g. Baisch et al. 2008; Hatch et al. 2020; Gao et al. 2021, Derode et al. 2023), a well-suited tool here, as the event locations are not reliable for a long time period. Moreover, identifying similar events based on their waveform similarity also provides an estimate of location uncertainties—according to the location differences inside a family of similar-waveform events (SWEs).

For the only permanent station in the Lacq area (URDF), all events between 2013 and 2021 are selected for cross-correlation study. The associated waveforms are filtered between 8 and 20 Hz; the waveforms are cut 1 s before the P-wave arrival, 10 s after—therefore encompassing the S-wave arrivals; and cross-correlated in the time domain following Derode et al. (2023). If waveforms from at least four events have their cross-correlation coefficients above 0.6, these events define a family with similar waveforms. Forty-nine families are found (Fig. 10a), encompassing 24 per cent of the catalogue.

We can note that the locations for the ‘June 2020’ SWE family (in red, Fig. 10a) are poorly constrained. Even if the waveforms are very similar (Fig. S2), with high correlation values and pointing clearly to close events, the locations can vary by several kilometres (std = 7 km). There are no well located Class-1 and Class-2 events in this family. For the other identified families, most of the SWEs are relatively well clusterized for a same family, usually spreading over a few kilometres (<2–3 km), and it became especially true when looking for the well located Class-1 and Class-2 events only (Fig. 10a). This coherency between similar waveform families and Class-1 and Class-2 event clusters is another clue of the accuracy of our final catalogue uncertainty estimation.

Interestingly, several similar-waveform events (SWEs) are found on the edges of the reservoir, especially at its deep northwestern part, where no permanent background seismicity (non-SWEs) is present. SWEs are also found around the injection wells (Fig. 10a). It is also surprising to note that most of these SWEs occur preferentially during the summer months (Fig. 10b), while the non-SWEs occur over all periods of the year (Fig. 11). We will discuss these SWEs and their annual cycle in the next section.

## 5 STUDYING THE LINK BETWEEN INJECTION AND SEISMICITY

### 5.1. Stable observed features of the seismicity over time

Despite higher uncertainties (84 per cent of the catalogue is composed of Class-0 events), our new catalogue allows us to affirm: (1) Seismicity has not decreased after 2013 (the end of the intensive extraction period), with earthquakes regularly reaching $M_L$ 3–$M_L$ 4 magnitudes until present (Fig. 9). We note a decrease of released seismic energy and of events with $M_L$ > 2.4 until around 2000, 13 yr before the end of the extraction period. Later, the seismicity even seems to increase again (Fig. 9, Fig. 12). (2) There is no obvious outward migration of seismicity over time (Fig. 13a). Thanks to SISLACQ2021, we even clearly highlight the same seismicity areas
<figure><img src="image_9.png" /><figcaption>**Figure 8.** (a) All events in the final catalogue. The colour scales the origin catalogue (see Table 2). Catalogues 15–10 ( derivative of the pink colour) are locations derived from the regional catalogues OMP & RENASS. Others catalogues are 1-SISLACQ2021 (black), 2-EPOS2018 (dark green), 3-Bardaine *et al.* (2008) relocation (violet), 4-Bardaine *et al.* (2008) original (red), 5-MAUPASSAC (green), 6-Pyrope (blue), 7-Dubos-Sallée *et al.* (2006) relocation (braun), 8- Dubos-Sallée *et al.* (2006) original (Orange) and 9- Géopetrol (Yellow). The black square represents the limit of the box represented in (b) where all Class-1 and Class-2 events are constrained. (b) Location of Class-1 (small dots) and Class-2 events (big dots) events. Panels (c) and (d) are SN cross sections of, respectively (a) and (b), depth in metres. The red line indicates the position of the well LA109.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** (a) All estimated magnitude (MLCEA) reported in this study. The colour scales with the catalogue. The black line is the $M_L$ 2.4 threshold. (b) Evolution of the number of detected events with a magnitude above $M_L$ 2.4, according to the date (for each day $X$, we count the number of events with magnitude above $M_L$ 2.4 detected between $X$ – 6 months and $X$ + 6 months, to smooth the variations).</figcaption></figure>

in 2021 as for the years before 2000, except for SWEs at the reservoir periphery. (3) Seismicity is concentrated in the few kilometres below the top of the reservoir, well defined by active seismic surveys and boreholes (Guyoton *et al.* 1992). The seismicity extends

mainly between 0 and 1.5 km from the top reservoir. There are very few constrained earthquakes above the reservoir (Fig. 13b). It is worth noting that the reservoir has a constant thickness of 450–500 m (Guyoton *et al.* 1992; CLSIC4000 2019). Therefore, most earth-
<figure><img src="image_11.png" /><figcaption>**Figure 10.** (a) Epicentral location of the 49 similar-waveform events (SWEs) families detected by cross-correlation during the period 2013–2021 (Full dots, the colour scales the year of each family). We identify well-constrained events, Class-1 and 2 events (large full dots) and less-constrained events, Class-0 events (small full dots) which are part of the 49 identified families. Black dashed lines are the reservoir isolines. (b) Number of SWEs detected each month (red histogram) and non-SWEs (black histogram) according to the date, from 2013 to 2021. Grey rectangles delimit the 4 winter months (November–February).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** (a) Location of SWEs (filled red dots) and non-SWEs (black dots) occurring in November and December for the years 2013–2021. The purple black dashed lines are the reservoir boundaries. Big red dots are Class-2 events and medium-size red dots are Class-1 events. (b) Local magnitude ($M_L$) for the events of (a), for each year from 2013 to 2021. (c) and (d) Same as (a) and (b) for January–February. (e) and (f) Same as (a) and (b) for March–April. (g) and (h) May–June (i) and (j) July–August (k) and (l) September–October.</figcaption></figure>

quakes are located within or less than 1 km below the reservoir (except for some peripheral deeper events, between 6 and 12 km, which are mostly SWEs and/or clustered events). (4) The strongest earthquakes seem to be concentrated a few kilometres around the

injection wells (Fig. 13c). This stability of the active zones over time (Fig. 13a) probably reflects an identical process driving the seismicity. Since the extraction stopped in 2013, while injection has always been carried out, it thus appears that injection could
<figure><img src="image_13.png" /><figcaption>**Figure 12.** (a) Seismic energy released per year (black dots), estimated from local magnitude converted to moment magnitude *M*<sub>w</sub> (see Supplementary S7). The smoothed averaged variation of this seismic energy (black line) is computed on a 5-yr running window. Dashed lines represent the period of the four phases of operations, as described by Grasso *et al.* (2021). (b) Yearly total injected volume in the deep gas reservoir (black thick line) and injected volume inside the well LA102 solely (in red) and in the well LA109 solely (in blue). From 1974 until 2000, the injections were performed only in the LA109 well. From 2000, the two wells are used (c) Reservoir pressure evolution over time (black line) and LA109 (blue triangle) and LA102 (red triangle) measured pressures.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** (a) and (b) Location of Class-1 and Class-2 events, coloured to date (1974–2021). Bigger dots are Class-2 events. The squares are the location of the injection wells and the dashed black lines are the isolines of the reservoir top. (c) and (d) Location of events with local magnitude above 3. The red big dots are Class-2 events, the orange moderate size dots are Class-1 events, and the small blue dots are Class-0 events. 11 of the 15 well-located events with $M_L > 3$ are part of the seismic cluster related to the injection wells La102 and La109 (The 4 others seem related to the northwestern cluster, Arthez-de-Béarn, cluster: -0.655, 43.46).</figcaption></figure>
be a major cause for the seismicity, as proposed by Grasso *et al.* (2021).

## 5.2. Observed relationship between injection and released seismic energy from 1969 to 2021

Only annual injection values are provided from the two wells LA109 and LA102 from 1974 until 2021 (CLSIC4000 2021). Comparison of these annual injected volumes with the seismic energy released over the entire reservoir shows very similar trends until 2012–2013 (Fig. 12a). We can note an increase of seismic energy and injection volumes until 1980 and then a simultaneous decrease until the 2000 s. Then, a period of increase in seismic energy associated with an increase in injected volumes until 2007, finally again a decrease of both volumes and seismicity until 2011–2012. At first order, between 1969 and 2012 this correlation tends to prove a major control of injections on the high magnitude earthquakes. From 2012 to 2013, injections remain moderate (on average <400 m³ d⁻¹ over the year) but the released seismic energy increases (Fig. 12a; an increase driven by the two M_L 4 events in 2013 and 2016). The injection-seismicity relationship seems to surprisingly come to a halt, and this roughly coincides with the end of the gas extraction period in 2013. Around 2013, several notable changes occurred that could explain changes in the injection-seismicity relationship. The amount of re-injected reservoir water (which used to be important, sometimes more than half of the total volume injected) dropped drastically in 2013 and filtered-water began to be injected instead, to keep the injectivity of the wells (CLSIC4000 2021, slide 5). The end of the intense extraction period also marks the beginning of a slow reservoir pressure rise following a period of fast pressure decrease (CLSIC4000 2019, see Fig. 12c). The pressure at the LA102 injection-well even seems to have stabilized as early as 2010 with the extraction operation decreasing (Fig. 12c). Do these changes in the reservoir properties affect the seismicity-injection relationship? The well LA102 shows the exact same pressure evolution as the global reservoir pressure (Fig. 12c). On the contrary, LA109 has a 5 times higher pressure (Fig. 12c), and its variation is delayed compared to the reservoir pressure evolution. This could mark a weaker connectivity of LA109 to the reservoir when LA102 is better connected. The addition of the use of LA102 in 2000 could have affected the relationship between seismicity and injection.

## 5.3. Relationship between injection, higher magnitude events, and SWEs from 2016 to 2021

We also focused on the recent period (2016–2021), as we had access to daily data on the volume of wastewater injected during this period (Cretace4000 Committee 2023). This highly detailed information on injection rate and volume (now daily instead of annual) enables a more refined study of the connection between seismic activity and injection in recent years.

First, we examined the correlation between the injected volume and earthquakes of higher magnitude ($M_L > 2.4$, i.e. $M_w > 2$). We observed that the number of days of ‘intense’ injection during 3 months (>390 m³ d⁻¹) correlates with the number of $M_L > 2.4+$ events for the same period (Figs 14a and c). Tests were conducted with other thresholds for the minimum injection rate (Fig. S6). The threshold of 390 m³ d⁻¹ showed the strongest correlation, but thresholds of 370 or 410 m³ d⁻¹ also

appeared acceptable, encompassing the annual variation in seismicity. We also observed an anticorrelation between the number of days with low injection rates (<300 m³) and seismic activity (Fig. S7).

Next, we identified a systematic time delay of a few months between a period of intensive injection and the activation of SWEs (Figs 14b and c). We have intense injection days, usually from November to March, associated with $M_L > 2.4+$ events (Fig. 14a), followed by the activation of SWEs, mainly from April/May to October (Fig. 11, Fig. 14). It appears that the number of days with high-intensity injections controls seismic activity, including both $M_L > 2.4+$ events and SWEs. On the contrary, we observed that the reduction in the total injected volume in 2018–2019, with fewer days of intense injections (>390 m³ d⁻¹), was accompanied by a period of reduced seismic activity, with only $2 M_L > 2.4$ events recorded from April 2018 to June 2020 (Fig. 14a).

If we now examine variations in injection rate from one day to another (i.e. injection rate derivative, Fig. 14d), a correlation between variations in daily injections (>135 m³) and activation of SWEs (1, 2, 3, 4; Fig. 14d) seems to be observed. However, the period of important injection variations observed in June 2020 started 3 d after the activation of SWEs and the two main earthquakes (magnitudes of 2.6 and 3.3). Injection operations are therefore not driving the seismicity in this case (but seismicity may have influenced injection operations—the information is not made available). For the other SWEs (1, 2, 3; Fig. 14), we do find injection variations preceding the appearance of SWEs. Injection rate variations (±135 m³ compared to the day before) may therefore also play a role in controlling SWEs in some cases.

Furthermore, the time lag of a few months between intensive injection days and the SWEs corresponds roughly to velocities of tens of metres per day (5–15 m d⁻¹ for the 2016-SWEs and 2018-SWEs, 1, 3, Fig. 14b; 60–80 m d⁻¹ for the 2017-SWEs, 2, Fig. 14b; 40–50 m d⁻¹ for the 2021-SWEs, 5, Fig. 14b), falling within the range of migration velocities observed for fluid pressure diffusion processes (with observed diffusivity coefficients ranging between 0.01 and 10 m s⁻², Shapiro *et al.* 2002; Lohman & McGuire 2007; De Barros *et al.* 2020), while aseismic processes would involve higher velocities on the order of tens of kilometres per day (Shapiro *et al.* 2002; Lohman & McGuire 2007; De Barros *et al.* 2020). With such a diffusion mechanism, it would be possible to explain the appearance of the June 2020 SWEs further (6–7 km) from the wells (4, Fig. 14b) by intense injection days at the end of 2019 and diffusion velocities around 30–40 m d⁻¹ (Fig. 14c). However, the range of observed velocities is quite large between the swarms (variations from 5 to 80 m d⁻¹). It could be due to errors on the well-swarm distance estimation, but it could also come to a more complex mechanism of triggering than a simple fluid pressure diffusion. Indeed, fluid diffusion and aseismic slip may also coexist (e.g. De Barros *et al.* 2020; Hatch *et al.* 2020). Aseismic slip could allow the activation of seismic clusters (e.g. Lohman & McGuire 2007; Danré *et al.* 2022; Yang *et al.* 2023), and/or the fluid pressure primarily induces aseismic slip which, in turn, migrates and triggers the seismicity (e.g. Danré *et al.* 2022; Yang *et al.* 2023). These different possible triggering mechanisms are difficult to study, as event localizations are not well constrained for the Lacq area, especially for the years 2016–2020. Obtaining a more precise quantification of this potential migration pattern of an order of a few tens of metres per day is therefore challenging. For the year 2021 solely, where better resolved event locations are available, a fluid diffusion process from each day with injection above 390 m³
<figure><img src="image_15.png" /><figcaption>**Figure 14.** (a) Number of higher magnitude events ($M_L$ > 2.4+, coloured curve, without SWEs) during 3 months. The colour scales the average distance of the events to the injection wells (in km) (b) Number of SWEs observed during 3 months (c) Number of days with injection above 390 m³, observed during 3 months. (d) During 3 months, the number of days where the injection increases or decreases more than 135 m³ compared to the day before (injection derivative above 135 m³). The blue light rectangles mark the winter periods (November–February). The dashed lines are the beginning of the 2016, 2017, 2018 and 2020 SWEs. The 1, 2, 3, 4 and marks are the 5 maxima that can be noted in the 4 plots.</figcaption></figure>

(mostly between January and April), with a diffusion coefficient of $D = 0.036 \text{ m}^2 \text{ s}^{-1}$ (i.e. diffusion velocity around 40–50 m d$^{-1}$), could be coherent with the activation of CEs/SWEs closer to the wells (2–3 km) for the January–June period, and further (8–9 km) from July to October (Fig. S8). Events which are not CE/SWEs do not seem to obey a diffusion pattern (Fig. S8). However, no definitive relationship can be established based on only 1 yr of analysis.

## 5.4. Relationship between injection and clustered events (CEs) identified by their spatio-temporal proximity (1969–2021)

At first order, an injection at the beginning of the year and a weak diffusion could explain the SWEs mainly observed in summer/autumn. But one could also argue that a seasonal external force such as the water table loading could also explain this seasonality, as observed by Johnson *et al.* (2017). We therefore looked at the seismicity before 2013, in order to quantify any possible seasonality in the seismicity rate. As SWEs cannot be detected from the cross-correlations before 2013 (waveforms are not available), clusters are now identified by their spatio-temporal proximity in our final catalogue. If at least 5 earthquakes are found within 2 km in 10 d, they are considered to form the same cluster of events (CEs). We found 8 CEs after 2013 (which are all part of SWEs) and 23 CEs between 1974 and 2013 (The Geopetrol catalogue was discarded for this CEs identification, because all the events are forced to be located close to

the wells, which does not allow to define a very satisfactory spatial proximity criterion to dissociate the earthquakes into different clusters). It is likely that some clusters have not been identified by this procedure and some may be poorly located. Nevertheless, it seems possible to affirm that these CEs are not systematically detected during the summer periods in the 1975–2013 period (Fig. 15c). This seasonality is only resolved in the 2013–2022 period (Fig. 15b). The precise date of appearance of this seasonality is not easy to estimate. There were several CEs in winter in 1995–1997 (Fig. 15c). We found one cluster in the summer 2005, and then, the clusters are in the spring/summer/autumn period from 2013 until 2021 (Fig. 15b). The tipping point is therefore somewhere between 1997 and 2013.

Seasonal external perturbations, such as the loading of the water table, have no reason to change behaviour between 1997 and 2013. This change seems rather to be explained either by the beginning of injections on LA102 in the 2000s, by changes in injection procedures/techniques/periods around the years 2013–2014 with a change of operator (from 2014, the new operator is the Geopetrol company), or related to the change of the reservoir pressure, which decreases until the years 2013–2014 and increases since then (Fig. 12c). However, details of injections before 2016 are missing to better study a link between this cluster seasonality and injections before 2016. More data on injections, on the distribution of injections between LA102 and LA109, on seismicity, are needed to make progress in understanding the origins of these SWEs/CEs and their apparent seasonality post 2013.
<figure><img src="image_16.png" /><figcaption>**Figure 15.** (a) Distance of Class-1 and Class-2 events, which are not clustered events (CEs) to injection wells (km), according to the day of the year (from 0 to 366). The black line is the average distance for 30 consecutive days. The Colour scale the year. (b) Same as (a) but for CE events from the year 2013, and removing Geopetrol events (see text). Large black triangles are Class-2 SWEs and small black triangles are Class-1 SWEs. (c) Same as (a) and (b) for CE events before the year 2013. (d) Number of injection days above 370 m³ (blue) and above 390 m³ (red) each month, from 2016 to 2021.</figcaption></figure>

## 6 DISCUSSIONS, CONCLUSIONS AND PERSPECTIVES

The review of all existing seismic catalogues for the Lacq area, their improvement and compilation, has revealed that the background seismicity is primarily concentrated within or just below the reservoir boundaries. By deploying a new seismic network in 2021 and analysing temporary networks deployed near Lacq, it was possible to compensate for the lack of instrumentation in the area since the 2000s. These efforts demonstrated that seismicity has consistently been concentrated in the same areas of the reservoir since the 1970s: (1) near the injection and production wells and (2) in a more diffuse band to the north of the reservoir towards Arthez-de-Béarn. These same zones remain active in 2021 (Figs 8 and 13). Few events are found clearly above the top of the reservoir for Class-2 events (Figs 8, 13 and 16). For the 2021 earthquakes, non-SWEs/non-CEs Class-2 events are mainly located either inside the reservoir (37 per cent, Fig. 16a) or 500 m below it (40 per cent, Fig. 16a), and the majority of event localizations consistently follow the deepening of the reservoir. However, a systematic shift of the depth estimation of a few hundred metres (due to velocity model uncertainties) cannot be discarded and an overestimation of the vertical extent of the seismogenic zone (90 per cent of events are confined within 1.2 km, Fig. 16a) is likely due to depth uncertainties around 1–2 km. The seismicity could be therefore slightly shallower and fully inside the 500 m thick reservoir, or fully below it, or overlapping between the reservoir and the lower layers. Therefore, it is likely that most of these events are located on faults extending in/between the reservoir

and the Upper Jurassic layers below it (e.g. Figs 1 and 7). The injection below the gas reservoir could also explain why the seismicity is on average slightly below the gas reservoir (Fig. 16).

The Lacq seismicity is also characterized by the presence of short duration clusters (SWEs and CEs), some around the injection wells, and some peripheral clusters. These SWEs/CEs seem to have similar depth than non-SWEs/non-CEs for the reservoir central part (top reservoir depths <4500 m), but become deeper (6–12 km) at its periphery (for top reservoir depths >4500 m, Fig. 16b). The deep peripheral SWEs/CEs seem therefore not linked to the reservoir faults. The August 2021 deep cluster (at around 7 km: red dots between 2 and 3 km below the top of the reservoir in Fig. 16b, see also Fig. 7), well constrained by the SISLACQ2021 network, is clearly deeper than the reservoir. One can imagine here that the CEs/SWEs are activated on external critically stressed faults, a few kilometres (<15 km) further from the injections and the reservoir, as a temporal correlation between injection and SWEs/CEs did exist. It could be on the North Pyrenean Frontal Thrust system, at the southern part of the reservoir (for the southern deep clusters, Figs 5, 6 and 8), and below the reservoir (for the northwestern deep clusters, Figs 6, 7 and 8). Such activation distances (maximum 15 km) have already been observed (e.g. The 2016 $M_w$ 5 Oklahoma earthquake triggered at more than 40 km from fluid disposal wells, Goebel et al. 2017), and, in several other cases of induced seismicity, events occurred below the injection zones, at depth greater than 10 km (e.g for Oklahoma, California, Texas induced seismicity; Keranen et al. 2013; Hornbach et al. 2015; Goebel et al. 2016; Hennings et al. 2019). It seems reasonable to think that these deep SWEs/CEs
<figure><img src="image_17.png" /><figcaption>**Figure 16.** (a) Depth difference between the top of the reservoir and each earthquake as a function of the depth of the top of the reservoir. The depth difference is represented for Class-2 events, without CEs/SWEs. The colour of the dots scales the catalogue origin (Table 2). The blue lines represent the reservoir top and bottom depths; the dashed blue lines represent the reservoir top and bottom depths with a 500 m variation. The black line represents the average variation of the earthquake depth to the top of the reservoir. (b) Same as (a) but for the CEs and SWEs events.</figcaption></figure>

could be a clue about the activation of external critically stressed faults. Such (re)activation of assumed inactive faults (Choukroune 1992) would impact the seismic hazard map for this Lacq area. It is therefore crucial to confirm such activation, especially if the shallower part of the North Pyrenean Frontal Thrust shows any activity, as suggested by the Dubos-Sallée *et al.* (2006) catalogue (but that we did not confirm during this study). Deploying a denser network, on a longer time period, potentially developing an improved velocity model, and the use of deep borehole sensors could help reduce depth uncertainties and, most importantly, correct any potential systematic location bias. These efforts would definitively address the issue of seismicity localization relative to the reservoir and external faults for the background permanent seismicity and the CE/SWEs.

This study also allowed the most complete view of the variations of seismic activity from 1969 until 2021 in the Lacq area. Although there was a peak in seismicity during the 1980s, the seismic activity has not decreased over time and even appears to be increasing since 2013, despite the end of the intense extraction period. Based on these observations, it appears that the same process has been responsible for this continuous and consistent seismicity, which has persisted for more than 50 yr: the injection of wastewater in the deep gas reservoir emerges as a promising factor, supported by the strong temporal relationship between injection, higher magnitude events, CEs and SWEs.

Examining the relationship between annual injection volumes and the seismic energy released in the reservoir, a positive correlation is observed until the 2010s. This supports the hypothesis that injections have influenced seismicity (more injections correlating with a higher frequency of high-magnitude earthquakes). However, around 2010–2013, a change in this injection-seismicity relationship was observed. Although injection volumes are lower, seismicity remains persistent. Starting in 2016, daily injection data became available, revealing a strong temporal link between recording earthquakes of magnitude greater than 2.4 and days with higher rate/volume injections (approximately above 400 m³ d⁻¹). Thus, even in the most recent years, higher magnitude seismicity appears

to be triggered by injection when it exceeds a certain threshold, as noted by Rayleigh *et al.* (1976).

Furthermore, our analysis reveals that SWEs occur predominantly during the spring, summer, and early autumn from around the years 2000–2013. This seasonality could be attributed to injection activities, as higher injection volumes are typically observed during winter, and a slow diffusion process (tens of metres per day) could explain the delayed activation of SWEs a few months later and at a large distance from the injection wells (<15 km). Other external repetitive factors, such as water table loading, could be considered if an explanation is found for the absence of CEs’ seasonality before 2000–2013 (potentially due to different initial stress conditions).

Studying the individual influence of the two injection wells LA109 and LA102 could provide new insights to refine these observations (seismic energy variations and CEs/SWEs seasonality). LA102 appears to have better connectivity with the entire reservoir, while LA109 remains at higher pressure levels, which supports the hypothesis of a compartmentalized reservoir (CLSIIC4000 2019). A change of porosity and permeability around LA102 could have impacted the seismicity (e.g Yang & Dunham 2021). It would be valuable to investigate whether the use of LA102, starting from 2002, enables better fluid diffusion into other parts of the reservoir and influences the injection-seismicity causality and the seasonality of SWEs.

These new findings raise questions about the relative influence of extraction activities, which have been traditionally considered the primary factor in driving and controlling seismicity. The depletion model fails to replicate the observed variations in seismic energy release over time, particularly the two periods of increasing energy observed after 2000, a period when gas extraction was continuously declining and nearly ceased in 2013.

When using a poro-elastic model, the depletion process gradually moves rock matrices outside the reservoir towards failure (e.g. Segall *et al.* 1994), and a time delay between the onset of depletion and the observed seismicity is expected (e.g. Grasso *et al.* 2021). The initial induced event in 1969, which occurred prior to deep injections, could therefore potentially align with the deep
reservoir depletion model (e.g. Grasso et al. 2021). It could also be linked to shallow injections in the oil reservoir conducted between 1955 and 1974 (phase I, Table 1). However, there is no evidence of seismicity around or inside this shallow oil reservoir, and the seismic energy released during these periods is minimal (Fig. 12). Therefore, a strong impact of shallow injections on seismicity is not expected, and we rather propose that the onset of seismicity is driven by the depletion process, as previously suggested (e.g. Grasso & Wittlinger 1990; Grasso et al. 2021). The Lacq case study may represent a case of induced seismicity where both production and disposal activities contribute to fault reactivation, as proposed by Grasso et al. (2021). It is also worth noting that once earthquakes have initiated, interactions between earthquakes themselves can influence the stress distribution on the fault (Wang et al. 2020). We propose that a portion of the background seismicity is associated with the long-term gradual reservoir depletion process, affecting the overall stress field around the reservoir rocks, while higher magnitude events that coincide with significant changes in injection rates are driven by injections into the deep gas reservoir.

One other future objective could be to better evaluate the exact triggering mechanisms driving the seismicity, from a possible mix of fluid diffusion, aseismic slip and long term depletion. One approach could be to rely on SISLACQ2021 data, covering a longer recording period spanning several years, with consistently reliable event localizations. Above all, a complete reservoir and geomechanical modelling would be necessary to investigate the link between injection, depletion and seismicity, and for more reliable earthquake forecasting, as proposed by Meyer et al. (2023).

## DATA AND RESSOURCES

PYROPE (2013): data obtained from Chevrot, S., Sylvander, M., & RESIF. (2017). *Seismic network X7: PYROPE PYREnean Observational Portable Experiment (RESIF-SISMOB)* [Data set]. RESIF—Réseau Sismologique et géodésique Français. https://doi.org/10.15778/RESIF.X72010.

**FDSN Datacenter access:** https://www.fdsn.org/networks/detail/X7.2010/

MAUPASSAC (2017): data obtained from Lehujeur et al. (2021). Three-dimensional shear velocity structure of the Mauléon and Arzacq Basins (Western Pyrenees). *Bulletin de la Société Géologique de France*, 192(1).

**FDSN Datacenter access:** https://www.fdsn.org/networks/detail/XD.2017/

OMP (2021): catalogue obtained from Sylvander, M., Rigo, A., Sénéchal, G., Battaglia, J., Benahmed, S., Calvet, M., . . . & Pauchet, H. (2021). Seismicity patterns in southwestern France. Comptes Rendus. Géoscience, 353(S1), 79-104.

RENASS (2021): catalogue obtained from the French National seismological network.

**Catalogue access:** https://renass.unistra.fr/fr/zones/

BCSF-Rénass (2022): Instrumental seismicity in mainland France. EOST UAR830, Université de Strasbourg, CNRS. (Collection). doi:10.25577/fv3f-sq09

CEA (2021): catalogue derived from Duverger, C., Mazet-Roux, G., Bollinger, L., Trilla, A. G., Vallage, A., Hernandez, B., & Cansi, Y. (2021). A decade of seismicity in metropolitan France (2010–2019): the CEA/LDG methodologies and observations. Une décennie de sismicité en France métropolitaine (2010–2019): les méthodes et observations du CEA/LDG. *Bulletin de la Société Géologique de France*, 192(1).

**Catalogue access:** https://www-dase.cea.fr/evenement/syntheses_resultat.php?lang=fr

SISLACQ2021 (2021): data obtained from Letort, J., Coton, F., Jacquemond, L., Grimaud, F., Pauchet, H. (2023). SISLACQ2021. GFZ Data Services. Other/Seismic Network. doi:10.14470/0Z082998

**FDSN Datacenter access:** https://www.fdsn.org/networks/detail/5M.2021/

SISLACQ, UPPA (2002-2006): data obtained from Dubos-Sallée, N., Bardainne T., Sénéchal, G. (2006). Lacq gas field seismicity: spatio-temporal evolution over 30 years. Geophysical Research Abstracts, Vol. 8, 03 299, 2006. SRef-ID: 1607-7962/gra/EGU06-A-03299

**Data access:** on demand by email at guy.senechal@univ-pau.fr

EPOS (2018): data from EPOS (2018) Episode: LACQ GAS FIELD,

**Data access:** at https://tcs.ah-epos.eu/#episode:LGF (last accessed December 2022).

Cretace4000 Committee (2023): detailed injection monitoring reports (daily injection data), from the injection operator company Geopetrol, and presented to the local injection monitoring committee of the Crétacé4000 project. These reports have been provided by the French Regional Directorate for Environment (DREAL) to the authors.

**Data access:** on demand by email at DREAL@developpement-durable.gouv.fr

Geopetrol catalogue (2010-2017): data from local ‘industrial’ network, 20 132 016 are available at CLSIC4000 2019, Comité local du suivi des injections Cretacé 4000, rapport quinquenal, 2013-2018, Mai 2019, 51 pp.

**Data access:** at http://www.pyrenees-atlantiques.gouv.fr/content/download/32 169/207 505/file/Rapport per cent20quinquennal per cent20C4000 per cent 202013-2018 per cent20v per cent20finale.pdf (last accessed September 2020).

## DECLARATION OF COMPETING INTERESTS

The authors acknowledge that there are no conflicts of interest.

## ACKNOWLEDGMENTS

We thank the French Regional Directorate for Environment (DREAL) for providing the useful daily injection data. We also thank Noelwenn Dubos-Sallée for careful manual picks made from the SISLACQ1 network that have allowed us to locate events for the 2002–2006 period. Finally, we also thank all the owners who have welcomed the temporary seismic SISLACQ2021 stations in their garden/domain. Thanks to them, we managed to provide crucial event localizations for the recent years. Finally, we deeply thank the 2 reviewers (Professor Heather DeShon & Anonymous) and the Editor (Dr Sidao Ni), as well as the editorial board, for their very useful comments that allowed us to improve a lot this paper.

## SUPPORTING INFORMATION

Supplementary data are available at GJIRAS online.

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by
the authors. Any queries (other than missing material) should be
directed to the corresponding author for the paper.

REFERENCES

Aochi, H. & Burnol, A.(2018). Mechanism of the ML 4.0 25th April 2016 earthquake in southwest of France in the vicinity of the Lacq gas field. *J. Seismol.*, **22**, 1139–1155.

Baisch, S., Ceranna, L. & Harjes, H.P. 2008. Earthquake cluster: what can we learn from waveform similarity? *Bull. seism. Soc. Am.*, **98**(6), 2806–2814.

Bardainne, T. 2005. Étude de la Sismicité De Lacq et Analyse Des Formes D'ondes Par Décomposition En Chirplets, *Doctoral dissertation*, Université de Pau et des Pays de l'Adour (in French).

Bardainne, T., Dubos-Sallée, N., Sénéchal, G., Gaillot, P. & Perroud, H. 2008. Analysis of the induced seismicity of the Lacq gas field (Southern France) and model of deformation, *Geophys. J. Int.*, **172**(3), 1151–1162.

Boyer, E., 1996. *Sismicité induite et production pétrolière*, Final Report Engineering Study, University of Grenoble I, France (in French).

Cara, M., Cansi, Y., Schlupp, A. et al. 2015. SI-Hex: a new catalogue of instrumental seismicity for metropolitan France, *Bull. Soc. géol. France*, **186**(1), 3–19.

Chen, X. & Shearer, P.M. 2011. Comprehensive analysis of earthquake source spectra and swarms in the Salton Trough, California, *J. geophys. Res.*, **116**(B9), doi:10.1029/2011JB008263.

Chen, X., Shearer, P.M. & Abercrombie, R.E. 2012. Spatial migration of earthquakes within seismic clusters in Southern California: evidence for fluid diffusion, *J geophys. Res.*, **117**, doi:10.1029/2011JB008973.

Chevrot, S. & Sylvander, M., & RESIF. 2017. Seismic network X7:PYROPE PYRenean Observational Portable Experiment (RESIFSISMOB) [Data set]. RESIF—Réseau Sismologique et géodésique Français, https://doi.org/10.15778/RESIF.X72010.

Choukroune, P. 1992. Tectonic evolution of the Pyrenees. *Annu. Rev. Earth planet. Sci.*, **20**(1), 143–158.

CLSIC4000 2012. Follow up committee 2012 Cretacé4000 (French only). Comité local du suivi des injections Cretacé 4000, rapport quinquenal, 2007-2012, Octobre 2013, 241 pp., available at https://www.pyrenees-atlantiques.gouv.fr/content/download/32173/207549/file/Rapport_quinquennal2007-2012.pdf (last accessed September 2020).

CLSIC4000 2019. Follow up committee 2019 Cretacé4000 (French only). Comité local du suivi des injections Cretacé 4000, rapport quinquenal, 2013-2018, Mai 2019, pp. 51, available at http://www.pyrenees-atlantiques.gouv.fr/content/download/32169/207505/file/Rapport%20quinquennal%20C4000% 202013-2018%20v%20finale.pdf (last accessed September 2020).

CLSIC4000 2021. Follow up committee 2021 Cretacé4000 (French only), available at https://www.pyrenees-atlantiques.gouv.fr/Actions-de-l-Etat/Cadre-de-vie-eau-environnement-et-risques-majeurs/Risques-technologiques/Cretace-4000/Comite-local-de-suivi-CLS-reunions-2021/Sea-nce-du-26-fevrier-2021

Danré, P., De Barros, L., Cappa, F. & Ampuero, J.P. 2022. Prevalence of aseismic slip linking fluid injection to natural and anthropogenic seismic swarms, *J. geophys. Res.*, **127**(12), e2022JB025571.

De Barros, L., Cappa, F., Deschamps, A. & Dublanchet, P. 2020. Imbricated aseismic slip and fluid diffusion drive aseismic swarm in the Corinth Gulf, Greece, *Geophys. Res. Lett.*, **47**, doi:10.1029/2020GL087142.

Derode, B., et al. 2023. Fluid-driven seismic swarms in the Gripp valley (Haute-Pyrénées, France), *Geophys. J. Int.*, **234**(3), 1903–1915.

Dubos-Sallée, N., Bardainne, T. & Sénéchal, G. 2006. Lacq gas field seismicity: spatio-temporal evolution over 30 years, *Geophys. Res. Abstr.*, **8**, SRef-ID: 1607-7962/gra/EGU06-A-03299.

Duverger, C., Mazet-Roux, G., Bollinger, L., Trilla, A.G., Vallage, A., Hernandez, B. & Cansi, Y. 2021. A decade of seismicity in metropolitan France (2010–2019): the CEA/LDG methodologies and observations, *Bulletin de la Société Géologique de France*, **192**(1), doi:10.1051/bsgf/2021014.

Foulger, G.R., Wilson, M.P., Gluyas, J.G., Julian, B.R. & Davies, R.J. 2018. Global review of human-induced earthquakes, *Earth Sci. Rev.*, **178**, 438–514.

Gao, D., Kao, H. & Wang, B. 2021. Misconception of waveform similarity in the identification of repeating earthquakes, *Geophys. Res. Lett.*, **48**(13), e2021GL092815.

Giardini, D. 2009. Geothermal quake risks must be faced, *Nature*, **462**(7275), 848–849.

Goebel, T.H., Hosseini, S.M., Cappa, F., Hauksson, E., Ampuero, J.P., Aminzadeh, F. & Saleeby, J.B. 2016. Wastewater disposal and earthquake swarm activity at the southern end of the Central Valley, California, *Geophys. Res. Lett.*, **43**(3), 1092–1099.

Goebel, T.H.W., Weingarten, M., Chen, X., Haffener, J. & Brodsky, E.E. 2017. The 2016 Mw5.1 Fairview, Oklahoma earthquakes: evidence for long-range poroelastic triggering at >40 km from fluid disposal wells, *Earth planet. Sci. Lett.*, **472**, 50–61.

Grasso, J.R. & Feignier, B. 1990. Seismicity induced by gas production: II. Lithology correlated events, induced stresses and deformation, *Pure appl. Geophys.*, **134**(3), 427–450.

Grasso, J.R. & Wittlinger, G. 1990. Ten years of seismic monitoring over a gas field, *Bull. seism. Soc. Am.*, **80**(2), 450–473.

Grasso, J.R., Amorese, D. & Karimov, A. 2021. Did wastewater disposal drive the longest seismic swarm triggered by fluid manipulations? Lacq, France, 1969–2016, *Bull. seism. Soc. Am.*, **111**(5), 2733–2752.

Grigoli, F. et al. 2017. Current challenges in monitoring, discrimination, and management of induced seismicity related to underground industrial activities: a European perspective, *Rev. Geophys.*, **55**, 310–340.

Guyoton, F., Grasso, J.R. & Volant, P. 1992. Interrelation between induced seismic instabilities and complex geological structure, *Geophys. Res. Lett.*, **19**(7), 705–708.

Hatch, R.L., Abercrombie, R.E., Ruhl, C.J. & Smith, K.D. 2020. Evidence of aseismic and fluid driven processes in a small complex seismic swarm near Virginia City, Nevada, *Geophys. Res. Lett.*, **47**(4), e2019GL085477.

Helmstetter, A., Nicolas, B., Comon, P. & Gay, M. 2015. Basal icequakes recorded beneath an Alpine glacier (Glacier d'Argentière, Mont Blanc, France): evidence for stick-slip motion?, *J. geophys. Res.*, **120**(3), 379–401.

Hennings, P.H., et al. 2019. Injection-induced seismicity and fault-slip potential in the Fort Worth Basin, Texas, *Bull. seism. Soc. Am.*, **109**(5), 1615–1634.

Hornbach, M.J. et al. 2015. Causal factors for seismicity near Azle, Texas, *Nat. Commun.*, **6**, doi:10.1038/ncomms7728.

Improta, L., Bagh, S., De Gori, P., Valoroso, L., Pastori, M., Piccinini, D. & Buttinelli, M. 2017. Reservoir structure and wastewater-induced seismicity at the Val d'Agri Oilfield (Italy) shown by three-dimensional vp and vp/vs local earthquake tomography, *J. geophys. Res.*, **122**(11), 9050–9082.

Johnson, C.W., Fu, Y. & Bürgmann, R. 2017. Seasonal water storage, stress modulation, and California seismicity, *Science*, **356**(6343), 1161–1164.

Keranen, K. & Weingarten, M. 2018. Induced seismicity, *Annu. Rev. Earth planet. Sci.*, **46**, 149–174.

Keranen, K.M., Savage, H.M., Abers, G.A. & Cochran, E.S. 2013. Potentially induced earthquakes in Oklahoma, USA: links between wastewater injection and the 2011 Mw 5.7 earthquake sequence, *Geology*, **41**(6), 699–702.

Lahaie, F. & Grasso, J. R. 1999. Loading rate impact on fracturing pattern: Lessons from hydrocarbon recovery, Lacq gas field, France, *J. Geophys. Res: Solid Earth*, **104**(B8), 17941–17954

Laporte, M. 2022. Contribution à l'amélioration de l'estimation de la profondeur hypocentrale à partir de réseaux régionaux ou globaux, *Phd thesis*.

Lehujeur, M., Chevrot, S., Villaseñor, A., Masini, E., Saspiturry, N., Lescoutre, R. & Sylvander, M. 2021. Three-dimensional shear velocity structure of the Mauléon and Arzacq Basins (Western Pyrenees), *Bulletin de la Société Géologique de France*, **192**(1), 47, doi:10.1051/bsgf/2021039.

Letort, J., Cotton, F., Jacquemond, L., Grimaud, F. & Pauchet, H. 2023. GFZ Data Services. Other/Seismic Network. doi:10.14470/Z082998.
Lohman, R.B. & McGuire, J.J. 2007. Earthquake swarms driven by aseis-
mic creep in the Salton Trough, California, *J. geophys. Res.*, **112**(B4),
doi:10.1029/2006JB004596.

Lomax, A., Virieux, J., Volant, P. & Berge-Thierry, C. 2000. Probabilistic earthquake location in 3D and layered models, in *Advances in Seismic Event Location. Modern Approaches in Geophysics*, Vol. **18**, pp. 101–134, eds Thurber, C.H. & Rabinowitz, N., Springer.

Maury, V.M.R., Grasso, J.R. & Wittlinger, G. 1992. Monitoring of subsi-
dence and induced seismicity in the Lacq gas field (France): the conse-
quences on gas production and field operation, *Eng. Geol.*, **32**(3), 123–
135.

Meyer, H., Smith, J.D., Bourne, S. & Avouac, J.P. 2023. An integrated framework for surface deformation modelling and induced seismicity forecasting due to reservoir operations, *Geol. Soc., Lond., Spec. Publ.*, **528**(1), 299–318.

Nadeau, R.M. & McEvilly, T.V. 1999. Fault slip rates at depth from re-
currence intervals of repeating microearthquakes, *Science*, **285**(5428),
718–721.

Nur, A. 1972. Dilatancy, pore fluids, and premonitory variations of ts/tp travel times, *Bull. seism. Soc. Am.*, **62**(5), 1217–1222.

Pauchet, H., Rigo, A., Rivera, L. & Souriau, A., 1999. A detailed analysis of the February 1996 aftershock sequence in the eastern Pyrenees, France, *Geophys. J. Int.*, **137**(1), 107–127.

Rayleigh, C.B., Healy, J.H. & Bredehoeft, J.D. 1976. An experiment in earthquake control at Rangely, Colorado, *Science*, **191**(4233), 1230–1237.

Rigo, A., et al. 2015. Present-day deformation of the Pyrenees revealed by GPS surveying and earthquake focal mechanisms until 2011, *Geophys. J. Int.*, **201**(2), 947–964.

Roecker, S.W. 1982. Velocity structure of the Pamir-HinduKush region: pos-
sible evidence of subducted crust, *J. geophys. Res.*, **87**, 945–959.

Roland, E. & McGuire, J.J. 2009. Earthquake swarms on transform faults, *Geophys. J. Int.*, **178**, 1677–1690.

Rothé, J.P. 1970. Seismes artificiels, *Tectonophysics*, **9**(2/3), 215–238 (in French).

Rothé, J.P. 1977. Seismes artificiels et exploitations petrolières: l’exemple de Lacq (France), *Ann. Geofis.*, **30**, 369–383 (in French).

Schmittbuhl, J. et al. 2021. Induced and triggered seismicity below the city of Strasbourg, France from November 2019 to January 2021, *Comp. Rend. Géosci.*, **353**(S1), 561–584.

Segall, P., Grasso, J.R. & Mossop, A. 1994. Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France. *J. geophys. Res.*, **99**(B8), 15 423–15 438.

Shapiro, S.A., Rothert, E., Rath, V. & Rindschwentner, J. 2002. Characteri-
zation of fluid transport properties of reservoirs using induced microseis-
micity, *Geophysics*, **67**, 212–220.

Smith, J.D., Avouac, J.P., White, R.S., Copley, A., Gualandi, A. & Bourne, S. 2019. Reconciling the long-term relationship between reservoir pore pressure depletion and compaction in the Groningen region, *J. geophys. Res.*, **124**(6), 6165–6178.

Sylvander, M., et al. 2021. Seismicity patterns in southwestern France, *Compt. Rend. Géosci.*, **353**(S1), 79–104.

Thelen, W., Malone, S. & West, M. 2011. Multiplets: their behavior and utility at dacitic and andesitic volcanic centers, *J. geophys. Res.*, **116**(B8), doi:10.1029/2010JB007924.

Tselentis, G.A., Martakis, N., Paraskevopoulos, P. & Lois, A. 2011. High-resolution passive seismic tomography for 3D velocity, Poisson’s ratio *v*, and P-wave quality QP in the Delvina hydrocarbon field, southern Albania, *Geophysics*, **76**(3), B89–B112.

Van Thienen-Visser, K. & Breunese, J.N. 2015. Induced seismicity of the Groningen gas field: history and recent developments, *Leading Edge*, **34**(6), 664–671.

Vidale, J.E. & Shearer, P.M. 2006. A survey of 71 earthquake bursts across southern California: exploring the role of pore fluid pressure fluctuations and aseismic slip as drivers, *J. geophys. Res.*, **111**(B5), doi:10.1029/2005JB004034.

Volant, P. 1993. Mécanisme des déformations et aspect fractal de la sismicité induite par l’exploitation d’un gisement d’hydrocarbures (Lacq, France), *Doctoral dissertation*, Université Joseph Fourier (in French).

Volant, P. & Grasso, J.R. 1994. The finite extension of fractal geometry and power-law distribution of shallow earthquakes: a geomechanical effect, *J. geophys. Res.*, **99**, 21 879–21 889.

Waldhauser, F. & Ellsworth, W.L. 2000. A double-difference earthquake location algorithm: method and application to the northern Hayward fault, *Bull. seism. Soc. Am.*, **90**, 1353–1368.

Wang, J., Li, T., Gu, Y.J., Schultz, R., Yusifbayov, J. & Zhang, M. 2020. Sequential fault reactivation and secondary triggering in the March 2019 Red Deer induced earthquake swarm, *Geophys. Res. Lett.*, **47**(22), e2020GL090219.

Wittlinger, G. 1980. Etude De la Sismicité En Champ Proche par Un Réseau Sismologique à Faible Ouverture: application au Frioul (Italie) et au Gisement De Lacq (France), *Doctoral dissertation*, University of Strasbourg, France (in French).

Wittlinger, G. 1980. Etude de la sismicité en champ proche par un réseau sis-
mologique à faible ouverture: application au Frioul (Italie) et au gisement
de Lacq (France). *France: University of Strasbourg*.

Yang, Y. & Dunham, E.M. 2021. Effect of porosity and permeability evolution on injection-induced aseismic slip, *J. geophys. Res.*, **126**(7), e2020JB021258.

Yang, Y., Yang, H. & Zi, J. 2023. Stress transfer outpaces injection-induced aseismic slip and triggers seismicity, **13**, doi:10.1038/s41598-023-43760-0.

Zaliapin, I. & Ben-Zion, Y. 2016. A global classification and characterization of earthquake clusters, *Geophys. J. Int.*, **207**(1), 608–634.

Zhai, G., Shirzaei, M., Manga, M. & Chen, X. 2019. Pore-pressure dif-
fusion, enhanced by poroelastic stresses, controls induced seismicity in
Oklahoma, *Proc. Natl. Acad. Sci.*, **116**(33), 16 228–16 233.