Building 1D and 3D Mechanical Earth Models for Underground Gas Storage—A Case Study from the Molasse Basin, Southern Germany

Abstract

Hydromechanical models of gas storage in porous media provide valuable information for various applications ranging from the prediction of ground surface displacements to the determination of maximum reservoir pressure and storage capacity to maintain fault stability and caprock integrity. A workflow to set up such models is presented and applied to a former gas field in southern Germany for which transformation to a gas storage site is considered. The workflow comprises 1D mechanical earth modeling (1D MEM) to calculate elastic properties as well as a first estimate for the vertical and horizontal stresses at well locations by using log data. This information is then used to populate a 3D finite element model (3D MEM) which has been built from seismic data and comprises not only the reservoir but the entire overburden up to the earth’s surface as well as part of the underburden. The size of this model is 30 × 24 × 5 km³. The pore pressure field has been derived from dynamic fluid flow simulation through history matching for the production and subsequent shut-in phase. The validated model is ready to be used for analyzing new wells for future field development and testing arbitrary injection-production schedules, among others.

1. Introduction

A thorough understanding of the pre-production stress state and its changes with time plays an important role for the safe and economic operation of hydrocarbon reservoirs and underground storage sites. The subsurface stresses affect various operational aspects such as wellbore stability, caprock integrity, stress induced material parameter changes, fault reactivation, reservoir compaction and subsidence as well as stimulation techniques such as hydraulic fracturing [1,2,3]. Some regional information about at least one of the components of the stress tensor, i.e., the orientation of the maximum horizontal stress, can be gathered from the World Stress Map (WSM) [4]. Magnitude information concerning the three principal stresses is generally sparse. In addition, on the scale of a reservoir both stress magnitudes and orientations can be quite variable due to lithological heterogeneities and structural complexities [5]. Thus, a tool is required which can provide a robust prognosis of the in situ stresses in space and time, i.e., the complete stress tensor and its changes during subsurface operations, honoring the specific conditions of the hydrocarbon reservoir and storage site, respectively. Numerical modeling and, in particular, mechanical earth modeling (MEM) has proven to be a very valuable tool to integrate various data sets and to investigate the hydromechanical response of a reservoir under different operational conditions throughout its life cycle [3,6,7,8,9].

The data base for such an MEM approach is typically derived from a wide range of geological, geophysical and engineering data including field measurements, core tests, well logs, drilling, and production data. Once the numerical model has been verified by calibration data, such as in situ stress measurements and historical production data, it can be utilized for testing future operational scenarios such as, for example, stability of new boreholes, optimal orientation of horizontal well trajectories for multiple hydraulic fracturing, and subsidence due to pore pressure reduction. With respect to underground storage, MEM can be used, among others, to assess the maximum safe storage capacity and maximum threshold pressure avoiding fault reactivation and maintain caprock integrity [8,10].

In this study, we present a worked example on how to set up and populate a 3D hydromechanically coupled model of an underground gas storage site. This case study is typical for the conversion of depleted reservoirs into underground gas storage reservoirs. In such cases the data base required for a MEM is rather old, i.e., from the exploration phase of the reservoir. Nevertheless, even in such cases the workflow outlined below allows us to generate a numerical model for testing of future injection-production schedules and assessment of maximum storage capacity. The validated coupled 3D geomechanical model provides the stress state at any location and along well path in order to plan a new well. Furthermore, the model is an open platform for short term (weekly) as well as long term (seasonal) gas storage scenario testing in future. This type of study is unique in Molasse Basin in South Germany and can be applied on underground gas storage, CO₂ sequestration, and geothermal projects in the area or anywhere in the world.

2. Methodology

The entire modeling workflow is illustrated in Figure 1. It starts with 1D modeling of individual wells to calculate mechanical properties from log data. The 3D model comprising the reservoir along with underburden, overburden, and sideburden formations is built using a structural/geological model which is typically based on seismic and borehole data. Upscaling of hydromechanical properties from the various 1D MEM’s to the 3D MEM can be achieved through geostatistical methods such as Gaussian and Kriging interpolation method [11]. The pore pressure field in the model domain is calculated by a fluid flow simulation which uses the same structural/geological model as the 3D MEM. Mechanical and hydraulic calculations are coupled via pore pressures and effective stresses leading to deformation which in turn affects porosity and permeability. Both the 3D MEM and the fluid flow model need to be calibrated by comparing calculation results to observational data, e.g., in situ stress measurements and historical production data. Once both approaches have been validated, the coupled model can be used to study potential future injection-production schedules such as the seasonal cycles typical for underground gas storage or more frequent cycles to cover short-term variations in demand and supply.

Modeling outlines below utilize the Schlumberger software suite, i.e., Petrel for geological modeling, Techlog for 1D MEM, Eclipse for fluid flow modeling and Petrel RG (Visage) for 3D MEM. Geomechanical modules of Techlog and Petrel RG provide the user interface for geomechanical analysis and various elastic properties calculation methods. Eclipse is the standard industry software to provide pore pressure evolution of the reservoir throughout its operating history. All these Schlumberger software suites can be coupled to perform hydromechanical simulations.

2.1. 1D Mechanical Earth Modeling (1D MEM)

1D mechanical earth modeling (1D MEM) is the numerical representation of rock mechanical properties and the state of in situ stresses along a borehole [12]. Therefore, well log data are used to obtain various material properties and estimates for pore pressure and in situ stress along a wellbore section assuming a poro-elastic approach. The quality of a 1D MEM model depends on the availability of well log data and the availability of laboratory tests on cores for calibration. The complete 1D MEM modeling scheme (step by step) is presented in Figure 2.

2.1.1. Vertical (Overburden) Stress

The 1D MEM process starts with the overburden stress calculation based on density information. If a density log is readily available, the vertical stress or overburden stress can be determined according to:

$$S_{\mathsf{\nu }}=g{{\displaystyle \int }}_0^{TVD}{\rho }_{b\left(z\right)dz}$$

(1)

where $S_{\mathsf{\nu }}$ is the vertical stress, ${\rho }_{b\left(z\right)}$ is the bulk density at depth z below the earth’s surface, and g is the gravitational acceleration. For absent or incomplete density log data, which frequently are the case especially for old wells, data gaps must be filled first with synthetic densities which can be obtained by using several empirical correlations [13]. The extrapolation density method can be used in order to fill the air gap between bulk density logs in the computing process of overburden stress (Equation (2)). The combination of both extrapolated (synthetic density) together with bulk densities integrate the splice density to compute the overburden stress.

$${\rho }_{extrapolated}={\rho }_{mudline}+ A_0\ast {\left(TVD-AirGap-WaterDepth\right)}^{\alpha }$$

(2)

where ${\mathsf{\rho }}_{extrapolated}$ is the extrapolated synthetic density, $A_0$ and α are the fitting parameters, ${\rho }_{mudline}$ is the density at the sea floor (offshore) or ground level (onshore), and TVD is the true vertical depth.

2.1.2. Pore Pressure

The pore pressure is a key parameter required to determine effective stresses in 1D MEM. Various methods and correlations have been put forward to get a pore pressure profile (see [13]), but the most popular one providing reasonable results is Eaton’s method [14]. Either compressional slowness from sonic logs (Equation (3)) or resistivity from corresponding electric logs (Equation (4)) can be used to calculate pore pressure by using Eaton’s method.

$$P_p=S_v-\left(S_v-H_g\right) \ast \alpha \ast {(\frac{\Delta t_n}{\Delta t})}^n$$

(3)

$$P_p=\left(S_v-H_g\right){(\frac{R_{log}}{R_n})}^n$$

(4)

where, $P_p$ is the pore pressure, $S_v$ is the overburden gradient, $H_g$ is hydrostatic gradient, $\Delta t_n$ and $R_{log}$ are compressional slowness in shales at normal pressure and resistivity log, respectively, corresponding to the normal compaction trend, $n$ is the Eaton exponent 1.2 when using the resistivity log and 3.0 for the velocity log, and α is the Eaton factor [14].

2.1.3. Elastic Properties

Static elastic properties required for geomechanical modeling such as Young’s modulus, Poisson’s ratio, and shear and bulk moduli as well as the Biot coefficient have to be obtained from core analysis and laboratory measurements, respectively. MEM requires information for the entire overburden and at least part of the underburden section. A continuous set of at least dynamic elastic properties can be computed from specific logging data. Though these dynamic properties do not yet reflect the static properties actually required, and they can be transferred by lithology-dependent empirical correlations whose applicability can be crosschecked at calibration points (if lab tests on cores are available) [15].

The following equations can be used to derive dynamic elastic properties such as dynamic Poisson’s ratio and dynamic Young’s modulus from sonic log data (Equations (5) and (6)) [12]. The static elastic properties may then be calculated by using various empirical correlations, e.g., the Bradford correlation [16].

$${\mathsf{\nu }}_{\mathrm{d}} =\frac{0.5{\left(\Delta t_s/\Delta t_c\right)}^2-1}{{\left(\Delta t_s/\Delta t_c\right)}^2-1}$$

(5)

$${\mathrm{E}}_{\mathrm{d}} =\frac{{\mathsf{\rho }}_{\mathrm{b}}}{\Delta t_c{}^2}\frac{\left(1-2{\mathsf{\nu }}_d\right)\left(1+{\mathsf{\nu }}_d\right)}{\left(1-{\mathsf{\nu }}_d\right)}$$

(6)

In these equations, ${\mathsf{\nu }}_d$ is the dynamic Poisson’s ratio, E_d is the dynamic Young’s modulus in Mpsi, $\Delta t_c$ is compressional wave interval transit time in μs/ft, $\Delta t_s$ is shear wave interval transit time in μs/ft and ${\rho }_b$ is the bulk density in g/cm³, and $\Delta t_s$ is the shear wave slowness in μs/ft [13].

In the absence of shear wave slowness, i.e., if only the standard compressional sonic log is available, the Greenberg and Castagna correlation [17] (Equation (7)) can be used accordingly:

$$V_s=0.8042 V_p-855.9$$

(7)

where $V_s$ is the shear wave velocity and $V_p$ is the compressional wave velocity, i.e., the inverse of the compressional wave slowness.

The Biot coefficient can be calculated from the static bulk moduli of the dry skeletal frame and the rock material (Equation (8)) [18].

$$\alpha =1-\frac{K_{dry}}{K_{solid}}$$

(8)

where $\alpha$ is the Biot coefficient, $K_{solid}$ is the static bulk modulus of the rock material, and $K_{dry}$ is the static bulk modulus of the dry skeletal frame.

2.1.4. Strength Parameters

The indirect estimation of strength parameters is possible by using various correlation methods [19].

Bradford [16] has proposed a correlation of unconfined compressive strength (UCS) with static Young’s modulus (Equation (9)) for soft formation (Young’s modulus < 12 GPa) that can be used to calculate UCS.

$$\mathrm{UCS}=2.28+4.1089E_s$$

(9)

where, UCS and $E_s$ are the unconfined compressive strength and static Young’s modulus, respectively.

Tensile strength is considered to be 1/10th of the UCS, if the Bradford correlation [16] has been used to compute UCS from static Young’s modulus [20]. Friction angle (FANG) can be estimated by using Gamma ray log and porosity; it uses the Plumb’s correlation (1994) (Equation (10)).

$$FANG=26.5-37.4\left(1-\Phi -V_{clay}\right) + 62.1{\left(1-\Phi -V_{clay}\right)}^2$$

(10)

where FANG is friction angle, $\Phi$ is porosity, and $V_{clay}$ is volume of clay. It uses gamma ray log to estimate V_clay. In the absence of porosity log, compressional slowness (DTCO) can be used to compute the porosity log by using Wyllie’s correlation [21] (Equation (11)).

$$\Phi =\frac{DTCO-47.8549}{21371.5414-47.8549}$$

(11)

2.1.5. Horizontal Stresses

The horizontal stresses can be estimated based on a poro-elastic approach (Equations (12) and (13)) [22].

$$S_{hmin}=\frac{\mathsf{\nu }}{1-\mathsf{\nu }}{S}_{\mathsf{\nu }}-\frac{\mathsf{\nu }}{1-\mathsf{\nu }} \mathsf{\alpha }{P}_p+ \mathsf{\alpha }{P}_p+\frac{E}{1-{\mathsf{\nu }}^2}{\epsilon }_h+\frac{\mathsf{\nu }E}{1-{\mathsf{\nu }}^2}{\epsilon }_h$$

(12)

$$S_{Hmax}=\frac{\mathsf{\nu }}{1-\mathsf{\nu }}{S}_v-\frac{\mathsf{\nu }}{1-\mathsf{\nu }} \mathsf{\alpha }{P}_p + \mathsf{\alpha }{P}_p + \frac{E}{1-{\mathsf{\nu }}^2}{\epsilon }_H + \frac{\mathsf{\nu }E}{1-{\mathsf{\nu }}^2}{\epsilon }_h$$

(13)

In the above equations, $S_{hmin}$ and $S_{Hmax}$ are minimum and maximum horizontal stress magnitudes, respectively, S_v is vertical stress, E is Young’s modulus, $\mathsf{\nu }$ is Poisson ratio, α is Biot coefficient, P_p is pore pressure, and ${\mathsf{\epsilon }}_{\mathrm{H}}$ and ${\mathsf{\epsilon }}_{\mathrm{h}}$ are the maximum and minimum horizontal strain magnitudes. Appropriate values for the two horizontal strains and, hence, the two horizontal stress magnitudes, can be found by fitting these values to observational data on in situ stress magnitudes, e.g., from extended leak-off tests and minifracs [22].

2.2. 3D Mechanical Earth Modeling (3D MEM)

3D MEM uses the workflow illustrated in Figure 3. The subsurface geometry consisting of faults and lithostratigraphic layers is adopted from the geological/structural model which in turn is usually based on seismic and well data. The 3D MEM is then populated with material properties calculated from the log-based 1D MEM approach described above. This upscaling of properties from the various 1D MEM’s of the individual wells to the 3D MEM can be achieved through geostatistical methods such as the Gaussian and Kriging interpolation method [11]. If available, information from seismic amplitude versus offset (AVO) analysis which captures lithological changes can be used to populate the 3D MEM [7]. The final 3D MEM is discretized, and finite element (FE) techniques are applied to calculate stresses and strains in the model domain. In an iterative process similar to the 1D MEM, the calculated stresses are compared to observed stress data now from wells of the entire model domain. If a satisfactory fit is achieved, this model is regarded as validated in mechanical terms and can be used for coupling with the fluid flow model.

2.3. Fluid Flow (Dynamic) Model

The flow simulation accounts for multiphase fluid flow in porous media. For the flow of two immiscible fluids in a porous medium, Darcy’s law and the corresponding mass balance can be written as [23],

$$v_{\pi i}=-\frac{K_{ij}{K}_{r\pi }}{{\mu }_{\mu }}\left(\frac{\partial }{\partial x_j}{P}_{\pi }+{\rho }_{\pi }gh\right)$$

(14)

$$-\frac{\partial }{\partial x_i}\left(\frac{1}{B_{\pi }}{v}_{\pi i}\right)=\frac{\partial }{{\partial }_t}\left(\frac{1}{B_{\pi }}\phi S_{\pi }\right) + q_{\pi }$$

(15)

where, $\pi$ is fluid phase (oil, gas, water), $v_{\pi i}$ is fluid velocity for phase π, $K_{ij}$ is permeability tensor, $K_{r\pi }$ is relative permeability for phase $\pi$, ${\mu }_{\mu }$ is viscosity for phase $\pi$, $P_{\pi }$ is fluid pressure for phase $\pi$, ${\rho }_{\pi }$ is density for phase $\pi$, g is gravitational acceleration, h is fluid height, $B_{\pi }$ is volume factor for phase $\pi$, $\phi$ is porosity, $S_{\pi }$ is saturation for phase $\pi$, and $q_{\pi }$ is fluid source or sink. The development of reservoir and fluid properties is modelled in a series of discrete steps through space and time using finite difference techniques [24]. Thereby, the behavior of the multiphase system can be described by complex pressure, volume, temperature (PVT) and special core analysis (SCAL) relations.

2.4. Coupled Hydromechanical (HM) Modeling

For processes such as storage of gas in a porous medium the flow simulation and the geomechanical simulation have to be coupled: the pore pressure controls the effective stresses and, hence, deformation which in turn changes rock porosity and permeability which again affects fluid flow. Thus, two sets of equations, one for fluid flow and one for equilibrium of forces have to be solved on the discretized model domain. There are different coupling approaches solving these sets of equations. These approaches could be one-way coupling, two-way coupling, and full coupling. The one-way coupling scheme transfers pore pressure from the flow simulator to the mechanical model at scheduled time steps to calculate stresses and rock deformation. No modeling results are passed back to the flow simulator. This scheme is usually sufficient for HM modeling of gas reservoirs as gas compressibility dominates the bulk rock compressibility and the mass balance is mainly controlled by gas pressure rather than by the stresses of solid rock [25]. In contrast, two-way coupling involves sending information on porosity and permeability changes due to rock deformation back to the flow simulator for updating the hydraulic properties there. This in turn affects fluid flow and the resulting pore pressure field which is transferred to the mechanical model. This loop is performed at each time step until convergence is reached. Full coupling or implicit coupling implicates that both simulators, i.e., the flow as well as geomechanical simulator calculate the coupling parameters by using a system of equations in which these coupling parameters would be taken as unknown assuring internal consistency [24,26]. One-way coupling has been used in this case study.

3. Case Study

The case study used to apply the workflow outlined above is a depleted gas field for which a potential conversion to an underground gas storage site is investigated. The data base was kindly provided by Uniper SE. The former gas reservoir is located about 65 km east of Munich/Germany in the northern foreland of the Alps, the so-called Molasse Basin (Figure 4a). It is located in a structural trap, i.e., an anticlinal structure bounded by a north-dipping normal fault (Figure 4a) at about 1800 m depth. The reservoir rock is of Late Oligocene age, i.e., the Chattian Hauptsand with three gas bearing layers having a total thickness of 85 m. The reservoir formation is highly porous with an average porosity of 23% and permeabilities in the range from 20 mD to 500 mD [27].

The field has been produced from 1958 till 1976 and has remained in a shut-in state ever since. The initial reservoir pressure at the time of discovery was 16 MPa and dropped to a minimum level of 3 MPa during the production phase. Due to inflow of water, the present-day pressure has increased to about 15.3 MPa again.

3.1. Geological Setting

As the 3D MEM does not only comprise the reservoir but has to incorporate the entire overburden up to the earth’s surface as well as part of the underburden, some information on the geological setting of the case study has to be gathered. The wedged-shaped fill of the Molasse basin consists of Cenozoic sediments that are underlain by Mesozoic strata (Triassic to Cretaceous limestones and sandstones) and crystalline basement rocks of the European plate (Figure 4b and Figure 5). The maximum thickness of the basin fill is about 6 km in the South and progressively decreases towards the North. The Molasse Basin proper developed from Late Eocene to Late Miocene times during the late stage of the Alpine orogeny. Major subsidence took place during Oligocene and Lower Miocene time which lead to rapid sedimentation in the basin [10]. The deposition of pelitic deposits occurred during Early Oligocene (Rupelian) due to shallow to deep marine settings. Then, a subsequent sea level fall resulted in the deposition of terrestrial to coastal facies which formed Aquitanian and Chattian sand during the Late Oligocene to Lower Miocene [28]. The sea level rise during Lower Miocene resumes the marine settings again in the Molasse basin, which resulted in the deposition of marine and tidal deposits (Burdigalian deposits). During Middle to Late Miocene, the terrestrial and fluvial environment replaced the marine settings leading to the formation of the Upper Freshwater Molasse [28]. The sedimentary rocks in the Molasse Basin provide a source and reservoir as well as seal rocks for hydrocarbon reservoirs [29]. However, only small oil and gas fields have been found in the German part of the Molasse basin. Most of these hydrocarbon fields have been depleted by now and some are considered for underground gas storage projects.

3.2. Data Set

3.2.1. Well Data

For eight wells (X1, X2, X3, X3a, X4, X5, X6, and X7) some log data are available (Figure 6). As most of the wells have been drilled before 1980, these data are incomplete and lack the modern logging available nowadays. However, this is the typical situation for old, depleted fields for which a conversion to storage sites is considered. In the present case, a density log is available only for well X6. It is used as a reference density log to calculate overburden and pore pressures for each of the wells. No gamma ray is available for any of the wells. Compressional sonic logs exist for seven wells; in addition the absent shear slowness for all seven wells was calculated from the empirical modeling described in Section 2.1.3.

3.2.2. Calibration Data

Some calibration data for the geomechanical as well as the fluid flow modeling are available. Laboratory tests such as triaxial tests were carried out on the core samples of X6 taken at reservoir level which lead to average values of 4.18 GPa for static Young’s modulus ($E_s$), 0.35 for static Poisson’s ratio (${\mathsf{\nu }}_s$) and 0.82 for the Biot coefficient α [27]. While the data can be used to calibrate the empirical correlations used to find mechanical properties in 1D MEM, stress magnitudes of 38.2 MPa and 36.7 MPa observed at reservoir level of well X6 for S_Hmax and S_hmin [33] are used to constrain the stresses calculated by the 3D MEM. For calibration of the fluid flow simulation, the evolution of the reservoir pressure at well X6 from 1969 to present is available.

3.2.3. Static Geological Model

The model geometry was provided by Uniper SE, which was built up using seismic and well data from which depth and thickness contour maps were generated for the entire model domain. The stratigraphic layers from the earth’s surface to the reservoir and even further down to the underburden formations have been extracted from drilling reports. The scheme of stratigraphic layers for reservoir (top and base of Chattian sandstone), overburden (Top Aquitansand, Bases Burdigal, Bases Helvet up to the earth surface) and underburden (Top Lithothamnienkalk, Top Cretaceous, Top Malm Zeta, Top Malm Gamma, Top basement down to the flat model base) formations is illustrated in Figure 7. The static reservoir model is shown in Figure 6. The model contains three faults with the main fault which cuts through the middle of the reservoir in the east-west direction. The volume of the reservoir model proper is about 8 × 4 × 0.3 km³ in X, Y and Z direction, respectively, and located at a depth of about 1200 m. The 3D grid is made up of 495,720 cells which are 162 × 85 × 36 cells in X Y and Z directions, respectively. The gridded part of the reservoir covers the area of 34 km² around the anticlinal structure of the reservoir. The reservoir geometry comprises an anticlinal structure with a north-dipping normal fault.

3.2.4. Dynamic Hydraulic Model

The reservoir model (Figure 6) was exported to the dynamic fluid flow simulator. The corresponding high-resolution reservoir model has 14 layers with individual petrophysical properties of each layer. The reservoir is highly heterogeneous with a porosity distribution of 0% to 30% and water saturation from 0 to 1. The porosity and permeability range is from 0% to 21% and 0 to 80 mD, respectively. The active cells cutoff is determined by a minimum pore volume of 200 m³.

4. Model Generation

4.1. 1D Model

The extrapolation density method (Equation (2)) is used to calculate the overburden stress (Equation (1)) in this case study. The pore pressure profile is acquired from the Eaton’s method (Equation (3)). Furthermore, the Plumb and Bradford [12] correlation is used to compute dynamic elastic properties and the same correlation is used to estimate static elastic properties as well as rock strength properties. These parameters are calibrated against the core material properties available [27]. The poro-elastic strain modeling approach is used to estimate horizontal stresses (Equations (12) and (13)). Although the average Biot coefficient of reservoir was calculated as 0.82 through 1D MEM, a large section of the log profile consists of overburden formations which are less consolidated relative to reservoir section, for this reason Biot coefficients close to 1 are used for these layers.

4.2. 3D MEM

The 3D MEM was built up using the static geological model described above, which comprises a high-resolution reservoir section and regions of lower resolution away from reservoir named as sideburden, overburden, and underburden. The higher resolution in the area of interest (reservoir) and the lower in the outside regions make up a grid which creates a balance between simulation precision and computational demand. Topography was extracted from the elevation maps of the ground level in order to include the top surface of the model. The horizons bounding the reservoir were used to generate overburden layers and underburden layers in the entire realm of the geomechanical model. The basal unit of the model comprises crystalline basement rocks at a depth of about 5 km. However, as none of the wells had actually reached this depth, this information is inferred from regional geological knowledge. The final 3D geomechanical model consists of 12 horizons and 11 lithostratigraphic units with dimensions of about 30 × 24 × 5 km³ in X, Y, and Z direction, respectively. The grid of the 3D MEM with the reservoir model embedded is shown in Figure 8. The calculated and calibrated log-derived properties including pore pressure, Young’s modulus, Poisson’s ratio, and density were upscaled from the well locations to the entire model domain. Thereby, the kriging method was used to populate the 3D geomechanical model [11]. The ranges of important material properties such as Young’s modulus, Poisson ratio, bulk density, and Biot coefficient of the reservoir as well as over- and underburden regions which are listed in

Table 1. Material properties of the reservoir along with overburden and underburden units of 3D geomechanical model.

Stratigraphic Unit	Young’s Modulus (GPa)	Poisson Ratio (-)	Density (g/cm3)	Biot Coefficient (-)
Quaternary, sweet-water molasse and Helvet	0.23–0.56	0.34–0.40	1.91–2.22	0.90–0.97
Burdigal	2.79–3.80	0.32–0.34	2.22–2.42	0.90–0.94
Top Aquitansand (Aquitan II) to Upper Aquitan down	1.38–4.2	0.32–0.36	2.43–2.47	0.89–0.95
Aquitansand (Aquitan I) + Chatt Hangende Tonmergel	2.71–3.70	0.32–0.35	2.45–2.49	0.90–0.92
Chatt Hauptsand	2.97–5.58	0.30–0.32	2.32–2.55	0.88–0.91
Chatt liegende Tonmergel, Rupel + Lattorf Fischschiefer	4.99–9.75	0.29–0.30	2.53–2.76	0.81–0.88
Lithothamnienkalk	3.58–16.09	0.29–0.34	2.77–2.85	0.77–0.90
Malm Zeta to Top Cretaceous	25–40	0.20–0.30	2.50–2.682	0.7–0.90
Malm Gamma to Top Malm Zeta	20–35	0.25–0.28	2.5–3.0	0.7–0.90
Top Basement to Top Malm Gamma	30–45	0.20–24	2.5–3.0	0.7–0.90
Basement	45	0.25	2.7	0.90

. The precision is of course decreasing with increasing distance to the wells, but the model fits well with overall trends. The initial pore pressures were also upscaled and interpolated from the 1D MEM’s, but in the coupling stage later-on pore pressures of different time steps are taken from the fluid flow model.

Boundary Conditions and Stress Orientations

Strain increments of ${\epsilon }_h$ = 0.0009 and ${\epsilon }_H$ = 0.0012 derived from 1D modeling have been applied on the lateral side of the model which provides a simulated regional stress state. Displacements in north-south direction lead to magnitude of $S_{Hmax}$ and displacements in west-east direction lead to $S_{hmin}$. The direction of the horizontal stresses can be estimated from the analysis of borehole breakouts on image logs and caliper data or earthquake focal mechanisms as well as from drilling induced tensile fractures. The orientations of horizontal stresses have been calibrated according to the World Stress Map (Figure 9) [4].

5. Results

5.1. 1D MEM

Figure 10 shows the results of Well X6 with calibration data as an example for the 1D MEM output. Shear velocity was calculated by using the Greenberg and Castagna correlation (Equation (7)) [17]. The density log along with compressional and shear slowness were used to calculate dynamic and static elastic properties by using the Bradford correlation. The average elastic Young’s modulus and Poisson’s ratio are 4 GPa and 0.35, respectively, at the X6 location at a depth of about 1786 m. The calculated material properties were calibrated using laboratory measurements. The pore pressure was calibrated first to get a gradient profile matching the calibration point at corresponding depth. Both calculated and measured pressures of 16 MPa coincide at a depth of 1786 m. The pore pressure profile follows the hydrostatic pressure until the top of the reservoir section, which shows an over-pressured zone slightly surpassing the hydrostatic pressure by 5 MPa. The elastic properties at the same depth were also calibrated against measured data. At last, the horizontal stresses were calculated using a poro-elastic approach and calibrated against measured stresses by using strain increments (Equations (14) and (15)) i.e., ${\epsilon }_h=$ 0.0009 and ${\epsilon }_H$ = 0.0012.

5.2. Fluid Flow (Dynamic) Model

The dynamic model incorporates the gas production history of the case study reservoir. Production started in 1958 and lasted until 1976. The shut-in phase then started and replenish the field pressure (up to 15.3 MPa). The corresponding history match has been performed and the calibrated production history is shown in the Figure 11.

5.3. 3D MEM

Petrophysical properties, such as porosity and permeability, were updated in the 3D geomechanical model from the dynamic model through history matching. The mechanical properties were upscaled and populated by the Kriging method [11]. The spatial distribution of Young’s modulus is presented in Figure 12. The mechanical properties indicate distinct vertical variations, but only slight lateral changes due to changes in lithology. The static Young’s modulus in the overburden ranges from 2 GPa to ~8 GPa near the bottom. The western and eastern part of the top reservoir layer shows higher Young’s moduli of about 10 to 12 GPa, which increase to the east where they reach up to 16 GPa. The average value of the Young’s modulus decreases up to 6 to 8 GPa in the middle part of the reservoir near wells X6, X6a, and X1. The higher Young’s modulus and lower Poisson ratio in the eastern part of the reservoir indicate more shaley lithologies in this area. In contrast, the decrease in Young’s modulus and increase in Poisson ratio in the middle and western part indicate the presence of more sandstones than in the eastern part of reservoir. The higher Young’s modulus in the underburden area may be due to higher compaction and presence of basement rocks close to the base of model.

The stress state in the 3D MEM was achieved by calibrating minimum and maximum horizontal stresses by iterative adjustment of strain increments (both minimum and maximum strains). The calibrated stresses are shown in Figure 13 for well X6. The reservoir zone lies in a normal faulting regime as the vertical stress $Sv$ (= S₁) is larger than the maximum horizontal stress $S_{Hmax}$ (= S₂) which in turn is larger than the minimum horizontal stress $S_{hmin}$ (= S₃). The calibration points match the calculated stresses, i.e., 38.2 MPa and 36.7 MPa for $S_{Hmax}$ and $S_{hmin}$, respectively [33]. Likewise, the calculated vertical stress $S_v\approx 42$ MPa in the reservoir zone fits the observed vertical stress of 42.2 MPa from [33]. However, the vertical stress in large parts of the underburden and overburden sections becomes the smallest among the principal stresses which means a thust faulting regime above and below the reservoir. The stress state in the section above 600 m (below sea level) to the surface is $S_{Hmax}$ > $S_v$ > $S_{hmin}$ making it a strike-slip region. The abrupt changes in the tectonic regime with depth are a consequence of the different mechanical properties in the model domain.

The magnitudes of the three stresses $S_{Hmax}$, $S_{hmin}$ and $S_v$ in the immediate reservoir area are presented in Figure 14. $S_{Hmax}$ ranges from 38 MPa to 45 MPa and tends to increase in the eastern part of the reservoir as effect of the deeper reservoir. Similarly, $S_{hmin}$ and $S_v$ are also showing this trend of increasing stresses towards the eastern side of the reservoir. $S_{hmin}$ is in the range of 36 MPa and 44 MPa with a minimum value at the well location X6 at 1786m depth. Moreover, the vertical stress ${\sigma }_v$ reaches up to 48 MPa to the lower eastern edge of the reservoir. The typical ratio of ${\mathrm{S}}_{Hmax}$ to $S_v$ is about 0.91 at the X6 location whereas the $S_{hmin}$ to $S_{Hmax}$ ratio is roughly 0.96. The ${\mathrm{S}}_{Hmax}$ orientation is similar to the data from the World Stress Map which is trending N-S in the wider region of the reservoir [4].

The validated 3D MEM gives opportunity to assess the drilling conditions (potential and risk) of a new well at any location within the model domain. Planning of well trajectory requires information on the complete stress tensor, i.e., the magnitude as well as the orientations of all three principal stresses. This is important especially close to faults and other complexities. Such a comprehensive stress prognosis can be provided by the 3D MEM as is exemplified by the stress profile for a randomly positioned well shown in Figure 15. This stress prediction method based on 3D MEM’s is universally applicable.

Another outcome of the 3D MEM is the displacements which are calculated for the Earth’s surface. The vertical displacement of the Earth’s surface above the reservoir was predicted maximum at the ultimate depletion time (t1) in 1978, i.e., nine years after production, and the upheaval was observed during the maximum replenishment time (t2), which took place in 1990, i.e., 12 years after the maximum depletion time (t1). The subsidence of the ground surface at t1 is due to the compaction processes during the depletion phase, while the slight upheaval during t2 is due to increase in pore pressure by the active aquifer below the gas-bearing layers [9]. Since the displacement is a function of all cumulative displacements of all underlying cells, the vertical displacement is greatest at the ground surface directly above the reservoir area [9]. The maximum subsidence of the ground surface above the reservoir is observed during t1 with −13 mm (negative sign stands for subsidence), which decreases laterally on each side and consequently becomes zero at the edges of the reservoir. During t2 an inverse surface behavior was observed, which led to a slight upheaval of 11 mm at the ground surface above the reservoir. The translation of the reservoir deformation to the ground surface is shown in Figure 16.

6. Conclusions

A workflow is presented that uses 1D mechanical earth models (1D MEM) to calculate static elastic properties and pore pressures based on well log data information and provide a first approximation for the stress distribution with depth. The 1D modeling results are then used to populate a 3D geomechanical model (3D MEM) with elastic properties by using the kriging interpolation method. The geomechanical model is coupled with a dynamic hydraulic model to get the pore pressure field for the reservoir for certain time steps in order to simulate total and effective stresses and related deformation inside and outside of the reservoir as well as surface displacement.

The workflow was successfully applied to a case study gas reservoir located in the Molasse Basin of southern Germany. Modeling results were calibrated against various observational data sets. The validated 3D geomechanical model provides the complete stress state at any location in the model domain. This information can be used, for example, for planning of future well trajectories. Furthermore, model results provide the starting point for future scenario testing, e.g., conversion of the former gas field to a storage site with short-term (weekly) injection-production schedules on top of long-term (seasonal) gas storage cycles. The modeling results also provide information about surface ground movement during depletion and replenishment times. The workflow outlined and tested in this study is not site-specific but generally applicable to any gas storage in a porous medium including methane, CO₂, and hydrogen.