Jean-Robert Grasso
email: jrgrasso@me.com

Daniel Amorese

Abror Karimov

Did Wastewater Disposal Drive the Longest Seismic Swarm Triggered by Fluid Manipulations? Lacq, France, 1969-2016

The activation of tectonics and anthropogenic swarms in time and space and size remains challenging for seismologists. One remarkably long swarm is the Lacq swarm. It has been ongoing since 1969 and is located in a compound oil-gas field with a complex fluid manipulation history. Based on the overlap between the volumes where poroelastic model predicts stresses buildup and those where earthquakes occur, gas reservoir depletion was proposed to control the Lacq seismic swarm. The 2016 M w 3.9, the largest event on the site, is located within a few kilometers downward the deep injection well. It questions the possible interactions between the 1955-2016 wastewater injections and the Lacq seismicity. Revisiting 60 yr of fluid manipulation history and seismicity indicates that the impacts of the wastewater injections on the Lacq seismicity were previously underevaluated. The main lines of evidence toward a wastewater injection cause are (1) cumulative injected volume enough in 1969 to trigger M w 3 events, onset of Lacq seismicity; (2) 1976 injection below the gas reservoir occurs only a few years before the sharp increase in seismicity. It matches the onset of deep seismicity (below the gas reservoir, at the injection depth); (3) the (2007-2010) 2-3 folds increase in injection rate precedes 2013, 2016 top largest events; and (4) 75% of the 2013-2016 events cluster within 4-8 km depths, that is, close to and downward the 4.5 km deep injection well. As quantified by changepoint analysis, our results suggest that timely overlaps between injection operations and seismicity patterns are as decisive as extraction operations to control the Lacq seismicity. The seismicity onset is contemporary to cumulative stress changes (induced by depletion and injection operations) in the 0.1-1 MPa range. The interrelation between injection and extraction is the most probable cause of the Lacq seismicity onset and is sustenance over time. The injected volume-largest magnitude pair for Lacq field is in the same range (90% confidence level) than wastewater volume-magnitude pairs reported worldwide, in a wide variety of tectonic settings.

KEY POINTS

• We highlight the role of wastewater disposal on seismicity in a context of gas reservoir depletion.

• The interrelation between injection and extraction is the most probable cause of the Lacq seismicity.

• The Lacq volume-magnitude pair is in the same range (90% confidence level) than pairs reported worldwide.

INTRODUCTION

Understanding and managing induced and triggered seismicity are a crucial challenge for many geo-resource applications. Not only the changes in the perception of societies and regulators, but also new kinds of applications, are required to improve the scientific understanding of induced seismicity and to validate frameworks for risk governance (e.g., [START_REF] Muntendam-Bos | A guideline for assessing seismic risk induced by gas extraction in the Netherlands[END_REF][START_REF] Langenbruch | Value at induced risk: Injection-induced seismic risk from lowprobability, high-impact events[END_REF]. Most of the current concerns related to geo-resource production are due to fluid manipulation. The dramatic rise in seismicity in intraplate regions of North America (e.g., Oklahoma, Western Canada) is an unintended consequence of fluid injection. These injections are mainly related to (1) the disposal of wastewater by injection below the productive reservoir (e.g., central United States, [START_REF] Ellsworth | Injection-induced earthquakes[END_REF]; and (2) hydraulic fracturing (e.g., increase in seismicity rate in the western Canada sedimentary basin, [START_REF] Atkinson | Hydraulic fracturing and seismicity in the Western Canada Sedimentary Basin[END_REF][START_REF] Bao | Fault activation by hydraulic fracturing in western Canada[END_REF]. For the largest earthquake in central U.S. context, high-resolution hypocenter locations show that the induced earthquakes principally activate previously unknown faults in the basement, well oriented for failure in the contemporary stress field. The time and space evolutions of earthquake swarms are difficult to deterministically relate to specific well history, due to the current inability to track fluid flow over space and time (e.g., [START_REF] Schoenball | A systematic assessment of the spatiotemporal evolution of fault activation through induced seismicity in Oklahoma and Southern Kansas[END_REF][START_REF] Schoenball | How faults wake up: The Guthrie-Langston, Oklahoma earthquakes[END_REF]. Our study focuses on a very long, still ongoing seismic swarm recorded in southwest France. The Lacq seismic swarm overlaps with a gas and oil field, with two wells used for deep wastewater disposal. It is now the longest ever seismic swarm related to hydrocarbon recovery worldwide. In the Lacq area, the most significant shock (M w 3.9, [START_REF] Aochi | Mechanism of the M L 4.0 25 April 2016 earthquake in southwest of France in the vicinity of the Lacq gas field[END_REF] occurred in 2016, that is, 60 yr after the fluid manipulation onset.

There are two vintages of probabilistic seismic hazard analysis for France: 2002 and 2020. Neither of them excluded induced events from their analysis. The first French national probabilistic hazard map [START_REF] Martin | Preliminary probabilistic seismic hazard assessment of France[END_REF] is the map that was translated into the current regulatory documents and building code. The Lacq field located in a low-seismicity area is called "Nord Pyrénéen" domain [START_REF] Martin | Preliminary probabilistic seismic hazard assessment of France[END_REF]. In this zone, the maximum observed magnitude (magnitude of the reference earthquake) is 5.2, and the maximum magnitude upper bound used for the probabilistic model is 5:2 1 6:2. In the Lacq area, seismic building code is required for most edifices, such as collective buildings and buildings essential for civil security. In 2020, about 20 "Seveso-classified" companies are listed in the region (in the European Union [EU] regulation, "Seveso-classified" sites are hazardous industrial facilities where dangerous substances are used or stored in large quantities; EU-Directive, 2012).

The Lacq seismic swarm is well isolated from the North Pyrenean fault seismicity; this latter clusters 30 km southward (e.g., [START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF][START_REF] Rigo | Present-day deformation of the Pyrenees revealed by GPS surveying and earthquake focal mechanisms until 2011[END_REF]. Deformation rate through the North Pyrenean fault is not well resolved (less than 1 mm yr -1 ) from geodetic surveys [START_REF] Rigo | Present-day deformation of the Pyrenees revealed by GPS surveying and earthquake focal mechanisms until 2011[END_REF]. Geodetic and seismic estimates support a trans-tension regime on the North Pyrenean fault seismicity [START_REF] Rigo | Present-day deformation of the Pyrenees revealed by GPS surveying and earthquake focal mechanisms until 2011[END_REF].

Figure 1 demonstrates that the Lacq seismic swarm is the only swarm that emerges in the whole North Pyrenees foreland as resolved by each of the regional seismicity surveys that deployed dense temporary network since the 1980s (e.g., [START_REF] Gallart | Relationship between deep structure and seismicity in the western Pyrenees[END_REF][START_REF] Chevrot | A preliminary catalog of moment tensors for the Pyrenees[END_REF][START_REF] Rigo | Present-day deformation of the Pyrenees revealed by GPS surveying and earthquake focal mechanisms until 2011[END_REF]. This pattern also exists in Figure 1b when using 1962-2016 (Si-Hex), French earthquake catalog [START_REF] Cara | SI-Hex: A new catalogue of instrumental seismicity for metropolitan France[END_REF]. The onset and the end and the M max value for this swarm remain open questions. Numerous studies (e.g., [START_REF] Rothé | Seismes artificiels[END_REF][START_REF] Rothé | Seismes artificiels et exploitations petrolieres: l'exemple de Lacq (France)[END_REF][START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF]) robustly demonstrate that there was no felt historical earthquake in the Lacq field area up to 1969. After the occurrence of a pair of events in 1969-1970, there was no significant (M w > 2:5) earthquake up to 1976. The drastic increase in M w 2.5+ events started in 1976, reaching a peak value of 8 events/yr in the 1990s. Since 2000, the yearly M w 2.5 event rate for the Lacq swarm decreases to 1-2 events=yr (Fig. 2). The largest energy release (M w 3.9, 2016;[START_REF] Aochi | Mechanism of the M L 4.0 25 April 2016 earthquake in southwest of France in the vicinity of the Lacq gas field[END_REF] occurs 4 yr after the peak of shallow and deep injection rates for wastewater disposal. It questions the relationship between seismicity and local fluid manipulations, in particular, with regard to the M max value and timing (Fig. 3). The new seismic and industrial data we analyze support that the history of wastewater disposal is as decisive as the depletion history to trigger and control the Lacq seismicity. Our analysis used a changepoint method to investigate the relationship between gas and oil production, shallow and deep injection and seismicity, respectively.

DATA

Earthquake catalogs, overview of the Lacq seismicity history, and past scientific results

Data from joined academic and industrial seismic networks that operated on the gas field during 1975-1997 were used to cross-analyze M w > 1:6 seismicity patterns over time and space with the hydrocarbon production history [START_REF] Lahaie | Loading rate impact on fracturing pattern: Lessons from hydrocarbon recovery, Lacq gas field, France[END_REF][START_REF] Bardainne | Etude de la sismicité de Lacq et analyse des formes d'ondes par décomposition en chirplets[END_REF][START_REF] Bardainne | Analysis of the induced seismicity of the Lacq gas field (Southwestern France) and model of deformation[END_REF]. Data from the local seismic network, as operated by industrial partners since 2000, are available during the 2013-2016 period solely (CLSIC4000, 2019, see Data and Resources). To robustly analyze the half-century duration of the ongoing Lacq seismic swarm, we are bound to use (Si-Hex) French national earthquake catalog [START_REF] Cara | SI-Hex: A new catalogue of instrumental seismicity for metropolitan France[END_REF]. The Si-Hex catalog merges all the available arrival times from regional seismic networks to extract the best estimates for the location and size of the 1962-2009 events [START_REF] Cara | SI-Hex: A new catalogue of instrumental seismicity for metropolitan France[END_REF][START_REF] Cara | Magnitude M w in metropolitan France[END_REF]. Respecting the same Si-Hex criteria, we extend the database in the vicinity of the Lacq gas field up to 2016. For the Si-Hex earthquake catalog, the completeness value for the magnitude distribution (M c ) is of the order of M c 2:2, as derived from the frequency-magnitude distribution (FMD; see Fig. 3 for M c definition). Because of network geometry variability back in the 1970s and for analysis robustness, we choose a conservative M c 2:5 value for this study (Fig. 3). The largest local event occurred in 2016 as an M w 3.9 event with a 5:3 0:5 km depth from local network data. The depth accuracies from the 2013-2016 local catalog are within 0.5 km, due to borehole seismometers, all within 300-570 m depth range (Fig. 4).

The Si-Hex catalog we use revisits and homogenizes the magnitude attribute for each event over time. For this catalog, the depth attributes are not accurate enough (the accuracy of the depth estimate is absent) to identify any possible depth changes. Alternatively, using local seismic network data, several studies report deepening of seismicity over time [START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF][START_REF] Guyoton | Interrelation between induced seismic instabilities and complex geological structure[END_REF][START_REF] Lahaie | Loading rate impact on fracturing pattern: Lessons from hydrocarbon recovery, Lacq gas field, France[END_REF][START_REF] Bardainne | Etude de la sismicité de Lacq et analyse des formes d'ondes par décomposition en chirplets[END_REF][START_REF] Bardainne | Analysis of the induced seismicity of the Lacq gas field (Southwestern France) and model of deformation[END_REF]. [START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF] and [START_REF] Guyoton | Interrelation between induced seismic instabilities and complex geological structure[END_REF] resolve the first deep earthquakes (below the reservoir level) that appeared during late December 1982. Since that time, the swarming style of the deeper seismicity contrasts with the more diffuse organization of shallow earthquakes [START_REF] Guyoton | Interrelation between induced seismic instabilities and complex geological structure[END_REF][START_REF] Volant | The finite extension of fractal geometry and power law distribution of shallow earthquakes: A geomechanical effect[END_REF][START_REF] Bardainne | Analysis of the induced seismicity of the Lacq gas field (Southwestern France) and model of deformation[END_REF]. The data from the local seismic network operated by the producer points on 70% of the 2007-2012 seismicity is located in between 3 and 5 km depth, the most significant event on this period occurring at 7 km depth (CLSIC4000, 2013, see Data and Resources).

-2˚00' -1˚30' -1˚00' -0˚30' 0˚00' 0˚30' 1˚00' 1˚30' 2˚00' 2˚30' 3˚00' 42˚00' 42˚30' 43˚00'
During 2013-2016, the depth histogram of the 810 events (-0:4 ≤ M w ≤ 4) located by the producer network (CLSIC4000, 2019, see Data and Resources) evidences that 75% of the seismicity is deeper than 4 km, with a mean depth value of 4.6 km (Fig. 4). Apart from the late occurrence of the largest 2016 event, the main stages of the Lacq seismicity history develop as follows:

1969 corresponds to the onset of the Lacq swarm as the onset of felt seismicity within a several century long aseismic zone [START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF]; 1976 corresponds to a drastic increase of M w 2.5 event rates for a peak value of 8 events/yr in the 1990s. Since 2000, the yearly M w 2.5 event rate for Lacq swarm decreases to 1-2 events=yr (Figs. 2 and3). When using the local seismic network data, possible migration for seismicity was suggested around the Lacq field, 1976[START_REF] Lanzante | Resistant, robust and non-parametric techniques for the analysis of climate data: Theory and examples, including applications to historical radiosonde station data[END_REF](e.g., Bardainne et al., 2008). Using the 1969-2016 data set, we did not resolve lateral migration. On 1976-1996, the key pattern resolved by [START_REF] Bardainne | Analysis of the induced seismicity of the Lacq gas field (Southwestern France) and model of deformation[END_REF] is the decrease of seismicity in the central part of the field, except for the seismic cluster close to the deep injection well. Since 2000, the seismic clustering close to the deep injection well is enhanced. It climaxes on 2013-2016 in which most of the earthquake cluster around the deep injection well (Fig. 4b).

History of Lacq hydrocarbon field operations and previous studies

The Lacq hydrocarbon fields (Fig. 4) are compound oil and gas reservoirs that are trapped within an anticline structure (e.g., [START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF][START_REF] Maury | Monitoring of subsidence and induced seismicity in the Lacq gas field (France): The consequences on gas production and field operation[END_REF][START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF]. The productions for the (0.7 km) shallow oil field and the (3-5 km) deep gas reservoir started in 1950 and 1957, respectively [START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF][START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF].

Oil field production: The oil extraction (cumulative 1950-2016 production 3 × 10 6 tons; standard condition) was balanced by water inflow toward the reservoir, this later being connected to a large aquifer [START_REF] Rothé | Seismes artificiels et exploitations petrolieres: l'exemple de Lacq (France)[END_REF]. It inhibits strong fluctuation of oil reservoir pressure over time. The initial 1950 6.1 MPa pressure changed in between a 1955 minimum 5.9 MPa value and a maximum 7.1 MPa value in 2010, during a reactivation phase of shallow injection in the oil reservoir.

Gas field production: The 1957-2012 gas field production drives a cumulative 64.5 MPa depletion within the gas reservoir rocks (Fig. 2). It corresponds to a cumulative production of 254:25 × 10 9 m 3 in standard condition (Fig. 5). The initial gas pressure was 66.1 MPa at 3.7 km depth (e.g., [START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF]. The production formally stopped in 2012, when the reservoir pressure was of the order of 1.5 MPa (CLSIC4000, 2019, see Data and Resources). There is still a vanishingly close to zero production value since 2013 for local chemical industry use only.

Wastewater disposal: On the Lacq hydrocarbon fields, wastewater disposals through shallow (primarily up to 1975) and 4.8 km deep (1974 onset) injection take place (e.g., [START_REF] Rothé | Seismes artificiels et exploitations petrolieres: l'exemple de Lacq (France)[END_REF][START_REF] Wittlinger | Etude de la sismicité en champ proche par un réseau sismologique à faible ouverture: Application au Frioul (Italie) et au gisement de Lacq[END_REF][START_REF] Grasso | Seismicity induced by gas production: II. Lithology correlated events, induced stresses and deformation[END_REF][START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF][START_REF] Maury | Monitoring of subsidence and induced seismicity in the Lacq gas field (France): The consequences on gas production and field operation[END_REF][START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF][START_REF] Bardainne | Etude de la sismicité de Lacq et analyse des formes d'ondes par décomposition en chirplets[END_REF]. The injection processes started in 1957, as wastewater disposal from four wells (Figs. 1 and4) within the formerly depleted shallow oil field (e.g., [START_REF] Wittlinger | Etude de la sismicité en champ proche par un réseau sismologique à faible ouverture: Application au Frioul (Italie) et au gisement de Lacq[END_REF]. Yearly injection rates in the shallow oil field (0.7 km depth) are in the 6 × 10 5 m 3 =yr magnitude range Yearly rate for pressure changes (i.e., derivative of (a)) and seismicity energy. Seismicity and seismic energy (gray bars) is M w ≥ 2:5 within 1L r distance from the gas reservoir. Peak values are 3.5 and 0:3 MPa=yr for yearly gas and oil reservoir pressures, respectively. Note that a 64.5 MPa gas reservoir depletion value induces a maximum 0.2 MPa poroelastic stress increase in the rock matrices outside the gas reservoir (e.g., [START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF]. The 0.6-1 MPa pore-pressure increase in the shallow oil field corresponds to a direct effective pore-pressure change. The color version of this figure is available only in the electronic edition.

during the 1957-1967 period, with an increase up to 9 × 10 5 m 3 =yr in the 1970-1975 period (Table 1; Fig. 5). Between 1957 and 1975, a cumulative 13 × 10 6 m 3 of wastewater are injected through the oil reservoir (Fig. 5). It causes the pressure to overcome the initial oil field pressure (0.6 MPa increase) to reach a 6.5 MPa peak value by 1975 (Fig. 2). Subsequently, the oil field pressure declined back to 6.2 MPa by 1989. A new phase of pressure increases from 6.6 to 7.1 MPa exists during 2007-2010. The decrease of the oil field pressure back to a 6.6 MPa value is reached by 2012. The corresponding injection rates are in the range of 900-500 m 3 =day in between 2007 and 2012, respectively (Fig. 5). Because of the impermeable cap rock that traps and builds up the gas reservoir, there is a weak probability for the shallow injection to reach the rock matrices in and below the gas reservoir rocks (e.g., [START_REF] Maury | Monitoring of subsidence and induced seismicity in the Lacq gas field (France): The consequences on gas production and field operation[END_REF]. The injected volumes we report in this study are the minimum values that correspond to controlled wastewater disposals over time. Several other production operations, including but not restricted to urgency responses to production constraints, involve fluid injections whose volumes and timings are not currently available (e.g., the 1996 injection during a well collapse; [START_REF] Bardainne | Characterization of seismic waveforms and classification of seismic events using chirplet atomic decomposition. Example from the Lacq gas field (Western Pyrenees, France)[END_REF]. Since 1975, the wastewater disposal is dominated by 4.8 km depth injections (Fig. 4), that is, below the gas reservoir rocks [START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF][START_REF] Maury | Monitoring of subsidence and induced seismicity in the Lacq gas field (France): The consequences on gas production and field operation[END_REF][START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF][START_REF] Bardainne | Etude de la sismicité de Lacq et analyse des formes d'ondes par décomposition en chirplets[END_REF]. Using two nearby reformed production wells, the yearly injection rates are in the 1-3 × 10 5 m 3 =yr magnitude range. The headand bottom-well pressures are of the order of 8 and 10 MPa (at 3700 m=Nm), respectively (CLSIC4000, 2019). For clarity, we define four fluid manipulation phases since the onset of fluid manipulation (Table 1). For all the I-III phases, there is a common gas extraction activity. 1955-1974, phase I corresponds to (0.7 km depth) wastewater injections; 1974-2006, phase II is characterized by pure (4.5 km) deep wastewater injection; 2006-2012, phase III corresponds to simultaneous increases of deep and shallow wastewater injection rates; 2013-2016, phase IV is a deep wastewater injection phase with negligible gas extraction (Fig. 5). The phase II-IV injection rates are three times smaller than during the phase I (Fig. 5b).

METHOD

There is no consensus on the a priori distance-time-size patterns to be selected to relate a given earthquake to a given geo-resource production style, including fluid extraction or injection (e.g., [START_REF] Grasso | Mechanics of seismic instabilities induced by the recovery of hydrocarbons[END_REF][START_REF] Grasso | Le role des fluides dans les instabilitées de la croute supérieure: l' exemple des exploitations d'hydrocarbures[END_REF][START_REF] Grasso | Triggered earthquakes as stress gauge: Implication for the uppercrust behavior in the Grenoble area, France[END_REF][START_REF] Davies | Induced seismicity and hydraulic fracturing for the recovery of hydrocarbons[END_REF][START_REF] Klose | Mechanical and statistical evidence of the causality of human-made mass shifts on the Earth's upper crust and the occurrence of earthquakes[END_REF]National Research Council [NRC], 2013;[START_REF] Grasso | Patterns of reservoir-triggered seismicity in a low-seismicity region of France[END_REF]. This is due to the fact that the understanding of triggering processes for earthquakes is hindered by the poor knowledge of fracture and fault networks. Several studies report the use of the absolute distances to capture earthquake-earthquake triggering patterns (as aftershockmainshock cascades) that hide the key properties of the earthquake interactions (e.g., [START_REF] Bak | Unified scaling law for earthquakes[END_REF][START_REF] Parsons | On near-source earthquake triggering[END_REF][START_REF] Tahir | The largest aftershock: How strong, how far away, how delayed?[END_REF][START_REF] De Arcangelis | Statistical physics approach to earthquake occurrence and forecasting[END_REF]. This way, -0.7°-0.65°-0.6°-0.55°-0.5°4 following [START_REF] Grasso | Patterns of reservoir-triggered seismicity in a low-seismicity region of France[END_REF][START_REF] Grasso | Anthropogenic seismicity as aftershocks for geo-resource production? Implications for M max estimates (reservoir impoundment cases)[END_REF] for the analysis of lake reservoirs seismicity, we analyze in this study the seismicity triggered around the Lacq hydrocarbon reservoirs using normalized distances to the gas reservoir as, 1L r , 3L r , and 10L r distances, L r being the representative reservoir length (Fig. 1). Such an approach is driven by tectonic earthquake analysis in which the static earthquake interactions exist over areas that scale with the earthquake size, that is, in the 1L-3L ranges for aftershocks triggering [START_REF] Parsons | On near-source earthquake triggering[END_REF][START_REF] Tahir | The largest aftershock: How strong, how far away, how delayed?[END_REF][START_REF] Tahir | Aftershock patterns of M s > 7 earthquakes in the India-Asia collision belt: Anomalous results from the Muzaffarabad earthquake sequence, Kashmir, 2005[END_REF]. For the tectonic earthquake analysis, (L) is the mainshock fault length as derived from the magnitude scale (e.g., [START_REF] Wells | New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement[END_REF].

To mimic the method used for tectonic earthquake interaction, we choose in this study (L r ) the gas reservoir length, the characteristic dimension that drives the volumetric stress-change pattern induced by reservoir depletion (e.g., [START_REF] Segall | Earthquakes triggered by fluid extraction[END_REF][START_REF] Grasso | Mechanics of seismic instabilities induced by the recovery of hydrocarbons[END_REF][START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF]. This way, the volumetric stress changes induced by the reservoir depletion mimics the stress change induced by each of the tectonic earthquakes. This is equivalent to consider the reservoir pore-pressure change as a slow earthquake, the production-induced seismicity being the corresponding aftershocks (e.g., [START_REF] Grasso | Patterns of reservoir-triggered seismicity in a low-seismicity region of France[END_REF][START_REF] Grasso | Anthropogenic seismicity as aftershocks for geo-resource production? Implications for M max estimates (reservoir impoundment cases)[END_REF]. We test 1L r -distance as near-field reservoir effect for triggering and 10L r -distance as the far-field distance with the null effect from reservoir-induced stress changes on seismicity, respectively (Figs. 1 andA2). Within these normalized distance ranges, we use changepoint analyses to compare and interrelate hydrocarbon recovery operations and seismicity. Such an approach has already been successfully applied to analyze the evolution of earthquake series in the context of wastewater injection [START_REF] Gupta | A Bayesian change point model to detect changes in event occurrence rates, with application to induced seismicity[END_REF][START_REF] Fiedler | Multiple change-point detection in spatiotemporal seismicity data[END_REF]. These authors implemented a Bayesian changepoint model to assess whether temporal features of observed earthquakes between 1974 and 2014 [START_REF] Gupta | A Bayesian change point model to detect changes in event occurrence rates, with application to induced seismicity[END_REF] or between 1980 and 2015 [START_REF] Fiedler | Multiple change-point detection in spatiotemporal seismicity data[END_REF] in Oklahoma support the hypothesis that a change in seismicity rate matches wastewater injection timing. Changepoint analysis (as other statistical tools) is not designed to test for causality. The only way to address causality is through a controlled experiment in which the observed response to a single imposed change in conditions is determined. Thus, for instance, from 1969 to 1973, variations in seismicity had been produced by controlled variations in the fluid pressure in Rangely, Colorado [START_REF] Raleigh | An experiment in earthquake control at Rangely, Colorado[END_REF]. Changepoint analysis is an interesting tool to analyze earthquake time-series analysis, because component simplification makes data more easily interpretable [START_REF] Lykou | Change point analysis on the Corinth Gulf (Greece) seismicity[END_REF]. It is under this aspect that we envisage the use of changepoint analysis in this study.

A profusion of methods have been introduced to detect the time points of possible existing changes in time series (see [START_REF] Horvath | Extensions of some classical methods in change point analysis[END_REF]. Some have been applied to earthquake sequences [START_REF] Pievatolo | Analysing the interevent time distribution to identify seismicity phases: A Bayesian nonparametric approach to the multiple-changepoint problem[END_REF][START_REF] Gupta | A Bayesian change point model to detect changes in event occurrence rates, with application to induced seismicity[END_REF][START_REF] Montoya-Noguera | Bayesian identification of multiple seismic change points and varying seismic rates caused by induced seismicity[END_REF][START_REF] Fiedler | Multiple change-point detection in spatiotemporal seismicity data[END_REF][START_REF] Lykou | Change point analysis on the Corinth Gulf (Greece) seismicity[END_REF]. For most of these authors, the approaches they use require that the event occurrences follow a Poisson distribution. Thus, the analyses carried out by [START_REF] Gupta | A Bayesian change point model to detect changes in event occurrence rates, with application to induced seismicity[END_REF] and [START_REF] Fiedler | Multiple change-point detection in spatiotemporal seismicity data[END_REF] on induced seismicity series are based on the method of [START_REF] Raftery | Bayesian analysis of a Poisson process with a change-point[END_REF], which is specifically designed for Poisson series. [START_REF] Raftery | Bayesian analysis of a Poisson process with a change-point[END_REF] assume that a change model is better than the constant rate of occurrences in model (i.e., a breakpoint is detected) when the Bayes factor falls below a given limit. The Bayes factor is the ratio of the likelihood of the constant rate model (hypothesis H1) to the changepoint one (hypothesis H2):

E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 1 ; 5 3 ; 4 3 2 B 12 prD | H 1 prD | H 2 -1 ; 1
in which prDjH k is the probability density of data D assumed to have arisen under H k . This density is an integral:

E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 2 ; 5 3 ; 3 6 7

prD | H k Z prD | θ k ; H k πθ k | H k dθ k ; 2
in which θ k is the parameter under H k , πθ k jH k is its prior density, and prDjθ k ; H k is the likelihood function of θ. For a narrow class of models (such as exponential distributions), the exact (and then accurate) evaluation of integral (equation 2) is possible. Both in [START_REF] Gupta | A Bayesian change point model to detect changes in event occurrence rates, with application to induced seismicity[END_REF] and [START_REF] Fiedler | Multiple change-point detection in spatiotemporal seismicity data[END_REF] studies, the assessment of equation ( 2) is tuned for a Poisson distribution of earthquakes. Actually, the method of [START_REF] Raftery | Bayesian analysis of a Poisson process with a change-point[END_REF] assumes that the investigated series can be generated by two independent Poisson processes with each one a different rate value. This assumption, however, does not hold for earthquakes where clustering of interevent times is the rule (as suggested for instance by Omori's law [START_REF] Omori | On aftershocks of earthquakes[END_REF], for a review De Arcangelis et al., 2016]). In this study, use of resistant, robust, and nonparametric statistical techniques has been applied successfully for the analysis of climate or seismological data [START_REF] Lanzante | Resistant, robust and non-parametric techniques for the analysis of climate data: Theory and examples, including applications to historical radiosonde station data[END_REF][START_REF] Lanzante | Temporal homogenization of monthly radiosonde temperature data. Part I: Methodology[END_REF][START_REF] Amorese | Applying a change-point detection method on frequency-magnitude distributions[END_REF][START_REF] Amorese | Changepoint analysis of geophysical time-series: Application to landslide displacement rate (Séchilienne rock avalanche, France)[END_REF]. This approach, unlike Raftery and Akman's method, does not require the data to obey any Poisson distribution. It follows the classic Wilcoxon-Mann-Whitney (WMW) nonparametric test, also referred to as Wilcoxon rank sum test [START_REF] Wilcoxon | Individual comparisons by ranking methods[END_REF] or the Mann-Whitney U test [START_REF] Mann | On a test of whether one of two random variables is stochastically larger than the other[END_REF]. Because this changepoint extraction technique uses rank sums, it is named as rank-sum multiple changepoint method (RSMCPM; [START_REF] Amorese | Changepoint analysis of geophysical time-series: Application to landslide displacement rate (Séchilienne rock avalanche, France)[END_REF]. This iterative method, fully described in [START_REF] Amorese | Changepoint analysis of geophysical time-series: Application to landslide displacement rate (Séchilienne rock avalanche, France)[END_REF], is designed to search for multiple changepoints in an arbitrary time series. The RSMCPM has already proven that it can successfully compete with Bayesian methods [START_REF] Amorese | Changepoint analysis of geophysical time-series: Application to landslide displacement rate (Séchilienne rock avalanche, France)[END_REF]. The iterative process of the RSMCPM can be described as follows: At each point i in the series of n points, the sum of the ranks (SRi) from the beginning of the series to that point is calculated. Because the sum of the ranks depends on the number of points, SRi is adjusted. Then, an adjusted value for the sum of the ranks, SAi, is as follows:

E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 3 ; 3 2 0 ; 5 6 2 SAi j2SRiin 1j: 3

In the right side of equation ( 3), in 1, is connected with After this, the WMW test is used to decide whether or not the null hypothesis (that there is no change in the sequence at n1) is rejected in favor of the alternative hypothesis. The chosen level of significance is 5% as the generally accepted and expected alpha (type I error rate) value in most disciplines for statistical tests. The RSMCPM is applied to a given series as long as the statistical significance of each new changepoint is less than the specified significance level. For each iteration, a list of N changepoints is delivered that defined N 1 segments. At each iteration, the series is adjusted by subtracting the median of its segment from each point.

E Q -T A R G E T ; t
In addition, a signal-to-noise ratio (SNR), which quantifies the magnitude of each discontinuity, is computed. For a given changepoint, this ratio appraises the variability associated with the shift in level between the adjacent segments relative to the variability within each segment. For more detailed calculations of this parameter, for the sake of brevity, we kindly refer readers to the article by [START_REF] Amorese | Changepoint analysis of geophysical time-series: Application to landslide displacement rate (Séchilienne rock avalanche, France)[END_REF].

The SNR can be used to discriminate between "weak" changepoints and "important" ones: "important" changepoints usually show SNR values of at least 0.05 or 0.1 [START_REF] Lanzante | Resistant, robust and non-parametric techniques for the analysis of climate data: Theory and examples, including applications to historical radiosonde station data[END_REF]. One must note that these SNRs differ from the ones computed from powers of signal and noise. The definition we used and applied is the one given by [START_REF] Lanzante | Resistant, robust and non-parametric techniques for the analysis of climate data: Theory and examples, including applications to historical radiosonde station data[END_REF], in which 0.05 or 0.1 values are definitively values for important breakpoints.

The RSMCPM offers the advantage of being relevant without concern for the shape (distribution) of the values in the populations providing the samples. Actually, the underlying assumption of the method is that the difference in the sums of ranks (adjusted for the sample size) mirrors the difference between statistical populations. Thus, the premise of homogeneity of variances (identically distributed data) is not required in this application of the WMW statistics. We believe that this approach is preferable, because it allows the analysis without the data first being distorted by any declustering operation. To complete the description of our calculations, let us add that, in this study, we used series differencing (the subtraction of two sequential values) to remove each series dependency on time structures such as seasonality or trends. This way, the breakpoint emergence is facilitated (e.g., [START_REF] Amorese | Changepoint analysis of geophysical time-series: Application to landslide displacement rate (Séchilienne rock avalanche, France)[END_REF].

The effectiveness of the method is illustrated in Figure A1 by its ability to produce results comparable to those of the multiple changepoint Bayesian detection method used by [START_REF] Fiedler | Multiple change-point detection in spatiotemporal seismicity data[END_REF] for the analysis of the Oklahoma induced seismicity. Fiedler et al.'s findings highlighted two changepoints in the declustered Oklahoma catalog: the first one in early 2009 and the second one in late 2013-early 2014. Applied to the raw (nondeclustered) catalog, the RSMCPM we use detects two discontinuities in 2008 and 2012, respectively (Fig. A1). We also resolve (2012, 2014, and 2016) changepoints in the injected volumes (Fig. A1).

RESULTS

The impact of wastewater disposals on Lacq seismicity through 0.7 km shallow (phases I and IV) and 4.8 km deep (phases I-III) injections was considered to be negligible, relative to the one of the gas extraction both from the involved volumes and from the pressure change approaches (e.g., [START_REF] Rothé | Seismes artificiels et exploitations petrolieres: l'exemple de Lacq (France)[END_REF][START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF][START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF][START_REF] Bardainne | Analysis of the induced seismicity of the Lacq gas field (Southwestern France) and model of deformation[END_REF].

The injection history we comprehensively report in this study allows to back-analyze the possible impact of Lacq fluid injections on seismicity in the context of the recent seismicity related to wastewater disposals. At the time of the first felt an earthquake on the Lacq field (M w 3, 24 November 1969), the injected volume in the shallow reef aquifer (through injection in the oil reservoir), is of the order of 7:5 × 10 6 m 3 (Fig. 6). This value is in the same magnitude range as the injected volumes reported to trigger M w 5 events in the central United States and Canada (e.g., for a review van der [START_REF] Van Der Elst | Enhanced remote earthquake triggering at fluid-injection sites in the midwestern United States[END_REF][START_REF] Van Der Elst | Induced earthquake magnitudes are as large as (statistically) expected[END_REF][START_REF] Mcgarr | Maximum magnitude earthquakes induced by fluid injection[END_REF][START_REF] Galis | Induced seismicity provides insight into why earthquake ruptures stop[END_REF]. This first local 1969 M w 3 event occurred after 12 yr of shallow injection. It is also contemporary to the peak value of injected volume rate (Fig. 5).

Since 1976, the phases (II-IV) of wastewater injection concentrate below the gas reservoir except for a short-lived shallow injection revival during 2007-2012, phase (III) (Figs. 2 and5). The shift from shallow to deep injection in 1976 precedes for a few years: (1) the sharp increase of seismicity rate and energy release for the Lacq seismicity (Fig. 5); (2) the onset of seismicity below the gas reservoir (e.g., [START_REF] Guyoton | Interrelation between induced seismic instabilities and complex geological structure[END_REF][START_REF] Bardainne | Etude de la sismicité de Lacq et analyse des formes d'ondes par décomposition en chirplets[END_REF][START_REF] Bardainne | Analysis of the induced seismicity of the Lacq gas field (Southwestern France) and model of deformation[END_REF]. The injection phase (III) increases the injection rates for both shallow and deep injection wells. It corresponds to a major increase in seismic energy release (2008 M w 3.0, 2013 M w 3.2, 2016 M w 3.9; Figs. 2 and5).

The 2008, 2013, and 2016 occurrences of these M w ≥ 3 events are concomitant with (1) a seismicity cluster a few kilometers from the deep injection well (Fig. 4) and ( 2) the resumption of shallow and deep injections during 2007-2012 phase (III). It points toward a cooperative effect of the long-lasting stress buildup in response to the gas reservoir depletion and the modification of the effective pressure due to the wastewater injections to drive the triggering of these late and relatively large earthquakes (e.g., Chang and Segall, 2016b).

The 2016 M w 3.9 earthquake is contemporary to a 2:44 × 10 7 m 3 of total injected volume (phases I-IV), including 0:85 × 10 7 m 3 for deep injection solely (Fig. 6). It matches the magnitude-volume pattern recently reviewed for seismicity driven by wastewater injections (e.g., McGarr 2014; [START_REF] Buijze | Moment partitioning for injection-induced seismicity: Case studies & insights from numerical modeling[END_REF][START_REF] Dieterich | Modeling injection-induced seismicity with the physics-based earthquake simulator RSQSim[END_REF].

To help identify interrelations between hydrocarbon recovery operations and seismicity, we further use a changepoint method (e.g., for Oklahoma case study, [START_REF] Gupta | A Bayesian change point model to detect changes in event occurrence rates, with application to induced seismicity[END_REF][START_REF] Fiedler | Multiple change-point detection in spatiotemporal seismicity data[END_REF]. The RSMCPM we use [START_REF] Amorese | Changepoint analysis of geophysical time-series: Application to landslide displacement rate (Séchilienne rock avalanche, France)[END_REF] does not identify any changepoint on the oil-gas production and gas depletion time series, for Lacq data. It supports relatively smooth gas reservoir operations over time. Changepoints for the differenced gas production are extracted in 1961, 1970, and 1977 (Table 1; Figs. 7 and8). Using the M w 2.5 as the completeness threshold for seismicity catalog, 1977, 1988, and 1997 changepoints emerge (Fig. 7a). The 1977 seismicity changepoint matches a main increase in seismicity rate (Figs. 2345). This seismicity changepoint is within a few years from changepoints on (1) the 1977 differenced gas production, (2) the 1974 deep injection volume (i.e., the onset of deep injection below the gas reservoir), and (3) the 1976 shallow injection volume, that is, end of injection phase (I) (Table 1; Figs. 7 and8). These results suggest the coupling between the increase of the gas production rate and the shift from shallow to deep wastewater injection (change from injection phase (I) to phase (II)), to control the significant 1977 change in seismicity rate we observe in the Lacq zone (Figs. 2 and5; Table 2).

DISCUSSION

Control parameters for the Lacq seismicity history [START_REF] Rothé | Seismes artificiels et exploitations petrolieres: l'exemple de Lacq (France)[END_REF] states that the comprehensive understanding of the Lacq seismicity is complex due to the double oil-gas extraction and the injections operations. When more seismicity were recorded over time, the oil field operations were rejected as driver of the seismicity because (1) the observed seismicity spans over a much larger area than the 6:8 km 2 of the oil reservoir, but overlaps with the 80 km 2 of gas field area; (2) most of the earthquakes occurred below the low-permeability marls that cap the gas reservoir and thus could not be connected to the shallow oil field and aquifer [START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF][START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF][START_REF] Bardainne | Analysis of the induced seismicity of the Lacq gas field (Southwestern France) and model of deformation[END_REF]. The only perturbation of the stress or strain fields that has the same spatial and temporal scale as the seismicity is the 30-60 MPa drop in gas pressure due to of gas production [START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF][START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF].

Previous studies point on that 80% of the seismicity in the 1975-1996 period is located in the area where the Coulomb stress changes, as estimated from a reservoir depletion model, are positive (e.g., [START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF]. This volume is still the one that bounds the 1969-2016 Lacq seismicity (Fig. 2), that is, within 1L r distance to the reservoir center. The poroelastic model we refer to is a mean field model that reproduces the global deformation pattern within a homogeneous elastic medium. It macroscopically fits the reported few centimeters subsidence [START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF], this latter being roughly proportional to the gas pressure drop [START_REF] Grasso | Mechanics of seismic instabilities induced by the recovery of hydrocarbons[END_REF][START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF]. On average, the stresschange model by [START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF] can be used as a proxy for the in situ stress change induced by the reservoir depletion. At a smaller scale, the Lacq field deformation is known (1) to be localized on discontinuities (e.g., [START_REF] Maury | Monitoring of subsidence and induced seismicity in the Lacq gas field (France): The consequences on gas production and field operation[END_REF][START_REF] Odonne | Abnormal reverse faulting above a depleting reservoir[END_REF];

(2) to correspond to both seismic and aseismic slip (e.g., on Lacq case study, the seismicity remains a few 1% of the observed surface displacement; [START_REF] Grasso | Seismicity induced by gas production: II. Lithology correlated events, induced stresses and deformation[END_REF][START_REF] Grasso | Mechanics of seismic instabilities induced by the recovery of hydrocarbons[END_REF].

In Lacq field, as in several gas fields worldwide [START_REF] Grasso | Mechanics of seismic instabilities induced by the recovery of hydrocarbons[END_REF][START_REF] Grasso | Testing self-organized criticality by induced seismicity[END_REF]NRC, 2013), the estimated maximum effective stress change that triggers the seismicity is of the order of 0.1 MPa. Specifically, the onset of the first local M w 3, 1969, in Lacq area coincides with (1) a 30 MPa gas reservoir depletion, that is, a 50% reservoir pressure drop and (2) the 3:5 MPa=yr peak value of the depletion rate (Fig. 2). For the Lacq gas reservoir, these values correspond to stress changes (as estimated through poroelastic model) in the seismic zone outside the reservoir rocks of the orders of 0.06 MPa and 0:01 MPa=yr, respectively (e.g., [START_REF] Segall | Earthquakes triggered by fluid extraction[END_REF][START_REF] Grasso | Mechanics of seismic instabilities induced by the recovery of hydrocarbons[END_REF][START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF]. These values are in the same order of magnitude as the estimated coseismic stress changes that are admitted to trigger aftershock sequences (e.g., [START_REF] King | Static stress changes and the triggering of earthquakes[END_REF][START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF][START_REF] Harris | Introduction to special section: Stress triggers, stress shadows, and implications for seismic hazard[END_REF]. The role of the rate of pressure variation (and production rate) to control the earthquake triggering is reported for other types of anthropogenic seismicity (e.g., for reservoir triggered seismicity, [START_REF] Simpson | Rate of change in Lake level and its impact on reservoir triggered seismicity[END_REF].

Contemporary to the poroelastic loading that is driven by the gas reservoir depletion, the onset of M w 3 events matches the climax value of the injected volume rate for wastewater in the shallow aquifer (Figs. 56789). These shallow injections increase the oil reservoir pressure by 0.6 MPa (Fig. 2). It is the same magnitude order as the effective 0.1 MPa poroelastic stressing outside the gas reservoir levels [START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF]. The shallow oil reservoir is known to connect to a vast water-saturated reef that tops the deep gas reservoir [START_REF] Grasso | Seismicity induced by gas production: II. Lithology correlated events, induced stresses and deformation[END_REF][START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF][START_REF] Maury | Monitoring of subsidence and induced seismicity in the Lacq gas field (France): The consequences on gas production and field operation[END_REF]. This connection suggests the wastewater injection to spread out within the rock matrices in between the 0.7 km depth of the oil reservoir and the 3.2 km roof of the gas reservoir (Fig. 4). This depth range is comparable to the depths where most seismicity was located during phase (I) of wastewater injection (e.g., [START_REF] Grasso | Ten years of seismic monitoring over a gas field[END_REF][START_REF] Guyoton | Interrelation between induced seismic instabilities and complex geological structure[END_REF]. Our analysis supports that three processes cumulate as trigger mechanisms for the Lacq M w 3 seismicity onset. The 1955-1969 gas reservoir depletion builds stresses in the rock matrices above and below the gas reservoir. This phase appears as a nonbrittle (aseismic) deformation, that is, no M w 2.5-3 earthquake recorded for shear stress increase of the order of 0.1 MPa. This gas reservoir depletion effect is re-enforced, above the gas reservoir, by phase (I) of wastewater injections in the shallow oil field (and the connected aquifer). The injection induces a 0.6 MPa increase of effective pore (Fig. 3). The peak values for both the depletion rate of the gas reservoir and the injection rate merge during the 1968-1969 period (Figs. 2 and5). Thus, the interrelation between injection and extraction is the most probable cause of the Lacq seismicity onset. The cumulative values for depletion and injection stress changes act as "static" trigger thresholds (e.g., [START_REF] Dieterich | Modeling injection-induced seismicity with the physics-based earthquake simulator RSQSim[END_REF]Chang and Segall, 2016b), with a second-order trigger as the loading rate peak values through injection and depletion rate, respectively.

The same mechanisms hold to understand the 1977 increase in seismicity rate. As detected by a changepoint method, the shift from the shallow phase (I) to deep phase (II) injections and the peak value for gas production are all within 1-2 yr before the 1977 seismicity increase (Table 2; Figs. 5, 7, and8).

As a final note, interactions between tectonics and anthropogenic seismicity may add complexity to the understanding of Lacq seismic swarm. The top largest M w 5.2, 5.0 Pyrenean earthquakes (1967 and[START_REF] Wittlinger | Etude de la sismicité en champ proche par un réseau sismologique à faible ouverture: Application au Frioul (Italie) et au gisement de Lacq[END_REF]Fig. 1b) lead by 2 yr (1) the 1969 onset of M w 3 and felt seismicity on Lacq field, and (2) the 1982 onset of deep seismicity below the gas reservoir. Grasso et al. (1992) suggested the observed time delay to be consistent with a viscoelastic model in which salt rock acts as a viscoelastic channel below an elastic bed (the regionally salt rock level being known as a "decollement" zone at the basement interface).

The late and largest 2016 M w 3.9 Lacq event: Implication for M max ?

There is a growing realization that the principal seismic hazard from earthquakes related to fluid manipulations comes from those associated with disposal of wastewater into deep strata or basement formations (NRC, 2013;[START_REF] Rubinstein | Myths and facts on wastewater injection, hydraulic fracturing, enhanced oil recovery, and induced seismicity[END_REF][START_REF] Eaton | What controls the maximum magnitude of injection-induced earthquakes?[END_REF]. During injection phases (III and IV), recurrent detections of seismicity from the operator seismic network point on most of the Lacq seismicity clusters in a 2 km zone in the southeast part of the Lacq field, in the immediate proximity of the deep injection well (Fig. 4b). The earthquake sequence (2008 M w 3.0, 2013 M w 3.2, and 2016 M w 3.9) is within a few years from the 2007-2012 injection, phase (III) in which simultaneous increases of deep and shallow injection volumes are reported (Fig. 5). We also point on the immediate proximity of the largest (2016 M w 3.9) event of the sequence to the deep injection well (Fig. 4).

Long-term and high-volume injection in deep wells clearly carries risk (e.g., [START_REF] Frohlich | Two-year survey comparing earthquake activity and injection-well locations in the Barnett Shale, Texas[END_REF][START_REF] Muntendam-Bos | A guideline for assessing seismic risk induced by gas extraction in the Netherlands[END_REF][START_REF] Rubinstein | Myths and facts on wastewater injection, hydraulic fracturing, enhanced oil recovery, and induced seismicity[END_REF][START_REF] Langenbruch | Value at induced risk: Injection-induced seismic risk from lowprobability, high-impact events[END_REF]. When most wastewater injection wells remain apparently aseismic worldwide (NRC, 2013), such a risk is magnified in the Lacq field context in which the gas reservoir depletion succeed to move the rockmass around the reservoir toward critical thresholds for failure [START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF][START_REF] Grasso | Testing self-organized criticality by induced seismicity[END_REF]. The cumulative effects of fluid extraction and injection apply to understand the late occurrence of the (Continued)

M w 3.9 event after the gas depletion stops (e.g., Chang and Segall, 2016b). Following the drastic reduction in gas production since 2013, larger events occurring post shut-in are thus not unexpected (e.g., [START_REF] Segall | Injection-induced seismicity: Poroelastic and earthquake nucleation effects[END_REF]. About 50 yr of gas reservoir depletion slowly increases the stresses above and below the reservoir, the same way the slow plates tectonic drives the natural seismicity. The Lacq anticline structure is now, on average, closer to instability than before the production started (e.g., [START_REF] Grasso | Testing self-organized criticality by induced seismicity[END_REF]Chang and Segall, 2016b).

To calibrate the possible impact of Lacq injection process on seismicity, we quantify that the injected volume and the recorded largest magnitude for Lacq field is in the same range (90% confidence level) than volume-magnitude pairs that are reported worldwide in a wide variety of tectonic settings (e.g., 2014; [START_REF] Galis | Induced seismicity provides insight into why earthquake ruptures stop[END_REF]. All these observations prompt to a re-evaluation of the local seismic hazard.

A second line of evidence that points toward a re-evaluation of the local seismic hazard emerges from the change in FMD for Lacq seismicity (Fig. 3). We observe that the break-in slope above M w 3 events for the FMD of Lacq seismicity in the 1L r distance range (valid up to 1996;[START_REF] Volant | The finite extension of fractal geometry and power law distribution of shallow earthquakes: A geomechanical effect[END_REF][START_REF] Lahaie | Loading rate impact on fracturing pattern: Lessons from hydrocarbon recovery, Lacq gas field, France[END_REF]) vanished on 1997-2016. This change of the FMD on is also evidenced by a decrease of b-value over time (Fig. 3). The reliable interpretation of Figure 3b for the 1969-1996 time period is delicate, because of the small number of points of the FMD that makes the b-value result unstable. There is, nevertheless, a noticeable difference with the b-value obtained for the [START_REF] Fourmaintraux | Utilisation de l'enregistrement sismique continu pour l'estimation du risque de sismicité associée aux réservoirs d'hydrocarbures et déclenchée par leur exploitation[END_REF]-2016time period, relatively to 1969-1996 (Fig. 3c) (Fig. 3c). We observed that the two patterns only exist when using the local seismicity within the 1L r box, not 10L r box size (Figs. 1 and3).

These statistical seismology patterns are contemporary to the deepening of seismicity. The depth of the recent significant shocks matches the basement rock matrices and faults (Fig. 4). During 2007-2012, injection phase (III), 70% of the seismicity is in between 3 and 5 km depth, the largest (M w 3) event on this period occurring at 7 km depth (CLSIC4000, 2013, see Data and Resources). Since 2013, phase (IV) most earthquakes cluster all around the deep injection well, 75% of the 2013-2016 events being at depth larger than 4 km (Fig. 4). Our analysis suggests that the depth of the Lacq triggered earthquakes increased with time to possibly reach the basement fault (Fig. 4). The three top largest events, 2008, 2013, and 2016 as M w 3, 3.2, and 3.9 have 7, 4.3, and 5.2 km depths, respectively. We surmise the activations of basement fault slips, below the gas reservoir, are possibly triggered by the coupled effect of poroelastic loading due to gas reservoir depletion and effective stress decrease driven by wastewater disposal (e.g., NRC, 2013;[START_REF] Dieterich | Modeling injection-induced seismicity with the physics-based earthquake simulator RSQSim[END_REF]Mutendam-Bos et al., 2015;Chang and Segall, 2016a;[START_REF] Maurer | Magnitudes of induced earthquakes in low-stress environments[END_REF]. This article's key target and result acknowledge the possible role of injection on seismicity. It is of key importance because the earthquakes related to wastewater disposal are reported to overpass the one related to pure extraction (e.g., [START_REF] Suckale | Induced seismicity in hydrocarbon fields[END_REF]NRC, 2013;[START_REF] Mcgarr | Maximum magnitude earthquakes induced by fluid injection[END_REF][START_REF] Foulger | Global review of human-induced earthquakes[END_REF].

From hazard assessment, this implies the possible earthquake on the Lacq area to be larger than the one expected from gas depletion processes solely. A larger and possibly shallowest earthquake relative to the tectonic earthquakes may induce local ground motion in the Lacq area that are not included when using tectonic earthquake rule (e.g., higher frequencies content in the near field evidenced by [START_REF] Fourmaintraux | Utilisation de l'enregistrement sismique continu pour l'estimation du risque de sismicité associée aux réservoirs d'hydrocarbures et déclenchée par leur exploitation[END_REF][START_REF] Klose | Human-triggered earthquakes and their impacts on human security[END_REF][START_REF] Bommer | Developing an application-specific ground-motion model for induced seismicity[END_REF][START_REF] Crowley | Probabilistic damage assessment of buildings due to induced seismicity[END_REF]. Next step studies need to improve the knowledge of the geological structure in the neighborhood of the injection wells, with a detailed knowledge of injection operations (e.g., [START_REF] Buttinelli | Inversion of inherited thrusts by wastewater injection induced seismicity at the Val d'Agri oilfield (Italy)[END_REF][START_REF] Vadacca | Slip tendency analysis, fault reactivation potential and induced seismicity in the Val d'Agri oilfield (Italy)[END_REF].

Conceptual model for Lacq seismic deformation

Following [START_REF] Volant | b-Value, aseismic deformation and brittle failure within an isolated geological object: Evidences from a dome structure loaded by fluid extraction[END_REF], we tentatively map the seismicity and deformation patterns we record on the Lacq field on a stress-strain path diagram that is used at the laboratory scale to characterize rock rheology (Fig. 9). As a first-order macroscopic diagram, the strain values for the Lacq anticline are directly estimated from the subsidence values as normalized by the reservoir depth [START_REF] Segall | Earthquakes triggered by fluid extraction[END_REF][START_REF] Grasso | Mechanics of seismic instabilities induced by the recovery of hydrocarbons[END_REF][START_REF] Maury | Monitoring of subsidence and induced seismicity in the Lacq gas field (France): The consequences on gas production and field operation[END_REF][START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF]. Surface subsidence is used as proxy for strain as subsidence over reservoir depth. Gas pressure drop works for stress outside the reservoir rocks as predicted by poroelasticity, Δσ α ΔP (e.g., [START_REF] Segall | Earthquakes triggered by fluid extraction[END_REF]; α is a function of reservoir rock properties of the order of 3 × 10 -2 for Lacq gas field (e.g., [START_REF] Grasso | Mechanics of seismic instabilities induced by the recovery of hydrocarbons[END_REF][START_REF] Grasso | Testing self-organized criticality by induced seismicity[END_REF]. The corresponding stress changes are controlled by the pore-pressure drop within the reservoir rock that transfers poroelastic shear stresses outside the reservoir rocks [START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF] and by the pore-pressure increase due to wastewater disposal. Both fluid perturbations correspond to effective stress changes outside the gas reservoir rocks that are of the same few 0.1 MPa magnitude order. In Figure 9, the first-order proportionality between the gas pressure drop and strain supports the poroelastic stress transfers that reproduce the gross deformation of the area [START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF]. On the 1976-1986 period, [START_REF] Grasso | Seismicity induced by gas production: II. Lithology correlated events, induced stresses and deformation[END_REF] estimated the brittle seismic deformation to account for only a few percent of the global subsidence. The M w 3.9 occurrence in 2016 does not change this ratio. Since 1969, the second-order deviation from the global elastic response may be tracked (Fig. 9). Such a departure from elasticity matches the onset of the Lacq seismic sequence. It corresponds to the onset of the brittle damage within a global elastic response (e.g., [START_REF] Grasso | Seismicity induced by gas production: II. Lithology correlated events, induced stresses and deformation[END_REF][START_REF] Grasso | Mechanics of seismic instabilities induced by the recovery of hydrocarbons[END_REF][START_REF] Volant | b-Value, aseismic deformation and brittle failure within an isolated geological object: Evidences from a dome structure loaded by fluid extraction[END_REF]. These patterns are similar to the ones observed from lab scale acoustic emission (AE) during compression tests on rock samples (e.g., for Lacq deep rock samples; [START_REF] Fabre | Mechanical behavior of deep rock core samples from a seismically active gas field[END_REF]. At the lab scale, such an increase in AE counts maps the divergence of brittle failure rate before macroscopic failure (e.g., [START_REF] Scholtz | The Mechanics of Earthquakes and Faulting[END_REF][START_REF] Amitrano | Brittle creep, damage, and time to failure in rocks[END_REF][START_REF] Davidsen | Scaling and universality in rock fracture[END_REF]. These paths toward failure are further enhanced when the coupling between poroelastic stress and fluid pressure increase exists (e.g., [START_REF] Amitrano | Fiber bundle model under fluid pressure[END_REF]Chang and Segall, 2016a). On such a basis, there is evidence for the Lacq seismicity to be neither deterministically predictable in duration nor in size for the M max value, that is, the maximum expected magnitude for macrofailure.

For now, we failed to identify relationships that would help anticipate the seismicity rate or the seismic energy release over time as a function of the local fluid manipulations, including the future development of this seismic swarm over space (Fig. 8). There is no French seismic regulation for oil and gas operations, except the use of a seismic network to monitor the earthquake activity. No traffic-light system either for the earthquake rate or for the magnitude is used. For Lacq case, the possible largest earthquake when related to waste disposal is larger than the one related to pure extraction operation (e.g., NRC, 2013;[START_REF] Foulger | Global review of human-induced earthquakes[END_REF]. It implies that the possible largest earthquake on the Lacq area to overpass an M w 4.5-5 event as previously expected from pure gas depletion hazard (e.g., [START_REF] Grasso | Mechanics of seismic instabilities induced by the recovery of hydrocarbons[END_REF][START_REF] Fourmaintraux | Utilisation de l'enregistrement sismique continu pour l'estimation du risque de sismicité associée aux réservoirs d'hydrocarbures et déclenchée par leur exploitation[END_REF].

Possible lessons for fluid-triggered seismic swarms

One important output of the analysis of the Lacq seismic swarm as a function of the fluid manipulations is its unexpected duration. When the sequence started in 1969, 15 yr of gas reservoir depletion and shallow wastewater injection had built up stresses that allow a long-lasting cascading process. These stresses build on through a slow anthropogenic loading rate. We estimate the stress rate of the order of 0:01 MPa=yr, as deduced from the peak value of the poroelastic stressing rate in 1969 (Fig. 2). An estimate for the strain rate derives from the few mm/yr value of the subsidence rate (e.g., [START_REF] Grasso | Seismicity induced by gas production: II. Lithology correlated events, induced stresses and deformation[END_REF][START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF]. This slow loading relatively to earthquake slip velocity favors the seismicity to self-organize as slip instabilities [START_REF] Grasso | Testing self-organized criticality by induced seismicity[END_REF][START_REF] Vespignani | How self-organized criticality works: A unified mean-field picture[END_REF]. Even when such a critical state, in term of phase transition, is reached, the stressed volume around the Lacq reservoirs remain insensitive to dynamic triggering by a nearby M w 5 earthquake (Fig. A2). The forecast of the Lacq swarm duration and how it depends on changes in fluid manipulations remain an open challenge. It questions the long-term pattern of seismicity around other depleted gas reservoirs (e.g., Groningen case study; [START_REF] Vlek | Induced earthquakes from long-term gas extraction in Groningen, the Netherlands: Statistical analysis and prognosis for acceptable-risk regulation[END_REF]. Even when we observe a short-term correlation of seismicity rate with gas production rate (e.g., for Lacq field, [START_REF] Lahaie | Loading rate impact on fracturing pattern: Lessons from hydrocarbon recovery, Lacq gas field, France[END_REF]Groningen field, Bourne and Oates, 2017;[START_REF] Bourne | The exponential rise of induced seismicity with increasing stress levels in the Groningen gas field and its implications for controlling seismic risk[END_REF][START_REF] Vlek | Induced earthquakes from long-term gas extraction in Groningen, the Netherlands: Statistical analysis and prognosis for acceptable-risk regulation[END_REF], our study points on the large uncertainties that rely on the occurrence and size of the maximum expected event magnitude (Fig. 8). Furthermore, here we evidence that the fluid injection cannot be rejected as a control parameter of the Lacq seismicity pattern. This result supports to not use the Lacq seismic swarm as a benchmark case for a pure depletioninduced seismicity.

Alternatively, the seismicity patterns we observe during the controlled fluid history around Lacq fields may help understand tectonic earthquake swarms in nonvolcanic context (e.g., [START_REF] Vidale | A survey of 71 earthquake bursts across southern California: Exploring the role of pore fluid pressure fluctuations and aseismic slip as drivers[END_REF][START_REF] Fischer | Intra-continental earthquake swarms in West-Bohemia and Vogtland: A review[END_REF][START_REF] De Barros | Fluid-induced swarms and coseismic stress transfer: A dual process highlighted in the aftershock sequence of the 7 April 2014 earthquake (M L 4.8[END_REF]. We learn from Lacq case study that leak off from overpressured lenses (as an analog, the Lacq gas reservoir onset of production) and fluid transfer can trigger seismicity with a decennial time lag (onset of Lacq seismicity 15 yr after production started). Multiple interactions between reservoir pressure changes and fluid flows (Fig. 8) possibly explain the complex pattern of either tectonic or anthropogenic seismic swarm, far beyond a simplistic hydraulic diffusion pattern (e.g., [START_REF] Vidale | A survey of 71 earthquake bursts across southern California: Exploring the role of pore fluid pressure fluctuations and aseismic slip as drivers[END_REF][START_REF] Fischer | Intra-continental earthquake swarms in West-Bohemia and Vogtland: A review[END_REF][START_REF] Schoenball | A systematic assessment of the spatiotemporal evolution of fault activation through induced seismicity in Oklahoma and Southern Kansas[END_REF][START_REF] Goebel | The spatial footprint of injection wells in a global compilation of induced earthquake sequences[END_REF][START_REF] De Barros | Fluid-induced swarms and coseismic stress transfer: A dual process highlighted in the aftershock sequence of the 7 April 2014 earthquake (M L 4.8[END_REF]. Similar to the aftershock zone size, our study supports that most of the triggered earthquakes remain in the near-field distance of the fluid reservoir. It may be used to size up the tectonic reservoir equivalent that drives seismic swarms in nonvolcano context worldwide.

CONCLUSION

In the 1969-2016 period, the Lacq seismicity has developed as a long-lasting seismic swarm. Revisiting the seismicity patterns and the fluid manipulations on Lacq site, we use seismicity rate and seismic energy for M w > 2:5 events. We propose that three processes may combine to trigger the Lacq seismicity onset. First, the 1955-1969 gas reservoir depletion builds stress in the rockmass above and below the gas reservoir. This phase appears as a nonbrittle (aseismic) deformation, that is, no M w > 2:5 earthquake recorded. In the same period, the stress change induced by the gas reservoir depletion is reinforced, above the gas reservoir, by 1955-1974 phase (I) injections in the shallow oil field. The 1969 onset of seismicity coincides with peak values of stress rate changes, as estimated by both the depletion rate of the gas reservoir and the shallow wastewater injection rate. The shift from shallow to deep injection in 1975 precedes by 2 yr the sharp increase in Lacq seismicity, these latter being also deeper than in the previous period. The interrelation between injection and extraction is the most probable cause of the Lacq seismicity onset and sustained activity. The same cumulative effects of fluid extraction and injection apply to trigger the late occurrence of the (2016 M w 3.9) largest (up to 2020) Lacq event. This earthquake occurred 5 yr after the largest injected volume rate on the past 30 yr. Our study points to the large uncertainties that rely on the occurrence and size of the maximum expected magnitude on the Lacq field. Primarily, our study points toward the underestimation of the possible impact of injection on the Lacq seismicity. In the context of reported magnitude-volume pairs for wastewater-induced seismicity, the amount of injected volumes on the Lacq site may substantiate triggering distances to the reservoir and maximum expected magnitude (M max ) to be larger than the ones estimated from the history of gas reservoir depletion solely.

Because of the possible role of injection on the current seismicity, we suggest that further studies need to improve the knowledge of the geological structure in the neighborhood of the injection wells, with a detailed knowledge of injection operations in parallel with warning system criterion to calibrate the occurrence of large events. . Lacq reservoir seismicity response to exogenous mainshock (M w 5 regional tectonic seismicity). Stars are M w 5.2 and 5.0 earthquakes that occurred in 1967 and 1980, respectively. Squares and circles are M w 2.5 aftershocks within 60 days from the mainshock for 1967 and 1980 events, respectively. Crosses are M w 2.5 events within 60 days before the mainshock for 1967 and 1980 events, respectively. Note that there is no triggering within the Lacq area after the two M w 5 Pyrenean earthquakes. Light gray squares and circles are earthquakes within 3L distance from each mainshock. Lacq reservoir contour line is gray (red). Boxes around M w 5 mainshocks are 1L and 3L distances from mainshock using M f L [START_REF] Wells | New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement[END_REF]. The color version of this figure is available only in the electronic edition.

Figure 1 .

 1 Figure 1. Regional and local seismicity map in Lacq area. (a) Pyrenees seismicity, 1989-2011. The western swarm located north of the Pyrenees mountain range overlaps with the Lacq gas field extension. Blue square is the Lacq swarm as seen in (b,c). Black lines are major faults, as NPFT, North Pyrenean frontal thrust; SPFT, South Pyrenean frontal thrust. In red, the current active faults as suggested by Lacan and Ortuño (2012) (modified from Rigo et al., 2015). (b) Regional seismicity map around the Lacq gas field, 1962-2016, M w 2.5. Box sizes represents 1L r , 3L r , and 10L r distances to the gas reservoir center, L r 10 km being the reservoir width. Largest

Figure 2 .

 2 Figure2. Yearly values for oil-gas field pressures and seismicity rate energy within 1L r distance from the gas reservoir. (a) Pressure history for oil and gas fields and seismicity rate energy. (b) Yearly rate for pressure changes (i.e., derivative of (a)) and seismicity energy. Seismicity and seismic energy (gray bars) is M w ≥ 2:5 within 1L r distance from the gas reservoir. Peak values are 3.5 and 0:3 MPa=yr for yearly gas and oil reservoir pressures, respectively. Note that a 64.5 MPa gas reservoir depletion value induces a maximum 0.2 MPa poroelastic stress increase in the rock matrices outside the gas reservoir (e.g.,[START_REF] Segall | Poroelastic stressing and induced seismicity near the Lacq gas field, southwestern France[END_REF]. The 0.6-1 MPa pore-pressure increase in the shallow oil field corresponds to a direct effective pore-pressure change. The color version of this figure is available only in the electronic edition.

Figure 3 .

 3 Figure 3. (a) Magnitude (M w ) values as a function of time. Black color is the near field seismicity (within 1L r distance from reservoir); gray color, far-field seismicity (10L r distance from reservoir). L r is the gas reservoir dimension. (b-e) Frequency-magnitude distributions. Open circle, cumulative distributions; triangle, discrete distributions. (b,c) Frequency-magnitude distribution for 1L r and (d,e) 10L r distances from the Lacq gas reservoir, respectively. (b,d) 1969-1996; (c,e) 1997-2016. The seismicity catalogs we used are open-access data at IS-EPOS (2018, see Data and Resources).

Figure 4 .

 4 Figure 4. (a) Schematic cross section of the Lacq anticline structure; (dark [green] zone), the shallow oil reservoir; (deep dark [red] zone), the gas reservoir. Thick vertical (purple) bars are the injection well locations (modified from Lahaie and Grasso, 1999); the blue star is the location of the largest Lacq event, 2016 M w 3.9; vertical axis is depth in km; (b) seismicity map, 2013-2016; square is the deep injection well location. Triangle is the seismic antenna location. Broken line is the gas reservoir contour as Figure 1b,c for comparison; (c) hypocenter depth histogram (2013-2016); (b,c) -0:4 ≤ M w ≤ 4:0, data from CLSIC4000 (2019, see Data and Resources). Note the local network estimates 2016 event as M w 4.0 while the regional catalog estimate is M w 3.9. The color version of this figure is available only in the electronic edition.

Figure 5 .

 5 Figure 5. Yearly rate for seismicity, gas production, shallow, and deep injection volumes. (a) Normalized values for seismic energy rate and fluid manipulations; (b) same as (a) for absolute values of injected volume and seismicity rate. Peak values are 11 × 10 9 m 3 =yr for gas extraction, 9:1 × 10 5 m 3 =yr and 3:5 × 10 5 m 3 =yr for shallow and deep injected volume, respectively. Seismicity is for M w ≥ 2:5 event. The color version of this figure is available only in the electronic edition.

  e m p : i n t r a l i n k -; d f 4 ; 3 2 0 ; 4 9 7EWi in 1=2;4 which is the expected value of the rank sum for the i first observed ranks out of a total of n points. The next step of the procedure is to find the maximum of SAi to divide the series into two segments. The point n1 being where the value of SAi is maximum, the following variables are defined:E Q -T A R G E T ; t e m p : i n t r a l i n k -; d f 5 T A R G E T ; t e m p : i n t r a l i n k -; d f 6 ; 3 2 0 ; 3 4 0 n2 n -n1: 6

Figure 6 .

 6 Figure 6. The maximum observed magnitude and the injected volume pairs on of 50 case studies, worldwide. Red line is the regression line through the wastewater injection-earthquake pairs using data from McGarr (2014) and Buijze et al. (2015). Blue line is the regression line through the hydrofrack-earthquake pairs using data from Atkinson et al. (2016) and Maxwell (2013). The transparency plot around each regression line shows the 90% prediction band. Dotted and black circles are the Lacq volume-magnitude pairs as 1969 M w 3.0 first felt earthquake and the maximum observed magnitudes 2016 M w 3.9, respectively. Injected volume values we used are 7:5 × 10 6 m 3 for the total injected volume at the time of the 1969 event (all injections are shallow injections) and 0:85 × 10 7 m 3 for the total of deep injected volume at the time of the 2016 deep event. The Lacq 2016 magnitude-volume pair is accepted to fit the worldwide pattern for wastewater injection seismicity at a 90% confidence level. The color version of this figure is available only in the electronic edition.

Figure 7 .

 7 Figure 7. Pattern recognition for changepoints in fluid manipulations and seismicity yearly time series. (a) Seismicity rate; (b) shallow injection rate; (c) deep injection rate; (d) differenced gas production; and (e) synthesis of changepoint detection on seismic data as plotted on rough gas production and gas reservoir pressure drop. Vertical (red) lines are seismicity changepoints from (a); the (red) dotted line is the Lacq seismicity onset as defined by the first 1969 M w 3. For (a-d), each insert reports the changepoint time, the signal-to-noise ratio (SNR), and the corresponding p-value for confidence level (p-value for the Wilcoxon-Mann-Whitney test, e.g., Amorese et al., 2018). The color version of this figure is available only in the electronic edition.

Figure 8 .

 8 Figure 8. Time series for seismicity and fluid manipulations. (a) Cumulative M w 2.5 seismicity, seismic energy, gas production, and wastewater injection, reservoir pressures for Lacq fields. Cumulative values are 254:25 × 10 9 m 3 for gas production, 244 × 10 5 m 3 for injected wastewater, 64.5 MPa for gas reservoir pressure drop, and 1 MPa for oil reservoir pressure increase. (b) Differentiated time series for yearly rates of fluid manipulations and seismicity: shallow and deep wastewater injections, oil and gas pressure, and seismic energy and event count (M c 2:5). Solid orange vertical lines are the seismicity changepoints as reported in Figure 7a. The dashed orange vertical line marks the 1969 onset of the Lacq swarm. The color version of this figure is available only in the electronic edition.(Continued)

Figure 8 .

 8 Figure 8. Continued

Figure 9 .

 9 Figure 9. Conceptual stress-strain diagram for brittle and total deformation on the Lacq gas field. Step (gray) curve is the cumulative number of earthquakes. Diagonal (blue) line is the gas pressure drop as a function of surface subsidence on Lacq field, as normalized to convert to strain and stress outside the reservoir rocks (see Conceptual Model for Lacq Seismic Deformation section for details). Diamond denotes the 2016 seismicity level. The color version of this figure is available only in the electronic edition.

Figure A2

  Figure A2. Lacq reservoir seismicity response to exogenous mainshock (M w 5 regional tectonic seismicity). Stars are M w 5.2 and 5.0 earthquakes that occurred in 1967 and 1980, respectively. Squares and circles are M w 2.5 aftershocks within 60 days from the mainshock for 1967 and 1980 events, respectively. Crosses are M w 2.5 events within 60 days before the mainshock for 1967 and 1980 events, respectively. Note that there is no triggering within the Lacq area after the two M w 5 Pyrenean earthquakes. Light gray squares and circles are earthquakes within 3L distance from each mainshock. Lacq reservoir contour line is gray (red). Boxes around M w 5 mainshocks are 1L and 3L distances from mainshock using M f L[START_REF] Wells | New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement[END_REF]. The color version of this figure is available only in the electronic edition.

TABLE 1

 1

	Fluid Manipulation Phases, 1955-2016, Lacq Field
	Lacq Fluid Manipulations (1955-2016)
	Fluid Manipulation Phase	Extraction Operation (Depth, Type)	Injection Operation (Depth, Type)
	Phase I (1955-1974)	3-5 km, Gas reservoir	0.7 km, Wastewater
	Phase II (1974-2006)	3-5 km, Gas reservoir	4.5 km, Wastewater
	Phase III (2006-2012)	3-5 km, Gas reservoir	0.7 and 4.5 km, Wastewater
	Phase IV (2013-2016)	Negligible	4.5 km, Wastewater

1955 is the onset of the gas extraction.

Volume XX Number XX -2021 www.bssaonline.org

Bulletin of the Seismological Society of America • 7

TABLE 2

 2 Changepoint Detection for Seismicity and Fluid Manipulation History

		Timing for Changepoint Analysis (RSMCPM)
	Seismicity		Gas Production	Injection	Injection
	(Observation)	Seismicity	(Differenced)	(Shallow)	(Deep)	Injection Phase

• Bulletin of the Seismological Society of America www.bssaonline.org Volume XX Number XX -2021 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

Volume XX Number XX -2021 www.bssaonline.org Bulletin of the Seismological Society of America • 3 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

• Bulletin of the Seismological Society of America www.bssaonline.org Volume XX Number XX -2021 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

Volume XX Number XX -2021 www.bssaonline.org Bulletin of the Seismological Society of America • 5 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

• Bulletin of the Seismological Society of America www.bssaonline.org Volume XX Number XX -2021 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

• Bulletin of the Seismological Society of America www.bssaonline.org Volume XX Number XX -2021 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

Volume XX Number XX -2021 www.bssaonline.org Bulletin of the Seismological Society of America • 9 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

• Bulletin of the Seismological Society of America www.bssaonline.org Volume XX Number XX -2021 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

Volume XX Number XX -2021 www.bssaonline.org Bulletin of the Seismological Society of America • 11 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

• Bulletin of the Seismological Society of America www.bssaonline.org Volume XX Number XX -2021 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

• Bulletin of the Seismological Society of America www.bssaonline.org Volume XX Number XX -2021 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

Volume XX Number XX -2021 www.bssaonline.org Bulletin of the Seismological Society of America • 15 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

• Bulletin of the Seismological Society of America www.bssaonline.org Volume XX Number XX -2021 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

Volume XX Number XX -2021 www.bssaonline.org Bulletin of the Seismological Society of America • 17 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

• Bulletin of the Seismological Society of America www.bssaonline.org Volume XX Number XX -2021 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by

Volume XX Number XX -2021 www.bssaonline.org Bulletin of the Seismological Society of America • 19 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

• Bulletin of the Seismological Society of America www.bssaonline.org Volume XX Number XX -2021 Downloaded from http://pubs.geoscienceworld.org/ssa/bssa/article-pdf/doi/10.1785/0120200359/5310981/bssa-2020359.1.pdf by 19225

ACKNOWLEDGMENTS

Abror Karimov is supported by EPOS-IP Anthropogenic Hazard Grant (EC H2020-EPOS IP: 676564). Jean-Robert Grasso's contribution is partially supported by EPOS-IP (Anthropogenic Hazard) and Fonds européen de développement régional (FEDER: SISM@LP-Swarm research project, POAI-PA0014885). The authors thank M. Galis, J. Letort, C. Martin, S. Maxwell, and A. Rigo, for sharing data set and information. Comments and suggestions from Associate Editor, D. Fitzenz, Editor-in-Chief, T. Pratt, and anonymous reviewers greatly improve previous versions of this article. Discussion with partners of the EPOS-Anthropogenic Hazard group motivated the study and strengthened the contents.

DATA AND RESOURCES

All the fluid manipulation data we report and the seismicity catalogs we used are open access data at Induced Seismicicity -European Plate Observing System (IS-EPOS) (2018), from Implementation Phase -European Plate Observing System (IP-EPOS) Platform as https://tcs .ah-epos.eu/#episode:LGF. Data should be cited as IS EPOS (2018), Episode: LACQ GAS FIELD, available at https://tcs.ah-epos.eu/ #episode:LGF (last accessed December 2019). Report on local seismicity is available at CLSIC4000, 2013, 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). Data from local "industrial" network, 2013-2016 are available at CLSIC4000, 2019, Comité local du suivi des injections Cretacé 4000, rapport quinquenal, 2013-2018, Mai 2019, 51 pp., 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). The rank-sum multiple change-point method (RSMCPM) algorithm and software for changepoint detection are fully available in the ACA R package (ACA: Abrupt Change-Point or Aberration Detection in Point Series), R package version 1.1, available at https://CRAN.R-project .org/package=ACA (last accessed September 2020).

DECLARATION OF COMPETING INTERESTS

The authors acknowledge that there are no conflicts of interest recorded.

APPENDIX

Figure A1 is the changepoints we identify when testing our method on the same seismicity data (Oklahoma, 1980(Oklahoma, -2018))


David Bourrier
email: bourrier@laas.fr

Ayako Mizushima

Naonobu Shimamoto

Etsuko Ota

Hugues Granier

Yoshio Mita

Clean

Clean room engineers' exchanges: a win/win international collaboration

à 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.

Clean room engineers' exchanges: a win/win international collaboration

David Bourrier 1* , Ayako Mizushima 2 , Naonobu Shimamoto 2 , Etsuko Ota 2 , Hugues Granier 1 and Yoshio Mita 2

1 LAAS-CNRS, Université de Toulouse, CNRS, Toulouse, France 2 The University of Tokyo; 7-3-1 Hongo, Bunkyo-city, Tokyo, Japan

Abstract:

In the world of research, international collaborations between researchers are common, well-structured, and benefit from a lot of funding. Exchanges between engineers/technical staff in general are rarer, it is difficult to find dedicated funding Through the opportunity offered to two engineers of LAAS-CNRS by Professor Mita of the University of Tokyo, supported by MEXT X-nics UTokyo Agile-X project, we will show that these exchanges are an opportunity for real scientific and technical collaborations; and that as such they deserve to be developed. The exchange focused on two aspects. The training and transfer of skills by the LAAS-CNRS in the field of electrolytic deposition, on equivalent equipment of the Yamomoto-MS company for the engineers of the Mita Lab, d.lab Platform Device Research Division. The discovery and the learning from the Mita Lab of a new technology of conformal seed layer deposition by supercritical CO 2 [START_REF] Usamil | IEEE 32 nd International Conference on Microelectronic Test Structures (ICMTS)[END_REF] deposition to study its feasibility and its transfer to the LAAS. We will present the very good preliminary results obtained in only 15 days as well as all the benefits brought by this exchange; as well on the technological, scientific and humanly. They will benefit the whole community of the two laboratories and will allow to accelerate new scientific results. These international exchanges allow us to confront and enrich different and complementary visions of how to implement technologies. Which is all the more enriching and motivating for the evolution of ways of working and exchanging. We hope that our presentation will create other opportunities, and show the real interest for our supervisory bodies to organize this type of opportunities beneficial to the whole scientific and technological community.

  a. 0.5 µm of copper electroplated on chip b. wafer Si 3 N 4 80 nm + Au 100nm c. Japanese-French electroplating Team

This work was supported by MEXT Initiative to Establish Next-generation Novel Integrated Circuits Centers (X-NICS) Grant Number JPJ011438.

This work was supported by LAAS-CNRS micro and nanotechnologies platform, a member of the Renatech french national network