Multidisciplinary Analysis of Ground Movements: An Underground Gas Storage Case Study

The paper presents a multi-physics investigation of the ground movements related to the cyclical and seasonal injection and withdrawal of natural gas in/from a depleted reservoir located in the Po Plain area, Italy. Interferometric Synthetic Aperture Radar (InSAr) data (from 2003) and Global Navigation Satellite System (GNSS) data (from 2008) provided a full and coherent panorama of almost two decades of ground movement in the monitored area (more extended than the field boundary). The analysis of the acquired millimetric-scale movements together with the detailed geological analysis, both at reservoir and at regional scale, represents the focal point for understanding the investigated phenomena. Based on this information, a fully integrated and multidisciplinary geological, fluid-flow and geomechanical numerical modeling approach was developed to reproduce the main geometrical and structural features of the involved formations together with the poromechanics processes induced by the storage operations. The main achievement of the adopted methodology is a deep knowledge of the system and the involved processes, which is mandatory for the safety of the urbanized areas and the effective management of the underground resources.


Introduction
The ground surface movements due to anthropogenic activities, such as underground fluid production/injection, have been deeply investigated by the scientific literature [1][2][3][4], among the others, particularly in potentially critical environments such as high-urbanized areas and costal zones. In the Italian panorama, a large number of hydrocarbon fields were discovered and produced since the early 1950s onward in the Po Plain area (e.g., [5]); due to its high degree of urbanization, the area has been the focus of numerous studies [6][7][8][9][10][11], among others. The monitoring, investigation and forecasting approaches adopted for the subsidence analyses have shown a constant and continuous improvement of the technologies and methodologies to fulfill the higher and higher standards of the system safety and the social acceptance. The ground movement monitoring has improved from Global Navigation Satellite System (GNSS) to Interferometric Synthetic Aperture Radar (InSAR) technologies and the analysis of the involved multiple physical processes has improved from analytical to 3D numerical approaches.

Materials and Methods
The aim of the surface monitoring plan was the determination of the planar and vertical components of movement and of the uplift/subsidence area ascribed to UGS operations via two complementary technologies: the satellite SAR interferometry and the GNSS techniques. SAR interferometry investigated an area larger than the field area providing relative measurements of movements whereas GNSS monitoring provided punctual and absolute measurements.
The SAR interferometry is widely used to detect and monitor slow terrain movements with millimeter accuracy along the satellite line of sight (LOS), covering long-term period with frequent updates and large areas with tens of thousands of measurements per square kilometer. The SAR image phase is dependent both on the distance between the radar antenna and the ground target, and on the radiometric characteristics of the transmission medium and of the target. From a stack of interferometric phases (i.e., phase difference between two SAR acquisitions), displacement information of the observed scene is extracted [14] for a set of sparse points (Persistent Scatterers-PSs). The PSs are physical targets characterized by temporally stable backscattering properties and they are commonly found in poorly vegetated and urbanized areas [15][16][17]. The present research was developed according to the Persistent Scatterer Pairs (PSP-IFSAR) approach; it is a multi-interferogram InSAR technique for deriving dense and reliable information in a set of sparse points (PSs). The adopted PSP-IFSAR method, which integrates both ascending and descending acquisition geometries, measures the same ground deformation along two-satellite lines of sight, and, consequently, the actual vertical and East-West planar components of the movements can be determined. The North-South component cannot be derived because the SAR sensors are insensitive to the North-South displacements.
The GNSS (formerly "GPS") is a space geodetic technique able to determine 3D coordinates of a monitoring permanent site with sub-centimeter accuracy and its over-time kinematic evolution (i.e., its Remote Sens. 2020, 12, 3487 3 of 19 velocity of movement along the planar and vertical components). This geolocation technique is based on more than one constellation of satellites (NAVSTAR for the GPS), numerous ground permanent stations receiving satellite signals and an International GNSS Service (IGS), which ensures high-quality GNSS data products (e.g., GNSS satellite ephemerides, Earth rotation parameters, velocities, etc.). For the present research, the GNSS data were processed using BERNESE software (developed by the Astronomical Institute, University of Berne, Switzerland) [18] through two different approaches ( Table 1). The Network Approach adopts about nine reference stations of known coordinates, in addition to the estimating network, in order to set the terrestrial reference frame; the process uses double-differences of the GNSS observables (then a network of stations acquired simultaneously is required). The Precise Point Positioning Approach uses one-way GNSS observables; it does not require a network of stations acquiring simultaneously, but it strictly requires: GNSS data collected by the station to be estimated, precise GNSS orbits, related Earth rotation parameters, corrections for the clocks bias and drift, and specific mapping function for the atmospheric delay modeling (neutral part of the atmosphere). The acquired and interpreted ground monitoring data were used during the back analysis of the geomechanical modeling, which represents the final step of a complete geological, fluid-flow and geomechanical numerical modeling workflow ( Figure 1). The 3D numerical modeling approach was based on a coherent and full dataset, including: 3D seismic survey, primary production and storage data, well-logs, petrophysical and fluid properties; in situ Image log and MDT stress tests, deformation and strength parameters from lab tests, logs and technical literature.
based on more than one constellation of satellites (NAVSTAR for the GPS), numerous ground permanent stations receiving satellite signals and an International GNSS Service (IGS), which ensures high-quality GNSS data products (e.g., GNSS satellite ephemerides, Earth rotation parameters, velocities, etc.). For the present research, the GNSS data were processed using BERNESE software (developed by the Astronomical Institute, University of Berne, Switzerland) [18] through two different approaches (Table 1). The Network Approach adopts about nine reference stations of known coordinates, in addition to the estimating network, in order to set the terrestrial reference frame; the process uses double-differences of the GNSS observables (then a network of stations acquired simultaneously is required). The Precise Point Positioning Approach uses one-way GNSS observables; it does not require a network of stations acquiring simultaneously, but it strictly requires: GNSS data collected by the station to be estimated, precise GNSS orbits, related Earth rotation parameters, corrections for the clocks bias and drift, and specific mapping function for the atmospheric delay modeling (neutral part of the atmosphere). The acquired and interpreted ground monitoring data were used during the back analysis of the geomechanical modeling, which represents the final step of a complete geological, fluid-flow and geomechanical numerical modeling workflow ( Figure 1). The 3D numerical modeling approach was based on a coherent and full dataset, including: 3D seismic survey, primary production and storage data, well-logs, petrophysical and fluid properties; in situ Image log and MDT stress tests, deformation and strength parameters from lab tests, logs and technical literature. Starting from the detailed reservoir static model, describing all the main stratigraphic, structural, lithological and petrophysical characteristics of the reservoir, an extended (regional-scale) geological Starting from the detailed reservoir static model, describing all the main stratigraphic, structural, lithological and petrophysical characteristics of the reservoir, an extended (regional-scale) geological model was set up (around 450 10 3 cells). The latter effectively reproduces the main stratigraphic and structural features of the investigation domain for the ground movement analysis, which includes the reservoir itself and its surrounding formations, up to the surface (Figure 2). A multiphase flow numerical model (FDM) was then set up based on the geological model and on all the available petrophysical and fluid properties together with well completion data. The reservoir consists in sand layers characterized by 18-23% porosity and 100-300 mD (milliDarcy) permeability ranging values; the permeability anisotropy at the reservoir scale reflects the occurrence of pelitic interlayers within the reservoir. The model was then initialized and calibrated to reproduce the pressure distribution in the reservoir during more than 60 years of production and storage history, considering data from 80+ wells. The pressure distribution in the reservoir is homogeneous as the reservoir has no internal compartmentalization; the cyclical pressure variations follow the cyclical gas storage volume variations that were compared with the ground movement measures during the InSAR and GNSS analyses. The dynamic model is periodically updated at the end of each injection/withdrawal phase as new pressure data are acquired from the monitoring wells.
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 19 model was set up (around 450 10 3 cells). The latter effectively reproduces the main stratigraphic and structural features of the investigation domain for the ground movement analysis, which includes the reservoir itself and its surrounding formations, up to the surface ( Figure 2). A multiphase flow numerical model (FDM) was then set up based on the geological model and on all the available petrophysical and fluid properties together with well completion data. The reservoir consists in sand layers characterized by 18-23% porosity and 100-300 mD (milliDarcy) permeability ranging values; the permeability anisotropy at the reservoir scale reflects the occurrence of pelitic interlayers within the reservoir. The model was then initialized and calibrated to reproduce the pressure distribution in the reservoir during more than 60 years of production and storage history, considering data from 80+ wells. The pressure distribution in the reservoir is homogeneous as the reservoir has no internal compartmentalization; the cyclical pressure variations follow the cyclical gas storage volume variations that were compared with the ground movement measures during the InSAR and GNSS analyses. The dynamic model is periodically updated at the end of each injection/withdrawal phase as new pressure data are acquired from the monitoring wells. Finally, the ground displacement analysis was addressed via a stress-strain finite element (FEM) model, set up based on the extended geological study. The model was characterized by mechanical properties derived from lab and log data ( Table 2). During the initialization phase, the initial pore pressure distribution and the original stress state were defined for all the investigated volumes according to the available data (i.e., initial static pressure, well logs and in situ tests). Dirichlet boundary conditions were imposed (i.e., zero displacement on the lateral and bottom boundaries, zero normal stress on the surface). Rock mechanics engineering analyses were developed according to the elastic-purely plastic constitutive law and adopting the Mohr-Coulomb failure criteria. As a matter of fact, different studies focused on ground subsidence induced by UGS for analogous case studies have shown the reliability of the elastic constitutive model in simulating the measured surface movement data [19][20][21]. It should be taken into consideration that the UGS-related pressure variations affect the formation both cyclically (following the withdrawal/injection phases) and in a very small time frame (5/6 moths): it clearly differs from a standard monotonic pressure decrease over some decades that typically characterizes the primary production, and the formations react accordingly. On the base of the preliminary sensitivity analysis results, the model of the case under study reacted within the elastic domain even at the end of primary production when the system experienced the maximum pressure variation. An adequate time stepping was set up accordingly for the UGS system analysis: the mechanical simulations were performed at the end of the primary Finally, the ground displacement analysis was addressed via a stress-strain finite element (FEM) model, set up based on the extended geological study. The model was characterized by mechanical properties derived from lab and log data ( Table 2). During the initialization phase, the initial pore pressure distribution and the original stress state were defined for all the investigated volumes according to the available data (i.e., initial static pressure, well logs and in situ tests). Dirichlet boundary conditions were imposed (i.e., zero displacement on the lateral and bottom boundaries, zero normal stress on the surface). Rock mechanics engineering analyses were developed according to the elastic-purely plastic constitutive law and adopting the Mohr-Coulomb failure criteria. As a matter of fact, different studies focused on ground subsidence induced by UGS for analogous case studies have shown the reliability of the elastic constitutive model in simulating the measured surface movement data [19][20][21]. It should be taken into consideration that the UGS-related pressure variations affect the formation both cyclically (following the withdrawal/injection phases) and in a very small time frame (5/6 moths): it clearly differs from a standard monotonic pressure decrease over some decades that typically characterizes the primary production, and the formations react accordingly. On the base Remote Sens. 2020, 12, 3487 5 of 19 of the preliminary sensitivity analysis results, the model of the case under study reacted within the elastic domain even at the end of primary production when the system experienced the maximum pressure variation. An adequate time stepping was set up accordingly for the UGS system analysis: the mechanical simulations were performed at the end of the primary production, and at the end of each production/injection period for all the InSAR monitoring time frames. Furthermore, a negligible variation of the petrophysical parameters due to pressure change was expected, in agreement with other UGS in analogous clastic formations studies [7,11,21]. Consequently, an uncoupled fluid-flow and stress-strain approach was adopted: the time and space pressure evolution represents the forcing function applied to the geomechanical model, whereas the petrophysical parameter values adopted in the fluid-flow model remain unchanged. The historical pressure evolution and the vertical ground displacement data were finally integrated in the numerical model via a back analysis approach: the adopted pseudo-elastic parameters of the reservoir and of the cap rock were fine-tuned to achieve a suitable match between simulated and measured movements. Previous experiences of the authors and the technical literature (please consider the abovementioned references) highlighted that the UGS-related ground movements for a clastic reservoir at about 1000-1500 m depth, placed in the specific geological contest under analysis (paragraphs 3.1 and 3.2), are mainly influenced by the stress-strain behavior of the gas bearing formations and, secondly, by the stress-strain behavior of the cap rock.

Regional Geological Setting
A deep knowledge of the geological characteristics of the examined formations is fundamental for the construction of accurate and representative 3D geological models and for the interpretation of ground monitoring data. Regarding the case study, a detailed geological analysis was developed at both reservoir and regional scales and it provided insight of the internal stratigraphic-structural architecture and of the heterogeneity of the investigated volume.
The studied UGS field is placed in the eastern Po Plain (North of Italy) ( Figure 3). The Po Plain-Northern Adriatic region was extensively studied through exploration well data and 2D-3D seismic surveys acquired over the last 40 years (e.g., [22][23][24][25]). During Jurassic-Triassic time the Remote Sens. 2020, 12, 3487 6 of 19 region experienced east-west directed extensional phases which led to the formation of fault-bounded basins (e.g., [23,26]). In response to the subsequent Eurasia-Adria plate convergence, the rift-related framework was overprinted by Cenozoic compressional structures that are currently buried below the Po Plain sedimentary infill. The present-day architecture of the Po Plain is mainly the result of the Late Messinian-Pleistocene tectonic evolution; the latter led to the development of the Po Plain-Adriatic basin as a shared, complex foredeep separating the converging Southern Alpine (to the north) and Northern Apennine (to the south) chains, which locally collided in the Po Plain subsurface (e.g., [27,28]).
The buried Northern Apennines is an arc-shaped fold-and-thrust belt subdivided into three adjoining arcs named Monferrato, Emilia and Ferrara-Romagna arcs (from the west to the east, respectively) [29]. Its subsurface architecture consists of synclines and ramp anticlines associated with blind thrusts formed in response to a series of successive tectonic pulses that have led to the progressive North-Northeastward migration of the outermost thrust fronts ( Figure 3). During the latest Miocene-Pliocene event, the eastern Po Plain has undergone a strong tectonic activity; the Late Pliocene-Pleistocene tectonic events led to the complete development of the outer fold and thrust system of Ferrara and to the uplift and tilting of the Bologna area [24]. The area of the UGS field is characterized by NE-verging, arcuate blind thrusts which involved the Mesozoic units, their pre-Mesozoic basement and the younger clastic infill of the Plio-Pleistocene foredeep ( Figure 3). In the field area there is no evidence of the presence of thrusts dislocating the Pleistocene sediments or with reactivations during the Late Pliocene-Pleistocene times: none of the thrusts that created significant dislocation of the pre-Pliocene sequences propagated across the base of the Early Pleistocene deposits (e.g., [25,30,31] The buried Northern Apennines is an arc-shaped fold-and-thrust belt subdivided into three adjoining arcs named Monferrato, Emilia and Ferrara-Romagna arcs (from the west to the east, respectively) [29]. Its subsurface architecture consists of synclines and ramp anticlines associated with blind thrusts formed in response to a series of successive tectonic pulses that have led to the progressive North-Northeastward migration of the outermost thrust fronts (Figure 3). During the latest Miocene-Pliocene event, the eastern Po Plain has undergone a strong tectonic activity; the Late Pliocene-Pleistocene tectonic events led to the complete development of the outer fold and thrust system of Ferrara and to the uplift and tilting of the Bologna area [24]. The area of the UGS field is characterized by NE-verging, arcuate blind thrusts which involved the Mesozoic units, their pre-Mesozoic basement and the younger clastic infill of the Plio-Pleistocene foredeep ( Figure 3). In the field area there is no evidence of the presence of thrusts dislocating the Pleistocene sediments or with reactivations during the Late Pliocene-Pleistocene times: none of the thrusts that created significant dislocation of the pre-Pliocene sequences propagated across the base of the Early Pleistocene deposits (e.g., [25,30,31]).  The Pliocene foredeep infill in the UGS field area includes the thick syntectonic, sand turbidites of the Porto Corsini Fm. (Early-Middle Pliocene) and Porto Garibaldi Fm. (Middle-Late Pliocene); meanwhile, the clay sequence of the Santerno Fm. (Pliocene-Early Pleistocene) was deposited above the submarine paleo-highs and part of the foreland and foreland ramp. Finally, the late Pliocene-Pleistocene sedimentary infill recorded the gradual transition from deep-marine to continental environment (i.e., Asti Group), through the development of the basin-scale Po Plain prograding complex as a result of the east-ward progradation of the Po River delta. The overall Plio-Quaternary clastic sequence reaches a thickness of about 7-8 km in the deepest depocentral areas of the basin [22,24,32].

UGS Gas Field Description
The investigated gas field is located approximately 20 km NE of the Bologna city ( Figure 3) at an average depth of 1300 m ssl. The gas primary production started in 1956 and, after two decades, the reservoir was converted into a gas storage system. Currently the system is managed at a maximum injection pressure close to the initial pressure of the reservoir with injection and production rates in the order of 3 to 25 MSm 3 /d and 5 to 40 MSm 3 /d, respectively; 50+ wells were drilled on the trap culmination for the storage operations and pressure monitoring.
The reservoir consists of a mixed (stratigraphic-structural) trap that combines an asymmetric anticline controlled by an NE-verging, blind thrust and the deposition of sandy levels onlapping on the structural high with local pinching-out features. The pinch-out closure over the structural high defines a non-deposition area in the SE part of the field. Minor, extensional faults, sub-parallel to the trap axis, affect the axial area. Neither the stratigraphic nor the structural features produced reservoir compartmentalization.
The reservoir sequence belongs to the Porto Garibaldi Fm. (Middle-Late Pliocene) and consists of an alternance of sand layers, meter-to-decameter thick, and clay layers, 1 to 10 m of average thicknesses. In 2011, a 3D seismic survey was acquired and processed covering the entire UGS area (Figure 3). The approximate extension of the acquisition was 11.5 km × 7 km and the investigation depth almost 12 km. The results were integrated and interpreted together with data from the existing wells. The cross-section derived from the 3D seismic volume shown in Figure 4 highlights the geometry of the trap and of the main thrust bounding the reservoir to the North. The three main marker surfaces here represented correspond to the Top reservoir (Late Pliocene), Top Santerno (Calabrian, i.e., Early Pleistocene) and Top Prograding (Middle Pleistocene). The seismic data show that the thrust that bounds the trap ends within the Santerno Fm. without reaching the Top Santerno horizon and thus not affecting the Quaternary sequence as well as none of the faults inside the seismic volume. According to the literature, it confirms the absence of evidence of recent reactivation for the Neogenic faults present in the investigated area.
here represented correspond to the Top reservoir (Late Pliocene), Top Santerno (Calabrian, i.e., Early Pleistocene) and Top Prograding (Middle Pleistocene). The seismic data show that the thrust that bounds the trap ends within the Santerno Fm. without reaching the Top Santerno horizon and thus not affecting the Quaternary sequence as well as none of the faults inside the seismic volume. According to the literature, it confirms the absence of evidence of recent reactivation for the Neogenic faults present in the investigated area.

SAR interferometry results
The huge archive of InSAR data (C-band) ( Table 3) was analyzed via the PSP-IFSAR technique. Radarsat-1/2 data have provided a very long series of ground displacement measurements starting from 2003 with a temporal resolution of 24 days in both observation geometries. Starting from 2014, Sentinel-1 data have provided more detailed information on the last 5 years with a new acquisition every 6 days in both observation geometries. The satellite data cover an area larger than the field area, as shown in Figure 3. Table 3. Main characteristics of the datasets. RSAT1: Radarsat-1; RSAT2: Radarsat-2; S1: Sentinel-1.

Satellite
Observation Based on ascending and descending results, for each sensor, the vertical and East-West planar movements were derived. The PSP measurements of the ground movement at each acquisition date allowed the evaluation of its mean velocity (as linear trend between the first and the last acquisition of each interferometric stack) and its mean amplitude (in case of sinusoidal behavior).
The mean velocity data show the overall stability of the field area with a slight subsidence trend; local zones with slight uplift or subsidence are evident in the vertical velocity map, with a notable subsidence phenomenon in the South-West area towards the city of Bologna ( Figure 5). The mean vertical velocity in the area above the field ranges between movements were derived. The PSP measurements of the ground movement at each acquisition date allowed the evaluation of its mean velocity (as linear trend between the first and the last acquisition of each interferometric stack) and its mean amplitude (in case of sinusoidal behavior).
The mean velocity data show the overall stability of the field area with a slight subsidence trend; local zones with slight uplift or subsidence are evident in the vertical velocity map, with a notable subsidence phenomenon in the South-West area towards the city of Bologna ( Figure 5). The mean vertical velocity in the area above the field ranges between −1.0 and −2.0 mm/year. The submillimetric difference between the Radarsat and Sentinel results can be attributed to the different time spans they investigated (2003-2019 vs. 2015-2019, respectively).  Figure 6 shows the Radarsat maps of the mean vertical and horizontal amplitude of the ground displacement. The well-defined area of maximum vertical amplitude is located above the center of the UGS field ( Figure 6a) and it coincides with the area of minimum horizontal amplitude ( Figure  6b). Instead, the areas of the maximum horizontal amplitude identify two isolated peaks in the East-West direction around the center of the reservoir; zero amplitude values are detected along the North-South direction because of the limit of the InSAR technique.   Six representative points were selected in the InSAR-monitored area within and outside the field boundary ( Figure 6a): their historical displacement along the vertical and the East-West directions were compared with the curve of the storage gas cumulative volumes of the whole field. The temporal evolution of the two ground displacement components for the P1, P2, P3 internal points shows a clear sinusoidal signal with a strong correlation with the trend of the gas cumulative volume curve in terms of amplitude and periodicity (Figures 7 and 8). The vertical displacement always shows the maximum and minimum peaks at the end of each injection and withdrawal period, respectively, with a delay of about 30 days, which could be attributed primarily to fluid-flow phenomena (i.e., time and space propagation of the pressure sink in the reservoir and surrounding aquifer) (Section 2. Materials and Methods). Six representative points were selected in the InSAR-monitored area within and outside the field boundary ( Figure 6a): their historical displacement along the vertical and the East-West directions were compared with the curve of the storage gas cumulative volumes of the whole field. The temporal evolution of the two ground displacement components for the P1, P2, P3 internal points shows a clear sinusoidal signal with a strong correlation with the trend of the gas cumulative volume curve in terms of amplitude and periodicity (Figures 7 and 8). The vertical displacement always shows the maximum and minimum peaks at the end of each injection and withdrawal period, respectively, with a delay of about 30 days, which could be attributed primarily to fluid-flow phenomena (i.e., time and space propagation of the pressure sink in the reservoir and surrounding aquifer) (Section 2. Materials and Methods). boundary (Figure 6a): their historical displacement along the vertical and the East-West directions were compared with the curve of the storage gas cumulative volumes of the whole field. The temporal evolution of the two ground displacement components for the P1, P2, P3 internal points shows a clear sinusoidal signal with a strong correlation with the trend of the gas cumulative volume curve in terms of amplitude and periodicity (Figures 7 and 8). The vertical displacement always shows the maximum and minimum peaks at the end of each injection and withdrawal period, respectively, with a delay of about 30 days, which could be attributed primarily to fluid-flow phenomena (i.e., time and space propagation of the pressure sink in the reservoir and surrounding aquifer) (Section 2. Materials and Methods).  The horizontal displacement shows alike clear sinusoidal signals: during the withdrawal periods P1-P2 points show positive trends (i.e., eastward direction) and P3 point shows negative trend (i.e., westward direction); vice versa, during the injection periods, P1-P2 points show negative trends and P3 point shows a positive one. The divergent responses can be attributed to the different locations of the points in respect to the area of maximum vertical displacement (i.e., reservoir center). The time series of the points outside the field area (P4, P5, P6) show displacement trends, oscillation characteristics and periodicity, both in vertical and horizontal directions, unrelated with the storage activity (Figures 9 and 10). The horizontal displacement shows alike clear sinusoidal signals: during the withdrawal periods P1-P2 points show positive trends (i.e., eastward direction) and P3 point shows negative trend (i.e., westward direction); vice versa, during the injection periods, P1-P2 points show negative trends and P3 point shows a positive one. The divergent responses can be attributed to the different locations of the points in respect to the area of maximum vertical displacement (i.e., reservoir center). The time series of the points outside the field area (P4, P5, P6) show displacement trends, oscillation characteristics and periodicity, both in vertical and horizontal directions, unrelated with the storage activity (Figures 9 and 10). westward direction); vice versa, during the injection periods, P1-P2 points show negative trends and P3 point shows a positive one. The divergent responses can be attributed to the different locations of the points in respect to the area of maximum vertical displacement (i.e., reservoir center). The time series of the points outside the field area (P4, P5, P6) show displacement trends, oscillation characteristics and periodicity, both in vertical and horizontal directions, unrelated with the storage activity (Figures 9 and 10).  Note that the Radarsat series are more sensitive to the long-term displacement trends whereas Sentinel series give a higher detail on the last monitoring years because of the different monitored westward direction); vice versa, during the injection periods, P1-P2 points show negative trends and P3 point shows a positive one. The divergent responses can be attributed to the different locations of the points in respect to the area of maximum vertical displacement (i.e., reservoir center). The time series of the points outside the field area (P4, P5, P6) show displacement trends, oscillation characteristics and periodicity, both in vertical and horizontal directions, unrelated with the storage activity (Figures 9 and 10).  Note that the Radarsat series are more sensitive to the long-term displacement trends whereas Sentinel series give a higher detail on the last monitoring years because of the different monitored Note that the Radarsat series are more sensitive to the long-term displacement trends whereas Sentinel series give a higher detail on the last monitoring years because of the different monitored time spans. The time series of the vertical displacement confirm a long-term trend of gentle surface downwarping in the field area.

GNSS results
The permanent GNSS station was installed at the end of 2008 in the area of the storage plant located at the south-eastern edge of the gas field (Figures 3 and 6). The station was equipped with an ASHTECH UZ-12 receiver since the installation date in 2013, when it was replaced by a LEICA GR10 receiver. The GNSS antenna is a TOPCON, TPSCR4 model, with TPSH radome, dual frequency GPS (L1/L2).
The daily RINEX data of the GNSS permanent station were analyzed via the TEQC software (http://facility.unavco.org/software/teqc), an international standard for the pre-processing phase of GPS data and for the evaluation of their quality. The data quality was stated by using the TEQC parameters MP1 (multipath on L1) and MP2 (multipath on L2). Based on the values of the IGS reference stations, the values of MP1 (<0.5 m) and MP2 (<0.75 m) indicate a good quality station ( Figure 11).
Remote Sens. 2020, 12, x FOR PEER REVIEW 11 of 19 time spans. The time series of the vertical displacement confirm a long-term trend of gentle surface downwarping in the field area.

GNSS results
The permanent GNSS station was installed at the end of 2008 in the area of the storage plant located at the south-eastern edge of the gas field (Figures 3 and 6). The station was equipped with an ASHTECH UZ-12 receiver since the installation date in 2013, when it was replaced by a LEICA GR10 receiver. The GNSS antenna is a TOPCON, TPSCR4 model, with TPSH radome, dual frequency GPS (L1/L2).
The daily RINEX data of the GNSS permanent station were analyzed via the TEQC software (http://facility.unavco.org/software/teqc), an international standard for the pre-processing phase of GPS data and for the evaluation of their quality. The data quality was stated by using the TEQC parameters MP1 (multipath on L1) and MP2 (multipath on L2). Based on the values of the IGS reference stations, the values of MP1 (<0.5 m) and MP2 (<0.75 m) indicate a good quality station ( Figure 11). The results of the two processing approaches, the Network Approach and the Precise Point Positioning, are reported in Figure 12 in terms of coordinates time series: their consistency attests the robustness of the results. The areal and vertical trends of displacement indicate a marked NE-ward movement and a gentle subsidence of the monitored site, respectively. The results of the two processing approaches, the Network Approach and the Precise Point Positioning, are reported in Figure 12 in terms of coordinates time series: their consistency attests the robustness of the results. The areal and vertical trends of displacement indicate a marked NE-ward movement and a gentle subsidence of the monitored site, respectively. In order to discriminate the effects of the UGS activities along the three components of displacement, the signals ascribed to other known phenomena were removed. The typical seasonal signals of the GNSS time series are essentially due to load variations caused by the redistribution of fluid masses on the Earth's surface, (e.g., [33]). It includes changes in air mass, which results in changes in surface air pressure, in ocean level due to terrestrial and solar tides, wind and pressure atmospheric, and, above all, in soil moisture (surface hydrological load). Once removed the longterm velocity from the GNSS time series, the residual time series were obtained. The latter were compared with the displacements due to changes in the surface hydrological load (the so-called "surface hydrological mass loading") by using the data of the EOST Loading Service (http://loading.u-strasbg.fr/) and the data estimated by MERRA2 [34]. The seasonal vertical movement recorded by the GNSS station is largely due to the variation in the hydrological load (the hydrological model well followed the analytical seasonal model) and does not show an appreciable displacement-UGS correlation. On the contrary, the hydrological model does not explain the seasonal shifts observed in the planar components, attributed to the UGS activity ( Figure 13). The two horizontal displacement components are seasonally dependent: during withdrawal periods the North-South component shows a positive trend (i.e., northward direction), and the East-West component shows a negative trend (i.e., westward direction), in agreement with the movement toward the peak of the maximum vertical displacement highlighted by InSAR data. On the contrary, during injection periods, the North-South component shows a negative trend, and the East-West component a positive trend. Furthermore, the main direction of planar movement is almost North-South (about N198°). In order to discriminate the effects of the UGS activities along the three components of displacement, the signals ascribed to other known phenomena were removed. The typical seasonal signals of the GNSS time series are essentially due to load variations caused by the redistribution of fluid masses on the Earth's surface, (e.g., [33]). It includes changes in air mass, which results in changes in surface air pressure, in ocean level due to terrestrial and solar tides, wind and pressure atmospheric, and, above all, in soil moisture (surface hydrological load). Once removed the long-term velocity from the GNSS time series, the residual time series were obtained. The latter were compared with the displacements due to changes in the surface hydrological load (the so-called "surface hydrological mass loading") by using the data of the EOST Loading Service (http://loading.u-strasbg.fr/) and the data estimated by MERRA2 [34]. The seasonal vertical movement recorded by the GNSS station is largely due to the variation in the hydrological load (the hydrological model well followed the analytical seasonal model) and does not show an appreciable displacement-UGS correlation. On the contrary, the hydrological model does not explain the seasonal shifts observed in the planar components, attributed to the UGS activity ( Figure 13). The two horizontal displacement components are seasonally dependent: during withdrawal periods the North-South component shows a positive trend (i.e., northward direction), and the East-West component shows a negative trend (i.e., westward direction), in agreement with the movement toward the peak of the maximum vertical displacement highlighted by InSAR data. On the contrary, during injection periods, the North-South component shows a negative trend, and the East-West component a positive trend. Furthermore, the main direction of planar movement is almost North-South (about N198 • ).

Geomechanical Simulation Results
In the back-analysis phase of the mechanical model, the calibration of the pseudo-elastic parameters of the reservoir and the cap rock allowed a suitable match between the simulated vertical ground movements and the measured values discussed in the previous paragraph. Radarsat data were considered because of their longer period of acquisition. Figure 14 shows the comparison between simulated and measured values, in terms of relative variations within a single withdrawal/injection cycle, for the three points placed at surfaces within the boundary of the field (location in Figure 6) during a monitored time frame. For P2 point, a descending trend recognized by InSAR analysis, and not ascribable to UGS operations, was superimposed to the simulation results. The analysis of the InSAR data allowed the identification also of the areal extension of the ground surface cyclically affected by the UGS activities: Figure 15 shows an example of comparison between measured and simulated delta displacement maps at surface level, for a historical withdrawal period during UGS activity.

Geomechanical Simulation Results
In the back-analysis phase of the mechanical model, the calibration of the pseudo-elastic parameters of the reservoir and the cap rock allowed a suitable match between the simulated vertical ground movements and the measured values discussed in the previous paragraph. Radarsat data were considered because of their longer period of acquisition. Figure 14 shows the comparison between simulated and measured values, in terms of relative variations within a single withdrawal/injection cycle, for the three points placed at surfaces within the boundary of the field (location in Figure 6) during a monitored time frame. For P2 point, a descending trend recognized by InSAR analysis, and not ascribable to UGS operations, was superimposed to the simulation results. The analysis of the InSAR data allowed the identification also of the areal extension of the ground surface cyclically affected by the UGS activities: Figure 15 shows an example of comparison between measured and simulated delta displacement maps at surface level, for a historical withdrawal period during UGS activity.

Ground Monitoring Data and Comparison with UGS
Since 1950s, the eastern Po Plain area, where the studied field is located, has been experiencing a widespread and strong subsidence mainly due to anthropic causes, with a minor contribution (of a few mm/year) from ongoing natural processes (e.g., tectonics, sediment consolidation) [8,35,36]. In particular, the subsidence evolution in this region mostly depended on the groundwater pumping [8]. This phenomenon was particularly relevant in the Bologna area, where a subsidence rate in the order of several cm/year was historically measured; subsequently, the impressive subsidence reduction monitored during recent decades was a direct result of the reduction in groundwater pumping activities [36].
The huge spatial and temporal extension of the InSAR monitoring allowed the investigation, not only of the UGS field area but also of the surroundings; it ensured the consistency of the measured displacements with respect to regional trends over a large time span. The analysis and interpretation of the InSAR data have showed that: (i) the field area is characterized by a general long-term, gentle subsidence trend with no uplift evidence above the blind thrust, potentially ascribable to the growth of the buried anticline; a pronounced subsidence occurs in the direction of the Bologna city, outside the field area; (ii) the UGS activity does not affect the mean horizontal and vertical displacement velocities in the gas field area, which are coherent with the velocity range within the entire monitored domain; (iii) a strong correlation exists between the curve of the storage gas cumulative volumes and the historical series of the ground displacement above the reservoir, considering both the vertical and East-West planar components; (iv) the UGS-related, short-term, cyclical subsidence/uplift is limited to the field area; it is maximum in the center of the area while it dissipates near the field boundary.
The analysis and interpretation of the GNSS data have showed that: (i) the GNSS site is characterized by a long-term, gentle subsidence trend, in agreement with InSAR data; (ii) the estimated areal velocities are consistent with the NE-ward vergence of the Northern Apennines driven by the movement of the Eurasiatic plate (e.g., [35]); (iii) The short-term, seasonal vertical displacement of the GNSS station (placed close to the field boundary) is largely due to the variation in hydrological load rather than the UGS activities, resulting in the absence of an appreciable correlation between the curve of the storage gas cumulative volumes and the residual sinusoidal trend. The GNSS response is coherent with InSAR data, which show the decrease in the UGS-related cyclical vertical displacement towards the field boundary; (iv) a strong correlation exists between the UGS activities and the short-term, seasonal horizontal displacements; in particular, the GNSS station is more sensitive to the North-South planar component of movement due to its location at the south-eastern edge of the gas field area.
In conclusion, the integration between GNSS and InSAR data provides a full and coherent panorama of the ground movement in the monitored area.
Note that the present analysis investigated the effects of the storage activities and did not deeply analyze other variables such as seasonal temperature variations, groundwater withdrawal or large-scale natural subsidence/uplifting.

Discussion about the Geomechanical Simulation Results
The pseudo-elastic parameters adopted for the mechanical characterization of the reservoir and the cap rock were derived from laboratory tests (triaxial compression tests under isotropically consolidated undrained conditions), from the interpretation of well logs (sonic and density logs) and corroborated by data from the technical literature. The obtained Young's modulus values differ by almost one order of magnitude: from few GPa (values from the lab data interpretation) to 8-10 GPa (values from the log data interpretation). As documented by analogous case studies reported in the technical literature [11,19], the mechanical behavior of clastic formations at a depth of around 1000-1500 m ssl could exhibit an important (non-linear) influence of the strain on the formation stiffness, especially in case of soft rocks [37,38]. Experimental values of the deformation parameters (i.e., Young's modulus) increase as the imposed strain levels decrease from small strain (i.e., the values obtained from the stress-strain measurements in a rock mechanical test) to very small strain (i.e., values obtained from the interpretation of acoustic acquisition) [39,40]. Furthermore, if acoustic acquisitions are available at wellbore level, the scale effect can also be taken into consideration.
The discrepancy in the experimental values could generate uncertainty in defining the real deformation behavior of the formations; it is trivial to say that high variability of input parameters leads to high uncertainty in model prediction.
For the investigated case study, the real deformation behavior of the formation was identified via the back-analysis of the mechanical model. The formation stiffness that resulted from the calibration process is around three times higher than the values from lab tests, in agreement with the lithologies under analysis and with the in-situ stress condition during UGS operations [11,19]. As a theoretical explanation, during primary production, the deformation behavior of the normally consolidated formations is controlled by the loading static elastic modulus E I (in the case presented in this paper, E I are the values obtained from the lab tests). In this phase, the stress-strain path follows the critical state line (CSL). Due to the formation compaction caused by storage activities, the formations could become slightly over-consolidated and their elastic deformation response during the storage phase is usually ruled by the unloading/reloading static elastic modulus E II , around 2-4 times higher than E I .
The good match obtained between measured and simulated ground displacements assuming a constant unloading/reloading elastic modulus over almost two decades of UGS operation denotes a sound and stable system response in time.
Once calibrated, the geomechanical model represents an effective tool for forecasting future scenarios.

Conclusions
The paper proposes a fully-integrated and multi-physics approach for the analysis of the ground movement related to a gas storage system (eastern Po Plain area). The gas field is hosted in a clastic Pliocene sequence involved in a mixed trap. The trap combines the onlap of sand levels and a structural high associated with an NE-verging blind thrust that ends within the Pliocene sequence without propagating up to the Quaternary clastic sequence. The InSAR data excludes surface evidence of an active growth of the buried anticline during the entire monitored period.
The monitoring of the ground movements was achieved via both the SAR interferometry (PSP-IFSAR) and the GNSS continuous technologies. The two complementary technologies provided quantitative information about vertical and planar (North-South and East-West) components of the surface movements (with millimetric accuracy) and the subsidence/uplift areal extension of the UGS influence. The InSAR data measured the East-West planar and vertical components; the robustness of the adopted methodology is attested by the use of two different sensors (Radarsat-1/2 and Sentinel-1) to get independent and independently processed measures, and by the consistency of the results. On the other hand, the GNSS station measured the punctual displacement component associated with the UGS activity along the North-South direction, missed by the InSAR acquisition. The information from GNSS and InSAR data provided a full and coherent panorama of the ground movement in the monitored area: in particular, the identification of a regional trend was used for interpreting the effects of UGS operation.
The reliable ground movement monitoring plan and the detailed geological analysis, both at reservoir and regional scales, are key information for the geological, fluid-flow and geomechanical numerical simulation process. The multidisciplinary simulation approach deeply investigated and reproduced the main geological and structural features of the system together with the fluid-flow and poro-mechanical processes induced by UGS activities. In particular, almost two decades of displacement data effectively supported the calibration process of the geomechanical model, which showed a sound and stable mechanical response of the system in terms of induced subsidence/uplift for a given variation of pressure.
Furthermore, the two calibrated dynamic and geomechanical models represent the most reliable tool in forecasting the performance and in verifying the safety of the UGS system under different management conditions.