The complex geological setting that characterizes the Italian peninsula is the result of different geodynamical processes closely acting in time and space; consequently, today the crust in this zone is characterized by a complex stress field with tightly spaced compressional and extensional regions (Amato and Montone, 1997). The most predominant geomorphologic features in this region are the Southern Alps and the Apennines mountain chains, which are characterized by thrust-and-fold belts originated from the interaction between the European and the Adriatic-African tectonics plates (see e.g., Calamita et al., 1994; Cello and Mazzoli, 1998; D’Agostino et al., 2001; Bertello et al., 2010; Handy et al., 2010; Cazzini et al., 2015; van Hinsbergen et al., 2020). Such an active and complex tectonic setting makes of Italy a seismically active region where on average, every year, more than 2,000 seismic events with magnitude ≥2.0 are located by the Italian national seismic network (see e.g., the Bollettino Sismico Italiano, Pagliuca et al., 2020). Moreover, Italian seismic sequences are generally very complex, often characterised by the occurrence of either foreshocks, multiple mainshocks, or strong aftershocks; this feature of the seismicity in this region has a strongly influence on the seismic hazard (Guidoboni and Valensise, 2015).

In recent years, the interest in the possible influence of anthropogenic activities on seismicity has significantly grown mainly because of a generalized public concern, which has stimulated the development of independent scientific research to support objective policy making (e.g., Ellsworth, 2013; Dahm et al., 2015; van der Voort and Vanclay, 2015; Garcia–Aristizabal et al., 2020). The Italian territory is the scenario of a wide number of underground industrial activities such as oil and gas extraction, geothermal energy production, and gas storage, many of which have been suspected to have direct or indirect causal links with some seismic events located nearby; however, to date there are no unambiguously documented reports of damaging seismic events associated with anthropogenic activities in the country (see e.g., Braun et al., 2018). Probably, the only clear cases of seismicity linked to underground geo-resource development in Italy are the low-magnitude seismicity occurrences recorded in connection with wastewater reinjection at the Costa Molina 2 well in the High Val d’Agri, southern Italy (Valoroso et al., 2009; Improta et al., 2015), and the seismicity recorded near geothermal power plants in Tuscany (Evans et al., 2012).

Among all the anthropic activities having the potential to stimulate earthquakes to occur, the effects of fluid injection or extraction from the crust are probably the processes arising more concern, in particular the activities related with oil and gas production. In Italy, hydrocarbons are found in several oil and gas provinces, most of which are located in the Po plain, the northern Adriatic sea, the southern Apennines, and in Sicily (e.g., Bertello et al., 2010); as a consequence, these are the areas hosting most of the development leases in the country (Figure 1). According to data published by the Italian Ministry of Economic Development (MISE), as of 2019 there were 193 development leases in the country, 127 of which are onshore and the other 66 offshore (UNMIG, 2020).

Discriminating natural from induced seismicity in seismically active regions is a particularly complex task. Early attempts to discriminate induced from natural seismicity were performed, for fluid injection operations, by Davis and Frohlich (1993), and for fluid withdrawal by Davis and Nyffenegger (1995); however, these approaches were mainly based on qualitative assessments.

More quantitative, physically-based and/or stochastic methods for discriminating natural from induced seismicity have recently been proposed in literature (a review can be found, e.g., in Grigoli et al., 2017). For example, Dahm et al. (2015) propose a quantitative probabilistic approach to discriminate induced, triggered, and natural earthquakes, calculating the probability that events have been anthropically triggered/induced from the modeling of Coulomb stress changes and a rate-and-state dependent seismicity model.

Schoenball et al. (2015) analyzed inter-event times, spatial distribution, and frequency-size distributions for natural and induced earthquakes around a geothermal field. Determining the distribution of nearest neighbor distances in a combined space-time-magnitude metric, they identify clear differences between both kinds of seismicity. For example, it is suggested that compared to natural earthquakes, induced earthquakes feature a larger population of background seismicity and nearest neighbors at large magnitude rescaled times and small magnitude rescaled distances. They argue that unlike tectonic processes, stress changes caused by anthropic underground operations occur on much smaller time scales and appear strong enough to drive small faults through several seismic cycles. As a result, it is likely to record seismicity close to previous hypocenters after short time periods.

Zhang et al. (2016) compared moment tensors of both natural and induced events in the Western Canadian Sedimentary Basin. These authors calculated full moment tensors and stress drop values for eight induced earthquakes (magnitudes between 3.2 and 4.4), as well as for a nearby M5.3 event considered as a natural earthquake. This study suggests that, first, it may be possible to discriminate between induced and natural seismicity considering region-specific attributes, as for example the focal depths (which they suggest as the most robust parameter since the induced events in their study area are significantly shallower than most of the intra-plate earthquakes in the Canadian Shield). Moreover, they found a non-negligible (>25%), non double couple component for most of the induced events studied.

Zaliapin and Ben–Zion (2016) analyzed statistical features of background and clustered subsets of earthquakes in California and in South Africa. These authors suggest that, compared to regular tectonic activity, induced seismicity in the analyzed data sets exhibit remarkable features as i) a higher rate of background events, ii) faster temporal offspring decay, iii) higher rate of repeating events (i.e., earthquakes located in the rupture area of some previous earthquakes, but that occur at times far exceeding the typical duration of an aftershock series), iv) larger proportion of small clusters, and v) larger spatial separation between parent and offspring.

Discriminating human-induced from natural seismicity is therefore a difficult problem for which virtually any solution can be a source of controversy. For this reason, the evaluation of possible interactions between seismicity and hydrocarbon production, if possible, should rely on multidisciplinary analyses such as e.g., detailed physically-based modelling complemented by sophisticated stochastic methods able to provide probabilistic assessments and to take uncertainties into account. It is worth noting however that the ways in which these interactions may occur are complex and their identification in a context characterized by high levels of naturally-occurring seismicity is not straightforward. Moreover, given the relatively high number of development leases active in Italy, performing such analyses at the national scale may be considered an intractable problem.

These reasons pushed us to explore the possibility to implement large-scale screening methods aiming at tracking measurable phenomena, such as e.g., changes in seismicity rates, that plausibly could occur if notable interactions between underground human operations and nearby seismicity sources are actually occurring in a given area. Spatial and temporal correlation between human activity and event rates is usually considered a key parameter to suspect possible relationships between seismicity and underground anthropic activity (e.g., Shapiro et al., 2007, 2010; Cesca et al., 2014; Leptokaropoulos et al., 2017; Garcia-Aristizabal, 2018; Molina et al., 2020); for example, significant changes in seismicity rates with respect to background seismicity, as well as the spatial and temporal correlation of gas injection operations and seismicity were analyzed by Cesca et al. (2014) to suggest a possible case of triggered/induced seismicity near an offshore platform used for gas storage in Spain (the Castor project).

However, underground human-induced perturbations (as e.g., pore pressure variations due to fluid injections) can produce changes at large distances and/or with large temporal delays, potentially causing earthquakes to occur several kilometers away as well as months/years after the industrial operations have stopped or reached the maximum peak (e.g., Mulargia and Bizzarri, 2014); likewise, natural seismicity may also occur within few kilometers from industrial sites. Therefore, in seismically active regions (such as in Italy), spatio-temporal correlations between industrial activity and significant changes in seismicity rates with respect to background activity by themselves usually do not provide irrefutable proofs of causal relationships between hydrocarbon production and seismic activity. Despite this, we argue that studying such correlations has a remarkable added value since it gives us the possibility of performing large-scale, systematic analyses of a huge amount of seismic and production data and, under the working hypotheses considered, to identify hotspot areas where it could be possible to perform, in a later stage, more detailed research to verify possible causal relationships.

The industrial data publicly available for this study consists of a time series of hydrocarbon production volumes; for this reason, our analyses are particularly focused on studying the possible effects on seismicity of stress perturbations caused by fluid withdrawal processes. The working hypothesis in this work therefore starts from assuming that fluid withdrawal from the crust may induce deformations in the surroundings of the host rock (especially in depleting reservoirs); the magnitude of such deformations will depend on multiple factors such as, for example, the litology, the structural geology, the volume and rate of fluid removed, the geomechanical features of the reservoir, and its behavior during the fluid withdrawal process, among others. The deformations may in turn alter the local stress field and, as a consequence, stimulate or inhibit seismicity in the surroundings.

We assume that in absence of other stress perturbation sources, the regional release of background seismicity is predominantly influenced by the regional tectonic stress field; in such a context, it can be expected that a steady stress field should tend to generate background seismicity with stationary rates; however, if other natural (e.g., hydrogeology: see Hainzl et al., 2006; Pintori et al., 2021) or man-made (e.g., pressurized fluid injections: see Shapiro et al., 2007; Garcia-Aristizabal, 2018) processes are able to perturb the local stress field, then it is possible that the rate at which seismicity is released in that specific area can be altered. In such a case, slight deviations from stationarity could possibly be measured.

In this work we are interested in identifying significant departures of background seismicity from a stationary behavior around productive oil and gas development leases in Italy. We are particularly interested in exploring cases in which the long-term hydrocarbon production may have left a measurable footprint on the release of background seismicity (as e.g., by processes related to reservoir depletion), and to test whether possible changes in the rate at which background seismicity is released in such areas are correlated with the main changes in the hydrocarbon production patterns.

The article is structured as follows: first we present the seismic and the hydrocarbon production data available for this study. Second, we present the methodological approach used in order to identify zones with possible anomalies in seismicity rates concomitant with significant changes in oil and gas production. Finally, we present the results and discuss the importance and limitations of these findings.

For this study we use a national-wide seismic catalog containing earthquake locations and magnitudes, as well as the most detailed public oil and gas production data from development leases in Italy. All the used data are freely accessible from public sources (Section 7 for details).

We perform a large-scale screening of the behavior of background seismicity (i.e., the events considered independent, as described in Section 3.3) in areas around oil and gas production sites in Italy. Our goal is to attempt i) to identify correlated changes between seismicity rates and hydrocarbon production, and ii) to test the significance of these seismicity rate changes.

First, a pre-processing step is performed in order to identify the areas around the target leases where a reasonable amount of both seismic and production data is available. Afterwards, we proceed with the proper data analyses, namely i) identification of independent background seismicity in the lease area, ii) selection of background events located within a distance δx from the production wells, and iii) test the significance of possible seismicity rate changes correlated with changes in hydrocarbon production. These steps are explained in the following paragraphs.

One of the main issues faced to perform the data analysis was to concentrate our efforts in areas where both seismic and production data were sufficiently representative to avoid, as much as possible, data-driven biases in our results. For this reason, and given the time span covered by the seismic and production data available for this study, we avoided analyzing any region in which the minimum data requirements defined in the Methods section were not accomplished. In particular, the minimum length of the time window of seismic data before and after t _m, the time at which the maximum production is reached, becomes one of the main constraints, forcing us to not considering about 56% of the leases for which production data are actually available (that is, 57 out of the 102 leases). In such cases, the peak production occurred too close to the end or the beginning of the seismic catalog, which prevents an accurate estimation of the pre– and post– t _m seismicity rate. On the remaining 45 leases (39 produce gas only, two produce oil only, and four produce both gas and oil), we applied the analyses described in Section 3. The spatial distribution of the analyzed leases cover different areas of the country, including the western Po plain, the northern and central Adriatic sea, southern Apennines and Sicily.

The completeness magnitudes and b values determined for the regional seismicity around each lease are summarized in Table 1 (for the oil-producing leases) and Table 2 (for the gas-producing leases). It is worth noting that, in general, the completeness magnitudes tend to be relatively high due to our interest in having seismic data sets as long as possible in time; in fact, the M _c has been selected for the seismicity around each lease so that completeness is ensured, at least, since 1985. The seismic catalogs, composed of nearby regional events with magnitudes above the M _c specifically determined for each study area, are then declustered using the NN clustering analysis technique (Section 3.3) to obtain a sample of background seismicity composed of independent events at regional scale.

It is worth noting that the declustering step deserves some specific considerations. We argue that in the context of this work, declustering is an important step given our particular interest in detecting seismicity rate changes induced by eventual stress changes caused by depleting reservoirs located in areas where regional stresses are assumed to generate stationary background seismicity. Therefore, the proposed approach can hardly be applied to cases where induced seismicity tends to be strongly clusterized in time and space, as for example cases of seismicity induced by pressurized fluid injections. Moreover, uncertainties related to the effects of declustering could overshadow the reliability of the seismicity rate variations identified with the proposed approach; to mitigate this risk, we also discuss the results that can be obtained by analyzing the full complete catalog (that is, considering events above M _c without declustering); such results, along with the data sets (both the full and the declustered catalogs) are available in the Supplementary Material.

After declustering, we identify background events located within δx = 5 and δx = 10 km from the production wells in a given lease, and select the events before and after the respective t _m. In this way we calculate the number of events before (n _pre) and after (n _post) the peak production time, t _m, as well as the respective time window lengths T _pre and T _post (see e.g., Figure 6D for reference). The same procedure is performed using the full complete catalog (Supplementary Material) to compare the results obtained with and without the declustering processing.

Using these data, we first calculate the seismicity rate (in terms of average number of events/year) observed before and after t _m for events located within a distance δx from the producing wells. Comparing the seismicity rates before and after the peak production, it is possible to highlight the areas with seismicity rate variations. For example, Figure 7A and Figure 8A show the location and estimate of seismicity rate variations around, respectively, the six oil-producing and 43 gas-producing leases analyzed in this study when considering seismicity within δx = 5 km from production wells; colors indicate the behavior of the seismicity rate (i.e., increase, decrease, unchanged) before and after t _m. Regarding the oil-producing leases, the seismicity rate around three sites results null and unchanged before and after t _m (green squares in Figure 7A; names can be seen in Table 1), in 2 cases (Masseria Verticchio and Val d’Agri) the seismicity rate before the peak production results higher than the seismicity rate in the time window after t _m (blue squares in Figure 7A), whereas in one case (Ragusa) we observe the opposite situation (that is, the seismicity rate before is lower than the seismicity rate after t _m, red squares in Figure 7A). Regarding the gas-producing leases, we observe that in 20 sites the seismicity rate does not change (18 of which have zero events before and after t _m). These sites are represented as green squares in Figure 8A. In the remaining 23 leases (names can be seen in Table 2), a seismicity rate variation has been detected: nine sites exhibit an increase (red squares in Figure 8A), whereas 14 sites exhibit a decrease in the rate (blue squares in Figure 8A). When considering δx = 10 km, one oil-producing and 10 gas-producing leases do not exhibit seismicity rate changes, in 4 and 22 leases (oil and gas, respectively) the seismicity rate before t _m results higher than the seismicity rate after the peak production, whereas the opposite behavior is observed in one oil- and 11 gas-producing leases.

We then evaluate the significance of these seismicity rate variations using the binomial test described in Section 3.4. For the areas exhibiting a reduction in seismicity rate after t _m, we calculate p ₁ (Eq. 3) to evaluate the significance of the observed seismicity rate reduction in the time interval T _post = [t _m, t ₂] with respect to the rate observed in the interval T _pre = [t ₁, t _m]. Likewise, for the areas exhibiting an increase in seismicity rate after t _m we calculate p ₂ (Eq. 4) to evaluate the significance of the observed increase in the time interval T _post with respect to the rate observed in the interval T _pre. The results of these calculations, considering seismicity located within both δ x = 5 and δ x = 10 km from production wells, are presented in Table 1 for the oil-producing and in Table 2 for the gas- producing leases.

The resulting probability values (p ₁ and p ₂) provide an indication about how likely it is to observe the number of events n _post under the null hypothesis (which basically reflects what would be expected in case of stationary seismicity in the whole interval [t ₁, t ₂], see Figure 6D for reference). Therefore, the lower the p value, the more unlikely is to observe n _post under such hypothesis (and therefore the more evidence in favor of a significant seismicity rate change). Considering for example a significance level s.l. = 0.05 we find that, on the one hand, none of the cases in which an increase in seismicity rate was observed after t _m is statistically significant (considering seismicity located within both 5 and 10 km from production wells). On the other hand, regarding the cases exhibiting a seismicity rate decrease after t _m, our results suggest that the observed seismicity rate change is statistically significant for one oil-producing lease (the Val d’Agri) and two gas-producing leases (Fiumetto and Rocca Cavallo). The location of all the analyzed leases, classified by the calculated p values considering as reference a s.l. = 0.05, are shown in Figure 7B for the oil-producing and Figure 8B for the gas-producing leases. The Val d’Agri lease is located in the Basilicata region, in the Southern Appenines, whereas Fiumetto and Rocca Cavallo leases are located in Sicily, near Etna Volcano.

The leases for which the seismicity rate change is considered statistically significant are identified using events within both δx = 5 and δx = 10 km (see Tables 1 and 2). It is worth to note that in most of the cases the rate change (i.e., increase or decrease) is coherent for seismicity selected considering both δx values. However, in a few cases of data selected around some gas-producing leases (5 sites out of 43, namely A.C 28. EA, Bronte-S. Nicola, Monte Morrone, Recovato, and S. Andrea), there is a change in the observed trend of seismicity rate variation (for example, a decrease observed for the seismicity located within 5 km contrasts with an increase observed when selecting events within a 10 km distance); discrepancies in such areas are due to instabilities mainly caused by either low seismic activity (such as, e.g., in the Po plain or offshore in the Adriatic sea) or sites close to seismically active sources (such as, e.g., very active tectonic or volcanic areas). It is worth noting however that in all these cases, as well as for cases in which no events are identified within the first 5 km, the observed rate changes are in any case not statistically significant.

In order to compare these results with the results that would be obtained if we do not use the background events only, the same analyses were performed using the full, complete catalog (i.e., considering events above M _c). It is worth noting that using the full catalog may somehow violate the methodological assumption considered for the statistical test (that is, assuming that in the absence of external stress perturbations, seismicity is expected to be stationary). Nevertheless, in this case such comparison is useful to demonstrate that the significant seismicity rate changes found in this work are not artificially created by the declustering procedure. These results are reported in Supplementary Tables S1 and S2 and in Supplementary Figures S1 and S2 in the Supplementary Material. When using the full catalog, we still find evidence of significant seismicity rate changes in the same three leases described in the previous paragraph (i.e., Val d’Agri, Fiumetto, and Rocca Cavallo). Nevertheless, it is worth observing that using the full catalog a few additional zones show some evidence of significant rate changes, especially when considering seismicity located at distances within 10 km from production wells; we consider however that such results are unreliable (for example, in the same area the rate may change e.g., from rate increase to rate decrease– when considering events at different distances, see for instance Tables S1 and S2 in the Supplementary Material); such unstable results are most probably due to the presence of clusters of aftershock events randomly located around some leases.

In this work we monitor significant departures of background seismicity from a stationary behavior around active oil and gas development leases in Italy. We collected oil and gas production data from 102 leases in the period 1980–2017, and we used the seismic data from the HORUS catalog (Lolli et al., 2020b). After a close and conservative inspection of the available data, it has been possible to implement the proposed analyses in six oil-producing and 43 gas-producing leases in the country (including about 44% of the leases from which production data was available). We identify statistically significant seismicity rate changes (considering a s.l. = 0.05) concomitant with outstanding changes in hydrocarbon production (i.e., before and after peak production) in one oil-producing lease, Val d’Agri (Figure 7B), and two gas-producing leases, Fiumetto and Rocca Cavallo (Figure 8B). In all these cases, the seismicity rate after the peak production results significantly lower with respect to the seismicity rate observed before it.

It is worth to remind that spatial and temporal correlations between background seismicity rates and industrial activity, as the one highlighted in these areas, do not constitute an absolute proof to establish a causal relationship between hydrocarbon production and seismic activity. Therefore, we stress that the main value of this finding is to highlight these areas as hotspot zones deserving further detailed analyses to verify (or discard) possible causal relationships. If the oil or gas production in the Val d’Agri, Fiumetto and Rocca Cavallo leases played a role to alter shallow earthquake occurrences in these areas, then more sophisticated, physically-based studies are required to understand if the observed changes in seismicity rates are actually an observable consequence of physical mechanisms able to generate such changes. However, such analyses require access to detailed, geological, structural and geomechanical information.

Val d’Agri, located in the Basilicata region, is a particularly interesting case (Figure 9). This area hosts the largest onshore oil and gas field in Europe, where possibly-induced seismicity has been detected in connection with wastewater reinjection at the Costa Molina 2 well (e.g., Stabile et al., 2014b; Improta et al., 2015; Buttinelli et al., 2016; Improta et al., 2017). Moreover, some authors reported clustered seismicity located to the south of the nearby Pertusillo artificial water reservoir, whose origin has been suggested to be induced by rapid water level changes of the Pertusillo impoundment (Valoroso et al., 2009, 2011; Stabile et al., 2014a).

Figure 9 shows the location of the Val d’Agri lease, the production wells, as well as the regional seismicity in the area (gray squares), the events with magnitude above M _c (all circles), the background seismicity identified using the NN clustering analysis technique (colored circles), and the selected background seismicity located within a distance δx = 5 km (Figure 9A, red circles) and δx = 10 km (Figure 9B, yellow circles) from production wells. Only a few selected events are located in the southern part of the lease, close to the cluster of seismicity located south of the Pertusillo lake. Figure 9C shows the time series of annual oil and gas production from the Val d’Agri, and the plot of event times and magnitudes of selected seismicity. What we actually observe is a clear reduction in the number of shallow (z ≤ 15 km) background seismic events just after the oil production in Val d’Agri reached its maximum in 2005.

A similar significant change in the seismicity rate has been detected in the surroundings of Fiumetto and Rocca Cavallo leases, which are located in Sicily, near Etna Volcano (Figure 10). These two leases are located in an area where other five gas-producing leases operate (namely Bronte-S. Nicola, Case Schillaci, Gagliano A, Gagliano B, and Samperi), forming a “cluster” of production leases distributed in a relatively small area (about 400 km², see e.g., Figure 10A). The results of the correlation analysis for Bronte-S. Nicola and Case Schillaci are also reported in Table 2 and Figure 8, whereas the last three leases were not included in this study because they did not meet the data requirements defined in Section 3.1. For completeness, we also reproduced the study using the full catalog of events (that is, without removing the most likely aftershock sequences identified using a declustering approach). Such a test allowed us to demonstrate that the observed reduction in seismicity rate after hydrocarbon peak production in the three production leases outlined in this study is not artificially introduced by the declustering procedure.

Figure 10 shows the geographical location and gas production for Fiumetto and Rocca Cavallo leases, and the background seismicity identified for this zone. We also plotted the location of all the other leases in the area, as well as the gas production of the other two cases included in this study (Bronte-S.Nicola and Case Schillaci). It is worth noting that Gagliano has been producing gas for much more longer time with respect to the other four, but it was not possible to include it in the analysis since the peak production precedes the time interval considered in this study. In these plots we highlight all the background events located within 5 and 10 km from the production wells in Fiumetto (Figures 10A,B) and Rocca Cavallo (Figures 10D,E) leases.

In this context, a correlation analysis between seismicity and hydrocarbon production for this area is particularly complex, probably requiring an integrated analysis considering all the leases together. Analyzing the results obtained considering each lease independently, similar to what is observed in Val d’Agri, we observe a decrease in the rate of shallow seismicity after the peak gas production in both Fiumetto and Rocca Cavallo (Figures 10C,F). A similar behavior is observed in the Case Schillaci site as well, but the p value in this case is not statistically significant at the s.l. adopted (neither considering δx = 5 km nor δx = 10 km). Finally, and as mentioned before, the results for Bronte-S. Nicola are contrasting when considering different δx values, as a probable effect of intense, shallow seismicity related to activity at the nearby Etna volcano. It is worth noting that the higher peaks in gas production in this area (in the time interval considered in this study) are reached in Fiumetto, Rocca Cavallo and Bronte-S. Nicola between ∼1997 and ∼2005, a period in which the reduction in the number of events starts to be evident in the zone. However, it should be kept in mind that this area is a particularly complex case due to different factors, such as the closeness of different active leases and the proximity of a seismically active volcano, therefore any further analysis probably requires taking into account the whole cluster of leases together (considering also the influence in this particular area of possible stress perturbations caused by activity at the nearby Etna volcano).

In the areas highlighted by low p values in this study we observe that the seismicity after t _m occurs at a lower rate with respect to the seismicity occurring before the peak production. In the case of a link with fluid withdrawal, such observation could result from either an increase in seismicity rate associated with the phase preceding t _m , when the fluid withdrawal rate had an increasing trend (see e.g., Figure 5E for reference), or because seismicity is somehow inhibited as the fluid withdrawal process goes on (which results particularly evident after t _m ). An analysis of the seismicity rate before and after the start of fluid withdrawal operations would be useful to understand which of the two scenarios dominates; however, the seismic data for times preceding the start of operations is too scarce for reliable and systematic analyses. Looking however at the seismicity and production data in Figures 9, 10, it is interesting to note that background seismic events are anyhow present before the start of production in the highlighted cases, and this observation contrasts with the apparent seismicity rate drop observed after the peak production. Moreover, some authors have reported evidences of reservoir depletion processes in the Val d’Agri field; for instance, Improta et al. (2017) reported lacking seismicity correlated with low Vp/Vs zones and interpreted these zones as possible depleted sectors of the Val d’Agri production reservoir. They conclude that these observations agree with increased frictional resistance on preexisting faults within the hydrocarbon reservoir caused by the increase of the effective normal stress resulting from pore pressure depletion induced by significant fluid withdrawal (Improta et al., 2017). These observations lead us to think that the change in rate would primarily be due to a decrease in the number of earthquakes after the peak.

In an attempt to better understand if deformative phenomena associated with fluid withdrawal from the crust may explain an apparent inhibition of shallow seismicity, we analyzed a hypothetical case of deformation related to reservoir depletion and calculate the Coulomb stress changes on fault planes optimally oriented for failure (see e.g., Lin and Stein, 2004; Toda et al., 2005). In practice, we assume the contraction of a sub-horizontal dike in an elastic halfspace. For this theoretical exercise we use approximate geomechanical properties of the Val d’Agri zone taken from literature. The modeled dike has an area that roughly covers the area outlined by the production wells active in the Val d’Agri lease (that we assume as an approximate proxy of the actual areal distribution of the reservoir), and is located at a depth of ∼4 km below surface (information about the depth of the reservoir in Val d’Agri can be found in e.g., Stabile et al., 2014b; Improta et al., 2017).

In agreement with the regional stress field inferred for this region, we assume an extensional tectonic environment (characteristic of normal faulting regimes) with S V > S H max > S H min , and where the minimum horizontal stress is oriented in the direction NNE-SSW (Cucci et al., 2004; Montone et al., 2012; Improta et al., 2017). The results of the Coulomb stress change resolved on optimally oriented planes for this setting is shown in Figure 11; cold (blue) colors in Figure 11 indicate zones with negative Coulomb stress, outlining areas in which the modeled deformation would tend to inhibit seismicity. Such results depict a simple example in which fluid withdrawal operations are allowed to induce a contraction (compaction) of the reservoir; in such a case, the resulting deformation would tend to inhibit seismicity at least in the shallow crust around the area of a depleting reservoir. Therefore, under some circumstances it seems feasible to observe a reduction in seismicity rates following the period of largest hydrocarbon production. Nevertheless, more accurate analyses using detailed geomechanical information and based on the results of a more adequate poroelastic model (as done e.g., by Segall, 1989; Segall et al., 1994) would be required for better capturing the nature of such processes.

In areas characterized by high levels of natural seismicity, the identification of human-induced seismicity is a difficult task for which, virtually, any result can be a source of controversy. In Italy, a relatively high number of oil- and gas-producing leases have been operating in the last decades, many of which are located in the surroundings of seismically active regions (e.g., Sicily, and the central and southern Apennines). Besides hydrocarbon production, underground waste water injection is also performed in different areas of the country, but these activities were not considered in this study because injection-related data are not publicly available.

Therefore, our analyses are focused on areas where fluids (oil and gas) are withdrawn from the crust. Performing detailed, physically-based analyses at the national scale to identity cases of anthropogenic seismicity in Italy may be considered an intractable problem. Consequently, we implemented a large-scale screening procedure aiming at tracking measurable phenomena, such as, e.g., changes in seismicity rates, that plausibly could occur if notable interactions between fluid withdrawal from the crust and nearby seismicity sources are actually occurring in a given area. We stress however that spatial-temporal correlations between proxies of industrial activity and significant changes in seismicity rates do not provide an absolute proof of causal relationships between hydrocarbon production and seismic activity; rather, studying such correlations gives us the possibility of performing large-scale, systematic analyses of a huge amount of seismic and production data, and to identify in this way hotspot areas where to focus more detailed research to verify, in a later stage, possible causal relationships.

In this context, we analyzed seismicity rates before and after peak production in six oil-producing and 43 gas-producing leases, and evaluate the significance of potential seismicity rate changes. Such cases are about the 44% of the development leases active in Italy (the other 56% were not analyzed due to data limitation constraints). The main findings resulting from this study can be summarized as follows:

• When considering the background seismicity located within 5 km from the production wells, about 50% of the oil-producing and 46% of the gas-producing leases analyzed in this study do not show any change in seismicity rate before and after the time at which the peak production is reached; in most of these cases the seismicity rate is zero (basically due to the short distance considered). The percentage of oil- and gas-producing leases with no observed change in the seismicity rate reduces respectively to about 17 and 23% when considering seismicity located within 10 km from production wells.

In conclusion, our analyses highlight areas near some oil- and gas-producing leases in Italy where the seismicity rate reduces after peak production is reached (as compared to the seismicity rate preceding it). We emphasize however that our results just put in evidence a correlated change between the rates of shallow seismicity and hydrocarbon production in these areas, and that assessing actual causal relationships between these two processes will require further detailed, physically-based research. Despite this, we argue that should a physical link exist between these processes, the observed seismicity rate reduction could either be due to increased seismicity during the progressive increase in production rate before t _m , or to actual seismicity rate reduction after t _m . This second scenario would put in evidence possible processes of seismicity inhibition. Considering that the occurrence of seismicity before the start of hydrocarbon production in the hotspot areas contrasts with the reduction of events observed after the peak production, we suspect that the seismicity inhibition is a plausible hypothesis in these cases. With a simple theoretical exercise based on modelling Coulomb stress changes we showed that, at least under some simplified conditions, inhibition of seismicity is actually possible; nevertheless, more accurate models (e.g., using poroelastic theory) are required for better understanding the nature of such processes.

Our observations draw attention to an interesting research problem: the characteristics and implications of increased seismicity caused by anthropogenic activities (e.g., pressurized fluid injection) have so far had a prominent role in research on induced seismicity; the implications, instead, of anthropic processes potentially capable of inhibiting seismicity—in particular, on seismic hazard in seismically active regions—have received less attention, or have been neglected. We consider that more efforts to study the mechanisms and the possible consequences of anthropogenically-driven seismicity inhibition, as well as to document possible cases of such phenomenon, are required.