Send Orders for Reprints to Reprints@benthamscience.ae How to Approach Subsidence Evaluation for Marginal Fields: a Case History

This paper presents the evaluation of the subsidence potentially induced by underground storage of natural gas in a marginal depleted field located in Southern Italy. The critical aspect of the study was the lack of data because economic and logistic reasons had restricted data acquisition at the regional scale to perform a geomechanical study. This limitation was overcome by accurately gathering the available data from public sources so that the geometry of a large-scale 3D model could be defined and the formations properly characterized for rock deformation analysis. Well logs, seismic data and subsidence surveys at the regional scale, available in open databases and in the technical literature, were integrated with the available geological and fluid-flow information at the reservoir scale. First of all, a 3D geological model, at the regional scale, incorporating the existing model of the reservoir was developed to describe the key features of a large subsurface volume while preserving the detail of the storage reservoir. Then, a regional geomechanical model was set up for coupled mechanic and fluid-flow analyses. The stress and strain evolution and the associated subsidence induced in the reservoir and surrounding formations by historical primary production as well as future gas storage activities were investigated. Eventually, the obtained results were validated against the measurements of ground surface movements available from the technical literature for the area of interest, thus corroborating the choice of the most critical geomechanical parameters and relevant deformation properties of the rocks affecting subsidence.


INTRODUCTION
Subsidence is a well-known phenomenon in Italy.Regions such as the Po Plain have been the focus of extensive research over the last decades due to the major social and economic impact of subsidence on highly urbanized areas [1][2][3].It is common knowledge that subsidence can be caused both by long-term natural processes and by anthropogenic activities, but the effects occur at different time and spatial scales.The main components of natural subsidence can be classified as: tectonic loading, sediment loading, sediment compaction and post-glacial rebound; anthropogenic subsidence is mainly due to underground fluid exploitation, namely diffuse massive water utilization and local hydrocarbon production.Underground gas storage mainly causes a cyclic ground movement (i.e.subsidence and rebound) related to seasonal gas withdrawal and injection operations.Generally speaking, ground surface movements induced by a storage cycle are much smaller than those induced by reservoir primary production and groundwater withdrawal.Nevertheless, the potential impact of underground storage activities on existing constructions and infrastructures needs to be assessed to comply with more and more stringent recommendations from governmental authorities and -just as importantly to gain social acceptance.Subsidence needs to be evaluated *Address correspondence to this author at the Dipartimento di Ingegneria dell'Ambiente, del Territorio e delle Infrastrutture (DIATI), Politecnico di Torino, Corso Duca degli Abruzzi 24, 10129 Torino, Italy; Tel: 00390110907610; Fax: 00390110907699; E-mail: vera.rocca@polito.it both in terms of magnitude of the vertical displacement and extension of the involved area.
When a depleted hydrocarbon reservoir is converted into gas storage, the geological, structural, petrophysical and fluid dynamic properties of the reservoir formations are usually well known because they were investigated during the primary production phase.Generally, static and dynamic studies are available, which standardly include geological and fluid flow numerical models.Even if this information is key in the construction of the geomechanical model, it must be extended at the regional scale and integrated with geomechanical characterization of both intact rock and fractures/faults.Because the reservoir and the undisturbed rock around it remain connected, pore pressure perturbations in the reservoir induced displacements at the reservoir boundaries that caused deformations and stress changes also in the rocks around the reservoir [4,5].As a consequence, the domain of the mechanical investigations is not limited to the reservoir but includes all the rock volume affected by the stress-strain variations induced by periodic gas withdrawal and injection.In case of subsidence analysis via a 3D numerical modelling approach the rock volume which has to be characterized and analyzed from a geological, fluid flow and mechanical standpoint comprises the formations surrounding the storage and the overburden up to the surface.However, the costs involved in a detailed characterization of such a large volume can compromise the economic feasibility of a storage project; furthermore, data acquisition can be difficult due to neighboring exploitation/storage permissions awarded to different oil companies.These critical issues in obtaining the data, needed to set up a reliable numerical model for geomechanical simulations, can be partially dealt with by gathering public information available in the technical literature and from national and local governmental agencies.In Italy, the database of the National Mining Office for Hydrocarbons and Georesources (UNMIG) comprising thousands of well logs and 2D seismic sections provides valuable data that can significantly assist regional scale modeling.
Through the case history of a marginal field located in Southern Italy, this paper shows how to perform a reliable subsidence analysis despite the lack of direct data acquired at the regional scale for this specific purpose.First of all, the geological background of the area under investigation was delineated.Then, the reservoir geological model was extended at the regional scale.To this end the workflow for a static model definition at the reservoir scale is described in this paper.Subsequently, the regional geomechanical model for coupled mechanic and fluid-flow analyses was defined.Eventually, the effects of stress and strain evolution due to historical primary production and future gas storage activities were investigated in terms of induced subsidence.

Geological Background
The studied area is located in the Basilicata region, Southern Italy, between the provinces of Matera and Potenza (Fig. 1).In the 60's, a large number of wells were drilled in the area and several seismic sections were acquired for hydrocarbon prospecting.In fact, an early exploration phase revealed the presence of a thick (>2000 m) sedimentary sequence constituting the filling of a basin later known as Bradanic trough (see Lucanian basin), that is the Pliocenepresent day southern Apennines foredeep.

Fig. (1).
Simplified geological map of Southern Apennines and location of the study area at the limit of the Apennine frontal thrust.The boundary of the regional geological model is depicted by the black dashed line.

Tectonic and Sedimentary Evolution
The tectono-sedimentary evolution of the area, related to the Apennines formation, began during the Early Pliocene, when the compressional tectonic phase started [6].As a consequence of the Apennines frontal thrusts progression, the basin is characterized by a general north-east migration of its axes of subsidence and related depressions [6].The basin is characterized by an instable internal margin with a tendency to be strongly uplifted and a subsiding external margin progressively involving the inner foreland zone.Two main sequences can be distinguished for both the internal and external margins, namely the 'Appeninic' and the 'Murgian (Apulian)'.The Apenninic sequence is characterized by the presence of both the Allocton, a chaotic geological complex of pre-pliocene age, interposed between the plio-pleistocene terrigenous deposits and consequent to the middle pliocene tectonic phase [7] and the transgressive deposits on the carbonatic substratum.The Murgian sequence is only characterized by transgressive plio-pleistocene deposits on the carbonatic substratum [8].The pre-pliocene substratum is constituted by miocene or older deposits [8].From the Eocene to the Lower Miocene, the area currently occupied by the basin acquired the characteristics of instable foreland [9].The Eocene is characterized by a carbonatic deposition associated to a distensive tectonic phase relevant to both the foredeep basin and the Apulian platform.It is assumed that a generalized emersion of the area occurred at the end of the Eocene, because low to no deposition took place during the Miocene, while the Oligocene is completely absent [10].The pliocene sequence lies transgressive on the pre-pliocene substratum with layers being younger towards the southeast [8].

Structural Framework and Stratigraphy
The main structural features characterizing the basin are presented in Fig. (2) and can be summarized as follows [10]: 1. a NW-SE elongation with asymmetric transversal profile; 2. a depocentral area affected by an extensional regime.The inner eastern margin is characterized by a regional highgradient structural ramp derived by the direct faulting of the meso-cenozoic limestones of the Apulian Platform, while the outer eastern low-gradient zone, also known as 'pre-murgian shelf' [11], is characterized by a horstgraben structure.The regional ramp is subject to subduction underneath the Apennines chain; 3. a western margin defined by the accretional frontal wedge of the Apenninic thrusts (Allocton) propagating towards NE and involving the pliocene deposits.
The stratigraphic setting of the Lucanian basin was deduced and largely studied on the basis of well and outcrop data [8,12].Three main stratigraphic units could be identified [12]: • the Pre-pliocene substratum • the Plio-pleistocene deposits

• the Allocton
Confirming the results of Crescenti [8] and Balduzzi [12], the following facies sequence can be distinguished in the Pre-pliocene substratum: Based on Balduzzi [12], the Plio-pleistocene deposits can be divided into three main stratigraphic intervals: 1. a top interval (post-turbiditic phase), characterized by shaly deposition on the lower part and sandy deposition at the top of the sequence and representing the filling up of the foredeep, dated Pleistocene and characterized by variable thickness up to 1 kilometer.Finally, the Allocton is represented by pre-pliocene and pliocene formations and it is made up of various lithology and chaotic texture.Its thickness increases towards WSW (at least 2 kilometers at well Masseria Rigirone 1) and it was reached by all the wells drilled on the western margin of the Lucanian basin [12].

Reservoir Description
The study presented in this paper was performed on the Grottole-Ferrandina reservoir and all the available data for reservoir modeling was provided by GEOGASTOCK S.p.A., holder of the concession.The reservoir lies below the broader area of the towns of Grottole and Ferrandina of the Matera province in southern Italy.The names of the drilled wells depend on their vicinity to the two nearby towns.
The reservoir is located at an average depth of 900 meters s.s. and it lies above the horst resulting from the Upper Pliocene tectonic regime [12].The reservoir is hosted inside the Pliocene and Pleistocene clastic formation known as 'Argille del Santerno'.This formation is mainly composed of shales and minor sandy intercalations 'draping upon the pre-pliocene substratum high' and formed in response to syn-sedimentary tectonics and/or differential compaction [13,14], (Fig. 3).These sand drapes constitute a peculiar kind of hydrocarbon (stratigraphic) trap.The Grottole-Ferrandina reservoir is characterized by both stratigraphic and structural features: a facies heterotopy from sand to shale characterizes the northern and eastern closures of the reservoir whereas direct faulting limits the reservoir along the southern and western boundaries as a result of tectonic activity.Faulting also occurs inside the reservoir, defining an internal compartmentalization of two blocks.Eventually, the cap rock of the reservoir is constituted by the main shaly portion of the 'Argille del Santerno' formation.

REGIONAL VERSUS RESERVOIR GEOLOGICAL MODELING
During the lifetime of a hydrocarbon reservoir 3D numerical static and dynamic reservoir models are built and periodically updated when new information (e.g.new well logs, well tests, static pressure surveys, production history) becomes available.As a consequence, a good knowledge of the system in terms of the geological structure and layering, the petrophysical and the fluid-flow properties is generally achieved at the end of reservoir life, yet this knowledge is typically limited to the reservoir volume.
When a depleted field is considered for conversion into gas storage, geomechanical analyses are often required in order to assess the potential impact of the storage activities on the ground surface.This implies that a geomechanical model incorporating all the main geological features should be set up, the dimensions of which are necessarily much larger than those of the reservoir volume.Thus a new 3D geological model has to be constructed for geomechanical simulation purposes; as a consequence, regional geological information should be gathered and incorporated into this new model.
In the following, both the reservoir and the regional geological models are described: the reservoir model that incorporates the reservoir characteristics in detail and the regional geological model that describes the area at a larger scale from both the structural and stratigraphic viewpoints.At the end, an integrated model, which will be used later on for the geomechanical simulations, is created combining all the information available at the reservoir and regional scale.

Dataset
Available data for regional modeling was collected from the online public national database (UNMIG database) [15].3).Schematic geological section across the North Bradanic Trough.Some oil and gas fields located in the area are projected [14].In particular, the Grottole-Ferrandina reservoir is shown to the north-east: the lenticular (lobate) shape of the sandy intercalation and the drapping upon the pre-pliocene substratum high can be observed.
The database comprises composite logs (scale 1:1000) and seismic sections acquired by oil and service companies since the beginning of exploration activities in Italy.
The retrieved data mostly contains composite logs with well information (i.e.: wellhead coordinates, deviation surveys, rotary table, etc.), lithologic logs (typically gamma ray or spontaneous potential registrations) and resistivity logs for formation fluid identification; in addition to the lithological information deriving from the cuttings, the main stratigraphic units are also reported.The lack of core data as well as of density and neutron logs was a main issue to finalize the lithological and petrophysical characterization of the regional model.In particular, the available dataset was not useful to distinguish terrigenous from carbonatic rocks.It is well known that gamma-ray logs can consistently identify sandstones (low gamma ray) and mudstones (high gamma ray), providing a measure of the clay volume [16].However, also carbonate rocks cause low gamma ray responses and they can be distinguished from sandstones only based on their greater density and acoustic velocity on porosity logs.
Seismic sections were available in raster format, occasionally containing structural and stratigraphical interpretation information.Well trajectories were in most cases overprinted, but as no velocity information was available, a rigorous time-depth conversion could not be performed.Generally, well and seismic data could only be correlated qualitatively or by considering the main tectono-stratigraphic events based on regional geological knowledge.
As the original data was collected from various sources (i.e.oil/service companies), which often adopt different standards and use different levels of detail in their reports, the UNMIG database is highly non-uniform.As a consequence, the need for a consistent dataset led to a significant reduction of the usable data.
Well data inside the Grottole-Ferrandina concession area were kindly provided by GEOGASTOCK S.p.A.

Well Data
Composite logs were available for all the 57 wells drilled in the Grottole-Ferrandina concession area.In particular, 18 wells in the Ferrandina area (FE in Fig. 4) and 39 wells in the Grottole area (GR in Fig. 4).Composite logs available from UNMIG database, to be used for (regional) modeling purposes, were selected on the basis of three requisites.The first one was the well location so as to adequately cover the area of interest.The second was the level of detail, which is essential in the case of lithological or litho-stratigraphic observations.The last one was the homogeneity of the final dataset, which is crucial when performing well correlation in order to avoid inconsistent results due to different approaches in the litho-stratigraphic description.It was found that well profiles coming from eni (formerly AGIP) oil company were preferable among all the others because they typically satisfied the preset requisites.In the end, the selected wells were: Castelluccio 2, Ferrandina 15, Grottole 2, Grottole 17, Masseria Caniuccio 1, Masseria Rigirone 1, Monte S. Vito 2, Pisticci 3, Pizzo Corvo 1, Pomarico 4, Pomarico 7, Salandra 2, Serra d'Olivo 1 (Fig. 5).

Seismic Lines
Various seismic sections available in the UNMIG database were acquired in the studied area.The available seismic sections and those selected for the model construction are shown in Fig. (6).After a preliminary quality check, the selection of the seismic sections was operated according to their proximity to the studied area, their orientation (only the ones perpendicular to the axis of the basin were considered) and their position with respect to the main geological features (such as the Allocton, the carbonates, faults, etc.).The  4).Regional and reservoir models created for this study.On the left, the regional model is presented and on the right, the reservoir model along with their average dimensions to emphasize the different size of the two grids.The wells located inside or very close to the reservoir area are also shown.

Fig. (5).
Map of the wells used in the construction of the regional geological model.The boundary of the regional model is shown by the black line.

Fig. (6). Map of the available seismic sections from UNMIG [15]
from the broader study area.The red lines represent the seismic lines used in this study.The yellow region in the middle of the figure denotes the extension of the reservoir in the area, while the black line the boundary of the regional geological model.sections aligned vertically to the buried Apennines frontal thrust show the geometrical relationship between the Allocton, the foredeep deposits and, partially, the foreland, providing an insight into the geological structure in the proximity of the area to be modeled.
The seismic sections within the model boundaries (Fig. 7) clearly show the characteristic geometry of the Allocton and its relation with the underlying carbonates.It can be observed that the Allocton front appears more advanced to the north (section 1, Fig. 7) and to the south (section 3, Fig. 7), while it appears retreated in the middle (section 2, Fig. 7).All these features were represented in the 3D regional geological model.Furthermore, the carbonatic platform is strongly faulted with some of the faults also penetrating the overlying geological layers -including the Allocton.
The seismic sections, namely MT-558-95 and PZ-692-94, acquired 20 km to the north and 10 km to the south of the boundaries of the studied area, respectively, were selected for a further fine-tuning and validation of the regional geological model, with a particular focus on the Allocton geometry and on the stratigraphic trend of the plio-pleistocene sequence.
The MT-558-95 profile (section 5, Fig. 8) was published with an overlaid interpretation [17].It is characterized by low resolution, thus only the main structural and  stratigraphic features can be recognized.The Allocton wedge is clearly identifiable (but the internal geometry cannot be deduced) as well as the main stratigraphic trend of the foredeep deposits.The typical faulting pattern of the underlying Mesozoic carbonates, caused by the load of the Allocton on the foreland platform, could not be identified.
The PZ-692-94 seismic profile (section 4, Fig. 8) shows a remarkably higher resolution than the other sections.The interpretation was performed according to Lazzari's approach [10] applied to the seismic profile PZ-445-81, located in the northern part of the basin at about 30 km from the northern boundary of the studied area and oriented almost perpendicularly to the axis of the basin.Along the PZ-692-94 profile the overall geometry of the Allocton as well as the main internal thrusts can be identified.A strong reflector lying below the thrusts system can be easily recognized.An attempt to date this surface was performed, based on the integrated analysis of the composite log stratigraphy and on the trajectory of well S. Chirico 3; the reflection could be caused by the lithological variation occurring in the last 25 meters of the sequence intercepted by well S. Chirico 3, from the thick pliocene shales to the underlying Eocene-Miocene carbonatic sequence.Thus the surface should be reasonably dated as top or near top of the Miocene.This would locate the Plio-Pleistocene sequence lying between this reflector and the Allocton.Because no relevant features could be identified below the reflector, the geometry of the underlying Mesozoic carbonate platform could not be defined.Furthermore, a shallow sequence sealing the outer (eastern) thrust fault can be identified on the upper portion of the profile.Conversely, in the inner (western) zone, the seismic signal appears not to change from the deeper zone up to the surface, hence it is not possible to identify any other element except for the Allocton thrust system.

Regional geological model
The modeled area, shown in Fig. ( 4), extends approximately 30 km in the NW-SE direction and 35 km in the SW-NE direction; in the vertical direction 2000 meters of sequence was considered.Geologically, the area encompasses all the main structural and stratigraphic features characterizing the Lucanian basin, as described in paragraph 2.2.

Structural Model
The geometry of the Allocton was defined based on the stratigraphic correlation between wells Masseria Rigirone 1, Masseria Caniuccio 1, Pisticci 3, Serra D'Olivo 1 and Salan-dra 2. The interpolation of the results led to the definition of the overall geometry of the 'Allocton' wedge (Fig. 9).The SW-NE sections of the model qualitatively match the interpretation of seismic profiles previously analyzed.With respect to the main NW-SE orientation of structural features [10], a relevant deviation was observed in the southern part of the model, where the main orientation was defined WNW-ESE.This deviation is due to the absence of the Allocton at well Pizzo Corvo 1.
The geometry and orientation of the horst-graben structure of the meso -Cenozoic carbonatic platform was defined, coherent with the Allocton frontal thrust geometry, by the analysis of the well data in the reservoir area.The identification of different blocks was performed via well correlation and the main faults were mapped by measuring the distance of each well from the top of the carbonate sequence.This approach led to the identification of four main blocks.From west (inner part of the basin) towards east, they are defined as follows (Fig. 10): • Block 1: intercepted by wells Masseria Caniuccio 1 and Pisticci 3.
• Block 3: intercepted by wells Pomarico 4, Pomarico 7, Ferrandina 15, Grottole 17 (public data) and by all the reservoir production wells.It corresponds to the Grottole-Ferrandina horst previously described by Balduzzi et al. [12].The integration between reservoir and regional data provided a detailed definition of this structure.
• Block 4: intercepted by wells Monte San Vito 2 and Grottole 2.  An additional block (Block 0) was modeled westward to Block 1 based on the resulting overall Allocton and substratum geometry.However, no well data is available in the area to confirm this hypothesis.
The lack of seismic data and the scarcity of well data did not allow identification of the fault system in the pliopleistocene sequence.

Stratigraphic Model
The stratigraphic modeling was performed on the basis of the well data and selected bibliography [12].The geological conceptual model, represented by the various seismic sections described in the selected bibliography, was integrated in the model.The well correlation was performed in order to reconstruct the geometrical relationship among the main stratigraphic and structural features characterizing the area.Integrating the results of Balduzzi et al. [12], all the main chronostratigraphical subdivisions of the sequence, identifying the main geological phases, were correlated and, eventually, the overall 3D stratigraphic setting was defined.
Two main stratigraphic sequences were observed on selected composite logs depending on their location in the basin (see paragraph "2.1.Tectonic and sedimentary evolution").The typical stratigraphy of the western wells (i.e.: Masseria Rigirone 1, Masseria Caniuccio 1, Pisticci 3, Serra D'Olivo 1, Salandra 2) is characterized by the presence of the Allocton.Considering the basin physiography, the sequence intercepted by these wells is representative of the inner or 'apenninic' margin of the basin.Conversely, the rest of the wells (i.e.: Castelluccio 2, Pomarico 4, Pizzo Corvo 1, Pomarico 7, Ferrandina 15, Grottole 2, Grottole 17, Monte S. Vito 2) located in the eastern part of the basin show the absence of the Allocton.Here the sedimentary sequence belongs to the 'murgian' margin.

Reservoir Geological Model
The reservoir model, shown in No seismic data was available for the reservoir area.The reservoir zone was intercepted by 24 wells.The spacing of the wells is quite regular, which offers a good level of reliability for the reconstruction of the stratigraphic surface geometry and for the petrophysical parameter estimation.

Structural Model
The main structural features identified by the analysis of the well data are the faulting systems, typical of the carbonatic pre-pliocene deposits and resulting in the horstgraben structure previously described, and the clastic pleistocene deposits, hosting the gas bearing formation.
On the basis of well Grottole 4, which did not encounter the gas-bearing formations, the north-eastern boundary of the reservoir was accurately identified.The west-eastern boundary was set based on the wells drilled along the margin.This reservoir closure is partially confirmed by the presence of a minor block, targeted by wells Grottole 31 and Grottole 32 separated from the reservoir block to the north.From the analysis of the well data (logs and production) the fault separating the reservoir block from this minor block also acts on the reservoir sands, deposited during Early Pleistocene.This is the only plio-pleistocene structural feature that could be identified from the available data.

Stratigraphic Model
The reservoir is hosted inside two clastic layers, characterized by a limited areal extension (7 km NW-SE, 2 km NE-SW).Lithologically, these two layers are constituted by se-quences of sand and shale, deposited during the Calabrian (Early Pleistocene), when a mainly shaly deposition occurred over the area [12].A lithostratigraphical correlation was performed to define the reservoir geometry and continuity.In particular, the two layers are two isolated sandy lenses separated by a continuous shaly interlayer acting as a permeability barrier.The layers clearly pinch out to the north, close to well Grottole 27.

INTEGRATED 3D REGIONAL AND RESERVOIR MODEL
The main faults (Fig. 11) were built based on the structural analysis results at both the regional and reservoir scale.Their inferred geometry was mainly constrained by well data thus their final geometrical accuracy is subject to well spacing.
The main stratigraphic surfaces were built on the basis of the well correlation results, a qualitative analysis of the seismic sections and a conceptual geological model [12].

3D Grid Construction
The 3D model grid was built so as to match the multiple purposes of the model, which were: • represent the geological setting at both the regional and reservoir scales; • simulate fluid flow at the reservoir scale; • assess subsidence through geomechanical simulations at the regional scale.
The geological reliability of the model at both regional and reservoir scales mainly lies on the capability of representing the geometry of the structural and stratigraphic features.Thus, special effort was made to define the structural setting of the meso-cenozoic carbonatic platform and the frontal thrust of the Allocton wedge as well as the stratigraphic setting of the plio-pleistocene filling of the basin (both its internal pattern and its relationship with the prepliocene substratum and with the Allocton wedge).
The detail of the geometrical reconstruction is strictly linked to the geometry, orientation and dimensions of the cells constituting the 3D grid.Furthermore, cell dimensions and regularity strongly affect the functionality of both the fluid flow and geomechanical modeling.The possibility of managing cell geometry and orientation is crucial for a realistic reconstruction of the structural features, as cells can be oriented and partially distorted to best fit the fault surfaces.Of course, the smaller the cell dimensions, the more realistic the model geometry.However, the cell dimensions are inversely proportional to the total number of cells, thus to the simulation run time.The adopted strategy for building the grid was based on the following: a) the application of a gridding technique which allows a different level of detail in the description of the reservoir zone and the surroundings, depending on the level of characterization and target of the study.
b) the constraint of the grid geometry to fault orientation so as to satisfy the need for a highly realistic-shaped fault system.This operation was carefully monitored insofar as only minor deformation to cells geometry should be accepted.
The final geological grid capturing the geometries identified in the seismic sections and obtained from bibliographic information is presented in Fig. (12).

GEOMECHANICAL MODEL
The geological and structural model previously defined at the regional scale represented the starting point for the geomechanical model construction.The grid geometry, i.e. extension and discretization, was further optimized for geomechanical analysis purposes.Input data was collected from the technical literature and UNMIG as well as from the performed geological analysis.The intact rock and fracture/fault systems were then characterized in terms of petrophysical and mechanical properties.The populated numerical model was initialized and then run to evaluate subsidence due to future gas storage activities, according to different operational scenarios.

Grid
A reliable geomechanical grid not only includes reservoir geometry, but also its side-, over-and under-burden formations.Usually, the investigation domain is identified by the reservoir, by its surrounding formation and, especially in the case of subsidence analysis, by the overburden up to the surface; in fact, this volume is the most affected by the stressstrain variations induced by production/injection operations.The extension to the under-burden has the main purpose of setting the boundary conditions far enough from the reservoir and of avoiding numerical stability issues: generally speaking, the numerical stability of a model can be significantly improved by adding stiff bottom cells [18].As a consequence, the optimization of the model grid has to address two critical aspects: global extension and discretization.An appropriate global model extension ensures undisturbed boundary conditions while maintaining acceptable computational time.According to preliminary analyses, a model areal extension of 33x25 km 2 with a vertical thickness of 5 km from the surface (Fig. 13) turned out to be suitable for a correct description of subsidence evolution, ensuring undisturbed boundary conditions.

Fig. (13). Areal and vertical extension of the geomechanical model.
A series of sensitivity analyses were performed to evaluate the potential impacts of the grid block dimensions on simulation results.The gridding of the reservoir zone was the same as in the calibrated dynamic model so as to facilitate data exchange between mechanic and dynamic models and to investigate with a high degree of reliability the volume subject to the largest stress gradient.Conversely, more coarse meshes were adopted for the description of the remaining model volume, marginally affected by the phenomena under analysis.In the reservoir region (10x5 km 2 ) a square 100x100 m 2 meshing was adopted; the grid was progressively coarsened towards the boundaries of the model according to a geometric progression (Fig. 14).The total number of cells was 117 x 69 (I and J directions, respectively).The vertical layering reflected the sedimentary sequence recognized from the geological and structural sections of the region.The cap rock and the overburden volume were further refined in order to adequately describe the geomechanical model response.The model consisted of 32 layers.
This gridding approach provided a good accuracy in the model volume affected by subsidence phenomena while limiting the total amount of cells.Furthermore, because reservoir dimensions are limited, the total number of cells in the geomechanical model and, consequently, the computational time do not represent a critical issue for the case under investigation.

Input Data
The required mechanical input for populating the model can be divided into: initialization parameters, deformation and strength parameters.
Due to the lack of laboratory tests and/or relevant log data, the deformation and strength parameters as well as the initialization data were defined on the basis of the authors' experience and of similar cases found in the literature.For example, the work developed by Montone et al. [19,20] for the Italian territory to improve and enrich the world stress map was a valuable source of information on the far field stress state of the zone under investigation.The deformation parameters were defined according to a basin-scale compressibility law derived by Teatini et al. [1] under the assumption that the system was isotropic.This empirical correlation was obtained from in situ deformation measurements via the radioactive markers technique.The investigated formations were silty to fine grained sandstones, which are similar to the lithologies investigated in this work.The law provides the vertical uniaxial compressibility, c M , as an exponential function of the vertical effective stress, v : c M = 1.3696 10 -2 v -1.1347 (1) According to this relation, the static elastic modulus, E s , adopted in the model characterization can be calculated as a function of depth, as shown in Fig. (15).Equation ( 1) holds true for rock compression in virgin loading conditions (first loading cycle), while rock expansion is controlled by the value of c M ' in unloading/reloading conditions (second loading cycle and the following).The first and second loading cycles can be assimilated to in situ loading conditions during primary production and subsequent storage activity, respectively.From in situ marker data and odometer tests, Baù et al. [21] and Ferronato et al. [22] estimated that on average the unloading/reloading c M ' values are 1.8-3.5 times the loading c M values in a depth range between 1000 and 6000 m under hydrostatic conditions.According to this information and to the authors' experience, the elastic modulus was more than doubled during the simulation of the field reconstitution and storage phases.
In order to assess the impact of the elastic parameters on subsidence phenomena, a series of sensitivity analyses were performed: the system's response was evaluated under both hypotheses that the static or the dynamic Young's modulus, E S and E D respectively, could apply.A ratio of E D ~2.5E S was assumed [23,24].In the regional scenario under analysis, the lithological change between the shaly-sandy upper formations and the Cretaceous carbonate bottom sequence deeply affected the elastic response of the system: this was reproduced by significantly increasing the Young's modulus values in the bottom layers, for both the static and dynamic conditions, as shown in Fig. (16).
The strength parameters are related to the adopted failure criterion.Given that the Mohr-Coulomb criterion was adopted in this study [25], the cohesion and friction angle were required.Their values were defined according to the formation lithology and depth.Cohesion and friction angle of clay formations were gathered from available geotechnical laboratory tests on similar formations [26,27], while the strength parameters of the sand formations and basal carbonate formations were determined on the basis of the technical literature [27] and on the authors' experience.

Model Characterization
The geomechanical characterization of the system was performed defining 12 different classes reflecting the identified rock lithotypes, mechanical properties and depth.Each class was assumed to exhibit a homogeneous and isotropous geomechanical behavior.The analysis of the composite logs (scale 1:1000) used in the model construction showed three prevailing lithotypes: shales, sands and carbonates.The terms "shales" and "sands" refer to formations located at 800-900 m depth or more and, therefore, naturally subject to horizontal stresses (medium and minimum stress) greater than 8 MPa.For this reason, the mechanical characteristics of these formations are more similar to those of rocks than to those of soils.
The geomechanical parameters assigned to each class are provided in Table 1.
The classes called "Soil 1", "Soil 2", "Soil 3" and "Soil 4" identify shaly-sand formations with different compaction degrees according to their depths."Allocton" represents the clastic chaotic sediments of the Allocton wedge."Reservoir" characterizes the gas bearing Calabrian sequence."Cap rock" and "Interlayers" identify the shale sequences at the top and inside the gas bearing formation, respectively."Carbonates" correspond to the upper Cretaceous sequence at the bottom of the geological model."Underburden" represents the extra layers added at the bottom of the model for numerical purposes.
Only the faults directly affected by the pressure variations caused by gas withdrawal and injection were imported into the geomechanical model, namely the two faults located inside the reservoir.The stress-strain equilibrium of all the other faults described in the regional model were not perturbed by fluid movement, thus they were not relevant in the geomechanical analysis.Fig. (17) shows faults 'Reservoir 1' and 'Reservoir 2' which confine and compartmentalize the reservoir, respectively.These two faults were characterized according to the values summarized in Table 2.

Initialization Process
During the initialization phase, the original stress field and the initial pore pressure distribution of the undisturbed formations were determined as a function of depth, formation characteristics and saturation fluids.
According to the measured initial static pressure of the reservoir, a substantially hydrostatic regime was assumed: the pressure of each water-saturated cell (below sea level) was determined according to the hydrostatic gradient, while the pressure of the gas-bearing formation cell was calculated according to the gas gradient and the original gas-water contact depth.
The analyzed south Apennine area between the Allocton and the foredeep is characterized by an active extensional regime, according to Montone et al. [19,28] and Maggi et al. [29].On the basis of the few fault-planes data available   in the analyzed region, it was not possible to clearly determine whether the stress field in the foredeep is normal (like in the Apennine belt) or strike slip (like in the foreland).A series of sensitivities showed that the system is hardly affected by the far-field stress orientation.The far-field stress of the model was then assumed to be normal, thus the vertical stress was set equal to the maximum principal stress 1 .
In order to calculate the original stress field, the model was initialized by assigning a gravitational stress state; an anisotropic stress field in the horizontal plane was assumed ac-cording to the following ratios: and where ' V is the effective vertical stress and ' H and ' h are the effective maximum and minimum horizontal stresses, respectively.The direction of the maximum horizontal stress was 40° N [19].

Operational Working Conditions
The system's behavior during future storage activities was investigated according for different scenarios.The maximum operational pressure was imposed to be equal to,  but also greater than, the initial reservoir pressure so as to assess the benefits provided by delta-pressuring the storage.By definition, delta-pressuring (or over-pressuring) means operating the storage at a maximum reservoir static pressure which exceeds the initial formation pressure; deltapressuring conditions are typically reached toward the end of the injection period, thus are maintained for a limited period of time.This is a common option to enhance the storage performance, especially the working gas, at much reduced costs and virtually no hazards for a safe storage management [30,31].
During seasonal withdrawals the minimum tubing head pressure was set at 25 barsa (obviously this option requires the use of compressors).

SUBSIDENCE ANALYSES
A coupled mechanic and fluid-flow approach was adopted for analyzing the stress/deformation evolution induced in the reservoir and surrounding formations by primary production and future gas storage activities and the associated subsidence phenomena.The coupled approach is based on the integration of both rock mechanics and petroleum engineering principles and it assesses the effects that fluid flow phenomena have on porous media deformation and vice versa [32].The technical literature offers several theoretic methodologies to model the formation behavior with different degrees of coupling between rock deformation and fluid flow: from partially coupled to fully coupled methods.The fully coupled approach is based on the simultaneous determination of all variables, i.e. fluid flow and displacement calculations are performed together.In the partially coupled approach, the basic equations for multiphase porous flow and rock deformation are solved separately and sequentially and the coupling terms are iterated at each time step.In the one-way coupling technique, only geomechanical parameters are updated at each time step according to fluid dynamic reservoir behavior (i.e., pressure and temperature variations) defined via conventional reservoir models.The fluid flow simulation, instead, is not affected by the geomechanical behavior of the formation [33][34][35][36].
Considering the physics of the case-study and the involved phenomena, the one-way coupling approach was a good solution for predicting the impact of reservoir depletion on formation displacements.In fact, the adoption of a more complex coupled modeling technique would involve the need for additional (and not available) input data and supplementary time, without an appreciable improvement in the understanding of the phenomena under investigation.
According to the adopted one-way coupling methodology (Fig. 18), a pressure map at the geomechanical model scale is determined for each time step of analysis: the pore pressure disturbance induced in the gas-bearing layers is calculated by dynamic simulation, while the pressure distribution in the undisturbed formation remains constant and equal to the initial hydrostatic value (Fig. 19).Then, according to the selected elasto-plastic constitutive law, the geomechanical simulator determines the stress-strain state evolution due to the imposed pressure variation.When new geomechanical equilibrium is reached, for each time step, the induced deformations and displacements are calculated and so the subsidence cone can be estimated both in terms of extension radius and maximum vertical displacement.In the case history under investigation, the pore pressure variations in the reservoir were provided by a calibrated dynamic numerical model, which described the historical field dynamic behavior during primary production.Subsidence evolution was evaluated adopting a time-stepping of 10 years during the historical primary production, i.e. from 1962 to 2008 (time steps t 1 to t 5 ).Then, three storage scenarios were forecast: the planned refilling strategy from 2008 to 2015 was the same for all the examined cases (time step t 6 ); sub-sequently, the storage scenarios from 2015 to 2018 were defined according to different injection pressures.Three maximum operational pressures were forecast.In case 1, the injection pressure was set equal to the field discovery pressure (p i ); in cases 2 and 3 field overpressure conditions were simulated, with a maximum operating pressure equal to 110% p i and 120% p i , respectively.During the forecasted storage activities, the cyclic seasonal ground surface movements, i.e. subsidence and rebound, were evaluated at the end of each withdrawal and injection period.Fig. (20) shows the average reservoir pressure evolution in time during primary production and for the three forecast scenarios.
As already mentioned, in order to assess the impact of the elastic parameters on subsidence phenomena, the response of the system was evaluated using both the static and the dynamic Young's moduli.

RESULTS AND DISCUSSION
During historical production, the pressure reduction caused a progressive subsidence until a maximum vertical displacement was reached at the end of primary production (time step t 5 @ 11/2008).Fig. (21) shows the average pressure distribution at the top of the reservoir at t 5 , and the resulting subsidence at ground surface (in terms of vertical displacement and areal extension) both under the assumptions of static and dynamic elastic moduli.Fig. (22) shows cross sections of the subsidence cone along the longitudinal direction (section AB), at the minimum reservoir pressure when either the static or the dynamic elastic parameters were assumed.Eventually, the average pressure evolution in time and consequent vertical displacement evolution for a reference point on the ground surface, in both the static and dynamic elastic cases, is shown in Fig. (23).The assumption of static elastic parameters is the most critical in terms of maximum vertical displacement (7.5 cm) and subsidence cone extension (km 7 along the AB section, considering a minimum vertical displacement of 1 cm).Conversely, under the hypothesis that dynamic elastic parameters apply, a maximum vertical displacement of 1.8 cm is obtained with a subsidence cone extension of 4 km along the AB section.
For the three forecasted storage scenarios, the average pressure evolution and associated vertical displacement evolution versus time for a reference point, under both the static and dynamic elasticity assumptions, are provided in Figs.(24)(25)(26).Vertical displacement values are referred to the original elevation of ground surface, i.e. before reservoir exploitation.In case 1 (injection pressure = p i ), the position of the ground surface due to the cyclic seasonal subsidence  and rebound ranges between 7 and 5 cm below initial elevation for the static case and between 1.7 and 1.2 cm for the dynamic case (Fig. 24).Assuming that a maximum injection pressure equal to 110% p i could be reached (case 2), the ground surface position ranges between 6.8 and 4.7 cm below original elevation for the static elastic case and between 1.6 and 1.1 cm for the dynamic elastic case (Fig. 25).Finally, in case 3 (injection pressure = 120% p i ), the ground surface position ranges between 6.6 and 4.3 cm below original elevation for the static elastic case and 1.6 and 1 cm for the dynamic elastic case (Fig. 26).
The maximum vertical displacements and the subsidence extension along the AB section (considering a minimum vertical displacement of 1 cm) for each analyzed case, assuming either elastic or dynamic elastic parameters are summarized in Table 3.
The results obtained were compared with the measurements of ground movements available for the provinces of Matera and Potenza (Basilicata region).The accessible Synthetic Aperture Radar (SAR) data acquired and analyzed in the period 1997-2008 allowed the definition of a reference trend of ground movement evolution in the area under  investigation (Fig. 27) [37][38][39].A part from local landslides and other slope instabilities related to specific geological/structural issues, the average vertical movements of the ground surface due to natural processes and anthropogenic activities including gas production from the subsoil ranges between [-2; +2] cm for the monitored area (Fig. 28).As a consequence, the evaluation of subsidence phenomena based on the assumption that dynamic elastic parameters should be applied is the most realistic, whereas the simulation results obtained when static parameters were adopted provide more conservative scenarios.Furthermore, from a theoretical standpoint [40], the stress-strain behavior of the soils and rocks constituting the reservoir and the surrounding formations is strongly nonlinear and depends on the magnitude of the induced deformation of the porous media.When the system is subject to extremely limited deformations (i.e. 10 -3 m/m), such as in the case of storage activities, the dynamic elastic moduli, obtained from ultrasonic measurements, are suitable for describing the elastic response of the porous medium.For larger deformations (up to 10 -2 m/m), the static elastic parameter values, estimated via lab analyses, are better suited to describe the system behavior.The Eurocode [41,42] sets 1/300 as maximum admissible value for differential structural yielding.In the worst analyzed case (static elastic moduli), the maximum induced gradient at the end of primary production was equal to ~1/93000.It would be trivial to point out that estimated ground movements due to future storage activities in the studied field are not expected to compromise the safety of any structures or infrastructures, even if considering deltapressuring scenarios with an injection pressure up to 120% of the initial value.

CONCLUSION
The conversion of a marginal field into a gas storage facility implies the development of a reliable fluid flow/geomechanical numerical model at the regional scale for defining safe operational conditions.This includes the evaluation of potentially induced subsidence phenomena.The construction of an integrated model for geomechanical simulation purposes requires the integration of geological, fluid flow and mechanical data both at the regional and reservoir scales.But the costs related to a detailed characteriza-tion of such a large rock volume can compromise the economic feasibility of the storage project.
Through a case history of a marginal field located in Southern Italy, the present work showed how a reliable subsidence analysis could be performed despite the lack of direct data at the regional scale.The results obtained were validated against the measurements of ground surface movements available from the technical literature for the area under investigation.Not only could the calculated magnitude of subsidence be confirmed but also the appropriateness of the selected deformation properties of the rocks.
Despite the fact that the expected overall geomechanical storage response could be reliably assessed, the developed model might be subsequently updated and upgraded based on the results of possible future acquisition campaigns.The availability of the time evolution of the ground movements in the investigated area via satellite acquisition [43][44][45][46][47] together with the reservoir production history might allow the numerical model calibration via a back analysis approach, where the rock deformation parameters might be fine-tuned.
Furthermore, it should be noted that the safety of any storage system must also be evaluated in terms of cap rock integrity and possible faults (re)activation, especially in the case of delta pressuring management.Thus a direct mechanical characterization of the reservoir and cap rock formations via lab analyses and well logs is still necessary, yet data acquisition can be targeted at the most critical portion of the system rather than extended to the entire volume to be modeled.

Fig. (
Fig. (3).Schematic geological section across the North Bradanic Trough.Some oil and gas fields located in the area are projected[14].In particular, the Grottole-Ferrandina reservoir is shown to the north-east: the lenticular (lobate) shape of the sandy intercalation and the drapping upon the pre-pliocene substratum high can be observed.

Fig. (
Fig. (4).Regional and reservoir models created for this study.On the left, the regional model is presented and on the right, the reservoir model along with their average dimensions to emphasize the different size of the two grids.The wells located inside or very close to the reservoir area are also shown.

Fig. ( 7
Fig. (7).Three characteristic seismic sections taken from the studied area where the geometry of the Allocton (yellow color) and the Carbonate platform (Blue color) is shown.The front of the Allocton wedge appears more advanced to the north and south of the studied area and retreated in the central part.The inset shows the locations of the three main sections and a graphical sketch of the Allocton front.

Fig. ( 8 )
Fig. (8).Two seismic sections lying outside the model boundary and used to control the geometry of the main features in a broader scale.The inset shows the locations of the two seismic sections and a graphical sketch of the Allocton front.

Fig. ( 9 ).
Fig.(9).Characteristic geometry of the Allocton wedge together with the locations of the regional wells.The retreat of the wedge close to well Pizzo Corvo 1 can be observed.

Fig. ( 10
Fig.(10).Characteristic geometry of the faulted Carbonate platform and the creation of a Horst-Graben geometry, together with the locations of the regional wells.
Fig. (4), extends for 10 km in the NW-SE direction and almost 2 km in the NE-SW direction.

Fig. ( 11 ).
Fig.(11).Main faults used to define the geometry of the carbonatic platform (Top Cretaceous).The fault affecting the reservoir area is also shown.

Fig. ( 12 ).
Fig. (12).Regional model representation where the main stratigraphic features are shown in different colors.

Fig. ( 15
Fig. (15).Rock compressibility as a function of vertical stress and Elastic Modulus as a function of depth.

Table 1 .
Characterization for each geomechanical class: deformation, initialization and strength parameters.

Fig. ( 21 ).
Fig. (21).Average pressure distribution at the top of the reservoir (a); subsidence cone on ground surface in static (b) and dynamic (c) cases.

Fig. ( 22
Fig.(22).Comparison between subsidence cone on ground surface in case of static and dynamic assumption along AB section.

Fig. ( 23 ).
Fig.(23).Evolution in time of mean reservoir pressure and vertical displacement of ground surface during primary production.

Fig. ( 24
Fig.(24).Evolution in time of mean reservoir pressure and vertical displacement of ground surface for case 1.

Fig. ( 25 ).
Fig. (25).Evolution in time of mean reservoir pressure and vertical displacement of ground surface for case 2.