Dynamic Fluid Flow and Geomechanical Coupling to Assess the CO 2 Storage Integrity in Faulted Structures

The SiteChar research on the Southern Adriatic Sea site focused on the investigation of the geomechanical and hydrodynamic behaviour of the storage complex in the case of CO2 injection in a reservoir consisting of fractured carbonate formations. Special attention was paid to the effects that natural faults and fractures might have on CO2 migration, and the effects that injection might have on the stability of faults. This assessment was originally performed via a hydro-geomechanical oneway coupling which relies on an adequate representation of faults in the model, allowing one to simulate fluid flow along the fault plane and inside faults as well as evolution of the stress state due to CO2 injection. The geological model was populated with petrophysical and geomechanical parameters derived either from laboratory measurements performed on samples from a reservoir analogue, or published literature. Since only sparse data were available, various scenarios were simulated to take into account the uncertainties in the fluid flow and geomechanical properties of the model: the different state of faults ( i.e., open or closed) and various in situ stress state, commonly named geostatic stresses as the earth’s crust deformation is assumed to be slow regarding the shortterm study. Various fluid flow parameters were also considered, although only one set of petrophysical data corresponding to the most realistic ones is considered here. Faults modeled as volumetric elements behave as flow pathways for fluids when they are conductive. The injected CO2 migrates inside and through the Rovesti fault, which is located near the injection well. The fluid flow also induces overpressure in the faults. The overpressure in the Rovesti fault reaches 2.2 MPa while it reaches 4.4 MPa at the bottom hole of the injector. Extending to about 30 km, the pore pressure field reaches the Gondola fault located at 15 km from the injection zone but the overpressure does not exceed 0.1 MPa at such a distance from the injection well. Using this overpressure as loading in the geomechanical model allows one to compute the effective stress variation in the whole geological model. The total effective stress is then computed by adding an estimation of the regional stress. Post-processing is performed to derive the likely damage of the faults according to the MohrCoulomb criterion. The results are illustrated on the Rovesti fault, which is located near the injection well and consequently the most likely to be reactivated. On the basis of available data, for all the modeled scenarios (various initial stress regimes, closed or open fault), no fault damage is observed, as the stress state stays below the Mohr-Coulomb criteria. Résumé—Couplage des modélisations hydrodynamique et géomécanique pour évaluer l’intégrité d’un stockage de CO2 dans des structures faillées— Les travaux de recherche conduits dans le projet SiteChar sur le site situé en mer Adriatique sud se sont concentrés sur le comportement géomécanique des formations carbonates fracturées de ce site dans un contexte d’une injection de CO2. Une attention particulière est portée sur les effets que les failles et fractures naturelles peuvent avoir sur la migration du Oil & Gas Science and Technology – Rev. IFP Energies nouvelles, Vol. 70 (2015), No. 4, pp. 729-751 A. Baroni et al., published by IFP Energies nouvelles, 2015 DOI: 10.2516/ogst/2015017 This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. CO2 et inversement les effets que l’injection de CO2 peut avoir sur la stabilité des failles. Cette étude a été réalisée via un couplage hydro-géomécanique « one-way » reposant sur une représentation appropriée des failles qui permet de simuler les écoulements de fluide le long du plan de faille et à l’intérieur de celle-ci ainsi que l’évolution du champ de contrainte résultant de l’injection de CO2. Le modèle géologique a été habillé avec des propriétés pétrophysiques et géomécaniques issues de mesures de laboratoire réalisées sur des échantillons provenant d’un analogue ou sur base d’informations provenant de la littérature. Au vu du peu de données disponibles, différents scénarios ont été simulés pour appréhender les incertitudes sur les propriétés d’écoulement et les propriétés géomécaniques du modèle. Différents comportements hydrauliques de la faille (ouverte ou fermée) et différents états de contrainte in situ, communément appelée contrainte géostatique puisque la déformation de la croûte terrestre est supposée faible sur l’échelle de temps de l’étude, ont donc été simulés. Différents paramètres d’écoulement ont aussi été considérés même si un seul jeu de données pétrophysiques correspondant aux données les plus réalistes est ici présenté. Les failles modélisées comme des éléments volumiques, lorsqu’elles sont conductrices se comportent comme des chemins d’écoulement pour les fluides. Le CO2 injecté migre ainsi à l’intérieur et à travers la faille Rovesti qui est située proche du puits d’injection. Les écoulements de fluides induisent aussi une surpression dans les failles. La surpression dans la faille Rovesti atteint 2.2 MPa alors qu’elle atteint 4.4 MPa au fond du puits d’injection. S’étendant sur une trentaine de kilomètres, le champ de pression de pore atteint la faille Gondola située à 15 km de l’injection mais la surpression ne dépasse pas 0.1 MPa à une telle distance du puits. L’utilisation de cette surpression comme chargement dans le modèle géomécanique permet de calculer la variation de la contrainte effective dans l’ensemble du modèle géologique. La contrainte effective totale est alors calculée en rajoutant une estimation de la contrainte régionale. Un post-traitement est réalisé pour déduire l’endommagement éventuel des failles en utilisant un critère de Mohr-Coulomb. Les résultats sont illustrés sur la faille Rovesti, dont la réactivation est la plus probable puisque la faille est située près du point d’injection. À partir des données disponibles et pour tous les scénarios simulés (différents états de contraintes initiaux, failles ouvertes ou fermées), on n’observe pas d’endommagement des failles, l’état de contrainte, pendant et après l’injection, restant sous le critère de Mohr-Coulomb.


INTRODUCTION
The SiteChar project has refined the methodology for characterisation of geological storage of carbon dioxide.Five potential CO 2 storage sites across Europe were selected as test sites for the research.One of these sites is a carbonate deep saline aquifer located offshore in the Southern Adriatic Sea (Fig. 1) that might represent an opportunity to launch the first CCS project in Italy, taking advantage of the vicinity of the storage area to one of the major Italian power plants where the energy company Enel started a pilot plant for CO 2 capture in April 2010.
CO 2 injection into deep saline aquifer increases the pore pressure in porous and fractured formations and changes the state of stress in the reservoir as well as in the surrounding rock.As a consequence, natural fracture systems and faults might be reactivated and new fractures might be created.These fractures and faults can thus act as migration paths for CO 2 or fluid and lead to loss of containment of injected CO 2 .Other effects might also occur: deformation of the storage reservoir (e.g.expansion) and the surrounding rock and ground surface movement (e.g.uplift) generated by stress changes, and micro-seismic events induced by sudden slip on pre-existing faults; they are not addressed in this paper.
Modeling these kinds of phenomena calls for coupled fluid flow and geomechanical simulations as frequently used in the oil and gas domain (Guy et al., 2011(Guy et al., , 2012;;Longuemare et al., 2002).Three methods can be considered to understand the geomechanical effects: the explicit sequential method, the iterative sequential method and the fully coupled method.In the fully coupled method, the equations of both geomechanical and hydrodynamic problems are solved simultaneously within the same simulator thus leading to consistent solutions.However most often current fully coupled simulators integrating geomechanics in fluid flow models rely on simplifications in each of these two domains.The sequentially coupled method presents the advantage of allowing the use of specialised mechanical and reservoir simulators that are coupled in an external way; advanced functionalities available in both software can thus be used.Another intrinsic advantage of sequential methods is that both methods of modeling can use different time discretisation.The main drawback is that most often dynamic fluid flow and geomechanical simulators do not share the same meshes, thus hampering consistent hydromechanical modeling of faults behaviour.In the fluid flow modeling, the fault is usually represented by a plane (Longuemare et al., 2002;Al-Busafi et al., 2005).Therefore, fluids can flow across the fault but not inside.Besides, the transverse flux between the two cells in contact with the fault is modeled by a transmissivity coefficient (Manzocchi et al., 1999;Cappa and Rutqvist, 2011).In this study, the main objective is to simulate the evolution of the stresses acting on the faults while fluids flow inside the faults.For that purpose, a fault is modeled by a volumetric element.
This paper presents (Sect.2) a one-way methodology that couples two simulators: a fluid flow model (1) , and a geomechanical one (2) .Fluid flow and geomechanical simulations are actually performed on the same grid, thus allowing correct modeling of the fluid flow in the reservoir, cap rock and faults.Correct and coherent estimates of the pore pressure variations in the faults induced by CO 2 injection in the reservoir are thus loaded in the geomechanical simulator to update the normal stresses acting on the fault planes.The one-way approach is sufficient for our objective, as we state that geomechanical effects do not influence fluid flow parameters, as long as the stress state is below the damage criterion.If damage is reached, following computation is no longer valid.The Mohr-Coulomb criterion is used in this paper as a damage criterion; if the stress state is over the criterion value, the material is assumed to have encountered irreversible deformation, but not necessarily failure.
The approach is applied to injection of CO 2 in the Southern Adriatic Sea site, which presents several faults likely to be reactivated.Fluid flow simulation of CO 2 injection is first presented (Sect.3) with specific emphasis on the extent of the CO 2 plume and the overpressure footprint.The overpressure term is used to refer to the initial pressure, assumed to be the hydrostatic pressure.The resulting pressure changes are then used as loading for the geomechanical modeling and interpreted in terms of potential fault damage (Sect.4).Conclusions are derived regarding the geomechanical stability of the fault.

GEOLOGICAL SETTING
The geological characterisation of the Southern Adriatic Sea (Volpi et al., 2015) led to three sites potentially suitable for CO 2 storage: the Grazia, Rovesti and Grifone structures (Fig. 2).Rovesti contains oil and gas, whereas the Grifone and Grazia reservoirs are represented by a pure saline formation.Among these three structures, the Grazia anticline was selected as the storage site for the hydrodynamical simulations because of its proximity to the potential CO 2 emitter (namely the Brindisi power plant), the water depth (i.e., around 180 m, which is the shallowest of the three sites) and the absence of interference with other usage, Rovesti being under permit for hydrocarbon exploration.
The 3D geological model is described in Volpi et al. (2015).It is composed of the following geological elements (Fig. 2): the Seabed and the Top of the Messinian sediments, both considered as unconformities; the Top and the Bottom of the Scaglia Formation, considered as horizons; the Jolly Fault, considered as the main fault of the structural system; the Gondola Fault; the Rovesti Lineament; the Grifone structure; the Rosaria and Cigno faults; and four main stratigraphic layers: -Plio-Pleistocene: from the "Seabed" to the "Top Messinian"; -Oligo-Miocene) (cap rock): from the "Top Messinian" to the "Top Scaglia"; -Scaglia (reservoir): from the "Top Scaglia" to the "Base Scaglia"; -Pre-Scaglia: below the reservoir layer.
For the three upper layers (Plio-Pleistocene, cap rock and reservoir), a further internal partition was performed to increase the vertical resolution of the grid.
The 3D volumetric grid contains cells of 500 m 9 500 m lateral extension which allows a very good match with the geological surfaces in the structural model.The final 3D grid model has about 2.4 million cells (356 cells in x direction, 300 cells in y direction, and 22 cells in z direction).

Facies Attribution
The 3D geological model is composed of four geological formations as described in Table 1.Each formation contains one or several facies among these five facies types: marl, shale, dolomite, basinal carbonate and platform carbonate.The Plio-Pleistocene formation is only composed of shale.It is the most sealing part of the cap rock.One additional facies per numerical layer is introduced in order to model its porosity evolution with depth.
Finally, a facies corresponding to fault rock is added in order to simulate the fluid flow and the geomechanical behaviour of the faults.

METHODOLOGY FOR COUPLED FLUID FLOW AND GEOMECHANICAL MODELING
The objective is to model the potential geomechanical deformations which result from the pore pressure changes in the reservoir and the cap rock induced by CO 2 injection.The one-way coupling between fluid flow and geomechanical modeling relies on the following process: accurate fluid flow modeling firstly computes the pore pressure changes induced by CO 2 injection and migration.These changes are secondly used to load the geomechanical simulator and derive the distribution of stresses.Reciprocally, geomechanical deformations might change the petrophysical properties, i.e., porosity and permeability of the rock, and thus impact the fluid flow.Changes in porosity and permeability due to geomechanical modeling should thus be taken into account in the fluid flow modeling.This issue, which calls for a two-way coupling, is not addressed here.

Workflow for the One-Way Coupling between Fluid Flow and Geomechanical Simulators
Figure 3 presents the workflow developed to simulate the geomechanical behaviour of a faulted medium.The pore pressure evolution is computed by the fluid flow simulator PumaFlow TM .Pore pressure distributions at given time steps are transferred to the geomechanical simulator Abaqus FEA, that computes the resulting stress and strain changes.It is a sequential process that relies on six main steps: -Step 1: Construction of the geological model grid in conformity with both PumaFlow TM(3) (fluid flow simulator using finite volume) and Abaqus FEA (geomechanical simulator using finite elements) requirements.Production of an "Eclipse format (4) " grid.
Usually the construction of a geological model does not take into account the specific requirements of both geomechanical and fluid flow modeling (Sect.2.2).In particular, fault cells are not always on the same side of the fault surface, which hampers the correct handling of the fault throw in the fluid flow/geomechanical models.A new list of fault cells adapted to these requirements has thus to be provided, either by the geomechanical model or by an independent process (Sect.2.2.1).-Step 2: Conversion of the grid from "Eclipse format" into a finite element model for the geomechanical model.The gridding of faults requires specific attention.Two facing groups of elements surrounding each fault are introduced to give a volume to the faults.These two facing groups define each side of the fault surface.In addition, one side of the fault is attributed with cohesive elements that are tied on the other fault side (Fig. 4).These cohesive elements are finite elements with a surfacic behaviour, that allow one to attribute a small width to According to the selected CO 2 injection scenario, the fluid flow simulator computes the fluid saturation and the pore pressure increase due to the CO 2 injection and migration at different time steps.-Step 5: Geomechanical simulations.
The pore pressure grid resulting from the fluid flow simulations is interpolated from cell centres to gridpoints as required by the geomechanical simulators.This is done at different simulation times in order to see the evolution of the geomechanical state of the fault.The resulting stresses are then derived.-Step 6: Post-processing to estimate the fault damage.
The resulting stresses are analysed in terms of potential fault damage.The initial stress state is added to the stress changes computed in Step 5 and the distance to the Mohr-Coulomb criterion in the stress space is calculated for the damage occurrence.Deploying this workflow raises a number of technical issues that are described in Section 2.2.

Fault Representation Compatible with Both Fluid Flow and Geomechanical Modeling
The main difficulty when applying a 'one-way' coupling comes from the geological gridding that does not provide a correct representation of faults and other subsurface discontinuities.Most often grids for fluid flow simulation differ from geomechanical meshes.A flow simulation grid is typically a corner-point grid generally restricted to the reservoir with faults represented by stair stepping (a polygonal surface, which does not really fit with the real slope of the fault).Geomechanical meshes are usually finite element grids that take into account the over-, under-and side burden and that must handle faults modeled with a polygonal surface.This surface has to be as close as possible to the real surface to have the right projection of the stress tensor on it.An important challenge is thus to build the same grid of the reservoir, over-, under-and side-burden that will allow one to run both the fluid flow and geomechanical simulations.Faults are actually complex elements that require various representations according to the investigated property and the software to be used (Fig. 4): in the geomodel, the fault is represented by a surface which is the interface between two facing groups of elements (Fig. 4a); in the fluid flow model, the fault is represented by adjacent volume elements located at one side of the fault surface (Fig. 4b); in the geomechanical model, the fault is represented by cohesive elements with geometric initial null width tied to the incoherent side of the fault (Fig. 4c).

Fault Cells
The faults are modeled based on the horizontal curved lineaments (Fig. 2) and vertically along lines (Fig. 4).
For the fluid flow simulation, the cells associated with the fault are all those located along the a priori defined side of the fault surface (Fig. 5).Due to the horizontal curved shape of some faults and the grid orientation in these areas, the connectivity between cells has to be checked: some fault cells have to be added to avoid gridding artefacts and enable fluid to flow along the whole fault, as illustrated in Figure 6.
In this study, the porous rock in the fault is assumed to be associated with a single facies: the fault rock.Its petrophysical and geomechanical properties are detailed in Tables 2, 3 and 4.

Equivalent Fault Properties
The model mesh is in general not sufficiently refined in the fault zone to correctly model the small thickness of the fault.In this study, the size of the fault cell is about 80 m.An average real fault thickness was estimated from Figure 7, which shows several fault thickness data versus the fault throw, also called displacement, collected in Wibberley et al. (2008).For an average displacement of 1 000 m along the fault plane (which is the average throw of faults in the geological model), the fault thickness is of the order of 10 m.
To mimic fluid flow in a 10-m-thick fault with cells of 80 m, an equivalent porosity and an equivalent permeability are introduced for fault facies.Porosity / is thus replaced by equivalent porosity / eq such that: where e is the real fault thickness and Á is the fault cell size.Fault thickness versus displacement (Wibberley et al., 2008).
Regarding the permeability, conservation of the flux per surface unit in the three directions is assumed.Therefore, only the flux through the fault is affected by the thickness and then only the permeability perpendicular to the fault surface has to be taken into account in the equivalent permeability which is thus defined as: where K is the permeability (m 2 ), e is the real fault thickness (m) and Á is the fault cell size (m).Permeability in the fault plane is considered homogeneous.However, this equivalent permeability model could be improved in future work by integrating a SGR (Shale Gouge Ratio) computation along the faults (Manzocchi et al., 1999) to compute a more realistic permeability.

Fault Modeling in the Geomechanical Model
For geomechanical purposes, the fault is considered as a surface, the "fault surface", which is meshed using cohesive elements.Cohesive elements allow one to model the behaviour of fractures considering fracture formation as a gradual phenomenon in which the separation of the surfaces in the crack takes place across an extended cohesive zone and is resisted by surfacic stiffness.A surfacic energy allows one to take into account the energy to open the fracture and the behaviour law in Figure 8 models the behaviour of the fracture zone.Even if with this kind of element the geometrical width is null, it can be considered as 10 metres for the geomechanical behaviour according to the value of parameter e (Fig. 8).

Model Gridding
In order to simplify the data exchange between the fluid flow model and the geomechanical model, a common refined grid was built.
The Tartan refinement (Fig. 9) was chosen for the horizontal refinement of the CO 2 injection area.It consists of progressively refining all the cells in the direction perpendicular to the fault.Three refinement levels were used in practice in the Grazia model to gradually reduce the horizontal cell size from 1 000 m outside the refinement zone to 80 m in the injection zone.Such a refinement reduces numerical diffusion and thus allows for a better simulation of the pressure gradients and the CO 2 plume migration.
Vertical refinement was performed by dividing the Plio-Pleistocene formation into four grid layers and the Oligo-Miocene into six.Moreover, to avoid artefacts at the bottom boundary of the computational domain, the Pre-Scaglia formation was associated with a thickness that is here fully arbitrary and that will not influence the results.This arbitrary thickness is taken to be 500 m according to the estimation of the thickness of the Pre-Scaglia formation provided by OGS (Istituto Nazionale di Oceanografia e di Geofisica Sperimentale) geologists of between 116 m and 732 m.This arbitrary layer was divided into 10 numerical layers but only the shallowest one is considered in the fluid flow model.A no-flux boundary condition is assumed at the bottom of the fluid flow model.
Figure 10 shows the final mesh that comprises 2.1 millions cells.

Remaining Issues
A first limitation of the approach is related to the meshing strategy, which makes all discontinuity surfaces discretised through a set of non-horizontal lines on which cohesive elements are tied.As a consequence, the real vertical curvature of the fault cannot be correctly described (Fig. 10 right).This has a major impact on the stress state, since the surface on which the stress is projected does not present the correct orientation.
Secondly, this approach requires specific data that are not always easily available.In particular, fluid flow modeling requires the distribution of petrophysical properties of the reservoir and seal rocks; the heterogeneity of these properties will influence the resulting pressure distribution and stress states.In deep saline aquifers, this heterogeneity is difficult to infer on the basis of the few data that are usually available.Even for depleted hydrocarbon fields, information on the seal rocks might be insufficient.A sensitivity analysis should be used to infer the impact of these uncertainties in a second step.
Thirdly, the petrophysical properties of faults are also required to correctly model the fluid flow within the faults.Such information is often only assumed from geological knowledge and a SGR estimation that predicts the fault rock types for simple fault zones developed in sedimentary sequences dominated by sandstones/carbonates and shales (Adeoti et al., 2009;Manzocchi et al., 1999;Yielding, 2002).This parameter is widely used in the oil and gas exploration and production industries to enable quantitative predictions of the hydrodynamic behaviour of faults.
Besides, geomechanical modeling requires the mechanical rock properties that can be derived from log data or laboratory experiments performed on cores.A mechanical damage criterion is also required to assess the potential damage of the rock.
Last but not least, the in situ stress state is very important to correctly initiate the geomechanical model.The regional stress state is available from the World Stress Map Project (Heidbach et al., 2008) but often it does not give a sufficiently accurate estimation of the in situ stress state.Even if the in situ stress state might be available at some wells when data from leak-off tests are available, its spatial distribution is not straightforward: it might vary significantly locally due to geological events such as fractures, fault zones or deformation induced by horizontal stress.

Fluid Flow Model
The fluid flow simulations were carried out using PumaFlow TM .This reservoir simulator is a three-phase flow model based on mass conservation equations for oil species and water, and Darcy laws for flow modeling coupled with thermodynamic equilibrium equations.A classical fully implicit numerical formulation as well as mixing implicit and explicit time discretisation methods is implemented.These equations are discretised in space with a finite volume scheme and linearized with a Newton-type iterative method.

Petrophysical Model
Due to the lack of data, it was not possible to use any geostatistical approach to populate the model.As a consequence, all facies properties were assumed constant per facies.This section describes the methodology for the estimation of petrophysical properties, i.e., porosity, permeability, pore compressibility, capillary pressure and relative permeability associated with the matrix and the fault.

Attribution of Petrophysical Properties to the Matrix
Porosity and permeability values are derived from different sources.Those of the dolomite and carbonate facies are from well logs.The core measurements (Volpi et al., 2013) were carried out to assess the marl porosity.
Regarding shales, no data were available so petrophysical properties were derived from the literature (listed below).As described in the geological setting section, the Plio-Pleistocene formation is refined in four numerical layers to take into account the porosity decrease with depth due to rock compaction.Figure 11 shows the porositydepth curve, which is used in the fluid flow model for shales.The trend is based on several published porositydepth data which have been collected for shale (Hermanrud, 1993).The curve matches the mean value of the porosity data (50%) that was estimated from well logs.From this porosity curve, the porosity values for the four additional facies were determined.They are displayed as red points for facies 11 to 14 and a blue point for facies 10.
Permeability K (m 2 ) for the marl and shale facies was determined by the Kozeny-Carman law (Kozeny, 1927;Carman, 1937Carman, , 1956)): where U is the porosity (fraction) and S o is the specific surface area defined as the quantity of pore surface area per unit volume of solid (m 2 /m 3 ).The value 5 9 10 7 m 2 /m 3 is from the TEMIS TM(5) data base for a facies composed of 90% shale and 10% sand.
Regarding the vertical permeability, the commonly used ratio between vertical and horizontal permeability of 0.1 (Fanchi, 2006) is taken for all facies in the absence of relevant data.
A Corey-type function is used for the relative permeability Kr n of phase n: Kr n ¼ S n r where S n is the saturation of phase n.For carbonates, faults and dolomites, r = 1, whereas r = 2 for marls and shales.In the case of no data measurement, a common value of 20% is used for the irreducible water saturation.
The hysteresis phenomenon on relative permeability and capillary pressure is not considered here, since this study focuses on the CO 2 injection period, which corresponds to a drainage phase.One of the characteristics of the imbibition phase is, however, taken into account by considering a residual gas saturation of 5%.
The threshold pressure Pc th (bar) is computed for gas by the Thomas, Katz and Tek correlation (Thomas et al., 1968): where / is the porosity (fraction).
All the petrophysical data are presented in Table 2. 3.1.1.2Attribution of Petrophysical Properties to the Fault Same functions or correlations are used to evaluate the relative permeability, the capillary pressure and the pore compressibility of the faults and the matrix.As there is no information regarding fault transmissivity, scenarios were run assuming faults to be either open or closed.Closed faults are assumed to be filled with marl or dolomite facies and open faults with platform carbonate facies.
Regarding the porosity, an arbitrary value of 20% is assumed for all the faults whatever their transmissivity.With 20% porosity corresponding to a 80-m-thick fault, the equivalent porosity is 2.5%, which is a reasonably low value for the porosity of a 10-m-thick fault.In the open fault case, the equivalent permeability is 960 9 10 À9 m 2 .For this equivalent porosity/permeability couple, no numerical instability was observed (as might occur for too low a porosity associated with too high a permeability).It might be interesting to infer the impact of the fault porosity uncertainty by conducting a sensitivity analysis of the fluid flow results.
Equivalent porosity and equivalent permeability are used for fault facies in order to correctly model the fault volume and the fluid flow through the fault.

Thermodynamic Model
Pure CO 2 is assumed to be injected in a saline aquifer with 0.035 kg/L salinity, which corresponds to the salinity measured at the Rovesti well in the Scaglia formation.Water density and viscosity are computed from internal correlation models (Schmidt, 1969) based on pressure, temperature and salinity.The Peng-Robinson equation of state is used to calculate the gas density, and the Lohrenz-Clark equation (Lohrenz and Bray, 1964) for the gas viscosity.A comparison with different equations of state (Peng-Robinson, Span-Wagner, Duan-Moller-Weare) and data (Kuhn and Munch, 2013) shows a good fit of the Peng-Robinson model in the temperature and pressure range of the studied storage site.
CO 2 is injected into the Scaglia formation through the Grazia well.At the top of the formation, i.e., 1 634 m subsea depth, the temperature is 40°C and pressure is about 16 MPa.At these temperature and pressure conditions, CO 2 is a supercritical fluid.Its viscosity (0.067 cp) is the viscosity of the gas phase and its density (767 kg/m 3 ) is the density of the liquid phase.However, in the case of open faults, CO 2 can migrate upwards and reach a depth where the temperature or pressure conditions are below the critical point (T = 31.1°Cand P = 73.8bar).In this situation, CO 2 is no longer supercritical.CO 2 dissolution in brine is not considered in this study because the large size of cells would lead to numerical diffusion and, as a consequence, to an overestimation of the amount of dissolved CO 2 in the first years of injection.

Heat Model
In this study, the heat transfer is not taken into account; only the geothermal gradient is considered.This was estimated from temperature values at different depths: Figure 12 displays the temperature data for three different wells and the estimated geothermal gradient of 0.011°C/m.

Initial State
A virgin oil field is located in the south-west part of the studied domain (Fig. 2).However, it is not considered in the fluid flow model because no in situ data were available.As a consequence, the computational domain is assumed to be a saline aquifer.The introduced error is likely negligible since the distance between the oil field and the injection zone exceeds 100 km.
At 2 385 m subsea depth, the pressure measured at the Rovesti well is 24 MPa.With only one piece of pressure data, it is assumed there is no regional circulation and a hydrostatic initial pressure.

Boundary Conditions
Regarding outer boundaries, analytical aquifers (i.e., constant pressure) are used to simulate the fluid flow at the outer boundaries.
The chosen injection well is an old exploration well.It is located 0.4 km from the Rovesti fault (Fig. 2).It is assumed to be vertical.It is assumed to be perforated in the deepest 50 metres of the Scaglia reservoir in order to move the highest overpressure away from the top of the Scaglia formation and thus decrease the risk of fracturing the cap rock.Due to lack of information, a common radius of 0.010 m is used for the injection well.
Preliminary work was performed to design a realistic CO 2 injection rate: overpressure management imposes an injection rate of 1 Mt/year for 10 years.In addition, two CO 2 injection wells each injecting 0.5 Mt/year for 10 years are considered to further reduce the induced overpressure at the wellbore.They will be spaced at about 1 km.

Flow Modeling of CO 2 Injection
Four scenarios were simulated in a 30-year period with a total injection rate of 1 Mt/year over a ten-year period.They correspond to different states of fault transmissivity (closed or open) and different numbers of injection wells (one or two).On average, each computation run took between 1 and 2 hours.Several results were analysed: the maximum overpressure (i.e., difference between the current and initial pore pressure) of each scenario in the matrix; the maximum overpressure in the Rovesti fault; and the CO 2 gas saturation.
Figure 13 gives the CO 2 gas saturation map of four scenarios at 30 years.The diameter of the CO 2 plume is about 2 km.The injected CO 2 is mainly located in the Scaglia formation, but a part also migrates through the lower part of the Oligo-Miocene, which contains carbonate facies.The threshold pressure of marls, which are the other facies composing the Oligo-Miocene, is, however, sufficiently high to prevent CO 2 gas from migrating throughout the marls.When faults are conductive, CO 2 migrates up to the top of the fault (which corresponds to the bottom of the Plio-Pleistocene formation), and slows down its vertical migration that might be completely stopped due to the lower conductivity of the Plio-Pleistocene formation.When faults are open, the extension of the CO 2 plume is smaller in the Scaglia and Oligo-Miocene formations than when faults are closed, since a part of the amount of CO 2 has migrated into the Rovesti fault.The small thickness of the fault results in a very long vertical extension of the CO 2 plume in the fault.When the same rate of CO 2 is injected into two wells instead of one, the CO 2 plume is more extended in the Scaglia formation, and consequently also in the Oligo-Miocene formation, and the overpressure is lower.
Figures 14 to 17 show the overpressure field (more precisely, the overpressure above 1 bar) at different time steps during injection and relaxation periods in two cross-sections: the Rovesti fault (or the neighbouring cross-section in the case of closed faults) and the surface perpendicular to this one.Four scenarios are displayed: closed/open fault scenarios with one or two injection wells.It can be observed that CO 2 injection induces an overpressure field which increases with injection time.Injected CO 2 spreads out around the injection zone, where the overpressure is maximum.Similarly to the CO 2 plume, the overpressure field is mainly located in the Scaglia formation since it is the most permeable formation.As the overpressure depends on both the reservoir capacity to allow fluids to flow (i.e., permeability) and the compressibility and viscosity of the fluids (i.e., thermodynamic properties), the induced overpressure strongly depends on the fault transmissivity: it increases up to 6.9 MPa when faults are closed (Fig. 14), whereas it reaches only 4.4 MPa when they are open (Fig. 15).When faults are conductive (Fig. 15, 17), CO 2 migrates up to the top of the fault, and slows down its vertical migration that might be completely stopped due to the lower conductivity of the Plio-Pleistocene formation.This slowdown/stop induces a secondary maximum overpressure zone.This area which is far from the Scaglia formation, thus shows the highest overpressure just after the end of the injection period.The maximum overpressure inside the Rovesti fault reaches 2 MPa with one well (Fig. 14) and 1.8 MPa with two wells (Fig. 16).Comparison of the one-well and two-well scenarios for the closed fault scenario (Fig. 14, 16) shows that the bottom-hole overpressure is smaller when CO 2 is injected via two wells (5.7 MPa) than when injected via a single one (6.9MPa).This difference (1.2 MPa) is smaller when faults are open (1 MPa, as can be observed in Fig. 15,17).
The overpressure field is mainly extended in the Scaglia formation.Over time, its maximum diameter extends to about 30 km for the open fault scenario (Fig. 15).In this case, the Gondola fault is impacted by the overpressure field in addition to the Rovesti fault, but here the overpressure does not exceed 0.1 MPa.It can be noted that the extension of the overpressure field is not sufficient to impact the boundaries of the fluid flow model, which confirms the assumption of a constant pressure at the boundaries.When the faults are closed (Fig. 14), the fluid flow is constrained by both the Jolly and Rovesti faults.In the first years, only the Rovesti fault is impacted.The overpressure field looks like half of a disk with a diameter of 36 km at the end of the injection.Once the Jolly fault is reached, the fluids flow in the north-west and south-east directions.The overpressure field extends roughly within a rectangle of 58 km and 19 km at 30 years.The Jolly fault is also impacted by the overpressure field but it does not exceed 0.1 MPa.
The maximum overpressure is located in the wellbore of the injection well and it occurs during the injection period.Figure 18 compares the overpressure of the injection wellbore with time for the different scenarios.It can be observed that the maximum overpressure occurs at the end of the injection period for the closed fault scenario, whereas it occurs at the beginning of the injection when faults are open, thus providing a migration pathway.The fluid flow through the Rovesti fault partly absorbs the injection flux.
In conclusion, the fluid flow simulations show the impact of the fault permeability on both the CO 2 migration and the overpressure field evolution.Considering fault as volumetric allows one to simulate potential scenarios of permeable fault behaviour: CO 2 migration and the associated pore pressure variation inside the fault and all along the fault.
These results have to be taken with caution because no comparison with available data was possible.In addition to the uncertainty in the geological model and the petrophysical properties previously mentioned, it should also be noted that discretisation of the fluid flow equations on a rather rough mesh induces some inaccuracy besides neglecting CO 2 dissolution.Therefore, it would be fruitful to perform some refinement in order to reduce the mesh impact on the simulation results.In particular, a key point for the fault modeling would be to refine the fault cells by using the Local Grid Refinement (LGR) methodology instead of the Tartan refinement.This was not possible to achieve within the research project; it should be considered in future work to reduce the uncertainty in the results.

Geomechanical Model
The geomechanical modeling computes stress changes resulting from this pore pressure increase (i.e., overpressure).In a post-processing step, the geomechanical initial stress state is added to the computed stress increment to obtain the total stress.This resulting stress state is then compared with a damage criterion, which is here the commonly used Mohr-Coulomb criterion, in order to evaluate the potential fault damage resulting from CO 2 injection.Since the model is structurally complex, the real initial stress state is also complex and cannot be satisfactorily estimated.As an alternative, various possible initial stress states were evaluated in term of potential fault reactivation.

Elastic Parameters
Geomechanical modeling requires an estimation of the initial stress state (i.e.before CO 2 injection), an estimation of the mechanical parameters of the model and the pore pressure variations resulting from CO 2 injection, which are computed by the fluid flow simulator as described before.
In the geomechanical model, each lithology is assumed to be a homogeneous material associated with a poroelastic behaviour.Faults also have a poroelastic behaviour, and are modeled with surface behaviour elements.Two sets of parameters are thus given for each facies: drained parameters are used when the facies is in a layer where fluids are allowed to flow, undrained parameters are used in a layer where the fluid is not allowed to flow.They are given in Table 4.In the case of closed faults, as pore pressure is not given inside the fault by the fluid flow simulation, the fluid is implicitly taken into account in the fault through undrained parameters.A correction is thus made to take into account the pressure variation due to the volume change: the pressure increment in faults for the closed scenarios is computed assuming where K f is the incompressibility of the saturating fluid (assuming the fluid is water); / is the fault porosity; and e V is the volumetric strain.

Boundary Conditions
As the model is embedded in the earth, each side of the geomechanical grid has to be constrained in displacement or in force.Commonly used boundary conditions are chosen: the top surface is a free surface and the bottom surface is assumed to have a null vertical displacement.The side surfaces are chosen far from the hydromechanical solicitation in order to avoid any mechanical impact in the area under interest.The side surfaces of the model are assumed to have constant regional stresses.

Initial Conditions -The Regional Stress State
The potential CO 2 storage site is located offshore in the Southern Adriatic Sea, closed to Bari (Fig. 1).The presentday geological setting of this part of the Italian peninsula is very complex due to the interaction of different geodynamic processes that have been acting closely in space and time.Both compression and extensional processes coexist (Amato and Montone, 1997;del Gaudio et al., 2007).As a consequence, two stress regimes might coexist: a normal fault stress regime (r v > r Hmax !r hmin ) in the South (del Gaudio et al., 2007); a strike-slip fault stress regime (r Hmax > r v > r hmin ) in the North; where r Hmax , r hmin and r v correspond, respectively, to the maximum total horizontal stress, the minimum total horizontal stress and the total vertical stress; -The World Stress Map Project provides very few data in the area under investigation (Fig. 19).Indeed, among the offshore wells examined in the area for breakout detection, some have been discarded and, for some other ones, no breakout was detected (Amato and Montone, 1997).Plausible explanations proposed by Amato and Montone (1997) are either a horizontally isotropic stress, or a high rock strength, which prevents the genesis of breakout.Beyond the studied area in the South in the foreland of the southern Apennines, r Hmax orientation is mainly NW-SE between N120°and N135° (Montone et al., 2004).Numerical modeling recently performed by Petricca et al. (2013) shows that the r Hmax orientation in the study area turns from N135 to N180 from the South to the North (Fig. 20).
The following hypotheses are derived from Amato and Montone (1997), del Gaudio et al. (2007) and Montone et al. (2004): for the North part, a strike-slip fault stress regime with r Hmax = 1.1 r v , oriented N135; r hmin = 0.9 r v ; for the South part, a normal stress regime with r Hmax = r hmin = 0.9 r v .
For the sake of simplicity, the regional initial stress state is assumed to be homogeneous in the whole model.To illustrate the influence of this initial state, two stress regimes are considered: -NF: Normal Fault regime with r Hmax = r hmin = 0.9 r v ; -SS135: Strike-Slip regime with r Hmax = 1.1 r v , oriented N135, and r hmin = 0.9 r v .

Geomechanical Modeling
The geomechanical simulator computes the stress variations resulting from the pore pressure variations induced by CO 2 injection.These pressure variations were computed for the 30 years period: to perform the geomechanical computation, 12 time steps were chosen for all the scenarios to analyse the evolution of the stress state over the injection and the relaxation periods.Therefore, total stress is obtained by adding the initial stress to the calculated stress variations.The stress which is referred to here is the effective stress, i.e., the stress to which the rock is submitted.With the compression negative convention, this effective stress is: where r total is the total stress, r effective is the effective stress, P is the pore pressure, b is the Biot coefficient (equal to 1 in this study) and I is the identity matrix.The results presented in this paper correspond to the Rovesti fault, since this fault is the closest to the injection well and thus experiences the highest pressure elevation Figure 21 shows an example of variation for shear stress and normal effective stress at the end of the injection period, for the open fault case.The stress changes are mainly normal and the shear stress variation remains very low.

Interpretation of Geomechanical Results in Terms of Potential Fault Reactivation
Assessment of geomechanical potential fault damage (possibly leading to reactivation) related to faults consists of comparing the stress state on the fault to a damage criterion/ envelope that defines the limit between elastic and inelastic behaviour.This envelope can be determined through laboratory tests: various stress paths are imposed on a given rock sample, and when damage occurs, the envelope is gradually determined.
The criterion chosen in this study is the Mohr-Coulomb envelope, commonly used for faults.This criterion is represented by a line with a slope, tan(/), governed by the friction /, and an intercept (governed by the cohesion), in a stress space defined by normal and shear stress.The common value for cohesion in faults is 0, which was chosen here.For friction, according to Byerlee's faulting theory (Byerlee, 1978), given a slope with a 0.6 value, a value of 30°for friction is chosen.
As illustrated in Figure 22 left, the stress vector variation due to pressure increment is composed of Dr n , which is the resulting increment of the effective normal component (normal to fault plane) and Ds elas which is the resulting increment shear component (contained in the fault plane).To compute the final stress (black point of Fig. 22b), the stress vector variation is added to the effective initial stress (green star of Fig. 22b).The final stress state projected on   r Hmax direction and stress regimes in the study area (Heidbach et al., 2008).the fault is compared with the Mohr-Coulomb criterion.If it is above, some damage has occurred on the fault.
To evaluate the impact of CO 2 injection on the potential fault damage, the distance of the stress state to the Mohr-Coulomb criterion is computed, as illustrated in Figure 23, with the following equation: where / is the friction angle and s the shear stress.
- The impact of CO 2 injection can be clearly observed on the fault plane for all scenarios.It should first be noted that the range of DCrit is roughly [À0.7 MPa; +1.2 MPa], whereas the range of Crit INI is around 9.8 MPa in the reservoir, near the injection point, so that CO 2 injection has a low impact on the distance to the Mohr-Coulomb criterion.This, of course, is strongly related to the injection scenario (position of injection wells, rate of injected CO 2 ).DCrit varies with time, increasing or decreasing according to the scenarios.It goes to zero when the injection is stopped because of the decrease in the pressure with time and the elastic behaviour assumption.The size of the DCrit footprint depends on the stress regime, but in any case its magnitude remains smaller than Crit INI which makes Crit END remain negative and the fault stable.
Figure 25 presents the results obtained with the SS135 initial stress regime for the open fault scenario at 2, 10, 12 and 30 years. the end of the relaxation period becomes similar to the state before CO 2 injection.Figure 27 presents the results obtained with the SS135 initial stress regime for the closed fault scenario at 2, 10, 12 and 30 years.Figure 28 presents the same scenario with the NF initial stress regime.For both initial stress regimes, CO 2 injection increases the distance to the Mohr-Coulomb criterion line at the reservoir level and slightly decreases it in the reservoir vicinity.It is observed that CO 2 injection makes the stress state in the fault closer to the Mohr-Coulomb criterion at the reservoir level, but here mechanical effects make the stress state above the reservoir level farther.The effect of CO 2 injection is here much higher for the SS135 initial stress regime than for the NF initial stress regime.These results illustrate the importance of the initial stress regime.
It is actually tricky to compare the closed and open faults scenarios.Even if the same amount of CO 2 has been injected, the variation of pore pressure due to CO 2 injection is different in the two scenarios due to the hydraulic fault conductivity difference.In the closed fault scenario, the pore pressure increases in the fault at the reservoir level but less than when faults are open.Indeed, fluids are not allowed to flow in the fault so that the increase in pore pressure in the fault is only due to a geomechanical effect because of its undrained behaviour.

CONCLUDING REMARKS
As part of an investigation of potential CO 2 storage in carbonate formations, the Grazia site was investigated to test a methodology of coupling a fluid flow model with a geomechanical simulator.The focus was placed on the assessment of possible fault reactivation induced by CO 2 injection.The coupling methodology between fluid flow modeling and geomechanical modeling is an iterative process that was performed in six steps: construction of the geological model, in conformity with both PumaFlow TM (fluid flow simulator) and Abaqus FEA (geomechanical simulator) requirements; conversion of the "Eclipse format" model into a finite element model for the geomechanical model using an IFPEN house-built programme; construction of the fluid flow model; -fluid flow modeling to compute overpressure inside each cell; geomechanical simulations to compute effective stress variation due to overpressure; post-processing to estimate the fault stress state distance to a failure line, while comparing the final stress state with the Mohr-Coulomb criterion for instance.Faults are complex elements that require various representations according to the investigated property and the modeling software to be used.To make the geometries required by the geomechanical and the fluid flow modeling compatible, faults are modeled with cohesive elements in the geomechanical model, and fault neighbouring cells were used to model fluid flow in the fault.Besides, to handle the uncertainties due to lack of data, several scenarios were simulated associated with different petrophysical properties, the number of injection wells and fault transmissivity.
The fluid flow simulations show that modeling faults as volumetric elements allows faults when conductive to behave as flow pathways for fluids.It is then observed that the injected CO 2 migrates inside and through the Rovesti fault, which is located near the injection well.The fluid flow also induces overpressure in the faults.The overpressure in the Rovesti fault reaches 2.2 MPa, while it reaches 4.4 MPa at the bottom hole of the injector.Extending to about 30 km, the pore pressure field reaches the Gondola fault located at 15 km from the injection zone, but the overpressure does not exceed 0.1 MPa at such a distance from the injection well.When non-conductive, faults behave as barriers and the overpressure field extends about 60 km along the Jolly and Gondola faults.There is no overpressure inside the faults, but the difference in pressure between both sides of the faults can potentially cause damage to the faults.Postprocessing of the geomechanical stresses was performed to compute the resulting distance of the stress state from the Mohr-Coulomb criterion line, according to various stress regimes, and assess the resulting risk of fault reactivation.On the basis of available data, all simulated scenarios show that the Rovesti fault, which is near the injection well, remains below the Mohr-Coulomb criterion line.However, the results have to be considered with care because of the few available data to inform the fluid flow and geomechanical models.Besides, in addition to the mesh refinement level, some assumptions, such as ignoring CO 2 dissolution, were made and impact the simulation results.Therefore, increasing the result accuracy requires both improving the geomechanical/fluid flow coupling and integrating more data.The focus is thus placed on the methodology rather than on the results themselves.
The next step could be to use a SGR estimation, adapted to geological context, to derive the fault permeability or a law describing the mechanical behaviour.Techniques of homogenisation could then be used to achieve this goal.This was out of the scope of this study but it should achieve a significant improvement in the approach and more realistic results.

Figure 1
Figure1Location of the potential storage site in the Italian peninsula (our model is shown by the grey rectangle) and the Brindisi power plant in the Southern Adriatic Sea site (Google Maps).

Fault
Figure 4 Various vertical representations of a fault in a) the geological model, b) the fluid flow model and c) the geomechanical model.
Figure 6Fault discretisation (horizontal view), illustrating the need for an additional fault cell.
Figure 8Behaviour law for cohesive element.
Figure 10 Final refined grid.Showing map (left) in the Scaglia reservoir (basinal and platform carbonates facies) and vertical section AA' (right).
th ¼ 7:37 Â K À0:43 Â 0:06895: CO 2 is a supercritical fluid in the Scaglia reservoir.However, if the faults are open, the CO 2 phase might change from supercritical fluid to gas as it migrates upwards.The InterFacial Tension (IFT) therefore increases from about 30 mN/m (supercritical CO 2 ) to about 70 mN/m (gas).As it is not possible to consider this IFT variation in the model, the correlation used in this study is adapted for upper layers but overestimates the threshold pressure of the reservoir and around.The expected effect is an overestimation of the pore pressure.The pore volume compressibility C (1/10 5 Pa) is estimated by the Hall correlation(Hall, 1953): Figure 11Shale porosity depth (m) for shales (facies 10 to 14), adapted from Hermanrud (1993).5 IFPEN basin simulator.
Figure 12Mean geothermal gradient from temperature data measured at the Grifone, Rovesti and Grazia wells (data provided by OGS).
Figure 13 CO 2 gas saturation at 30 years for the (top) 1 injection well and (bottom) 2 injection wells in the case of (left) closed and (right) open scenarios.
Figure 14 3D overpressure for the closed fault scenario with one injection well (top) during the injection period (2 and 10 years) and (bottom) during the relaxation period (12 and 30 years).
Figure 15 3D overpressure for the open fault scenario with one injection well (top) during the injection period (2 and 10 years) and (bottom) during the relaxation period (12 and 30 years).
Figure 16 3D overpressure for two injection wells in the case of the closed fault scenario (top) during the injection period (2 and 10 years) and (bottom) during the relaxation period (12 and 30 years).
Figure 17 3D overpressure for two injection wells in the case of the open fault scenario (top) during the injection period (2 and 10 years) and (bottom) during the relaxation period (12 and 30 years).
Figure 18 Wellbore (bottom-hole) overpressure with time for 1 or 2 wells: open/closed fault scenarios.

r
Figure 21 Stress variation on the Rovesti fault in the case of an open fault after 10 years of injection.These two panels are a zoom of the Rovesti fault around the injection well (lateral extent of 26 km; vertical extent of 5.6 km from the seabed), located by a yellow rectangle.
Figure 23 Schematic representation of the criterion, DCrit = Crit END À Crit INI , used to assess the potential damage of the fault.A strike-slip stress regime is chosen.

Figure 27 Figure 28 Figure 26 Figure 25
Figure 27DCrit parameter on the Rovesti fault (seen from SW) for the closed fault scenario and the SS135 initial stress regime at 2, 10, 12 and 30 years.These two panels are a zoom of the Rovesti fault around the injection well (lateral extent of 26 km; vertical extent of 5.6 km from the seabed), located by a yellow rectangle.