About geomechanical safety for UGS activities in faulted reservoirs

A critical issue concerning geomechanical safety for UGS (underground gas storage) in compartmentalized reservoirs is fault reactivation. Indeed, the displacement (land subsidence, land upheaval) and the stress fields caused by the seasonal injection and production of CH4 into and from deep reservoirs is peculiar. The need of improving our understanding of compartmentalized reservoir behavior and to define safe bounds for the pressure fluctuation in order to prevent undesired land movements and induced seismicity is becoming even more important. This also in view of the expected energy transition when large amount of green energy will potentially be stored and recovered through UGS of compressed air or hydrogen. In this framework, an indepth modelling investigation has been carried out for the typical UGS geological setting and operations in The Netherlands. The specific goals of the study are the following: (i) explaining the possible mechanisms responsible for seismic events unexpectedly recorded during UGS phases; (ii) understanding which are the critical factors (e.g. the geological configuration, the geomechanical properties, and the reservoir operations) that increase the probability of fault reactivation during the various UGS stages; and (iii) advancing possible guidelines for safe UGS operations. This contribution summarizes the main outcomes obtained by the modelling simulations: the combinations of factors causing fault reactivation during primary production (PP) are also more prone to generate fault failure during cushion gas injection (CG) and UGS. In fact, fault activation during PP leads to a stress redistribution and a new (deformed) “equilibrated” configuration that is newly loaded, in the opposite direction, when the pressure variation changes the sign because of CG and/or UGS. Finally, the various combinations have been ranked to highlight the conditions where the fault system is most likely reactivated during CG and UGS operations: the initial stress regime of the system, the geomechanical properties of the fault, and dislocation of the reservoir compartments are the major influencing drivers to fault instability.


Introduction
Induced seismicity is become a major issue in fluid production from and injection into deep formations. Apart from fracking, where low-permeability formations are intentionally fractured to increase productivity, thus (micro-) seismicity are an inevitable consequence of a successful project (e.g., Farahbod et al., 2015), several seismic events caused by "conventional" removal or injection of diverse fluids are listed in specific databases, such as HiQuake (https: //inducedearthquakes.org/, last access: 1 March 2020), and related publications (e.g., Foulger et al., 2018).
Seismicity induced by fluid removal from faulted formations is somehow "expected" when the extracted volumes or the associated pressure decline overcome a certain bound (Fig. 1a). McGarr (1991) suggested that net extraction of oil and water reduced the average density of the upper crust, thus decreasing the normal stress that prevents slip of faults. In compartmentalized gas reservoirs, seismicity has been reported when large pressure declines cause significant differential compaction between produced and unproduced blocks, i.e. large increase of the shear stress on the sealing faults between the compartments (van Eck et al., 2006). Nicholson Figure 1. Sketch of "expected" and "unexpected" seismic events. "Expected" seismicity can occur during (a) hydrocarbon primary production (PP) and (b) during fluid injection when the pore pressure depletion or increase, respectively, overcomes a certain threshold. (c) "Unexpected" seismicity can occur during cushion gas injection (CG) and underground gas storage/injection (UGS) UGS when the pore pressure is comprised between the value in initial undisturbed condition and the minimum at the end of PP.
and Wesson (1992) suggested that an earthquake might occur in response to larger stresses imposed by fluids migrating into the mid-to-lower crust. The change in pressure resulting from withdrawal of oil induces fluid migration within nonsealing faults that brought the faults closer to failure because of the reduced effective stress normal to the discontinuity surfaces and the deterioration of the fault mechanical properties (i.e., a decrease of the friction angle).
According to the Mohr-Coulomb criterion, also the introduction of fluids into faulted reservoirs is expected to encourage fault failure (Fig. 1b). A number of cases of seismic activity has been associated to geologic CO 2 sequestration (e.g., Verdon et al., 2013) and wastewater disposal (e.g., Ake et al., 2005). The explanation of these seismic occurrences has been provided by Vilarrasa et al. (2019): the overpressurized fluid flows within the faults, decreasing the nor-mal stress acting on the discontinuities that are consequently re-activated.
Induced seismicity cannot be explained by "simply" invoking the Mohr-Coulomb criterion when an event develops in correspondence to a fluid pressure distribution already experienced by the geologic system. This is the typical condition of UGS reservoirs where the pore pressure fluctuates seasonally between a maximum P max and a minimum P min values due to CH 4 injection during summer and withdrawal during in winter (Fig. 1c). A certain seismic activity has been recorded in a few faulted UGS fields in The Netherlands, toward the end of the primary production (PP), at the end of the cushion gas (CG) injection (when the pressure returns approximately to the natural value because a volume of gas is injected in an underground storage reservoir to provide the necessary pressure to deliver working gas volumes to customers), or during injection of gas within the UGS cycles (TNO, 2015;NAM, 2016).
In this framework, the main goal of the study is (i) to understand which the mechanisms are responsible for this "unexpected" seismicity, and (ii) consequently provide some considerations about the possibility of safe UGS activities, i.e. seasonally storing gas into and producing gas from UGS fields reducing the risk of re-activating the faults bounding and/or crossing the reservoir.

Conceptual model
A set of numerical simulations has been carried out to investigate the processes of interest. The geological setting, geometric features, and production history typical of UGS reservoirs in The Netherlands have been used (e.g., Fokker et al., 2016;Wassing et al., 2017).
The reservoir is made of two adjacent square blocks, 2000×2000 m wide, confined laterally by four faults, namely F1, F2, F4 and F5 (Fig. 2). Another fault, F3, subdivides the reservoir in two compartments, so that the pore pressure changes P 1 and P 2 in the two compartments can differ. The reservoir is embedded in a 30 km wide square domain. Faults F4 and F5 are vertical, F1 and F2 are inclined with a dip angle equal to ±10 • . The dip angle of fault F3 can vary. The reservoir is 200 m-thick and 2000 m deep. The bottom of the model is 5000 m deep and the land surface has an elevation of 0 m. The faults extend from −3000 to −1500 m depth, i.e. they terminate within the caprock sealing the reservoir (Zeichestein formation). Block 2 can be offset in the vertical direction of 100 and 200 m, corresponding to half the thickness and the entire thickness of the reservoir (Fig. 2).

Set-up of the numerical model
A 3D hexahedral finite element (FE) -interface elements (IE) mesh is developed (Fig. 2). The mesh consists of Table 1. Scenarios addressed in the sensitivity analysis. M 1 and M 2 are the ratios between the minimum and maximum horizontal principal component of the natural effective stress regime and the vertical principal component, respectively.
the intra-field fault F3 is characterized by a dip = 65 • (instead of being vertical) 3b the intra-field fault F3 is characterized by a dip = −65 • (instead of being vertical) 3c offset in the vertical direction = 100 m between the two reservoir blocks (instead of zero offset) 3d offset in the vertical direction = 200 m between the two reservoir blocks (instead of zero offset) 4a the natural stress regime in the horizontal plane is rotated by 90 • 4b reservoir stiffness E = 20 GPa (instead of 11 GPa) 7a uneven pressure changes in the two reservoir blocks during UGS: P 1 = −100 bar and P 2 = 0 bar (instead of P 1 = P 2 = −100 bar) 7b uneven pressure changes in the two reservoir blocks during UGS: P 1 = −100 bar and P 2 = −200 bar (instead of P 1 = P 2 = −100 bar) 8 pressure changes during UGS P 1 = P 2 = −150 bar (instead of −100 bar) 9 viscous caprock (instead of elastic) 253 165 nodes and 236 208 FEs with a finer discretization in the reservoir layers, i.e., at depth between 1800 and 2200 m. The element size within the reservoir is 100 × 100 × 20 m. IEs, which are zero-thickness FEs composed of two surfaces that can slide or move apart one from each other when a failure criterion is overcome, have been integrated into the 3D FE mesh. Figure 2b shows the fault system embedded in the continuous 3D grid. An overall number of 5215 IEs is used to discretize the five faults. Stress and strain fields associated to UGS have been simulated using the M3E_GEPS3D (Geomechanical visco-Elasto-Plastic Simulator -3D) simulator, an in-house developed code (Spiezia et al., 2017;Isotton et al., 2019). The code simulates the possible activation of pre-existing faults with a quasi-static approach. The discontinuity surfaces are modelled according to the principles of contact mechanics as inner boundaries embedded in the continuous body. IEs are implemented for fault discretization (Franceschini et al., 2016. The rupture activation is governed by the Mohr-Coulomb failure criterion (Labuz and Zang, 2012), i.e. a fault reactivation occurs whenever the shear stress exceeds the limiting value τ L : with σ n and τ the normal and shear stresses, respectively, acting on the fissure surfaces, φ and c the fault friction angle and cohesion, respectively. The principle of "Maximum Plastic Dissipation" is used to define the direction of the limiting shear stress. The pore pressure variation within the faults is the average between the values experienced by the two adjacent portions separated by the discontinuity. Standard conditions with zero displacements on the outer and bottom boundaries are prescribed, and the land surface is a no-stress boundary.

Simulated scenarios
The model has been initially applied to the so-called reference scenario, i.e. a scenario based on the typical characteristics of the UGS reservoirs in the Netherlands. The two reservoir blocks experience a pressure change P amounting to −200 bar over a 10-year PP phase. PP is followed by a 2-year CG injection phase when the pressure recovers to the initial (undisturbed) value P i and then UGS cycles characterized by a 6-month production phase, during which P = −100 bar, and a 6-month injection period when the pressure returns to P i . For a detailed description of the reference scenario, see Teatini et al. (2019).
Then, a sensitivity analysis has been developed to understand the geometric features, geomechanical parameters, and pressure distributions that make the fault system likely to be reactivated during CG and UGS. Firstly, each parameter characterizing the system has been varied at a time. In a second stage, simultaneous modifications of various parameters of the reference set-up have been addressed. The various scenarios addressed by the sensitivity analysis are summarized in Table 1. Notice that only likely configurations have been tested. The aim of the study is not to analyse "extreme" conditions.

Outcome of the numerical model
The model results provide the 3D displacement and stress fields on each nodes of the 3D grid. In this study, the main interest is focused on the IE solution and, specifically, on the stress conditions and eventually sliding of the discontinuity surfaces.
Two parameters are used to quantify the fault state in relation to their possible reactivation. Firstly, the criticality index χ = τ/τ L , which represents the ratio between the actual tangential stress τ acting on a fault and the limit tangential stress τ L according with the Mohr-Coulomb criterion. When χ = 1.0 the fault starts sliding. The second one is t 80 , i.e. the fault thickness with a criticality index χ > 0.8. An example of |τ | and χ distribution at year 10 for the scenario 3d is provided in Fig. 3. Consistent with previous modelling outcomes (e.g., Wassing et al., 2017), fault reactivation develops mainly at the top and bottom of the reservoir formation and propagates over a short distance into the caprock.

Processes responsible for "unexpected" seismicity
The modelling outcome reveals that the scenarios causing fault reactivation during PP are more prone to fault failure during CG injection and UGS. Indeed, reversing the sign of the pressure change when only elastic deformations developed (i.e., no fault has been activated) causes the system to be unloaded, returning to the original stress regime. Conversely, fault activation during primary production (generally at the end of this development phase in the scenarios addressed in this study) leads to a stress redistribution and a new (deformed) "equilibrated" configuration that is newly loaded, in the opposite direction, when the pressure increases due to CG injection.
The faults approach a critical stress state during CG and UGS generally at the end of the injection/production phases when the cumulative pressure change -both decreasing and increasing -reaches the largest value. Inspection of the results allows pointing out that the initial stress regime can play a major role: a significant decrease of the horizontal principal components, as tested in scenario 4b, favours an early fault reactivation. Other main factors yielding the faults close to failure are a reduced friction angle (scenarios 5b and 5c), a large dislocation between producing compartments (scenario 3d), a significant difference between the reservoir stiffness and that characterizing the caprock, sideburden, and underburden (scenario 6a), and an uneven pressure change within adjacent compartments (scenarios 7a and 7b).

Ranking the conditions prone to "unexpected" seismicity
Post-processing of the modelling outcomes has been aimed at defining a methodology to rank the mechanisms, geological settings, geomechanical and production parameters in relation to their potentiality of inducing "unexpected" seismic events. The sensitivity scenarios have been ranked following these nested criteria: 1. (a) χ max during UGS; or (b) χ * = element (t e ·χ | χ>0.8 normalized over max (χ * ) | scenarios where t e is the element area. χ * is indicative of the fault extent where activation is close to occur; 2. maximum value of average sliding (δ avg ), evaluated on active elements only; 3. loading step of activation, i.e. the minimum pressure change causing fault activation.
The analysis has been carried out for the central fault F3 and fault F2 that is representative of the discontinuities bounding the reservoir. Notice that fault F3 is inactive in some of the investigated scenarios because of symmetry in the geological setting and driving forces. The results of the two ranking procedures are reported in Tables 2 and 3 for fault F2, and in Tables 4 and 5 for fault F3. It is interesting to note that the ranks obtained with the two criteria are similar which means that, when 0.8τ L < τ < τ L , fault reactivating is very likely to occur. As expected, the critical factors influencing fault activation during CG and UGS cycles arranges differently for the boundary and central faults. On fault F2, the stability is mainly jeopardized by the initial stress regime of the system, the geomechanical properties of the system (e.g., reservoir stiffness) and the characteristics of the faults (cohesion, static friction angle, presence of fault weakening). Due to symmetry conditions of fault F3 for most of the scenarios, the major influencing drivers to fault instability are given by geometrical parameters characterizing the fault/reservoir system (Tables 4 and 5). Dislocation of the reservoir compartments, non-vertical fault plane, different pressure changes in the two compartments, and fault dip are the features threating the stability of the fault F3.
It is worth noting the results of the scenario with a viscous caprock (scenario 9). Conversely to results obtained by previous modelling studies (e.g., Wassing et al., 2017), Tables 2 and 3 reveal that the presence of a salt viscous formation sealing the top of the reservoir is not a key parameter in favoring fault reactivation. This difference with previous outcomes, which were focused on production reservoirs, is likely due to the short-term (seasonal) fluctuation of the pressure change in the underlying UGS reservoir that limits the full development of the viscosity effect. However, it must be specified that the caprock discretization used in this 3D mod- elling study (20 m along the vertical direction) is one order of magnitude larger than those used in 2D investigations specifically devoted to the caprock behavior. The relatively large mesh could smooth the viscous effects.

Conclusions
The numerous scenarios investigated within the study have clearly revealed that fault reactivation can occur during CG (cushion gas) injection and UGS (underground gas storage), and is more likely if seismicity has been recorded during PP (primary production). Although the results are qualitative because of the theoretical/general framework of the modelling applications, preliminary "suggestions" can be sketched to avoid potential induced seismicity during CG and UGS.
A limitation on the maximum pressure P max (below the initial pore pressure) at the end of CG and UGS injection should be prescribed depending on the geometry of the fault/reservoir: the presence of sloped faults, dislocation of the reservoir compartments, differential pore pressure between adjacent reservoir compartments and within each reservoir block are critical factors to be accounted for to define a safe P max bound. A large difference between the reservoir and caprock stiffness is also a criticality factor that can threaten the fault stability. In these conditions, limitations on the operational pressure fluctuation during a UGS cycle and on the rate of pressure recovery during CG injection should be properly prescribed depending on the specific reservoir features.
Specific investigations are also fundamental to characterize the reservoir setting and thereby to avoid or reduce the risk during operations. Importance should be given to improve the knowledge on the following aspects (from the most to the least importance): the initial stress regime of the faulting system; the reservoir stiffness; the geomechanical parameters of the faults (failure criterion); the reservoir permeability (or, equivalently, the pore pressure distribution within the reservoir compartment).
Data availability. Data are available by writing to the corresponding author.
Author contributions. All the authors contributed to design the modelling study, develop the numerical analyses, and write the paper.
Competing interests. The authors declare that they have no conflict of interest.

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