Influence of hysteretic stress path behavior on seal integrity during gas storage operation in a depleted reservoir

Pierre Jeanne, Yingqi Zhang, Jonny Rutqvist

Lawrence Berkeley National Laboratory, Energy Geosciences Division, Berkeley, CA, 94720, USA

ABSTRACT In this study, we numerically investigate the influence of hysteretic stress path behavior on the seal integrity during underground gas storage operations in a depleted reservoir.Our study area is the Honor Rancho Underground Storage Facility in Los Angeles County (California, USA), which was converted into an underground gas storage facility in 1975 after 20 years of oil and gas production. In our simulations,the geomechanical behavior of the sand reservoir is modeled using two models:(1)a linear elastic model(non-hysteretic stress path) that does not take into consideration irreversible deformation, and (2) a plastic cap mechanical model which considers changes in rock elastic properties due to irreversible deformations caused by plastic reservoir compaction (hysteretic stress path). It shows that the irreversible compaction of the geological layer over geologic time and during the reservoir depletion can have important consequences on stress tensor orientation and magnitude. Ignoring depletion-induced irreversible compaction can lead to an over-estimation of the calculation of the maximum working reservoir pressure. Moreover, this irreversible compaction may bring the nearby faults closer to reactivation. However, regardless of the two models applied, the geomechanical analysis shows that for the estimated stress conditions applied in this study, the Honor Rancho Underground Storage Facility is being safely operated at pressures much below what would be required to compromise the seal integrity.? 2020 Institute of Rock and Soil Mechanics, Chinese Academy of Sciences. Production and hosting by Elsevier B.V. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/).

Keywords:Gas storage Stress path Geomechanical simulation Caprock integrity Fault stability Modified cam-clay model Honor rancho underground storage facility

1. Introduction

Although natural gas is used for a wide variety of purposes,it is mostly used to generate electricity and heat buildings, which results in a higher demand during the winter (e.g. CCST, 2018). In order to meet high demands for heating in the winter, natural gas can be stored during periods of low demand. The most common method of storing natural gas is underground storage in depleted oil and gas reservoirs.As an example,in USA,329 of the 415 natural gas storage facilities operating (79%) utilize depleted oil and gas fields(Ground Water Protection Council and Interstate Oil and Gas Compact Commission, 2017). Sites with such reservoirs involve an extensive network of distribution and gathering lines for injection and withdrawal operations and provide a large reservoir volume in the pore space with an effective caprock seal as demonstrated by the accumulation and production of hydrocarbons.

Two important characteristics of an underground storage reservoir are the amount of natural gas that can be stored and its capacity to hold it for future use. Indeed, the integrity of the reservoir or gas storage zone and associated seals are of paramount importance and should not be compromised.For this reason,there are limits to the pressure required to inject intended gas volumes,particularly at total inventory. This maximum operating pressure depends on the in situ state of stress that determines the fracture gradient(a function of in situ stress and rock tensile strength)and the stability of nearby faults (a function of stress, fault orientation and frictional strength). A starting point for maximum operating pressure analysis is a field’s discovery pressure(original formation pressure prior to hydrocarbon production and reservoir depletion).If an operator chooses to operate the facility above the discovery pressure, then additional reservoir and caprock analyses are required to be able to justify safe operations at that pressure level(Ground Water Protection Council and Interstate Oil and Gas Compact Commission, 2017). Even if the operator decides to keep reservoir pressure at or below discovery pressure, field evidence has suggested that caprock can still be compromised (Lynch et al.,2013). The reason, as pointed out by Santarelli et al. (1998, 2008),is that the stress state usually follows a different path during the production-induced reservoir depletion from the path during subsequent gas storage operations,and this hysteresis in the stress path may lead to erroneous calculation of the maximum reservoir pressure,which could potentially compromise the sealing capacity(Lynch et al., 2013). In general, sealing capacity can be compromised by fracturing of the caprock which could lead to subsurface leakage, or sealing capacity can be compromised by induced earthquakes (reactivation of existing faults) that could severely impact the caprock or well integrity.

Fig. 1. (a) 3D geomechanical model showing the different geological formations. (b) Layer forming the Wayside 13 sand formation divided into six domains to consider the anisotropic distribution of the hydraulic properties within the turbidities deposit and faults bounding the reservoir.

In this paper, we studied the geomechanical effects associated with natural gas storage in a depleted reservoir. The study area is the Honor Rancho Underground Storage Facility in Los Angeles County(California,USA).From 1955 to 1975,23 wells were drilled for oil and gas production in this field. In 1975, the field was converted to gas storage and 38 wells were completed to the storage zone(Wayside 13 sand)(Southern California Gas Company, 2009).The Honor Rancho Underground Storage Facility is the subject of a project led by Lawrence Berkeley National Laboratory to develop a risk management system for underground gas storage infrastructure. This risk management system relies on process models to understand physical processes occurring at the storage site.In this context, we developed a three-dimensional (3D) geomechanical model of the Honor Rancho gas storage field to evaluate stress conditions during 20 years of production-induced reservoir depletion followed by 40 years of gas storage operations. Specifically,we focus on the potential impact on the caprock integrity and on the stability of the reservoir bounding faults,considering impact of irreversible geomechanical behavior during the initial reservoir depletion and subsequent pressure cycling. To investigate the impact of irreversible behavior, the reservoir is modeled using a linear elastic model as a reference case for non-hysteretic stress path(ignoring irreversible behavior),and a plastic cap mechanical model for hysteretic stress path(considering irreversible behavior).

2. Three-dimensional geomechanical model

The structure of the Honor Rancho Underground Storage Facility is a southwest-dipping monocline(Schlaefer,1978).The gas storage reservoir is within the Wayside 13 sand,a coarse-grained turbidite deposit, and corresponds to the basal unit in the upper Miocene and lower Pliocene Towsley formation. The reservoir is capped by the shale beds of the Towsley formation.The formations above the caprock are the Pico formation which consists of marine sandstone and conglomerate (Schlaefer, 1978) and the Saugus formation composed of coarse to fine grained, unconsolidated to loosely consolidated sand and silt (Schlaefer,1978). The upper part of the Pico formation is called the Yule sand and is utilized for wastewater disposal. We used the surface elevations of these different rock formations(SoCalGas,2016)to build our 3D geomechanical model(Fig.1a).

The 3D geomechanical model extends to a depth of 6.3 km,and laterally ～10 km in the eastern and ～11 km in the northern direction(Fig.1a).The model is composed of 53,130 cells,with 4830 cells for the Wayside 13 sand formation. We have represented the Wayside 13 sand formation with a constant thickness of 170 m and with two layers(2415 cells per layer).The reservoir is bounded on the north by a normal fault oriented N090°-to-N100°dipping 80°N and on the south by a reverse fault oriented N090°dipping 58°N.These faults do not intersect the geological formations above the caprock. On the east and west, the reservoir boundaries are represented by a change in the hydraulic properties (lower porosity and lower permeability). Also on the east, the San Gabriel Fault is present. It is a subvertical right-lateral strike slip fault oriented～N140°(Schlaefer, 1978) (Fig. 1b). In our simulations,these faults are 30 m wide and are represented with different hydraulic and mechanical properties (Table 1).

3. Numerical approach

Numerical simulations were done using the iTOUGH-FLAC simulator (Blanco-Martin et al., 2016), which links the iTOUGH2(finite volume) multiphase flow and heat transport simulator(Finsterle,2004)and the FLAC3D(finite-difference)geomechanical code (Itasca, 2012) for coupled thermo-hydro-mechanical analysis under multiphase flow conditions. Simulations are carried out to investigate the influence of a hysteretic stress path behavior on the seal integrity during gas storage operation in a depleted reservoir.

3.1. Fluid flow simulation

Our primary goal is to estimate changes in stresses induced by changes in reservoir pressure without consideration of the nature of the produced fluids. Hence, for the sake of simplicity, we consider the reservoir formation to be fully water-saturated before reservoir depletion, and the hydrocarbon fluid production is simulated by water extraction. When the liquid is extracted from the pore space, it has to be replaced by something else (water, air,etc.). In our case, because we are interested in storage and production of methane and to simplify the calculation of multiphase and multicomponent fluid flow, the extracted liquid is directly replaced by methane. It should be noted that the quantity of methane created by this approach is very low as demonstrated by the very low reservoir pressure (～2.5 MPa) at the end of the reservoir depletion.To simulate this process,we used the EOS-CH4module (Finsterle, 2004) which provides the thermophysical properties of fluid mixtures (water and methane) needed for assembling the governing mass- and energy-balance equations in iTOUGH2.

We set the initial(pre-production)conditions considering linear pore pressure and temperature gradients(9.81 MPa/km and 22°C/km, respectively). This results in an initial pore pressure of～30 MPa and a temperature of ～66°C at reservoir depth(～2400 m depth below sea level). Laterally, we set no-flow boundaries (no mass or heat flow) except on the west where we impose an open flow boundary. The top of the model is of a constant atmospheric pressure. The simulations are conducted in an isothermal mode, such that the initial geothermal temperature gradient is maintained during simulation.

Porosity and permeability measurements made in the Honor Rancho Gas Storage Field show porosities ranging from 3.2% (in WEZU-C2A well) to 42.4% (in WEZU-5 well) and permeability ranging from 5 × 10-15m2(in WEZU-30 RD well) to 5 × 10-12m2(in WEZU-16A well).However,most of the porosity measurements range between 6% and 26%, and permeability is around 10-14m2(DOGGR,1991;GeoMechanics Technologies,2016;SoCalGas,2016).In our simulation, to consider this complex distribution of the hydraulic properties within the reservoir, we create six reservoir domains with various hydraulic properties (Fig. 1b). The size and shape of these reservoir domains are re-constructed based on maps of porosity and permeability provided by field operator(Davis and Kuncir, 2004). Hydraulic properties used for the different rock formations are summarized in Table 1.

To simulate 20 years of reservoir depletion from 1955 to 1975,we impose a stepwise declining bottom hole pressure for the 21 withdrawal wells (Fig. 2a) by imposing a progressive decrease in pressure and a progressive increase in methane saturation(Fig.2b).The frequency of the change is one year, and within a year, a constant pressure and a constant gas saturation are used.From 1955 to 1975,the reservoir pressure declines progressively from 30 MPa to 2.5 MPa in the 21 cells representing the 21 withdrawal wells in activity during this period. Then, to simulate the gas storage operation from 1977 to 2017 (no data are available from 1975 to 1977),we imposed bottom hole pressures calculated from monthly wellhead pressure data as reported in the California State Department of Conservation, Division of Oil, Gas, and Geothermal Resources (DOGGR) database (https://maps.conservation.ca.gov/doggr/wellfinder). To calculate the bottom hole pressure from the wellhead pressure, we used the following relation which was developed for dry gas wells (Lapeyrouse, 2002):

wherePbhis the bottom hole pressure,Pwhis the wellhead pressure,Sgis the specific gravity of gas,Ris a gas constant for the American Petroleum Institute(API)standard condition air(53.36(ft lb)/lb),His the true vertical depth of a well, andTavis the average temperature (in Rankin). In 2017, bottom hole pressures were updated every month in 36 wells. Fig. 2c shows the location of the 39 cells where pressure and gas saturation were imposed from 1977 to 2017(because some wells were abandoned and others drilled at the end of 2017, pressure and gas saturations were imposed in 37 wells).

When we impose the bottom-hole pressure in one cell,we also calculate the total flow rate through the cell. The reservoir hydraulic properties are calibrated by making sure the simulatedmonthly flow rate is roughly consistent with the “measured flow rate”(same order of magnitude)(Fig.3).The main difficulty is that the DOGGR database provides one value of wellhead pressure for the whole month,but the operator may not operate with a constant pressure.Whereas in our simulation,we calculate an average flow rate for each month by assuming the operation is the same each month.

Table 1Hydraulic and mechanical properties used during the geomechanical simulations.

Fig.2. Map of the reservoir domains with the locations of(a)the 21 production wells(white stars)used during reservoir depletion and(c)the 39 wells(red stars)used during gas storage operations.(b)Evolution of the imposed pressure at the 21 production wells during reservoir depletion(from 1955 to 1975)and(d)example showing the evolution of the imposed bottom hole pressure at three wells during gas storage operation (1977—2017). The blue, red and green triangles show the location of the three control points used to monitor the mechanical reservoir compaction simulated with the modified Cam-Clay model.

3.2. Geomechanical simulation

Geomechanical simulations were performed to analyze (1) the likelihood of a failure of the caprock by examining the critical pressures that could induce hydraulic fracturing (PFRAC) or induce slip along existing fractures critically oriented for shear reactivation(PSHEAR), and (2) the stability of the faults bounding the reservoir.The geomechanical behavior of the reservoir is modeled using both a linear elastic model which results in non-hysteretic stress and strain changes, and a plastic cap mechanical model, possibly leading to hysteretic stress and strain changes.However,during the initialization of the geomechanical model using a plastic cap mechanical model, irreversible plastic reservoir compaction occurs that changes the initial state of stress. Thus, an exact comparison between these two models is not possible. However, despite this difference, we still try to investigate how changes in elastic properties within the gas storage zone influence the seal formation integrity.

Fig. 3. Comparison between the simulated monthly flow rates and the measured monthly flow rates.

Modeling of reservoir compaction associated with pressure depletion is usually done with mechanical models originally developed for unconsolidated soil e.g. the Drucker-Prager yield criterion(Drucker and Prager,1952),and the Pelessone smooth cap model (Pelessone,1989; Foster et al., 2005; Motamedi and Foster,2015). Here, we adapt the modified Cam-Clay model as a plastic cap model for simulating the hysteretic stress path due to plastic void compaction. The modified Cam-Clay model has been successfully applied in petroleum geomechanics for modeling irreversible reservoir compaction behavior (Zoback, 2007). A more recent study shows that the compaction behavior of Permian Rotliegend sandstones from reservoirs in the Netherlands (at comparable depth as the reservoir in this study) measured in triaxial experiments can be described by modified Cam-Clay plasticity(Pijnenburg et al., 2019). It should be also noted that reservoir deformations are not monitored during reservoir depletion and during most of the gas storage operation,thus we cannot compare the deformation associated with these two phases,and also we cannot discuss which of the two models (the linear elastic model or the plastic cap mechanical model)is better to describe the mechanical behavior of the Wayside 13 sand.

Studies made by Townend and Zoback (2004) and Wilde and Stock (1997) suggest that the stress regime at Honor Rancho is a thrust faulting stress regime, withSv≤Shmin≤SHmax, whereSvis the vertical stress,Shminis the minimum horizontal stress, andSHmaxis the maximum horizontal stress. Wilde and Stock (1997)analyzed 14 wells in the East Ventura Basin west of the Honor Rancho Gas Storage Field.The dominant breakout azimuth is to the northwest with an average direction of N44°W, implying a minimum compressive stress direction of N44°W, and a maximum compressive stress direction of N46°E. Therefore, based on these results, we simulate a case of a thrust faulting stress regime withShmin = 1.1SvandSHmax= 1.2Svand withSHmaxorientated N046°(Wilde and Stock,1997).The stresses follow the lithostatic gradient and we assumeSv= ρgh, where ρ is the rock wet density (here assumed to be ρ=2600 kg/m3),gis the acceleration of gravity andhis the depth.Roller boundaries conditions(no displacement in the normal direction to boundary)are imposed on the side boundaries whereas the ground surface is free to move and null displacement is imposed at the bottom of the model.Mechanical properties used for the different rock formations are summarized in Table 1.

3.2.1. Sealing integrity

PFRACis estimated by assuming that a hydraulic fracture could develop as soon as the fluid pressure exceeds the minimum compressive principal stress, σ3(Jaeger et al., 2007), thus, with compressive stress positive,we have

A lower limitPSHEARis estimated by assuming that a fault of arbitrary orientation could exist anywhere in the caprock.For such a case,the Coulomb failure criterion can be written in the following form (Jaeger et al., 2007):

whereS0and φ are the coefficient of internal cohesion and angle of internal friction, respectively; τmand σmare the two-dimensional maximum shear and mean stresses in the σ1-σ3plane, respectively, defined as

Here,S0is considered negligible because the fault cohesive strength contributes little compared to the compressive field stresses at a depth of several kilometers(Zoback et al.,2003)and φ is considered equal to 36.86°(Goodman,1989). Under these conditions and according to Eq. (3), we have

We also use the Coulomb failure criterion to estimate the changes in proximity to failure (Δσc) for the faults bounding the reservoir:

where Δτ is the change in shear stress on the plane in the expected rake (slip) direction on the target fault, Δσnis the change in normal stress (compression is positive), μ is the coefficient of friction, and ΔPis the change in pore fluid pressure (Harris and Simpson, 1992). Positive values of Δσcindicate that the fault plane in question is closer to failure, whereas negative values indicate stress changes that move the stress state on the fault farther from failure.

3.2.2. The modified Cam-Clay model

The modified Cam-Clay model (Roscoe and Burland, 1968)represents a particular form of nonlinear elasticity and hardening/softening behavior governed by volumetric plastic strain.The model originally is developed and applied on clay assuming that when the porous formation is slowly compressed under isotropic stress conditions and under perfectly drained conditions, voids will compact and cause irreversible deformation changing the rock elastic properties when the failure envelope is reached.

The failure envelope is an ellipse in the von Mises equivalent stress (q) versus mean effective stress (p) plane. The flow rule is associated and the yield function is given by

whereMis the slope of the critical state lineq=Mp;andpcis the consolidation pressure defined as the mean effective stress at which yielding occurs in hydrostatic conditions (q= 0). Thus,pccontrols the size of the yield surface. Upon failure,pccan increase(hardening and associated plastic compaction),decrease(softening and associated plastic dilation) or remain constant (no change in plastic volumetric deformation). Therefore, the failure surface is defined by a series of ellipses in the (p, ν,q) space, where ν is the specific volume, defined as

whereVis the total volume,Vsis the volume of the non-porous(solid) space, φ is the porosity, andeis the void ratio (ν, φ andeare dimensionless).

Specific volume and volumetric strain (ε) are related according to

Changes in consolidation pressure are given by

where εplis the incremental plastic volumetric strain; andkand λ are the material parameters corresponding, respectively, to the slopes of the elastic unloading-reloading (swelling) lines and the normal consolidation line, as shown in Fig. 4 (Itasca, 2012; Roscoe and Burland, 1968). As the normal consolidation pressure,p, increases,the specific volume,ν,of the material decreases.The point representing the state of the material moves along the normal consolidation line defined by

Fig. 4. Normal consolidation line and unloading-reloading (swelling) line for an isotropic compression test (Itasca, 2012).

where νkis the specific volume corresponding top1on the consolidation line.

An unloading-reloading excursion, from pointAorBin Fig. 4,will move the point along an elastic swelling line of slopek,back to the normal consolidation line where the path will resume. The equation of the swelling lines has the following form:

where νλis the specific volume corresponding top1on the swelling line.

Finally, given the nonlinear relationship betweenpand ν, the tangent bulk modulus is not constant, but evolves following Eq.(14). Then a new shear modulus (G) is calculated from Eq. (15),considering a constant Poisson’s ratio,ν, which is equal to 0.25:

The reservoir at Honor Rancho is a southwest-dipping monocline,meaning that the initial effective stress strongly varies within the reservoir. The initial mechanical parameters are established to allow compaction to occur during reservoir depletion both in the shallowest and deepest parts of the reservoir while keeping the calculated porosity (Eq. (9)) close to that measured on rock samples. To reach these objectives, we use a sloping normal consolidation line,λ=0.1,a pre-consolidation pressure,pc0=5×107Pa,a reference pressure,p1 = 4 × 105Pa, and a specific volume at reference pressure,p1,on the normal consolidation line,νλ=1.3,1.5,1.6,1.7,1.8 and 2.1 for reservoir domains 1 to 6,respectively.We also use the slope of elastic swelling line,k=0.003(Eq.(14))to calculate realistic bulk modulus for porous sandstone and a frictional constant for sandstoneM= 1.4 (Rutter and Glover, 2012). Because of these different parameters and the variation in effective stress across the reservoir,the initial compaction is not evenly distributed within the porous sandstone. As a result, within the reservoir, the bulk modulus varies from 4 GPa to 11.7 GPa(with an average value of 6.6 GPa and a standard deviation of 1.2 GPa), and the shear modulus varies from 2.3 GPa to 7 GPa (with an average value of 3.9 GPa and a standard deviation of 0.7 GPa).These values are used as input parameters for the linear elastic model.

4. Simulation results

In our simulations, the geomechanical behavior of the sand reservoir is modeled using(1)a linear elastic model(meaning that there is no plastic deformation, no changes in porosity, and no changes in mechanical properties which result in a non-hysteretic stress patch from reservoir depletion to the gas storage operation);and(2)a plastic cap mechanical model to consider changes in rock elastic properties due to irreversible deformations caused by plastic reservoir compaction (hysteretic stress patch). In the following,we present(1)how the irreversible deformations caused by plastic reservoir compaction simulated with the modified Cam-Clay model influence the reservoir mechanical properties (Section 4.1), and (Section 2) the comparison between the linear elastic model and the plastic cap mechanical model with the modified Cam-Clay model to highlight the influence of reservoir compaction on caprock integrity (Section 4.2) and the stability of bounding faults (Section 4.3).

4.1. Mechanical reservoir compaction simulated with the modified Cam-Clay model

Figure 5 shows the simulation results of the mechanical behavior of the reservoir formation during oil production (from 1955 to 1975)and the gas storage operations(from 1975 to 2017)at three control points located at depths of -2957.4 m, -2817 m,and -2756 m (blue, red and green curves, respectively, in Fig. 5)obtained with the modified Cam-Clay model.The locations of these three control points are shown in Fig. 2 (blue, red and green triangles corresponding to the blue, red and green curves, respectively). The pressure drop associated with reservoir depletion causes an increase in effective stress(Fig.5a),resulting in reservoir compaction (highlighted by the negative values of the total volumetric strain in Fig. 5b), and by the decrease in porosity (highlighted by the decrease in specific volume(Eq.(9))in Fig.5c).This compaction also causes an increase in bulk (Fig. 5g) and shear moduli (Fig. 5d). In the simulation, an increase in mean effective stress from 24 MPa to 36 MPa results in an increase in bulk modulus from 9 GPa to 14 GPa (blue curve in Fig. 5a and g). This is in good agreement with laboratory tests made on porous sandstone as reported in the literature (Li and Fj?r, 2012).

It can also be observed that both during the reservoir depletion and the gas storage operations, the mean effective stress reaches the consolidation pressure several times (Fig. 5e), causing irreversible porosity compaction(highlighted by the increase in plastic volumetric strain in Fig. 5f).

4.2. Influence of reservoir compaction on caprock integrity

Figures 6 and 7 show the distribution and evolution ofPSHEARandPFRACin the caprock (just above the reservoir) before oil production(1955),after 20 years of oil production and 43 years of gas storage operations(2017),as well as the changes between 1955 and 2017. These results are shown both for the case of a reservoir mechanical behavior modeled with the linear elastic model and with the modified Cam-Clay model.

Common to all the simulations, the simulation results show:

(1) The lowest values ofPFRACandPSHEARare calculated towards the northeast in the shallowest part of the caprock and the highest values are found towards the southwest in the deepest part of the caprock (Fig. 6a—d and 7a-d,respectively);

(2) A drop inPFRAC(here the vertical stress) (Fig. 6g and h) andPSHEAR(Fig. 7g and h) during the reservoir depletion;

Fig.5. Evolution at three control points within the reservoir of(a)mean effective stress,(b)accumulated total volumetric strain,(c)change in specific volume,(d)shear modulus,(d) pre-consolidation pressure, (f) accumulated plastic volumetric strain, (g) bulk modulus, and (h) mean deviatoric stress.

Fig. 6. Distribution in the caprock of PFRAC (PFRAC in 2017 minus PFRAC in 1955) (a, b) before reservoir depletion and (c, d) after 20 years of production and 43 years of gas storage operations(2017)calculated with the modified Cam-Clay model and with the linear elastic model.(e,f)Distribution in the caprock of ΔPFRAC for the two models.(g,h)Evolution of ΔPFRAC at two control points (located in Fig. 6e) calculated with the modified Cam-Clay model and with the linear elastic model, respectively.

(3) An increase inPFRACandPSHEARtoward their original values in the area where the gas injection and withdrawal wells are located (Figs. 6h and 7h). This area also corresponds to the zone of the highest porosity and permeability in our reservoir simulation (Fig. 2); and

(4) No or very small stress recovery in areas located farther from the injection and withdrawal wells(Figs.6g and 7g)in parts of the reservoir with lower porosity and permeability(Fig.2).

The pressure drop (from ～30 MPa to 2.5 MPa) and associated reservoir compaction cause these changes in stress. The compaction causes the overburden to subside,but vertical support from the sideburden (the rock just outside the permeable storage reservoir)prevents the overburden from subsiding at the same rate as the reservoir compaction causing a slight increase(from 0.4 MPa to 2 MPa in Figs. 8a and 10a) in the maximum horizontal stress above the reservoir formations. As a result, the vertical distance from the surface to the reservoir increases. This stretching of the overburden is responsible for the decrease in vertical stress (up to -8 MPa in Fig. 6) that is leading to the calculated decrease inPFRAC. In the case of the reservoir mechanical behavior modeled with the modified Cam-Clay model, the calculated reservoir compaction during reservoir depletion is larger due to the plastic volumetric compaction strain. This additional compaction increases the stretching of the overburden and results in a higher vertical stress reduction in the caprock(Figs.6g and 7g).However,it can be noted that locally we observe the opposite behavior with a smaller decrease inPFRAC(here the vertical stress) during the reservoir depletion in case of the modified Cam-Clay model(Figs. 6h and 7h). It is because at this location, the initial bulk modulus is slightly higher than that in the elastic simulation.

Fig.7. Distribution in the caprock of PSHEAR and ΔPSHEAR(PFRAC in 2017 minus PFRAC in 1955)before reservoir depletion and after 20 years of production and 43 years of gas storage operations(2018)calculated with the modified Cam-Clay model(a,c,e)and with the linear elastic model(b,d,f).(g,h)Evolution of ΔPSHEAR at two control points(located in Fig.6e)calculated with the modified Cam-Clay model and with the linear elastic model, respectively.

4.3. Influence of reservoir compaction on fault stability

The model simulation shows that the reservoir exploitation and the irreversible compaction have resulted in complex changes in the stress tensor below,within and above the reservoir.To illustrate these changes and to understand how they influence the stability of the surrounding faults,we present the variation in the stress tensor at locations far (～2.5 km from the center of the injection and withdraw zone, Fig. 8), near (～0.9 km from the center of the injection and withdraw zone, Fig. 9) and very close (～0.3 km from the center of the injection and withdraw zone,Fig.10)to the center of the production zone from 1955 to 2017. At these locations,stresses are monitored ～50 m above,within and ～50 m below the reservoir.

Fig. 8. Variations in (a) σ1,σ2 and σ3 and in (b) their orientations relative to the east (angle α), to the north (angle β) and to the horizontal (angle θ); and (c) the variations in the effective normal stress the shear stress(Δτ)and Δσc at three control points located along the fault N140—90N above,within and below the reservoir,for the case of the linear elastic model (dark blue, blue and light blues curves) and for the modified Cam-Clay model (red, orange and yellow curves).

However, before looking at the stress tensor evolution and its influence on the nearby faults, we present in Table 2 the initial σcvalues calculated by Coulomb failure criterion(Eq.(7))for both the linear elastic model and the modified Cam-Clay model. These σcvalues are calculated at the beginning of the simulation (before depletion) at three control points along the faults N140—90N,N090—58N and N090—80N and located above, within and below the reservoir. The σcvalues are given for a reservoir mechanical behavior simulated with the linear elastic model and the modified Cam-Clay model.It can be observed that these values are very high due to the initial stress tensor(Shmin=1.1SvandSHmax=1.2Sv)that result in a very low deviatoric stress and thus a very low probability of shear failure(high σcvalues).Nevertheless,we can also observe that the lower values of σc(in green in Table 2)calculated at control points located below the reservoir and close to the withdrawal zone are obtained for the modified Cam-Clay model along the fault N090—58N(south boundary)and along the fault N090—80N(north boundary).At the other control points,faults are closer to failure for the linear elastic model. The reason these differences occur between the two models is that the initial reservoir compaction occurring during the initial reservoir depletion phase causes the stresses to be redistributed.

Fig. 9. Variations in (a) σ1,σ2 and σ3 and in (b) their orientations relative to the east (angle α), to the north (angle β) and to the horizontal (angle θ); and (c) the variations in the effective normal stress the shear stress(Δτ)and Δσc at three control points located along the fault N090—58N above,within and below the reservoir,for the case of the linear elastic model (dark blue, blue and light blues curves) and for the modified Cam-Clay model (red, orange and yellow curves).

Figures 8—10 show, at control points along faults N140—90N,N090—58N and N090—80N as well as at locations above,within and below the reservoir, the variations in (a) the three principal stress components (Δσ1, Δσ2, Δσ3); (b) their orientations relative to the east(angle α),to the north(angle β)and to the horizontal(angle θ);and(c)the variations in the effective normal stressthe shear stress(Δτ)and the changes in proximity to failureResults are shown for the case of the linear elastic model(dark blue,blue and light blues curves) and for the modified Cam-Clay model (red, orange and yellow curves).

4.3.1. Far from the withdrawal zone: fault N140—90N

At the beginning of the simulation,the orientations of the stress tensor are different when the reservoir mechanical behavior is simulated with the linear elastic model or with the modified Cam-Clay model (Fig. 8b). This difference is caused by the initial compaction within the reservoir during the initialization of the numerical model. However, despite these initial differences, the changes in the stress tensors are quite similar.

Fig.10. Variations in(a) σ1,σ2 and σ3 and in(b) their orientations relative to the east (angle α), to the north(angle β) and to the horizontal (angle θ);and (c) the variations in the effective normal stress the shear stress(Δτ)and Δσc at three control points located along the fault N090—80N above,within and below the reservoir,for the case of the linear elastic model (dark blue, blue and light blues curves) and for the modified Cam-Clay model (red, orange and yellow curves).

Table 2Initial σc values calculated by Coulomb failure criterion for three control points above, within and below the reservoir along the faults N140—90N, N090—80N and N090—58N for both the linear elastic model and the modified Cam-Clay model.

Far from the withdrawal zone, these changes in the stress tensors are small.The reservoir depletion causes a stress drop of a few MPa (except a small increased in σ1above the reservoir up to～0.3 MPa)and a small rotation of the stress tensor of a few degrees(Fig.8a).The σ2rotates towards the west(increase in angle β)and σ1towards the north(decrease in angle α)(it is a counterclockwise rotation) (Fig. 8b).

Within and below the reservoir, these small changes in the stress tensor have little effect on the fault stability.Theincreases up to ～2 MPa during the depletion but returns to its original value during gas storage operation,and τ slightly increases and the fault becomes slightly closer to instability (increase in Δσc). Above the reservoir, stresses stay lower than their initial magnitudes and do not rotate back to their original orientation due to the drainage of the caprock during reservoir depletion.Because of this behavior,stays high and thus the fault stays farther to instability during the gas storage operation (decrease in Δσc).

4.3.2. Near the production zone:fault N090—58N(south boundary)

As observed before,the reservoir depletion causes a stress drop(up to 10 MPa,Fig.9a)and a stress rotation:σ1rotates towards the north and σ2towards the west (～4°—10°counterclockwise rotation,Fig.9b).In case of the linear elastic model,these rotations are associated with a rotation of σ3from vertical to the horizontal plane up to ～20°.

In case of the modified Cam-Clay model,the initial compaction calculated during the initialization of the numerical model has caused an initial rotation of the stress tensor.Indeed,in this case,σ3is initially tilted of 20°,40°and 10°from the vertical above,within and below the reservoir (θ = 70°,θ= 50°,θ= 80°, respectively, in Fig.9b).The consequence is that during the reservoir depletion,the stretching of the overburden causes the vertical stress to decrease and σ3to rotate back towards the vertical axis (increase in θ).

During the reservoir depletion, these small changes in stress tensor and pore pressure cause an increase inand a decrease in τ that brings the fault farther from instability(decrease in Δσc).Then,during the gas storage operation,slightly decreases above the reservoir, whereas it reaches its initial value within and below the reservoir. The τ decreases at the three control points both during the reservoir depletion and during the gas storage operation.However, in case of the modified Cam-Clay model, because of the relatively low fault angle (58°) and the rotation of σ3towards the vertical axis, the decrease in τ is less important than the one simulated in the linear elastic model. These changes cause Δσcto(1) decrease during depletion and stay constant during the gas storage operation above the reservoir, and (2) decrease during depletion and become slightly higher or very close to its original value during gas storage operation within and below the reservoir.

In general, the calculated variations in the Coulomb criterion during the simulation with the modified Cam-Clay model keep the fault closer to instability(either increase in Δσcor smaller decrease in Δσc).

4.3.3. Very close to the withdrawal zone: fault N090—80N (north boundary)

In the case of the linear elastic model and during the reservoir depletion, the stress variations are very similar to what observed previously:a stress depletion(up to ～10 MPa,Fig.10a)and a ～15°counterclockwise stress rotation above, within and below the reservoir (Fig. 10b). These stress changes cause an increase inand a decrease in τ that brings the fault farther from failure(decrease in Δσc) (Fig. 10c). During the gas storage operation,returns to its initial value whereas τ increases and reaches values above the initial one. These changes cause Δσcto (1) decrease during depletion and stay constant during gas storage operation above the reservoir,and(2)decrease during depletion and become slightly higher or very close to its original value during gas storage operation within and below the reservoir.

In case of the modified Cam-Clay model, the influence of the initial reservoir compaction calculated during the initialization of the numerical model is very important. The stress redistribution associated with this compaction has caused a rotation of σ3from vertical towards the horizontal with an angle of θ= 18°above the reservoir,θ=51°within the reservoir and up to θ=73°below the reservoir causing a change from a thrust stress regime to a strikeslip stress regime (green rectangle in Fig.10b).

This stress variation causes τ to decrease during reservoir depletion and to increase during gas storage operation.As a result,the fault stress state is closer to failure during reservoir depletion than it is in the linear elastic model. We observe the opposite during the gas storage operation with the fault stress state farther to failure compared to the case of the linear elastic model.

5. Conclusions

In this study, we numerically investigate the influence of a hysteretic stress path behavior on the seal integrity during gas storage operations in a depleted reservoir. Our study area is the Honor Rancho Underground Storage Facility in the Los Angeles County (California), which is converted into an underground gas storage facility after 20 years of oil and gas production.

It is well known that during the production-induced reservoir depletion, as pore pressure dropped, the reservoir compacted and overburden subsided, but support from the sideburden prevented it to subside at the same rate as the reservoir compaction.This leads to the formation of a stress arch, with a decrease of the vertical stress and an increase of the horizontal stress above the reservoir and the opposite in the pillars of the stress arch.Then it is assumed that during the gas storage operation, the pore pressure increase should return the stresses back to their initial stress state before hydrocarbon production. This is why if an operator chooses to operate the facility above the discovery pressure(which is not the case at the Honor Rancho Underground Storage Facility), then additional reservoir and caprock analyses are required to be able to justify safe operations at that pressure level.

Here we show that irreversible deformations caused by plastic reservoir compaction can have important consequences for the calculation of the maximum working reservoir pressure on the seal integrity. The calculation of the maximum reservoir pressure to avoid the creation of a hydrofracture or the reactivation of fractures favorably oriented for shear reactivation in the seal formation tends to be overestimated when irreversible deformations are ignored(Fig.6c versus Fig.6d).In our case,considering a hysteric stress path behavior, the calculated changes inPSHEARandPFRACat the end of the reservoir depletion are few MPa higher than the calculated changes for a non-hysteretic stress patch.

These differences in the simulations results using the two models are caused by two phenomena:

(1) An initial compaction of the turbidite deposit. This initial compaction may reflect the compaction occurring during the compaction of the geological layers within a heterogeneous geological layer (as in a turbidite deposit) made of softer zones embedded in a more rigid zone.This phenomenon can cause locally a strong stress rotation and even a change in stress regime that may influence the stability of the nearby faults.However,it is very likely that stresses equilibrate over time due to creep, yielding and other slow time-dependent processes resulting in a more uniform stress distribution within the heterogeneous formation.

(2) An irreversible compaction of the reservoir during reservoir depletion that increases the stretching of the overburden causes the vertical stress to decrease and the stress tensor to rotate.

Nevertheless, regardless of the two models applied, the geomechanical analysis shows that for the best estimated stress conditions applied in this study, the Honor Rancho Underground Storage Facility is being safely operated at reservoir pressure much below what would be required to compromise seal integrity.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

The authors are deeply grateful to their project partner Southern California Gas(SoCalGas)company for their support.This work was conducted with funding provided by the California Energy Commission under the contract PIR-16-027 for Research on Risk Management Framework for Underground Natural Gas infrastructure in California. We are also grateful to the Special Issue Guest Editor Antonio Pio Rinaldi and to the two anonymous reviewers of this paper.