Language selection

Search

Patent 2735038 Summary

Third-party information liability

Some of the information on this Web page has been provided by external sources. The Government of Canada is not responsible for the accuracy, reliability or currency of the information supplied by external sources. Users wishing to rely upon this information should consult directly with the source of the information. Content provided by external sources is not subject to official languages, privacy and accessibility requirements.

Claims and Abstract availability

Any discrepancies in the text and image of the Claims and Abstract are due to differing posting times. Text of the Claims and Abstract are posted:

  • At the time the application is open to public inspection;
  • At the time of issue of the patent (grant).
(12) Patent: (11) CA 2735038
(54) English Title: STRESS AND FRACTURE MODELING USING THE PRINCIPLE OF SUPERPOSITION
(54) French Title: MODELISATION DE CONTRAINTE ET DE FRACTURE UTILISANT LE PRINCIPE DE LA SUPERPOSITION
Status: Granted and Issued
Bibliographic Data
(51) International Patent Classification (IPC):
  • G01V 09/00 (2006.01)
  • G06F 17/10 (2006.01)
(72) Inventors :
  • MAERTEN, FRANTZ (France)
  • MAERTEN, LAURENT (France)
(73) Owners :
  • SCHLUMBERGER CANADA LIMITED
(71) Applicants :
  • SCHLUMBERGER CANADA LIMITED (Canada)
(74) Agent: SMART & BIGGAR LP
(74) Associate agent:
(45) Issued: 2014-08-19
(22) Filed Date: 2011-03-23
(41) Open to Public Inspection: 2011-09-25
Examination requested: 2011-03-23
Availability of licence: N/A
Dedicated to the Public: N/A
(25) Language of filing: English

Patent Cooperation Treaty (PCT): No

(30) Application Priority Data:
Application No. Country/Territory Date
13/052,327 (United States of America) 2011-03-21
61/317,412 (United States of America) 2010-03-25

Abstracts

English Abstract

Stress and fracture modeling using the principal of superposition is provided. A system simulates linearly independent far field stress models for a subsurface earth volume, computing stress, strain, and displacement values based on superposition of independent stress tensors. Based on the precomputed values, the system generates real-time recovery of paleostress values, or, stress, strain, and displacement parameters for any point in the subsurface volume as the user varies far field stress values. The system recovers one or more tectonic events, or a stress tensor represented by a ratio of principal magnitudes and associated orientation, using fault geometry, well bore data (fracture orientation and secondary fault plane data), GPS, InSAR, folded and faulted horizons, tiltmeters, slip and slikenlines on faults. The system uses different geologic data from seismic interpretation, well bore readings, and field observation to provide numerous results, such as predicted fracture propagation based on perturbed stress field.


French Abstract

La modélisation de contrainte et de fracture utilisant le principe de la superposition est présentée. Un système simule des modèles de contrainte en zone éloignée indépendants linéairement pour un volume terrestre souterrain, une contrainte calculée, la pression et les valeurs de déplacement fondées sur la superposition des tenseurs de contrainte indépendants. En fonction des valeurs précalculées, le système génère la récupération en temps réel des valeurs de paléocontrainte, ou les paramètres de contrainte, de déformation et de déplacement de tout point dans le volume souterrain au moment où l'utilisateur varie les valeurs de contrainte en zone éloignée. Le système récupère un ou plusieurs événements tectoniques ou un tenseur de contrainte représenté par un rapport entre les magnitudes principales et l'orientation associée, à l'aide de la géométrie de faille, des données du puits de forage (orientation de la fracture et données de plan de faille secondaire), des données GPS, des données InSAR, des horizons de plis et de failles, des inclinomètres, du glissement et des lignes de batture sur les failles. Le système emploie différentes données géologiques issues de l'interprétation sismique, des lectures de puits de forage et d'observation sur le terrain pour fournir plusieurs résultats, comme la propagation de fracture prédite fondée sur le champ de pression perturbé.

Claims

Note: Claims are shown in the official language in which they were submitted.


CLAIMS:
1. A computer-readable storage medium storing an executable program,
the executable program containing instructions, which when
executed by a computer, perform a process, comprising:
simulating three linearly independent far field stress models for a
subsurface earth volume;
computing stress, strain, or displacement values for data points in the
subsurface earth volume, based on a principle of superposition applied to the
three linearly independent far field stress models; and
predicting a stress attribute of the subsurface earth volume in real-time,
based, on the precomputed stress, strain, and displacement values.
2. The computer-readable storage medium as recited in claim 1, wherein
the three linearly independent
far field stress models are based on different data sets, each data set
capable
of being weighted differently and each data set comprising one of fault
geometry data, fracture orientation data, stylolites orientation data,
secondary
fault plane data, fault throw data, slickenline data, global positioning
system
(GPS) data, interferometric synthetic aperture radar (InSAR) data, laser
ranging data, tilt-meter data, displacement data for a geologic fault, or
stress
magnitude data for the geologic fault.
3. The computer-readable
storage medium as recited in claim 1, wherein the predicted stress attribute
comprises one of a stress inversion, a stress field, a far field stress value,
a
stress interpolation in a complex faulted reservoir, a perturbed stress field,
a
stress ratio and associated orientation, one or more tectonic events, a
displacement discontinuity of a fault, a fault slip, an estimated
displacement, a
perturbed strain, a slip distribution on faults, quality control on
interpreted
38

faults, fracture prediction, prediction of fracture propagation according to
perturbed stress field, real-time computation of perturbed stress and
displacement fields while performing interactive parameters estimation, or
discernment of an induced fracture from a preexisting fracture.
4. A computer-executable method, comprising:
simulating linearly independent far field stress tensors in constant time
based on multiple types of geologic data associated with a subsurface
volume;
computing precomputed strain, stress, or displacement values based on
a principle of superposition applied to the simulations; and
iteratively minimizing a cost function to determine optimization
parameters for predicting a stress attribute of the subsurface volume in real-
time based on the precomputed strain, stress, and displacement values.
The computer-executable method of claim 4, wherein the
predicted optimization parameters are applied during one of:
fitting a local perturbed stress field to a far field stress value in the
subsurface volume, in real-time; or
computing a perturbed stress, a displacement field, and a displacement
on faults around a complex faulted area, in real-time.
6. The computer-executable method of claim 4, wherein simulating
linearly independent far field stress tensors in constant time enables the
principle of superposition to be applied to the simulations regardless of an
individual complexity of each model associated with a given simulation.
7. The computer-executable method of claim 4, wherein the multiple
types of geologic data are weighted differently and wherein each datum is
capable of supporting a weight for each coordinate.
39

8. The computer-executable method of claim 4, wherein the multiple
types of geologic data are obtained from multiple sources, including one of
seismic interpretation data, well bore data, or field observation data; and
wherein the multiple types of geologic data include one of fault
geometry data, fracture orientation data, stylolites orientation data,
secondary
fault plane data, fault throw data, slickenline data, global positioning
system
(GPS) data, interferometric synthetic aperture radar (InSAR) data, laser
ranging data, tilt-meter data, displacement data for a geologic fault, or
stress
magnitude data for the geologic fault.
9. The computer-executable method of claim 4, wherein the
predicted stress attribute further comprises one of a stress inversion, a
stress
field, a far field stress value, a stress interpolation in a complex faulted
reservoir, a perturbed stress field, a stress ratio and associated
orientation,
one or more tectonic events, a displacement discontinuity of a fault, a fault
slip, an estimated displacement, a perturbed strain, a slip distribution on
faults,
quality control on interpreted faults, fracture prediction, prediction of
fracture
propagation according to perturbed stress field, real-time computation of
perturbed stress and displacement fields while performing interactive
parameters estimation, or discernment of an induced fracture from a
preexisting fracture
10. The computer-executable method of claim 4, further comprising
applying a Monte Carlo technique for iteratively minimizing the cost function.
11. The computer-executable method of claim 4, further comprising:
applying fracture orientation data from well bore data to recover a far
field stress and displacement continuity on active faults;

deforming initially flattened horizons with a displacement field from the
far field stress and the displacement continuity;
comparing a geometry of the deformed horizons to an interpretation of
the horizons; and
detecting mismatches between the deformed horizons and the
interpreted horizons as a quality control.
12. The
computer-executable method of claim 4, wherein the multiple
types of geologic data include only fracture orientation data and/or secondary
fault plane data from well bores, and the predicted stress attribute comprises
one or more recovered tectonic events, wherein a recovered stress tensor is
given by an orientation and ratio of principal magnitudes.
13. A computer-readable storage medium tangibly containing
instructions, which when executed by a computer, perform a process,
including:
receiving fault geometry data for one or more faults in a subsurface
volume;
receiving at least one data set containing data points of collected data
including fault orientation data, displacement data, and/or stress magnitude
data in the subsurface volume;
simulating three linearly independent far field stress tensor models in
constant time to generate three sets of component strain, stress, and/or
displacement values for each of the data points;
computing strain, stress, and/or displacement values for each data point
based on a principle of superposition applied to the three stress tensor
models;
iteratively determining optimization parameters to fit the precomputed
strain, stress, and displacement values to a far field stress value for each
data
point in the subsurface volume; and
41

applying the optimization parameters to the precomputed strain, stress,
and/or displacement values to predict stress and displacement parameters in
the subsurface volume in real-time.
14. The computer-readable storage medium of claim 13, further
comprising instructions for receiving a user-selected tectonic stress for the
subsurface volume; and calculating the resulting stress and displacement
parameters in the subsurface volume in real-time.
15. The computer-readable storage medium of claim 13, further
comprising instructions for representing faults in the subsurface volume as a
model comprising triangulated surfaces with discontinuous displacements.
16. The computer-readable storage medium of claim 13, wherein
iteratively determining the optimization parameters further includes:
selecting a set of optimization parameters;
scaling the strain, stress, and displacement values by the optimization
parameters;
assessing a cost of the scaled strain, stress, and displacement values;
wherein when the cost is not substantially minimized during an
assessment iteration, then selecting a next set of optimization parameters,
rescaling the strain, stress, and displacement values according to the next
set
of optimization parameters, and reassessing the cost in a subsequent
assessment iteration; and
when the cost is substantially minimized, then applying the scaled
strain, stress, and displacement values to predict a stress result in the
subsurface volume.
42

17. The computer-readable storage medium of claim 13, further
comprising instructions to predict a tectonic stress in real-time based on the
precomputed strain, stress, and displacement values.
18. The computer-readable storage medium of claim 13, wherein at
least one data set is weighted and wherein at least one data set comprises
one of fault geometry data, fracture orientation data, stylolites orientation
data,
secondary fault plane data, fault throw data, slickenline data, global
positioning system (GPS) data, interferometric synthetic aperture radar
(InSAR) data, laser ranging data, tilt-meter data, displacement data for a
geologic fault, or stress magnitude data for the geologic fault.
19. The computer-readable storage medium of claim 13, wherein the
three linearly independent far field stress tensor models are each based on
independent data sources, each data source obtained from one of seismic
interpretation data, well bore data, or field observation data.
20. The computer-readable storage medium of claim 13, wherein the
predicted stress and displacement parameters comprise one of a stress
inversion, a stress field, a far field stress value, a stress interpolation in
a
complex faulted reservoir, a perturbed stress field, a stress ratio and
associated orientation, one or more tectonic events, a displacement
discontinuity of a fault, a fault slip, an estimated displacement, a perturbed
strain, a slip distribution on faults, a quality control on interpreted
faults,
fracture prediction, prediction of fracture propagation according to perturbed
stress field, real-time computation of perturbed stress and displacement
fields
while performing interactive parameters estimation, or discernment- of an
induced fracture from a preexisting fracture.
43

Description

Note: Descriptions are shown in the official language in which they were submitted.


CA 02735038 2014-03-24
= ,
50866-123
STRESS AND FRACTURE MODELING USING THE
PRINCIPLE OF SUPERPOSITION
[0001]
BACKGROUND
[0002] A fault may be considered a finite complex three-dimensional
surface discontinuity in a volume of earth or rock. Fractures, such as joints,
veins, dikes, pressure solution seams with stylolites, and so forth, may be
propagated intentionally, to increase permeability in formations such as
shale,
in which optimizing the number, placement, and size of fractures in the
formation increases the yield of resources like shale gas.
[0003] Stress, in continuum mechanics, may be considered a measure
of the internal forces acting within a volume. Such stress may be defined as a
measure of the average force per unit area at a surface within the volume on
which internal forces act. The internal forces are typically produced between
the particles in the volume as a reaction to external forces applied on the
volume.
[0004] Understanding the origin and evolution of faults and the
tectonic
history of faulted regions can be accomplished by relating fault orientation,
slip

CA 02735038 2011-03-23
direction, geologic and geodetic data to the state of stress in the earth's
crust.
In conventional inverse problems, the directions of the remote principal
stresses and a ratio of their magnitudes are constrained by analyzing field
data on fault orientations and slip directions as inferred from artifacts such
as
striations on exposed fault surfaces.
[0005] Conventional methods for stress inversion, using measured
striations and/or throw on faults, are mainly based on the assumptions that
the
stress field is uniform within the rock mass embedding the faults (assuming no
perturbed stress field), and that the slip on faults has the same direction
and
sense as the resolved far field stress on the fault plane. However, it has
been
shown that slip directions are affected by: anisotropy in fault compliance
caused by irregular tip-line geometry; anisotropy in fault friction (surface
corrugations); heterogeneity in host rock stiffness; and perturbation of the
local
stress field mainly due to mechanical interactions of adjacent faults.
Mechanical interactions due to complex faults geometry in heterogeneous
media should be taken into account while doing the stress inversion. By doing
so, determining the parameters of such paleostress in the presence of multiple
interacting faults requires running numerous simulations, and therefore an
enormous amount of computation time in order to fit the observed data. The
conventional parameters space has to be scanned for all possibilities, and for
each simulation the model and the post-processes need to be recomputed.
[0006] Motion equations are typically not invoked while using
conventional methods, and perturbations of the local stress field by fault
slip
are ignored. The mechanical role played by the faults in tectonic deformation
is not explicitly included in such analyses. Still, a relatively full
mechanical
treatment is applied for conventional paleostress inversion. However, the
results might be greatly improved if multiple types of data could be used to
better constrain the inversion.
2

CA 02735038 2014-03-24
50866-123
SUMMARY
[0006a] According to an aspect of the present invention, there is
provided a
computer-readable storage medium storing an executable program, the executable
program containing instructions, which when executed by a computer, perform a
process, comprising: simulating three linearly independent far field stress
models for
a subsurface earth volume; computing stress, strain, or displacement values
for data
points in the subsurface earth volume, based on a principle of superposition
applied
to the three linearly independent far field stress models; and predicting a
stress
attribute of the subsurface earth volume in real-time, based on the
precomputed
stress, strain, and displacement values.
[0006b] According to another aspect of the present invention, there is
provided
a computer-executable method, comprising: simulating linearly independent far
field
stress tensors in constant time based on multiple types of geologic data
associated
with a subsurface volume; computing precomputed strain, stress, or
displacement
values based on a principle of superposition applied to the simulations; and
iteratively
minimizing a cost function to determine optimization parameters for predicting
a
stress attribute of the subsurface volume in real-time based on the
precomputed
strain, stress, and displacement values.
[0006c] According to another aspect of the present invention, there is
provided
a computer-readable storage medium tangibly containing instructions, which
when
executed by a computer, perform a process, including: receiving fault geometry
data
for one or more faults in a subsurface volume; receiving at least one data set
containing data points of collected data including fault orientation data,
displacement
data, and/or stress magnitude data in the subsurface volume; simulating three
linearly independent far field stress tensor models in constant time to
generate three
sets of component strain, stress, and/or displacement values for each of the
data
points; computing strain, stress, and/or displacement values for each data
point
based on a principle of superposition applied to the three stress tensor
models;
3

CA 02735038 2014-03-24
50866-123
iteratively determining optimization parameters to fit the precomputed strain,
stress,
and displacement values to a far field stress value for each data point in the
subsurface volume; and applying the optimization parameters to the precomputed
strain, stress, and/or displacement values to predict stress and displacement
parameters in the subsurface volume in real-time.
[0007] In some embodiments, stress and fracture modeling using the
principal
of superposition is provided. An example system simulates linearly independent
far
field stress models for a subsurface earth volume, computing stress, strain,
and/or
displacement values based on superposition of independent stress tensors.
Based
on the precomputed values, the system can generate real-time results, such as
recovery of paleostress values, or, stress, strain, and displacement
parameters for
any point in the subsurface volume as the user varies the far field stress
value. The
system may recover one or more tectonic events, or a stress tensor represented
by a
ratio of principal magnitudes and associated orientation, using only fault
geometry,
well bore data (including fracture orientation and secondary fault plane
data), GPS,
InSAR, folded and faulted horizons, tiltmeters, slip and slikenlines on
faults. The
system can use different types of geologic data from seismic interpretation,
well bore
readings, and field observation to provide numerous results, such as predicted
fracture propagation based on perturbed stress field.
[0008] This summary section is not intended to give a full description of
stress
and fracture modeling using the principle of superposition, or to provide a
comprehensive list of features and elements. A detailed description with
example
implementations follows.
3a

CA 02735038 2011-03-23
BRIEF DESCRIPTION OF THE DRAWINGS
[0009] FIG. 1 is a diagram of an example stress and fracture modeling
system.
[0010] FIG. 2 is a block diagram of an example computing environment
for performing stress and fracture modeling using the principle of
superposition.
[0051] FIG. 3 is a block diagram of an example stress and fracture
modeling engine.
[0052] FIG. 4 is a block diagram comparing techniques for recovering
paleostress.
[0053] FIG. 5 is a diagram of an example method applied to fracture and
conjugate fault planes using data sets without magnitude information.
[0054] FIG. 6 is diagram comparing example stress inversion for normal
fault and thrust fault configurations.
[0011] FIG. 7 is a diagram of cost functions for the stress inversions of
Fig. 6.
[0012] FIG. 8 is a diagram of an example model configuration showing
InSAR data points and a fault surface.
[0013] FIG. 9 is a diagram comparing fringes from an original InSAR grid
and an InSAR grid recovered by an example stress and fracture modeling
method using the principle of superposition.
[0014] Fig. 10 is a diagram showing a top view and a front perspective
view of a plot of the cost surface for an example method applied to the InSAR
case of Figs. 8-9, with cost solution indicated.
[0015] FIG. 11 is a diagram showing results of an example method
applied to a flattened horizon scenario.
[0016] Fig. 12 is a flow diagram of an example method of stress and
fracture modeling using the principle of superposition.
4

CA 02735038 2011-03-23
,
[0017]
Fig. 13 is a flow diagram of an example method of stress and
fracture modeling using the principle of superposition and a cost function.

CA 02735038 2011-03-23
DETAILED DESCRIPTION
Overview
[0018] This
disclosure describes stress and fracture modeling using the
principle of superposition. Given diverse input data, such as faults geometry,
and selectable or optional data sets or data measures, including fault throw,
dip-slip or slickenline directions, stress measurements, fracture data,
secondary fault plane orientations, global positioning system (GPS) data,
interferometric synthetic aperture radar (InSAR) data, geodetic data from
surface tilt-meters, laser ranging, etc., the example system can quickly
generate or recover numerous types of results. The systems and methods
described herein apply the principle of superposition to fault surfaces with
complex geometry in 3D (not only planar), and the faults are, by nature, of
finite
dimension and not infinite or semi-infinite. The results are often rendered in
real-
time and may include, for example, real-time stress, strain, and/or
displacement parameters in response to a user query or an updated
parameter; remote stress states for multiple tectonic events; prediction of
intended future fracturing; differentiation of preexisting fractures from
induced
fractures; and so forth. The diverse input data can be derived from well bore
data, seismic interpretation, field observation, etc.
[0019] The
example systems and methods described below are
applicable to many different reservoir and subsurface operations, including
exploration and production operations for natural gas and other hydrocarbons,
storage of natural gas, hydraulic fracturing and matrix stimulation to
increase
reservoir production, water resource management including development and
environmental protection of aquifers and other water resources, capture and
underground storage of carbon dioxide (CO2), and so forth.
[0020] In an
example implementation, a system applies a 3-dimensional
(3D) boundary element technique using the principle of superposition that
applies to linear elasticity for heterogeneous, isotropic whole- or half-space
6

CA 02735038 2011-03-23
media. Based on precomputed values, the example system can assess a
cost function to generate real-time results, such as stress, strain, and
displacement parameters for any point in a subsurface volume as the user
varies the far field stress value. In one implementation, the system uses only
fault geometry and well bore data, including, e.g., fracture orientation,
secondary fault plane data, and/or in-situ stress measurement by hydraulic
fracture, to recover one or more tectonic events or a stress tensor
represented
by a ratio of principal magnitudes and the associated orientation. The system
can use many different types of geologic data from seismic interpretation,
well
bore readings, and field observation to provide a variety of results, such as
predicted fracture propagation based on perturbed stress field.
[0021] Fig. 1 shows an example stress and fracture modeling system
100. The example system 100 can solve both common and unusual
geomechanical problems. The faults geometry is often known (and optionally,
measured fault throw and imposed inequality constraints such as normal,
thrust, etc., may be knowns). The user may typically have access to data from
well bores (fracture orientation, in-situ stress measurements, secondary fault
planes), geodetic data (InSAR, GPS, and tilt-meter), as well as interpreted
horizons. An example stress and fracture modeling engine 102 and/or
corresponding example methods can recover the remote stress state and
tectonic regime for relevant tectonic event(s), as well as displacement
discontinuity on faults, and estimate the displacement and perturbed strain
and stress fields anywhere within the system, for example.
[0022] Using the principle of superposition, the example stress and
fracture modeling system 100 or engine 102 may perform each of three
linearly independent simulations of stress tensor models in constant time
regardless of the complexity of each underlying model. Each model does not
have to be recomputed. Then, as introduced above, applications for the
example system 100 may include stress interpolation and fracture modeling,
recovery of tectonic event(s), quality control on interpreted faults, real-
time
7

CA 02735038 2011-03-23
computation of perturbed stress and displacement fields when the user is
performing parameters estimation, prediction of fracture propagation,
distinguishing preexisting fractures from induced fractures, and numerous
other applications.
Example Environment
[0023] Fig. 2 shows the example system 100 of Fig. 1 in the context of a
computing environment, in which stress and fracture modeling using the
principle of superposition can be performed.
[0024] In the shown implementation, a computing device 200
implements a component, such as the example stress and fracture modeling
engine 102. The example stress and fracture modeling engine 102 is
illustrated as software, but can be implemented as hardware or as a
combination of hardware and software instructions.
[0025] In the illustrated example, the computing device 200 is
communicatively coupled via sensory devices (and control devices) with a
real-world setting, for example, an actual subsurface earth volume 202,
reservoir 204, depositional basin, seabed, etc., and associated wells 206 for
producing a petroleum resource, for water resource management, or for
carbon services, and so forth.
[0026] The computing device 200 may be a computer, computer
network, or other device that has a processor 208, memory 210, data storage
212, and other associated hardware such as a network interface 214 and a
media drive 216 for reading and writing a removable storage medium 218.
The removable storage medium 218 may be, for example, a compact disk
(CD); digital versatile disk / digital video disk (DVD); flash drive, etc. The
removable storage medium 218 contains instructions, which when executed
by the computing device 200, cause the computing device 200 to perform one
or more example methods described herein. Thus, the removable storage
medium 218 may include instructions for implementing and executing the
8

CA 02735038 2011-03-23
example stress and fracture modeling engine 102. At least some parts of the
example stress and fracture modeling engine 102 can be stored as
instructions on a given instance of the removable storage medium 218,
removable device, or in local data storage 212, to be loaded into memory 210
for execution by the processor 208.
[0027] Although the illustrated example stress and fracture modeling
engine 102 is depicted as a program residing in memory 210, a stress and
fracture modeling engine 102 may be implemented as hardware, such as an
application specific integrated circuit (ASIC) or as a combination of hardware
and software.
[0028] In this example system, the computing device 200 receives
incoming data 220, such as faults geometry and many other kinds of data,
from multiple sources, such as well bore measurements 222, field observation
224, and seismic interpretation 226. The computing device 200 can receive
many types of data sets 220 via the network interface 214, which may also
receive data from the Internet 228, such as GPS data and InSAR data.
[0029] The computing device 200 may compute and compile modeling
results, simulator results, and control results, and a display controller 230
may
output geological model images and simulation images and data to a display
232. The images may be 2D or 3D simulation 234 of stress and fracture
results using the principle of superposition. The example stress and fracture
modeling engine 102 may also generate one or more visual user interfaces
(U1s) for input and display of data.
[0030] The example stress and fracture modeling engine 102 may also
generate or ultimately produce control signals to be used via control devices,
e.g., such as drilling and exploration equipment, or well control injectors
and
valves, in real-world control of the reservoir 204, transport and delivery
network, surface facility, and so forth.
[0031] Thus, an example system 100 may include a computing device
200 and interactive graphics display unit 232. The computing system as a
9

CA 02735038 2011-03-23
whole may constitute simulators, models, and the example stress and fracture
modeling engine 102.
Example Engine
[0032] Fig. 3
shows the example stress and fracture modeling engine
102 in greater detail than in Fig. 1 and Fig. 2. The illustrated
implementation
is only one example configuration for the sake of description, to introduce
features and components of an engine that performs the example stress and
fracture modeling using the principle of superposition. The
illustrated
components are only examples. Different configurations or combinations of
components than those shown may be used to perform the stress and fracture
modeling functions, and different or additional components may also be used.
As introduced above, the example stress and fracture modeling engine 102
can be implemented in hardware, or in combinations of hardware and
software. Illustrated components are communicatively coupled with each
other for communication as needed. Arrows are shown only to suggest
process flow or data flow, since the components can communicate with each
other as needed.
[0033] The
example stress and fracture modeling engine 102 illustrated
in Fig. 3 includes a buffer for data sets 302 or at least access to the data
sets
302, an initialization engine 304, stress model simulators 306 or at least an
access to the stress model simulators 306, a optimization parameters selector
308, a cost assessment engine 310, and a buffer or output for results 312.
These components are shown for descriptive purposes. Other components,
or other arrangement of the components, can enable various implementations
of the example stress and fracture modeling engine 102. The functionality of
the example stress and fracture modeling engine 102 will be described next.

CA 02735038 2011-03-23
Operation of the Example System and Engine
[0034] Fig. 4
shows various methods for recovering paleostress. Fig.
4(a) shows a conventional technique that uses only slip distributions on
faults
as data input, but provides unstable results. The technique of Fig. 4(a)
cannot
augment the slip distributions with alternative forms of data input, such as
GPS data, and so forth. Fig. 4(b) depicts a conventional (2D) Monte Carlo
method, but without the optimization using the principle of superposition.
Fig.
4(c) shows an example method as described herein, including techniques that
utilize the principle of superposition and drastically reduce the complexity
of
the model. The example method shown in Fig. 4(c) may be implemented by
the example stress and fracture modeling engine 102. For example, in one
implementation the initialization engine 304, through the stress model
simulators 306, generates three precomputed models of the far field stress
associated with a subsurface volume 202. For each of the three models, the
initialization engine 304 precomputes, for example, displacement, strain,
and/or stress values. The optimization parameters selector 308 iteratively
scales the displacement, strain, and/or stress values for each superpositioned
model to minimize a cost at the cost assessment engine 310. The values thus
optimized in real-time are used to generate particular results 312.
[0035] In one
implementation, the three linearly independent far field
stress parameters are: (i) orientation toward the North, and (ii) & (iii) the
two
principal magnitudes. These far field stress parameters are modeled and
simulated to generate a set of the variables for each of the three simulated
models. Four variables can be used:
displacement on a fault, the
displacement field at any data point or observation point, a strain tensor at
each observation point, and the tectonic stress. The optimization parameters
selector 308 selects an alpha for each simulation, i.e., a set of "alphas" for
the
three simulated stress models to act as changeable optimization parameters
for iteratively converging on values for these variables in order to minimize
one or more cost functions, to be described below. In one implementation, the
11

CA 02735038 2011-03-23
optimization parameters selector 308 selects optimization parameters at
random to begin converging the scaled strain, stress, and or displacement
parameters to lowest cost. When the scaled (substantially optimized)
parameters are assessed to have the lowest cost, the scaled strain, stress,
and/or displacement parameters can be applied to predict a result, such as a
new tectonic stress.
[0036] Because the example method of Fig. 4(c) uses precomputed
values superpositioned from their respective simulations, the example method
or the stress and fracture modeling engine 102 can provide results quickly,
even in real-time. As introduced above, the stress and fracture modeling
engine 102 can quickly recover multiple tectonic events responsible for
present conditions of the subsurface volume 202, or more quickly discern
induced fracturing from preexisting fracturing than conventional techniques,
or
provide real-time parameter estimation as the user varies a stress parameter,
or can rapidly predict fracturing, and so forth.
[0037] While conventional paleostress inversion may apply a full
mechanical scenario, the stress and fracture modeling engine 102 improves
upon conventional techniques by using multiple types of data. Data sets 302
to be used are generally of two types: those which provide only orientation
information (such as fractures, secondary fault planes with internal friction
angle, and fault striations, etc.), and those which provide magnitude
information (such as fault slip, GPS data, InSAR data, etc.). Conventionally,
paleostress inversion has been computed using slip measurements on fault
planes.
[0038] The conventional technique shown in Fig. 4(a) executes an
operation that proceeds in two stages: (i) resolving the initial remote stress
tensor ap onto the fault elements that have no relative displacement data and
solving for the unknown relative displacements (bj in Fig. 4(a)); and, (ii)
using
the computed and known relative displacements to solve for ciR (Fig. 4(a)). An
12

CA 02735038 2011-03-23
iterative solver is conventionally used, that cycles between stages (i) and
(ii)
until convergence.
[0039] The conventional technique shown in Fig. 4(b) is based on a
Monte-Carlo algorithm. However, this conventional technique proves to be
unfeasible since it requires a long computation time, for which the complexity
is C grza + where n and p are a number of triangular elements modeling the
faults and the number of data points, respectively. For a given simulation, a
random far field stress aR is chosen, and the corresponding displacement
discontinuity b on faults is computed. Then, as a post-process at data points,
and depending on the type of measurements, cost functions are computed
using either the displacement, strain, or stress field. In this conventional
scenario, for hundreds of thousands of simulations, the best cost (close to
zero) is retained as a solution.
[0040] On the other hand, the example method diagrammed in Fig. 4(c),
which can be executed by the stress and fracture modeling engine 102,
extends inversion for numerous kinds of data, and provides a much faster
modeling engine 102. For example, a resulting fast and reliable stress
inversion is described below. The different types of data can be weighted and
combined together. The stress and fracture modeling engine 102 can quickly
recover the tectonic event(s) as well as displacement discontinuity on faults
using diverse data sets and sources, and then obtain an estimate of the
displacement and perturbed strain and stress field anywhere within the
medium, using data available from seismic interpretation, well bores, and
field
observations. Applying the principle of superposition allows a user to execute
parameters estimation in a very fast manner.
[0041] A numerical technique for performing the example methods is
described next. Then, a reduced remote tensor used for simulation is
described, and then the principle of superposition itself is described. An
estimate of the complexity is also described.
13

CA 02735038 2011-03-23
[0042] In one
implementation, a formulation applied by the stress and
fracture modeling engine 102 can be executed using IBEM3D, a successor of
POLY3D (POLY3D is described by F. Maerten, P. G. Resor, D. D. Pollard, and
L. Maerten, Inverting for slip on three-dimensional fault surfaces using
angular
dislocations, Bulletin of the Seismological Society of America, 95:1654-1665,
2005, and by A. L. Thomas, Poly3D: a three-dimensional, polygonal element,
displacement discontinuity boundary element computer program with
applications to fractures, faults, and cavities in the earth's crust, Master's
thesis, Stanford University, 1995.) IBEM3D is a boundary element code based
on the analytical solution of an angular dislocation in a homogeneous or
inhomogeneous elastic whole- or half-space. An iterative solver is employed
for speed considerations and for parallelization on multi-core architectures.
(See, for example, F. Maerten, L. Maerten, and M. Cooke, Solving 3d
boundary element problems using constrained iterative approach,
Computational Geosciences, 2009.) However, inequality constraints cannot
be used as they are nonlinear and the principle of superposition does not
apply. In the selected code, faults are represented by triangulated surfaces
with discontinuous displacement. The advantage is that three-dimensional
fault surfaces more closely approximate curvi-planar surfaces and curved tip-
lines without introducing overlaps or gaps.
[0043] Mixed
boundary conditions may be prescribed, and when traction
boundary conditions are specified, the initialization engine 304 solves for
unknown Burgers's components. After the system is solved, it is possible to
compute anywhere, within the whole- or half-space, displacement, strain or
stress at observation points, as a post-process. Specifically, the stress
field at
any observation point is given by the perturbed stress field due to slipping
faults plus the contribution of the remote stress. Consequently, obtaining
only
the perturbed stress field due to the slip on faults is not enough. Moreover,
the estimation of fault slip from seismic interpretation is only given along
the
dip-direction. Nothing
is known along the strike-direction, and a full
14

CA 02735038 2011-03-23
mechanical scenario is necessary to recover the unknown components of the
slip vector as it will impact the perturbed stress field. Changing the imposed
far field stress (orientation and or relative magnitudes) modifies the slip
distribution and consequently the perturbed stress field. In general, a code
such as 1BEM3D is well suited for computing the full displacement vectors on
faults, and has been intensively optimized using an It -matrix technique. The
unknown for purposes of modeling remains the estimation of the far field
stress that has to be imposed as boundary conditions.
[0044] In an example method, which may be implemented by the stress
and fracture modeling engine 102, a model composed of multiple fault
surfaces is subjected to a constant far field stress tensor up defined in the
global coordinate system by Equation (1):
la11
0-R = 22 a2,l
(1)
Assuming a sub-horizontal far field stress (but the present methodology is not
restricted to that case), Equation (1) simplifies into Equation (2):
al, ai2 0
aR 0
a3õ (2)
Since the addition of a hydrostatic stress does not change aR , the far field
stress tensor aR can be written as in Equation (3):
R = (111 23; al2a3, 00 = a1
a221
0 k2 (3)
[0055] Consequently, a definition of a far field stress with three
unknowns is obtained, namely (flit, (T'22.11-121.

CA 02735038 2011-03-23
[0056] The far field stress tensor, as defined in Equation (3) can be
computed using only two parameters instead of the three ail .571. ti-10 .
Using
spectral decomposition of the reduced R, Equation (4) may be obtained:
o-R = 40-Ro (4)
where, as in Equation (5):
o_R Ea, o
I. (71 (5)
the far field stress tensor up, is the matrix of principal values, and in
Equation
(6):
R e = [cps 6 ¨ sin el
kin cos I (6)
is the rotation matrix around the global z -axis (since a sub-horizontal
stress
tensor is assumed).
[0057] By writing, in Equation (7):
02 = kai (7)
Equation (4) then transforms into Equation (8):
C7R (0, = GiR70. 101 ]Re
(8)
16

CA 02735038 2011-03-23
[0058] Omitting the scaling parameter cfi due to Property 1 discussed
below (when oi = a in Property 1), aR can be expressed as a functional of
two parameters 6 and k , as in Equation (9):
0- R (0 , 1C) ¨ 1776; Eol ?AR e
(9)
[0059] These two parameters are naturally bounded by Equations (10):
77
2
¨10:5_ k < io
1 it 2
(10)
assuming that uniaxial remote stress starts to occur when k 10 . For k = 1, a
hydrostatic stress tensor is found, which has no effect on the model. Also,
using a lithostatic far field stress tensor (which is therefore a function of
depth
Z) does not invalidate the presented technique, and Equation (9) transforms
into Equation (11):
o-R (0,k,z) = zR78-[01 11R a (11)
which is linearly dependent on z. The simplified tensor definition given by
Equation (9) is used in the coming sections to determine (6,k), or
equivalently
fall, (1-2 f2a 1J, according to measurements.
[0060] It is worth noting that even when 2-dimensional parameter-space
is used for the Monte-Carlo simulation using (0- k ), three components are
still
used for the far field stress, specified by the parameters (a, a2. a. ). The
conversions are given by Expression (12):
(6, k)--= W.:Q. cro., 0-12) --. (cri, a,. a 2) (12)
17

CA 02735038 2011-03-23
where, in Equations (13):
goo
I = cosz 6 + k sin2 6
goi. =-
(k¨ 1)cos 6 sin 0
o-i it = cosz 6 + k sin2 6 (13)
and (al. a2. az) are given by Equation (20), further below.
Principle of Superposition
[0061] The example stress and fracture modeling engine 102 uses the
principle of superposition, a well-known principle in the physics of linear
elasticity, to recover the displacement, strain and stress at any observation
point P using the precomputed specific values from linearly independent
simulations. The principle of superposition stipulates that a given value I
can
be determined by a linear combination of specific solutions.
[0062] In the stress and fracture modeling engine 102, recovering a far
field stress implies recovering the three parameters (ai, rt2. a.). Therefore,
the
number of linearly independent solutions used is three. In other words, in
Equation (14):
f = F(a) = F(crio-(1. + cr.20-41) + cracol
= a1F(o-(1)) +a2F(cr12)) + a3Fign
. ail', + aõ,f,.+ a.f,., (14)
where (a, . a:, a2 ) are real numbers, and g 1 to 3 " (for i = ) are
three linearly
independent remote stress tensors. If F is selected to be the strain, stress
or
displacement Green's functions, then the resulting values E , Cr , and u , at
P
can be expressed as a combination of three specific solutions, as shown
below. Thus, the strain, stress and displacement field for a tectonic loading
18

CA 02735038 2011-03-23
are a linear combination of the three specific solutions, and are given by
Equation (15):
1 (1) (2) 7(3)
Ep = aiEp 4- azEi, +
(1) (2)
up = a,up + t.r.,o-i, + a,..o-i; i
up = 170(1 cr.,u:2) aiu(3) P) + P + P (15)
[0063] Similarly, using (ava2, a2) allows recovery of the displacement
discontinuities on the faults, as in Equation (16):
,
.,
1.õ, '
b, = rtikill + cr2b, + a:do 3); (16)
and any far field stress is also given as a combination of the three
parameters,
as in Equation (17):
,,
(1) 1 1
.- i a )
UR = aia- + ct,o-R . +C(;5ji (17)
Complexity Estimate
[0064] Changing crR usually requires recomputing the entire model in
order to determine the corresponding unknown displacement discontinuities.
Then, at any observation point P, the stress is determined as a
superimposition of the far field stress aR and the perturbed stress field due
to
slipping elements.
[0065] For a model made of n triangular discontinuous elements,
computing the stress state at point P first requires solving for the unknown
displacement discontinuities on triangular elements (for which the complexity
is 0(712), and then performing approximately 350n multiplications using the
standard method. By using the principle of superposition, on the other hand,
the unknown displacement discontinuities on triangular elements do not have
19

CA 02735038 2011-03-23
to be recomputed, and only 18 multiplications need be performed by the
example stress and fracture modeling engine 102. The complexity is
independent of the number of triangular elements within the model, and is
constant in time.
[0066] Some direct applications of the example methods will now be
described, such as real-time evaluation of deformation and the perturbed
stress field while a user varies a far field stress parameter. Paleostress
estimation using different data sets 302 is also presented further below, as
is
an example method to recover multiple tectonic phases, and a description of
how the example method can be used for quality control during fault
interpretation.
Real-time Computation
[0067] Before describing the paleostress inversion method, another
method is described first, for real-time computing displacement discontinuity
on faults, and the displacement, strain, and stress fields at observation
points
while the orientation and/or magnitude of the far field stress is varied.
[0068] If the tectonic stress CFR is given and three independent solutions
are known, there exists a unique triple (ai. a2. az) for which Equation (17)
is
honored, and Equations (15) and (16) can be applied.
[0069] In matrix form, Equation (17) is written in the format shown in
Equation (18):
,
I(ii
2 i 3i
cid a CTO a 06 , ' a ,TR
-0o
a( 1) a (2 ! a (3 :i 41 = a&
01 0 i 01 zy_
(1) (2 i (3"1
all ail ail ' (18)
or, in compact form, as in Equation (19):

CA 02735038 2011-03-23
(19)
[0070] Since the three particular solutions 0-01 are linearly independent,
the system can be inverted, which gives Equation (20):
a = '4;1 o- R (20)
[0071] In Equation (20), Aal is precomputed by the initialization engine
304. Given a user-selected remote stress, (TR , the stress and fracture
modeling engine 102 recovers the three parameters (a, a2. ), then the fault
slip and the displacement, strain and stress field are computed in real-time
using Equations (16) and (15), respectively. To do so, the three particular
solutions of the displacement, strain and stress are stored at initialization
at
each observation point, as well as the displacement discontinuity on the
faults.
In one implementation, the example stress and fracture modeling engine 102
enables the user to vary the orientation and magnitude of aR , and to
interactively display the associated deformation and perturbed stress field.
Paleostress inversion using data sets
[0072] As seen above, the main unknowns while doing forward modeling
for the estimation of the slip distribution on faults, and consequently the
associated perturbed stress field, are the orientation and relative magnitudes
of the far field stress 0R.
[0073] If field measurements are known at some given observation
points (e.g., displacement, strain and/or stress, fractures orientation,
secondary fault planes that formed in the vicinity of major faults, etc.),
then it is
possible to recover the triple (eri- 7.. fr.), and therefore also recover the
tectonic stress un and the corresponding tectonic regime (see Appendix A,
below). The next section describes the method of resolution and the cost
functions used to minimize cost for different types of data sets 302.
21

CA 02735038 2011-03-23
, .
Method of resolution
[0074] Applying a Monte Carlo technique allows the parameters
(ai. a:, a, ) to be found, which minimize the cost functions when given three
independent far field stresses (see Equation 15). However, even if (cr1,c(2,
ag )
imply a 3-dimensional parameters-space, this space can be reduced to two
dimensions (namely, to the parameters 6 and k), the conversion being given
by Equation (20) and (6,0--.(0-o0.0o1,0-11)--. (al, a:. CIO , where, in
Equations
(21):
1 uoc
= cos2 64 + ksjn2 0
aoi ¨ (k¨ i)cos 0 sin 0
ail = cosZ 6 + ksin2 6 (21)
(see also Fig. 4(c) and Algorithm 1 for a detailed description). Consequently,
the searching method (the search for optimized parameters) is accelerated by
reducing the parameters-space by one dimension.
[0075] A simple sampling method can be performed by considering
a
two-dimensional rectangular domain for which the axes correspond to 6 and
k. The 2D-domain is sampled randomly with 111%, points, and the associated
cost function (to be defined in the coming sections) is used to determine the
point of minimum cost. A refinement is then created around the selected point
and the procedure is repeated with a smaller domain. Algorithm (1) depicts a
simplified version of the example procedure, for which there is no refinement.
The example sampling method presented here can be greatly optimized by
various techniques.
ALGORITHM (1): Paleostress Estimation
Input: Faults geometry
22

CA 02735038 2011-03-23
Input: Data set
Output: o-R // the estimated paleostress
Initialize: Compute 3 simulations using 3 linearly independent
a1 (1 3) and the faults geometry. Store the resulting displacement and
stress fields at data points, and the displacement discontinuity at faults, if
necessary.
Let C = 1 II initial cost
Let a (0,0,0,) // initial (a) solution
for = 1 to np do
0 E [- 71 -1
Randomly generate 2, 2
Randomly generate k e [-10,10]
// Convert (6,k) c 9tt(2 ) to (a) E 9iT(3 )
Compute R6 using Equation (6)
Compute a R(0 .k) using Equation (9)
Compute a using Equation (20)
// Compute the corresponding cost
d = cost(a , data set)
if d c then
= cr
e = d
end
end
oR =
f=t
Data Sets
23

CA 02735038 2011-03-23
[0076] The particularity of this method lies in a fact that many different
kinds of data sets 302 can be used to constrain the inversion. Two groups of
data are presented in the following sections: the first one includes only
orientation information and the second includes displacement and/or stress
magnitude information.
Without magnitude information from the data set
[0077] For opening fractures (e.g., joints, veins, dikes) the orientation
of
the normal to the fracture plane indicates the direction of the least
compressive stress direction in 3D (72). Similarly, the normals to pressure
solution seams and stylolites indicate the direction of the most compressive
stress (0-0. Using measurements of the orientations of fractures, pressure
solution seams, and stylolites allows the stress and fracture modeling engine
102 to recover the tectonic regime which generated such features.
[0078] At any observation point P, the local perturbed stress field can be
determined easily from a numerical point of view by using three linearly
independent simulations. Fig. 5 shows fracture and conjugate fault planes.
Fig. 5(a) shows orientation of 0-2 relative to an opening fracture (joints,
veins,
dikes) given by its normal i in 3D. Fig. 5(b) shows the same as Fig. 5(a)
except for an orientation of u. relative to a joint given by its projected
normal
(e.g., trace on outcrop). Fig. 5(c) and Fig. 5(d) show the same as Fig. 5(a)
and Fig. 5(b) except shown for a stylolite. Fig. 5(e) shows orientation of 0-2
and oi relative to conjugate fault-planes given by one of the normal Ti in 3D
and the internal friction angle 6 . The goal is to determine the best fit of
the far
field stress 0-p, and therefore parameters (ai.a..a.), given some orientations
of opening fracture planes for which the normals coincide with the directions
of
the least compressive stress of at P, or equivalently for which the plane of
the
fracture contains the most compressive stress ( Crl as in Fig. 5(a) and Fig.
5(b).
24

CA 02735038 2011-03-23
[0079]By varying the
state of stress at any observation point
(ai. a,. a,),
P can be computed quickly using the three precomputed models. The cost
function to minimize is given in Equation (22):
ffrac(ai= az, az) ¨ 2mX(R- ne) +
(22)
where "." is the dot-product, is the normal to a fracture plane, and in is the
number of observation points.
The minimization of a function of the three parameters is expressed by
Equation (23):
Fir ac = a inAnaµ ffrac cto}
(23)
[0080]
Similarly, for pressure solution seams and stylolites, the cost
function is defined as in Equation (22) using the least compressive stress c72
as in Equation (24) (see Fig. 5(c) and Fig. 5(d)):
2 ^.
P p
fstvi(ai= elvaz) = E(K-TiP) ))
p (24)
Using secondary fault planes
[0081] The
orientation of secondary fault planes that develop in the
vicinity of larger active faults may be estimated using a Coulomb failure
criteria, defined by Equation (25):
1
tan(2.8) = ¨
(25)

CA 02735038 2011-03-23
where 6 is the angle of the failure planes to the maximum principal
compressive stress al and P is the coefficient of internal friction. Two
conjugate failure planes intersect along cr; and the fault orientation is
influenced only by the orientation of the principal stresses and the value of
the
friction.
[0082] The cost function is therefore defined by Equation (26):
1
isfn -it'?. _ 0)2r
ffuaL(a.. (I:- az) = 27-n - P)2 (
(26)
where ui is the direction of the most compressive stress and O is the
direction of the intermediate principal stress. The first term of the right
hand
side in Equation (26) maintains an orthogonality between the computed (T2
and the normal of the fault plane, whereas the second term ensures that the
angle between the computed al and the fault plane is close to 6 (see Fig.
5(e)).
Example 1: Normal and thrust fault
[0083] Fig. 6 shows a synthetic example using an inclined planar fault as
one of two conjugate fault planes selected randomly. The inclination of the
two conjugate fault planes is presented for a normal fault configuration (Fig.
6(a)) and a thrust fault configuration (Fig. 6(b)). Dip-azimuth and dip-angle
of
each conjugate fault plane are used to perform the inversion and the internal
friction angle is 6 = 30. The main active fault is represented by the inclined
rectangular plane 602.
[0084] Initially, the model is constrained by a far field stress at some
observation points 604, where the two conjugate planes are computed using
an internal friction angle of 30 degrees. Then, for each observation point
604,
26

CA 02735038 2011-03-23
one of the conjugate fault planes is chosen randomly and used as input data
for the stress inversion.
[0085] Fig. 7 shows the cost function for the synthetic example from Fig.
6. Fig. 7(a) shows the cost function for the normal fault, and Fig. 7(b) shows
the cost function for the thrust fault. In both cases, the recovered regional
stress tensor, displacement on fault, and predicted conjugate fault planes
provide a good match with the initial synthetic model.
Using fault striations
[0086] In the case of fault striations, the cost function is defined as in
Equation (27):
1
N2
(1
m
(27)
where di; and dr represent the normalized slip vector from a simulation and
the measured slip vector, respectively.
Data sets containing magnitude information
[0087] The magnitude of displacements may be used to determine not
only the stress orientation, but also the magnitude of the remote stress
tensor,
instead of just the principal stress ratio.
[0088] To do so, the procedure is similar to that described previously.
However, given Equations (15) and (16), it is evident that there exists a
parameter 8 for which the computed displacement discontinuity on faults and
the displacement, strain and stress fields at observation points scale
linearly
with the imposed far field stress. In other words, as in Equation (28):
27

CA 02735038 2011-03-23
1
6o- R 6b,,
tio-R57.1p
8CfRoEp
o'ciR6o-p (28)
[0089] This leads to the following property:
[0090] Property 1: Scaling the far field stress by (3 F It scales the
displacement discontinuity on faults as well as the displacement, strain, and
stress fields at observation points by 6 .
[0091] Using this property, all measurements at data points are globally
normalized before any computation and the scaling factor is noted (the
simulations are also normalized, but the scaling factor is irrelevant). After
the
system is solved, the recovered far field stress, displacement and stress
fields
are scaled back by a factor of 67-7-t .
Using GPS data
[0092] In the case of a GPS data set, the cost function is defined in
Equation (29):
1
- -.
lin: ii c \T 2 2
fgpslits- a2.(12) = ¨X I - .-- -... 1." III')I -111-
U."
P /
211 P 11111c2111uT11/a 111411, (29)
where TIT is the globally normalized measured elevation changed at point P
---c
from the horizon and up is the globally normalized computed elevation
change for a given set of parameters (ai=cr2.cra). The first term on the right
hand side in Equation (29) represents a minimization of the angle between the
28

CA 02735038 2011-03-23
two displacement vectors, whereas the second term represents a minimization
of the difference of the norm.
Using InSAR data
[0093] When using an InSAR data set, there are two possibilities. Either
the global displacement vectors of the measures are computed using the
displacement u along the direction of the satellite line of sight in which
case Equation (30) is used:
= riznsar = (30)
and the same procedure that is used for the GPS data set (above) is applied
with the computed 4- or, the computed displacement vectors are computed
along the satellite line of sight, in which case Equation (31) is used:
licp' =11.1 (31)
where "." is the dot product. The cost function is consequently given by
Equation (32):
1
2
= 71-1 1 ¨
11 p-
(32)
Example
[0094] Figs. 8, 9, and 10 present a synthetic example using an InSAR
data set. Fig. 8 shows a model configuration showing the InSAR data points
802 as well as the fault surface 804. Fig. 9 shows a comparison of the fringes
from the original InSAR grid 902 and the recovered InSAR grid 904. Fig. 10
shows a plot of the cost surface, as a function of 6 (x-axis) and k (y-axis).
29

CA 02735038 2011-03-23
On the left is a top view 1002 of the plot, and on the right is a front
perspective
view 1004 of the plot. The minimized cost solution 1006 in each view is
marked by a small white circle (1006).
[0095] To utilize an InSAR data set, a forward model is run using one
fault plane 804 and one observation grid (Fig. 8) at the surface of the half-
space (see surface holding the data points 802). A satellite direction is
selected, and for each observation point 802, the displacement along the
satellite line of sight is computed. Then, an example stress and fracture
modeling method described herein is applied using the second form of the
InSAR cost function given in Equation (32). Fig. 9 compares the original
interferogram 902 (left) to the recovered interferogram 904 (right). Fig. 10
shows how complex the cost surface can be, even for a simple synthetic
model. In one implementation, the cost surface was sampled with 500,000
data points 802 (number of simulations), and took 18 seconds on an average
laptop computer with a 2GHz processor and with 2GB of RAM running on
Linux Ubuntu version 8.10, 32 bits.
Using a flattened horizon
[0096] Using the mean plane of a given seismic horizon (flattened
horizon), the stress and fracture modeling engine 102 first computes the
change in elevation for each point making the horizon. Then, the GPS cost
function is used, for which only the uz component is provided, giving Equation
(33):
1
, 1
fixonzon(ava2. az) = -
m p Urpn
(33)
If pre- or post-folding of the area is observed, the mean plane can no longer
be used as a proxy. Therefore, a smooth and continuous fitting surface has to

CA 02735038 2011-03-23
. .
be constructed which must remove the faulting deformations while keeping the
folds. Then, the same procedure as for the mean plane is used to estimate
the paleostress. In some circumstances and prior to defining the continuous
fitting surface, it is sometimes necessary to filter the input horizon from
noises
possessing high frequencies, such as corrugations and bumps, while
preserving natural deformations.
Example
[0097] Fig. 11 shows results when applying an example stress and
fracture modeling method to a synthetic example using a flattened horizon.
Fig. 11(a) presents a model configuration showing the horizon 1102 and the
fault surface 1104. Fig. 11(b) shows a comparison of the original dip-slip
1106
(left) and the recovered dip-slip 1108 (right). Fig. 11(c) shows a comparison
of
the original strike-slip 1110 (left) and the recovered strike-slip 1112
(right). Fig.
11(d) shows original vertical displacement 1114 (left) from the flattened
horizon
(left) and recovered vertical displacement 1116 (right) from the flattened
horizon.
[0098] In the example shown in Fig. 11, a complex shaped fault
is initially
constrained by a far field stress, and consequently slips to accommodate the
remote stress. At each point of an observation plane cross-cutting the fault,
the stress and fracture modeling engine 102 computes the resulting
displacement vector and deforms the grid accordingly. Then, inversion takes
place using the fault geometry. After flattening the deformed grid, the change
in elevation for each point is used to constrain the inversion and to recover
the
previously imposed far field stress as well as the fault slip and the
displacement field. The comparison of the original and inverted dip-slip (Fig.
11(b)) and strike-slip (Fig. 11(c)) show that they match well (same scale). A
good match is also observed for the displacement field at the observation grid
(Fig. 11(d)).
31

CA 02735038 2011-03-23
Using dip-slip information
[0099] When dip-slip data is used, the cost function is defined as in
Equation (34):
=
b:
fds(ct al.. cza) = ¨ 1 I / 1
m
P (34)
where V is the measured dip-slip magnitude for a triangular element e, and
b5 is the computed dip-slip magnitude.
Using all available information
[00100] The example stress and fracture modeling engine 102 can
combine the previously described cost functions to better constrain stress
inversion using all available data (e.g., fault and fracture plane orientation
data, GPS data, InSAR data, flattened horizons data, dip-slip measurements
from seismic reflection, fault striations, etc.). Furthermore, data can be
weighted differently, and each datum can also support a weight for each
coordinate.
Multiple tectonic events
[00101] For multiple tectonic events, it is possible to recover the major
ones, e.g., those for which the tectonic regime and/or the orientation and/or
magnitude are noticeably different. Algorithm (2), below, presents a way to
determine different events from fracture orientation (joints, stylolites,
conjugate
fault planes) measured along well bores.
[00102] After doing a first simulation, a cost is attached at each
observation point which shows the confidence of the recovered tectonic stress
relative to the data attached to that observation point. A cost of zero means
a
good confidence, while a cost of one means a bad confidence. See Fig. 7 for
an example plot of the cost. By selecting only data points that are under a
32

CA 02735038 2011-03-23
given threshold value and running another simulation with these points, it is
possible to extract a more precise paleostress value. Then, the remaining
data points above the threshold value are used to run another simulation with
the paleostress state to recover another tectonic event. If the graph of the
new cost shows disparities, the example method above is repeated until
satisfactory results are achieved. During the determination of the tectonic
phases, the observation points are classified in their respective tectonic
event.
However, the chronology of the tectonic phases remains undetermined.
ALGORITHM (2): Detecting Multiple Tectonic Events
Input: E as a user threshold in [0, 1]
Input: S = all fractures from all wells
Let: T = 0
while S * 0 do
Simulations: Compute the cost for each fracture in S
if max(cost) <E then
Detect a tectonic event CFR for fractures in S
if T = 0 then
Terminate
else
S=T
T=
continue
end
= fractures below E
T + = fractures above E
end
33

CA 02735038 2011-03-23
Seismic interpretation quality control
[00103] It is useful to have a method for quality control (QC) for
interpreted faults geometries from seismic interpretation. The basic idea is
to
use the fracture orientations from well bores to recover the far field stress
and
the displacement discontinuities on active faults. Then, the computed
displacement field is used to deform the initially flattened horizons. The
geometry of the resulting deformed horizons can be compared with the
interpreted ones. If some mismatches are clearly identified (e.g., interpreted
uplift and computed subsidence), then the fault interpretation is possibly
false.
For example, an interpreted fault might dip in the wrong direction. Note that
an unfolded horizon can be approximated by its mean plane, as described
above in relation to flattened horizons.
Conclusion and perspectives
[001041 The example stress and fracture modeling engine 102 applies the
property of superposition that is inherent to linear elasticity to execute
real-
time computation of the perturbed stress and displacement field around a
complexly faulted area, as well as the displacement discontinuity on faults.
Furthermore, the formulation executed by the example stress and fracture
modeling engine 102 enables rapid paleostress inversion using multiple types
of data such as fracture orientation, secondary fault planes, GPS, InSAR,
fault
throw, and fault slickenlines. In one implementation, using only fracture
orientation and/or secondary fault planes from well bores, the stress and
fracture modeling engine 102 recovers one or more tectonic events, the
recovered stress tensor being given by the orientation and ratio of the
principal
magnitudes. The example stress and fracture modeling engine 102 and
34

CA 02735038 2011-03-23
associated methods can be applied across a broad range of applications,
including stress interpolation in a complexly faulted reservoir, fractures
prediction, quality control on interpreted faults, real-time computation of
perturbed stress and displacement fields while doing interactive parameter
estimation, fracture prediction, discernment of induced fracturing from
preexisting fractures, and so forth.
[00105] In a variation, another important application of the stress and
fracture modeling engine 102 and associated methods is evaluation of the
perturbed stress field (and therefore the tectonic event(s)) for recovering
"shale gas." Since shale has low matrix permeability, gas production in
commercial quantities requires fractures to provide permeability. This is
typically done by hydraulic fracturing to create extensive artificial
fractures
around well bores, and therefore requires a good understanding of how
fractures will propagate according to the perturbed stress field.
Example Methods
[00106] FIG. 12 shows an example method 1200 of stress and fracture
modeling using the principle of superposition. In the flow diagram, the
operations are summarized in individual blocks. The example method 1200
may be performed by hardware or combinations of hardware and software, for
example, by the example stress and fracture modeling engine.
[00107] At block 1202, linearly independent stress models for a
subsurface volume are simulated.
[00108] At block 1204, stress, strain, and/or displacement parameters for
the subsurface volume are computed, based on a superposition of the linearly
independent stress models.

CA 02735038 2011-03-23
[00109] At block 1206, an attribute of the subsurface volume is iteratively
predicted based on the precomputed stress, strain, and/or displacement
values.
[00110] FIG. 13 shows an example method 1300 of stress and fracture
modeling using the principle of superposition. In the flow diagram, the
operations are summarized in individual blocks. The example method 1300
may be performed by hardware or combinations of hardware and software, for
example, by the example stress and fracture modeling engine.
[00111] At block 1302, faults geometry for a subsurface earth volume is
received.
[00112] At block 1304, at least one data set associated with the
subsurface earth volume is also received.
[00113] At block 1306, three linearly independent far field stress tensor
models are simulated in constant time to generate strain, stress, and/or
displacement values.
[00114] At block 1308, a superposition of the three linearly independent
far field stress tensor models is computed to provide precomputed strain,
stress, and/or displacement values.
[00115] At block 1310, a post-process segment of the method
commences, that can compute numerous real-time results based on the
principle of superposition.
[00116] At block 1312, optimization parameters for each of the linearly
independent far field stress tensor models are selected.
[00117] At block 1314, the precomputed stress, strain, and/or
displacement values are scaled by the optimization parameters.
[00118] At block 1316, a cost associated with the scaled precomputed
stress, strain, and/or displacement values is evaluated. If the cost is not
satisfactory, then the method loops back to block 1312 to select new
optimization parameters. If the cost is satisfactory, then the method
continues
to block 1318.
36

CA 02735038 2014-03-24
50866-123
[00119] At block 1318, the scaled strain, stress, and/or displacement .
values are applied to the subsurface earth volume, e.g., with respect to a
query about the subsurface earth volume or in response to an updated
parameter about the subsurface earth volume.
[00120] At block 1320, a query or updated parameter regarding the
subsurface earth volume is received, that seeds or initiates generation of the
post-process results in the real-time results section (1310) of the method
1300.
Conclusion
[00121] Exemplary systems and methods have been described in language
specific to structural features and/or methodological acts. It is to be
understood that
the specific features and acts are disclosed as exemplary forms of
implementing the
claimed systems, methods, and structures.
=
=
37

CA 02735038 2011-03-23
APPENDIX A
For a given far field stress tensor, the tectonic regime can be easily
determined using the principal values Gri and vectors Vi , (ordered in
decreasing order). Since one axis is vertical, the tectonic regime is given by
Equation (35):
{ Trust fault regime: if di is vertical
Tectonic regime = StrikeEislip faukregime: if ii2 is vertical
Normal fault regime: if (13 is vertical (35)
37a

Representative Drawing
A single figure which represents the drawing illustrating the invention.
Administrative Status

2024-08-01:As part of the Next Generation Patents (NGP) transition, the Canadian Patents Database (CPD) now contains a more detailed Event History, which replicates the Event Log of our new back-office solution.

Please note that "Inactive:" events refers to events no longer in use in our new back-office solution.

For a clearer understanding of the status of the application/patent presented on this page, the site Disclaimer , as well as the definitions for Patent , Event History , Maintenance Fee  and Payment History  should be consulted.

Event History

Description Date
Common Representative Appointed 2019-10-30
Common Representative Appointed 2019-10-30
Change of Address or Method of Correspondence Request Received 2018-03-28
Grant by Issuance 2014-08-19
Inactive: Cover page published 2014-08-18
Inactive: Final fee received 2014-04-29
Pre-grant 2014-04-29
Notice of Allowance is Issued 2014-04-10
Letter Sent 2014-04-10
Notice of Allowance is Issued 2014-04-10
Inactive: Q2 passed 2014-04-07
Inactive: Approved for allowance (AFA) 2014-04-07
Letter Sent 2014-03-27
Reinstatement Request Received 2014-03-24
Amendment Received - Voluntary Amendment 2014-03-24
Reinstatement Requirements Deemed Compliant for All Abandonment Reasons 2014-03-24
Reinstatement Requirements Deemed Compliant for All Abandonment Reasons 2014-03-24
Inactive: Abandoned - No reply to s.30(2) Rules requisition 2013-11-01
Inactive: Abandoned - No reply to s.29 Rules requisition 2013-11-01
Inactive: S.30(2) Rules - Examiner requisition 2013-05-01
Inactive: S.29 Rules - Examiner requisition 2013-05-01
Application Published (Open to Public Inspection) 2011-09-25
Inactive: Cover page published 2011-09-25
Letter Sent 2011-06-22
Inactive: Single transfer 2011-05-30
Inactive: IPC assigned 2011-05-18
Inactive: First IPC assigned 2011-05-18
Inactive: IPC assigned 2011-05-18
Inactive: Filing certificate - RFE (English) 2011-04-07
Application Received - Regular National 2011-04-07
Letter Sent 2011-04-07
Request for Examination Requirements Determined Compliant 2011-03-23
All Requirements for Examination Determined Compliant 2011-03-23

Abandonment History

Abandonment Date Reason Reinstatement Date
2014-03-24

Maintenance Fee

The last payment was received on 2014-02-11

Note : If the full payment has not been received on or before the date indicated, a further fee may be required which may be one of the following

  • the reinstatement fee;
  • the late payment fee; or
  • additional fee to reverse deemed expiry.

Patent fees are adjusted on the 1st of January every year. The amounts above are the current amounts if received by December 31 of the current year.
Please refer to the CIPO Patent Fees web page to see all current fee amounts.

Owners on Record

Note: Records showing the ownership history in alphabetical order.

Current Owners on Record
SCHLUMBERGER CANADA LIMITED
Past Owners on Record
FRANTZ MAERTEN
LAURENT MAERTEN
Past Owners that do not appear in the "Owners on Record" listing will appear in other documentation within the application.
Documents

To view selected files, please enter reCAPTCHA code :



To view images, click a link in the Document Description column. To download the documents, select one or more checkboxes in the first column and then click the "Download Selected in PDF format (Zip Archive)" or the "Download Selected as Single PDF" button.

List of published and non-published patent-specific documents on the CPD .

If you have any difficulty accessing content, you can call the Client Service Centre at 1-866-997-1936 or send them an e-mail at CIPO Client Service Centre.


Document
Description 
Date
(yyyy-mm-dd) 
Number of pages   Size of Image (KB) 
Abstract 2011-03-22 1 30
Description 2011-03-22 38 1,366
Drawings 2011-03-22 13 394
Claims 2011-03-22 6 246
Representative drawing 2011-09-06 1 36
Description 2014-03-23 39 1,422
Claims 2014-03-23 6 256
Acknowledgement of Request for Examination 2011-04-06 1 189
Filing Certificate (English) 2011-04-06 1 166
Courtesy - Certificate of registration (related document(s)) 2011-06-21 1 104
Reminder of maintenance fee due 2012-11-25 1 111
Courtesy - Abandonment Letter (R30(2)) 2013-12-29 1 164
Courtesy - Abandonment Letter (R29) 2013-12-29 1 164
Notice of Reinstatement 2014-03-26 1 170
Commissioner's Notice - Application Found Allowable 2014-04-09 1 161
Correspondence 2014-04-28 2 75