The Birth of an Underground Gas Storage in a Depleted Gas Reservoir—Results From Integrated Seismic and Ground Deformation Monitoring

Abstract

In the Po Valley (Italy), near Lodi, a depleted methane reservoir was recently converted into an underground gas storage (UGS) facility. We describe the new monitoring infrastructure that integrates seismic and ground deformation capabilities. We also present results obtained before and after the start of UGS operation, namely: (a) the so‐called “background,” for seismicity and deformation estimated over several decades before the UGS; (b) the “seismic baseline,” assessed using the new monitoring network over nearly 2 years before the gas injection began; and (c) the seismicity and deformation measured over the first 3 years of UGS operation. In practice, we observe the phenomena associated with the “birth” of a UGS with remarkable instrumental capabilities. Following three cycles of injection/extraction, about 30 events consistent with natural, tectonically related seismicity were located within 30 km of the UGS. Moreover, the observed uplift of about 2 cm is consistent with theoretical expectations of ground deformation. Our study confirms that UGS in depleted gas reservoirs, if well managed, can result in negligible, if any, human‐induced seismicity and limited ground surface deformation. It also shows the importance of measuring the undisturbed conditions prior to industrial activities over a sufficiently long period of time, to correctly interpret the phenomena observed later.

Introduction

Underground gas storage (UGS) facilities are strategic components of a gas supply infrastructure. They are needed for several reasons: to temporarily store natural gas imported from other places or countries; to compensate for seasonal fluctuations in gas demand; to ensure the so-called strategic gas reserve, which is drawn upon in case of extraordinary events (e.g., geopolitical crises, adverse weather periods, accidents). UGS operations involve reinjection and withdrawal of gas into underground reservoirs; this is usually done in annual cycles—extraction during the cold season and injection during the warm season—although some modern technologies allow much faster reversal of regimes, even on a daily basis. Deep aquifers, depleted gas fields, salt caverns, and deep mines are the most commonly used solutions for UGS.

All underground industrial activities—UGS is one of them—can change the stress field of the Earth's shallow crust and thus trigger, induce or generate earthquakes. In the case of existing faults, failures occur when the shear stress acting on the fault plane exceeds the Coulomb shear strength. Monitoring of microseismicity and ground surface deformation plays a key role in better understanding the physical mechanisms controlling induced seismicity and is a fundamental tool used by decision makers to decide whether to suspend, reduce or continue industrial operations. Recently, a comprehensive study on how faults can be reactivated by the operation of UGS has been published (KEM-01, 2021) with the aim of establishing a protocol to manage/limit induced seismicity. For more details on the physical mechanisms that can trigger seismicity and strategies for monitoring human activity underground, see also Grigoli et al. (2017) and the references therein.

As said, injection of fluids at seismogenic depths has the potential to perturb stress conditions and reactivate critically stressed faults (Ellsworth, 2013). However, the effect of reinjecting gas into a reservoir is much weaker than that of water or other incompressible fluids, which hydraulically transfer all overpressure into the rock. UGS in depleted fields is prevalent, accounting for about 80% of the world's working gas volume (Cornot-Gandolph, 2021), and scientific literature reports very few cases of felt-to-damaging earthquakes occurred during UGS operations in depleted reservoirs, over a total of nearly 700 plants (American Petroleum Institute (API), 2015; Cornot-Gandolph, 2021; Evans, 2008; Pipeline and Hazardous Materials Safety Administration (PHMSA), 2020; The American Petroleum Institute (API) et al., 2016; U. S. Department of Energy, 2016). Therefore, as far as seismicity is concerned, depleted gas fields are probably the safest option compared to other storage options (e.g., deep aquifers, salt caverns, deep mines), provided that proper preparatory studies are carried out prior to their implementation. Nano- or microseismicity can only be detected at the reservoir scale using downhole instrumentation or local networks designed specifically for this purpose (Cesca et al., 2014; Muntendam-Bos et al., 2022). Only a few UGS cases document induced seismicity and are briefly discussed here.

The Castor Project planned to reuse as a UGS facility the depleted Amposta oil field in the Gulf of Valencia, Spain, about 21 km offshore, in which the oil was initially trapped at a pressure of about 190 bar (Juanes et al., 2017). Gas injection began in September 2013 with an increasing injection pressure, which according to Juanes et al. (2017) started at around 177 bar and reached a maximum of around 185 bar. Shortly after the start of injection work, more than 1,000 weak events were recorded near the injection platform in about 40 days (Grigoli et al., 2017). UGS activity was halted due to public concerns. However, seismic activity continued for weeks after the end of the injections, culminating in a moment magnitude M_w 4.3 event on 4 October 2013, which proved to be the strongest induced earthquake ever associated with gas storage operations (Cesca et al., 2014, 2021). Subsequent studies showed that the seismic sequence was triggered by pore pressure diffusion from the injection sites, which activated a secondary fault of the Amposta fault system directly beneath the reservoir (Tang et al., 2018).

The Hutubi gas field is located in the southern Junggar Basin, adjacent to the seismically active region of the Tien Shan Mountains (China). The reservoir had an original gas pressure of 339.5 bar at a central depth of 3.5 km, which dropped to 137 bar after full exploitation and before conversion to UGS (Qiao et al., 2018). The UGS has been in operation since 2013 and it is the largest gas storage field in China with a storage capacity of 10.7 bcm (billion cubic meters). The UGS injection pressure varied from 13 to 28 MPa during its operation from June 2013 to December 2016 (Qiao et al., 2018; Zhou et al., 2019). About 2 years after UGS was put into operation, several earthquake sequences with a maximum local magnitude M_L 3.6 occurred in the region within 10 km of the gas field. The earthquakes followed an abrupt change in gas injection rate during the first two cycles of injection (Zhou et al., 2019) and occurred on faults not thought to be hydraulically connected to the reservoir formation (Jiang et al., 2020). According to Tang et al. (2018) and Jiang et al. (2020), these earthquakes were triggered by poroelastic stress disturbances caused by the gas injection.

The Gazli (Uzbekistan) case could also be included in this group (Suckale, 2010), but there is too much doubt whether these earthquakes were actually induced (National Research Council (NRC), 2013).

The case of the Bergemeer UGS is also worth mentioning, although only microseismicity was associated with UGS activity. The Bergermeer field is a natural gas deposit in the north-western part of the Netherlands. The field was in production from 1970 to 2007 and the reservoir pressure was originally 238 bar (Kraaijpoel et al., 2013). During the production phase, four seismic events were recorded by the national network: two events in 1994 (M_L 3.0 and M_L 3.2) and two events in 2001 (M_L 3.5 and M_L 3.2), which were attributed to the central fault (Van Eck et al., 2006). Muntendam-Bos et al. (2008) questioned the accuracy of the earthquake location, especially for the 1994 events. They argued that, because of the uncertainty in the earthquake location, the seismic events cannot be placed unambiguously on the internal fault but that it is equally possible that they occurred on a surrounding fault. After exploitation, the field was converted to a UGS facility: cushion gas injection started in 2010 and the storage system has been fully operational since spring 2015 with a working gas capacity of ∼6 bcm. The system is monitored with a new 6-level borehole geophone string. The monitoring array recorded hundreds of induced events that occurred either within or above the reservoir and had a maximum magnitude of M_L 0.7 (TAQA, 2015; Muntendam-Bos et al., 2022 and references therein). Seismic activity on the central fault ended in late 2014.

Another experience worth reporting is that of Silverii et al. (2021); the authors demonstrate the potential of using established software to simulate temporal variations in pore pressure and examine their relationship to recorded microseismicity (about 500 events over a two-and-a-half-year period). It is worth noting that the UGS in this case uses a saline aquifer and the seismicity background is not known: the USG location, although generally indicated to be in southeastern France, remains hidden, probably for reasons of industrial confidentiality. They conclude that “a part of the observed seismic patterns could be related to the UGS activity” but also that “at least a portion of the observed seismicity is most likely due to tectonic processes rather than human activities” and “more accurate models—and longer data time-series—are needed to confirm or reject this hypothesis.”

To address public concern about potential induced seismicity associated with the energy industry, government agencies in several countries have commissioned specific state-of-the-art studies (e.g., Evans, 2008; NRC, 2013), published new protocols (e.g., the Italian Guidelines for Monitoring Underground Activities Associated with Hydrocarbon Production and Storage, MiSE-UNMIG, 2014, hereafter ILG), or updated existing ones (Grigoli et al., 2017). These protocols are intended to take into account seismicity, ground deformation, and pore pressure in the vicinity of storage facilities. The key element of the existing protocols is the implementation of appropriate monitoring infrastructures tailored to the type of activities.

Today, induced seismicity is treated just as the natural one to the point that in the United States, the USGS has incorporated the hazard of induced seismicity into the national seismic hazard model (Petersen et al., 2016), while several European governments (e.g., the Netherlands, Switzerland, and Italy) require mining companies to install appropriate microseismic and ground deformation monitoring infrastructures. In Italy, for example, seismic monitoring is strongly recommended, although not yet mandatory, for new reservoir concessions according to regulations set by the Ministry of Environment.

In the Netherlands, for example, in addition to the previously reported case of UGS in Bergemeer, a large microseismic network with more than 60 borehole stations has been set up to monitor the industrial activities of gas production in the Groeningen field. In Italy, the Collalto UGS site (northern Italy) is monitored by a dedicated microseismic network (Priolo et al., 2015) consisting of 10 stations equipped with borehole seismometers. Another experience of gas storage monitoring in Italy is that conducted in the depleted gas field of Minerbio, where the local microseismic network consists of 8 stations, 4 of which are equipped with borehole sensors at a depth of 100–150 m (Carannante et al., 2020, 2021). Another monitoring experience described by Kwiatek et al. (2019) for the Enhanced Geothermal System near Helsinki (Finland) is not directly related to UGS but is relevant in the context of the present study. Fluid stimulation of deep geothermal wells near Helsinki is monitored in real-time by a telemetric network of 24 borehole seismographs, and near-real-time seismic monitoring is successfully used to modify stimulation parameters to keep the magnitude of induced earthquakes below M_w 2.0.

In this paper, we describe the new monitoring infrastructure set up at the Cornegliano Stoccaggio UGS, that integrates seismic and deformation measurements at the ground surface. Together with the gas production and storage fields in the Netherlands, this is one of the few cases in which, to our knowledge, monitoring was planned before the operational phase of the gas storage facilities. We will present both the detailed picture of seismicity and ground surface deformation prior to the start of gas injection—the so-called “background” for seismicity and deformation analysis, and the “baseline” for seismicity recorded by the new monitoring network before the start of storage operations—and what has been measured so far after the start of UGS operations. We will highlight the importance of estimating the seismic baseline as a reference for evaluating seismicity variations due to storage activities. Considering that the evidence of seismicity induced by UGS activities in the literature is limited to the first years of UGS operations (e.g., Cesca et al., 2021; Tang et al., 2018; Zhou et al., 2019), we are witnessing the phenomenology associated with the “birth” of a UGS with remarkable observational capabilities.

The UGS and the Integrated Monitoring Infrastructure

The Cornegliano Stoccaggio UGS

Cornegliano Stoccaggio is one of the newest UGS plants developed in Italy. It is operated by Ital Gas Storage S.p.A. and is located near the village of Cornegliano Laudense, a few kilometers from the city of Lodi and about 30 km south of Milan (Figure 1).

[IMAGE OMITTED. SEE PDF]

This part of the Po Valley is one of the less seismic areas in Italy (Figure 1a); most of the tectonic deformation occurs further eastward, in the Friuli Plain and the Eastern Alps, and southward, along the Apennine chain (Devoti et al., 2010, 2017; Nocquet, 2012; Rovida et al., 2021). Nevertheless, seismogenic faults are mapped beneath the alluvial sediments; they are due to the still ongoing crustal shortening between the southern Alpine front and the northern Apennine front, which are moving toward each other at an estimated average rate of 1–3 mm/year in a roughly N-S direction (Figures 1b–1d, from Martelli et al., 2016). The thickness of alluvial sediments deposited by the Po River and its tributaries reaches a maximum of about 600 m in the eastern part of the Po Plain toward the Adriatic Sea and is approximately 150 m near the UGS (Bigi et al., 1990; Carcano and Piccin, 2002; Ogniben et al., 1975). This cover prevents the identification of active faults by surface geological-geomorphological surveys and makes it difficult to determine the sources of major historical earthquakes. It also has strong effects on seismic wave propagation with known phenomena such as: amplification of ground motion due to seismic wave trapping and high-velocity contrasts; lower resonant frequency and longer earthquake duration; scattering and generation of surface waves due to the presence of buried structures (Barnaba et al., 2014; Bragato et al., 2011; Klin et al., 2019; Laurenzano et al., 2017; Mascandola et al., 2017; Tondi et al., 2019). This was the case for the 1951 Caviaga earthquakes, originally interpreted as events triggered by oil and gas production and recently revised as deep earthquakes of tectonic origin (Caciagli et al., 2015; Peruzza et al., 2021).

Similar to the other dozen UGSs operated in the Po Valley, Cornegliano Stoccaggio uses a depleted gas reservoir that was discovered in 1951 and used for gas production until 1995. The reservoir is located at the top of an anticline at a depth of 1.3–1.4 km; the concession covers an area of about 24 km² (MiSE-DGISSEG-UNMIG, 2021). The reservoir and sealing rocks are Pliocene sands (Caviaga Sand Formation) and Lower Pleistocene clays (Santerno Clay), respectively (Boccaletti et al., 2011; Fantoni & Franciosi, 2010).

The UGS facility consists of two main sites, called “clusters” (see Figure 1c), located near the town of Cornegliano Laudense and separated by about 1.5 km. Cluster A houses the headquarters, the main processing facilities, and has six active wells reaching the reservoir; in Cluster B, four wells were active until 30 September 2020, when two more wells were opened. The natural gas reservoir has an estimated storage capacity of approximately 1.8 billion Sm3—Sm3 stands for Standard cubic meters, that means the amount of gas that would occupy a volume of one cubic meter under standard conditions, that is, at temperature of 15 degrees Celsius and under an absolute pressure of 1.01325 bar—at the original confining pressure of approximately 160 bar, measured at depth of approximately 1,300 m. In 2018, the pressure within the reservoir was approximately 114 bar at a depth of approximately 1,400 m, with a gradient of 0.081 bar/m. Storage operations started in December 2018.

During the design phase of the UGS facility, 3D geomechanical modeling was commissioned to assess the potential response of the reservoir to the storage operations and to guide the OGS operator toward optimal management of the storage operations (Ital Gas Storage, pers. comm.). Some outcomes of the modeling results are given in the discussion chapter and compared with the results of our observations.

Following the ILG for monitoring underground activities (MiSE-UNMIG, 2014), the company decided to carry out integrated seismic and geodetic monitoring on a voluntary basis. Two Italian public research institutes were entrusted with the development and management of the monitoring system, namely OGS-CRS (National Institute of Oceanography and Applied Geophysics—Center for Seismological Research) and CNR-IREA (National Research Council—Institute for Electromagnetic Sensing of Environment). OGS-CRS is responsible for seismic monitoring and coordinates all monitoring activities, while CNR-IREA is responsible for geodetic monitoring. As recommended in the Ministry's guidelines, seismic activity should be analyzed over a period of at least 1 year prior to the start of storage operations in order to define a baseline for seismicity.

The New Monitoring Infrastructure

The Cornegliano Laudense Seismic Network (RSCL) was designed with the objective of monitoring and investigating both microseismicity potentially induced by storage activities within or near the reservoir and natural seismicity in a wider, surrounding area. In accordance with the ILG requirements, the network combines both the instrumental characteristics of high dynamic range suitable to record moderate to strong earthquakes from nearby seismogenic sources (known or suspected) and the high sensitivity required to monitor the gas reservoir and surrounding geologic structures with the best-resolving power.

RSCL (registered with code OL at FDSN-International Federation of Digital Seismograph Networks, ) currently consists of 10 seismic stations (Figures 1b and 1c, Table 1, and Text S1 in Supporting Information S1). The network geometry is pivoted on the most productive part of the reservoir (Figures 1b and 1c): an internal core of six stations (OL01-OL06) is arranged in a pentagon with one station in the center, and the maximum distance between stations is 4–5 km. A group of four stations (OL07-OL10) is arranged in an outer quadrilateral with E-W and N-S aligned axes at distances between 10 and 20 km from the central core; they provide a seamless transition toward the national and regional seismic networks (RSNC, the Italian national seismic network managed by the National Institute of Geophysics and Volcanology (INGV) Seismological Data Centre, 2006; the RSNI-Regional Seismic Network of Northwest Italy managed by the Genova University; and the OGS network in Northeast Italy, Bragato et al., 2021; Priolo et al., 2005). All the OL stations are connected in real-time to the OGS acquisition center. RSCL was originally designed with nine stations; the tenth station OL10 was added later and put into operation in 2022. Therefore, station OL10 does not contribute to the results described below. Table 1 summarizes the information on RSCL station locations.

Table 1 Location of the Rete Sismica di Cornegliano Laudense Stations

Code	Lat	Lon	Elevation above sea level (m)	Location (municipality, province code)
OL01	45.291133	9.464850	72	UGS Cluster A (Cornegliano Laudense, LO)
OL02	45.266917	9.485002	71	Massalengo (LO)
OL03	45.266723	9.453194	70	Cascina Castagna (Pieve Fassiraga, LO)
OL04	45.290477	9.435672	78	Cascina Taietta (Lodi Vecchio, LO)
OL05	45.293630	9.526859	75	Cascina Cavrigo (Lodi, LO)
OL06	45.309150	9.464777	77	Cascina Bracca (Lodi, LO)
OL07	45.291950	9.325655	77	Gugnano (Bascapè, PV)
OL08	45.362077	9.538975	75	Dovera (CR)
OL09	45.169950	9.440460	70	Miradolo Terme (PV)
OL10	45.273900	9.604020	70	Cavenago D’Adda (LO)

All stations are equipped with the same instruments manufactured by Guralp Systems Ltd. consisting of a Radian digital broadband borehole velocimeter (at a depth of about 75 m), a Fortis compact broadband accelerometer (on the surface) and a Minimus 24-bit digitizer. The acquired seismic data is transmitted in real-time via the SeedLink protocol to the processing and storage system. The real-time acquisition also includes some stations managed by OGS to monitor Northeastern Italy. Continuous waveform data are permanently stored in the OASIS system (), the OGS archive for instrumental seismological data (Priolo et al., 2012); here, all waveform data are freely available with a delay of one day from their acquisition.

Data acquisition, processing, and analysis procedures are based on the Antelope software system developed by BRTT (), supplemented by some procedures and functions developed by OGS staff. The monitoring system has both real-time and off-line data processing and uses similar procedures to those adopted by OGS for monitoring the Collalto UGS, as described in Priolo et al. (2015).

For event localization we use the software Hypo 71 (Lee & Lahr, 1975) and Hypoellipse (Lahr, 1999). The former is used to ensure compatibility with the localizations of the regional OGS network and to check the quality of the localizations obtained in the arrival time selection phase. Hypoellipse is used to relocate all events once the arrival times of phases P and S have been confirmed. The local magnitude (M_L) is calculated using the attenuation formula of Bragato and Tento (2005).

Further information on the seismic monitoring can be found in Text S1.1 in Supporting Information S1.

Geodetic monitoring is carried out through two components, one based on an in-situ detection system with Global Navigation Satellite System (GNSS) stations and the other on satellite remote sensing (i.e., Differential Interferometry Synthetic Aperture Radar—DInSAR technology). The low tectonic deformation of the region where the UGS site is located allows the detection of possible deformation episodes related to the UGS activity, since the amplitude of the detectable surface deformation exceeds the one of the background noise and can be well distinguished from the signal caused by seasonal variations.

A GNSS station was installed at site OL01, near the homonymous seismic station, to obtain direct ground displacement data on potential anomalous deformation phenomena in the area above the UGS. It was installed in early 2017 and it consists of a Topcon NET-G5 receiver and a Topcon CR G5 choke ring antenna. The antenna is coupled to a stainless steel mast with a SCIGN GNSS adapter located at the top of an 8 m long iron pipe. The iron pipe, with a diameter of 88.9 mm and a thickness of 8 mm, descends 6 m deep into the ground and is anchored on the surface to a small cylindrical concrete pillar about 80 cm high. All power and data transmission systems are shared with the seismic station.

Since our goal is to monitor the potential deformation of the UGS, we establish a local reference frame using four stations of the interregional satellite positioning service SPIN3 GNSS (), namely CREA, CREM, MIL2 (which replaced MILA in 2019), and PAVI (see Figure 1b). All these four reference stations are located within ∼30 km of OL01 and have been in operation for more than a decade. The position of the “target” station OL01 is calculated with respect to this reference frame. The coordinates of the GNSS stations used in this study are listed in Text S7 and Table S7 in Supporting Information S1. Data from permanent GNSS stations are processed using the GAMIT/GLOBK software package (Herring et al., 2018).

Since the time-interval before the start of the UGS activities is short (less than 2 years), it is not possible to estimate OL01 velocity from the time series, ruling out the eventual effect of the storage activities. Hence, we assumed that OL01 velocity can be approximated to the average velocity of the reference stations. We calculated it, and then we subtracted it for each component day by day from the OL01 time series. In this way, the residual signal of OL01 represents the local behavior of the station.

Further information on the GNSS monitoring can be found in Text S1.2 in Supporting Information S1.

The DInSAR technology has been applied to spaceborne SAR data to measure surface deformations estimated with millimeter accuracy in the area of the Cornegliano Laudense UGS. Both archived and current SAR data are used to analyze ground displacements before and after the start of storage activities. The surface deformation analysis covers an area of approximately 190 × 70 km², which includes the UGS site.

In particular, the estimation of ground deformations prior to the construction of the UGS is performed using the entire SAR data archives collected from both ascending and descending orbits by the European Space Agency (ESA) ERS-1/2 and ENVISAT sensors, during the 1993–2010 time period. Moreover, for the most recent operational geodetic monitoring, we analyze the SAR data acquired by the European Union's Copernicus Earth Observation Program's Sentinel-1 (S-1) constellation, collected between March 2015 and December 2021 from both ascending and descending orbits and with a repeat time of 6 days. The DInSAR processing has been performed using the parallel computing resources of the Amazon Web Services cloud computing environment, which is appropriately equipped and configured.

Specifically, the DInSAR data were processed using the SBAS approach (Berardino et al., 2002), which is known as Parallel SBAS (P-SBAS), in its parallel version, and can effectively exploit distributed computing infrastructures (cluster, grid, cloud) (Casu et al., 2014; Zinno et al., 2020). We also remark that SBAS is an algorithm allowing the analysis of the spatio-temporal characteristics of the observed deformation phenomena by generating maps of surface deformation and time series with medium and full spatial resolution (a few tens of meters and meters, respectively), with a high spatial density of measurement points and a high temporal frequency given by the repetition time of the SAR sensor. The SBAS technique can be used to generate time series of deformation and the corresponding maps of average velocity with an accuracy of about 5–10 mm and 1–2 mm/year, respectively (Bonano et al., 2013; Casu et al., 2006; Manunta et al., 2019).

In addition to the possibility to carry out analysis based on the SAR data acquired by the currently operating SAR sensors (as for the mentioned Sentinel-1 constellation), one very relevant capability of the SBAS (and P-SBAS) technique is the generation of very long historical deformation series, lasting about 20 years. This is possible thanks to the exploitation of SAR data relevant to the same ground scene, acquired by different but geometrically compatible sensors, as in the case of the ESA ERS-1, ERS-2, and ENVISAT sensors (Bonano et al., 2012). Moreover, the SBAS multi-sensor approach is also able to generate products with different spatial resolutions (Lanari et al., 2004). Accordingly, the SBAS approach can be applied to study the currently ongoing slow deformation phenomena at both regional and local scales and also to reconstruct their past evolution (baseline) by using the large archives of SAR data collected by the ERS-ENVISAT radar systems between 1992 and 2010 (Bonano et al., 2012; Pepe et al., 2005).

We further remark that only the SAR pixels where the time series of deformation and the corresponding average velocities are reliably retrieved are considered; these pixels are referred to as coherent points and are properly geocoded. Moreover, the mean deformation information is also represented through false color maps corresponding to ground displacements and, following the geocoding step, they can be overlaid with other maps of the study area.

Further information on the DInSAR monitoring can be found in Text S1.3 in Supporting Information S1.

Three spatial domains to be monitored have been defined, with geometry derived from that of the reservoir. We represent the most productive part of the reservoir as a thin, ellipsoidal polygon that is about 2 km by 4 km and 1.4 km deep (dashed area in Figures 1b and 1c). The monitoring domains are the volumes surrounding the reservoir and are defined as follows (Figure 1):

• the Inner Domain (ID) extends 3 km from the reservoir both laterally and at depth, and extends vertically to the ground surface (see also Text S2 and Figure S5 in Supporting Information S1);

• the Extended Domain (ED) extends from ID for an additional 10 km, in the same manner as ID defined above;

• the Outer Domain (OD), the largest domain, is defined as a spherical volume with a radius of 30 km, centered on a point within the reservoir below the gas storage facility (near station OL01).

Both ID and ED are defined in accordance with the ILG (MiSE-UNMIG, 2014; see also Text S2 and Figure S5 in Supporting Information S1) and have a direct bearing on the interpretation of the detected seismicity. In particular, ID defines “… the volume within which induced seismicity and ground deformation could be potentially caused by anthropic activities. It represents the reference volume within which seismicity and ground deformation will be monitored, analyzed and, when possible, identified with maximum sensitivity.”, while ED “… is used to better constrain monitoring and to help the interpretation of the measured quantities (i.e.,: seismicity, deformation, and pore pressure) within the existing structural and geological background.” We introduced OD as a third, larger monitoring domain, to capture the natural deep seismicity that characterizes this part of the Po Valley.

Comprehensive information on the integrated seismic and geodetic monitoring system, and public access to monitoring reports and data, is provided through a dedicated website (). Scientific reports are released every 6 months and sent to the concession holder and controlling authorities at the completion of each injection/extraction cycle.

Results of the Integrated Monitoring

The integrated monitoring aims to provide detailed information on both real-time and/or quasi-real-time seismicity in predefined areas surrounding the UGS and surface deformation in the area above and around the UGS, with a delay within a couple of weeks. According to ILG, a baseline was established prior to the start of storage operations, using the same instrumental infrastructure as for operational monitoring. Such a baseline extends back as far as possible and includes at least 1 year of records with the final monitoring layout.

We describe all the results of the integrated monitoring of the Cornegliano Laudense UGS according to the following scheme (a graphical representation is given in Figure S9 in Supporting Information S1):

• the seismicity before the start of UGS operations, which includes:

  • the instrumental seismicity recorded before the RSCL deployment, critically reviewed and commented together with the major historical earthquakes—we refer to it as “background seismicity”;

  • the seismicity recorded by RSCL, that is, the same seismic network developed for operational monitoring, before the start of storage operations (i.e., for almost 2 years, from 1 January 2017, to 11 December 2018)—we refer to it as the “baseline seismicity”;

• the seismicity recorded by RSCL during the operational monitoring, that is, for 3 years after the start of storage operation (period 12 December 2018, to 31 December 2021)—we refer to it as “current seismicity”;

• the ground surface deformation before the start of UGS operations, which includes:

  • the interferometric analysis of archived SAR data collected by the constellations no longer in operation, the so-called “first-generation” SAR, ERS-1/2 and ENVISAT, for the 18-year period 1993–2010—we refer to it as “background surface deformation”;

  • the interferometric analysis of data collected by the currently operating SAR satellite constellation, namely Sentinel-1, before the start of UGS operations (i.e., for the nearly 4-year period 2015–2018)—we refer to it as “baseline surface deformation”;

  • c.the analysis of GNSS data from station OL01 for a period of about 2 years before the start of UGS operations (i.e., for the period from 1 January 2017, to 11 December 2018);

• the ground surface deformation estimated during operational monitoring, that is, for 3 years after the start of UGS operations (period from 12 December 2018, to 31 December 2021)—we refer to it as “current surface deformation,” which includes:

  • the interferometric analysis of the Sentinel-1 satellite data;

  • the analysis of the GNSS data from the OL01 station.

For each field of investigation (i.e., either seismicity or surface deformation), specific study areas are defined that include the integrated monitoring domains, so that local observations can be constrained at the regional scale.

Seismic Monitoring

Background Seismicity

The analysis of background seismicity focuses on a square area about 70 km wide around the UGS, corresponding approximately to OD. It considers the major historical earthquakes and the entire instrumental period up to the start of the reservoir operations. Peruzza et al. (2021) recently published a critical analysis of the instrumental data, and uniformly relocated the seismicity in the study area from 1951 to 2019, that is, for about 70 years of nonuniform observations. Figure 2 shows the major historical earthquakes as reported by CPTI15 (Rovida et al., 2021) and the relocated data set of Peruzza et al. (2021) from 1980 up to 2017.

[IMAGE OMITTED. SEE PDF]

Damaging historical earthquakes in the study area mainly affected the northern sector, with the most important one occurring in 1802 (M_w = 5.6). The events closest to Lodi are those of 1642, 1781, 1786, and 1951, the latter also known as the Caviaga earthquake. The 1951 earthquake sequence (main events on 15 May at 22:54 UTC, M_w = 5.4, and on 16 May at 2:27 UTC, M_w = 4.6, see Caloi et al., 1956; Caciagli et al., 2015) is extremely important from the point of view of induced seismicity because it occurred in an area rich in hydrocarbon reserves that were intensively exploited after World War II and that was considered aseismic at the time of the earthquake. Only these two moderate earthquakes occurred during the instrumental period within a distance of 30 km from the Cornegliano Laudense reservoir, and have been considered for decades to be caused by hydrocarbon production. The 15 May 1951 earthquake is still included in worldwide compilations of induced seismicity and is listed as the only M ≥ 5 earthquake triggered by gas production in Europe (Foulger et al., 2018). Caciagli et al. (2015) relocated the two main events using phase readings reported in the International Seismological Summary (ISS, now ISC, 2022; event identifiers nos. 894347 and 894349 for the main and largest aftershocks, respectively; note that different parameterizations are given in each quoted source).

Peruzza et al. (2021) also relocated these two events (using the original phase readings published by Caloi et al., 1956) and calculated the focal mechanisms using the first motion polarities. Although large localization errors exist and some ambiguity remains regarding the stress regime and causative faults, these events are most likely deep events representing active sources in the lower crust where the northern Apennines and southern Alps meet beneath the Po Valley.

Since the 1980s, very few earthquakes have been identified and located near the reservoir; only two events of the newly located subset are in the ED (Figure 2, gray-colored events; the fact that three rather than two events appear within ED in Figure 2d is only a perspective effect of the angles of view taken), namely the earthquakes of 24 December 1996, 6:02 UTC, and 2 February 2002, 3:43 UTC, with recalculated depths of 9.1 and 4.9 km, respectively. The magnitudes reported by the national authoritative agency INGV are duration magnitude Md 3.0 and Md 2.7, respectively (ISC, 2022; event identifiers Nos. 1003254 and 3861366, respectively); these values were lowered in the latest revised data sets published by INGV. Note, however, that the statistical errors in localization may not be representative of the actual uncertainties because the nearest recording stations are more than 50 km away in both cases.

Baseline Seismicity

RSCL seismic monitoring became fully operational on 1 January 2017, and for nearly 2 years (until 12 December 2018) the recordings have been considered representative of natural seismicity, prior to the start of UGS filling. During this period, RSCL has located 12 events within the monitored area (Table 2, Figure 2, colored dots), with local magnitude (M_L) ranging from 0.9 to 2.2. Four events are deep (30–40 km) and are most likely due to natural tectonic deformation. Two events with a depth of 5–10 km are located in ED. The remaining six events are very shallow (1–4 km depth), and four of them are located within ID. The analysis of baseline seismicity with the new network, which is more sensitive than the existing ones, confirms that the RSCL target area has a very weak seismicity, probably due to some existing shallow active tectonic structures.

Table 2 Earthquakes Located by Rete Sismica di Cornegliano Laudense (RSCL) From 1 January 2017 to 31 December 2021

Id	Area	Date	Time	Lat	Lon	Depth	ML	Q	Gap	No	Ns	errh1	errh2	errz	rms	Distance	Location
1	ED	2017/03/06	20:12:37.93	45.370	9.371	9.1	2.2	D	251	7	0	1.77	4.16	99.00	0.41	12.6	Dresano
2	OD	2017/05/23	04:16:32.54	45.448	9.246	38.0	1.9	A	175	16	7	0.44	1.05	0.65	0.15	43.1	Milano
3	ID	2017/07/11	22:25:03.54	45.306	9.422	1.4	0.9	D	294	8	3	0.11	0.17	99.00	0.07	1.9	Lodi Vecchio
4	ID	2017/07/12	00:17:53.87	45.304	9.419	2.2	1.1	A	156	15	6	0.14	0.20	0.14	0.18	2.2	Lodi Vecchio
5	ID	2017/07/12	00:18:30.95	45.304	9.435	1.3	–	A	256	8	3	0.05	0.25	0.04	0.04	1.1	Lodi Vecchio
6	ED	2017/10/20	21:37:06.65	45.325	9.544	3.2	1.2	C	208	5	2	0.12	0.71	4.97	0.10	6.2	Lodi
7	ID	2018/01/27	04:01:44.44	45.311	9.435	1.5	0.8	A	278	11	5	0.06	0.09	1.16	0.11	1.8	Tavazzano con V.
8	ED	2018/05/27	20:01:35.35	45.257	9.456	5.7	0.7	A	150	13	6	0.14	0.19	0.19	0.14	4.9	Pieve Fissiraga
9	OD	2018/06/18	21:45:20.87	45.276	9.253	33.3	1.2	A	191	18	9	0.40	0.65	0.47	0.37	35.3	Vidigulfo
10	OD	2018/06/20	00:49:07.48	45.276	9.238	34.0	1.3	A	301	15	7	0.59	0.72	0.37	0.23	36.4	Vidigulfo
11	OD	2018/10/30	05:03:37.61	45.383	9.572	34.0	1.9	A	153	22	11	0.39	0.79	0.72	0.17	34.9	Palazzo Pignano
12	ED	2018/11/21	03:26:44.47	45.308	9.532	3.1	2.1	A	284	12	6	0.18	0.22	0.38	0.19	4.4	Lodi
13	OD	2018/12/23	15:06:12.55	45.440	9.583	35.6	1.6	B	197	14	7	0.29	2.29	0.93	0.00	38.8	Vailate
14	OD	2019/04/25	23:11:31.35	45.355	9.605	32.0	1.5	A	255	21	9	0.36	0.59	0.49	0.10	32.3	Bagnolo Cremasco
15	ID	2019/06/12	10:24:02.07	45.302	9.502	1.7	1.0	B	191	10	5	0.11	0.17	1.50	0.16	1.8	Lodi
16	ED	2019/06/21	15:22:39.47	45.308	9.535	2.2	1.0	A	290	10	5	0.10	0.22	0.08	0.21	4.3	Lodi
17	ED	2019/06/21	15:49:01.16	45.307	9.529	2.3	0.9	A	179	10	5	0.08	0.16	0.17	0.16	3.9	Lodi
18	OD	2019/06/30	12:01:20.47	45.332	9.312	35.0	2.0	A	115	30	14	0.20	0.46	0.42	0.15	35.3	Cerro al Lambro
19	ED	2019/07/13	03:01:35.97	45.300	9.528	2.5	1.2	A	180	16	8	0.10	0.15	0.06	0.16	3.5	Lodi
20	ED	2019/08/17	09:40:07.14	45.305	9.649	4.2	2.0	A	222	17	7	0.71	0.85	1.02	0.23	13.0	Credera Rubbiano
21	OD	2020/10/28	09:28:30.08	45.453	9.467	12.2	1.7	D	177	5	2	0.58	3.84	99.00	0.21	20.5	Truccazzano
22	OD	2020/12/17	15:59:23.86	45.445	9.169	41.6	3.3	A	224	25	9	0.83	0.97	1.08	0.24	48.5	Milano
23	OD	2021/01/17	10:27:25.98	45.473	9.390	37.0	2.7	A	124	31	15	0.18	0.20	0.22	0.27	40.9	Settala
24	OD	2021/02/21	08:54:21.16	45.212	9.112	31.8	1.7	B	257	15	7	0.72	1.45	0.71	0.22	41.0	Pavia
25	ID	2021/08/14	06:26:48.34	45.303	9.436	1.2	0.2	A	252	5	2	0.08	0.48	0.10	0.01	1.0	Lodi Vecchio
26	ID	2021/08/14	06:27:10.49	45.304	9.427	1.9	0.5	A	259	11	5	0.17	0.33	0.49	0.14	1.6	Lodi Vecchio
27	OD	2021/10/11	23:18:21.81	45.306	9.246	23.9	2.9	A	179	27	12	0.14	0.23	0.20	0.47	27.3	Landriano
28	OD	2021/11/06	11:00:57.57	45.496	9.581	34.95	1.7	A	113	20	9	0.27	0.41	0.57	0.14	41.29	Casirate d'Adda
29	ED	2021/12/21	23:37:57.32	45.313	9.536	2.5	1.6	A	195	16	8	0.14	0.46	0.18	0.25	4.7	Lodi
30	ED	2021/12/22	05:17:15.26	45.312	9.533	2.4	2.2	A	189	16	8	0.09	0.17	0.07	0.14	4.5	Lodi

During this initial operational phase, we also verified that the RSCL event detection procedures were working as expected. A remarkable situation arose when we compared the events detected by RSCL with those listed in the INGV national seismic bulletins. During the same period, the national seismic network detected four events within 32 km of the reservoir. These events have M_L between 1.2 and 1.9, which is lower than the completeness magnitude (Mc) of about 2.2–2.5 estimated by Schorlemmer et al. (2010) in this area. For two events, the depth was set at 6 and 10 km, while for the others it was estimated at 32 and 45 km. Only one of these four events was detected by RSCL, namely event M_L 1.9 on 23 May 2017 at 04:16:32. A more detailed analysis of the data used for the locations of the national network confirms that the three events not detected by RSCL were fake locations, due to incorrect phase readings or association, for stations that are either too far from the site or with poor azimuthal coverage. A more detailed analysis of this problem can be found in Text S3, Tables S1 and S2, and Figure S6 in Supporting Information S1.

Current Seismicity

Since the start of storage operations on 12 December 2018, we have located 18 seismic events (see Table 2 and Figure 2) with M_L in the range of 0.2–3.3, with the strongest one occurred near Milan on 17 December 2020. As in the pre-storage period, some events occurred within the lower crust (20–42 km) and outside ID. One event is characterized by a depth of about 12 km, but is located in the OD. The remaining 3 events are located in the ID with depths in the 1–2.5 km range. The 6 events within ED are at shallow depths (1–5 km). No events occurred within the reservoir volume.

During the same period (see Table S2 in Supporting Information S1), the Italian National Seismic Network recorded 9 events (about half of the 17 detected by RSCL) with M_L ranging from 1.7 to 3.0, all located also by RSCL.

Surface Deformation Monitoring

GNSS analysis results are time series of stations' position in a selected reference frame or, equivalently, displacement with respect to a reference position, based on a daily data set sampled every 30 s. Figure S7 (Text S4) in Supporting Information S1 shows the time series we obtained for the period since 2015: they show coherent motions among themselves in agreement with previous work (e.g., Devoti et al., 2017). The results of the analysis are presented below.

Analysis of surface deformation in the vicinity of the Cornegliano Laudense UGS was performed using SAR data from both decommissioned and currently operating spaceborne radar systems to provide an estimate of relative ground displacements for both a comparatively long period before the storage activities begun and to monitor the current trend from the start of storage activities onwards. Deformation time-series retrieved from satellite data are compared with those estimated by some GNSS stations located in the analyzed area to verify the accuracy of the SAR measurements.

Background Surface Deformation: ERS-1/2 and ENVISAT SAR Data

To estimate the surface deformation of the ground before the onset of storage activities (the so-called background surface deformation), we use the total archive of SAR data collected by the ERS-1/2 and ENVISAT sensors of the ESA, from ascending and descending orbits over the area of interest. These data cover a period of more than 15 years, from 1993 to 2010, and our analysis covers an area of about 95 × 60 km² encompassing the UGS site (Figure 3). The data set obtained from the descending orbit (track 208) consists of 141 satellite images covering the period from May 1993 to September 2010; the data set obtained from the ascending orbit consists of 76 satellite images covering the period from June 1995 to October 2010. Table 3 summarizes the main characteristics of the ERS-1/2 and ENVISAT SAR data sets for the two orbits.

[IMAGE OMITTED. SEE PDF]

Table 3 Main Common Features and Cardinality of the Evaluated Synthetic Aperture Radar Data Sets

	ERS-1/2 and ENVISAT	S-1
Wavelength	5.6 cm	5.56 cm
Nominal spatial resolution (azimuth, range)	∼4 m × 20 m	∼20 m × 4 m
Spatial resolution of the interferometric products (lat, lon)	∼80 m × 80 m	∼30 m × 30 m
Spatial extent (lat, lon)	∼95 km × 60 km	∼190 km × 70 km
Observation period descending/ascending	May 1993–September 2010/June 1995–October 2010	March 2015–December 2021
No. of acquisitions descending/ascending	141/76	359/344
No. of generated interferograms descending/ascending	412/213	1,006/982

Baseline and Current Surface Deformation: Sentinel-1 SAR Data

The Sentinel-1 (S-1) satellite constellation was developed as part of the European Copernicus program. Interferometric analysis uses data from SAR, collected from both ascending and descending orbits, and has been used for both baseline estimation and operational monitoring. S-1 consists of two twin sensors, S-1A and S-1B, launched in April 2014 and April 2016, respectively, which allow data to be collected over the area of interest every 6 days. In this paper, we present the results for the period from March 2015 to December 2021, covering both baseline and current surface deformation, obtained by analyzing 344 S-1 images from the ascending orbit and 359 from the descending orbit. The DInSAR analysis covers an area of about 190 × 70 km² encompassing the UGS site (Figure 4). Table 3 summarizes the main characteristics of the S-1 SAR data sets for the two orbits. Figure 4 also shows the comparisons between the vertical and east-west components of the GNSS and DInSAR time series at three selected sites far from the UGS; the standard deviations of the differences between the two measurements are less than 0.5 cm on average.

[IMAGE OMITTED. SEE PDF]

Discussion

The measurement obtained so far by the new RSCL monitoring infrastructure shows a very scarce seismicity before and after the start of storage operations. Figure 5 compares the gas storage activity with the recorded seismicity. Panel (a) shows the timeline of the seismicity recorded by RSCL with respect to the hypocentral distance of the events from the reservoir. The events shown are those listed in Table 2. The size of the symbols is proportional to the local magnitude of the events. The other two panels show the recorded seismicity compared to (b) the average pressure of gas in the reservoir measured at the Cluster A wellhead and (c) the normalized volume of gas stored in the reservoir during UGS activity, respectively. The volume values are intentionally omitted because they are confidential, commercially sensitive data. The activity data are average daily values, although Ital Gas Storage provides OGS with hourly sampling data.

[IMAGE OMITTED. SEE PDF]

The diagrams show a sequence of alternating phases during UGS initial filling aimed at reducing local stress accumulation. The maximum pressure at the wellheads does not exceed the original confining pressure at the reservoir level (160 bar). The volume curve shows about two and a half storage cycles—three injections from spring to late fall and withdrawals in winter—and a gradually increasing amount of gas stored in the reservoir. We recall that the initial gas volume mainly contributes to the formation of cushion gas, that is, the quantity of gas permanently stored in a natural gas storage reservoir, which is necessary for the efficiency of injections and withdrawals.

Basic statistical analyses of the earthquake data presented in Figure 5 also show that there is no significant change in earthquake intertimes, and distances from the reservoir before and after the start of UGS operations (baseline duration 716 days; storage activities, 1,109 days, see Text S9 and Figure S9 in Supporting Information S1). The very weak seismicity recorded in the ID is statistically at the same distances from the reservoir during the baseline and injection/extraction periods (95% confidence level respectively 4.39 ± 2.43 km before 18 December 2018 and 4.22 ± 2.47 after). Similar considerations arise from the time series (details in Text S8 in Supporting Information S1). However, due to the very scarce local seismicity, we acknowledge that we do not have sufficient statistical samples to produce reliable statistics or to calculate meaningful seismological parameters such as the b-value and completeness magnitude, and to perform more sophisticated analyses (e.g., declustering). Unfortunately, the few events recorded in the vicinity of the deposit cannot be further interpreted seismologically, for example, with regard to the fault mechanism or the estimation of the moment tensor, as these events were very weak and were only recorded by a few stations. Moreover, it should be noted that most of the weak events detected by RSCL were not detected by the surrounding networks; conversely, some events reported by the National Seismic Network in the vicinity of Cornegliano Laudense, which do not appear in the RSCL signals, turned out to be false events.

In summary, from the analysis of seismicity:

• -

  Number, magnitude, depth, and distribution after December 2018 events are consistent with the general characteristics of background seismicity and are similar to the baseline measured by RSCL;

• -

  A few events occur in the ID during UGS activity, but they are consistent with the main characteristics of the previously recorded seismicity;

• -

  No events are located within the UGS reservoir.

These observations are in agreement with the results obtained by the 3D geomechanical modeling committed by the concession holder during the design phase (Ital Gas Storage, pers. comm.). The UGS operations were simulated over 25 years to allow the best storage management. The modeling suggests that, if the maximum operating pressure does not exceed the original confining pressure (i.e., approximately 160 bar): (a) no fractures occur in the caprock or within the reservoir, during the injection/withdrawal cycles; (b) the existing faults are excited well below the cracking conditions. The enhanced seismic monitoring confirms that no ruptures occurred during the first storage cycles within the reservoir or in the caprock.

It is quite difficult to draw a comparison between the UGS Cornegliano Laudense and some of the few documented cases, such as the UGS Bergemeer. In Bergermeer, the production from the gas field ceased in 2007 and in the following years preparations were made to convert the field to a UGS system. The injection of cushion gas began in 2010 and the storage system has been fully operational since spring 2015. In Cornegliano Laudense, production of the gas field ended in 1995 and injection began in 2018, so a first major difference is that in Cornegliano Laudense UGS the stress field in the reservoir had a period of about 25 years to rebalance after production. The analysis of seismicity during the first 5 years of UGS operation at Bergermeer (TAQA, 2016), which includes the period of cushion gas formation, reports about three hundred events with a maximum magnitude M_W 0.9. It further states that relocation based on waveform similarity and the use of master events is required to determine accurate hypocenters and identify clusters. Finally, seismicity during refill is strongly clustered in space and time and is associated with some of the faults that were already seismically active during production. In Cornegliano Laudense, in the first three storage cycles, we have only observed 3 weak events with maximum magnitude M_L 1.0 close to the deposit, and no clustering on known or unknown faults.

For the surface deformation analysis, we first consider the GNSS signal. The residual signal from OL01, which corresponds to the behavior of the local station, is shown in Figure 6. Residuals were obtained by subtracting the trend from each reference station time series (calculated in the ITRF system). The relative displacement of station OL01 (purple dots in the foreground) was determined by subtracting the average trend of the four reference stations (black dots) from the absolute measurement. To highlight the movement of OL01 with respect to the reference frame, the mean value of OL01 calculated from January 2017 to the UGS start of operation was set to zero. It can be seen that the displacement of OL01 with respect to the assumed reference frame is consistent with those of the reference stations in the two measurement years before the start of UGS operation. Station OL01 departs from the trend of the other stations only in the second half of 2019, recording a horizontal shift to the NW and a vertical upward movement. Since it is not possible to estimate OL01 velocity directly from the time series ruling out the eventual effect of the storage activities, we cannot exclude that there is a limited effect due to the error in the estimation of the trend characterizing the time series before the start of the UGS activities.

[IMAGE OMITTED. SEE PDF]

We estimate the velocity of OL01 from the start of UGS operation with respect to the European Reference Frame EURA and we compare it with the velocities of reference stations. The velocity of OL01 is: 1.2 mm/yr on the north component (between 0.2 mm/yr and 0.6 for reference stations), 0.0 mm/yr on the east component (between 0.4 and 0.7 mm/yr for reference stations) and 8.1 mm/yr for the vertical component (between 0.2 and 0.5 mm/yr for the reference stations). Therefore OL01 moves in north-west and up direction with respect to the reference stations after the UGS start of operation. This deformation trend seems to be related to the inflation of the UGS due to gas injection.

The DInSAR analysis allows us to obtain dense maps of surface deformation velocities and the corresponding time series of displacement for the entire area, as shown in Figure 4. In this figure the comparisons between the deformation time series of three GNSS stations and the corresponding DInSAR measurements are also shown. It is worth noting that it is unfortunately not possible to compare the GNSS and the DInSAR deformation time series at the OL01 station because it is located in an area without coherent electromagnetic signals, which means that it is not possible to have sufficiently close DInSAR measurements points to allow a reliable comparison as in the case of the other GNSS stations considered in Figure 4.

In particular, Figures 7a and 7b show the zoom of the mean deformation velocity maps for the vertical and the east-west components of the surface deformation in the area around the UGS indicated by the white box in Figure 4. The time series of displacement (units in cm) of some neighboring points are also shown in the boxes below the map. The results presented show that there are no significant deformations in the areas of Cornegliano Laudense marked by points P1 and P3 during the period from March 2015 to December 2021. In contrast, the area corresponding to point P2 shows a clear deformation pattern with a significant vertical component, exhibiting uplift behavior since the beginning of 2019 and characterized by an average deformation rate that varies between 0.5 and 1 cm per year. The same panel of P2 also shows the total volume of gas contained in the natural reservoir during UGS operations (red line), normalized to the maximum volume value. The correlation between the two trends is easily seen. Consistent with the uplift phenomenon, the area near the UGS also exhibits a deformation trend in the east-west direction, as shown by the time series of displacement of points P2 and P4 in the middle panel of Figure 7b, which shows a total displacement of about 1 cm to the west in the area of point P2 and of about 1.5 cm to the east in the area of point P4.

[IMAGE OMITTED. SEE PDF]

Figure 4, and the corresponding zoom of Figure 7, also show a negative vertical displacement trend (i.e., subsidence) with a deformation rate of about −1 cm/yr in the area of point P4, which is located in the Turano Lodigiano area, about 15 km from the UGS. The vertical deformation rate at this point is about 1 cm/yr. Such a displacement is independent of the storage activities, as already shown by the results of the background analysis (see Figure 3a).

The observed behavior of surface deformation reflects the pattern of the gas volume contained in the natural reservoir during UGS operations (Figure 5c) and is consistent with other studies conducted at various sites with storage activities (Jiang et al., 2020; Qiao et al., 2018; Teatini et al., 2011; Zhou et al., 2019). In fact, the ground displacement initially follows the reservoir loading phase and then oscillates around a uniformly increasing mean value according to the periodicity of the process. Moreover, in Lanari et al. (2020), several long-term deformation time series relevant to other gas storage sites in Northern Italy are shown, which are retrieved identifying their periodic behavior with respect to time, by computing the correlation coefficient of each pixel's detrended time series with a sinusoid with a 1-year period. The surface deformation phenomena due to gas injection/storage operations are easily distinguishable from those relevant to water withdrawal activities because of the different periodicity of the oscillations for the time window.

Even though the surface deformation above the UGS is not uniformly sampled by DInSAR, it exhibits a low spatial gradient. Similarly to the seismicity case, the observed deformation is consistent with some results of the 3D geomechanical modeling performed during the design phase (Ital Gas Storage, pers. comm.), which predicts that: (a) the volume behaves predominantly elastically during the injection/withdrawal cycles both outside and inside the reservoir; (b) the deformation is reversible; and (c) the ground uplift calculated during the injection phases is non-uniformly distributed.

This work helps to address the lack of studies in the scientific literature on seismicity due to UGS in depleted gas reservoirs. It confirms that the operations of gas storage facilities in such reservoirs can result in negligible, if any, seismicity and limited surface deformation; more importantly, it emphasizes the importance of measuring undisturbed conditions prior to industrial activities for a sufficiently long period of time, using instruments with adequate sensitivity and nearly uniform in time. Without a dedicated seismic network deployed for 2 years prior to storage activities, we would have classified the surface events that occurred near the UGS site after fill began as “induced” events.

There is a general consensus in the Netherlands that induced seismicity may occur during the early storage cycles required to produce cushion gas and working gas; in any case, the magnitudes are likely to be far below those observed during reservoir depletion (Muntendam-Bos et al., 2022). Local seismic networks or deep borehole seismometers improve the detection of micro or nano events; for example, in the Bergemeer UGS, where monitoring is done with a 6-level borehole geophone string, reported events range from M_L −1.5 to 0.7. Although such a high detection limit is of great scientific value for understanding the dynamics within the reservoir, it goes far beyond what is required to ensure safety. In Italy, no induced events have been reported for the 15 operating UGS facilities, although we cannot say with certainty whether or not microseismicity associated with UGS facility operations has been properly recorded, since seismic monitoring data are publicly available for only some of these facilities (Augliera, 2019; Priolo et al., 2005; Romano et al., 2019).

The few cases of induced seismicity reported for this type of UGS, where the human perception threshold was reached, concern two large reservoirs located in seismically active regions: first, Gazli, where the actual attribution of seismicity to gas storage activities remains questionable (Suckale, 2010), and second, Hutubi, where the results of the analysis in favor of storage activities as the cause of the events that occurred near the UGS are consistent with physics-based modeling, although the observation of seismicity before the UGS, made with the same seismic network, is missing to exclude beyond doubt the possibility of a natural origin for the same earthquakes. Another case is reported by Silverii et al. (2021), but it is a different type of UGS, namely in a deep aquifer. These have some relevant differences from UGS in depleted reservoirs, in terms of their effective tightness when gas injection significantly increases the existing natural pressure. In addition, the study by Silverii et al. (2021) examined the seismicity detected when the UGS had already been in operation for 3 years, and the seismic baseline before the start of gas storage was not estimated. In all cases, a seismic baseline estimated for a reasonable period of time (e.g., two or 3 years) prior to the start of storage activities with the same monitoring infrastructure used during UGS operations would have provided important additional information.

Depleted oil fields behave very differently, and their exploitation as UGS requires much more care, especially at the outset; the case of the Castor UGS represents well the dynamics of the dramatic consequences that can occur. Because the oil is much denser and more viscous than the gas, there is a serious possibility that the Castor reservoir was not able to secure the gas sealing within the original oil field (Salò-Salgado, 2016). Cesca et al. (2021) believe that the sequence resulted from progressive fault failure and unlocking, with seismicity occurring mainly on a secondary fault beneath the deposit and triggered first by pore pressure diffusion and then by stress transfer. This interpretation is disputed by Vilarrasa et al. (2022), who argue that the locations used were not accurate and believe that seismicity in extensional stress zones such as the Valencia Trough was more likely triggered by vertical stress changes caused by buoyancy. In any case, the microseismicity resulting from this migration was not detected in time due to the poor monitoring system (Cesca et al., 2014; Vilarrasa et al., 2022). All these factors, together with the slow decision-making and the “prudence” in interrupting ongoing activities for business reasons, let us think that the observed seismicity was recognized too late as being caused by reservoir activities.

A final comment on our monitoring experience concerns the importance of transparency of data and information to the public, local governments and agencies, and the scientific community. In addition to providing services and scientific support for safer management of the UGS, high-quality seismic monitoring provides detailed data and, where possible, clear scientific answers to the questions and doubts of various stakeholders. The integrated monitoring project of the Cornegliano Laudense UGS was presented to the public in the presence of local administrations and, above all, in a meeting organized by the Prefetto di Lodi—the Prefetto is the authority representing the Italian government at the provincial level and is responsible for public order and safety—at the end of 2019, which was attended by all local authorities and representatives of public interests. We have experienced two opposite situations: on the one hand, very reasonable points of view and requests that are easy to answer, and on the other hand, misbeliefs and confusion that require long-lasting and intense scientific communication campaigns to confront and clarify them. Both situations should be met with transparency, and the different actors (public administrations, energy companies, scientific institutions, etc.) should respect the established rules and roles. The presentation of the entire monitoring system already in operation and the direct access to various information through the website—in this context, the significant contribution of the concessionaire to the creation of the website is appreciated—was a key element in reassuring the public that the activities are carried out safely and according to the highest standards.

Conclusions

We described the new monitoring infrastructure of the Cornegliano Stoccaggio UGS, which integrates seismic and deformation measurements at the ground surface.

In addition, we presented both the data collected and the analyses performed so far for the integrated monitoring, namely:

• the detailed picture of seismicity and ground surface deformation for a period of several decades before the start of gas injection— the so-called “background” analysis;

• the “baseline” for seismicity and ground deformation, recorded by the new monitoring network for 2 years before the operational phase of gas storage;

• the “current” measurements of seismicity and ground surface deformation for 3 years (i.e., for the period 2019–2021) after the start of UGS operation.

This study shows that the local seismicity is very weak and has not changed after the commissioning of the UGS, even considering the limits of the sparsely populated earthquake catalog. During the studied period, some deep tectonic earthquakes occurred, similar to the past seismicity that was initially interpreted as anthropogenic, and later classified as natural (Caviaga earthquakes in 1951).

Thus, our study shows that the first cycles of UGS operation, commonly considered particularly critical, at the Cornegliano Laudense site do not change the main characteristics of the previously recorded seismicity from a statistical point of view. We recall that the injection cycles were characterized by a gradual increase of the total gas volume in the reservoir, with the maximum operating pressure not exceeding the original confining pressure of about 160 bar.

The new seismic network and ground deformation monitoring greatly improve the detection of seismicity and ground deformation throughout the area surrounding the Cornegliano Laudense UGS. As a result of UGS activity, the approximately 5 km × 2 km area over the vicinity of the UGS site, which includes the town of Cornegliano Laudense, exhibits a uniformly distributed deformation pattern consistent with a slight uplift typical of inflation beginning in 2019, characterized by an average deformation rate and displacement in the range of 0.5–1 cm/year and 1–1.5 cm, respectively, and an extremely low spatial gradient.

Our observations are consistent with the results of the 3D geomechanical modeling conducted outside of this study (Ital Gas Storage, pers. comm.), which predicted no fractures in the caprock or within the reservoir during the injection/withdrawal cycles, nor reactivation of existing faults, and an elastic behavior of the volume during the injection/withdrawal cycles with reversible deformation. However, as the surface deformation above the UGS is not uniformly sampled by DInSAR, we propose to establish additional measurement points either with persistent backscattering points or low-cost GPS instruments.

Some general conclusions can be drawn from this study that go beyond the specific case studied here. The first is that gas storage in depleted gas reservoirs can result in negligible seismicity if carefully managed, such as not exceeding the initial reservoir pressure or using slow ramp-up injection phases, as suggested by Juanes et al. (2017). This is an objective fact that results from high-quality monitoring. Second, it is important to measure undisturbed conditions for a few years before industrial activities—the so-called baseline analysis—in order to correctly interpret the phenomena observed afterward. While the Italian guidelines recommend a seismic baseline established using the final “operational” configuration of the network with a monitoring period of at least 1 year, prestational and not only prescriptive rules for monitoring would be desirable to capture human-induced perturbations. For example, a sufficiently long baseline analysis for seismicity should be that which allows statistically meaningful activity rates to be obtained, such as the values of the a and b coefficients of the Gutenberg-Richter relationship. Therefore, a global risk assessment is needed to determine what financial exposure would be acceptable due to the delay of the project in order to obtain reliable baselines. Third, the proposed monitoring approach is fully in line with the Italian guidelines and underlines the effectiveness of these recommendations.

Acknowledgments

The integrated monitoring of the Cornegliano Laudense UGS is carried out jointly by OGS and IREA-CNR on behalf of Ital Gas Storage S.p.A. The authors thank this company for its full cooperation in fulfilling all operational and scientific tasks. Any opinions, findings, conclusions, or recommendations expressed in this paper are those of the authors and do not necessarily reflect the point of view of the concession holder. The phase/magnitude readings of the earthquakes cited in the text with an event identifier code (EVENTID) are published in the Online Bulletin of the International Seismological Centre (), and can be accessed through the following link: . RingServer software produced by IRIS—Incorporated Research Institutions for Seismology () is used to acquire data from the seismic stations. The BRTT Antelope system (; last accessed February 2022) is used to manage and process seismological data. The graphics were created using: ArcGIS from ESRI (), with contributions from OpenStreetMap 2021, for Figure 1; Python (pandas, matplotlib, and numpy), GMT (, Wessel et al., 2013), for Figures 2 and 6 and Figure S7 in Supporting Information S1; Matlab v. R2018b Update 4 (9.5.0.1067069) The MathWorks Inc., Natick, Massachusetts, for Figure 5; GoogleEarth for Figure 6a; Figures 3, 4, and 7 were created using a GIS software developed by IREA, with contributions from GoogleMaps and IDL 8.8.0, Harris Geospatial Solutions, Inc. 2020.

Data Availability Statement

All data used for this study are freely available and their distribution is licensed under a Creative Commons Attribution-ShareAlike 2.0 Generic License. Specific links on integrated seismic and ground deformation monitoring of the Cornegliano Laudense UGS can be found on the RSCL website () under the “Data and Documentation” section. Specifically for the different types of data:

Seismic data—Data used to locate events in this study (i.e., velocity models, the list of additional stations used to locate each event, and phase arrival times) can be found in Text S6 in Supporting Information S1. Continuous seismic waveforms and station information for the OL network are available through the OGS Archive System of Instrumental Seismology (OASIS, ); a link is also available on the RSCL website under the “Data and Documentation/Seismic Data/Access to Seismic Data” section. Continuous seismic data for the IV and GU networks were obtained from the EIDA database (). Seismic data were then processed using the BRTT Antelope system () and HYPOELLIPSE () to detect seismic events and locate earthquakes, respectively.

GNSS data—GNSS station OL01 (named LODI to adequate it to international repository, like ) data through the service provided by the OGS FReDNet geodetic network at or via anonymous ftp at . A link is also available on the RSCL website under the “Data and Documentation/GNSS Data” section. GNSS data from the SPIN network, that is, the CREA, CREM, MILA, MIL2, and PAVI stations, are also available at through the SPIN3 GNSS service. The GNSS time series shown in Figure 6b were created using the files in the Res folder and computed from the time series in the ITRF08 folder using the procedure described in Text S1.2 in Supporting Information S1. The data used to generate the solid red line are the same as in Figure 5. The GNSS time series shown in Figure S7 in Supporting Information S1 are taken from the files in the EURA and ITRF08 folders.

DInSAR data—The ERS-1/2 and ENVISAT SAR images used in this paper can be downloaded from the on-line ESA catalogs available at (for ERS-1/2) and (for ENVISAT), respectively.

The Sentinel-1 SAR images are freely available on-line and they can be downloaded from the on-line Copernicus Open Access Hub () or from the Alaska Satellite Facility (ASF)—Distributed Active Archive Center ().

The SAR images were processed by using the P-SBAS processing chain (Casu et al., 2014; Manunta et al., 2019) developed at IREA-CNR and available as a processing service on the ESA Geohazards Exploitation Platform (GEP) () and through the European Plate Observing System (EPOS) portal .

The DInSAR results shown in the paper are relevant to the generated time series (and the related mean deformation velocity estimates), which are made available by the authors in the repository . A special link is also on the RSCL website under the heading “Data and Documentation/ERS -1/2 and EVINSAT Data.” This repository also contains a document describing the format of the provided DInSAR products.

References

American Petroleum Institute (API). (2015). API recommended practice 1171—Functional integrity of natural gas storage in depleted hydrocarbon reservoirs and aquifer reservoirs.

Augliera, P. (2019). Concessione di stoccaggio di gas naturale “Minerbio stoccaggio” (BO), Struttura Preposta al Monitoraggio, Relazione Finale (p. [eLocator: 125]). INGV‐Istituto Nazionale di Geofisica e Vulcanologia. Retrieved from http://cms.ingv.it/sperimentazioni/minerbio

Barnaba, C., Laurenzano, G., Moratto, L., Sugan, M., Vuan, A., Priolo, E., et al. (2014). Strong‐motion observations from the OGS temporary seismic network during the 2012 Emilia sequence in northern Italy. Bulletin of Earthquake Engineering, 12(5), 2165–2178. [DOI: https://dx.doi.org/10.1007/s10518-014-9610-4]