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

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³ d⁻¹) and both clustered events and higher magnitude earthquakes (greater than 2.4).

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 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 ductile 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 1990 s (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 seismicty 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 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 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 experiment (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 improve 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 ± 0.1 for local stations in the Lacq area, and 1.8 ± 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 University of Strasbourg (IPGS) installed a network of seven, 1-component stations and one, 3-component station (Wittlinger 1980). The network 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 processed to build up two earthquake catalogues; (i) the Volant catalogues, 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 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 NonLinLoc 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 tP*0.95 + 0.2 s, or below tP*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 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 (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.

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 of117 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 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.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, 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, 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:

• NLL horizontal uncertainty lower than 1.5 km.

• NLL vertical uncertainty lower than 2 km.

• At least 6 phases used.

• RMS value lower than 0.2 s.

• Azimutal 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.

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:

• NLL horizontal uncertainty lower than 2 km.

• NLL vertical uncertainty lower than 3 km.

• At least 6 phases used.

• RMS value lower than 0.3 s.

• Azimutal 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 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 earthquakes 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 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³ (mostly between January and April), with a diffusion coefficient of D = 0.036 m² s⁻¹ (i.e. diffusion velocity around 40–50 m d⁻¹), 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.

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 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 (CLSIC4000 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 obtain 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., Cotton, 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.

Acknowledgement

We thank the French Regional Directorate for Environment (DREAL) for providing the useful daily injection data. We also thank Noalwenn 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.

References

Author notes