3-D geomechanical modelling of a gas reservoir in the North German Basin: workflow for model building and calibration

Introduction Conclusions References


Introduction
The tectonic stress field strongly affects the optimal exploitation of conventional and unconventional hydrocarbon reservoirs.Among others, wellbore stability, orientation of hydraulically induced fractures and -particularly in fractured reservoirs -permeability anisotropies depend on the recent 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., Introduction Conclusions References Tables Figures

Back Close
Full 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, respectively (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 regional scale 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 utilizing FE techniques to build and calibrate geomechanical reservoir models.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 non-linear 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 represents the data Introduction

Conclusions References
Tables Figures

Back Close
Full basis used for the build-up of the geomechanical model.Calibration data is a separate, independent dataset and includes various stress and fracture measurements, which are exclusively used to compare the modelling outcome to these field observations.Input data comprises 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 boundary condition for the numerical model.
Following the geometrical build-up of the model inside the FE code, the fault block volumes are discretized 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 re-calculation 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 Introduction

Conclusions References
Tables Figures

Back Close
Full reservoir model that is geometrically consistent with all available data, e.g. a Petrel ® project.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 had to be found.
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 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 (fault-block size) submodels of a reservoir.Not only does the extensive effort preclude this approach to be used for field-scale models, but also their inevitable larger cell size prevents the 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 re-generated 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.Introduction

Conclusions References
Tables Figures

Back Close
Full 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 also used in the build-up of the case study model described below.
The re-generated 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 discretizised into a finite number of elements.The smaller the elements are, the higher becomes the spatial resolution of the results.However, the total amount of elements and thus the computational effort are 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.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 also 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 Introduction

Conclusions References
Tables Figures

Back Close
Full sides of the pre-assigned 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 allows for elastic and plastic deformation described by Hooke's and the Mohr-Coulomb laws, respectively, but also for temperature-and/or strain rate -dependent 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, frequency-dependent 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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, respectively.
The boundary conditions representing the ambient stress field complete the description of the model.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).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, drilling induced tensile fractures and ultrasonic wave velocity analyses (Ljunggren et al., 2003).Borehole breakouts can be determined by caliper logs, while breakouts and drilling induced fractures can both be recognized on image logs of the wellbore.The Introduction

Conclusions References
Tables Figures

Back Close
Full • (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 is increasing 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.
A model is regarded as validated as soon as all reliable calibration data is 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.Introduction

Conclusions References
Tables Figures

Back Close
Full However, it is important to note that at any stage, all data regarded for input and calibration inherently comprises inaccuracies and errors.This includes seismic data and its 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 localizations 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 to 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 recent stress conditions.

Conclusions References
Tables Figures

Back Close
Full

Case study
The workflow outline 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 most economic gas reservoirs throughout the Southern Permian Basin area (Doornenbal and Stevenson, 2010;Glennie, 2007).
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. 2b).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 distant and hypothetical faults are left out.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 s wave 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 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 (high-performance computing) devices to achieve calculation times of less than a day.Two systems are 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) parallelization 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 focus-Introduction

Conclusions References
Tables Figures

Back Close
Full ing 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, but the horizontal stress magnitudes as well.
The final modelling outcome of the calibrated geomechanical model includes, among others, 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 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 were 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 summarizes 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 is 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 Introduction

Conclusions References
Tables Figures

Back Close
Full

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 kind 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 is successfully applied to a Rotliegend gas field in the North German Basin.A geomechanical reservoir model, which was 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 referred to critically stressed faults and eventually to spatial probabilities of seismicity.Introduction

Conclusions References
Tables Figures

Back Close
Full Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Screen / Esc Printer-friendly Version Interactive Discussion Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | determination of stress orientations by wave velocity analyses requires the extraction of the oriented drill core.All those measurement methods comprise errors ranging between 10-15 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | stress magnitudes largely suppresses strong perturbations in stress orientation.Only close to faults or in overstep regions between faults pronounced orientation changes occur.

Fig. 1 .Fig. 2 .Fig. 3 .Fig. 4 .Fig. 5 .Fig. 6 .
Fig. 1.Summary of the workflow used to build 3-D geomechanical FE models.Data that is taken as input must not be used for calibration purposes to avoid circular reasoning.After calibration the validated geomechanical model can be used for stress prediction in the inter-well space or undrilled parts of the reservoir.