An integrated approach to discretized 3D modeling of geomechanical properties for unconventional mature field appraisal in the western Canadian sedimentary basin


ORIGINAL PAPER - EXPLORATION GEOLOGY

An integrated approach to discretized 3D modeling
of geomechanical properties for unconventional mature field
appraisal in the western Canadian sedimentary basin

Kalyan Saikia1 • Clara E. Ikuku2 • Bhabesh C. Sarkar3

Received: 14 September 2016 / Accepted: 20 November 2017 / Published online: 30 November 2017

� The Author(s) 2017. This article is an open access publication

Abstract In mature field appraisal and development, dis-

cretized 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 geome-

chanical models. However, elastic properties measured this

way are restricted to the well location and interpolation

across the reservoir may not be always appropriate. To

overcome these challenges, this paper describes an inte-

grated 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 dis-

cretization 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 inno-

vative way of determining effective horizontal stress for

the entire reservoir through 3D distribution of elastic

properties.

Keywords Reservoir � Geomechanics � Basin model �
Unconventional � Geostatistics

Introduction

In unconventional tight reservoirs, where hydraulic frac-

turing is the key technique for enhanced oil production, a

thorough knowledge of distribution of geomechanical

properties helps to maximize return on investment. Hori-

zontal 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 min-

imum horizontal stress as it guides the stimulation design

which will in turn control the fracture propagation. Accu-

racy of effective minimum horizontal stress model is pri-

marily 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 con-

ventional laboratory methods which may not be repre-

sentative of the entire elastic regime existing in the

& Kalyan Saikia
kalyansaikia@gmail.com

1
Halliburton, 24335 Piazza Drive, Richmond, TX 77406, USA

2
Clenik Petrotech, 104 Tremblant Way SW, Calgary, AB

T3H0A7, Canada

3
Department of Applied Geology, Indian Institute of

Technology (Indian School of Mines), Dhanbad 826004,

India

123

J Petrol Explor Prod Technol (2018) 8:417–429

https://doi.org/10.1007/s13202-017-0406-3

http://orcid.org/0000-0002-8500-9690
http://crossmark.crossref.org/dialog/?doi=10.1007/s13202-017-0406-3&domain=pdf
http://crossmark.crossref.org/dialog/?doi=10.1007/s13202-017-0406-3&domain=pdf
https://doi.org/10.1007/s13202-017-0406-3


reservoir as in most cases these small samples are taken

from consolidated rocks and used to measure geome-

chanical 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 geome-

chanical 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 suffi-

cient enough to describe spatial distribution pattern of

elastic properties for entire reservoir. Havens (2012) car-

ried out similar approach for Bakken petroleum system to

establish the spatial behavior elastic properties and con-

cluded 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 geom-

etry 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 geome-

chanical 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 proper-

ties of the reservoir rock, explaining their response to

mechanical stresses at reservoir condition. Such relation-

ship is very important for defining trends in the reservoir

description. The current approach also incorporates wire-

line 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 geome-

chanical 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 geocel-

lular 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 lithofa-

cies within the reservoir layers employing advanced geo-

statistical algorithms. Lastly, spatial distribution model of

petrophysical property in terms of porosity, permeability,

saturation etc., is generated and constrained to the litho-

facies model. This helps in capturing geological uncer-

tainty 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 continu-

ous 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 vari-

ance 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 ran-

dom variables Z(ua), a = 1,…, N defined at N locations ua.
The aim is to generate L joint realizations {z(l)(ua),_ = 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 con-

ditional distributions:

418 J Petrol Explor Prod Technol (2018) 8:417–429

123



F(u1, …, uN; z1, …, zN|(n)) = F (uN; zN|(n ?
N - 1)) 9 F (uN - 1; zN - 1|(n ? N - 2)) 9 ��� 9
F(u2; z2|(n ? 1)) 9 F(u1; z1|(n)),

Where F (uN; zN|(n ? N - 1)) = Prob

{Z(uN) B 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)(ua), a = 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 analyti-

cally available (Deutsch 2002). The conditional probability

distribution in a Gaussian model is determined using sim-

ple 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 conve-

niently used in SGS making it one of the popular algo-

rithms 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 archi-

tecture 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 recon-

struction, 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 under-

stand effect of physical and chemical processes in a sedi-

mentary 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 spa-

tial grids. This computation results in many outputs of

different rock properties including properties such as vit-

rinite reflectance, temperature, effective stress, and poros-

ity. 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

J Petrol Explor Prod Technol (2018) 8:417–429 419

123



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, tempera-

ture 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 prob-

lems, 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 reser-

voir 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 relation-

ships 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.

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 Sas-

katchewan, 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 forma-

tion 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). 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 con-

sisting 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 pro-

duction in the Bakken depends on accurate placement of

horizontal wells with multistage fracture stimulations. The

effective minimum horizontal stress is a primary control-

ling factor for fracture growth. Thus, knowledge of distri-

bution 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 (Al-

berta 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, geo-

logical maps, etc. No seismic data are available covering

the study area to enable understanding of lateral conti-

nuity 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) repre-

senting three most important geological units, i.e., Man-

nville, 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.

420 J Petrol Explor Prod Technol (2018) 8:417–429

123



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

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

J Petrol Explor Prod Technol (2018) 8:417–429 421

123



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 thick-

nesses and their relationship to each other helped in

establishing the layering schemes in the structural frame-

work. 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 9 150 m with vertical layer thickness of 0.5 m

(Fig. 4b, d).

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 cap-

tured 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 geo-

cellular 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 inter-

polation, 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.

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

Fig. 3 Stratigraphic column of study area

422 J Petrol Explor Prod Technol (2018) 8:417–429

123



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 simu-

lation, 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

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

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

J Petrol Explor Prod Technol (2018) 8:417–429 423

123



these properties. While creating the porosity model, pre-

viously created facies model was used to constrain or bias

reducing geological uncertainty of modeling of petro-

physical properties. Distribution of facies and associated

porosity in the reservoir is shown in Fig. 6.

Building 1D basin model

The concept of basin modeling in general refers to inte-

gration of geological, physical, and geochemical processes

in a sedimentary basin to simulate sediment burial, com-

paction, heat transfer, hydrocarbon generation, and

migration and entrapment processes within it. Basin mod-

eling reconstructs the pressure regime of the basin during

sediment deposition and predicts zone of overpressure.

This process also delineates hydrocarbon migration path-

ways 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 con-

struction 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 infor-

mation. 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. In other words, it shows the time versus thickness

relationship for all the layers including the layers that have

been eroded or periods of non-deposition (Hiatus). Based

on the lithology information of the 1D model inputs,

depth–porosity relationship is also calculated while estab-

lishing burial history. The depth–porosity relationship is

then converted to calculate the effective stress regime and

temperature profile of each layer through the geological

time. The effective stress thus derived for each layer of the

basin provides crucial information for calculation of

geomechanical properties of the reservoir rock.

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

424 J Petrol Explor Prod Technol (2018) 8:417–429

123



Fig. 7 Output of 1D basin modeling process in terms of reservoir elastic properties and pressure regim

Fig. 8 Geohistory plot (burial history)—output from 1D basin modeling process

J Petrol Explor Prod Technol (2018) 8:417–429 425

123



Model integration

The final stage of the current integrated approach is to

combine the result of classical 1D basin model (burial

history) with geostatistical 3D geocellular model to gen-

erate discretized geomechanical model of the study area.

Lithology identifier of each cell from 3D facies model

provides the link of this integration. Each element of burial

history model has also been assigned a lithological iden-

tifier to define (1) mechanical compaction, i.e., porosity

versus depth relationships, (2) pore pressure relationships,

i.e., porosity versus effective stress, (3) fluid transport

properties, i.e., porosity versus permeability and capillary

threshold pressure, and (4) matrix properties such as ther-

mal conductivity. The lithology identifiers available in both

the models enable populating layer-wise distribution of

effective stress for each cell in the 3D geocellular grid

(Fig. 9). The effective stress estimation is obtained during

the process of burial history modeling of the sedimentary

basin. The burial (or geo) history model provides infor-

mation on rate of sedimentation including depth and age of

deposition. While calculating effective stress regime of the

basin, we assume that sedimentation and subsidence is the

primary cause of pressure and stress; stresses generally

increase with depth and rock stress generally interacts

negatively with porosity and compaction (Hantschel and

Kauerauf 2009). In 1D basin modeling approach, Terza-

ghi’s (1923) effective stress concept is widely used as the

fundamental relationship for effective stress modeling

(Hubert and Rubey 1959).

Effective stress ðr0Þ ¼ Total stress ðrÞ�Pore pressure pð Þ:
ð1Þ

In this relationship, lithostatic pressure (equaling the

sediment weight) is considered as total stress and

additional stresses due to compaction or extensional

forces are neglected (Hantschel and Kauerauf 2009). Biot

(1941) established the theory of poro-elasticity for rocks

and introduced Biot’s coefficient in the definition of the

effective stress (Schneider et al. 1996).

r0 ¼ r � a p; where a is Biot’s coefficient: ð2Þ

In these models, pore pressure of formation is related to

incomplete mechanical compaction and a fixed relation

between porosity reduction and sediment compaction is

assumed (Hantschel and Kauerauf 2009). Studies by

various other researchers established a number of

empirical relationships of porosity versus depth versus

compaction. These relationships show that changes in

porosity with increasing depth is a function of lithology of

the sediments. Thus, in basin modeling process these

lithology-wise empirical relationships are used to satisfy

the depth–porosity function and pore pressure calculation.

The complete mathematical explanations of these

relationships and numerical solutions to complex

equations can be obtained from Hantschel and Kauerauf

(2009) and Schneider et al. (1996).

Once the effective stress distribution model is generated

for the entire reservoir geocellular grid, rock mechanical

properties such as Vp, Vs, UCS, BRI, OCR, Young’s

modulus, Poison’s ratio, and bulk modulus are calculated at

each cell locations of geocellular grid using empirical

formula as provided below (Eberhert-Phillips et al. 1989;

Ingram and Urai 1999; Horsrud 2001; Yarus and Car-

ruthers 2014):

Vp ¼ 5:77 � 6:94/ � 1:73
p
C þ 0:446 Pe � e�16:7Pe

� �
;

ð3Þ

Vs ¼ 3:70 � 4:94/ � 1:57
p
C þ 0:361 Pe � e�16:7Pe

� �
;

ð4Þ
Log UCS ¼ �6:36 þ 2:45 log 0:86 Vp � 1172ð Þ; ð5Þ
OCR ¼ Max Pe=Pe; ð6Þ
BRIT ¼ UCS=UCSNC; ð7Þ

E ¼ qV2s 3V
2
p � 4V

2
s

� �
= V2p � V

2
s

� �h i
; ð8Þ

m ¼ 0:5 V2p � 2V
2
s

� �
= V2p � V

2
s

� �h i
; ð9Þ

K ¼ q V2p � 4=3V
2
s

� �
; ð10Þ

where Vp = P-wave velocity, Vs = S-wave velocity, /
= Porosity, C = Clay content, Pe = Effective stress,

UCS = Unconfined compressive strength, OCR = Over

consolidation ratio, BRIT = Brittleness, NC = Uncon-

fined, q = Density, E = Young’s modulus, m = Poisson’s
ratio, K = Bulk modulus.

For each of the calculations lithology, porosity and

effective stress serve as main input parameter. Figure 9

displays a representative view of the geocellular grid with

discretized geomechanical property derived from the inte-

gration of result of basin and geocellular modeling.

Results and discussion

Well coverage in the study area is not adequate as revealed

in Fig. 4c. The wells are mostly clustered in the central and

north-eastern part. Furthermore, in the study area conven-

tional cores were collected in only five wells out of 27

wells made available. Thus, physical measurement and

interpolation of geomechanical properties of these core

samples were not considered as reliable source of infor-

mation for future drilling and investment decision. Addi-

tionally, due to non-availability of 3D seismic data,

426 J Petrol Explor Prod Technol (2018) 8:417–429

123



interpolation of well-based geomechanical properties

beyond well control could not be validated through seismic

attribute correlation and integration as secondary input. As

a result, attempts to use the conventional approach of using

well core sample analysis-based geomechanical property

interpolation in the study area provided inappropriate

information for full field appraisal and development plan-

ning. However, use of current geostatistics-based modeling

approach explained in this study helped to nullify these

restrictions and provided a full field geomechanical prop-

erty distribution model with the integration of parameters

from basin modeling approach. Such integrated models are

highly flexible and appropriate for full field development

with reliable distribution of geomechanical property in 3D

space of unconventional reservoir.

The model developed in the current study clearly indi-

cates that elastic moduli decrease with increasing porosity;

therefore, such models would aid in the understanding of

the distribution of brittle versus ductile zones in the

reservoir. Statistical analysis in terms of correlation

coefficient of Young’s modulus and porosity reveals a

value of - 0.86. Figure 10 displays a plot of Young’s

modulus and porosity at K-level 23 of 3D reservoir grid.

Fig. 9 Distribution of
geomechanical properties in the

final integrated reservoir grid

Fig. 10 Plot showing correlation of Young’s modulus and porosity
for K-level 23 of 3D reservoir grid

J Petrol Explor Prod Technol (2018) 8:417–429 427

123



Based on the spatial distribution pattern of geomechanical

properties, the lower part of the Bakken formation in the

study area has been identified as zone of prime importance

with favorable values of geomechanical properties to

maximize hydrocarbon production. In the absence of 3D

seismic coverage in the study area, such information aids

operators tremendously for effective design of depletion

strategies with minimal geological risks and associated

costs. Furthermore, for a variety of other subsurface

operations, including flow simulation, history matching,

fluid typing, and hydrocarbon production, these models

serve as main input.

Understanding the burial history of the reservoir is a step

forward in modeling unconventional resources (Yarus and

Carruthers 2014). As part of full field development project,

a realistic picture of the depositional history of the reser-

voir through geological time provides a more solid foun-

dation to build a knowledge base for economic forecasting

and maximum return of investment. In the unconventional

upstream oil and gas industry, spatial distribution modeling

of rock mechanical properties of matured reservoir has

always been a major challenge. However, 3D spatial dis-

tribution of model of geomechanical properties created in

the current study provides an easy way of knowing pre-

vailing rock elastic properties and stress regimes of the

reservoir rocks.

Conclusions

Field development planning, in the current study area,

involves the use of horizontal wells with multistage

hydraulic stimulation in shale sweet spots as the primary

production strategies. Knowledge of in situ stress condi-

tions and facies is therefore very important for these

operations as it guides the stimulation design which will in

turn control the fracture propagation. Earlier studies on

Bakken formation clearly elaborated the use of velocities

and strains measured on samples from the core which are

not sufficient enough to optimize field development strat-

egy using hydraulic fracturing. This is because small

sample plugs or slabs may not be always sufficient to

describe elastic properties for the entire reservoir. For

instance, Vp measured from one core sample may not be

very diagnostic of lithology as it cannot be used to map

shales and sandstones across a reservoir with accuracy.

The current integrated approach incorporates wireline

logs in the form of hard data covering entire sediment

sequence. This provides a true representation of the sub-

surface geological condition of the reservoir even if there

are no 3D seismic data available. In this approach, use of

geostatistical techniques results into highly accurate 3D

facies model resembling real geological condition of

deposition of sediments. The geostatistical techniques via

spatial covariance models of sample values helped in

quantifying and mapping of geological heterogeneity in

sediment deposition in unconventional reservoirs. Addi-

tionally, conditional simulation algorithms used in geo-

statistical modeling approach honor the sample values at

location of measurements and produce minimum biasness

of interpolation of rock properties from well data. Thus,

results generated from final integrated model are not only

reliable, but also capable of providing accurate local error

of estimation (i.e., heterogeneity) within them. Similar

studies made elsewhere (Yarus et al. 2016; Yarus and

Yarus 2014; Yarus and Carruthers 2014) indicated how

combining geostatistics with basin modeling techniques

accurately identifies sweet spots in unconventional reser-

voirs. Verification and comparison of results obtained from

such integrated models with actual and blind wells pro-

vided high accuracy and reliability score with minimum

deviation of results. Such convincing results are also true

validation of consistency of result of current integrated

technique over conventional approach.

Traditional modeling methods for identifying good and

poor reservoir quality in shales have not been successful as

reported by many authors. It is mainly due to lack of

geospatial understanding of shales properties and difficulty

of quantifying their heterogeneity. However, the current

integrated geostatistics-based geomechanical models with

higher degree of confidence show how the velocity and

elastic moduli are decreasing with increasing porosity

throughout the reservoir. As a result, in the current study

area the concerned field team was able to identify target

locations or sweet spots with high degree of accuracy

within the shale for fracturing and completion. Further-

more, these models can be used alongside other reservoir

characterization methods in shale play workflows to aid in

the planning of unconventional well bores for enhanced

production.

Open Access This article is distributed under the terms of the
Creative Commons Attribution 4.0 International License (http://

creativecommons.org/licenses/by/4.0/), which permits unrestricted

use, distribution, and reproduction in any medium, provided you give

appropriate credit to the original author(s) and the source, provide a

link to the Creative Commons license, and indicate if changes were

made.

References

Alberta Geological Survey. http://ags.aer.ca/. Accessed July 2016

AlKawai W (2014) Integrated basin modeling with seismic technol-

ogy and rock physics. M.Sc. Thesis, Stanford University,

Stanford

Allard D, D’Or D, Biver P, Froidevaux R (2012) Non-parametric

diagrams for plurigaussian simulations of lithofacies. In: 9th

international geostatistical congress, Oslo, Norway, pp 11–15

428 J Petrol Explor Prod Technol (2018) 8:417–429

123

http://creativecommons.org/licenses/by/4.0/
http://creativecommons.org/licenses/by/4.0/
http://ags.aer.ca/


Armstrong A, Galli AG, Loch GL, Geoffroy F, Eschard R (2003)

Plurigaussian simulations in geosciences. Springer, Berlin, p 144

Arpat GB (2005) Sequential simulation with patterns. PhD Thesis.

Stanford University, Stanford

Biot MA (1941) General theory of three-dimensional consolidation.

J Appl Phys 12:155–164

Cannon S (2015) Petrophysics: a practical guide. Wiley, London,

pp 16–161

Deng H, Aguilera R, Settari T (2011) An integrated workflow for

reservoir modeling and flow simulation of the nikanassin tight

gas reservoir in the western Canada sedimentary basin. SPE

Annual Technical Conference and Exhibition, Denver, Colorado,

USA. SPE # 146953-MS. doi:10.2118/146953-MS

Deutsch CV (2002) Geostatistical reservoir modeling. Oxford

University Press, Oxford

Deutsch CV, Journel AG (1998) GSLIB: geostatistical software

library and user’s guide, 2nd edn. Oxford University Press,

Oxford

Dimitrakopoulos R, Luo X (2004) Generalized sequential gaussian

simulation on group size m and screen-effect approximations for
large field simulations. Math Geol 36(5):567–591

Eberhert-Phillips D, Han D-H, Zoback MD (1989) Empirical

relationships among seismic velocity, effective pressure, poros-

ity, and clay content in sandstone. J Geophys 54(1):82–89

Ferdous H, Bergman K, Qing H (2015) Viking lowstand deposits in

west central Saskatchewan: depositional model for the reservoir

units in Dodsland-Hoosier Area, Saskatchewan, Canada. AAPG

Search and Discovery Article #51076

Galli A, Beucher H, Le Loc’h G, Doligez B (1994) The pros and cons

of the truncated gaussian method. In: Armstrong M, Dowd PA

(eds) Geostatistical Simulations. Kluwer, Dordrecht, pp 217–233

Gonzalez L, Nasreldin G, Rivero J, Welsh P, Aguilera R (2014) 3D

modeling of multistage hydraulic fractures and two-way-cou-

pling geomechanics/fluid-flow simulation of a horizontal well in

the nikanassin tight gas formation, Western Canada sedimentary

basin. SPE Reserv Eval Eng 17(2):257–270

Government of Saskatchewan. http://www.ir.gov.sk.ca/ Accessed

July 2016

Hantschel T, Kauerauf AI (2009) Fundamentals of basin and

petroleum systems modeling. Springer, Berlin

Havens J (2012) Mechanical properties of the Bakken Formation.

M.S. Thesis, Colorado School of Mines, Golden

Heck TJ (1979) Depositional environments and diagenesis of the

Mississippian Bottineau interval (Lodgepole) in North Dakota:

University of North Dakota Unpublished M.S. Thesis

Horsrud P (2001) Estimating mechanical properties of shale from

empirical correlations. SPE Drilling and Completion conference,

pp 68–73 (SPE # 00056017)

Hubert MK, Rubey WM (1959) Role of fluid pressure in mechanics of

overthrust faulting. Bull Geol Soc Am 70:115–166

Ingram GM, Urai JL (1999) Top-seal leakage through faults and

fractures: the role of mudrock properties. Geol Soc Lond Spec

Publ 158:125–135

Lavoie D, Séjourné S (2016) Geomechanical properties of the upper

ordovician Macasty shale and its caprock, Anticosti Island: a

regional evaluation for a promising tight oil play. AAPG Search

and Discovery Article #80542

Matheron G, Beucher H, de Fouquet C, Galli A, Guerulot D, Ravenne

C (1987) Conditional simulation of the geometry of fluvio-

deltaic reservoirs. SPE annual technical conference and exhibi-

tion, Dallas, Texas. SPE# 16753

Meissner FF (1978) Petroleum geology of the Bakken Formation,

Williston Basin, North Dakota and Montana. In: Proceedings of

1978 Williston basin symposium, Montana Geological Society,

pp 207–227

O’Connell SC, Wood JM, Robinson SP (2000) An ancient tidal sand

ridge complex in the Bakken Formation (Fammenian-Tour-

naisian) of southwestern Saskatchewan. GeoCanada 2000 core

workshop, Calgary, Alberta. http://webapp1.dlib.indiana.edu/

virtual_disk_library/index.cgi/2870166/FID3366/PDF/1228.PDF

Peters KE (2009) Getting started in Basin and petroleum system

modeling. AAPG CD-ROM #16, AAPG Datapages

Porter JW (1955) Madison complex of Southeastern Saskatchewan

and Southwestern Manitoba. J Alberta Soc Petrol Geol

3(8):126–130

Press WH, Teukolsky SA, Vetterling WT, Flannery BP (2002)

Numerical recipes in C??, 2nd edn. Cambridge University

Press, Cambridge

Price LC, Ging T, Daws T, Love A, Pawlewicz M, Anders D (1984)

Organic metamorphism in the Mississippian-Devonian Bakken

shale, North Dakota portion of the Williston Basin. In: Wood-

ward J, Meisner FF, Clayton JL (eds) Hydrocarbon source rocks

of the greater Rocky Mountain region. Rocky Mountains

Association-Geologists, Denver, pp 83–133

Schneider F, Potdevin JL, Wolf S, Faille I (1996) Mechanical and

chemical compaction model for sedimentary basin simulators.

Tectonophysics 263:307–317

Shepherd M (2009) Oil field production geology: AAPG Memoir 91,

pp 182–183

Szydlik T, Helgesen HK, Brevik I, Prisco GD, Clark SA, Leirfall OK,

Duffaut K, Stadtler C, Cogan M (2015) Geophysical basin

modeling: methodology and application in deepwater Gulf of

Mexico. Interpretation 3(3):SZ49–SZ58

Terzaghi K (1923) Die Berechnung der durchlssigkeitsziffer des tones

aus dem verlauf des hydrodynamishen spannungserscheinungen.

Szgber Akad Wiss Vieu Math Naturwiss Klasse Ila

132(3–4):125–138

Tong A, Chi G, Pedersen PK (2005) Sequence stratigraphy and

preliminary diagenetic study of the lower cretaceous viking

formation, Hoosier Area, West-Central Saskatchewan. Summ

Investig Sask Geol Surv 1:2005

Wheeler HE (1963) Post-Sauk and Pre-Absaroka Paleozoic strati-

graphic patterns in North America. AAPG Bull 46:1497–1526

White WI (1969) Geology and petroleum accumulations of the North

Hoosier area, west central Saskatchewan. Department of Mineral

Resources, Geological Science Branch, Sedimentary Geology

Division, Report No. 133

Yarus JM, Carruthers D (2014) Modeling unconventional resources

with geostatistics and basin modeling techniques: characterizing

sweet-spots using petrophysical, mechanical, and static fluid

properties. AAPG Search and Discovery Article #41362

Yarus JM, Yarus JM (2014) Modeling with geostatistics, basin

modeling techniques, E & P Magazine. http://www.epmag.com/

modeling-geostatistics-basin-modeling-techniques-758646#p=

full

Yarus JM, Mohsenian E, Yarus JM (2016) Combing geostatistical

reservoir modeling with classical basin modeling techniques to

identify sweet-spots in unconventional reservoirs. AAPG Data-

pages/Search and Discovery Article #90254

Zhang L, Buatois LA, Knaust D (2016) Sedimentology, ichnology

and sequence stratigraphy of the upper devonian-lower Missis-

sippian Bakken formation in eastern Saskatchewan. Bull Can Pet

Geol 64(3):415–437

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in

published maps and institutional affiliations.

J Petrol Explor Prod Technol (2018) 8:417–429 429

123

https://doi.org/10.2118/146953-MS
http://www.ir.gov.sk.ca/
http://webapp1.dlib.indiana.edu/virtual_disk_library/index.cgi/2870166/FID3366/PDF/1228.PDF
http://webapp1.dlib.indiana.edu/virtual_disk_library/index.cgi/2870166/FID3366/PDF/1228.PDF
http://www.epmag.com/modeling-geostatistics-basin-modeling-techniques-758646%23p%3dfull
http://www.epmag.com/modeling-geostatistics-basin-modeling-techniques-758646%23p%3dfull
http://www.epmag.com/modeling-geostatistics-basin-modeling-techniques-758646%23p%3dfull

	An integrated approach to discretized 3D modeling of geomechanical properties for unconventional mature field appraisal in the western Canadian sedimentary basin
	Abstract
	Introduction
	Integrated modeling approach
	Step 1
	Step 2
	Step 3
	Geological background of the study area
	Data collection and conditioning

	3D geocellular modeling
	Building 1D basin model
	Model integration

	Results and discussion
	Conclusions
	Open Access
	References