## Abstract

In mature field appraisal and development, discretized geomechanical property models play a vital role in providing information on in situ stress regime as a guide for placement of directional wells. Laboratory methods of measuring these properties, in most cases, take only small samples from consolidated rocks. These isolated samples may not be representative of the entire elastic regime existing in the reservoir owing to sample size. In general, geomechanical studies are performed on a well-by-well basis and then these measurements are used as calibration points to convert 3D seismic data (if available) to geomechanical models. However, elastic properties measured this way are restricted to the well location and interpolation across the reservoir may not be always appropriate.

#### Author

Kalyan Saikia^{1}, Clara E. Ikuku^{2}, Bhabesh C. Sarkar^{3}

^{1}.Halliburton Richmond USA. ^{2}.Clenik Petrotech Calgary Canada. ^{3}.Department of Applied Geology Indian Institute of Technology (Indian School of Mines) Dhanbad India.

Received: 14 September 2016 / Accepted: 20 November 2017

© The Author(s) 2017

To overcome these challenges, this paper describes an integrated approach for deriving 3D geomechanical models of the reservoir by combining results of 3D geocellular models and basin models. The basin model reconstructs the geologic history (i.e., burial history) of the reservoir by back-stripping it to the original depositional thickness. Through this reconstruction, the mechanical compaction, pore pressures, effective stress, and porosity versus depth relationships are established. Next, these mechanical properties are discretized into 3D geocellular grid using empirical formulas via lithofacies model even if no 3D seismic data are available for the reservoir. The discretization of elastic properties into 3D grids results in a better understanding of the prevailing stress regimes and helping in design of hydraulic fracturing operations with minimal risks and costs. This approach provides an innovative way of determining effective horizontal stress for the entire reservoir through 3D distribution of elastic properties.

## Introduction

In unconventional tight reservoirs, where hydraulic fracturing is the key technique for enhanced oil production, a thorough knowledge of distribution of geomechanical properties helps to maximize return on investment. Horizontal wells with multistage hydraulic stimulation are the primary production strategies during development of the unconventional tight reservoirs. It is therefore important to understand in situ stress conditions such as effective minimum horizontal stress as it guides the stimulation design which will in turn control the fracture propagation. Accuracy of effective minimum horizontal stress model is primarily dependent on knowledge of the spatial distribution of reservoir’s elastic property.

In general, for upstream reservoir modeling processes, reservoir geomechanical properties are determined either by (1) analysis of seismic data (S-wave and P-wave velocities) or (2) by laboratory measurement of core plugs or (3) a combination of both techniques. These are conventional laboratory methods which may not be representative of the entire elastic regime existing in the reservoir as in most cases these small samples are taken from consolidated rocks and used to measure geomechanical properties. Furthermore, these measurements are performed on a well-by-well basis on core samples. Measurements taken at the core samples are often used in geomechanical analysis as calibration points to convert the 3D seismic data (if available) to 3D volumes of geomechanical properties.

Previous studies on the gross effect of mechanical properties on productivity of the unconventional reservoir clearly elaborated that the use of only seismic velocities and elastic properties measured on samples are not sufficient enough to describe spatial distribution pattern of elastic properties for entire reservoir. Havens (2012) carried out similar approach for Bakken petroleum system to establish the spatial behavior elastic properties and concluded that conventional approach of consideration of small core samples plugs in formulating 3D reservoir geomechanical model is not a practical solution. During fracture modeling of unconventional tight gas reservoirs, Gonzalez et al. (2014) and Deng et al. (2011) demonstrated that incorporation of 3D geomechanical model gives a more realistic representation of the orientation and geometry of hydraulic fractures compared to traditional way of semi-analytical calculations of fracture geomechanical properties. Other published studies on similar topics in the western Canadian sedimentary basin reveal that spatial distribution of rock mechanical properties is highly heterogeneous in nature and is result of complex basin-forming (depositional) process (Ferdous et al. 2015; Lavoie and Séjourné 2016; Tong et al. 2005).

To overcome these challenges, this paper presents a practical approach for deriving and discretizing geomechanical and other elastic properties in unconventional tight reservoir by integrating results of 3D geocellular and basin models. This technique provides a relationship between petrophysical properties and elastic rock properties of the reservoir rock, explaining their response to mechanical stresses at reservoir condition. Such relationship is very important for defining trends in the reservoir description. The current approach also incorporates wireline logs representing the entire sediment sequence which can be used as input for effective minimum horizontal stress modeling for the entire reservoir. The case study demonstrated here assumes that elastic properties measured conventionally are restricted to the well locations and cannot be interpolated across the reservoir if no 3D seismic data are available.

## Integrated modeling approach

The integrated approach of deriving discretized geomechanical model of the reservoir can be divided into three main steps. While the first two steps are highly dependent on the input data and selection of appropriate modeling algorithms, the third step is an integration of results of the first two steps.

### Step 1

This step involves creating a high-resolution 3D geocellular reservoir model of the petrophysical properties and depositional environment of the reservoir. Initially existing geological, geophysical, and well data are interpreted and incorporated to create a structural framework and empty geocellular grid. Next, facies curves derived from the well logs are used to generate cell-wise distribution of lithofacies within the reservoir layers employing advanced geostatistical algorithms. Lastly, spatial distribution model of petrophysical property in terms of porosity, permeability, saturation etc., is generated and constrained to the lithofacies model. This helps in capturing geological uncertainty while creating 3D petrophysical property distribution model for the reservoir.

While carrying out 3D geocellular modeling, mainly two geostatistical simulation algorithms, viz. sequential Gaussian simulation (SGS) and Plurigaussian simulation (PGS) were employed. Sequential Gaussian simulation (SGS) (Deutsch and Journel 1998) is one of the popular geostatistical algorithms used for interpolation of continuous variable in reservoir modeling. It was originally introduced as a solution to the smoothing problem of interpolation by kriging. Sequential Gaussian simulation algorithm reproduces a globally structured two-point statistics using variogram model, whereas kriging provides a best estimate at each location with minimum error variance and ignores estimates made at other locations. Thus, with regards to modeling of hydrocarbon reservoirs, SGS algorithms provide more relevant 3D reservoir models honoring global spatial variations in sample points.

The implementation of sequential simulation consists of reproducing the desired spatial properties through the sequential use of conditional distributions (Arpat 2005; Dimitrakopoulos and Luo 2004). Consider a set of N random variables Z(uα), α = 1,…, N defined at N locations uα. The aim is to generate L joint realizations {z(l)(uα),_ = 1, …, N} with l = 1, …, L of the N, conditional to n available data and reproducing the properties of a given multivariate distribution. To achieve this goal, the N-point multivariate distribution is decomposed into a set of N univariate conditional distributions:

F(u1, …, uN; z1, …, zN|(n)) = F (uN; zN|(n + N − 1)) × F (uN − 1; zN − 1|(n + N − 2)) × ··· × F(u2; z2|(n + 1)) × F(u1; z1|(n)),

Where F (uN; zN|(n + N − 1)) = Prob {Z(uN) ≤ zN|(n + N − 1)} is the conditional cumulative distribution function of Z(uN) given the set of n original data values and the (N − 1) realizations z(l)(uα), α = 1, …, N − 1 of the previously simulated values.

This decomposition allows generating a realization by sequentially visiting each node on the simulation grid. Multivariate Gaussian distribution model is the only model for which the above-mentioned decomposition is analytically available (Deutsch 2002). The conditional probability distribution in a Gaussian model is determined using simple kriging of each unknown value at any node uj given the original data and previously simulated values (n + j − 1) (Deutsch and Journel 1998). This concept is very conveniently used in SGS making it one of the popular algorithms in reservoir modeling. The whole process of SGS used in current study can be summarized in five steps:

1. Generate a random path through the grid nodes.

2. Visit the first node in that path and use kriging to estimate a mean and standard deviation at that node based on surrounding data (i.e., local conditional probability distribution).

3. Select at random a value from the distribution and set the node value to that number.

4. Include the newly simulated value as part of the conditioning data.

5. Repeat steps 1–4 until all grid nodes have a simulated value.

Truncated Gaussian simulation (TGS) (Matheron et al. 1987) and Plurigaussian simulation (PGS) (Galli et al. 1994) technique are two popular methods of reservoir modeling workflow for simulating lithotype and facies of sedimentary rocks. TGS technique is best suitable for reservoirs where the lithotypes occur in a sequential order and it helps in defining the geometry and internal architecture of the reservoir rocks (Armstrong et al. 2003). PGS is a natural extension of TGS and capable of handling complex geological situations when the lithotypes occur in complicated relationship rather than simple sequence.

The underlying principle of TGS and PGS simulations is to set up one or more simulations of Gaussian random fields at every grid point and then use some rules to convert these Gaussian values into lithotype indicators (Allard et al. 2012). In case of PGS simulation, generally two Gaussian random fields are used to define the lithofacies structure. Armstrong et al. (2003) distinguishes four steps in a PGS approach as:

1. Choosing the rock-type rules depending on the types of relations among lithotypes.

2. Estimating the parameter values. Two key factors control PGS including the thresholds at which the different Gaussians are truncated and the variogram model of the underlying Gaussian variable. The proportion of each facies, the “rock-type” rule, and the correlation between the two underlying Gaussian random functions determines the thresholds. Once the thresholds have been worked out, the variograms and cross-variograms of the indicators are calculated experimentally.

3. Generating Gaussian values at the conditioning points.

4. Simulating Gaussian values at each grid node with a usual Gaussian simulation algorithm.

### Step 2

In this step, a basin modeling study is initiated covering the reservoir area. The basin model reconstructs the geologic history (i.e., burial history) by back-stripping the reservoir to its original depositional condition. Through this reconstruction, the mechanical compaction, pore pressures, effective stress, and porosity versus depth relationships are established for the entire reservoir. Well data in the form of log curves constitute the main input for basin modeling process in association with other basin-related information such as regional temperature and pressure profile.

Basin modeling is a technique which allows to understand effect of physical and chemical processes in a sedimentary basin to generate hydrocarbon. Basin modeling techniques help to reconstruct the burial and temperature history of the basin through time and to understand source rock maturation and subsequently hydrocarbon expulsion and migration (Hanstschel and Kauerauf 2009; Peters 2009). The basin modeling process dynamically models changes in rock properties through geologic time by numerically solving coupled partial differential equations with moving boundaries on discretized temporal and spatial grids. This computation results in many outputs of different rock properties including properties such as vitrinite reflectance, temperature, effective stress, and porosity. In the final stage, the basin models need to be calibrated with existing data such as vitrinite reflectance, temperature, and porosity (AlKawai 2014) to make it more realistic.

Overall implementation of basin modeling process contains a wide range of mathematical algorithms and methods each of them appropriate for each “sub-model,” and detailed discussion of these goes beyond the scope of this paper. Additional information and overview of these numerical methods are provided by Press et al. (2002) and Hantschel and Kauerauf (2009). In basin modeling, numerical solution of differential equations is fundamental requirement and most challenging, complex and costly process to accomplish. In most of the scenarios, temperature and pressure are modeled with parabolic diffusion equations. Multi-dimensional differential equations with complicated geometries are usually solved with the finite element method, whereas finite differences are often used in special cases of approximately one-dimensional problems, such as simplified crustal layer models (Hantschel and Kauerauf 2009).

### Step 3

It is the integration step where results from basin modeling study (i.e., distribution of mechanical properties of reservoir rocks over geological time) are correlated with each cell of geocellular model using empirical relationships (Eberhert-Phillips et al. 1989; Ingram and Urai 1999; Horsrud 2001; Yarus and Carruthers 2014). Use of these empirical relationships provided a practical solution in discretizing geomechanical properties in the current study area. Other published evidences of use of these relationships in basin modeling process include Alkawai (2014), Szydlik et al. (2015), and Schneider et al. (1996). The output of this process is a discretized geomechanical model of the reservoir along with other petrophysical properties.

**Fig. 1** Diagrammatic representation of geocellular model and basin model integration

Thus, the methodology adopted in this study is a true integration of results and parameters of both geocellular and basin models. Diagrammatic representation of the integrated approach is provided in Fig. 1. A step-by-step description of the implementation of the entire modeling process for appraisal of unconventional mature field in western Canadian sedimentary basin is explained in the subsequent section of the paper.

## Geological background of the study area

The study area is located at the western boundary of Saskatchewan, Canada, close to the Alberta border in Hoosier field (Fig. 2). Previous studies by White (1969) identified two oil pools in the field. The study area falls in the North Hoosier oil pool at the west central Saskatchewan. In this field, the Bakken shale formation is the important formation due to its high hydrocarbon potential. The Bakken formation was deposited during the geological age of Late Devonian and Early Mississippian. It lies unconformably over the Big Valley formation (Upper Devonian) and is conformably overlain by the Madison group as shown in Fig. 3 (Zhang et al. 2016). Deposition of Bakken formation occurred through a series of onlap–offlap cycles during the Tamaroa sequence (Wheeler 1963).

**Fig. 2** Location Map of study area showing the Hoosier field. Modified after O’Connell et al. (2000).

Time-equivalent shale units of Bakken shale include the Exshaw/Banff in the Alberta Basin, the Woodford shale in the Anadarko Basin, the Chattanooga in the Southern Appalachian Basin, and the Antrim in the Michigan Basin (Meissner 1978). The Madison Group which overlies the Bakken formation is a thick carbonate sequence of Mississippian age. The rocks of Madison were deposited in a generally shallow marine setting as indicated by richly fossiliferous rocks and favor for accumulation of petroleum resources (Porter 1955). It has been divided into three formations, viz. Lodgepole, Mission Canyon, and Charles. The oldest Lodgepole consisting limestone and dolomites conformably overlies the Bakken formation is a major hydrocarbon-producing horizon (Heck 1979). Abnormally high pore pressures in the Bakken shales have been attributed to hydrocarbon generation associated with the thermal anomaly in some parts of the basin (Price et al. 1984). Hydrocarbon production in the Bakken depends on accurate placement of horizontal wells with multistage fracture stimulations.

**Fig. 3** Stratigraphic column of study area.

The effective minimum horizontal stress is a primary controlling factor for fracture growth. Thus, knowledge of distribution of elastic properties of the Bakken shales is very important to accurately determine the effective minimum horizontal stress (Havens 2012). In terms of wireline log data interpretation, Bakken shale is easily recognizable because of the strong contrast in lithology (Havens 2012). The upper and lower shales have unusually high gamma ray readings and high resistivity, while the middle member has a signature similar to clastic and carbonate rocks.

## Data collection and conditioning

For the current study, geologic data on Hoosier field were sourced and collected from public domain sources (Alberta Geological Survey; Government of Saskatchewan). The data gathered for the current study mainly comprised of well location information, interpreted formation tops demarcating geological boundaries, wireline logs, geological maps, etc. No seismic data are available covering the study area to enable understanding of lateral continuity of lithological properties within reservoir layers. Thus, facies logs interpreted and derived from available wireline log data were heavily relied upon for establishing spatial continuity of depositional model of the reservoir. Three main geological markers (formation tops) representing three most important geological units, i.e., Mannville, Detrital, and Bakken, were identified from the wireline log interpretations of 27 well data. Figure 3 shows a stratigraphic column with the formations of interest.

## 3D geocellular modeling

The interpreted formation tops namely Mannville, Detrital, and Bakken were first used to generate depth structure grids (i.e., maps). These depth structure grids (Fig. 4a) serve as input to create the structural framework of the reservoir with top and base of the framework as Mannville and Bakken formations. Internal stratigraphic unit thicknesses and their relationship to each other helped in establishing the layering schemes in the structural framework. The structural framework constitutes the basic grid skeleton for discretization of reservoir properties. The 3D structural framework is then divided into blocks of size 150 m × 150 m with vertical layer thickness of 0.5 m (Fig. 4b, d).

**Fig. 4** a Depth structural map from tops, b 3D grid layers, c Location of wells in the study area and d 3D block configuration.

A detailed multi-well lithostratigraphic correlation was carried out to establish the continuity of subsurface lithology. Multivariate statistical analysis technique is employed to group the log curves into electrofacies groups. Each of the electrofacies classes has been assigned a rock class based on the sedimentological interpretation of the wells. Lithological classes thus assigned are correlated with each of the wells covering the entire study area. The assigned lithology classes formed the percentage of rock values (Fig. 4c) that would be used to fill in the grid during facies modeling. The well-wise rock interpretation is next combined with petrophysical information to establish its effect on porosity, permeability, and saturation distribution controlling dynamic fluid behavior.

Log measurement data from 27 wells in terms of petrophysical properties are converted to point sets and then blocked to the grid. A point set is collections of generic points that has 3D spatial location information, i.e., easting (x), northing (y), and depth (z) and rock property measurement values such a porosity and permeability corresponding to each point locations. While converting a log curve to point set, first measured depth (MD) (z) values and corresponding rock property measurements from log curves (such as porosity) are captured in columns in a data file. Next, easting (x) and northing (y) information is captured from well head and added to each data points as two separate columns. Thus, the resultant point set data file has series of columns starting with x, y, z and followed by rock property values. Generally, each of these data points is taken at the scale of actual log measurement, i.e., 0.5 feet (approximately 15 cm). After creation of the point set, the next step is to block the point set to the grid resolution (Shepherd 2009).

Usually, the vertical dimensions of geocellular grids are larger than the vertical sampling interval of log curves (or the point set). As a result, each geocells passing through well locations carry multiple points. In order to use the point information for modeling and interpolation, the multiple points within each cell need to convert to one point per grid cells using an averaging technique such as arithmetic, geometric, and harmonic for continuous properties. In case of discrete properties (such as facies or lithology), the averaging approach is “most of,” i.e., most commonly occurring (mode) lithology within the points occurring in a single geocells (Cannon 2015). This process is referred as blocking of point sets to the geocellular grid. Lithotype assignment and vertical proportion curves in terms of three main lithotypes, i.e., dolomitic sandstone, siltstone, and shale are generated. These lithotype proportion curves resulted into lithotype proportion maps (Fig. 5) which are then incorporated in generation of facies distribution model for the reservoir.

**Fig. 5** Vertical proportion curve and lithotype proportion map showing vertical and areal distribution of rock types.

Variogram models have been constructed to establish data stationarity and spatial correlation of input data points over the study area. Establishing stationarity is expedient for an implicit assumption of a geostatistical calculation in the current study. The absence of a nugget effect in the variograms showed that there is no randomness in the sampled well population. The spatial position of well locations produced a large spatial correlation distance for each of the variables. This correlation allows to progress with a geostatistical approach to extrapolate the measured data for the entire grid cells.

Plurigaussian simulation (PGS) algorithm (Armstrong et al. 2003) was employed with two variogram models for creating 3D lithofacies distribution model. In PGS simulation, the controlling component is lithotype rule, which defines the depositional environment (Armstrong et al. 2003). In addition to lithotype rules, the lithotype proportion maps resulted from lithotype curves provided a background guideline for cell-wise facies distribution. PGS is one of the best techniques in geologically complex depositional environment modeling. Next, sequential Gaussian simulation (SGS) (Deutsch and Journel 1998; Deutsch 2002) was used for creating 3D spatial distribution model of petrophysical properties (such as porosity and permeability) employing variograms created for each of these properties. While creating the porosity model, previously created facies model was used to constrain or bias reducing geological uncertainty of modeling of petrophysical properties. Distribution of facies and associated porosity in the reservoir is shown in Fig. 6.

**Fig. 6** Distribution of facies (rock types) and porosity (facies-constrained).

## Building 1D basin model

The concept of basin modeling in general refers to integration of geological, physical, and geochemical processes in a sedimentary basin to simulate sediment burial, compaction, heat transfer, hydrocarbon generation, and migration and entrapment processes within it. Basin modeling reconstructs the pressure regime of the basin during sediment deposition and predicts zone of overpressure. This process also delineates hydrocarbon migration pathways and identifies structural and stratigraphic traps.

Basin modeling can be performed in 1D (one-dimensional), 2D (two-dimensional) or 3D (three-dimensional) to simulate the burial and thermal histories of the basin with increasing complexities. If the area is relatively small, such that thermal, pressure, and other effects are uniform over the study area, then a 1D basin model might be sufficient (Yarus and Carruthers 2014). If the basin development processes vary significantly over the area, a 3D study is warranted. In the current study, the necessity of 2D/3D basin model is not felt due to relatively uniform thermal and pressure effect on the basin-forming process in the study area. Comparison 1D basin model results from few places within the study area also revealed that important parameters controlling the basin-forming process do not vary much. Thus, the current study was restricted to construction of 1D basin model and uses the results for derivation of geomechanical properties.

The 1D basin model uses rock and fluid information to evaluate source rocks and thermal histories of basin along a single direction, i.e., vertical. In this process, layers in the 1D model are sequentially back-striped to its original depositional thickness resulting into effective stress–porosity relationships for all lithologic layers. The lithologies are generally assigned to the model layers by zones using petrophysical logs and other geological information. These lithologies provide thermal conductivity and heat capacity for each layer.

In the current study, the input parameters used in developing 1D basin model of the reservoir include:

- Geologic framework in terms of layer thicknesses and ages
- Rock properties (i.e., lithology) for each framework layers
- Well data information such as temperature and pressure
- Identification of source layers with assignment of kerogen types
- Regional geothermal gradient of the basin

Figure 7 shows the output of basin modeling process in terms of variation of elastic properties and pressure regim while depositing the sediments over geological periods. Modeling of the parameters is governed by the rates of pressure generation (e.g., compaction disequilibrium) and pressure dissipation (e.g., as controlled by permeability) during basin-forming process. Another important output of 1D basin modeling process is burial history or geohistory plot of the reservoir as shown in Fig. 8. Geohistory plots indicate how sediments were deposited over geological time in the basin providing thickness variation information of layers.