On the Source Parameters and Genesis of the 2017, Mw 4 Montesano Earthquake in the Outer Border of the Val d’Agri Oilfield (Italy)

On October 27, 2017, an Mw 4 earthquake occurred close to the municipality of Montesano sulla Marcellana, less than 10 km external to the concession of the largest European onshore hydrocarbon reservoir—the Val d’Agri oilfield (Southern Italy). Being a weak event located outside the extended monitoring domain of the industrial concession, the relevance of this earthquake and the possible links with the hydrocarbon exploitation were not extensively discussed. Actually, the analysis of shallow seismic events close to subsurface exploitation domains plays a significant role in the definition of key parameters in order to discriminate between natural, triggered, and induced seismicity, especially in tectonically active regions. The study of weak-to-moderate earthquakes can improve the characterization of the potentially destructive seismic hazard of this particular area, already struck by M > 6.5 episodes in the past. In this work, we analyze the source parameters of this Mw 4 earthquake by applying advanced seismological techniques to estimate the uncertainties derived from the moment tensor inversion and identify plausible directivity effects. The moment tensor is dominated by a NW–SE oriented normal faulting with a centroid depth of 14 km. A single ML 2.1 aftershock was recorded and used as the empirical Green’s function to calculate the apparent source time function for the mainshock. Apparent durations (in the range 0.11–0.21 s, obtained from S-waves) define an azimuthal pattern, which reveals an asymmetric bilateral rupture with 70% of the rupture propagation in the N310°W direction, suggesting a rupture plane dipping to the SW. Our results tally with the activation of a deeper fault segment associated with the Eastern Agri Fault System close to the basement as the origin of the Montesano earthquake. Finally, the Coulomb stress rate induced by depletion of the oilfield is calculated to quantify the trigger potential estimated for the Montesano earthquake yielding relatively low probabilities below 10%. Our analyses point toward the conclusion that the Mw 4 event was more likely due to the local natural tectonic stress, rather than induced or triggered by the long-term hydrocarbon extraction in the Val d’Agri oilfield.


INTRODUCTION
Induced seismicity can nowadays be considered as a relevant and pressing problem of increasing interest for the general public. Underground operations such as shale gas and oil exploitation, mining, and other energy technologies could generate seismic activity (McGarr et al., 2002;Grigoli et al., 2017;Braun et al., 2018;Foulger et al., 2018). Some directly induced earthquakes reached moderate moment magnitude (e.g., M w > 5) and were strong enough to cause damages, and even casualties (Grigoli et al., 2018). In some cases, moderate-to-strong earthquakes may also be triggered at distances larger than 10 km from the underground operations (Keranen et al., 2014;Goebel et al., 2017), with an impact on the regional seismic hazard related to such activities. Usually, a well-constrained hypocentral location with small uncertainties as well as the temporal correlation between geotechnical activities and seismicity patterns play an important role in distinguishing the origin of seismic events (Healy et al., 1968;Raleigh et al., 1976). However, a detailed and robust analysis about the seismic source of potentially induced events could reveal and clarify the discrimination between induced, triggered, or natural earthquakes. Generally, the term "induced seismicity" is referred to anthropogenic seismic events in a wide sense encompassing both pure induced and triggered seismicity; however, different specific definitions have been given for these two seismic processes. According to McGarr and Simpson (1997), induced earthquakes are events where most of the stress changes release during the rupture produced by human action, while triggered events release a substantial amount of tectonic stress. More precisely, according to Dahm et al. (2015), an induced event is when the rupture is driven by the induced stress over the full rupture plane, whereas for triggered events only the rupture initiation is caused by the stress rate at the hypocenter. In other words, induced seismic events are entirely controlled by stress changes caused by human operations and would not have occurred without anthropic intervention; on the other hand, in triggered seismicity, the tectonic stress plays a primary role, with a stress variation on favorably oriented faults in agreement with the existing regional or local background stress field and geological structure.
In the particular case of Italy, a framework for the management of the geophysical monitoring associated with hydrocarbon reservoir exploitation and gas storage was established in 2014 by the publication of the "Guidelines for Monitoring Seismicity, Ground Deformation and Pore Pressure in Subsurface Industrial Activities" (ILG, Dialuce et al., 2014). The ILGs are currently experimentally applied for selected pilot sites , including the largest European onshore hydrocarbon reservoir, the Val d'Agri (VA) oilfield (Southern Italy), which will be the focus of our study (Figure 1). With the purpose of distinguishing between natural seismicity from possibly induced or triggered one, the ILGs define limits, domains, and technical requirements intended to localize the seismicity in a volume surrounding the area where human activities take place. More specifically, the ILGs define: (i) the Internal Monitoring Domain (DI) as the volume that includes the mineralized zone (oilfield) extended by an additional 5 km wide volume and (ii) the Extended Monitoring Domain (DE) beyond the DI to a neighborhood of 5-10 km ( Figure 1b). Moreover, the ILGs also recommend the application of a Traffic Light System (TLS) protocol as guidance for the seismic surveillance management (Bommer et al., 2006;Bosman et al., 2016), which is based on the estimation of seismic parameters (hypocenter, magnitude, peak ground velocity, and peak ground acceleration) for events occurring within the DI close to operational reinjection wells (Dialuce et al., 2014;Braun et al., 2020).
The present-time tectonics of the VA is mainly characterized by NE-SW extension, which is accommodated by two main quaternary fault structures, the Monti della Maddalena (MMFS) and the Eastern Agri Fault System (EAFS), bounding the VA on the SW edge and the NE edge, respectively (Maschio et al., 2005). This active deformation is expressed seismically through mostly high-angle, normal-faulting earthquakes that occur on the two main NW-SE-trending fault systems. The occurrence of sparse background seismicity is well known, moderate seismic events in VA are significant, and historical seismicity suggests that large and destructive events can occur ( Figure 1A). The most destructive earthquake, with an equivalent magnitude M w 7.0, occurred on December 16th, 1857, causing extreme damage (I max XI, EMS-98) and accounting for more than 11,000 victims (Guidoboni et al., 2019;Rovida et al., 2019;Rovida et al., 2020). A conclusive characterization of the event source is still under debate: Benedetti et al. (1998), Cello et al. (2003, and Barchi et al. (2006) ascribe the main 1857 M7 event to the EAFS; more recently, Burrato and Valensise (2008) indicate a ∼40 km-long complex rupture of at least two segments of the MMFS ( Figure 1A). Comprehensive and historical reports about the effects of the 1857 M7 earthquake were performed by Mallet (1862) and later re-edited by Guidoboni and Ferrari (1987). While the building damage and number of victims due to the 1857 M7 event were definitely higher in VA, Mallet (1862) reported significant surface effects for adjacent areas; for example, in the area of Montesano sulla Marcellana (Vallo del Diano), located in the western part of the MMFS, Mallet (1862) reports, "real surface ruptures in old rock layers for an extension of 200 feet". Even if our study does not investigate the hypocenter of the "great Neapolitan earthquake" (Mallet, 1862), it is however important to keep in mind that the exact location of the 1857 M7 event is still discussed.
According to the earthquake bulletins of ENI (the Italian governmental energy company assigned to the VA concession) and INGV (Istituto Nazionale di Geofisica e Vulcanologia), since 1990 in VA only few seismic events with a magnitude of M > 3 have been reported . Note that two significant clusters linked to anthropogenic operations have been well-described in previous works ( Figure 1A). The main swarm (depth range between 1 and 5 km and magnitudes smaller than Ml 3) occurred a few kilometers SW of the artificial Pertusillo lake in the southern termination of the MMFS and was interpreted as water reservoir-induced seismicity, due to an increased impoundment during the late winter/spring rainfalls (Valoroso et al., 2009;Valoroso et al., 2011;Stabile et al., 2014)  ; this is actually the only confirmed case of waste-water induced seismicity in Italy.
In this study, we investigate the source parameters and genesis of a weak-to-moderate earthquake that occurred at 22:38 UTC on October 26th, 2017, and was located 4 kmN of the municipality of Montesano sulla Marcellana (hereon referred to as Montesano earthquake). Its magnitude of M w 3.8 makes this event the largest earthquake in the Monti della Maddalena since the last 2 decades. Being located in the proximity of the external margin of the VA oilfield's extended monitoring domain, it is generally of interest to investigate whether the seismic source of this event may be natural, induced, or triggered. Due to the hypocentral depth of more than 10 km (10.9 and 14 km from ENI and INGV catalogues, respectively), the national monitoring service classified it as a natural tectonic event, excluding any anthropogenic origin. Taking advantage of a wide coverage of seismic stations deployed in the VA region, we describe a detailed and robust analysis of the source parameters for the Montesano earthquake in order to decipher its genesis and apply a probabilistic approach to address the question of discriminating between natural, triggered, and induced origin. A combination of advanced seismological techniques is applied to estimate the uncertainties derived from the moment tensor inversion and identify plausible directivity effects from Apparent Source Time Functions (ASTFs), revealing preferred rupture directions and potential asymmetries in the rupture propagation as well as to resolve the fault plane ambiguities. Finally, we analyze the origin of the source of the Montesano earthquake through a probabilistic approach based on the modeling of depletion-induced stress changes and the previously calculated seismological source parameters.

Geological Setting, Val d'Agri Oilfield and Seismic Monitoring
The formation of the hydrocarbon reservoir in VA is the result of a succession of complex multiple tectonic phases. Below the Quaternary deposits of the VA Basin, upper Messinian and/or Pliocene terrigenous marine deposits are found (Doglioni et al., 1999;Scrocca et al., 2005), underlain by 6-7 km thick allochthonous units of Mesozoic-early Tertiary carbonate sequence, which form the Apulian Platform (AP) and the Lagonegro units (LU) (e.g., Cello and Mazzoli, 1999). The collision of the African and the Eurasian Plates, which started in the late Cretaceous, caused a shortening in anti-Apenninic direction, which was accommodated by a series of thrust and reverse faults. This process led to a displacement of the LU and the overlying Apennine platform toward the east, overriding thus the AP (D' Argenio et al., 1975;Cello and Mazzoli, 1999;Mazzoli et al., 2001;Butler et al., 2004;Valoroso et al., 2009). The internal deformation of the AP by thrust faults and related folds, sealed on top by the marine deposits, formed the main structural traps for the oilfields in the region (e.g., Mazzoli et al., 2001;Butler et al., 2004;Shiner et al., 2004).
The giant deposit of VA, discovered in 1988, is actually the most productive oilfield onshore in Europe. The structural trap is represented by a large-scale pop-up bounded by high-angle reverse faults oriented in the NW-SE Apennine direction (Bertello et al., 2011) hosted in the extensional homonym Quaternary basin. The reservoir rocks for the oil/gas production are hosted inside the Cretaceous to Miocene limestone and dolostone of the Apulian Platform. The carbonate reservoirs include vuggy intervals, as well as extensively developed fracture systems, which provide flow rates in the range of 500-2000 Sm 3 /day (Standard cubic meters per day) from oil columns that vary in thickness from 600 to 1000 m. The entire producing complex is sealed by the overlying Neogene sequences (Wavrek and Mosca, 2004). Depending on the type of depositional facies, the hydrocarbon column of the VA reservoir shows significant lateral heterogeneities in density (in the range of 835-1013 kg/m 3 ), in petrophysical characteristics and in compositional grading (Wavrek and Mosca, 2004). The reservoir includes three main productive regions corresponding to the three geographical culminations named Cerro Falcone (CF), Monte Enoc (ME), and Monte Alpi (MA) ( Figure 1B). After a long initial test period, production started at MA in 1996. Subsequently, the number of production wells increased to the current number of 24 (ranging from a depth of 1.8-3.5 km below sea level), the crude oil treatment lines were improved, the CM2 wastewater disposal well came into operation, and the entire VA field achieved full production in 2010. The wastewater reinjection well CM2 was initially drilled for exploitation purposes in 1982-83; due to its limited production potential, it was converted into a reinjection well and became operational in 2006.
ENI has continuously operated a seismic network in the VA domain since 2001 (Eni Spa, 2001), when the first eight seismic stations were installed, providing a significant refinement in the location capacity and sensitivity of the national seismic monitoring system for the area. Schorlemmer et al. (2010) estimated that by 2008, in the broader VA region including Lake Pertusillo, the magnitude of completeness for the INGV network was between 1.5 ≤ M C ≤ 2.0. In 2012, the local ENI seismic network was upgraded to 15 stations . For the southern sector of the oilfield, this new configuration led to a lowering of the detection threshold of local seismic events to 1.0 ≤ M C ≤ 1.2 . In the framework of specific research projects that focused on the observation of the microseismicity potentially induced by the CM2-well and the Lake Pertusillo, between 2005 and 2012, temporary seismic stations were deployed in the VA area to complement the permanent ENI and INGV networks, thus providing a significant improvement for the integration of local seismic catalogs (Valoroso et al., 2009;Stabile et al., 2014).

MOMENT TENSOR INVERSION
We performed a deviatoric moment tensor inversion for the mainshock based on broadband seismological data using the probabilistic seismic source inversion algorithm Grond (Heimann et al., 2018;Kuḧn et al., 2020). Grond allows for simultaneously fitting the available seismic data in different frequency ranges, using different velocity models, and implementing different misfit definitions, both in the time and frequency domains. The inversion algorithm optimization starts with a random search of the parameter space and then progressively scans the parameter space closer to a range of best-fitting solutions. Parameter uncertainties are estimated by a bootstrap approach, simulating many different data configurations.
Here, we chose to simultaneously model full waveforms and amplitude spectra using data from the permanent INGV network (IV, INGV Seismological Data Centre 2006) and the local network (VA, Eni Spa, 2001). For near stations, located within 30 km from the epicenter, we build synthetic seismograms using the local 1D velocity model of Improta et al. (2017), while for stations at distances of 30-200 km we use a regional model from the Crust 2.0 database (Bassin et al., 2000). Both models are provided in the Supplementary Material S1. The inversion is setup to simultaneously fit different datasets, using in all cases 3-components data, as follows: 1) full waveform (60 s time windows) amplitude spectra in the frequency band 0.04-0.20 Hz for near stations, 2) full waveform (180 s time windows) amplitude spectra in the frequency band 0.04-0.08 Hz for far stations, 3) full waveform (60 s time windows) displacement traces in the frequency band 0.04-0.20 Hz for near stations, 4) cross-correlation of body wave displacement signals (from 2 s before P phases to 4 s after S phases) in the frequency band 0.5-2.0 Hz for near stations, and 5) cross-correlation of full waveform (180 s time windows) displacements in the frequency band 0.04-0.08 Hz for far stations. The choice of the frequency bands has been based on a number of considerations. At regional distances, modeling full waveforms can only be successfully achieved for relatively low frequencies, as high frequency waveforms cannot be reproduced at these distances with a simplified 1D model; below 0.04 Hz, however, the signal-to-noise ratio becomes too low, due to the moderate magnitude of the earthquake. At local distances, and when modeling body wave pulses, the seismic signals can be well reproduced to higher frequencies using the local velocity model. Results (Figure 2 and Table 1) show a moment tensor dominated by a NW-SE oriented normal faulting, plus an additional negative compensated linear vector dipole (CLVD) component, the presence of which is not robustly resolved (see uncertainties in Supplementary Figure S2). The resolved focal mechanism is in a good agreement with the INGV solution, but we estimate a slightly larger moment magnitude of M w 4.0 ± 0.2 and a slightly deeper depth of 14.0 ± 2.8 km (see Table 1).

ASTFs and Rupture Directivity Analysis
Earthquake sources can be isolated from seismograms applying empirical Green's functions (EGFs) techniques (Hartzell, 1978). For this purpose, a detailed knowledge of the earth structure between the target earthquake and receivers through which the waves propagated is required. Waveforms of fore-and aftershocks, sharing a common path with the mainshock, are then used as EGFs to accurately model these propagation effects.
A point source representing a delta function is assumed for the EGFs (typically, one to two magnitude units smaller than the mainshock), which must also show similar hypocentral locations and focal mechanism implying similar waveforms with the mainshock. Under these conditions, ASTFs can be retrieved through a deconvolution procedure revealing the relative moment rate functions of the target event observed at different receivers (e.g., López-Comino et al., 2016;Abercrombie et al., 2017;Stich et al., 2020). The shape and duration of these ASTFs can be then modeled providing constraints to identify preferred rupture directions and source complexities (e.g., López-Comino and Wu et al., 2019).
Following the M w four Montesano earthquake, a single M L 2.1 aftershock was recorded 13 days later at 09:51 UTC on November 8th, 2017, located less than 2 km from the mainshock at a depth of 15 km (Bolletino Sismico Italiano, INGV), with a similar hypocentral location (just 1 km deeper than the mainshock). Both events show similar polarities in the P-wave arrivals (Supplementary Figure S3) as well as similar waveforms ( Figure 3) revealing a common focal mechanism (Supplementary Figure S4); therefore, we use this aftershock as an EGF to calculate the ASTFs of the Montesano earthquake. We perform the deconvolution in the frequency domain through spectral division using the seismic recordings of the mainshock and the EGF (López-Comino et al., 2012). To avoid numerical instabilities derived from the deconvolution procedure, a Gaussian low pass with a pulse width of ∼0.075 s and a water level of 0.01 of the maximum spectral amplitude are used. Successful deconvolution results are obtained for S-wave windows (length of 6 s starting 1 s before the S-wave arrival) using the horizontal components ( Figure 3). The resulting ASTF at each station is obtained by stacking both components (east-west and north-south), reaching a good azimuthal coverage for 21 near-regional and local stations up to a distance of 41 km (Supplementary Figure S5). ASTFs are normalized to unit area according to the total seismic moment of the mainshock and negative values derived from the deconvolution procedure are removed. Apparent durations are identified manually by picking onset and end of each ASTF at the intersection of the initial and final slopes with the baseline assuming uncertainties of about 0.01s. One single and consistent pulse can be clearly identified for the Montesano earthquake with slightly shorter apparent durations at NW azimuths, and slightly larger durations for stations located at the opposite directions ( Figure 4). Apparent durations range from 0.11 up to 0.21 s.
Assuming a line source after Haskell (1964), theoretical predictions for unilateral and asymmetric bilateral ruptures are considered to adjust the azimuthal pattern of the apparent durations identified by the ASTFs (for technical details, see Cesca et al., 2011a andLópez-Comino et al., 2016). Some fixed parameters are considered to solve the nonlinear functions: local S-wave velocity of 3.2 km/s (from Improta et al., 2017) and a rise time of 0.09 s. Our results reveal an asymmetric bilateral rupture with 69 ± 8% of the rupture propagation toward N310°W ± 12°, rupture duration of 0.16 ± 0.03 s and rupture length of 0.29 ± 0.12 km ( Figure 4C and Table 2). This model yields lower L1 misfit of 0.01 s in comparison with pure unilateral rupture. The largest uncertainties are observed for the rupture velocity because intrinsic trade-offs among rise time, rupture length, and velocity remain in these inversions; and, we obtain a reasonable value of 2.9 km/s, corresponding to 90% of the local S-wave velocity. Resolving the fault plane ambiguities in a normal faulting focal mechanism is challenging due to the similar strike values in both planes and its uncertainties around 50°estimated from the moment tensor analysis ( Table 2). However, the directivity direction is well constrained from this analysis and small differences in the strike values are observed from the best solution of the DC moment tensor (strike1 131°(311°); strike2 327°), which could allow to identify the plane dipping to the SW as our preferred solution ( Figure 4C).

PROBABILISTIC DISCRIMINATION
The discrimination problem between a natural, triggered, or induced event is addressed through a probabilistic scheme to quantify the likelihood of a correlation with the depletioninduced stress perturbations (Dahm et al., 2015). This methodology based on physical-statistical seismicity models considers an earthquake occurring close to an oil-or gasfield, that has continuously produced over a period of several years and is depleted; this is the case of the Montesano earthquake. The result of this analysis is the probability that the target event has been triggered by the stressing rate induced by the depletion of the oilfield. The scheme can account for the target event source parameters and their uncertainties, which were addressed above.
In a first step, the steady-state tectonic stress rate is estimated from the background seismicity before the production started, using (Hainzl et al., 2010;Dahm et al., 2015): Variables in this equation are the a and b values of the Gutenberg Richter frequency-magnitude distribution, the maximum observed magnitude (M max ), the area of the seismogenic zone (A), and the seismogenic width (D). These parameters are considered for the seismogenic zone 927 defined in the ZS9 model proposed by Meletti et al. (2008), which contains our study region. Assuming the same seismogenic zone, slightly different M min , M max , N (annual earthquake rate of M ≥ M min ), and b values were reported by different authors (Convertito et al., 2009;Iervolino et al., 2011). Then, two plausible values of tectonic stress rate of 6.1 and 2.3 kPa/yr are estimated for the VA region (Table 3). These results show a much higher level of background seismicity compared to other studied areas in Italy, such as the Po Plain Adriatic front in the region of Emilia earthquake (Dahm et al., 2015) and the Central Apennines graben system in the region of Aquila (Catalli et al., 2008), where values around 0.7 kPa/yr were calculated.
In a second step, the stress rate induced from the depletion of the oilfield is estimated. The cumulative annual production in the VA oilfield has reached an almost constant value since 2001. Before this date, it was relatively minor. Often, the pore pressure reduction, and thus the depletion of the reservoir layer, correlates more or less linearly with the production (Cesca et al., 2011b). The VA oilfield is more complex, as the reservoir is extremely thick and variable in vertical length, and the porosity and permeability are unusually small. Pore pressure and oil flow is mainly controlled by localized individual fractures rather than a distributed pore space. The permeability can change with varying pore pressure in the fractures. In addition, structural compartments between the CF, ME, and MA are present. Possibly because of this complexity, the measured pore pressure reduction varies between different wells, and the correlation with cumulative production is not strictly linear. Direct measures of the bottom-hole pressure since the beginning of production (late 1990s) are protected by proprietary rights; nevertheless, Improta et al. (2017) stated that a cumulative depletion of about 3MPa involved the southern productive region of the VA oilfield until 2017. The variation of the pressure depletion can be large for thick reservoirs with structural compartments. For a conservative estimate, and to simplify the problem, we assume a pressure drop in the range between 4.5 and 18 MPa (5% and 95% percentiles) through a time span of 16 years over the whole field. The uncertainties of the tectonic stress were taken in the range _ τ T [0.2 kPa/yr, 0.6 kPa/yr]. The induced stress from field pressure depletion is estimated using the nucleus of the strain approach, implemented in a 3D boundary element method (see Dahm et al., 2015). The outline of the field is adapted to the water-oil contact isoline, and the simplified model field was associated with a depth of 3 km, where a regular gridding of 2 km was used. Other parameters used for modeling were a reservoir thickness of 800 m, a Biot's constant of 0.1, and a Poisson ratio of 0.25. The strike, dip, and rake of our target source mechanism were taken from the best MT DC solution ( Table 1). A friction coefficient of 0.6 was used to calculate Coulomb stress changes, which is at the lower range of typical sedimentary rocks under high confining stress (e.g., Byerlee, 1978). The Coulomb stress rate is calculated on a grid between 4 and 17 km depth with 1 km spacing. Figure 5A shows the simulated stress rate induced by depletion for a target fault 14  km depth. At this depth, the largest stress changes induced by the field depletion are expected below the center of the field. The Montesano earthquake is located just at the border of the high stress rate. If the target depth is smaller and closer to the field, the Coulomb stress change increases and the pattern becomes more complex. In any case, the estimated stressing at the hypocenter of the Montesano is in the range of only _ τ D ∼ 0.25 kPa/yr.
The trigger potential (p D ) for an earthquake at a given location, that is, the probability that the rupture nucleation was triggered by field depletion, is described by (see Dahm et al., 2015, Eq. 4): H is the Heaviside function. Values for p D can vary in the range [0,1] with 0 meaning pure tectonic trigger potential and one meaning pure induced trigger potential. The estimated trigger potential for the target event, calculated at a depth of 14 km, is shown in Figure 5B, where we follow the procedure described in Dahm et al. (2015). The pattern shows a smooth field, with diffuse small potential below 0.2. The probability that the Montesano earthquake was triggered by the depletion of the VA oilfield is most likely only ∼6%, where the uncertainties of the tectonic stressing, pore pressure depletion, and earthquake location have been bootstrapped. If the auxiliary plane is assumed to have ruptured, the trigger probability is similar. Note a slight increase in the trigger probability is observed if we consider a shallower earthquake in the same epicentral location of our target event; still, low probabilities are reached, for instance, values around 15% at a depth of 4 km (Supplementary Figure S6).

DISCUSSION AND CONCLUSION
Standard procedures of location and magnitude calculation are carried out automatically by different seismological institutions, at the local or national scale as well as the global scale, including also moment tensor inversion for the largest events (e.g., M w > 4.0). Rupture processes for smaller events are typically neglected and a point source approximation is assumed. However, designed near-regional and local seismic networks, providing important datasets of microseismicity and weak earthquakes, can also reveal important information about the earthquake rupture. The dense network of permanent and temporal seismic stations in the VA region provides us a great opportunity to decipher the genesis of weak events, such as the 2017, M w 4 Montesano earthquake located outside the outer border of the DE. Other than the natural seismicity that is expected in this area, induced seismicity has also been well reported, which was previously associated with variations in the water level of the Lake Pertusillo and derived from wastewater reinjection operations in the disposal well CM2. Moment tensor solutions from natural and induced earthquakes do not reveal significant differences to discriminate the source origin. Typically, the main pattern of normal faulting and some strike-slip focal mechanisms are identified ( Figure 1A), which are associated with the main tectonic regime of NE-SW extension accommodated in the MMFS and EAFS fault systems. Therefore, further seismological and modeling analyses must be conducted to discriminate the rupture plane and assess cases of natural, triggered, or induced seismicity. Generally, an accurate relocation of aftershocks showing some alignments helps to identify the fault plane activated by the mainshock. However, we cannot proceed with this analysis for the Montesano earthquake because only one aftershock was observed. Despite this, rupture directivity has been inferred for the mainshock using this single aftershock as the EGF. Apparent durations, ranging from 0.11 to 0.21 s, have been obtained from S-waves and define an azimuthal pattern, which reveals an asymmetric bilateral rupture with 70% of the rupture propagation in the NW direction and a fault length of ∼0.3 km. The rupture directivity of N310°W is well constrained with small uncertainties and matches the N311°W strike of the fault plane identified in the best double couple of the moment tensor solution, suggesting the rupture plane dipping to the SW. A possible continuation of the EAFS NE-SW dipping normal faults below the allochthone down to the Apulia Platform is still debated (Maschio et al., 2005;Buttinelli et al., 2016 and references therein). Nevertheless, our results could exclude an activation of the MMFS dipping to the NE, and reveal the EAFS as the plausible structure hosting the Montesano earthquake. Further support comes from the observation that the hypocentral location of the target event at ∼14 km depth is consistent with the activation of a deep fault segment associated with the EAFS, which would extend to the basement as reported in previous geological profiles (Menardi Noguera and Rea, 2000;Candela et al., 2015; Figure 1C).
Although hydrocarbon extraction is performed at a shallower layer around 1.8-3.5 km depth, stress perturbation and pore pressure changes may induce or trigger deeper events on preexisting faults previously identified in our seismological analysis. In fact, potential destructive events with an M w larger than 6.5 that occurred in the past (Burrato and Valensise, 2008) could be reactivated by anthropogenic activities, increasing the seismic hazard of the VA region. We have calculated the Coulomb stress rate induced by the depletion of the oilfield to quantify the trigger potential estimated for the Montesano earthquake; our results yield relatively low probabilities below 10%. Note that our analysis does not include the role of the fluids in the reservoir and the surrounding rocks due to the reinjection of water in CM2, which would require specific modeling beyond the scope of the present paper; in particular, hydraulic connections with the VA oilfield should be considered. For instance, some extreme cases in Oklahoma reported triggering mechanisms for earthquakes in the far-field at distances larger than 40 km attributable to fluid disposal wells (Goebel et al., 2017); however, a larger number of injection wells, as well as larger injection volumes, and shallower earthquakes were involved in comparison with our study area. Therefore, we conclude that it is highly unlikely that our target event had originated from the depletion-induced stress rate of the VA oilfield and, rather, a natural cause controlled by the tectonic stress on preexisting faults can be assumed.
On the other hand, we can also apply other qualitative discrimination approaches based on a series of YES-NO questions (Davis and Frohlich, 1993;Davis et al., 1995;Frohlich et al., 2016), which remain used today and have been recently improved (Verdon et al., 2019). A new framework has been proposed including numerical scores to each question considering also some uncertainties, where positive and negative points are assigned depending on whether the answers indicate an induced or a natural cause. Following this scheme (Verdon et al., 2019), we obtain an induced assessment ratio (IAR) of +1.8%, which quantifies whether the overall assessment indicates a natural (−100%) or an induced cause (+100%), considering an evidence strength ratio (ESR) of 98.2% (ranging between 0 and 100%) describing quality and quantity of information used in the assessment (Supplementary Figure S7 and S8). Our seismological analysis allows for reaching a high ESR and evidences a very low IAR score; this suggests that we are unlikely to discriminate whether the Montesano earthquake was originated by an induced or a natural cause. In such a case, the probabilistic discrimination approach described in the present study has more relevance and brings us closer to a better understanding of the genesis of such seismic sources.
In conclusion, we describe a detailed seismological procedure to discriminate between induced, triggered, and natural earthquakes in the VA oilfield, which should be applied together with the previous TLS protocol proposed by the ILG. The relatively large magnitude (M w 4) of the Montesano earthquake and its location close to the external margin of the DE should require the automatic implementation of such an advanced seismological analysis in order to clarify and identify the activation of preexisting or unknown faults. Our results conclude that the Montesano Frontiers in Earth Science | www.frontiersin.org March 2021 | Volume 8 | Article 617794 earthquake activated a deeper fault segment associated with the EAFS close to the basement. The relative low trigger potential based on depletion-induced stress changes discards an induced or triggered event due to the longterm hydrocarbon extraction in the VA oilfield, and it rather suggests a natural cause due to the local tectonic stress.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
JL-C has realized the rupture directivity analysis. TB and SD compiled the data, and framed the geological and monitoring context. TD has realized the probabilistic discrimination analysis. SC has realized the moment tensor analysis. JL-C prepared the initial draft, with contributions and editing from all authors. All authors have contributed to the interpretation and discussion.