Combined approach of poroelastic and earthquake nucleation applied to the reservoir-induced seismic activity in the Val d’Agri area, Italy

Antonio P. Rinldi, Luigi Improt, Sestin Hinzl, Flmini Ctlli,d, Lu Urpi,Stefn Wiemer

a Swiss Seismological Service, ETH, Zurich, Switzerland

b Istituto Nazionale di Geofisica e Vulcanologia, INGV, Rome, Italy

c Deutsches GeoForshung — Zentrum, GFZ, Potsdam, Germany

d RatePAY GmbH, Berlin, Germany

ABSTRACT In this work, an approach is developed to study the seismicity associated with the impoundment and level changes of a water reservoir (reservoir induced seismicity — RIS). The proposed methodology features a combination of a semi-analytical poroelastic model with an earthquake nucleation approach based on rate-and-state frictional law. The combined approach was applied to the case of the Pertusillo Lake, located in the Val d’Agri area (Italy), whose large seasonal water level changes are believed to induce protracted micro-seismicity (local magnitude ML ＜3). Results show that the lake impoundment in 1962 could have produced up to 0.5 bar (1 bar = 100 kPa) changes in Coulomb failure stress (ΔCFS),while the seasonal water level variation is responsible for variation up to 0.05 bar.Modeling results of the seismicity rates in 2001-2014 show that the observed earthquakes are well correlated with the modeled ΔCFS. Finally, the reason that the seismicity is only observed at southwest of the Pertusillo Lake is provided, which is likely attributed to different rock lithologies and depletion caused by significant hydrocarbon exploitation in the northeastern sector of the lake.

Keywords:Reservoir induced seismicity (RIS)Poroelasticity Rate-and-state frictional law Pertusillo lake

1. Introduction

A water reservoir affects the underlying crustal stress state through the poroelastic response to the weight of the water volume stored and by the consequent fluid propagation (Talwani, 1997;Gupta, 2002). The perturbation of crustal stress state has been in some cases associated with small-to-large seismic events, with maximum seismicity magnitude up toM= 6.3 recorded in the largest confirmed case of reservoir-induced seismicity(RIS),which took place at the Konya reservoir in India in 1967. In recent years,understanding of the physical mechanism of RIS has emerged as pivotal case of discussion in the seismological community given the increasing reports of damaging (M> 5) earthquakes worldwide(Chen and Talwani, 1999; Barros et al., 2018; Grasso et al., 2018;Huang et al., 2018). Water level variation in a lake causes both hydraulic and mechanical effects on the rock volume.The weight of the water column affects volumetric and shear strains, which in turn could affect pore pressure.At the same time,the lake bottom constitutes a source of fluid altering the conditions at depth,as the pore pressure increases.Its mechanical and hydraulic effects could lead to an increase of pore pressure and changes in Coulomb failure stress (ΔCFS) with elapsed time. Fig. 1 summarizes the possible reservoir effects on pore pressure at depth in response to water level changes.

Fig.1. Schematic of coupled processes occurring at depth as the water level in a reservoir changes with elapsed time and affecting the pore pressure in the seismogenic volume.

RIS is generally investigated by means of seismological and hydrogeological techniques, which are essential in both assessing possible basic mechanisms of fault reactivation and identifying the anthropogenic origin of seismicity.The combined approach is even more advantageous when RIS activities occur in seismic active regions, possibly coexisting with other forms of anthropogenic seismicity (e.g. mine exploitation, oil production, and wastewater disposal).This context can be found in the Val d’Agri area,located in the southern Apennines seismic belt,where RIS and fluid-injection induced seismicity is associated to the production of a large oilfield(Valoroso et al., 2009; Stabile et al., 2014; Improta et al., 2015;Buttinelli et al., 2016). The artificial lake Pertusillo in Val d’Agri(Italy) had shown protracted seismic activities for several years after the initial filling in 1963. More than 800 small-magnitude events occurring between 2001 and 2004 (local magnitudeML＜3,and magnitude of completenessMc=1.1)were accurately located by local permanent and temporary networks (see Fig. 2).This constitutes a unique case of monitored RIS, with high-quality seismic logs providing key information for better interpretation of observed seismic activities (Valoroso et al., 2011; Stabile et al.,2014; Improta et al., 2017). Seismicity concentrates on waterfilled fractured Mesozoic limestone sealed by Cenozoic clayey sequences, at depths of 2-5 km to the south of the lake (Improta et al., 2017). During the same period, the lake water level fluctuated on average tens of meters between summer and winter periods.The observed seismicity rate positively correlates with these seasonal water level oscillations(Valoroso et al.,2009;Stabile et al.,2014). HighVp/Vs(VpandVsare P- and S-wave velocities, respectively) anomalies, resolved by high-resolution 3D/4D (three- and four-dimensional)local earthquake tomography (LET), suggest the essential role of pressurized fluid behind the seismicity (Valoroso et al., 2011; Improta et al., 2017). The state-of-the-art, however,fails in providing a consistent physical modeling for the RIS in the region. An accurate modeling is needed to properly discriminate the complex observed patterns.Moreover,given that the Pertusillo Lake micro-seismicity on active normal faults is capable of producing earthquakes on surface faults(M>6)(Improta et al.,2010),modeling can help in quantifying the potential hazards posed by the lake itself. In this work, the first modeling is attempted to understand the physical mechanism behind the seismicity at Pertusillo Lake.The approach accounts for a semi-analytical model with Green’s function solution given for a homogeneous, poroelastic half-space, and considers the decoupled approximation when solving the governing partial differential equations (i.e. elastic stress influences the pore pressure, but not vice versa). Then the calculated effective stresses are used to compute seismicity rate changes through a rate-and-state nucleation model. Stress and strain are calculated for the observed transient evolution of the water reservoir level, and the calculation allows to compute ΔCFSand to identify the mechanism of failure as a function of time.The current approach for stress calculation is only valid for static loading (i.e. water dam); while the nucleation approach is more generic and could be used with any ΔCFScalculation given the orientation of potential seismogenic structure. The current modeling approach is not meant to be a general tool for any hydromechanical analysis of underground processes.It only provides some insights on the impoundment and seasonal variation of water level in an artificial lake.Albeit simplified,the proposed model can easily explain the temporal evolution of seismicity in agreement with the reactivated structures at the south of the lake, and pose some important points of discussion:(i)If a temporal correlation is confirmed, why does the seismicity concentrate mostly on the south of the lake given that active normal faults optimally orient to slip located on both sides of the lake?and(ii)Is the oil production going to interfere with the spatio-temporal distribution of RIS? A conceptual model is proposed to explain the absence of seismicity and conclude that the oil production could have hindered the protracted seismic activity on the northeast of the lake.

Fig.2. (a)Observed Val d’Agri seismicity in 2001-2014.Highlights in figure also show two bordering fault systems,the oil field,the production wells,and the CM2 well that caused injection-induced seismicity.High-quality focal mechanisms for four earthquakes of the southwestern cluster are also presented(Improta et al.,2017).(b)The Pertusillo Lake water level (m) and the monthly observed number of earthquakes in southwestern region (longitude 15.87°-15.95°; latitude 40.2°-40.29° - 531 earthquakes, with ML ≥ 1.0 and depth ＜ 7 km in 2001-2014,showing a correlation behavior.(c)Schematic geologic vertical cross-section across the Val d’Agri basin(see Improta et al.,2017;Mazzoli et al.,2013)(1)Mesozoic-Tertiary fractured carbonates of the Apulian platform(hydrocarbon/aqueous reservoir);(2)Mesozoic basin rocks;(3)Mesozoic limestone of the western platform;(4)Tectonic melange formed by ductile,low-permeable Mio-Pliocene sediments(caprock);(5)Mio-Pliocene thrust-sheet-top sediments;and(6)Quaternary deposits of the Val d’Agri basin. Earthquakes within 3 km from the section are projected.

2. Seismicity at the Val d’Agri basin

The Val d’Agri Quaternary basin is bordered north by the southwest-dipping east Agri fault system (EAFS) and south by the northeast-dipping Monti della Maddalena fault system (MMFS);both fault systems are characterized by a normal-fault kinematic and show evidence of Late Pleistocene—Holocene activity(Improta et al., 2017).

Previous studies mainly focused on the correlation between Pertusillo water level oscillations and seismicity in the Val d’Agri area (Fig. 2a) from 2001 onwards. Prior to the installation of local monitoring networks in 2001, no permanent stations were operational in the Val d’Agri area. The first permanent station was installed in 2004,about 12 km to the south of the lake,followed by three other stations in 2006 along the MMFS western ridge. As a consequence,records of seismicity prior to 2001 are not reliable for the study of RIS at Pertusillo Lake, due to high magnitude of completeness (i.e.Mc> 2.5) and high uncertainties on events’location. Based on instrumental logs, it can be concluded that the reservoir impoundment,as well as the following seasonal loading/unloading stages, induced large seismic events. Since 1962, the closest moderate earthquake was aM= 4.7 strike slip event in 1971 at a distance of about 15 km to the northwest of the lake(Gasparini et al.,1985;Cucci et al.,2004).Valoroso et al.(2009)first observed a temporal correlation between an intense swarm-type micro-seismicity (ML＜2.8) in the southwest of the Pertusillo lake and the significant water level changes (～15 m), suggesting that this seismicity was likely reservoir-induced. Their study focuses on a time window spanning from May 2005 to June 2006,which was covered by a dense passive seismic survey.

Stabile et al. (2014) and Telesca et al. (2015) analyzed the correlation between the water level changes and seismicity extending the observation from 2005 to 2012 using data from a monitoring network located within the oil field.Fig.2b shows a comparison of the earthquake rate with the water level variation in 2001-2014(i.e. since the deployment of the local monitoring network). The earthquakes distribution (Fig. 2a) and the analysis of focal mechanisms highlight the presence of several clustered events (Improta et al., 2017). The northeastern part features reactivation, since 2006, of inherited blind thrust correlated with the wastewater injection at the disposal well Costa Molina 2 (CM2, see Fig. 2a—c) as highlighted by Improta et al.(2015)and Buttinelli et al.(2016).The injection well is located at marginal portion of the Val d’Agri hydrocarbon field,which has been in production since the late 1990s and caused a depletion of the carbonate hydrocarbon reservoir by few tens of bar in about 15 years (Improta et al., 2017). The southeastern clusters in Fig. 2a—c suggest the reactivation of northwest-trending, northeast-dipping splays of the MMFS along the border of the Val d’Agri basin,optimally oriented in the present state of stress, characterized by NE-trending extension (Improta et al., 2017). Interestingly, the southwest-dipping EAFS showed no seismicity in the monitored period, although being located at the same distance from the lake.

3. Modeling approach

3.1. Poroelastic pressurization

A Green’s function approach was used to model (i) the elastic stress changes induced by a static surface loading and(ii)the pore pressure distribution in a poroelastic half-space and its time evolution.The stress changes can be estimated using solutions derived for the classical elastic theory.Here we refer to the 3D Boussinesq’s solution for a vertical point loadfin an elastic half-space (Jaeger et al., 2007). The stress tensor σ due to a finite surface load is then obtained by integrating over a number of point forcesfidistributed over the surfaceA, which represents the water reservoir:

whereGBis the specific Green’s function tensor for the Boussinesq’s solution,Ais a dimensionless fault constitutive friction parameter usually estimated as ～0.01(Dieterich,1994;Dieterich et al.,2000).The poroelastic semi-analytical model,under the approximation of a decoupled problem (Roeloffs, 1988), assumes that the elastic stresses influence pore pressure but not vice versa, and the governing partial differential equation is reduced to the following inhomogeneous diffusion equation:

wherecis the diffusivity, and θ =B(σxx+ σyy+ σzz)/3 represents the undrained response of the medium (pore pressure changes due to compression)withBbeing the Skempton’s coefficient,pis the pore pressure, andtis the time. The solution for the pore pressurepis made of two parts: (1) the solution of the homogeneous diffusion equation, and (2) the contribution of the inhomogeneous term for the decoupled approximation (Kalpna and Chander, 2000):

wherepBLis the pressure caused by homogeneous diffusion from the bottom of the load/lake (flow boundary condition), which can be expressed through Green’s functionGDand the weight of the lakeWL, assuming that the medium is initial not pressurized, as

The parameterpCis the pressure caused by diffusion of the undrained response, and it is expressed through Green’s function with the following time evolution for discrete time intervalti:

where θiis computed by solving the Boussinesq’s problem for a surface load changes at timeti.

In order to understand the potential for inducing earthquake,the changes in ΔCFScan be computed as

where Δτ is the shear stress variation in the direction of slip,μ is the coefficient of friction, and Δσnis the normal stress changes (positive for compression) (King and Cocco, 2001). Δτ and Δσnare calculated by rotating the full stress sensor (see Eq. (1)) into the orientation of the fault system,andpis obtained by solving Eq.(2).A positive value of ΔCFScorresponds to a stress perturbation on the fault favoring destabilization.

3.2. Earthquake nucleation

The well-known earthquake generation model introduced by Dieterich(1994)is applied,which accounts for ΔCFSand rate-andstate dependent frictional law as observed in experimental data.The main assumptions of this model are that a large number of potential nucleation sites exist in any volume and that earthquakes are nucleating independently from each other. In this model, the earthquake nucleation rateRdepends on the state variable γ, the constant tectonic background stressing rate ˙S, and the tectonic background seismicity rateraccording to

Fig. 3. ΔCFS values simulated at depths of 2 km, 4 km, and 6 km at the time of impoundment and assuming a receiver orientation as the focal mechanisms in Fig.1.Points represent the projected hypocenters within the corresponding depth interval for the real data (Improta et al., 2017).

Fig.4. Comparison of results for different focal mechanisms:(a)ΔCFS value calculated at depth of 4 km for receiver oriented as the MMFS;(b)ΔCFS value calculated at depth of 4 km for receiver oriented as the EAFS; and (c) Pore pressure changes caused by Boussinesq stress variation as well as diffusion from the reservoir bottom.

The evolution of the state variable is governed by

where dCFSis the infinitesimal variation of ΔCFSwhich is calculated by Eq. (6) but using an effective friction coefficient μeff= μ- α.Here,α is a dimensionless constant accounting for the dependence on normal stress variations, Δσn, which has typical laboratory values in the range of 0.25-0.5.

For an arbitrary stressing history consisting of transient stress changes ΔCFS(t) in addition to the constant tectonic loading/stressing rate ˙S, the evolution of γ can be tracked by considering sufficiently small time steps leading to stress increments of ΔCFS(t)during time interval of Δt. Implementing the stress-step in the center of the time step Δt,the state variable is iterated according to Eq.(9)starting from the background level,that is,γ(0)=1/˙S(Hainzl et al.,2010).Tectonic forces alone would lead to a continuous stress change with constant stressing rate and thus to a constant background seismicity rate. However, the additional reservoir induced stresses lead to temporal variations of the model rate.

Fig. 5. Temporal variations: (a) Pertusillo water level (m) in 2001-2014; and (b) Pore pressure changes and (c) ΔCFS values in 2001-2014 for different values of the intrinsic permeability.

The seismicity rate calculation depends on the model parametersr,Aσ, and. These parameters are not well constrained by independent observations.In particular,σ and˙can vary significantly.The tectonic background seismicity rate (r) could be in principle estimated by declustering the seismicity observed prior to 2001.However,because of missing local high-quality records in the past and problematic declustering procedures,an independent estimate ofris difficult.Thus all three values are used as model parameters that can be calibrated to reproduce observations. To find the best model parameters, the maximum likelihood approach is used to optimize the model fit for the considered time period[ts,te]and in the spatial volume [x1,x2][y1,y2][z1,z2]. Given the model rate of seismicity (R) andNobserved earthquakes occurring at the spacetime pointsti,xi,yi,zi(i=1,…,N),the log-likelihood value(referred to as LL-value hereinafter)is given as a function of the parametersAσ,andrby(Daley and Vere-Jones, 2003)

Using a grid-search forAσ andthe corresponding value ofrwhich maximizes the LL-value can be analytically determined in each case by

A similar approach was used to model the reservoir effect in the Gujarat region, India (Hainzl et al., 2015). The proposed approach could potentially be used to predict the seismicity given a water level variation plan and accounting for some realistic parameters from literature,even if no seismic data are available.In the context,the approach aims at understanding the potential link between water level variation and seismicity.Building a predictive tool is far beyond the goal of the current work.

Fig.6. Results of the maximum LL-search.(a)The dependence of the maximum LL-value as function of μeff;(b)The maximum LL-value as function of the permeability;and(c)For μeff =0.3 and κ=5×10-12 m2,color-coded results in dependence of Aσ and ta=Aσ/˙S,where the colors refer to as the difference to the maximum LL-value and the blue dot marks the position of the maximum.

3.3. Model setup

The simulated period for the case of the Pertusillo Lake started from the reservoir impoundment in October 1962 (see Fig. A1 in the Appendix). The full history of water level followed the data from Stabile et al.(2014),assuming a linear increase in water level for the first 4.5 years since reservoir impoundment. Data from 2001 to 2013 were taken from an online database by a national agency (www.adb.basilicata.it). In a base case simulation, the intrinsic permeability (κ) was assumed to be 5 × 10-13m2, i.e.similar to estimates for the region (Improta et al., 2015). Water density and viscosity,and storativity were assumed to be constant(i.e.ρw=1000 kg/m3;ηw=10-3Pa s;Se= 1/M+α2/(K+4G/3),with Biot modulusM=22 GPa,Biot coefficient α=1,drained bulk modulusK= 25 GPa, and shear modulusG= 11.54 GPa). Such values correspond to a hydraulic diffusivityD= 7 m2/s and Skempton coefficientB= 0.47. A comparison between the semianalytical pressure calculation and a commercial simulator,FLAC3D (ITASCA, 2017), for a simplified load can be found in thesupplementary material (Fig. A2). The use of a semi-analytical approach provides a more flexible tool for quick calculation. The poroelastic stress and pressure variations were calculated in daily steps at grid nodes (xj, yj, zj) with 1.3 km lateral spacing and at three different depths of 2 km, 4 km, and 6 km ΔCFSvalues were calculated for the main focal mechanisms (FM1) of the MMFS as highlighted in Fig. 1. The receivers for ΔCFSare oriented with strike φf=330°,dip δ=60°,and rake λ=-80°,i.e.dipping toward east.For comparison,the ΔCFSvalues for receivers oriented as the EAFS in the oil production region were also calculated (FM2:φf= 120°, δ = 60°, and λ = -80°, dipping toward the Pertusillo Lake). A sensitivity study was performed for the permeability in the range of 5×10-15—5×10-12m2(diffusivityD=0.07-70 m2/s).Records of seismicity prior to local network installation in 2001 are not reliable due to high uncertainties on events location and high magnitude of completeness in the current study (i.e.Mc>2.5). The selected time period corresponds to time betweents= 14,190 d andte= 18,719 d after reservoir impoundment.

Table 1Optimized model parameters. Highlighted in bold is the best fit.

4. Results

A cumulative positive value of ΔCFSis found at the location of the recorded seismicity when accounting for the time of impoundment and assuming an orientation similar to the observed focal mechanisms (FM1, Fig. 3).

The seismicity in different depth intervals has been projected to the depth of 2 km, 4 km, and 6 km for comparison with stress calculation. All seismic events are located in the positive region of ΔCFS,despite the relative small value(below 0.5 bar).The value of ΔCFSslightly decreases with depth,with a value higher than 0.1 bar up to 6 km depth and 15 km away from the center of the Pertusillo Lake. Such value has been often reported as a threshold to induce seismic events (e.g. Cochran et al., 2004).

Fig.7. Comparison of the predictions and observations.(a)Observed M ≥1 earthquake locations(red points)in comparison to the color-coded predicted spatial density of events;(b) Cross-correlation between observed seismicity (Gaussian smoothed with standard deviation of 30 d), reservoir levels(blue) and model rates (red);and (c) The total predicted earthquake rates (cumulated in the whole seismogenic volume) in comparison to the observed earthquake activity in the same region.

It shows that the calculated ΔCFSvalue strongly depends on the orientation of the receiver fault(King and Cocco,2001).Calculation with the same orientation as the EAFS(FM2)also results in positive variation to the north of the lake at the location of the Val d’Agri oil field (orange dot in Fig. 3). Fig. 4 shows the comparison between the ΔCFSvalue calculated for the orientation of the MMFS (FM1;Fig.4a)and of the EAFS(FM2;Fig.4b)at depth of 4 km,respectively.As a rule of thumb,if a normal fault system dips toward the water reservoir (common geological condition for a basin), reactivation will favor these faults,given the similar rock properties.It is worth mentioning that, for both FM1 and FM2, the pore pressure is identical. The pressurization caused by both rock compression(Boussinesq/poroelasticity) and fluid diffusion from the bottom of the lake (boundary) is homogeneous at depth, slightly elongated given the profile of the Pertusillo Lake (Fig. 4c). The pressurized area extends up to a radius of 15 km from the center of the lake.

Fig. 8. Conceptual models of stress changes for both the MMFS (southwest, seismic region) and the EAFS (northwest, oil production region).

Fig.5 shows a comparison of the temporal variation in Coulomb stress changes and pore pressure for different values of permeability.The changes are calculated from the imposed water level of the Pertusillo Lake since impoundment (Fig. 5a). Overall, the pore pressure shows a seasonal water level variation in the order of 10-2bar,and the permeability value affects both the magnitude of the oscillation as well as the time delay with respect to the maximum(minimum)water level(Fig.5b).For a change in 30 m in water level(e.g.around October 2001),the simulated pore pressure varies about 0.03 bar.Interestingly,little changes were observed in terms of ΔCFSwhen the permeability varied.All considered values result in a similar trend and magnitude of variation:up to 0.05 bar for a 30 m water level change.The cross correlation of the pressure variation with the water level changes shows a delay up to 16 d for the case of low permeability,while the ΔCFSis always in phase for all permeability values (i.e. with a cross-correlation lags less than 1 d).

The results for the maximum LL-value are presented in Fig. 6 and the optimized parameter sets are listed in Table 1, together with the resulting correlation coefficients (C) and LL-value. The model best-fit observation is the one with the largest permeability(5 ×10-12m2), although the highest temporal correlation is achieved at a lower permeability(5×10-13m2).Both models feature a quick variation of pressure at depth when the water level changes in the reservoir (Fig. 5).

The maximum LL-value could be of even higher permeability(Fig.6b).However,a larger value of permeability would be unlikely for the depth considered (2—6 km). The model requires a highly sensitive frictional response (smallAσ-value of 1 kPa approximately) to explain the observed variations. In the case that fluids are confined in fault zones (e.g. fracture network), the expected pore pressure variations would be significantly higher and thus theAσ-value would increase accordingly and a larger effective permeability would be possible.

The poroelastic model with FM1 (an orientation as EAFS—FM2 cannot explain the observations and is thus not discussed herein)can well explain the spatial location of the cluster (Fig. 7a). The cross-correlations for the best-fit model (κ = 5 × 10-12m2) are presented in Fig. 7b as a function of the time delay between the model forecast and the seismicity. The correlation coefficients between model rates and seismicity are significantly higher than the correlations between the lake-level variations and seismicity.Finally, the corresponding time-dependent model predictions are shown in Fig. 7c in comparison with the observations.

5. Discussion

The proposed model can well explain the seismicity located at the southwest of the lake,with a ΔCFSvalue up to 0.5 bar.A similar value was estimated for the filling of the Zipingu dam,with a filling of 200 m of water causing ΔCFSvalue between 0.1 bar and 0.5 bar at depths of 10-20 km, depending on the modeling approach (Ge et al., 2009; Lei, 2011). Seismicity has also been reported to associate with pressure changes, for example the case of wastewater injection (e.g. 0.7 bar for Oklahoma — Keranen et al., 2014), hydraulic fractures for seasonal groundwater recharge (0.1 bar at depth of 4.5 km below Mt.Hood)(Saar and Manga,2003),for stress changes during hydraulic fracturing(0.1 bar,1—4 km from injection point at Crooked Lake) (Deng et al., 2016), and even for stress changes associated with Earth tides (as low as 0.01 bar) (Cochran et al., 2004). For the Apennines, a recent study (D’Agostino et al.,2018) was reported on seasonal changes in the rate of seismicity occurring at depths of 3-8 km in the Irpinia region,located about 80 km to the northwest of the Val d’Agri area,with very small stress changes (ΔCFSvalue of approximately 0.1—0.15 bar) caused by intense groundwater recharge of a giant,shallow karst aquifer.

Following the results of the ΔCFSmodel,it can be concluded that not only the MMFS but also the EAFS (Fig. 4) is favorably oriented with respect to water loading for receiving positive stress change to excite seismicity. Stabile et al. (2015) explained the reason for the absence of RIS to the northeast of the Pertusillo Lake,which may be attributed to the effect of the presumed lower permeability as a feature of this area. However, oil extraction data indicate that the northeastern sector is actually characterized by a high permeable fractured carbonate reservoir(κ=10-12-10-13m2),as highlighted by the elevated fluid flow associated production wells and also inferred by migration of the seismic diffusion front in the injectioninduced cluster (Improta et al., 2015). Structural cross-sections,constrained by subsurface exploration data, indicates that the Val d’Agri region features a caprock melange,composed of ductile,low permeable mudstones with high formation pressure and relatively high Poisson’s ratio (Improta et al., 2017). The thickness of this highly deformable caprock formation that overlies the carbonate reservoir strongly varies laterally across the lake, being thinner in the southwestern sector(few hundred meters)and much thicker in the northeastern sector(more than 1500 m)(see Fig.2c)(Catalano et al., 2004). Being highly ductile, the melange layer acts as an effective mechanical decoupling layer(D’Adda et al.,2017)and can accommodate larger strain caused by the weight of the lake. The thick caprock in the northeastern sector could prevent such strain to propagate into the seismogenic depths(2—6 km).Besides,there is an important factor that must be considered to explain the lack of RIS events below the oil field. As highlighted by Improta et al.(2017), the area to the northeast of the lake corresponds to a zone of the oilfield significantly exploited by oil production,causing a depletion of 30 bar in the carbonate reservoir since the late 1990s.Due to the reservoir compartmentalization, pore pressure depletion cannot propagate from the productive zones of the hydrocarbon carbonate reservoir,located to the northeast of the lake,to the southwestern zones where the observed RIS nucleates,outside the oilfield and featuring an aqueous carbonate reservoir.Thus, it is proposed that a combination of thick caprock and fluid production could have hindered the seismicity. Fig. 8 shows a conceptual model to summarize the discussion. The faults located at the southwest of the Pertusillo Lake are subjected to a certain variation of shear stress, normal stress, and pore-pressure by compression and diffusion(red line in Fig.8).Such variation is large enough to reach critical stress for reactivation(black line in Fig.8).The variation induced by the same water level change at the northwest of the Pertusillo Lake could be smaller because of different rock properties(blue line in Fig.8).The depletion caused by oil production could have further shifted the variation toward larger effective shear and normal stresses(Zbinden et al., 2017).

6. Conclusions

An approach was developed by combining poroelastic stress evolution with a rate-and-state earthquake nucleation. The approach was applied for the case of the Pertusillo Lake,Val d’Agri,Italy. This location represents a unique case, where the natural seismicity and reservoir induced seismicity, together with oil production,and wastewater injection reduced seismicity are observed.Overall, a combination of poroelastic stress variation and an earthquake nucleation model well explains the timing of the RIS at the MMFS,southwest of the Pertusillo Lake.The calculated ΔCFSis in the order of 0.5 bar after impoundment,with seasonal variation of 0.05 bar. Such value is at the limit of proposed threshold for triggering seismicity. In terms of temporal evolution, the model features a significantly higher correlation with seismicity compared to water level only.The model would also predict seismicity at the location of the oil production(EAFS—northwest of Pertusillo Lake).Geological conditions(e.g.different lithologies and rheologies)and the depletion at the Val d’Agri oil field explain the absence of seismicity.

The combined approach is obviously simplified, in particular with respect to the geological complexity and rheological and hydrogeological properties. However, the current work represents the first attempt to model the RIS observed at the Pertusillo Lake in Val d’Agri, not only with calculation of static ΔCFSbut also with coupling to rate-and-state model that allows for a better temporal evolution correlation.While previously,the link between RIS at the Pertusillo Lake and water level variation was only based on the temporal correlation, the current results clearly show a physical mechanism explaining the observed seismicity and temporal trend.Finally, the proposed approach represents a tool for a first-order analysis and discrimination of induced and natural seismicities,particularly relevant in case like the Val d’Agri where injectioninduced seismicity as well as RIS and natural seismicity can occur.

Declaration of Competing Interest

The authors wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.

Acknowledgments

This work was jointly funded by a research agreement between the Swiss Seismological Service(SED)and the Istituto Nazionale di Geofisica e Vulcanologia and by the RISE project under the European Union’s Horizon 2020 research and innovation programme(Grant No. 821115). A.P. Rinaldi was also partly financed by a SNSF Ambizione Energy grant (PZENP2 160555). Technical review comments by Antonio Petruccelli at SED are greatly appreciated.Comments from two anonymous reviewers helped improving the manuscript.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://doi.org/10.1016/j.jrmge.2020.04.003.