Finite Element Modelling of Coupled Fluid-Flow and Geomechanical Aspects for the Sustainable Exploitation of Reservoirs: The Case Study of the Cavone Reservoir

The study of the seismogenic mechanical effects induced by oil & gas activities is a socially impacting issue of environmental engineering as well as a challenging task in computational geomechanics. It requires the solution of a coupled problem governed by poroelastic and fluid flow equations in a faulted domain in the presence of in situ stress fields. As a viable alternative to state-of-the-art academical computational models, the present study contributes a simplified methodology based on a commercial Finite Element multiphysics software. The focus is on the evaluation of the link between the oil & gas activities of the Cavone oilfield reservoir, located in North Italy and adjacent to the Mirandola fault, and the recent seismic sequence that struck Emilia in May 2012. An operational coupled fluid-geomechanical procedure is developed where the Cavone reservoir is subjected to the typical in situ stresses, and the nearby Mirandola fault is modelled as an impervious thin layer.


Introduction
Since the early 20th century, almost everywhere, intensive oil-extraction activities have been carried out within tectonically active areas. Seismicity associated with injection of large volumes of fluids has been repetitively documented, for instance, in the Unites States [1], Asia [2] and Europe [3]. Therefore, research highlighting the causal link between earthquakes and oil production are of primary relevance to increase societal safety and reduce the inherent economical losses. Besides being a socially-impacting issue, the study of the seismogenic effects induced by oil & gas activities is a challenging task in environmental engineering and computational geomechanics, as it requires the solution of a coupled problem governed by poroelastic and fluid flow equations in a faulted domain. The present contribution aims to provide a reliable operative methodology based on a commercial multiphysics Finite Element software. As a representative case study, attention is focused on the Cavone reservoir in the presence of the Mirandola fault, hit by the May 2012 earthquake occurred in Emilia, a Northern Italy region.
It has been suggested that most aftershocks of the 2012 Emilia seismic sequence were triggered by previous seismic events [4]. For instance, a seismic event can lead to a slow or rapid change of the stress state and locally trigger earthquakes by static stress transfer [5]. Besides such coseismically generated earthquakes, anthropogenic activities can cause stress changes and seismically activate locally exposed faults. Fluid injection within a reservoir is one such activity [6]: fault slip can occur due to the reduction of the fault strength by shear stress accumulations and fluid pressures variations, to an extent depending on the initial stress state, the elastic moduli of the porous medium and the fault's frictional properties, among others.
Commissioned from the Italian Government to some of the most representative international scientists, the Cavone report [7] is the most extensive reference document for the evaluation of the sismogenic risk inherent to Cavone oilfield activities. A journal paper [8] with a synthesis of these results is also available. The Cavone report assessed stress and pressure variations on the Mirandola fault at the hypocentre of the May 29th aftershock, and investigated the possible role of oilfield injection and production activities on the seismic May 29th event. To this purpose, the reservoir was analyzed as a three-dimensional model where injection and production zones were implemented as point sources sunk in a half-space overlain by an impermeable layer. More recently, Albano et al. [9] compared the stress changes from water injection at the Cavone reservoir and its potential effects on 2012 Emilia earthquake by means of a two-dimensional finite element code.
Taking into account the coupling between fluid flow and geomechanics is key to gaining an insight into the possible interplay between oilfield activities and seismicity [10,11]. Especially exploited for analysing carbon dioxide CO 2 storage [6,12,13], some of the computational tools discussed in the relevant literature result from the coupling of different simulators. For instance, the Tough-Flac code, a simulator designed for modeling coupled multiphase flow and geomechanical problems, combines the TOUGH2 R multi-phase fluid flow simulator [14] with the FLAC3D R geomechanical simulator [15]. In the Cavone report [7,8], the General Purpose Research Simulator (GPRS) [16,17] for flows was coupled with the mechanical simulator PyLith [18,19]. Other prominent references in the field were contributed by Juanes and Jha and collaborators [8,[20][21][22]. For a more exhaustive list of coupled geomechanical models, reference can be made to Jha's Ph.D. thesis [20].
The choice of using a commercial multiphysics software may seem a viable alternative to the use of the aforementioned academical simulation codes. Accordingly, in the present study, the Finite Element (FE) solver Comsol Multiphysics c has been exploited to perform a coupled poroelastic simulation and estimate stress changes and pore pressure perturbations generated by fluid estraction/injection for the Cavone case study.
The original contribution of the present study is the proposal of an operational methodology consisting in a sequence of increasingly complex steps. As a starting point, the adopted simplified geometry is validated by comparing the pore pressure values computed in the reservoir without the fault and in situ stresses. For this purpose, in Section 3, a simplified model of the Cavone reservoir not including the Mirandola fault is constructed with wells modelled as lines [23]. Then, the Coulomb stress is computed in Section 4 considering both the oilfield activities at the Cavone reservoir and the Mirandola fault, the latter being modeled as a fault zone of different permeability. In a first stage, the in situ stress field has not been considered. In a second stage, in Section 5, the Coulomb stress change is determined when the fault, considered either emergent or blind, is subjected to in situ stresses. In particular, the Coulomb stress change is computed by substraction of the stress obtained in the case where the oil injection activities lack, from those computed in the circumstance of oil injection. In all the investigated examples, obtained results are compared with those of the Cavone report [7,8].
The main findings of this study concern the quantitative evaluation of the Coulomb stress changes associated with the injection/extraction activities of the Cavone oilfield, the assessment of the role of the fault, considered either blind or emergent, the importance of considering in situ stresses, and the results' sensitivity to the permeability of the bottom support of the reservoir. The latter finding is particularly addressed here for the first time in the context of the present case study.
The present study is addressed to technicians and scientists who are not necessarily experts in coupled fluid-geomechanical academical softwares, but are put in the position to rapidly assess the possible seismogenic effect induced by oil & gas activities.

Basic Statements and Problem Setting
According to the Coulomb failure criterion, the critical condition is reached when the magnitude of the shear stress computed in the slip direction exceeds the critical value [5,24] where µ s is the static friction coefficient, σ n is the magnitude of the normal stress on the rupture plane and p f is the pore fluid pressure within the fault. Let the slip on a fault occur when the applied stress change dCFF = dτ s + µ s (dσ n + dp f ) overcomes a threshold, where dτ s is positive if evaluated in the slip direction, and the fault normal stress change dσ n is positive for extension. Under undrained conditions, the pore pressure change resulting from a change in stress is [25][26][27] where B is Skempton's coefficient. The stress state induced by earthquakes is governed by the interplay between hydrologic and mechanical properties of the lithosphere, and changes of the pore pressure influence the coseismic stress distribution. The lithosphere response can be explained by adopting the poroelastic constitutive law relating the components of the stress tensor σ ij to those of the strain tensor ε ij [25] where G = E 2(1+ν) and K = 2(1+ν)G 3(1−2ν) are the shear and the drained bulk modulus, respectively, and α B is Biot's coefficient.
After replacement of the constitutive Equation (4) and the compatibility equation ε ij = 1 2 (u i,j + u j,i ) in the equilibrium equation where gravity has been neglected, the set of Navier's partial differential equations is obtained. It can be shown that the fluid mass change dm can be expressed as [24,26] where K u is the undrained bulk modulus. Furthermore, the fluid mass conservation equation for a three-dimensional flow in a porous medium with permeability k and dynamic viscosity µ f , whose velocity field q is ruled by Darcy's law writes Substitution of Equation (7) into Equation (9) leads to the following partial differential equation for The three partial differential Equation (6) and the partial differential Equation (10) form the searched coupled system of equations to be solved. For an extension of the above equations to the case in which gravity is considered, reference can be made to [27].

Preliminary Assessment of the Pressure at the Cavone Oilfield
As a starting point, in the present section, the pressure at Cavone at the time of the May 2012 earthquake has been computed by means of a poroelastic analysis without taking into account the Mirandola fault and the in situ stresses. The computed pressure has been compared with that indicated by the Cavone report [7].
The present simulation relies on the production and injection volumes and time intervals available from the Cavone report [7,28].The oilfield injection and extraction data are illustrated in Section 3.2.

Finite Element Model of the Reservoir
Cavone reservoir is delimited by an anticline which developed above a reverse fault structure. A sketch of the geometry can be drawn from Figure 1. The Cavone structure consists of a fold on a steeply south-dipping reverse fault, the Mirandola fault [7]. The geometry of the Cavone structure has been drawn from that deduced in the Cavone report based on seismic reflection profiles and well data furnished by the operators of the Mirandola concession. No data on the point sources depth or the thickness of the impermeable top layer are available. Two simple reservoir models have been implemented in the FE solver Comsol Multiphysics c : one with two and another one with three layers, reproducing the cases of unsupported and supported reservoir, respectively. In particular, the reservoir has been modelled as the rectangular parallelepiped of 30 km × 8 km × 4.75 km shown in Figure 2. Its thickness was established based on the log and wells data available from Assomineraria [28]. The reservoir thickness is an average value of the thickness of the layers containing oil or oil traces. In a first stage, the block is considered made of two layers, the upper caprock and the reservoir layers, with a thickness of 2.75 km and 2 km, respectively. The reservoir's bottom is in fact confined from an aquifer. Thus, the presence of the bottom support provided by the aquifer is simulated by subdividing the block into three layers: the upper caprock layer of the same thickness as before, the reservoir layer of thickness 1 km, and a lower caprock layer of thickness 2 km, as illustrated in Figure 3.  The mechanical properties of the lower caprock layer are the same as those of the upper one. The model data are in Table 1; whenever available, the same material properties used in the Cavone report [7] have been adopted.
The maximum and minimum element mesh sizes are 1050 m and 45 m for the unsupported and supported reservoir models, respectively.
Rollers have been applied to each external surface except the top surface, which is free.
In the present contribution, for the sake of simplicity, the fluid in the reservoir has been modelled as a single oil phase with density 970 Kg/m 3 , dynamic viscosity 0.4 MPa s, and compressibility 1.3 × 10 −9 Pa −1 .
Colormaps are consistent with the engineering notation where compression stresses are negative and tension stresses positive. This is the opposite of the notation used in geology.
Young's modulus and Poisson coefficient of the rock were assumed being E = 55 GPa and ν = 0.28. The fluid in the bulk is assumed to be water with density 1000 Kg/m 3 , viscosity 0.001 Pa s, bulk modulus 2.2 × 10 −9 Pa, and compressibility 45.8 × 10 −11 Pa −1 . We assumed the pores being fully saturated and set Skempton's coefficient B = 1.

Production Activity
Wells have been located at the World-Geodetic-System-1984 coordinates furnished by Assomineraria [28]. To this purpose, injection and production volume data have been gathered, and transformed into mass fluxes' values suited to modeling the wells as lines.
In the Cavone report [7], the injection and production zones were simulated as point sources in a half-space overlain by an impermeable layer. However, the point sources' depth and the thickness of the impermeable layer were not given. In a first attempt, wells were introduced in COMSOL Multiphysics as small 3D cylinders in a large layered subsurface domain subjected to suitable boundary conditions. This was not very efficient as it required meshing wells with a radius of 0.2 m, namely very small with respect to the surrounding domain. Therefore, each well has been replaced by a source line along its axis [23]. To this purpose, the "well boundary-condition" option offered by COMSOL Multiphysics 5.3 has been exploited. It is based on the Thiem's solution [29] for homogeneous media written for vanishing well radius and subsequent extension to heterogeneous aquifers through definition of an equivalent conductivity [23]. With the well boundary condition, cylinders are replaced by edges, thus reducing the computational burden associated with fine meshing. This feature was shown to provide accurate solutions compared to modelling the wells as three-dimensional cylinders. Additional data such as the wells' depth have been deduced from the ICHESE report [30]: the unknown depth values have been assumed equal to an hypothetical but realistic value of 3 km.
An average flux rate has been implemented for each well. All the injection-production rates have been entered in the model as the time-dependent average flux rates shown in Figure 4; mass fluxes produced and injected have negative and positive signs, respectively.
A time dependent simulation has been run over a time interval of 12,510 days, from March 1980 to June 2014.

Results
According to the available data [7], values of the reservoir permeability k res ranging from 0.3 mD to 1 mD were used. Major pressure change were however associated with smaller values of k res .
Thus, we choose to show in the current section the results obtained for the lowest data-based value of k res = 0.3 mD.
The results of the poroelastic analysis are shown for both supported and unsupported reservoir models. More specifically, Figure 5 displays the pressure on the top surface of the reservoir, the shaded surface in Figure 6. In both cases, the major pressure's increase has been detected near well C014, the main injector of the field, while the pressure's decreases have been measured near wells C007, C009, and C017. The prominent influence on pressure changes induced by the production activity of well C014 was also fully recognized and investigated by Albano et al. [9]. Though we represent wells as production-injection lines and not as point sources, the obtained pressure values fit with those of the Cavone report [7].   Figure 7 displays the vertical-displacement pattern obtained for both the two-layer and the three-layer models. In fact, the simulated injection-production activities induce subsidence, highlighted by the cold colormap, especially nearby the oilfield, while, expectedly, the entity of the subsidence-induce vertical displacements is sensitive to the layers' mechanical properties.

Modeling the Mirandola Fault and the Cavone Oilfield
The poroelastic response of a multi-layer block including both the Mirandola fault and the Cavone oilfield has been simulated. Two different subcases have been considered: with and without in situ stresses.

Model Geometry
A three-dimensional block domain whose dimensions are 30 km × 30 km × 20 km has been modelled, and the model properties of Table 1 have been assigned. The FE model consists of the three-dimensional block in Figure 8 subdivided into five different layers: the upper, the upper caprock, the reservoir, the lower caprock, and, finally, the lower layer.
Different permeability values for each of the layers displayed in Figure 8b have been adopted. The upper caprock layer has very low permeability to avoid communication of the reservoir with the layers above [7], while a lower permeability ranging from 0.1 to 0.001 mD has been assigned to the bottom caprock layer. The reservoir permeability has been set equal to 1 mD for the sake of establishing a comparison with the coupled geomechanical analysis [7], where this value was used. In the Cavone reservoir, the fault is a sealing stratum for the reservoir. Therefore, the permeability value of the fault is comparatively low, equal to 1 × 10 −9 mD. This approach is conservative, as low fault permeability favors slip [31].
Considering that typically the friction coefficient µ s varies in the interval 0.6-0.85, and that the choice of low µ s values is conservative [31], µ s has been assumed equal to 0.6 [32]. Because no information about Skempton's and Biot's coefficients were available, the values of 1 and 0.8, respectively, were set [9].
The Mirandola fault has a curved shape, concave to the South, as shown in Figure 1, where the coordinate system is such that the x-axis is easting, the y-axis is northing, and the z-axis is elevation in meters. However, the Mirandola fault has been implemented as a plane to simplify the analysis. When faults are modelled as finite thickness layers with mechanical properties softer than those of the surrounding domain, after imposition of the reverse fault regime by application of in situ stresses, spurious high stress peaks arise near the fault zone and propagate to other parts of the domain. Hence, it has been decided to model the fault zone as a comparatively thin layer of thickness 100 m with the same mechanical properties of the surroundings except for the permeability.
Though blind, the fault has been first modelled as emergent according to the Cavone report [7]. Eventually, it has been modelled as blind to see how results changed.

Production Activities
We drive the dynamic simulation with well rates as described above. A production-rate history more accurate than that of Section 3 has been implemented for the major producers and all three injectors according to the Cavone report [7]. The production-rate histories of the other wells are not available from the literature. In particular, the rates history of five producer wells and three injector wells out of nineteen wells are illustrated in Figures 9 and 10, respectively. For the other wells, we exploit the same rate history used in Section 3.
The time interval chosen for the poroelastic simulation covers 11,994 days, from the 1 March 1980 to 31 December 2012. For the time dependent study, a time step of 250 days was used, which is quite small compared to the 11,994 days, to account for the rate fluctuation shown in Figures 9 and 10.

Mirandola Fault and Cavone Oilfield without In Situ Stress
In a first stage, the changes in effective normal stress, Coulomb Failure Function and pressure have been evaluated in the presence of the oilfield activities without considering the in situ stress state; this extends the case treated in Section 3 by adding the Mirandola fault. These results were obtained on setting the permeability of the lower caprock layer equal to 0.1 mD. Figure 11 illustrates the color-maps of the changes of shear (a), normal (b) and Coulomb (c) stresses. A very similar colormap, hereby reproduced in Figure 12, was obtained in the Cavone report [7] assuming the reservoir in communication with the underlying aquifer. It is worth noting that the region of influence extends above and below the reservoir: the reservoir expansion due to injection near well C014 increases the effective normal compression above and below the reservoir, in blue on the contour map. Where production is the main activity, the consequent reservoir contraction leads to a drop in the effective compressive normal traction, highlighted in red in the colourmap. In Figure 11c, in the production areas, redthe Coulomb Failure Function change dCFF (2) is positive below the reservoir and negative above it, while colours invert in the injection areas, especially near the C014 well. As suggested in the Cavone report [7], contraction leads to up-dip and down-dip shear traction on the layers below and above the reservoir, respectively. However, the C014 injector induces expansion of the reservoir, and, thus, in its surroundings, the changes are opposite. Furthermore, reservoir contraction leads to a drop in the effective compressive normal traction (red color) in the reservoir depth interval as well as below the reservoir, except near C014 where injection-induced expansion causes an increase in the effective normal compression above and below the reservoir (blue in Figure).   Figure 13, the present simulation does not furnish values of dCFF comparable with those of the Cavone report. Hence, a more refined analysis accounting for the in situ stress state has been performed, as discussed in Section 5.

Coulomb Stress Changes
For each value of permeability of the lower caprock layer, two different models were set up. The former, called model (0), can be regarded as an "initial" condition; it takes into account the in situ stress and neglects the oilfield activities. The latter model, coined model (1), takes into account both in situ stress and wells activities. Assuming that application of the effects superposition is allowed, the changes of normal, shear and Coulomb stresses have been computed by subtracting the stress values obtained for model (0) from those established in model (1).

In Situ Stress Field
In situ stresses are prescribed to be lithostatic in the vertical direction, twice the lithostatic in north-south, and 1.5 times the lithostatic in the east-west direction. The lithostatic gradient is calculated with the bulk density assumed to be 2600 kg/m 3 . Rollers were applied to each external surface of the three-dimensional domain except the top surface, which is traction-free, and the northern face, on which an external load equal to the in situ stress in the y-direction was applied [7]. All boundaries are no-flow boundaries. Figure 14 displays the adopted boundary conditions. Let σ 1 , σ 2 and σ 3 denote the principal stresses in order of descending magnitude. The typical in situ stress field of Mirandola fault region for z starting from the bottom is [7] where, according to the Andersonian notation for reverse fault conditions of Figure 15, σ v is the vertical stress, while σ H and σ h are the maximum and minimum horizontal stresses, respectively. We have verified on different FE software besides Comsol that, to get corrected results, the in situ horizontal stresses have to be imposed as internal stresses, and the in situ vertical stress should not be prescribed as internal stresses but through the specific gravity tool of the adopted software.

Variable Lower Caprock Permeability
The results are plotted at points located within the hanging-wall block, where all the Cavone wells are placed. Stresses have been plotted close to, but not exactly at, the hypocentre, namely 1 km in the negative y-direction from the hypocentre, in order to avoid incurring stress oscillations on the fault. Figure 16 illustrates the changes in Coulomb failure function dCFF (a) and effective stress dσ n (b), respectively. It is noteworthy that, as the pore pressure increases, the Coulomb Failure Function and the effective normal stress increase. The results are plotted for values of permeability of the lower caprock layer ranging from 0.1 mD to 0.001 mD. Remarkably, the displayed pressure changes are displayed following the Cavone report sign convention. In particular, as the permeability of the lower caprock decreases, the pore pressure increases. Simultaneously, the Coulomb Failure Function CFF increases, favoring fault reactivation. Interestingly, the maximum dCFF attained at the end of the simulation, namely when the 2012 earthquakes occurred, is around 0.03 bar.  Figure 16a, namely the dCFF that can trigger earthquakes, are quantitatively similar to those of the Cavone report [7].

Results at Different Distances from the Hypocentre
The spatial variability of the stresses has been checked. In particular, at the locations indicated in Figure 17 having different distances from the hypocentre, the results have been plotted along the negative y-axis in Figure 18, and keeping a constant value of k lower equal to 0.1 mD. It can be drawn that the Coulomb stress change dCFF decreases for increasing distances from the fault on the hanging-wall block.   Figure 17).
To check how results change if we move towards the oilfield, the locations indicated in Figure 19 have been selected, where points from A to E gradually move from the hypocentre westward and then towards the Cavone well C014. The corresponding results have been plotted in Figure 20. As can be inferred from Figure 20a, the Coulomb stress change dCFF below the reservoir decreases when approaching well C014. This conclusion is consistent with those established in the Cavone report [7] and displayed in Figure 12. Once again, the positive values of dCCF, namely those that are assumed to induce seismogenic effects, are consistent with those reported in the Cavone report. On the other hand, the computed negative peaks of dCFF are smaller and are delayed with respect to those of the Cavone report, the difference decreasing when approaching the wells. The negative peaks occur at the peak of the injector well C014, as it can be drawn from Figure 10b. Hence, the present model is particularly sensitive to the activity injection of well C014. Note that, here, the wells have been modelled as line sources and not as point sources as in the Cavone report [7].  Figure 20. Emergent Mirandola fault and Cavone oilfield modeling: results at different points from the hypocentre for permeability of the lower caprock layer k lower = 0.1 mD; Coulomb stress change dCFF (bar) (a) and effective normal stress change dσ n (bar) (b) are plotted for different distances from the hypocentre, on the hanging-wall block (see Figure 19).

Blind Fault: Results at the Hypocentre for Variable Lower-Caprock's Permeability
As shown in Figure 1, according to Cavone seismic profiles, the Mirandola fault does not reach the Earth surface. The authors of the Cavone report [7] considered the fault cutting the entire domain only to facilitate the geometry-mesh building. During the simulation, the upper part of the fault was locked in order to avoid effects on the results. In the present study, to get a more realistic fault representation, the additional model of the blind fault of Figure 21 has been implemented. The material properties of the model and the boundary conditions are the same as in the case with the emergent fault. The same parametric analysis of the emergent fault case has been repeated. Remarkably, the results of Figure 22 are quite similar to those of the previous cases of non-emergent fault, except that the maximum dCFF slightly decreases, from 0.03 bar to 0.025 bar, suggesting that the fault emergence might not significantly affect the seismicity induced by the Cavone oilfield.

Discussion
The minimum stress change that is often recognised to potentially trigger seismicity by static stress transfer is around 0.1 bar [33]. For the Emilia seismic sequence, Pezzo et al. [34] computed a Coulomb stress change induced by the May 20th earthquake on the May 29th fault plane of about 6 bar, and suggested that the May 29th aftershock can be considered as the result of static stress changes from the earlier event.
The present results have shown that the positive change in Coulomb stress near the hypocentre of May 29th as a result of injection activities is very small compared to the Coulomb stress changes calculated at the May 29th hypocentre from the May 20th mainshock. In particular, the positive Coulomb stress change presently computed does not exceed a value of 0.03 bar, against the 0.02 bar established by the authors of the Cavone report [7], anyway, a value largely smaller than 6 bar.
It is commonly expected that earthquakes triggered by fluid injection occur in a radius of few kilometers around the injection well with the highest pressure. The obtained results suggest that fluid injection in the Cavone well C014 may have not caused earthquakes migration, even if injection activities have affected the Cavone field for more than 20 years. Fluid injection at well C014 generated a significant positive Coulomb stress change [9]. However, this change is almost localised in the vicinity of the well C014, namely sufficiently far from the May 29th hypocentre location. The same conclusion was also drawn by Albano et al. [9].
The present simulations have highlighted the role played by the choice of the permeability of the lower caprock layer on the determination of the maximum dCFF reached at the May 29th hypocenter. Remarkably, this conclusion is linked to the idea illustrated in the Cavone report that assuming or not the reservoir as supported by the bottom aquifer is key to a realistic evaluation of the stress field.
The possible emergence of the fault has been shown to slightly smooth the dCFF increments at the hypocenter.
The present study does not answer the question of whether pressure changes may have indirectly altered the time of occurrence of the 29 May 2012 event. Assuming that the main shocks of May 2012 had a tectonic origin, and the May 20th earthquake triggered the May 29th by Coulomb stress transfer, it should be established whether anthropogenic activity and its consequent perturbation of the stress or fluid pressure might have advanced the time of the second seismic event near the Cavone oilfield.
Finally, the FE poroelastic simulation can certainly be further improved, for instance, by simulating the real thrust fold conformation, and modelling the fluid in the reservoir as a multi-phase fluid [35], following, for instance Jha's high-resolution numerical scheme for the stable simulation of the viscous fingering process [20] and mixing process in viscously unstable laminar flows. However, it should be said that these limitations have not jeopardised the possibility of computing reliable values of potentially seismogenic Coulomb stress changes.

Conclusions
Proceeding through subsequent, increasingly complex problems, the present study has illustrated a simplified, coupled fluid-flow-poroelastic analysis for evaluating the seismogenic potential of oil extraction activities in the present case-study. The main finding is that the Coulomb stress change associated with Cavone oilfield activities is two orders of magnitude smaller than that computed on the Mirandola fault in the area where the May 29th event nucleated for the May 20th earthquake. Hence, the present analysis suggests that the seismic events of May 2012 seem not to have been directly triggered by the activity at the Cavone oilfield. Funding: The support of the FIR 2017 fund of the University of Ferrara and the financial support of the Agenzia Regionale per la Protezione Civile-Regione Emilia Romagna is gratefully acknowledged.

Conflicts of Interest:
The authors declare no conflict of interest.