Elasto-viscoplastic modeling of subsidence above gas fields in the Adriatic Sea

From the analysis of GPS monitoring data collected above gas fields in the Adriatic Sea, in a few cases subsidence responses have been observed not to directly correlate with the production trend. Such behavior, already described in the literature, may be due to several physical phenomena, ranging from simple delayed aquifer depletion to a much more complex time-dependent mechanical response of subsurface geomaterials to fluid withdrawal. In order to accurately reproduce it and therefore to be able to provide reliable forecasts, in the last years Eni has enriched its 3D finite element geomechanical modeling workflow by adopting an advanced constitutive model (Vermeer and Neher, 1999), which also considers the viscous component of the deformation. While the numerical implementation of such methodology has already been validated at laboratory scale and tested on synthetic hydrocarbon fields, the work herein presents its first application to a real gas field in the Adriatic Sea where the phenomenon has been observed. The results show that the model is capable to reproduce very accurately both GPS data and other available measurements. It is worth remarking that initial runs, characterized by the use of model parameter values directly obtained from the interpretation of mechanical laboratory tests, already provided very good results and only minor tuning operations have been required to perfect the model outcomes. Ongoing R&D projects are focused on a regional scale characterization of the Adriatic Sea basin in the framework of the Vermeer and Neher model approach.


Introduction
In the last two decades Eni S.p.A. has been developing a robust modeling methodology for production-induced subsidence (Capasso and Mantica, 2006). It is based on 1way hydro-mechanical coupling (Gambolati et al., 2005) and elasto-plastic modified Cam-Clay model (MCCM, Roscoe and Burland, 1968). Though taking into account inelastic strains is instrumental for describing compaction in clastic reservoirs (Pijnenburg et al., 2019) and the MCCM keeps providing accurate reproduction of monitoring data gathered from almost all the gas fields in the Adriatic Sea (e.g. Gemelli et al., 2015), a further constitutive modeling effort has been recently required for a few of them showing a certain delay between production trend and GPS data -a phenomenon broadly described in the literature (e.g. Hettema et al., 2002).
The modeling approach has been enhanced by adopting the elasto-viscoplastic model proposed by Vermeer and Ne-her (1999, VNM), capable to describe the viscous response of reservoir sands and derived from the extended overstress theory (Olszak and Perzyna, 1970;Yin et al., 2010).
Then, having the implementation already been validated at laboratory scale and tested at reservoir scale on synthetic hydrocarbon fields (Volonté et al., 2017;Musso et al., 2020), this paper presents a first application of the enhanced approach to the production-induced subsidence analysis of a real gas field in the Adriatic Sea, the GPS data of which exhibit a delay of about 1.5 year (Fig. 1).
Herein, because of confidentiality issues, field data have been anonymized and analysis results normalized.
Published by Copernicus Publications on behalf of the International Association of Hydrological Sciences.

Field and production
The off-shore gas field studied herein is located in the Adriatic Sea, at about 60 km from the Italian coastline, where the average water depth is around 60 m.
The sandy reservoir layers lie from 900 to 1800 m s.s.l. and are produced by 28 wells, connected to platforms A and B.
According to the Intersect ® fluid-dynamic model, the gas volume originally in place is approximately 30 GSm 3 and the recovery factor expected at forecast end is about 50 %.

Geomechanical modeling
The subsidence analysis has been performed by means of a 3D FE model built with the commercial code Abaqus ® .
Input data about geometry, geology and petrophysics, of both reservoir layers and hydraulically connected aquifers, have been provided by the fluid-dynamic model, the same for pore pressure distribution and time evolution.
The domain has been discretized with about 5.5×10 5 finite elements. The model has around 2 × 10 6 degrees of freedom.
For the 6 input parameters of the VNM (Volonté et al., 2017;Musso et al., 2020), a preliminary estimate has been   directly obtained from an experimental campaign of tailored laboratory oedometric compression tests, characterized by a creep phase and carried out on samples of bottom hole cores from the same field. An initial geomechanical simulation has been performed with these preliminary values of the param-eters. Then, in order to accurately reproduce the GPS data recorded at platform B, a calibration operation has been performed (Fig. 2), obtaining the final values of the parameters (Table 1). This step has required only a very minor tuning of the creep index µ * (less than 4 % variation), that is compatible with the uncertainty associated to interpretation of laboratory tests. All other parameters have been left unchanged.

Results
While GPS data from platform B have been used to calibrate the VNM parameters, other available data from the same platform are useful for evaluating the capability of the geomechanical model to accurately simulate the hydromechanical response of the field to gas withdrawal. To this purpose, Figs. 3 and 4 present comparisons in terms of reservoir compaction and iso-subsidence lines, respectively, for platform B.
In particular, Fig. 3 shows the cumulative compaction observed at reservoir depth along a well of platform B, where almost yearly a special logging tool is run for monitoring the distance between the radioactive markers. The comparison with the corresponding model estimates is very satisfactory: in fact, except for the measure acquired at year 11.7, which is out of trend, all the others are reproduced within the error bar or slightly overestimated.    Figure 4 shows the good agreement in terms of subsidence developed in the 7 year interval time elapsed between 2 bathymetric surveys performed around platform B, at year 3 and 10, respectively, after the production start. Figure 5 shows the expected subsidence evolution at platform A. Even if calibrated on data from platform B, the model reproduces properly the subsidence rate recorded at the GPS station installed on platform A. Figure 6 shows time evolution of the 2 cm iso-subsidence line, plus values of minimum distance from the coastline and maximum extent at simulation-end, which is usually set 30 years after the production-end.
Finally, Fig. 7 shows subsidence estimates provided by both VNM and MCCM, this latter with parameter values from lab tests and 1.2 overconsolidation ratio.

Concluding remarks
The subsidence analysis presented herein is the first application to a real gas field of the Eni's enhanced 3D finite element geomechanical workflow. The results show a very accurate reproduction of the monitoring data and the significative improvement obtained by adopting the VNM. Data availability. Data sets have been used but can not be published for a matter of confidentiality.
Author contributions. FG performed and supervised the whole work. AC provided support on numerical modelling and implementation of constitutive laws, GV on rock mechanics laboratory testing, GV and SM on constitutive modelling, SM also on subsidence monitoring data assessment and MDS on fluid-dynamic modelling.
Competing interests. The authors declare that they have no conflict of interest.

Special issue statement.
This article is part of the special issue "TISOLS: the Tenth International Symposium On Land Subsidence -living with subsidence". It is a result of the Tenth International Symposium on Land Subsidence, Delft, the Netherlands, 17-21 May 2021.