A workflow for building and calibrating 3-D geomechanical models – a case study for a gas reservoir in the North German Basin

Abstract. The optimal use of conventional and unconventional hydrocarbon reservoirs depends, amongst other things, on the local tectonic stress field. For example, wellbore stability, orientation of hydraulically induced fractures and – especially in fractured reservoirs – permeability anisotropies are controlled by the present-day in situ stresses. Faults and lithological changes can lead to stress perturbations and produce local stresses that can significantly deviate from the regional stress field. Geomechanical reservoir models aim for a robust, ideally "pre-drilling" prediction of the local variations in stress magnitude and orientation. This requires a numerical modelling approach that is capable to incorporate the specific geometry and mechanical properties of the subsurface reservoir. The workflow presented in this paper can be used to build 3-D geomechanical models based on the finite element (FE) method and ranging from field-scale models to smaller, detailed submodels of individual fault blocks. The approach is successfully applied to an intensively faulted gas reservoir in the North German Basin. The in situ stresses predicted by the geomechanical FE model were calibrated against stress data actually observed, e.g. borehole breakouts and extended leak-off tests. Such a validated model can provide insights into the stress perturbations in the inter-well space and undrilled parts of the reservoir. In addition, the tendency of the existing fault network to slip or dilate in the present-day stress regime can be addressed.


Introduction
The tectonic stress field strongly affects the optimal exploitation of conventional and unconventional hydrocarbon reservoirs.Among other things, wellbore stability, orientation of hydraulically induced fractures and -particularly in fractured reservoirs -permeability anisotropies depend on the present-day in situ stresses.Information on the regional stress orientations can be derived from large-scale data collections like, for example, the world stress map project (Zoback, 1992;Sperner et al., 2003).However, stress magnitudes and orientations are frequently not homogeneous on a reservoir scale, but can be substantially modified by the presence of faults as well as lithological changes and contrasts in rock mechanical properties (Fjaer et al., 2008;Zoback, 2007).In some fault-controlled reservoirs, local stress reorientations of up to 90 • relative to the regional trend have been reported (Maerten et al., 2002;Yale, 2003).In such cases, inference of local in situ stress orientations from regionalscale maps would inevitably lead to an incorrect pre-drilling prediction.Therefore, any robust prognosis has to incorporate the specific 3-D geological reservoir structure including faults as well as the specific rock mechanical behaviour of the reservoir under concern.Such complexities can only be treated adequately by a numerical modelling approach.3-D geomechanical reservoir models based on the finite element (FE) method have been proven to be valuable tools to gain a quantitative understanding of the in situ stresses in a reservoir (van Wees et al., 2003;Henk, 2009Henk, , 2010)).In this paper, we present a detailed workflow utilising FE techniques to build and calibrate geomechanical reservoir models (Fig. 1).Subsequently, the workflow is applied to a gas reservoir in Northern Germany.The reservoir rock is formed by a faulted aeolian sandstone from the Upper Rotliegend (Doornenbal and Stevenson, 2010).It makes an ideal case study for geomechanical modelling, as several data sets are available to set up the numerical simulation and compare model predictions to stress data actually observed.

General workflow
Finite element (FE) techniques were used to gain quantitative insights into the stress and strain distribution of reservoirs, as this numerical approach allows for robust simulations of heterogeneous structures with complex geometries and nonlinear material behaviour.Geomechanical modelling is based on the FE software ANSYS ® (Ansys Inc., Houston, USA).It requires two independent data sets -input and calibration data -and involves several work steps (Fig. 1) (Henk, 2005).Input data represent the data basis used for the build-up of the geomechanical model.Calibration data are a separate, independent data sets and include various stress and fracture measurements, which are exclusively used to compare the modelling outcome to these field observations.
Input data comprise information on the specific reservoir geometry, rock mechanical parameters as well as the regional stress field.The reservoir geometry includes faults and lithostratigraphic boundaries that can be derived from corresponding maps or 3-D geological models.Material parameters should be as reservoir-specific as possible and describe the mechanical properties of all participating lithologies and faults.The regional stress field of the reservoir area serves as a boundary condition for the numerical model.
Following the geometrical build-up of the model inside the FE code, the fault block volumes are discretised and populated with the respective material properties.After the initial calculation, modelling results need validation and hence stresses calculated by the model are compared to calibration data, i.e. stress orientations and magnitudes actually observed in well data.If required, poorly constrained parameters, such as the friction coefficient of faults, can be iteratively adjusted within geologically reasonable limits until a satisfactory fit between modelling results and measurements is obtained.However, each parameter change and corresponding iteration requires a recalculation of the entire FE model.
The final results of the calibrated geomechanical model include the complete 3-D stress and strain tensor for any location in the model, as well as shear and normal stress acting on the faults.Based on that, further stress quantities like, for example, mean stress and differential stress as well as the slip and dilation tendency of the faults can be calculated (Morris et al., 1996).

Geometry transfer
The most labour-intensive step in the build-up of a geomechanical FE model addresses the reservoir geometry.Lithostratigraphic horizon and fault surfaces are commonly derived from interpretation of 3-D seismics and are ideally available as a geological reservoir model that is geometrically consistent with all available data, e.g. a Petrel ® project.The most accurate option uses a different approach, i.e. a high-resolution point cloud and reverse engineering techniques for surface reproduction.This approach is by far the most labour-intensive approach and not suitable for field-scale models.
The finite difference grid frequently used for property modelling and flow simulations cannot be used for geomechanical FE modelling.The FE mesh has to fulfil special requirements, for instance dual grid nodes at fault faces (see below) and non-orthogonal material boundaries.Thus, the reservoir geometry has to be rebuilt in the FE modelling software to maintain the geometrical complexities of the subsurface reservoir structure.However, there is no direct way to export and import horizons or faults as surfaces from the geological modelling software to the FE program.Due to this lack of direct interfaces, other ways for transferring the reservoir geometry are required.
Three options were developed that differ in effort and in the accuracy of the regenerated surfaces (Fig. 2).One option is the extraction of high-resolution point clouds from the geological reservoir model describing, for instance, lithological boundaries.These point clouds are triangulated and processed to freeform surfaces of horizons using 3-D CAD (Computer Aided Design) and reverse engineering software.
This approach yields the highest accuracy, as basically no geometrical information is lost.However, it is also the most labour-intensive way and is suitable only for detailed (faultblock size) submodels of a reservoir.Not only does the extensive effort preclude this approach from being used for fieldscale models, but also their inevitable larger cell size prevents accurate consideration of the highly detailed horizon geometries.
Another option for geometry transfer uses the extraction of points only along horizon lines, i.e. the intersection lines between lithological and fault surfaces (Fig. 2).Based on these points, 3-D splines can be regenerated in the FE software following exactly the horizon lines.These splines are then used to create a so-called Coon's patch (Barnhill, 1982;Goodman and O'Rourke, 2004), which is a surface whose topology is interpolated by the slope of the bounding lines.This approach represents the fastest way to transfer the reservoir geometry to the FE software.However, the larger the created Coon's patch is, the more information on the internal topology of the surface is lost.This option is best suited for field-scale models with no pronounced horizon topology.
The third option represents a modification of the latter approach.In addition to the regular horizon lines, a network of auxiliary lines is created in the geological model.They follow the surface topology and are transferred to the FE software in the same way as the horizon lines.This yields significantly smaller Coon's patches and, thus, the internal topology of any face can be preserved much more accurately (Fig. 2).Especially regarding field-scale geomechanical reservoir models, this option offers the best compromise between effort and accuracy, and was used in the build-up of the case study model described below.
The regenerated surfaces representing lithostratigraphic boundaries and faults are subsequently combined to volumes of individual fault blocks comprising multiple mechanical layers (van Wees et al., 2003).These volumes are then discretised into a finite number of elements.The smaller the elements are, the higher the spatial resolution of the results.However, the total amount of elements and thus the computational effort is increased as well.

Incorporation of faults
Faults can be considered in the geomechanical model in two ways.In the first approach, the entire model is meshed continuously and those elements that are cut by fault faces receive significantly lower mechanical parameters than the surrounding area.Faults are thus represented by zones of weakness inside the numerical model (Herwanger and Koutsabeloulis, 2011).However, this approach limits potential displacements along faults due to the distortion of the mesh.Furthermore, no shear or normal stresses acting on the faults can be obtained from the modified 3-D elements, as these stresses demand, by definition, a surface to refer to.The lack of a fault surface implementation thus limits the results of this approach to the general 3-D stress tensor of the weak elements.
The second approach, which is used in the case study, implements existing faults as distinct planes of weakness cutting the finite element model.These planes are described by so-called contact elements.Contact elements are defined at opposing sides of the preassigned faults.Thus, dual grid nodes are required along the fault faces.Stiffness values of the contact similar to the Young's moduli of the rocks are used to enforce compatibility between adjacent fault surfaces.In geological terms, this contact stiffness can be related to a fault zone thickness, or damage zone of a fault, and can be regarded as the reciprocal value of the so-called fault zone compliance (van der Neut et al., 2008).
The contact elements can transmit shear and normal stress and are capable of describing frictional sliding and non-linear behaviour.During calculation, the state of contact thus can change from sticking to sliding in case stresses overcome the cohesive and frictional strength.The use of contact elements allows for large movements between the different, individually meshed parts of the model, but does not describe fault propagation itself.The latter would require automatic mesh modifications, i.e. mesh split-up and assignment of additional contact elements in areas where peak stress exceeds strength criteria.However, such an approach is beyond the capabilities of classic FE codes.
The usage of contact elements for fault description significantly enhances the classical continuum approach of the finite element method.Fault-specific shear and normal stresses can be retrieved and used to calculate slip and dilation tendencies of the faults that indicate their movement behaviour.

Material parameters and boundary conditions
The FE method not only allows for elastic and plastic deformation described by Hooke's and the Mohr-Coulomb laws, respectively, but also for temperature-and/or strain ratedependent creep.Depending on the material behaviour to be modelled, the corresponding reservoir-specific material parameters for all participating lithologies have to be obtained, for instance, from geomechanical logs and rock mechanical tests on drill cores.For numerical modelling, frequencydependent dynamic properties that are based on sonic wave propagation must be converted to their static equivalents, i.e. static Young's moduli and static Poisson's ratios.These static values are mandatory in numerical simulations to describe deformation at low strain rates.The conversion can be done by using reservoir-specific or published empirical correlations (e.g.Eissa and Kazi, 1989).
Assigning the required material parameters to the respective layers and elements establishes a mechanical stratigraphy.However, a more detailed mechanical stratigraphy has to be acknowledged by a corresponding reduction in element size and refined spatial resolution.
The boundary conditions representing the ambient stress field complete the description of the model.Information on the regional stress orientations can be derived from largescale data collections like, for example, the world stress map project (Zoback, 1992;Sperner et al., 2003).Boundary conditions are applied as nodal displacements to a rectangular frame oriented parallel and perpendicular to the orientations of the minimum and maximum horizontal stress, respectively.The displacements are calibrated to generate the desired stress magnitudes for both horizontal stresses.While the bottom of the model is fixed in vertical direction, a lithostatic pressure load is applied to the model's top, representing the load of the overburden at a distinct depth level.This load can also be varied laterally to consider, for instance, lateral changes in the density distribution of the overburden.Finally, gravity acts as a body force on the entire model.

Model calibration
Before any conclusions can be drawn from the geomechanical model, it has to be calibrated, i.e. the modelling outcome has to be compared to stress data actually observed.Subsequently, the fit is improved by iteratively adjusting poorly constrained parameters, such as the friction coefficient of the faults and the magnitude of the maximum horizontal stress, within geologically reasonable limits.
The types of data that can be used for calibrating the geomechanical model are manifold.Stress orientations can be derived, for instance, from borehole breakouts, drillinginduced tensile fractures and ultrasonic wave velocity analyses (Ljunggren et al., 2003).Borehole breakouts can be determined by caliper logs, while breakouts and drillinginduced fractures can both be recognised on image logs of the wellbore.The determination of stress orientations by wave velocity analysis requires extraction of the oriented drill core.All these measurements comprise errors ranging from 10 • -15 • (Grote, 1998).With respect to stress magnitudes, the least principal stress magnitude is commonly inferred from extended leak-off tests, hydraulic fracturing, etc. (Ljunggren et al., 2003).The magnitudes of both minimum and maximum horizontal stresses can be determined from ultrasonic wave velocity analyses and anelastic strain recovery (Ljunggren et al., 2003).While contour and vector plots of stress quantities provide good overviews, the comparison of the modelling outcome with measured well data described above requires the knowledge of the exact subsurface location, especially in cases of deviated or horizontal wells.The consideration of precise well trajectories becomes even more important for wells located in the vicinity of faults, as the intensity of perturbations in stress magnitudes and orientations increases towards the fault inducing them.Therefore, the coordinates of all well paths in the reservoir are imported into the FE program and all elements cut by the path of a specific well are selected (Fig. 3).Magnitudes and direction cosines of the principal stresses of those elements are then extracted from the geomechanical modelling results.While magnitudes can be directly compared to their measured counterparts, the direction cosines are used to calculate an azimuth value of the maximum horizontal stress, for instance.

Solid
A model is regarded as validated as soon as all reliable calibration data are reproduced within the assumed ranges of measurement errors.Following this calibration process, the validated geomechanical model can then be used for stress and strain predictions in the inter-well space and undrilled parts of the reservoir.As it is the case with all numerical models, the predictive quality of a geomechanical reservoir model depends on the availability and quality of input and calibration data.Thus, in the early exploration stage when well data will be sparse, modelling results have a higher uncertainty.Any new data coming up during the subsequent appraisal and development stages can be incorporated in order to improve the model.
However, it is important to note that at any stage, all data regarded for input and calibration inherently comprise inaccuracies and errors.This includes seismic data and their interpretation, providing the reservoir geometry, the derivation of material parameters, and all stress measurements.Therefore the uncertainty of all data must be kept in mind at all times during modelling and subsequent interpretation.

Post-processing
The final results of the calibrated geomechanical model comprise displacements as well as the complete 3-D stress and strain tensors for any subsurface location inside the FE model.Further stress and strain quantities can be calculated on demand, such as differential and mean stress.Stress magnitude distributions can be displayed by contour plots, while stress orientations are commonly shown by vectors.This comprehensive modelling outcome can be used for a variety of real-world applications.Strain localisations can be referred to increased fracture intensities in these areas and indicate locally increased probability of sweet spots.The provided stress orientations can help to ensure wellbore stability of newly drilled wells, while the magnitude of the least principal stress in combination with its orientation improves the planning of optimal oriented hydraulic fracs.Mean stress distributions may yield insights into potential hydrocarbon migration directions throughout the reservoir.
The contact elements used for fault representation provide shear and normal stresses.They can be used to calculate slip and dilation tendencies indicating the fault's behaviour to move and open.This helps, for example, to identify critically stressed and more permeable fault segments (Townend and Zoback, 2000) under present-day stress conditions.

Case study
The workflow outline presented above is applied to a mature gas field in the North German Basin, which is part of the much larger Southern Permian Basin.The reservoir rock is formed by an aeolian sandstone from the Upper Rotliegend (Early Wuchiapingian) (Fig. 4).This aeolian facies provides  The objective was to build a field-scale geomechanical reservoir model encompassing the entire reservoir area with numerous faults.Lithostratigraphic horizons and faults interpreted from 3-D seismic data were available in a geological reservoir model, in this case a Petrel ® project.The geometry was transferred to the FE code using the approach in which a network of Coon's patches is applied (see Fig. 2, B).In this way, the topology of the reservoir horizon was preserved accurately with reasonable effort.The reservoir horizon was embedded between thick over-and underburden layers to avoid boundary effects on the area of interest.An unfavourable boundary effect could result, for example, from the fixed bottom of the model if it is only a short distance away from the base of the reservoir unit under concern.A sufficiently thick underburden layer ensures that elastic rebound from beneath the reservoir is taken into account.
Besides the multiple stratigraphic layers, more than 80 faults are incorporated in the geomechanical model and simulated by contact elements.This represents the complete fault network in the production area of the reservoir as it is interpreted from 3-D seismic data.Only very small conceptual faults inferred from reservoir simulation are left out and those in great distance to the production area.By incorporating the entire fault network, the complex interaction of stress perturbations resulting from those faults can be addressed in the model.
Reservoir-specific material properties for all relevant lithologies were adopted from rock mechanical tests on drill cores as well as geomechanical log data providing p-and swave velocities.The frequency-dependent dynamic properties were transformed to reservoir-specific static values using the results of the mechanical tests.Vertical boundary conditions comprise a lithostatic pressure load derived from the integrated density of the overburden and incorporate lateral variations due to salt structures.Magnitudes and orientations of the regional maximum and minimum horizontal stresses were taken from Röckel and Lempp (2003) and applied to the model by calibrated displacements yielding the respective stress magnitudes.The bottom of the model is fixed with respect to vertical movements and gravity is applied as body force.The underburden layer comprises more than seven times the thickness of the reservoir layer and thus prevents any negative impact of the bottom fixation of the model on the in situ stress inside the reservoir.This justifies the assumption of a fixed bottom, which is numerically mandatory.
The final geomechanical reservoir model covers an area of more than 400 km 2 and comprises about 4 million elements.Spatial resolution and grid cell size inside the reservoir region is 100 m × 100 m × 25 m (length × width × depth).The total amount of elements and the complexity of the geomechanical model, especially regarding the faults, required the use of massive parallel computing techniques on HPC (highperformance computing) devices to achieve calculation times of less than a day.Two systems were used for calculation: the compute cluster of the bwGRiD1 at the Albert-Ludwigs-Universität Freiburg and an in-house server of the group of engineering geology at the Technische Universität Darmstadt.SMP (symmetric multi-processing) parallelisation was used in both cases, as well as computer servers comprising 12-32 cores and 192-512 GB of main memory.
Modelled stresses were compared to data from field measurements.The calibration data available for this case study included stress orientations from borehole breakouts and ultrasonic wave velocity analyses (Grote, 1998).Least principal stress magnitudes were provided by hydraulic frac reports.Ten different configurations of the model focusing on adjustments of poorly constrained parameters were iteratively calculated until a configuration yielded a satisfactory match between modelling outcome and field observation.This best fit model comprises a lowered fault friction coefficient and adapted lateral heterogeneities in the overburden load due to salt bodies.Specific densities were provided by gravimetric modelling of the entire area.The resulting local reduction of the overburden pressure affected the vertical stress as well as the horizontal stress magnitudes.
The final modelling outcome of the calibrated geomechanical model includes, among other things, the 3-D stress tensor throughout the model, from which various other stress quantities can be derived.Figure 5 shows the distribution of the minimum horizontal stress magnitude on the reservoir level in combination with the shear stress acting on the fault faces.This plot elucidates the spatial extent, intensity and origin of stress perturbations.Fault tips, strong curvatures and the interaction of neighbouring faults are the most important sources of perturbations.However, the juxtaposition of mechanically different layers at vertical displacements along faults can also lead to major variations in the stress distribution.
The modelled magnitudes of the least principal stress are compared to the reported magnitudes from hydraulic fracturing, while the magnitudes of the second principal stress are compared to those derived from ultrasonic measurements on drill cores.In the normal faulting regime prevailing throughout most of the reservoir area, these stresses represent the minimum and maximum horizontal stress.Figure 6 summarises this comparison by plotting all available and reliable measurements of maximum and minimum horizontal stress against the modelled stress magnitudes at the specific subsurface locations in the calibrated model.The fit is already satisfactory regarding uncertainties in input and calibration data, but could be further improved, for example, by incorporation of lateral variations in mechanical parameters.In addition to the magnitudes, the modelled orientations of the maximum horizontal stress are also compared to published data (Grote, 1998).Most of the model predictions, as well as the borehole observations, follow the regional NNW-SSE trend, as the large difference in horizontal stress magnitudes largely suppresses strong perturbations in stress orientation.Only close to faults or in overstep regions between faults do pronounced orientation changes occur.

Conclusions
The workflow presented above can be used to build 3-D geomechanical FE models for various types of reservoirs and ranging from field-scale models to smaller, highly detailed submodels of specific fault blocks.The workflow is applicable to all kinds of stress-sensitive reservoirs including conventional and unconventional hydrocarbon as well as geothermal reservoirs.All steps were elucidated in detail, including the transfer of reservoir geometry, incorporation of faults, assignment of material properties and boundary conditions, and the final calibration of the model.This workflow was successfully applied to a Rotliegend gas field in the North German Basin.
A geomechanical reservoir model that has been calibrated against well data can be used for stress predictions in the inter-well space and undrilled parts of the reservoir.This information on the in situ stress distribution can be used for an improved planning of numerous reservoir applications and may also provide further insights into the drainage pattern of the reservoir, since particularly in fractured reservoirs the tectonic stress causes fluid flow anisotropies (Heffer and Koutsabeloulis, 1995).In addition, the tendency of the existing fault network to slip or dilate in the present-day stress regime can be addressed.In very detailed models, this information can be applied to critically stressed faults and eventually to spatial probabilities of seismicity.

Fig. 1 .
Fig. 1.Summary of the workflow used to build 3-D geomechanical FE models.Data that are applied as input must not be used for calibration purposes to avoid circular reasoning.

Fig. 2 .
Fig. 2. Illustration of the three options for the transfer of an arbitrary surface from a geological subsurface model to the finite element model.(A) If bounding lines are used exclusively to create a Coon's patch, only little internal topological information is preserved in the regenerated surface.(B) By adding a network of auxiliary lines, the resulting Coon's patches are significantly smaller and the internal topology can be preserved more accurately.(C) The most accurate option uses a different approach, i.e. a high-resolution point cloud and reverse engineering techniques for surface reproduction.This approach is by far the most labour-intensive approach and not suitable for field-scale models.

Fig. 3 .
Fig. 3. Example of the selection of elements and modelling results along the drilling path of a highly deviated well.Elements are coloured according to their layer (top) and magnitude of minimum horizontal stress (Shmin) (bottom).

Fig. 4 .
Fig. 4. Oblique view of the depth-coloured top reservoir horizon of the case study (10 × vertical exaggeration).Black faces represent the faults.For geomechanical modelling the reservoir horizon is embedded between over-and underburden layers, thus forming a regular block in the subsurface with planar top and base.In total, more than 80 faults were considered.

Fig. 5 .
Fig. 5. Contour plot combining the reservoir-wide distribution of the magnitude of the minimum horizontal stress (Shmin) and the shear stress acting on the fault faces (both in MPa at mid-reservoir level).The plot illustrates that stress magnitude perturbations can range at least several hundreds of metres from the inducing fault.Some perturbations are induced by mechanically different layers juxtaposed to each other across the fault.

Fig. 6 .
Fig. 6.Correlation between measured and modelled magnitudes of the minimum and maximum horizontal stress.Measured stresses were derived from hydraulic fracturing and ultrasonic wave velocity analyses on drill cores.For a perfect fit, all points would lie on the bisecting line of the diagram.