Evaluation and Characterization of a Potential CO 2 Storage Site in the South Adriatic Offshore

The Southern Adriatic Sea is one of the five prospective areas for CO2 storage being evaluated under the FP7 European Sitechar project. The potential reservoir identified in the investigated area is represented by a carbonate formation (Scaglia Formation – Late Cretaceous). This paper shows the site characterization applied to one of the structures identified in the carbonate storage system of the South Adriatic offshore. The interpretation and analysis of seismic and borehole data allowed the construction of a 3D geological static model on both regional and local scales. Dynamic modeling was applied, adopting a sensitivity approach ( i.e. fault transmissivity, petrophysical properties of the caprock and reservoir, and different stress regimes). Coupled fluid flow and geomechanical simulation was applied to investigate the potential risk of leakage induced by mechanical solicitation on the faults occurring in the investigated area. Résumé — Évaluation et caractérisation d’un site de stockage potentiel de CO2 au sud de la Mer Adriatique— Le sud de la Mer Adriatique est l’une des cinq zones potentielles pour le stockage du CO2 évaluées dans le cadre du projet FP7 européen SiteChar. Le réservoir potentiel identifié dans la zone étudiée est constitué d’une formation de carbonate (Formation Scaglia – Crétacé supérieur). Cet article présente la caractérisation de site de stockage de CO2 appliquée à l’une des structures identifiées dans les formations carbonatées du sud de la Mer Adriatique. L’interprétation et l’analyse des données sismiques et des données de forage ont permis la construction d’un modèle statique géologique 3D à des échelles régionale et locale. Une modélisation dynamique a été effectuée, intégrant une analyse de sensibilité (sur la transmittivité des failles, les propriétés pétro-physiques des roches couverture et réservoir, et différents régimes de contraintes). Une simulation couplée géomécanique et d’écoulement de fluide a été menée pour étudier le risque potentiel de fuite provoqué par les contraintes mécaniques sur les failles de la zone étudiée. Oil & Gas Science and Technology – Rev. IFP Energies nouvelles, Vol. 70 (2015), No. 4, pp. 695-712 V. Volpi et al., published by IFP Energies nouvelles, 2015 DOI: 10.2516/ogst/2015011 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.


INTRODUCTION
Studies conducted on a regional scale on the Italian territory have identified several areas, both onshore and offshore, potentially suitable for CO 2 geological storage in deep saline aquifers within siliciclastic (Donda et al., 2011;Buttinelli et al., 2011) and carbonate successions (Civile et al., 2013).In this framework, the Southern Adriatic Sea revealed good potential conditions of porosity and permeability in some carbonate formations (Volpi et al., 2015).In this area, on the Puglia coast close to Brindisi, the main Italian CO 2 source emission is located, Enel's Federico II power plant, with about 16 Mt/yr (Fig. 1a).The Federico II is a coal-fired power plant which operates four 660 MW coal-fired units with a total capacity of 2 640 MW, where on March 2011, Enel launched Italy's first carbon capture pilot.The vicinity of the power plant to the area suitable for CO 2 storage makes the Southern Adriatic Sea a strategic area to perform the whole chain of the site characterization (static and dynamic modeling), also facing some particular issues related to storage in carbonate formations.
The sites potentially suitable for CO 2 geological storage are structural highs, probably generated by Jurassic salt tectonics, involving the Upper Triassic evaporites, and locally by episodic Cretaceous to Quaternary strike-slip tectonics developed along E-W-and NNE-trending tectonic lineaments (De Dominicis and Mazzoldi, 1987;De Alteriis and Aiello, 1993) (Fig. 1b).The Southern Adriatic offshore is considered as relatively stable; moderate tectonic activity is concentrated along the front of the circum-Adriatic orogenic chains (Southern Apennines and Ellenides-Albanides) and the E-W-trending Gondola seismogenic strike-slip fault zone located to the south of the Gargano promontory (Ridente et al., 2008;Di Bucci et al., 2009) (Fig. 1b).
The reservoir for the geological storage is represented by the Upper Cretaceous-Paleogene carbonate slope succession of the Scaglia Formation.This succession consists of mudstones to bioclastic packstones saturated by brine, with a variable thickness ranging from 10 to 100 m.The top of the Scaglia Formation, lying between 1 000 and 2 000 m depth, corresponds to the top of the platform succession representing the foreland (Apulian foreland) of the Southern Apennine chain.The potential caprock is represented by Upper Oligocene to Quaternary predominantly siliciclastic succession, several hundred meters thick (Volpi et al., 2015).After a first evaluation of the storage potential, considering the occurrence of suitable geological conditions and the safety of the storage system in terms of seismic risk, in this paper, we present the results of dynamic modeling consisting of fluid flow and coupled fluid flow/geomechanical simulations performed on the most suitable structure (Grazia Anticline) identified in the Southern Adriatic offshore (Fig. 1a, 2).The regional tectonic lineaments recognized in the study area were assessed in the modeling in order to infer the potential risk of CO 2 leakage produced by changes in the stress state.This required modeling the coupled fluid flow and geomechanical reservoir and seal behavior.

Seismic and Well Log Data
The South Adriatic offshore is covered by a large amount of data acquired by oil companies for hydrocarbon exploration and made available by the Ministry of Economic Development in the framework of the project "Visibility of Petroleum Exploration Data in Italy (ViDEPI)".About 1 500 km of 2D multichannel seismic lines and 15 sets of well data were analyzed to characterize the study area (Fig. 3).The public well data consist of "composite log charts" that contain the following information: lithology derived from cuttings; geological formation name (following the nomenclature provided by wells, which may vary among different areas and in time); formation age; depth, lithostratigraphy; occurrence and type of fluid; depositional environment; biostratigraphy; geophysical log curves (commonly resistivity, spontaneous potential, sonic and gamma ray) (Fig. 4).Where available, pressure and temperature values at depth, as well as porosity and permeability data from cores were also reported.

Petrophysical and Geomechanical Properties
Coupled fluid flow-geomechanical simulation, which allows one to compute the pressure and stress fields and simulate the CO 2 plume migration, requires petrophysical and geomechanical characterization of the storage system to associate property data with the matrix rock and the faults.The only available porosity and permeability values of the reservoir formations, ranging between 2.5 and 24% and between 2.6 and 2.7 mD, respectively, come from Grifone well, which is located about 60 km from the analyzed Grazia  structure (which lies below the homonymous well) (Fig. 1a).Thus, petrophysical properties were derived from different sources: well log analysis, measurements performed in the laboratory on core samples taken from outcropping analogs or the scientific literature.
The different methods for calculating the petrophysical and mechanical properties show a degree of uncertainty: the estimated permeability of the caprock obtained from well logs is too high to ensure the sealing characteristics; the measurements derived from outcropping analogs remain challenging because the burial history and diagenesis of the outcrop and the subsurface sample will probably be different.However, considering these uncertainties, four matrix scenarios were used in the dynamic simulation, integrating the petrophysical datasets obtained by the different methods.

Matrix Scenarios
The different petrophysical datasets considered in the dynamic modeling are: -LOGS dataset (Tab.1): corresponds to data derived from well log data.Caprock porosity was computed by applying a velocity-porosity relation proposed by Han (1986) for clastic sediments.A sonic-derived porosity (SPHI) was computed for the reservoir rocks (Scaglia formation) by applying the classical Wyllie time-average relation (Wyllie et al., 1956).Permeability values for the carbonate reservoir were computed starting from the rock-fabric petrophysical classification proposed by Lucia (1995), considering the empirical relationship based on the correlation between permeability and fractional interparticle porosity.For the clastic sediment of the caprock, the permeability was estimated using the relationship proposed by Xiangui et al. (2009).Density values for the caprock and reservoir were derived using empirical formulas reported in the literature (Gardner et al., 1974;Castagna et al., 1993;Brocher, 2005).
Regarding the mechanical properties, the Poisson ratio was computed by the relation proposed by Brocher (2005) and the elastic moduli (Young Modulus, Bulk Modulus and Shear Modulus) required for the geomechanical modeling were then obtained from the p-wave velocity formula, knowing the density and the Poisson ratio r: -BIBLIO dataset (Tab.2): this matrix scenario is based on the LOGS dataset except for the petrophysical data of caprock facies, which are replaced by values derived from published shale porosity-depth trends collected by Hermanrud (1993).A mean porosity trend was estimated from well logs and compared with the empirical porosity curves.The porosity values for the facies located at a certain depth were deduced using the porosity curve that fits the mean porosity curve derived from well logs.Permeability K for the caprock facies was determined by the model proposed by Temis TM (IFPEN basin simulator).It is based on Kozeny-Carman's law (Kozeny, 1927;Carman, 1937Carman, , 1956)): where U is the porosity and S 0 is the specific surface area defined as the quantity of pore surface area per unit volume of solid; -CORE dataset (Tab.3): this scenario is based on the LOGS dataset, but the petrophysical data of the caprock and the Example of a "composite log chart".
reservoir formations were updated using core measurements.The core samples were taken from an onshore analog succession outcropping in the area of Vieste-Peschici (Gargano promontory) (Fig. 5).The samples, collected from three different successions, consisted of limestones (reservoir analog) of the Turonian-Senonian Scaglia Formation, and marls (caprock analog) of the Lower Aptian Fucoidi Marls Formation.The Fucoidi Marls have a lithology and depositional environment very similar to the Scaglia Cinerea Formation that represents the lower portion of the caprock succession in the study area (Fig. 5).The following properties were measured by laboratory analysis: porosity, irreducible water saturation, threshold pressure, permeability, Young modulus, Poisson ratio, Coulomb parameters (cohesion and internal friction angle) and tensile strength; -OAb dataset (Tab.4): the dataset of this scenario is a synthesis of the three previous cases, where the shale facies properties of the caprock come from the BIBLIO scenario; the marl facies of the caprock come from the CORE scenario.The properties of the platform, basinal carbonates and dolomite facies of the reservoir come from the LOGS scenario.
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.In addition, since the Plio-Pleistocene Formation has been refined in 4 layers, it was decided to define one shale facies per layer in order to take into account the porosity decrease with depth due to rock compaction.

Faults
The same functions or correlations are used to evaluate the relative permeability, capillary pressure and pore compressibility of the faults as for the matrix.Transmissivity: as there is no information regarding fault transmissivity, three fault scenarios are studied, assuming faults to be either open, closed or in an intermediary case.
Permeability: closed faults are assumed to be filled with marl or dolomite facies, faults associated with an intermediate status with basinal carbonate facies and open faults with platform carbonate facies.
Porosity: an arbitrary value of 20% is assumed for all the faults whatever their transmissivity.This value was chosen in order to maximize the equivalent porosity, and consequently to limit numerical issues.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.

Static Model Building
The objective of the static modeling is to create a 3D computerized representation of the storage complex based on geophysical and geological data.The geometry of the storage complex was reconstructed by integrating the seismic data interpretation with borehole information.The 3D structural model is composed of the identified key seismic surfaces, including the top and bottom of the reservoir, caprock and faults.The volumetric grid is obtained from the structural model by splitting it up into a finite number of boxes/cells which can be populated by lithological, petrophysical and geomechanical properties.This final model is the "static model", which is a representation of the storage complex that does not consider any time-dependent variability of its internal properties.
The workflow for the static model building consists of six main stages (Fig. 6): data input; data analysis: identify key seismic horizons and tectonic features in order to build a 3D structural model; creation of the structural model; creation of a volumetric grid; facies modeling; matrix attribution: which consists of creating the petrophysical model (which requires the matrix attribution of petrophysical properties).

Structural Model
The structural model is a 3D representation on a regional scale of the four identified key seismic horizons: the Seabed, Top Miocene succession, Top and Bottom Reservoir/ Scaglia Formation and base of the Pre-Scaglia, which represents an artificial surface that defines the base of the model (Fig. 2) and the major tectonic lineaments (Jolly, Gondola, Rovesti, Grifone, Rosaria and Cigno tectonic lineaments) affecting the study area (Fig. 7a, b).These lineaments are sealed by Plio-Pleistocene successions, except for Grifone and Gondola, which seem to deform the sea floor.Structural interpretation of available seismic and borehole data was performed in a Two-Way Time (TWT) domain and successively converted to depth.The depth conversion required the creation of a velocity model calculated using velocity functions (i.e.sonic log) available at five wells (Grazia 1, Rovesti 1, Cigno Mare, Gondola 1 and Jolly 1) (Fig. 3 for locations).Grid maps of the average velocity for each horizon were created using the interpolation algorithm Ordinary Kriging.On the basis of this velocity model, the 3D structural-stratigraphic model was converted to depth and any possible mismatch in the depth of the horizons was cross-checked with the well composite logs.

Volumetric Gridding
The following four stratigraphic successions, bounded by the interpreted key seismic horizons, were defined in the volumetric grid from top to bottom: the Pliocene-Pleistocene, the Oligocene-Miocene (caprock), the upper Cretaceous-Eocene (reservoir -Scaglia Formation) and the Pre-Scaglia deposits (Fig. 7c, d).The maximum number of cells was connected to the minimum size of the features to be resolved.In this study, the vertical refinement was performed by dividing the Plio-Pleistocene and Oligo-Miocene successions into four and six "numerical" layers, respectively.Moreover, the potential development of artifacts in the computational domain along the bottom boundary of the model was avoided by assigning an arbitrary thickness of 500 m to the Pre-Scaglia succession.This arbitrary layer was subdivided into ten numerical layers, but only the shallowest one was considered in the fluid flow dynamic modeling.The horizontal grid cell dimensions were uniformly equal to 1 9 1 km, providing a final model grid with 318 600 cells (Fig. 7d).
The final 3D volumetric grid on a regional scale was the result of different adjustments to satisfy the computational constraints of dynamic modeling, while at the same time providing a good fit to the geological surfaces interpreted using seismic and borehole data.

Lithological Modeling
Five lithologies were identified in the study area and associated with the four sedimentary successions included in the The location of the three samples collected on the Gargano Promontory on the right, the composite basinal succession: the samples collected for the laboratory test are marked.
3D geological model (Tab.5): shales, marls, platform carbonate deposits, carbonate pelagic deposits and dolomites.Dolomites are predominantly in the Pre-Scaglia succession.Carbonate pelagic deposits include pelagic mudstoneswackestones deposited in a deep basinal environment and the slope deposits containing bioclastic packstones.This lithology was associated with the Upper Cretaceous-Paleogene Scaglia Formation, which is composed of 10-100-m-thick mudstones to bioclastic packstones sedimented in a carbonate slope depositional context.Shale and marls were associated with the Oligocene-Miocene caprock succession.
To create a distribution of the lithologies in the final 3D grid, a discrete log, representing the lithological distribution in depth, was created for all available wells.These logs were then upscaled to the size of the grid cells: the discrete value associated with the most represented lithology within the intersected grid cell was selected.In order to perform the stochastic lithological modeling, the initial proportion of each facies within each zone (derived from upscaled well data) was considered and, starting from this initial proportion, the vertical proportion curves, that express the probability that a certain facies exists in a certain layer, were derived.The vertical proportion curve shows how the proportion of each facies varies vertically in a stratigraphic zone (Fig. 8).Finally, the experimental variograms were obtained for each facies within each zone in order to define the variogram models applied using the truncated Gaussian simulation algorithm.This algorithm was chosen because it produces reliable results in systems where a natural transition through a sequence of facies is recognized (typical examples include carbonate depositional environments).
The results obtained are in agreement with the conceptual geological model of the area.

Facies-Based Petrophysical Property Attribution
The lack of available petrophysical data did not allow us to populate the model with petrophysical properties using a geostatistical approach.These properties were considered constant for all the identified facies.The uncertainties were handled considering several scenarios in the dynamic modeling, combining the petrophysical properties derived from well log analysis, core measurements and the scientific literature, as described in Section 1.2 (Tab.1-4).However, each scenario presents some critical issues: -LOGS dataset (Tab.1): the critical aspect is the calculated caprock permeability, which is too high to seal the injected CO 2 ; -BIBLIO dataset (Tab.2): in this scenario only parameters from the literature were considered, no direct measurements were included; -CORE dataset (Tab.3): this scenario presents better seal characteristics of the caprock marls (with lower permeability than that of the LOG scenario) but also very low permeability of the reservoir formation; -OAb dataset (Tab.4): in this scenario the reservoir is characterized by the highest permeability and the caprock with the best sealing capacity, thus representing an ideal situation.
Despite the uncertainty induced by the use of analogue core data or limited log data, these scenarios are still needed to initially characterize the storage site in the absence of essential field data, as is the case with the characterization of the Brindisi site.Uncertainty of the simulation results could be reduced by the use of more field static data; even a geostatistic model used to interpolate these local field data would also induce uncertainty in the simulation results.During the storage life, numerous geological, fluid flow and geomechanical models will be built by history matching, integrating more and more field static data and above all, dynamic data.These data would be very useful because they are the result of 4D coupled phenomena.Thus, field static data are always necessary to initially constrain the models but then, they have to be updated during the history matching process.

Dynamic Modeling
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 discretization methods are implemented.These equations are discretized in space with a finite volume scheme and linearized with a Newton-type iterative method.(IFPEN in-house Puma-Flow TM , trademark of Beicip Franlab.http://www.openflowsuite.com/).Workflow chart applied to build the geological static model.
The simulations were carried out using a two-phase model.This model is characterized by: the petrophysical properties defined in Section 2.1.4,a thermodynamic model based on an equation of state to simulate the different fluid states of the CO 2 at different temperatures and pressures, a heat model, an initial state, and boundary conditions such as the injector.
All these points are described in Baroni et al. (2015).This model was applied on a refined mesh (Baroni et al., 2015) built from the geological gridding (Sect.2.1.2).In addition to a vertical refinement, a horizontal Tartan refinement was chosen for the horizontal refinement of the CO 2 injection area to reduce numerical diffusion and thus to simulate the pressure gradients and the CO 2 plume migration better.The final mesh comprises 2.1 million cells.
The dynamic simulation demonstrated that the storage assessment is not capable of containing 10 Mt/years for 40 years of injection, which is the case where we would be able to store the total amount of CO 2 emitted at Brindisi power plant.The best CCS scenario for a pilot project was then evaluated considering a mean injectivity of 1 Mt CO 2 /year, injecting in 10 years a total of 10 Mtons of   The petrophysical properties of each scenario are given in Tables 1-4.
There is a large uncertainty on the reservoir injectivity, which could range from 0.01 Mt CO 2 /year to 1 Mt CO 2 /year.No demo or a pilot site are planned yet to give insight.In this context, 1 Mt/year for 10 years was assumed as reasonable, even if it is probably much lower than the estimated storage capacity of the reservoir.
Simulations performed on the CORE matrix scenario led to an overpressure exceeding 350 bar with a CO 2 injection rate of 1 Mt/year whatever the fault transmissivity.This is explained by the very low carbonate permeability of 0.19 mD.Therefore, it was decided not to simulate the fluid flow and geomechanical coupling for this matrix scenario.
Fluid flow simulations were run for the three other scenarios and the following results were analyzed: -CO 2 gas saturation, maximum overpressure of each scenario in the matrix, maximum overpressure at the Rovesti fault.
As the results of the intermediary fault scenarios are similar to those of the open fault scenarios, they are not presented here.

Fluid Flow Modeling
Figure 9 shows the CO 2 gas saturation map of the six scenarios at 30 years (i.e., three matrix scenarios and two fault scenarios).The shape of the CO 2 plume depends on the fault transmissivity, the caprock sealing capacity and the number of injection wells.Injected CO 2 remains mainly in the Scaglia Formation but a part also migrates through the lower part of the Oligo-Miocene caprock succession which contains Miocene carbonate deposits.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 caprock.
When faults are conductive (Fig. 9 right), CO 2 migrates up to the top of the fault.Here, the vertical migration of CO 2 slows down without stopping when the Pliocene-Pleistocene succession is conductive (1-3 mD for the LOGS case).It is completely stopped in the case of the very low conductivity of the Pliocene-Pleistocene succession in the OAb scenario.This slowdown/stop induces a small overpressure zone.This zone, which is far from the Scaglia Formation, thus shows the highest overpressure that occurs just after the end of the injection period.For all the simulations, it was observed that CO 2 injection induces an overpressure field that increases with time during the CO 2 injection period.Its diameter is about 14 km for the LOGS -open fault scenario but it extends to about 30 km for the OAb-open fault scenario (Fig. 10) and reaches the Gondola fault.The overpressure actually depends on the reservoir capacity to allow fluids to flow (i.e., permeability), and the compressibility and viscosity of the fluids (i.e., thermodynamic properties).The maximum overpressure is located at the wellbore of the injection well and occurs during the injection period.As observed in Figure 11, the OAb model with the closed faults induces the highest wellbore overpressure (67 bar) among the three matrix scenarios (48 bar for LOGS, 57 bar for BIBLIO).These differences come from the caprock sealing capacity.The sealing capacity of the caprock is actually quite poor for the LOGS scenario, with 1-3 mD permeability, whereas the best sealing capacity is observed for the OAb scenario (permeability of 0.0004 mD).
The BIBLIO case presents a sealing capacity between the two previous cases, i.e., 0.004 mD. Figure 12 shows that the overpressure zone simulated with the LOGS scenario crosses the caprock and reaches the sea bottom.It presents a vertical shape.Thus, in this scenario, the caprock appears as a migration pathway for fluids and this requires taking measures to reduce the overpressure.The overpressure obtained for the OAb scenario penetrates only the lower part of the caprock due to zones composed of carbonate facies.In this case, the overpressure zone has a mainly horizontal shape.All these results of course strongly depend on the geological model and the associated petrophysical properties.As they were determined with sparse available data,  the sealing capacity of the caprock remains quite uncertain and there is a need for a better characterization so as to gain confidence in the simulation results.
Comparison of the wellbore overpressure resulting from the OAb-closed fault scenario (Fig. 11) with that resulting from the OAb-open fault scenario shows that the fault conductivity allows fluids to migrate inside the fault, cross the fault and finally migrate to the other side (Fig. 12).This leads to a decrease in the overpressure from 67 bar to 41 bar.The other matrix scenarios show a smaller decrease than for OAb: 18 bar for BIBLIO and 12 bar for LOGS.Comparison of scenarios with one and two wells for the OAb-closed fault scenario (Fig. 11,12) shows that the wellbore overpressure is smaller when CO 2 is injected through two wells (55 bar) instead of a single one (67 bar).This effect (12 bar) is of course smaller when faults are open (9 bar).Lastly, the overpressure inside the open faults (Fig. 12) reaches 11 bar with the LOGS scenario and 20 bar with the OAb scenario.With the OAb scenario and two injection wells, it reaches 19 bar.

Geomechanical Modeling
The geomechanical simulations were performed with the Abaqus/Standard software.It is a general-purpose Finite Element software with an implicit integration scheme, developed by Dassault Systemes (http://www.3ds.com/products-services/simulia/products/abaqus).
The aim of the geomechanical modeling is to simulate the geomechanical response after CO 2 injection and investigate any potential damage of faults.Geomechanical modeling computes the stress changes that result from the pressure variations (or overpressure) calculated by the fluid flow simulation.The resulting stress, the sum of the stress changes and the initial stress, is then compared with a damage criterion, which is here the commonly used Mohr-Coulomb criterion.No feedback from the geomechanical modeling, e.g.changes in petrophysical properties, is given to the fluid flow simulator.We thus here consider a so-called "one-way coupling".
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 pressure variations resulting from CO 2 injection, which are computed by the fluid flow simulator as described above.The initial stress state is most often poorly known.The World Stress Map gives an idea of the regional stress and some assumptions are made to derive local variations.Based on the work of Petricca et al. (2013), two stress regimes were assumed and tested: the "SS135" regime, which corresponds to a strike-slip fault stress regime with r Hmax = 1.1 r v , oriented N135 and r vhmin = 0.9 r v ; the "NF" regime, which corresponds to a normal stress regime with r Hmax = r hmin = 0.9 r v .where r Hmax , r hmin and r v correspond, respectively, to the maximum horizontal stress, the minimum horizontal stress and the vertical stress.Only the "NF" regime is discussed in this paper.The reader can refer to Baroni et al. (2015) for results obtained with the SS135 stress regime.
The geomechanical model used in this study assumes that each layer is filled with a material associated with a poroelastic behavior.Faults also have a poroelastic behavior, and they are modeled with surfacic elements (Baroni et al., 2015).In layers where fluids are allowed to flow, drained parameters are needed, and the overpressure resulting from CO 2 injection is given by the fluid flow simulations.Where this is not the case (in closed faults or seals), undrained parameters are used.Fluid flow was not modeled in such layers and the geomechanical simulator is then loaded with a null overpressure.
Identical geomechanical properties were chosen for all fluid flow scenarios according to data derived from the available well logs, lab tests and literature (Tab.6).Two sets of parameters are 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 when the facies is in a layer where fluids are not allowed to flow.
Commonly used boundary conditions were applied at each side of the model: the top surface is a free surface, the bottom surface is assumed to have a null vertical displacement, and side surfaces (located far from the perturbation solicitation so as not to influence the results) are assumed to be subjected to constant regional stresses: in the modeling process, side surfaces are thus assumed to be free surfaces, as the initial effective regional stress is added in the post-processing as described below.
Fault integrity, versus fault potential reactivation, is assessed by comparing the fault stress state with a stress envelope that defines the limit between elastic and inelastic behavior.This is done in two steps: 1. the overpressure resulting from fluid flow simulation is used by the geomechanical simulator to compute the stress variations induced by CO 2 injection; 2. these stress variations are added to the initial stress state to derive the current stress state.This stress state is compared with a Mohr-Coulomb criterion in order to assess the potential damage of the fault.More precisely, the oriented distance of the stress state to the Mohr-Coulomb criterion is computed (Fig. 13): -Crit INI is the oriented distance of the initial stress state to the damage criterion; if negative, the initial stress state is below the damage criterion and the fault is stable; -Crit END is the oriented distance of the final stress state (after CO 2 injection) to the damage criterion; if negative, the final stress state is below the damage criterion; -DCrit is the difference between Crit END and Crit INI .
A positive value of DCrit indicates that the stress state is closer to the criterion after injection.A negative value of DCrit indicates that the stress state is farther from the criterion after injection.
Crit INI and Crit END were computed on the faults of the South Adriatic model and it was observed that these damage parameters remained negative for all tested configurations and time steps.This means that no damage is expected from CO 2 injection.
Since the effect of CO 2 injection is very small on the distance to the Mohr-Coulomb criterion, Crit END and Crit INI appear almost the same and the differences due to the different scenarios and the different fault states are hardly visible on Crit END (Fig. 14).Instead, Figure 15 presents, for different scenarios, the DCrit value, which is the evolution of the distance to the Mohr-Coulomb criterion after 10 years of injection.These results are displayed on the plane of the Rovesti fault, which is the closest one to the injection well, and consequently the most likely fault to be reactivated.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 [À8 bar; 15 bar], whereas the range of Crit INI and Crit END is around [À4.5 bar; À377 bar] so that CO 2 injection has a low impact on the potential damage based on the Mohr-Coulomb criterion.This, of course, is strongly related to the injection scenario (in particular, the position of injection wells and rate of injected CO 2 ).
When faults are considered open (Fig. 15 right), CO 2 injection leads to a pore pressure increase both in the reservoir and in the fault.This increase calculated by the fluid  flow simulator is used as mechanical loading in the geomechanical simulator that computes the resulting effective stress state.The value of DCrit is positive, which indicates that CO 2 injection makes the effective stress state closer to the Mohr-Coulomb criterion.Nevertheless, the damage remains below the Mohr-Coulomb criterion, which means that there is no risk of fault reactivation.It can be observed that the LOGS scenario leads to a smaller impact on the distance to the Mohr-Coulomb criterion than the OAb scenario.This is due to the poor quality of the seal in the LOGS scenario, that leads to a smaller increase in the pressure, and consequently a smaller effect on the stress state than those calculated for the OAb scenario.The two injection wells scenario shows a DCrit footprint quite close to the single well scenario due to the small difference in overpressure in the fault between both scenarios.However, the DCrit value for the two injection wells scenario is larger: this is due to the position of the wells, that are aligned perpendicularly to the fault, so that their interaction enhances the overpressure magnitude along the fault.Other well positions (such as parallel to the fault or more widely spaced) would have led to different overpressure in the fault, and so, a different impact on the fault behavior.
The closed fault scenario (Fig. 15 left) shows DCrit values quite different from the open fault scenario.It is actually tricky to compare the open and closed faults scenarios.Even if the same amount of CO 2 is injected, the variation of the pore pressure due to CO 2 injection is actually quite different in the two scenarios due to the possibility of fluids flowing or not in the fault.When faults are closed, the pore pressure increase in the reservoir is smaller than in the case of open faults.In the faults, fluids are not allowed to flow so that the increase in pore pressure in the fault is only due to a geomechanical effect associated with its undrained behavior.
Contrary to the open fault scenario, CO 2 injection induces a compression of the fault at the reservoir level.This makes the distance to the Mohr-Coulomb criterion greater at the reservoir level.
It is noted that the initial stress state is another influential parameter of the results that impacts both the normal and shear stress on the fault plane.The reader can refer to Baroni et al. (2015) for further details.injection scenarios.All simulations show that the fault behavior remains below the Mohr-Coulomb criterion, which confirms the integrity of the faults.However, these results have to be handled with care due to the large uncertainties that result from the very sparse available data.Some improvement should also be applied to the coupled fluid flow and geomechanical modeling.In particular, faults were considered here as a homogeneous material; however, some heterogeneities should be introduced by informing the relevant properties according to the shale gouge ratio (Adeoti et al., 2009;Yielding, 2002), for instance.In addition, a fully coupled modeling should be performed, in particular to validate the difference in the results observed when comparing closed fault and open fault scenarios.
The site characterization applied to the Southern Adriatic offshore is an example of the approach that could be adopted, in a first stage, when investigating this kind of reservoir formation where sparse information is available.This study also represents the first attempt to apply the full chain in the CO 2 storage site characterization in the Southern Adriatic offshore where major Italian CO 2 point sources are located.The study area represents a good opportunity to apply CCS on the scale of a demo project considering, however, that a real feasibility study would be required to access additional data or new data acquisition.
Figure 1 a) Location map of the investigated area; b) structural scheme of South Italy with the location of the main seismogenic tectonic lineaments (modified after Di Bucci et al., 2009; Angeloni, 2013).
Figure 2 a) Seismic view of the Grazia Anticline; b) interpreted lithological section at Grazia well.

Figure 3
Figure 3Location map of seismic lines and well sites.

Figure 6
Figure 6 Figure 7 Structural model: a) regional faults affecting the top reservoir; b) internal view of the structural model; c) 3D grid, with the number of cells in x (cell size in I direction is 500 m), y (cell size in J direction is 500 m) and z (cell size in K direction ranges from 1 m to 1 300 m) directions; d) vertical layering: 4 main zones delimited by the interpreted seismic horizons.

CO 2 .
Fourteen scenarios were run simulating a period of 30 years, including 10 years of CO 2 injection, and corresponding to: -LOGS, CORE, BIBLIO and OAb matrix scenarios with either closed, intermediary or open faults and one injection well; an OAb matrix scenario with either closed or open faults and two injection wells, each one with an injection rate half that of the single scenario one.
Figure 8Vertical proportion curves for the zone corresponding to the Oligo-Miocene succession of the caprock in the volumetric grid.
Figure 11 Wellbore overpressure (bar)evolution with time for the LOGS, BIBLIO, OAb (one well) and OAb (two wells) matrix scenarios: open/closed fault cases.

Figure 12 3D
Figure 12 3D overpressure (bar) at 10 years for the LOGS and OAb cases one or two wellsopen or closed fault scenarios at the end of the injection period (10 years): sections in the Rovesti fault plane and perpendicular to it.

CONCLUSIONSA
"one-way" coupling methodology between fluid flow modeling and geomechanical modeling was tested on the Grazia site, to assess the integrity of the faults in the context of CO 2 injection.Faults are integrated in the model by means of cohesive finite elements that allow one to take into account stress inside them.Post-processing was performed to assess potential damage on the faults according to the Mohr-Coulomb criterion.Different scenarios were run, according to different petrophysical parameters, different assumptions on the fault behavior and different Figure 14 Final distance to the Mohr-Coulomb criterion of stress state on the Rovesti Fault for the normal faulting initial stress regime.Different matrix scenarios are considered: a) LOGS, b) OAb considering one injection well and c) OAb considering two injection wells.Different states of faults are also considered: (left) closed fault scenario, (right) open fault scenario.Zoom around the injection well (lateral extend of 26 km; vertical extend of 5.6 km from the seabed).
Figure 15 Evolution of the Mohr-Coulomb criterion, Crit, on the Rovesti fault, as the difference between the Mohr-Coulomb criterion, Crit INI , for the initial stress regime (assumed as normal faulting) and the stress regime calculated after 10 years of CO 2 injection for different matrix scenarios: a) LOGS, b) OAb considering one injection well and c) OAb considering two injection wells.Different states of faults are considered: (left) closed fault scenario, (right) open fault scenario.Zoom around the injection well (lateral extend of 26 km; vertical extend of 5.6 km from the seabed).

TABLE 5
Geological formations and associated facies in the geological model