ProAZ 2014 11 DH
Uploaded by
Juan Francisco VillacresesProAZ 2014 11 DH
Uploaded by
Juan Francisco VillacresesProAZ Course
Azimuthal Attributes & Analysis
November, 2014
Table of Contents
Introduction to fracture analysis and anisotropy
Exercise 1 : Haynesville example – Setting up the project
Observing seismic anisotropy
Exercise 2 : Azimuthal gathers
Anisotropy and fractured rock physics
Exercise 3 : Azimuthal modeling
AVAz analysis & the Rüger equation
Exercise 4 : Azimuthal AVO (AVAz) analysis
Azimuthal Fourier Coefficients
Exercise 5 : Calculating Azimuthal Fourier Coefficients
Azimuthal Inversion
Acquisition and Processing considerations
Exercise 6 : Rose Diagram & Visualization
Velocity versus azimuth (VVAz) analysis
Azimuthal Case Study
References & Appendices
2 November, 2014 (DH)
Theme
Remotely detecting information about fractures and the
stress field is an important objective in the development
of unconventional and tight hydrocarbon reservoirs.
The direct detection of fractures is beyond the resolution
of traditional poststack seismic.
However, fractures and stress may cause the earth to
become anisotropic which is seismically observable.
By observing the prestack azimuthal amplitude and
traveltime variations it is possible to infer the presence
of fractures and their orientation.
ProAZ is new software from Hampson Russell to help
the interpreter observe and quantify these azimuthal
3
variations. November, 2014 (DH)
Objective
In this course we will learn how to predict fractures from P-wave seismic
azimuthal amplitude data using the Hampson-Russell ProAZ software.
In the figure below, azimuthal P-wave seismic data was used to predict the
fracture intensity and direction.
4 November, 2014 (DH)
The Goal – predict fractures and stress field remotely using seismic
These azimuthal
attributes convey both
magnitude and
orientation information.
The magnitude and
direction can be
simultaneously
displayed using vector
displays.
In this example, both
the color and the size
of the platelets convey
fracture intensity
information.
5 November, 2014 (DH)
Azimuthal Introduction to fracture analysis
In this introduction we show
1) How fractures and stress introduce
velocity anisotropy into the earth
2) How this fracture induced anisotropy
introduces
azimuthal traveltime variations in seismic data
Azimuthal amplitude variations in seismic data
Our focus in this course is on P-wave
seismic data, but we need also to
consider S-waves.
6 November, 2014 (DH)
P and S-Wave Velocities
First, let us note that there are two types of waves that travel in the earth with
different velocities:
P-wave, or compressional wave velocity, in which the direction of particle
motion is in the same direction as the wave movement.
S-wave, or shear wave velocity, in which the direction of particle motion is
at right angles to the wave movement.
x
y P-waves S-waves
7 November, 2014 (DH)
Polarized S-Wave Velocities
Because the direction of particle motion in an S-wave is at right angles
to the wave movement, note that there are two possible directions of
movement, or polarizations, of the S-wave.
In the case shown below, one is in the x-direction and one is in the y-
direction (shown by the red dot, coming out of the plane).
These are often called SV
(for vertical) and SH (for
horizontal) polarizations. x
However, this is only for
propagation along natural y S-waves
coordinates. z
For anisotropic rocks, we Polarization
really should talk about S1 in y direction
and S2, or Sfast and Sslow,
polarizations, which depend
on the type of anisotropy.
8 November, 2014 (DH)
Fractures and anisotropy
Slow wave
Fast wave
To see how fractures and anisotropy are related, consider this map view of a set of
fractures. Seismic waves at right angles to the fractures propagate slowly.
Seismic waves parallel to fractures propagate more quickly.
November, 2014 (DH)
9
CMP Bin, Azimuth and Offset
x x x x x x x x x x x x x x x
Slow wave
x x x x x x x x x x x x x x x
x x x x x x x x x x x x x x x
Bin
x x x x x x x x x x x x x x x
Fast wave
x x x x x x x x x x x x x x x
To measure the amplitude and traveltime effects as a function of azimuth, we
record a series of shot-receiver pairs with different offsets and azimuths. These
are gathered into the same bin. November, 2014 (DH)
10
Stress, fractures and anisotropy
Stress can also introduce
anisotropy
Stress induced anisotropy and
fracture induced anisotropy are
closely related.
Stress models typically assume
an isotropic rock with randomly
oriented and distributed cracks
under ambient stress.
Under differential principle
stresses these cracks can
change their shapes making the
rock anisotropic.
11 November, 2014 (DH)
Common Mid Point (CMP) Gather
The CMP gather (and stack) is made up of a number
of different source – receiver combinations all sharing
the same common mid point.
S4 S3 S2 S1R1 R2 R3 R4
vp1
Depth Reflector
vp2
12 November, 2014 (DH)
Common Mid Point (CMP) Gather
The velocity can be measured by observing how the
traveltime varies as a function of Source Receiver offset
Source – Receiver
S4 S3 S2 S1R1 R2 R3 R4 Offset
distance
time
VP1 vp1
time
Depth Reflector
vp2 By measuring the traveltime
can calculate the velocity
13 November, 2014 (DH)
Fractures change both the P-wave and S-wave velocities
P-wave Velocity S-wave Velocity
Layer 1
Top of
fractured zone
Slow velocity Fast velocity
at 90 degrees parallel to Anisotropic Layer 2
to fractures fractures
Base of
fractured zone
Layer 3
14 November, 2014 (DH)
Traveltimes perpendicular to the fracture are slower
than parallel to the fractures
N-S
0 degrees
E-W
90 degrees
15 November, 2014 (DH)
Traveltimes perpendicular to the fracture are slower than
parallel to the fractures
N-S
E-W
Top of
fractured zone
(same moveout)
Fast Slow
Base of
fractured zone
(different moveout)
16 November, 2014 (DH)
After NMO correcting the data with an average velocity
Fast Slow
Top of
fractured zone
Over-corrected Under-corrected
Base of
fractured zone
17 November, 2014 (DH)
Data recorded for 8 azimuth sectors
(with more azimuths easier to see azimuthal variations)
N-S
0 degrees
E-W
90 degrees
18 November, 2014 (DH)
Data recorded for 8 azimuth sectors
(with more azimuths easier to see azimuthal variations)
19 November, 2014 (DH)
Sort by offset, then azimuth
Common Offset / Common Azimuth (COCA) Gathers
Primarily
sorted by
offset
Top of
fractured zone
Offset
Azimuths
Fast Slow
Base of
fractured zone
20 November, 2014 (DH)
Azimuthal Velocities
The sinusoidal variation observed in the azimuthal residual
traveltimes can be modeled using elliptical velocities, which
are characterized by 3 parameters:
Isotropic Major axis: Vfast
media
Minor axis: Vslow
Φ, the azimuth of Vfast
relative to the North
• The orientation of Vfast is parallel to maximum stress.
• Vfast and Vslow are perpendicular.
21 November, 2014 (DH)
What causes the AVO Effect?
Surface
q2 q1
q3
r1 VP1 VS1 Reflector
r2 VP2 VS2
The traces in a seismic gather reflect from the subsurface at increasing angles of
incidence q. For isotropic media, The first order approximation angle dependent
reflection coefficient is given by:
R(q ) A Biso sin 2 q
where the isotropic gradient Biso term produces the AVO effect. It is dependent on
changes in density, r, P-wave velocity, VP, and S-wave velocity, VS.
22 November, 2014 (DH)
What causes the AVAz Effect?
Surface
q2 q1
q3
r1 VP1 VS1 Reflector
r2 VP2 VS2
Since the velocity changes as a function of azimuth the AVO relation also
changes as function of azimuth. This is often approximated by allowing the
gradient to change as a function of azimuth:
R(q , f ) A ( Biso Bani sin2 (f fiso )) sin2 q
Where Bani is the anisotropic gradient term and fiso is the azimuth of the
fracture strike. When the source-receiver azimuth f is identical to fiso it
reduces to the isotropic case.
23 November, 2014 (DH)
Azimuthal AVO 24
Recall the simple fractured model we started with.
Azimuthal AVO
Primarily
sorted by
offset
Top of
fractured zone
Offset
Azimuths
Fast Slow
Base of
fractured zone
25 November, 2014 (DH)
Need to flatten the data prior to AVAz analysis
Application of Azimuthal NMO or Trim Statics
Top of
fractured zone
Azimuths
Base of
fractured zone
26 November, 2014 (DH)
Analyze amplitude as a function of azimuth
27 November, 2014 (DH)
Summary
In this introduction we have seen how fracture and stress
introduce velocity anisotropy into the earth.
This anisotropy can be observed as azimuthal traveltime and
amplitude variations in the P-wave seismic data.
In subsequent sections we will
– more rigorously develop the link between anisotropy and fractures
– model the azimuthal traveltime and amplitude response due to fractures
– Quantify the azimuthal amplitude variations with two sets of Amplitude
Versus Azimuth (AVAz) attributes
– Near-offset Rüger equation
– Azimuthal Fourier Coefficients
– Quantify the azimuthal traveltime variations with Velocity Versus Azimuth
(VVAz) attributes
Before doing this, we set up a project in which we read in an
azimuthal pre-stack seismic dataset.
28 November, 2014 (DH)
Workflow for Azimuthal AVO
Exercise 1 Exercise 2B
AVO compliant, Azimuthally
COCA and CACA
preserved, NMO-corrected,
gathers, not used in
CDP gathers.
subsequent steps,
(The gathers contain both
but useful for visual
amplitude and time-delay effects
analysis
due to azimuthal variations)
Exercise 2A
Trim Statics Stack for horizon
(to correct time-delay effects) picking and well
correlation, if needed
Exercises 4 & 5
Azimuthal AVO analysis
Using either Ruger or Fourier
Coefficient approach
Exercise 6
Visualization of results
29 November, 2014 (DH)
Exercise 1: Setting up the project
As a starting point, we load a seismic dataset into the
fracture detection program that has a good
distribution of both offsets and azimuths.
This dataset focuses on the Haynesville Shale in
Texas.
We also read a well from this area.
30 November, 2014 (DH)
Exercise 1: Setting up the project
The Haynesville Formation is
a black, organic-rich shale of
Upper Jurassic age that
underlies much of the Gulf
Coast area of the United
States. The Haynesville
Formation overlies the
Smackover Formation and is
overlain by rocks of the Cotton
Valley Group. Its economic
viability has primarily been a
result of advances in
horizontal well drilling and
hydraulic fracturing.
31 November, 2014 (DH)
Exercise 1: Setting up the project
Before getting to the
theory, let us start the
project.
First start the Geoview
program by clicking the
Geoview icon on your
desktop. When you
launch Geoview, the
first window that you see
contains a list of any
projects previously
opened in Geoview, as
shown here.
32 November, 2014 (DH)
Exercise 1: Setting up the project
For this exercise, we start a new project. Before doing that, you can set all the
data paths to point to the location where you have stored the tutorial data. To do
that, click the Settings and Path tabs:
Now you see the default locations for the Data Directory, Project Directory, and
Database Directory. We change these to point to the directory where the
workshop data is stored.
33 November, 2014 (DH)
Exercise 1: Setting up the project
To change all of the directories to the same location, click the Set all default
directories to button and then click the 3-dot button to the right:
Then, in the File Selection dialog, select the folder which contains the tutorial data
(note that your view will probably be different from what we show below):
34 November, 2014 (DH)
Exercise 1: Setting up the project
After setting all three paths,
the Geoview window now
shows the selected
directories:
When you have finished
setting all the paths, click
Apply to store these paths:
35 November, 2014 (DH)
Exercise 1: Setting up the project
Now click the Projects
tab and choose the
option to create a New
Project.
36 November, 2014 (DH)
Exercise 1: Setting up the project
A dialog appears, where we set
the project name. We will call it
Haynesville, as shown here.
Enter the project name and click
OK on that dialog:
37 November, 2014 (DH)
Exercise 1: Setting up the project
Instead of using the database
below Specify a database.
Select Haynesville_well.wdb
and then click OK.
38 November, 2014 (DH)
Exercise 1: Setting up the project
The Geoview Window now looks like something like this:
39 November, 2014 (DH)
Exercise 1: Loading seismic data
Now that the wells have been loaded, our next
step is to load the seismic volumes.
On the far left side of the Geoview window, click
the Seismic tab:
The window to the right of this tab shows all
seismic data loaded so far.
Go to the bottom of the window (on the left of the
workspace) and click the Import Seismic button.
On the pull-down menu, select From SEG-Y File:
40 November, 2014 (DH)
Exercise 1: Loading seismic data
On the dialog that appears,
select the folder Seismic Data
and then the file
selected_gathers_before_azi_
vel_analysis.sgy and click
Select.
Then, click Next.
41 November, 2014 (DH)
Exercise 1: Loading seismic data
As shown above, click
Next to accept the
Geometry Type.
Then, as shown on the
right, click Next to
accept the General
Information.
42 November, 2014 (DH)
Exercise 1: Loading seismic data
On this dialog, fill in the
options as shown here to
read the azimuths from
the trace headers. Note
that the Data type for
Azimuth has to be
changed from Float to
Integer. Also note that
Calculate under Azimuth
Angles calculates the
azimuth from the
source/receiver
coordinates but these are
meaningless after
migration. Choose
Azimuth Angles Map
North. Click Next when
you have entered the
Map North: azimuth measured clockwise from north options.
Cartesian coordinates: azimuth measured counterclockwise from x-axis (east)
43 November, 2014 (DH)
Exercise 1: Loading seismic data
Now the following warning message appears because the program is about to
scan the entire SEG-Y file:
Click Yes to begin the scanning process.
When the scanning has finished, the Geometry Grid page appears, as
shown in the next slide.
44 November, 2014 (DH)
Exercise 1: Loading seismic data
The geometry as shown
here is correct.
Note, in this workshop
we are only reading 5
inlines out of the 3D
volume due to space
constraints.
Click OK to complete
the process of loading
the seismic data.
45 November, 2014 (DH)
Exercise 1: Loading seismic data
After building the geometry files, a new window appears, showing how the wells
are mapped into this seismic volume:
Since the mapping is correct, click OK to accept the locations shown on this window.
46 November, 2014 (DH)
Exercise 1: Loading seismic data
Now the seismic data appears within the Geoview window (move the slider down
to get the same display):
Before leaving this exercise, we will look at a few annotation options.
47 November, 2014 (DH)
Exercise 1: Loading seismic data
The first thing we would like to do is superimpose a well log curve at its
correct location.
To do this, click the well log icon on
the right side of the menu bar above
the seismic data.
Click the down arrow beside the icon
and select one of the wells to display
(here we select well-A).
48 November, 2014 (DH)
Exercise 1: Loading seismic data
A new display is then created at the location of the well (Inline 190, crossline 51)
with the sonic log from the well inserted. Scroll to the left to create the display
shown below:
49 November, 2014 (DH)
Exercise 1: Loading seismic data
Next, we turn on the
azimuth display.
First, click on the down
arrow next to the “eyeball”
icon and select View 1
Attributes, as shown
below:
Next, on the new dialog,
click Scale Annotation,
and turn on the Azimuth
option, as shown to the
right. Finally, click OK:
50 November, 2014 (DH)
Exercise 1: Loading seismic data
The azimuth values are
annotated at the top of the
plot. Data sharing common
azimuths are shaded in
alternating yellow and gray
colors. This helps
emphasize azimuthal
variations in the traveltime
or amplitude.
This data was migrated
using a series of azimuthal
migrations so looks like a
Common Offset / Common
Azimuth (COCA) gather,
however the offset sampling
is irregular.
51 November, 2014 (DH)
Exercise 1: Rose Diagram Analysis
The Rose Diagram
Analysis can be used to
determine if the dataset
has adequate azimuthal
and offset sampling to
perform the Azimuthal
analysis.
On the Processes Tab
under the Project
Manager and Azimuthal
AVO select the Rose
Diagram Analysis. This
initiates the window
shown to the right.
52 November, 2014 (DH)
Exercise 1: Rose Diagram Analysis
You may get the following message, telling you that your license is not active. If
you do, click OK and do the following steps: Under the Start tab select License,
and then click in the box to the right of HRS-proaz to turn on the AVAz option.
53 November, 2014 (DH)
Exercise 1: Rose Diagram Analysis
A Rose Diagram is a circular histogram which
plots directional data. It is an azimuthal
histogram plotting the frequency of each
azimuth for some bin size.
The Rose Diagram Analysis supports three
types of Analysis.
Well log
Trace Header
Seismic attributes
In this exercise we are interested in the
azimuth distribution of the input data so select
the second option – Seismic Volume…
Then select the input SEG-Y data
Selected_gathers_before_azi_vel_analysis
and then Scan the data.
54 November, 2014 (DH)
Exercise 1: Rose Diagram Analysis
A display similar to the right should be
displayed. The size of the red triangles
show the number of traces in each
azimuth bin.
This dataset is azimuthally sectored
migrated and as such is nicely
organized as a function of azimuth.
There are azimuths at 0, 30, 60, 90,
120 and 150 degrees. There are more
traces at 30 degrees azimuth than for
the 90 degrees azimuth.
The red histogram conveys no
information about the offsets. This
information is conveyed by the 1D
histogram to the bottom right outlined in
blue. In order to see how the offset
distribution varies with azimuth, select
the Histogram Parameters.
55 November, 2014 (DH)
Exercise 1: Rose Diagram Analysis
It is also possible to display a 2D polar
histogram based on two attributes such
as azimuth and offset.
By selecting the Histogram Parameters
tab the parameters specifying the
histogram can be modified. The user
has control over the number of
azimuths to display and their range.
By deselecting Rose the 2D polar
histogram becomes more obvious. The
offset radiates from the origin. For
each offset and azimuth the grey
shading indicates the number of traces
(fold) in this bin.
This display makes it clear that the
maximum offset in the East – West
direction is much smaller than in the
North-South direction.
56 End of exercise 1 November, 2014 (DH)
The CoffCaz display
The initial gathers that were displayed were a function of both
offset and azimuth and may be irregularly sampled.
When we add the extra dimension of azimuth, it is difficult to
visualize both the offset and azimuth dimensions at once.
For that reason, the next exercise introduces the concept of
the common offset, common azimuth, or CoffCaz, display.
This leads us to velocity analysis versus azimuth (VVAz) and
amplitude analysis versus azimuth (AVAz).
The CoffCaz display involves displaying subsets of increasing
azimuth within increasing offset panels, using a rolling mix to
improve the signal-to-noise ratio.
This is shown next, and we can notice anisotropic effects.
57 November, 2014 (DH)
Observations of anisotropy in pre-stack seismic data
180o
Az
0o
14000’ Increasing Offset 2000’
58 November, 2014 (DH)
After Gray et al. (2004)
Observations of anisotropy in pre-stack seismic data
180o
Az
0o
14000’ Increasing Offset 2000’
40ms
Fast
Slow
59 November, 2014 (DH)
After Gray et al. (2004)
Azimuthal traveltime variations
As seen in the previous slide, we observe anomalous normal
moveout (NMO) patterns within the azimuthally varying
component of the gathers.
This NMO can be used in several ways:
To analyze the data for azimuthally anisotropic velocity.
To use these velocities to correct the gathers.
To use these velocities for fracture analysis.
In the next few slides, we show how to correct the gathers
using azimuthal velocity analysis.
Later in the course we discuss how to use this velocity
analysis for VVAz to determine fracture patterns.
60 November, 2014 (DH)
Azimuthal NMO
In order to account for azimuthal traveltime variations the near-offset
NMO equation can be modified for elliptical azimuthal velocity
variations (Tsvankin, 2001):
N
x2
t t0 2 , where :
2 2 f
Vnmo
1 cos2 f sin 2 f
2
2
2
E
Vnmo Vslow V fast
We fit these elliptical NMO velocities by analyzing the data as a
function of offset and azimuth, as shown in the next slides.
61 November, 2014 (DH)
Correcting the gathers
Our first step is to use these velocity
calculations to correct the input data, which
involves:
Calculate continuous time variant statics
to estimate residual traveltimes.
For each analysis time, fit an elliptical
residual NMO surface to residual
traveltimes
From this calculate the elliptical stacking
velocities and azimuth:
– Fast stacking velocity: Vfast
– Slow stacking velocity: Vslow
– Azimuth of slow stacking velocity:
Apply elliptical stacking velocities to
correct the data.
An example is shown in the next figures, first
showing the application of the corrections on
the gathers and then showing the stacks.
62 November, 2014 (DH)
COCA Prior Azimuthal NMO Correction (shale gas example)
OFFSET
Az
N
Prior
info Vfast
Vslow Az. NMO
Final
Velocity
63 November, 2014 (DH) 63
COCA After Azimuthal NMO Correction (shale gas example)
OFFSET
Az
N
Prior
info Vfast
Vslow
Final
Velocity
64 November, 2014 (DH)
NMO STACK
65 November, 2014 (DH)
Az NMO STACK
66 November, 2014 (DH) 66
Correcting the gathers
In our software, an azimuthal NMO correction has not yet
been implemented.
However, the correction can be done quite effectively using
the trim static correction.
In the next exercise, we will flatten the gathers to prepare
them for Azimuthal NMO analysis.
We will create Common Offset, Common azimuth (COCA)
stacks before and after this correction for comparison to
model data.
67 November, 2014 (DH)
Exercise 2: Trim statics, Stacking and Gathering
We now return to
the Haynesville
project and perform
time variant trim
statics to flatten the
data.
From the
Processes tab,
select Trim Statics
under Seismic
Processing -
Utilities.
On the dialog to the
right, select
selected_gathers_
before_azi_vel_an
alysis as input.
68 November, 2014 (DH)
Exercise 2: Trim statics
Trim Statics is the process which corrects for
residual moveout errors and aligns the events on
the gathers.
In the Trim Statics process, a pilot trace is formed
by stacking each CDP gather.
Then, each gather trace is correlated with the pilot
trace, using a series of sliding windows.
The cross correlations are used to calculate an
optimal time shift for that window.
Finally, the shifts for the windows are interpolated
to produce a time-variant stretch of the trace.
The result is to align events with the pilot trace.
69 November, 2014 (DH)
Exercise 2: Trim statics
The dialog which appears shows
that the default is to apply the
process to the entire volume and
create a new output volume
called trim_statics:
We will use a series of windows
of size 100 ms, window step of
50 ms, with a maximum shift of
20 ms:
When you have filled in the
menu as shown above, click
Run to start the process.
70 November, 2014 (DH)
Exercise 2: Trim statics
The result shows a much improved alignment at the target zone around 2280 ms:
71 November, 2014 (DH)
Exercise 2: Stacking
Now we create a
CDP stack from
which we can pick a
Horizon.
From the
Processes tab,
select CDP Stack,
as shown
On the dialog to the
right, select
trim_statics as the
Input and click OK.
72 November, 2014 (DH)
Exercise 2: Stacking
The resulting stack is
then displayed in View
2, as shown on the
right.
73 November, 2014 (DH)
Exercise 2: Stacking
Next, we pick a horizon for later use. To do this, first select Pick Horizons from
the Horizon tab at the top of the seismic display.
Then, use Select View: to select View 2 and click OK.
When prompted whether to display a Map View, click Yes.
74 November, 2014 (DH)
Exercise 2: Stacking
Two things now happen. First, you see the map view in a separate window:
75 November, 2014 (DH)
Exercise 2: Stacking
Second, the picking dialog appears below the seismic window. Change the
Mode to Left & Right Repeat and keep the search as Peak. Then click the
peak around 2240 ms, as seen below:
The peak on this event on the stack is auto-picked with orange plus signs. The
yellow dashed line that appears on the gathers to the left shows the constant
picks from the stack but indicates that pre-stack picking has not been done.
76 November, 2014 (DH)
Exercise 2: Stacking
Next, select Automatic
Picking from the
Options pull-down, and
click OK:
77 November, 2014 (DH)
Exercise 2: Stacking
The auto-picked time structure map is then filled in, as shown here. Click the x
in the upper right part of the map to close this window.
78 November, 2014 (DH)
Exercise 2: Stacking
Then, click OK on the picking dialog to close this dialog and save the picks:
79 November, 2014 (DH)
Exercise 2: Stacking
The resulting
display now looks
like this, where the
picked event is in
blue.
80 November, 2014 (DH)
Exercise 2: Snail Display
For COV migrated
data, the data may be
sorted in manner
where it is not easy to
observe azimuthal
variations. The snail
display sorts the data 7 8 9 10
from the origin in a 6 1 2 11
spiral like a snail.
5 4 3
The number of traces
going into the
process is the same
as output, so the
snail display is just
for display purposes.
81 November, 2014 (DH)
Exercise 2: Snail Gather
Our dataset is well sorted but we will run the snail
gathers for illustration purposes.
This is done by selecting the snail gather process
from the Processes Menu Snail Gathers.
Choose the following the parameters and then run
the process.
82 November, 2014 (DH)
Exercise 2: Snail Gather
The resulting display
shows the input to
the left and the snail
gather to the right.
Change the display
so the gathers are
displayed vertically.
The snail gather is
now displayed on the
bottom for easy
comparison.
Later we will use the
snail display in the
azimuthal analysis.
83 November, 2014 (DH)
Exercise 2: Common Offset/Azimuth Gathers
Our next analysis step is to produce the
Common Offset / Common Azimuth
(CoffCaz) Gather. This is not actually used
in any further analysis step, but is an
important display to analyze the data.
To do this, double-click on Common
Offset/Azimuth Gathers from the
Azimuthal AVO (ProAZ) option.
84 November, 2014 (DH)
Exercise 2: Gathering
To save time and disk space we will run only a
Single Inline, 190 is to be run.
The CoffCaz gather organizes the data by
offset and azimuth To limit the size of the
display, only display to a maximum offset of
15000 ft, with 15 offset bins.
Choose the Number of Azimuth Ranges to be
6 starting at 0.0 degrees to match the azimuth
distribution observed in the Rose diagram
analysis. The resulting azimuth bins are
echoed back. By invoking reciprocity only
azimuths from 0 to 180 degrees are displayed
instead of 0 to 360 degrees.
In this run, we primarily sort by offset, then by
azimuth.
Click OK to run the process:
85 November, 2014 (DH)
Exercise 2: Gathering
The CoffCaz gather is
then displayed in View 2,
as shown:
The regularization makes
azimuthal variations more
obvious. It also
emphasizes the near
offset sampling issues.
The 3X3 super binning
also improves the S/N
ratio.
This result will be
compared to the
modeling later on. For
the actual analysis it is
better to use the raw trim
gathers.
86 November, 2014 (DH)
Exercise 2 – Angle of Incidence
Prior to performing
the azimuthal
analysis we need
to understand the
range of the angles
of incidence within
the dataset.
Turn off view 1 in
order to expand
the CoffCaz gather
to the full view.
87 November, 2014 (DH)
Exercise 2 – Angle of Incidence
Right click on the
gather, select
Color Data
Volume, Incident
Angle
This brings up the
Velocity Model Dialog.
Choose New, specify
type Single Well: P-
wave Curve, then
choose the velocity
curve EDIT_DTC_corr
and modify the vertical
smoother to 300 ft. and
then press OK.
88 November, 2014 (DH)
Exercise 2 – Angle of Incidence
At the zone
of interest
(Horizon 1)
there are
angles up to
30 degrees.
89
End of exercise 2 November, 2014 (DH)
Anisotropy
Anisotropy may arise due to a number of reasons
•Intrinsic anisotropy : minerals and crystals exhibit
anisotropy due to their regular molecular structure. Clay
minerals are anisotropic.
• Fine layers: anisotropy occurs when the seismic
wavelength is several orders of magnitude larger than the
thickness of individual strata so the seismic sees an
average or effective layer property (e.g. Backus average).
• Fracture induced anisotropy: Fractures can be thought of
as weak layers in an otherwise homogenous material again
giving rise to layer based anisotropy.
• Stress induced anisotropy: for an isotropic rock
composed of randomly orientated microcracks the stress
field can preferentially close certain fractures giving rise to
anisotropy (Nur and Simmons, 1969).
90 November, 2014 (DH)
Generalized stress-strain relation
All of our theory in this course 1
c11 c12 c13 c14 c15 c16 e1
follows from the relationship
2
c c22 c23 c24 c25 c26 e2
between stress and strain, 12
called Hooke’s Law. 3
c13 c23 c33 c34 c35 c36 e3
,
The stress-strain relation can 4 c14 c24 c34 c44 c45 c46 e4
be written in Voigt matrix 5 c15 c25 c35 c45 c55 c56 e5
vector notation =Ce.
6 c16 c26 c36 c46 c56 c66 e6
The stiffness's are a
symmetric matrix while the 1 11 e1 11
stress and strains are vectors. e
2 22 2 22
The stress and strain tensors 3 33 e3 33
in Voigt notation are , .
4 23 e4 2 23
5 13 e5 213
6 12 e6 212
91 November, 2014 (DH)
Stress and strain in 3D
The nine components of stress (ij) are shown here on a
rock cube (red = compressive, blue, orange, grey = shear):
However, note that the
shear components ij =
ji, or the cube would
rotate.
Thus, there are only
three unique shears
(12,13, and 23) and
three compressions
(11,22, and 33).
There are also six unique
strains ij.
92 November, 2014 (DH)
The stiffness matrix
The stiffness matrix C shown c11 c12 c13 c14 c15 c16
c c c c c c
here is from the most general 12 22 23 24 25 26
anisotropic case, triclinic. c13 c23 c33 c34 c35 c36
C
The coefficients cij allow us to c c c c c
14 24 34 44 45 46 c
calculate both seismic c15 c25 c35 c45 c55 c56
velocities and amplitudes.
c16 c26 c36 c46 c56 c66
In most of the course we consider three simpler cases derived
from the general stiffness matrix: isotropy, vertical transverse
isotropy (VTI) and horizontal transverse isotropy (HTI).
In each case, we show how the velocity is derived from the
stiffness matrix.
We also discuss Thomsen’s parameters, which are useful in
weakly anisotropic materials.
93 November, 2014 (DH)
Velocity related to stiffness matrix
The diagonal elements of the stiffness matrix define the P and S-wave
velocities traveling along the coordinate axes.
The velocity Vij is propagating in the ith direction and is polarized in the jth
direction.
P-wave velocities S-wave velocities
c11 c44
V11 V23 V32
r r
c22 c55
V22 V13 V31
r r
c33 c66
V33 V12 V21
r r
The S-wave may be either a SV or SH wave depending on the direction the
wave is propagating in.
94 November, 2014 (DH)
Isotropic Sandstone
In an isotropic x1
material, like this x2
q
sandstone sample,
x3 V(90 o
)
parameters such as
velocity are the
same in all
directions. Note that
the x2 direction is o
outward from this V(0 )
plane, and we only
consider velocity in
the x1 and x3 Isotropic sandstone, from Google images.
directions.
In this case V(0o) = V(90o), where the x3
o
95
direction has an angle q 0 .
November, 2014 (DH)
Isotropic stiffness
For isotropic rocks like our sandstone, there are only two
independent coefficients in the stiffness matrix, c11 and c44:
c11 c11 2c44 c11 2c44 0 0 0
c 2c c c 2 c 0 0 0
11 44 11 11 44
c11 2c44 c11 2c44 c11 0 0 0
C
0 0 0 c44 0 0
0 0 0 0 c44 0
0 0 0 0 0 c44
These two terms can also be written as combinations of
the Lamé coefficients l and m, where c11 = l + 2m, c44 = m
and c11 - 2c44 = l.
96 November, 2014 (DH)
Velocity related to stiffness matrix
As discussed earlier, the diagonal elements of the stiffness matrix
define the P and S-wave velocities traveling along the coordinate
axes.
We also saw that, for anisotropic rocks, the S-wave may be either
a SV or SH (Sfast or Sslow) wave, depending on which symmetry
plane it is travelling in.
However, for isotropic rocks where:
c11 = c22 = c33 = l + 2m, and c44 = c55 = c66 = m,
there are only two velocities:
l 2m
P - wave velocity VP ,
r
m
and S - wave velocity VS .
r
97 November, 2014 (DH)
Anisotropic models
We next look at two different types of anisotropy: Vertical transverse
isotropy (VTI) and horizontal transverse isotropy (HTI).
The figure below illustrates the difference between the two types.
Notice that the VTI model
can be transformed to the
HTI model simply by
rotating from the z
direction to the x
direction.
Rüger, 2002
The VTI model consists of horizontal layers and can be extrinsic, caused by fine
layering of the earth, or intrinsic, caused by particle alignment as in a shale.
The HTI model consists of vertical layers and is caused by parallel vertical
fractures or steeply dipping shale.
98 November, 2014 (DH)
Anisotropic Shale
x1 Consider this
x2 anisotropic shale
outcrop. This is an
x3
example of a VTI
V(90o) material.
The velocity is faster
in the direction parallel
1m to layering, so we find
that V(90o) > V(0o).
V(0o)
As in the isotropic
case, the velocity is
the same in all
Utica shale in outcrop, horizontal directions.
Courtesy Google images.
99 November, 2014 (DH)
VTI case: P-wave velocities
For the vertical transversely isotropic (VTI) case, as in our
shale, there are five independent stiffness coefficients:
c11 c11 2c66 c13 0 0 0
c 2c c c 0 0 0
11 66 11 13
c13 c13 c33 0 0 0
C
0 0 0 c44 0 0
0 0 0 0 c44 0
0 0 0 0 0 c66
This gives rise to the following P-wave velocities traveling along
the two principal directions (vertical and horizontal):
c11 c33
VP (0o ) ,VP (90o )
r r
100 November, 2014 (DH)
VTI case: S-wave velocities
There are two polarizations for the S-wave velocities in the two principal
directions (horizontal and vertical):
c55 c66
VSV , and VSH .
r r
As discussed in Appendix 1, we can use the Kelvin-Christoffel
equations to compute the P and S-wave phase velocities in arbitrary
directions.
However, a more intuitive way of understanding the velocities was
given by Thomsen (1986).
101 November, 2014 (DH)
Thomsen parameters (VTI media)
For VTI media Thomsen (1986) described a more intuitive
parameterization based on the perturbations , d and g,(now called the
Thomsen parameters) given by:
c11 c33 VP (90o ) VP (0o )
o
,
2c33 VP (0 )
c66 c44 VSH (90o ) VSH (0o )
g o
, and
2c44 VSH (0 )
d
c13 c55 c33 c55
2 2
VP ( 45o ) VP (0o )
4 .
2c33 c33 c55 o
VP (0 )
Notice that we can think of and g as deviations of the horizontal
and vertical P and SH velocities, and d as including the 45o P-wave
deviation.
102 November, 2014 (DH)
Velocities for the VTI case
Here are Thomsen approximations of the P, SV and SH phase
velocities as a function of his parameters:
VP q VP 0 1 d sin2 q cos2 q sin4 q ,
VP20
VSV (q ) VS 0 1 2 ( d ) sin q cos q ,
2 2
VS 0
VSH q VS 0 1 g sin2 q , where :
c33 c55
VP 0 and VS 0 .
r r
The next figures shows how close Thomsen’s approximation is to the
correct phase velocities found using the Kelvin-Christoffel equations
in Appendix 1.
103 November, 2014 (DH)
Typical Values for , d, and g
Typical values for , d, and g were given by Thomsen (1986).
Here are some representative values from his table (we
model the highlighted mudshale):
Lithology VP(m/s) VS(m/s) rho(g/cc) epsilon delta gamma
sandstone_1 3368 1829 2.50 0.110 -0.035 0.255
sandstone_2 4869 2911 2.50 0.033 0.040 -0.019
calcareous sandstone 5460 3219 2.69 0.000 -0.264 -0.007
immature sandstone 4099 2346 2.45 0.077 0.010 0.066
shale_1 3383 2438 2.35 0.065 0.059 0.071
shale_2 3901 2682 2.64 0.137 -0.012 0.026
mudshale 4529 2703 2.52 0.034 0.211 0.046
clayshale 3794 2074 2.56 0.189 0.204 0.175
silty limestone 4972 2899 2.63 0.056 -0.003 0.067
laminated siltstone 4449 2585 2.57 0.091 0.565 0.046
104 November, 2014 (DH)
Mesaverde mudshale
Note that the observed parameters for the Mesaverde
mudshale on the previous table where VP0 = 4.529 km/s,
VS0 = 2.703 km/s, r = 2.52 g/cc, = 0.034, g = 0.046, and
d = 0.211. The equivalent stiffness matrix is given by:
55.21 14.99 24.41 0 0 0
14.99 55.21 24.41 0 0 0
24.41 24.41 51.69 0 0 0
C
0 0 0 18.41 0 0
0 0 0 0 18.41 0
0 0 0 0 0 20.11
105 November, 2014 (DH)
Mesaverde mudshale
Here is the P-wave phase velocity as a function of angle q,
where the true value is shown as the solid line and
Thomsen’s approximation is shown as the dotted line:
106 November, 2014 (DH)
Mesaverde mudshale
Here is the SH-wave phase velocity as a function of angle
q, where the true value is shown as the solid line and
Thomsen’s approximation is shown as the dotted line:
107 November, 2014 (DH)
Mesaverde mudshale
Here is the SV-wave phase velocity as a function of angle
q, where the true value is shown as the solid line and
Thomsen’s approximation is shown as the dotted line:
108 November, 2014 (DH)
Vertically Fractured Limestone
V(90o)
V(0o)
(Reches, 1998)
Here is a vertically fractured limestone in vertical cross-
section view, where we would expect the fast velocity to
parallel the fractures, so that:
o o
V(0 ) > V(90 ).
This type of anisotropy is called horizontally transverse
isotropy,
109
or HTI. November, 2014 (DH)
Horizontal Transverse Isotropy
A better way to look at this rock is in an idealized three-dimensional view,
as shown below, where the vertical “sheets” represent the fractures:
After Rüger (1997)
We define the symmetry-axis plane at right-angles to the
fractures (f = 0o) and the isotropy-plane as parallel to
fractures (f = 90o).
110 November, 2014 (DH)
Horizontal Transverse Isotropy
Horizontal transverse isotropy (HTI) is described by 5 independent
parameters:
c11 c13 c13 0 0 0
c c c 2 c 0 0 0
13 33 33 44
c13 c33 2c44 c33 0 0 0
CHTI
0 0 0 c44 0 0
0 0 0 0 c55 0
0 0 0 0 0 c55
As shown in Appendix 2, an HTI
medium can be thought of as a VTI
medium rotated so that the symmetry
axis is in the x direction.
Thus, we can use the same
Thomsen parameters, but write them
as (V), d(V) and g(V) to indicate the
vertical transformation. Rüger, 2002
111 111 November, 2014 (DH)
Second fractured limestone
This second fractured
limestone displays
fractures in both the
horizontal and
vertical direction.
This type of fracturing
leads to orthorhombic
anisotropy, which is
discussed in
Appendix 3.
Fractured limestone outcrop, Musandam,
Oman: Courtesy Google images
112 November, 2014 (DH)
Orthorhombic media
Orthorhombic media are described by 9 independent parameters
c11 c12 c13 0 0 0
c c c 0 0 0
12 22 23
c13 c23 c33 0 0 0
C ,
0 0 0 c44 0 0
0 0 0 0 c55 0
0 0 0 0 0 c66
These media are characterized by 3 mutually orthogonal planes of symmetry
x2 x2
x1 x1
Orthogonal Identical Vertical fractures embedded in a
fractures fractures horizontally-layered shale
113 November, 2014 (DH)
Monoclinic media
There are more complex forms of anisotropy such as
monoclinic media. In this case there is only 1 mirror symmetry
plane. Monoclinic media might arise due to 2 sets of unequal
vertical fractures. Monoclinic media are described by 13
parameters.
c11 c12 c13 0 0 c16
c c22 c23 0 0 c26
12
c13 c23 c33 0 0 c36
Cmono
0 0 0 c44 c45 0
0 0 0 c45 c55 0
2 sets of vertical fractures c16 c26 c36 0 0 c66
114 November, 2014 (DH)
Fracture induced anisotropy
– The following is short summary of the rock physics of fractured rocks
designed to familiarize the interpreting geophysicist with the concepts
needed to perform forward modeling. More extensive material may be
found in Appendices 5 & 6.
– We model the anisotropy introduced by the fractures using linear slip
deformation theory (Schoenberg, 1980). The theory is written in terms
of compliance S which is the inverse of stiffness C.
– In this theory the fractures perturb the background rock acting as a
source of extra strain.
– Mathematically the total compliance of the fractured rock is the sum of
the background compliance plus the extra compliance due to the
fractures
S Sb S f
where :
Sb is the isotropic background compliance matrix
Sf is the excess compliance due to fractures
115 November, 2014 (DH)
LS fracture compliances
For parallel rotationally invariant fractures the
form of the fracture compliance matrix is
particularly simple.
It is described by two fracture compliance
relationships and parameters
N BN N and T BT T , where :
BN normal fracture compliance , T
T
BT tangential fracture compliance .
N N
Verdon and Wüstefeld, (2013) summarize the T T
published BN/BT ratio’s (See Appendix 6)
The median BN/BT ratio is 0.65
Vertical fractures lead to HTI symmetry
116 November, 2014 (DH)
LS fracture weaknesses
Rather than working with fracture compliances, it is easier to work the
dimensionless fracture weakness parameters.
These are fractional values that range from 0 to 1.
• At 0 the fracture has no influence on the background rock.
• At 1 the moduli of the resultant rock goes to zero.
The tangential fracture compliance becomes the tangential fracture weakness.
• Function of the crack density
mBT
• Can estimate this from gamma dT , d T [0,1)
1 mBT
The normal fracture compliance becomes the normal fracture weakness.
Function of the crack density and fluid fill
dN
l 2m BN , d [0,1)
1 l 2m BN
N
Penny shaped cracks can be written in terms of fracture weaknesses.
117 November, 2014 (DH)
Fracture rock physics: Vertical fracture
This gives rise to the HTI stiffness matrix
M b 1 d N lb 1 d N lb 1 d N 0 0 0
l 1 d M 1 2
dN lb 1 bd N 0 0 0
b N b b
l 1 d N lb 1 bd N M b 1 b d N
2
0 0 0
C S 1 b
0 0 0 mb 0 0
0 0 0 0 mb 1 d T 0
0 0 0 0 0 mb 1 d T
where : M b lb 2 mb ,
b lb / M b
118 November, 2014 (DH)
LS stiffness matrix
Although the LSD stiffness matrix looks quite daunting, several
simple observations can be made about it.
First, note that if dN and dT equal zero, it reduces to the isotropic
case (this is for the dry case, the saturated case is more complex).
Second, note from the form of the normal and tangential
weaknesses, that dN and dT equal zero when BN and BT equal zero
and approach 1 as BN and BT become very large.
Finally, since an HTI medium can be thought of as a rotated VTI
medium, these parameters can be related to Thomsen’s
parameters.
119 November, 2014 (DH)
Relationship between Thomsen parameters
and normal weaknesses
For vertical rotationally invariant vertical fractures Bakulin et al.
(2000) derived the following linearized relations between the
Thomsen parameters and normal and tangential weaknesses
V 2 g 1 g d N
d V 2 g d T d N
V 1
g dT
2
where g Vs / Vp and 1 - 2g
2
Note there is a simple relationship between g and dT so the two
variables can be used interchangeably.
There is an advantage to this, since g can be determined from
sonic scanner logs, thus providing an objective source for the
fracture parameters.
120 November, 2014 (DH)
Modeling - Elastic Model
How to construct anisotropic
stiffness matrix?
– For AVO modeling we need three pieces of
information: density, P-wave and S-wave
velocity.
– For HTI stiffness matrix we need seven
pieces of information: five stiffnesses,
density and azimuth.
Haynesville
121 November, 2014 (DH)
Elastic Model
Introduce tangential and normal
fracture weakness.
– Tangential fracture weakness is 2X
S-wave anisotropy (gamma)
– Borehole sonic log
– Initial guess of 0.2
– Normal weakness from BN/BT ratio
– Verdon and Wüstefeld, (2013)
summarized the published BN/BT
ratio’s
– Fluid filled fractures have low
BN/BT ratio
– Initial guess of BN/BT =0.1
Haynesville
From this calculate the stiffness
matrix for each layer
122 November, 2014 (DH)
Elastic Model
Haynesville
From this we can calculate the stiffness matrix and Thomsen Parameters
123 November, 2014 (DH)
Azimuthal Modeling
Once the elastic parameters
are specified then the forward
modeling can be run
The azimuthal modeling is
based on convolutional
modeling
– The reflectivity is calculated using a
linearized approximation of the
anisotropic Zoeppritz equation.
– The reflectivity is then convolved
with a user defined wavelet similar
to AVO modeling.
– The synthetic can be generated
using a variety of different
geometries.
– The synthetic is output as a seismic
gather which can then be compared
to the seismic in the AVAz analyzer.
124 November, 2014 (DH)
Azimuthal Modeling
The synthetic can be generated with (left) or without residual azimuthal traveltime
variations. The fractured carbonate introduces azimuthal traveltime variations for all
subsequent times.
125 November, 2014 (DH)
AVO and AVAz analysis
In this next theory section, we discuss Amplitude versus
Offset (AVO) and Amplitude versus Azimuth (AVAz)
analysis and show how AVAz leads naturally to fracture
analysis.
For isotropic and vertically transverse isotropic (VTI)
media, the data is a function of incident angle q alone.
However, for horizontally transverse isotropic (HTI)
media the data is a function of both incident angle q and
azimuth angle f.
Although this section focusses mainly on the AVAz
analysis of HTI media using Rüger’s equation, we start
with a discussion of AVO in both isotropic and VTI
media.
126 November, 2014 (DH)
AVO in isotropic and anisotropic media
This figure shows what happens when an incident P-wave strikes the boundary
between two VTI layers. When the d and coefficients equal zero, this simplifies
to the isotropic case.
127 November, 2014 (DH)
The Zoeppritz Equations
Zoeppritz derived the amplitudes of the reflected and
transmitted waves using the conservation of stress and
displacement across the layer boundary, which gives four
equations with four unknowns. Inverting the matrix form of the
Zoeppritz equations gives us the exact amplitudes as a function
of angle for isotropic media:
1
sin q1 cos f1 sin q 2 cos f2
RPP cos q sin f1 cos q 2 sin f2 sin q1
R 1 cos q
PS sin 2q r 2VS 22VP1 r 2VS 2VP1
cos 2f2
VP1 1
cos 2f1 cos 2f1
TPP 1
VS1 r1VS1VP 2
2
r1VS1 2 sin 2q1
r 2VP 2 r 2VS 2
TPS cos 2f1
VS 1
sin 2f1 cos 2f2 sin 2f2 cos 2f 1
VP1 r1VP1 r1VP1
128 November, 2014 (DH)
Anisotropic Zoeppritz equation
Daley and Hron (1977, 1979) published the first papers
on AVO in transversely isotropic media. Their second
paper treats ellipsoidally anisotropic media, where d .
Schoenberg and Protázio (1992) extended the
Zoeppritz equations to orthorhombic anisotropy (see
Appendix 4).
The Zoeppritz equations are nonlinear and therefore
difficult to understand.
We work with a linearization of Zoeppritz equations,
starting with the isotropic case.
129 November, 2014 (DH)
The isotropic Aki-Richards Equation
The Aki-Richards equation is a linearized approximation to the Zoeppritz
equations, and was rewritten by Wiggens et al (1984) in a form better
suited to reflection seismology. (Note that the velocity and density terms
are averages and differences, the angle is also the average, and this is
for isotropy, as indicated by the subscript iso):
Riso (q ) Aiso Biso sin 2 q Ciso tan 2 q sin 2 q
where:
1 VP r
Aiso
2 V p r
2 2
1 VP V VS V r
Biso 4 S 2 S
2 Vp VP VS VP r
1 VP
Ciso
2 Vp
130 November, 2014 (DH)
Isotropic AVO
To solve this equation, we pick q
q1 qN
the amplitudes of an N-trace
600
angle gather at each time t, as
Time t
shown on the right. For angles (ms)
less than about 30o, we can 650
drop the third term, to get:
R(qi ) Aiso Biso sin 2 qi , i 1, , N , where :
R(qi ) the reflection coefficien t at incident angle qi ,
A the intercept, and B the gradient.
As shown in Appendix 7, we can write this set of N linear
equations in matrix form and solve for Aiso and Biso.
131 November, 2014 (DH)
Linearized AVO in VTI media
Thomsen (1993) showed that VTI terms could be added to the Aki-Richards
equation using his weak anisotropic parameters d and , where R(q ) is the
AVO response and Riso(q ) is the isotropic AVO response.
Rüger (2002) gave the following form of Thomsen’s original equation:
d
R(q ) Riso (q ) sin q
2
sin 2 q tan 2 q ,
2 2
where : d d 2 d1 , and 2 1.
d 2 2
or : R(q ) Aiso Biso sin q Ciso sin q tan q
2
2 2
Note there are 5 unknowns with 3 basis functions so the problem is
underdetermined. Thus, only the following average parameters can
be solved uniquely:
d
Aiso , B Biso , C Ciso
2 2
132 November, 2014 (DH)
AVO and VTI
Blangy (1997) computed the effect of anisotropy on VTI models of the three
Rutherford-Williams type. Blangy’s models are shown below, but since he used
Thomsen’s formulation for the linearized approximation, his figures have been
recomputed in the next slide for the wet and gas cases using Rüger’s formulation.
The slide after that shows our example.
133 November, 2014 (DH)
VTI – AVO Effects
Class 1
Class 1
d = -0.15
= -0.3
Class 2
Class 2
Class 3
Class 3
Isotropic
--- Anisotropic
(a) Gas sandstone case: Note (b) Wet sandstone case:
that the effect of d and is Note that the effect of d and
to increase the AVO effects. is to create apparent AVO
decreases.
134 November, 2014 (DH)
In the Haynesville example the synthetic and data
have different AVO trends
135 November, 2014 (DH)
VTI background
Similar to Lin and Thomsen (2013), VTI anisotropy AVO Response
was introduced into the background shale in the
model to try and match the AVO trend observed in
the seismic.
• In order to match the AVO trend the following VTI
parameters where chosen for the shale background
d = 0.2 key variable used to match AVO trend
= 0.3 using constraint that 𝜀 > 𝛿, (Vernik and
Liu, 1997)
g = 0.3 using constraint that 𝛾 ≈ 𝜀 (Wang,
2002)
This leads to orthorhombic anisotropy in the
fractured reservoir.
136 November, 2014 (DH)
Linearized AVO in HTI media
Next we discuss AVO in HTI media, and AVAZ (Amplitude versus
Azimuth). Rüger (1998) linearized the Zoeppritz equation for the case
of an isotropic half-space over an HTI half-space.
The linearization can be extended to the case of two HTI half-spaces
provided their symmetry axes are aligned.
Recall, the symmetry-axis is perpendicular to the fractures and the
isotropy plane is parallel to the fractures:
After Rüger (1997)
137 November, 2014 (DH)
Azimuth angle
In addition to the raypath angle q, we now introduce an
azimuth angle f, which is defined with respect to the
Isotropy plane: z
x (North)
q f
y
Note that the isotropy plane has the same orientation as the
fractures, thus the azimuth angle f is equal to 0 degrees
parallel to the fracture and 90 degrees perpendicular to the
fractures.
138 November, 2014 (DH)
Linearized AVO in HTI media
With this definition of azimuth angle, we can derive the following
linearized modeling equation for AVO in HTI media:
R(q , f ) Aiso ( Biso Bani sin 2 f ) sin 2 q
(Ciso
1
2
d V cos2 f V sin 2 f sin 2 f ) sin 2 q tan 2 q ,
where Aiso , Biso , and Ciso are the isotropic AVO terms,
1 VS
2
Bani d 8 g
(V ) (V )
2 VP
is the anisotropic gradient, with :
(V ) Thomsen' s parameter defined with respect to vertical,
d (V ) Thomsen' s d parameter defined with respect to vertical,
g (V ) Thomsen' s g parameter defined with respect to vertical,
139
q incidence angle, and f azimuth angle.
November, 2014 (DH)
Linearized AVO in HTI media
In the VTI and HTI AVO expressions given by Rüger (2002),
he rewrites the gradient term B in the AVO equation in terms
of VP and m, as shown below:
Regular form:
1 VP
2
VS VS 2
VS r 1 VP VS 2 VS r
2
Biso 4 2 4
2 Vp VP VS VP r 2 Vp
VP VS r
Rüger form:
2
2
1
1 V V 2 VS r 1 VP V m 2
Biso P 4 S 4 S 2 ln ln r
2 Vp VP VS r 2 Vp VP r
1 VP VS
2
1 VP VS m
2
4 ln m 4
2 Vp P
V 2 p
V P
V m
140
November, 2014 (DH)
Linearized AVO in HTI media
To show the effects of HTI, Rüger (2002) created the
following four models:
Model Vp/Vp A m/m d (V) (V) g (V)
A 0.1 0.1 0.2 0 0 -0.1
B 0.1 0.1 0.2 -0.1 0 0
C 0.1 0.1 0.2 0 -0.1 0
D 0.1 0.1 0.2 -0.05 -0.05 -0.15
m 2 VS r
Note :
m VS r
The results of these models are shown on the next four
slides, first as a function of polar and then azimuth angle.
141 November, 2014 (DH)
Model A
HTI Model A AVO HTI Model A AVAZ
0.1 0.1
Azimuth 0 Angle 10
0.09 Azimuth 30 0.09 Angle 20
Azimuth 60 Angle 30
Reflection Coefficient
Reflection Coefficient
0.08 Azimuth 90 0.08 Angle 40
0.07 0.07
0.06 0.06
0.05 0.05
0.04 0.04
0.03 0.03
0.02 0.02
0.01 0.01
0 0
0 5 10 15 20 25 30 35 40 45 0 20 40 60 80 100 120 140 160 180
Incidence Angle Azimuth
The reflection coefficients for The reflection coefficients for Model
Model A (change in g 0.1) as a B (change in d 0.1) as a function
function of polar angle for 0, 30, 60 of polar angle for 0, 30, 60 and 90
and 90 degrees azimuth. degrees azimuth.
142 November, 2014 (DH)
Model B
HTI Model B AVO HTI Model B AVAZ
0.1 0.1
Azimuth 0 Angle 10
0.09 Azimuth 30 0.09 Angle 20
Azimuth 60 Angle 30
Reflection Coefficient
Reflection Coefficient
0.08 Azimuth 90 0.08 Angle 40
0.07 0.07
0.06 0.06
0.05 0.05
0.04 0.04
0.03 0.03
0.02 0.02
0.01 0.01
0 0
0 5 10 15 20 25 30 35 40 45 0 20 40 60 80 100 120 140 160 180
Incidence Angle Azimuth
The reflection coefficients for The reflection coefficients for Model
Model A (change in d 0.1) as a A (change in d 0.1) as a function
function of polar angle for 0, 30, 60 of azimuth angle for 10, 20, 30 and
and 90 degrees azimuth. 40 degrees polar.
143 November, 2014 (DH)
Model C
HTI Model C AVO HTI Model C AVAZ
0.1 0.1
Azimuth 0 Angle 10
0.09 Azimuth 30 0.09 Angle 20
Azimuth 60 Angle 30
Reflection Coefficient
Reflection Coefficient
0.08 Azimuth 90 0.08 Angle 40
0.07 0.07
0.06 0.06
0.05 0.05
0.04 0.04
0.03 0.03
0.02 0.02
0.01 0.01
0 0
0 5 10 15 20 25 30 35 40 45 0 20 40 60 80 100 120 140 160 180
Incidence Angle Azimuth
The reflection coefficients for The reflection coefficients for Model
Model C (change in 0.1) as a C (change in 0.1) as a function
function of polar angle for 0, 30, 60 of azimuth angle for 10, 20, 30 and
and 90 degrees azimuth. 40 degrees polar.
144 November, 2014 (DH)
Model D
HTI Model D AVO HTI Model D AVAZ
0.1 0.1
Azimuth 0 Angle 10
0.09 Azimuth 30 0.09 Angle 20
Azimuth 60 Angle 30
Reflection Coefficient
Reflection Coefficient
0.08 Azimuth 90 0.08 Angle 40
0.07 0.07
0.06 0.06
0.05 0.05
0.04 0.04
0.03 0.03
0.02 0.02
0.01 0.01
0 0
0 5 10 15 20 25 30 35 40 45 0 20 40 60 80 100 120 140 160 180
Incidence Angle Azimuth
The reflection coefficients for The reflection coefficients for Model
Model D (change in g 0.15, D (change in g 0.15, 0.05, and
0.05, and d 0.05) as a function of d 0.05) as a function of azimuth
polar angle for 0, 30, 60 and 90 angle for 10, 20, 30 and 40 degrees
degrees azimuth. polar.
145 November, 2014 (DH)
Model D versus polar and azimuth angle
Note that we can think HTI Model D Reflection Coef.
0.11
45
of the previous sets of
curves as “cuts” from a 40 0.1
three-dimensional 35 0.09
Incidence Angle
surface as a function of 30 0.08
polar and azimuth
25
angle. 0.07
20
0.06
The reflection 15
0.05
coefficient surface for
Model D (change in g
10
0.04
0.15, 0.05, and d 5
0.03
0.05) is shown here. 50 100 150
Azimuth w.r.t. isotropy plane
146 November, 2014 (DH)
Arbitrary Fracture Azimuth
In order to let the fracture have any orientation fiso with respect to
North, perform the following rotation N
f fSR fiso
fiso
fSR
f f
R(q , f ) Aiso ( Biso Bani sin 2 f ) sin 2 q
(Ciso d V cos 2 f V sin 2 f sin 2 f ) sin 2 q tan2 q ,
1
to 2
R(q , f ) Aiso ( Biso Bani sin 2 fSR fiso ) sin 2 q
(Ciso
1
2
d V cos 2 fSR fiso V sin 2 fSR fiso sin 2 fSR fiso ) sin 2 q tan2 q ,
147 November, 2014 (DH)
Rüger equation limitations
The Rüger equation was derived for the min
ISO
case of an isotropic half-space over and
HTI half-space. It can be extended to two
Rpp(θ,φ)
HTI half-spaces provided their symmetry
axes are aligned.
HTI(+45°) max
Rüger Zoeppritz Vavryčuk and Pšenčik
-90° -90° -90°
-45° -45° -45°
Azimuth
0° 0° 0°
45° 45° 45°
90° 90° 90°
0° 45° 0° 45° 0° 45°
Angle of incidence Angle of incidence Angle of incidence
148 November, 2014 (DH)
Rüger equation limitations
In the case the symmetry axes are not min
HTI(0°)
aligned the Rüger equation models the
reflectivity incorrectly. A more general
Rpp(θ,φ)
linearization such as Vavryčuk and Pšenčik
(1998) needs to be used in this case.
HTI(+45°) max
Rüger Zoeppritz Vavryčuk and Pšenčik
-90° -90° -90°
-45° -45° -45°
Azimuth
0° 0° 0°
45° 45° 45°
90° 90° 90°
0° 45° 0° 45° 0° 45°
Angle of incidence Angle of incidence Angle of incidence
149 November, 2014 (DH)
Azimuthal AVO (AVAz)
As discussed in a previous section, azimuthal velocity
analysis is also used to correct for the azimuthally-
distorted NMO effects in the COCA gathers.
The amplitudes of these corrected gathers can then be
used to perform azimuthal AVO analysis, or AVAz, using
Rüger’s equation.
The NMO corrected gathers from the previous theory
section are shown on the next slide and the AVAz effects
are as indicated.
150 November, 2014 (DH)
Azimuthal AVO (AVAz)
Az AVO
151 November, 2014 (DH)
November, 2013 After Gray et al. (2004)
Inverting the Rüger equation
To extract the AVAz terms for HTI media, we drop the third term of Rüger equation:
R(qi , f j ) Aiso [ Biso Bani sin 2 (f j fiso )] sin 2 qi ,
where :
R(qi , f j ) the reflection coeff. at incident angle qi
and azimuth angle f j , where i 1,, N , and j 1,, M ,
The Rüger equation is the extension of the standard AVO analysis equation
from isotropic to anisotropic media and solves for four values at each time
sample:
A : the standard intercept,
Biso: the isotropic gradient,
Bani : the anisotropic gradient, a measure of
fracture density.
fiso : the direction of the isotropy plane, which is the same as the
fracture strike.
The equation is non-linear but can be reparameterized in a linear form, which
is discussed next.
152 November, 2014 (DH)
Linearizing the Rüger equation
As shown in Appendix 7, we can re-write the Rüger equation in
linearized form as:
R(qi , f j ) A [ B C cos 2f j D sin 2f j ] sin 2 qi , where :
Bani B cos(2fiso ) B sin(2fiso )
A Aiso , B Biso , C ani , and D ani .
2 2 2
This allows us to solve for the coefficients A, B, C and D using matrix
inversion and solve for A, B, Bani and fiso as follows:
tan 1 ( D / C )
Aiso A, fiso (since : tan 2fiso D / C ),
2
B
Bani 2 C 2 D 2 and Biso B ani B C 2 D 2 ,
2
cos 2 (2fiso ) sin 2 (2fiso )
2
1 Bani
since C 2 D 2 2
Bani .
153
2 November, 2014 (DH)
2
Interpreting the Anisotropic Gradient
0.1
Assuming penny-shaped cracks (Hudson, 1981) gas
0.09
the anisotropic gradient Bani is proportional to hudson wet
Gassmann wet
0.08
the crack density.
0.07
With this interpretation Bani should be positive. 0.06
Bani
0.05
However, Bani is an interface property
0.04
• Upon entering a fractured media Bani should 0.03
be positive 0.02
• Upon exiting a fractured media Bani should 0.01
0
be negative 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
crack density
This interpretation of Bani is relatively simplistic,
a more nuanced discussion occurs in the rock Penny shaped crack theory
physics Appendix.
154 November, 2014 (DH)
Interpreting the Azimuthal Gradient
As just seen, gradient changes
as function of azimuth
If there is no anisotropy the
gradient does not change as a
function of azimuth (e.g. circle
with radius Biso)
In the isotropy plane the AVO
has a gradient of Biso
In the symmetry plane the AVO
has a gradient of Biso +Bani
There is also a sign ambiguity
implicit in the Rüger method,
which will be discussed later.
155 November, 2014 (DH)
Fracture orientation ambiguity
The least squares solution solves for the magnitude
of Bani. Because we don’t know the sign of Bani
there exists two equally valid solutions. Azimuthal reflectivity for a
The figure to the right shows the two isotropic constant angle of incidence
solutions in red and the other in black. 0.04
Both solutions represent the original background 0.03 +Bani
isotropic solution which is then perturbed by the -Bani
0.02
anisotropic gradient.
0.01 fsym fsym
The positive Bani perturbs the red background
isotropic reflectivity to the blue ellipse. In this case 0
the symmetry axis azimuth is defined by red line
-0.01
segment.
The negative Bani perturbs the black background -0.02
isotropic reflectivity to the blue ellipse. In this case -0.03
the symmetry axis azimuth is defined by black line
segment. -0.04
-0.04 -0.03 -0.02 -0.01 0 0.01 0.02 0.03 0.04
Note there is a 90 degree ambiguity between the two
isotropy plane orientations, both potentially
valid.
156 November, 2014 (DH)
Constrained AVAz inversion
In least squares inversion the optimal solution is determined by
minimizing the squared misfit between the model and the data.
Sometimes the problem is ill-conditioned (poor S/N for the amount of
data) or undetermined (not enough data to find a unique solution) so
that the inversion becomes unstable or unsolvable. In these cases we
can use constraints to find an optimal solution which is consistent with
our a priori knowledge.
Constraints have been implemented for the near-offset Rüger
equations. Constraints should help in areas with low fold and poor
azimuthal distributions such as the edge of the survey.
The constraints are calculated from well log statistics using the process
Calculate AVAz Constraints.
The process Azimuthal AVO Volume then optionally reads these
constraints and uses them to run constrained inversion.
Appendix 9 describes the mechanics of calculating and applying
constraints.
The next slide shows an example of the application of constraints.
157 November, 2014 (DH)
Anisotropic gradient: Impact of constraints and weighting on solution
Original With With constraints
constraints and IRLS
and SLS
Constraints can be used to help stabilize the solution. In the above example the
original (unconstrained) inversion of the anisotropic gradient becomes unstable
158
at the left edge of the 3D. Applying constraints and weights improve the solution
November, 2014 (DH)
Rüger equation: potential issues
It is a near offset approximation (ignoring the 3rd term) – potentially introducing
bias.
Assumes that azimuth cannot change as a function of layer.
Bani can be difficult to interpret.
Actually a weighted difference between the tangential and normal fracture
weakness parameters
Bani g ( d T 1 2 g d N )
Ambiguity of what this difference means
Certain Vp/Vs ratios for which two parameters annihilate each other
Only can estimate the magnitude of Bani introducing an ambiguity into the
interpretation of the azimuth estimate.
A new technique which gets around many of these issues, azimuthal Fourier
coefficients (Downton, 2011) is discussed later.
159 November, 2014 (DH)
Exercise 4 – Azimuthal analysis
Now we apply Azimuthal analysis using the Rüger equation.
As discussed in the theory section, the Rüger equation is the extension
of the standard AVO analysis equation from isotropic to anisotropic
media.
It solves for four values at every time sample of every CDP:
A : the standard intercept,
Biso: the isotropic gradient,
Bani : the anisotropic gradient, a measure of
fracture density.
fiso : the direction of the isotropy plane, which should be normal to
the fractures.
160 November, 2014 (DH)
Exercise 4 – Azimuthal AVO Analysis
The workflow is to first interactively to
perform the Rüger inversion in the Azimuthal
AVO analysis tool optimizing the parameters.
Then, once one is satisfied with the results,
perform the inversion on the whole dataset
using the Azimuthal AVO Volume.
To start the azimuthal AVO analysis, open
the Project Manager by clicking on it, go to
the Processes tab and move down the list
until you find the Azimuthal AVO Analysis
option under Azimuthal AVO (ProAZ), as
shown below on the right. Double-click that
option.
Then collapse the Project Manger window to
create more space.
161 November, 2014 (DH)
Exercise 4 – Azimuthal AVO Analysis
In the Azimuthal AVO Analysis dialog,
change the input to trim_statics, Choose
the inline and xline by selecting the well
log icon and selecting well-A.
The Type of Analysis set to Near Offset
Ruger, Perform AVO on Angles from 0 to
30 degrees, and calculate on Angles from
0 to 30, in 2 degrees steps as shown on
the right. The last set of parameters
control the number of AVAz curves
displayed in the AVAz tab.
Make sure the Snail Display is unchecked.
Once finished, click OK.
162 November, 2014 (DH)
Exercise 4 – Azimuthal AVO Analysis
This brings up the
AVO Analysis tool.
Drag the seismic
window so it is
centered around the
base of the fractured
zone at 2244 and
click on the peak.
The AVAz Curves
display show the
Amplitude versus
Azimuth. The data is
color coded by angle
of incidence. The
curve shows the
predicted data for 8
angles of incidence
from 0 to 40 degrees.
The AVAz is similar to
that observed earlier
in the modeling.
163 November, 2014 (DH)
Exercise 4 – Azimuthal AVO Analysis
It is possible to increase the size of graph by dragging the
164 edge to the left November, 2014 (DH)
Exercise 4 – Azimuthal AVO Analysis
Before doing any analysis it is good to look at the input geometry of the
data being viewed. Click on the Data Density Tab and then zoom in
around the gather.
165 November, 2014 (DH)
Exercise 4 – AVAZ Analysis (Data Density Display)
The Data Density display shows the
interpreter the input data density
providing information about how well
the data is sampled in azimuth and
offset.
This dataset was migrated using
azimuthally sectored data hence has
good sampling.
One can see that the maximum
offset in the X-direction is smaller
than the Y-direction. This is an
artifact of rectangular patch that the
data was acquired with.
The diagonal azimuths contain the
largest offsets. These offsets are
not useable because of the poor
azimuthal sampling.
Note the color bar and legend were
modified by right clicking the cursor
on the map.
166 November, 2014 (DH)
Exercise 4 – Azimuthal AVO Analysis
By clicking on the AVO tab one can observe
the Amplitude versus Offset view.
The average AVO based on the near offset
Rüger equation is shown as a solid line. The
best fit A, B parameters are shown at the top of
the screen.
There is quite a bit of near offset noise in this
dataset. Only angles beyond 10 degrees look
believable.
The data is color coded by azimuth using the
phase spectrum color bar. Note the systematic
shift in color about the AVO trend line indicating
the amplitude is also varying as a function of
azimuth in a predictable manner.
By default the data is displayed as a function of
angle of incidence. The user can optionally
toggle between offset and 𝑠𝑖𝑛2 𝜃. Only data less
than 15000 feet have good azimuthal
coverage.
167 November, 2014 (DH)
Exercise 4 – Azimuthal AVO Analysis
By clicking on the AVAz tab one can once
again observe the Amplitude versus
Azimuth display.
The data is color coded by angle of
incidence.
This data has been modeled using the
near-offset Rüger equation with the
estimated parameters shown in the yellow
box at the top of the display. The curves
are modeled over the angle range
specified by the input parameters; from 0
to 30 degrees every 2 degrees. These
ranges are shown under the Parameters
Tab.
168 November, 2014 (DH)
Exercise 4 – Azimuthal AVO Analysis
The next analysis view is the
Domain view. Domain is short for
time or depth Domain. Click on the
Domain tab to see this view.
In this view the amplitude is plotted
in color as a function of azimuth and
angle (/offset). The color of the
data points displays the amplitude at
a particular data location. The
background color shows the model
fitted to this data.
This window can be quite helpful to
see actual sampling of the data. It is
possible to see the data is poorly
sampled at near angles. In addition,
one can see beyond 40 degrees the
azimuthal distribution is no longer
sufficient to perform an AVAz
analysis.
169 November, 2014 (DH)
Exercise 4 – Azimuthal AVO Analysis
Perhaps more intuitively this display
can be shown as a Polar plot
170 November, 2014 (DH)
Exercise 4 – Azimuthal AVO Analysis
The Summary displays the
amplitude as a function of
offset and azimuth. It is
possible to cut and paste
these values into another
program such as Excel.
Note that the source receiver
azimuth is unwrapped so that
it displays between
-180 and 180 degrees
-90 and 90 degrees assuming
reciprocity
171 November, 2014 (DH)
Exercise 4 – Azimuthal AVO Analysis
At any time the
inversion parameters
can be optimized by
going to the
parameters tab.
It is possible to test
the influence of the
number of traces
going into the
supergather, the
maximum offset, the
angle range, the
number of terms
used to fit the AVO
trend, the AVAz
analysis type, etc.
172 November, 2014 (DH)
Exercise 4 – Azimuthal AVO Analysis
It is also possible to separate the graph from the main window by dragging its
title bar to another place on the screen. To attach it again drag the title bar
back to its original location.
173 November, 2014 (DH)
Exercise 4 – Azimuthal AVO Analysis
Using this facility it is possible to display all the views simultaneously
174 November, 2014 (DH)
Exercise 4 – Azimuthal analysis
Once you are satisfied with the
Azimuthal AVO analysis results,
Remove “Azimuthal AVO
Analysis” Scene and Delete,
and then run the inversion on
the whole dataset by running
the Azimuthal AVO volume.
Open up the Project Manager
window by clicking on Project
Manger.
Then go to the Processes tab
and move down the list until you
find the Azimuthal AVO
Volume option under
Azimuthal AVO (AVAz), as
shown below on the right.
Double-click that option.
175 November, 2014 (DH)
Exercise 4 – Azimuthal analysis
Choose the output name Ruger for the output.
The different attributes are output with the following
suffixes
_Biso: Isotropic gradient
_Bani: Anisotropic gradient
_Az: Isotropy plane azimuth
Note the intercept has no suffix.
Change the Trace Range to single inline 190.
Normally, one would run the full 3D but due to
space and time constraints we will run only one
inline. The process is run with 3X3 rolling window
to increase the fold.
Use the parameters as shown to the right. Near
Offset Ruger using angles from 0 to 30 degrees.
Constraints are sometimes helpfult at the edge of
surveys. More information about constraints
appears in Appendix 9.
Once finished, click OK.
176 November, 2014 (DH)
Exercise 4 – Azimuthal analysis
When the
processing is
finished, the
Azimuthal AVO
result is displayed
in View 2, as
shown here.
The default
display shows the
Intercept as the
both the color and
wiggle display.
You will have to
switch to a
horizontal display
And center on the
well.
177 November, 2014 (DH)
Exercise 4 – Azimuthal analysis
Let’s now display the Anisotropic Gradient
1) Open up the
Project
Manager
window, go to
the Project
Data and the
seismic Tab
and drag
Ruger_Bani to
view 1.
November, 2014 (DH)
178
Exercise 4 – Azimuthal analysis
The resultant display
should look similar
to the right. You will
need to optimize the
range. Right click
on the Bani color bar
and adjust the
range. Choose a
maximum value of
0.35 and click OK.
179 November, 2014 (DH)
Exercise 4 – Azimuthal analysis
The resultant display
shows an anomaly
at the top and Base
of the Haynesville.
The Horizon marks
the base of
Haynesville while
the top is 20 ms
above.
Note the anomaly on
the anisotropic
gradient at the both
top and base of the
Haynesville. The
bottom anomaly is
clear if the horizon is
turned off.
180 November, 2014 (DH)
Exercise 4 – Azimuthal analysis
To display the
azimuth in view 2,
right click on the
stack in view2 and
select Ruger_AZ
as the Color Data
Volume
181 November, 2014 (DH)
Exercise 4 – Azimuthal analysis
Resulting in the
Anisotropic
gradient and
isotropy plane
azimuth being
shown on the top
and bottom.
The azimuth only
makes sense if
the media is
anisotropic. If the
media is isotropic
then the azimuth
is meaningless
and should be
excluded.
182 November, 2014 (DH)
Azimuthal Fourier coefficients
It is possible to just measure the amplitude versus azimuth and use
this attribute to predict fractures.
One way to do this is to decompose the AVAz at a particular angle
of incidence into a series of sine waves characterized by their
magnitude and phase.
2nd FC
DC 4th FC
AVAz analysis centered on 35 degrees
183 November, 2014 (DH)
Azimuthal Fourier coefficients
Ikelle (1996) showed that it is possible to write the azimuthal
reflectivity calculated from the anisotropic Zoeppritz equation as a
Fourier series.
Azimuthal Fourier coefficients (FCs) are descriptive – measure of
AVAz.
Azimuthal FCs allow for the separation of the amplitude versus
azimuth (AVAZ) from the amplitude versus offset (AVO) analysis.
Useful information about fracture behaviour can be determined from
just Amplitude versus Azimuth information.
Assuming a single set of aligned vertical fractures, the azimuthal
FCs may be nonlinearly inverted for more fundamental fracture
parameters and an unambiguous estimate of the isotropy plane
azimuth.
184 November, 2014 (DH)
Calculating azimuthal FCs
The workflow for calculating the Fourier coefficients depends on
geometry of the data. The simplest case is for regularly sampled
data in azimuth, e.g. output of an azimuthally sectored migration.
In this case:
– For each azimuth sector, transform the seismic data from offset to
angle of incidence
– Create angle stacks at the angles you want to perform the
analysis on, e.g. 5, 15, 25, 35 degrees.
– For each angle of incidence calculate the Fourier coefficients
using a Fourier Transform.
185 November, 2014 (DH)
Interpretation of Azimuthal FCs
In the case of general anisotropy the linearized Zoeppritz
equations (Pšenčik, I. and J. L. Martins, 2001) can be written
as the azimuthal Fourier series:
R pp (f , q ) r0 r2 cos(2(f f2 )) r4 cos(4(f f4 ))
Where
r0: is the D.C. term
r2: is the magnitude of the n=2 cosine wave
f2: is the phase of the n=2 cosine wave
r4: is the magnitude of the n=4 cosine wave
f4: is the phase of the n=4 cosine wave
186 November, 2014 (DH)
Interpretation of Azimuthal FCs
In the case of HTI media f2 = f4 = fsym, the Rüger equation (2002)
can be written in terms of azimuthal Fourier coefficients:
Rpp (f ,q ) r0 r2 cos(2(f fsym )) r4 cos(4(f fsym ))
where r0 A0 B0 sin 2 q C0 sin 2 q tan 2 q
1 1
r2 Bani sin 2 q sin 2 q tan 2 q
2 4
1
r4 h sin 2 q tan 2 q
16
The D.C. term is the classic 3-term AVO expression where A0, B0, and C0 are
modified to account for the anisotropy.
Ignoring the far offset term the 2nd FC is proportional to the anisotropic gradient Bani.
The 4th FC is proportional to h d, and is a measure of anellipticity.
187 November, 2014 (DH)
Interpretation of Azimuthal FCs
For the special case of a single set of aligned rotationally, invariant
vertical fractures, Downton (2011) shows that:
Rpp (f ,q ) r0 r2 cos(2(f fsym )) r4 cos(4(f fsym ))
where r0 A0 B0 sin 2 q C0 sin 2 q tan 2 q
r2
1
Bani sin 2 q
1 g B h sin 2 q tan 2 q )
41 3g
ani
2
1
r4 h sin 2 q tan 2 q
16
In this special case, the number of parameters and unknowns are reduced by one.
Both Bani and h can be written as functions of the linear slip parameters as shown
on the next slide.
188 November, 2014 (DH)
What are Bani and h?
Bani and h are both weighted differences of the normal and tangential fracture
weakness parameters dN and dT
Bani g (dT d N )
where h 2 g ( d T gd N )
– dN normal fracture weakness:
– function of fluid, crack density and aperture
– fractional parameter d N [0,1)
– dT tangential weakness:
– function of crack density
– fractional parameter dT [0,1)
and
g (Vs /V p ) 2 , and 1 2 g
Note
1) the parameters are defined as the change in layer properties.
2) the fractures for both layers must be aligned.
189 November, 2014
2013 (DH)
2nd azimuthal Fourier coefficient
For small incidence angles the magnitude
of the 2nd Fourier coefficient is R2 (35)
1
r2 (q ) Bani sin 2 q
2 high
Magnitude of the 2nd FC is a scaled version
of the anisotropic gradient
2
Bani r2 (q )
sin q
2
Rüger Bani
• Approximate relationship 0 – 40 degrees
• Similar to Rüger equation, biased estimate
of anisotropic gradient.
• Uncertainty grows as the inverse of sin q .
2
Symmetry axis azimuth wraps every 90 low
degrees
November, 2014 (DH)
190
4th azimuthal Fourier coefficient
The 4th Fourier coefficient is
R2 (35)
1
r4 h sin 2 q tan 2 q
16 high
Exact relationship
Magnitude of the 4th FC is a scaled
version of k
1
h 16 r4
sin q tan q
2 2
R4 (35)
Uncertainty grows as the inverse of
sin 2 q tan 2 q .
low
For n=4 phase (azimuth) wraps every
45 degrees
191 November, 2014 (DH)
2nd FC (Nordegg example) compared to Rüger Bani
Well A Rüger Bani
25 – 35 degrees 30 – 40 degrees 0 – 40 degrees
high
Well B R=0.598 R=0.725 R=0.678
low
192 Correlation coefficient is between attribute and all wells November, 2014 (DH)
Azimuthal inversion
This last section on azimuthal amplitude techniques discusses Azimuthal Inversion,
in particular an approach published by Downton & Roure (2010).
Previous AVAz techniques discussed solve for interface parameters rather than
layer parameters, the Azimuthal Inversion solves for layer properties.
• Easier to interpret layer properties.
Previous methods are band-limited:
• Loss of low frequency.
• Side lobes.
The azimuthal inversion has certain theoretical advantages:
• Includes wavelet in the formulation of the problem.
• Allows the symmetry axis to change as a function of layer.
• No symmetry axis ambiguity.
• Coupling the layers improves the stability.
193 November, 2014 (DH)
Azimuthal elastic inversion is an extension of simultaneous
isotropic inversion (Coulon et al., 2006)
Angle stacks Initial model Final model
q1 q2 q3
TWT, Vp, Vs, ρ
Inversion minimizes a three term cost function using simulated annealing:
194
Data misfit Lateral continuity Prior model
November, 2014 (DH)
Azimuthal elastic Inversion (Downton & Roure, 2010)
f3
q2,3
q1,3
q1,1 q2,1 q3,1
f1
q3,2
TWT, Ip, Is, ρ dN : normal weakness
q2,2
q1,2 dN, dT, Φ dT : tangential weakness
Φ: symmetry axis
f2
Inversion still minimizes a three term cost function:
• Data misfit using the forward modeling explained previously
• Lateral continuity + vertical continuity on the symmetry axis
• Prior model
195 November, 2014 (DH)
Reflectivity
For each layer construct stiffness matrix
based on fractured media
parameterization
Calculate azimuthal reflectivity response
for each interface
Reflectivity calculated using
– Anisotropic Zoeppritz equations
– Schoenberg and Protázio (1992)
– Linearized Zoeppritz equations
– Vavryčuk and Pšenčik (1998)
– Pšenčik and Martins (2001)
The symmetry axis changing as a
function of layer causes a rotation of the
azimuthal reflectivity as a function of
incidence angle.
196 November, 2014 (DH)
P-impedance TWT (ms) Tangential weakness
1000
1100
1200
1300
1400 low high
197 November, 2014 (DH)
Azimuthal AVO: map of tangential weakness (hot=high) to confirm areas of
elevated stress. Areas of increased stress anisotropy are focused at the North
half pad 1 along the strike-slip fault.
(Goodway et al., The Leading Edge, December 2012)
198 November, 2014 (DH)
Summary
Near offset Azimuthal Az. FC Azimuthal
Rüger FCs inversion Inversion
Parameterization Model based/ Descriptive Model based Model based
“Descriptive”
# of fracture parameters 2 4 3 3
Theoretical model and Single set of Descriptive Single set of Single set of
assumptions (+/-) aligned vertical aligned vertical
fractures vertical fractures,
fractures azimuth may
vary
Parameter type Band-limited Band-limited Band-limited layer
interface interface interface
Fracture azimuth ambiguity Yes Yes No No
ease to interpret
Robustness
Maturity
Ease of calculation More Simple More Most
199
Involved Involved
November, 2014 (DH)
Involved
Exercise 6 – Visualizing Vector data
This last exercise demonstrates the use of the interactive
Rose Diagram Analysis and Visualization to analyze and
display the magnitude and azimuth information together.
Due to space and time constraints we only ran one inline in
the previous exercise. This exercise starts by reading in a
project where the Rüger inversion was run on the whole
dataset.
After reading the project, we then create maps at the base on
the Haynesville on the anisotropic gradient and azimuth.
These maps and volumes are then used by
– the interactive Rose Diagram to display the azimuth distribution for a
particular zone on a map, for example centered around a well location.
– the visualization software to display vector data. The results of the azimuthal
inversion contain both magnitude (anisotropic gradient) and direction
(Azimuth) information.
200 November, 2014 (DH)
Exercise 6 – Visualizing Vector data
Start up the Hampson-
Russell suite again.
Now click the Projects
tab and choose the
option to Find Project.
201 November, 2014 (DH)
Exercise 6 – Visualizing Vector data
A dialog appears, where we set
the project name. Here we select
the Visualizion.prj, as shown
here.
Enter the project name and click
OK on that dialog:
202 November, 2014 (DH)
Exercise 6 – Visualizing Vector data
This project contains the results of inverting the full volume using the Ruger
inversion. The anisotropic gradient is called Full_Bani and the isotropic plane
azimuth is called Full_az. First, let’s extract the amplitude of the anisotropic
gradient at the base of the Haynesville. Go to the Processes tab and select
Create Data Slice.
Selcted as Input
Full_Bani, then select
the Event HorizonA
this corresponds to the
base of Haynesville
and the peak we picked
previously. We will use
the default of a window
centered on the horizon
using a fixed window
size of 10 ms to extract
the amplitude. Name
the output BaniSurf
and then click OK.
203 November, 2014 (DH)
Exercise 6 – Visualizing Vector data
This creates a map of
the azimuth calculated
over a 10 ms interval at
the Base of the
Haynesville.
204 November, 2014 (DH)
Exercise 6 – Visualizing Vector data
Rerun the Create Data Slice
Process again, but this time
specifying as Input: Full_Az.
Specify as output AzSurf,
and click OK.
205 November, 2014 (DH)
Exercise 6 – Visualizing Vector data
This creates a map of
the azimuth calculated
over a 10 ms interval at
the Base of the
Haynesville.
206 November, 2014 (DH)
Exercise 6 – Rose Diagram Analysis
Initiate the Rose Diagram Analysis under the
Processes Tab.
In order to interactively determine the
azimuth distribution for a particular zone on
the map we use the Data Slice option.
207 November, 2014 (DH)
Exercise 6 – Rose Diagram Analysis
The primary attribute is azimuth data,
choose AzSurf.
Select BaniSurf as the secondary
attribute. The secondary attribute is
used both for display and to filter the
primary attribute.
Azimuths associated with small Bani
values are not meaningful, they are
often related to noise. These values
can be excluded by specifying a
minimum threshold.
It is typically easier to pick the zones to
analyze on the anisotropic gradient
attribute. Either attribute can be
displayed by selecting the Bring to
Front under the appropriate attribute.
In this case Bring to Front the BaniSurf.
Lastly, click on Scan to initiate the
analysis.
208 November, 2014 (DH)
Exercise 6 – Rose Diagram Analysis
This brings up the Rose
Diagram using all the
azimuth values.
To filter the values
associated with the small
anisotropic values select
the Histogram Parameters
tab.
209 November, 2014 (DH)
Exercise 6 – Rose Diagram Analysis
The minimum threshold based
on the anisotropic gradient is
selected by clicking on the
From text box.
This brings up a slider. Using
the slider specify a minimum
value of 39000 to exclude the
azimuths associated with small
anisotropic gradient values and
then click Scan.
The Azimuth distribution is now
more uniform.
Try out a number of thresholds
to see its influence.
210 November, 2014 (DH)
Exercise 6 – Rose Diagram Analysis
To interactively determine the azimuths
at a particular location on a map, select
the Zones on Map tap and the option
Draw temporary zones on the map.
The select the Interactive
Rose Diagram which brings
up the interactive menu
which allows you to choose
the ZOI using a Rectangle.
Interactively select the area
around well A.
Then click Scan.
211 November, 2014 (DH)
Exercise 6 – Rose Diagram Analysis
You should see a display similar to the one below. Most of the fractures
strike East West similar to the well control and the stress in the area.
212 November, 2014 (DH)
Exercise 6 – Visualization
Another way to view
the anisotropic
gradient with the
azimuth information
is to use the
Visualization
Software. To do this,
go to the View3D tab
on the right part of
the window and
select the Project
Data tab on the left.
Select the Full_Bani
volume and drag it
over to the View3D
tab.
213 November, 2014 (DH)
Exercise 6 – Visualization
Click OK on the
next dialog to
move the whole
volume over.
214 November, 2014 (DH)
Exercise 6 – Visualization
Next, select the
Full_az volume
and drag it over to
the View3D tab,
again clicking OK
on the next dialog.
Now there are two
volumes in the
viewing window,
but only the last is
being displayed in
color, as seen at
the top.
215 November, 2014 (DH)
Exercise 6 – Visualization
Next, open the
Horizon tab
and select
Horizon.
Drag it over to
the View3D
window.
216 November, 2014 (DH)
Exercise 6 – Visualization
Then, select well
well-A from the
Well data. Drag
it over to the
View3D window.
217 November, 2014 (DH)
217
Exercise 6 – Visualization
Finally, select
BaniSurf from
the slice data.
Drag it over to
the View3D
window.
218 November, 2014 (DH)
218
Exercise 6 – Visualization
Then, make more space by closing the two side bars:
219 November, 2014 (DH)
Exercise 6 – Visualization
Now, turn on
more viewing
planes. This
allows you to
slide through the
volume a slice at
a time in the X, Y
and Z directions.
You can also
switch between
viewing Azimuth
and Anisotropic
Gradient using
the Color
Volume button
above the view.
220 November, 2014 (DH)
Exercise 6 – Visualization
However, there is another
option for visualization, in
which we can view both the
Azimuth and Anisotropic
Gradient simultaneously in
a vector-type display.
To use this option, right-
click in the black space and
select Data Config from
the Vector option.
221 November, 2014 (DH)
Exercise 6 – Visualization
In the data configuration option,
select Azimuth in the upper box
and Anisotropic_Gradient in the
lower box. Click OK on the first
dialog to bring up the second
dialog.
Set the Min. Threshold to 0 and
the Max threshold to 150000,
make sure the Size Ratio to Min
Threshold is unchecked, click on
Horizon Based to attach glyphs
to surface, then click then OK.
Note these parameters can be
updated later.
222 November, 2014 (DH)
Exercise 6 – Visualization
Zoom a few times.
Note that as you slide
the horizontal slice
through the volume,
“glyphs” appear,
whose direction
indicates the azimuth
and color represents
the gradient.
You can use the arrow
keys move and center
the volume.
223 November, 2014 (DH)
Exercise 6 – Visualization
Turn off x, y and z so only
the surface is left
Right click screen and
select Surface View
Parameters
224 November, 2014 (DH)
Exercise 6 – Visualization
Then, click on
the BaniSurf,
modify the range
from 0 to
100,000 and use
the AVO
Envelope color
scheme.
The background
color is set using
Domain Data.
In this case we
set it to TWT.
Click on x to exit
the View
Parameters
dialog
225 November, 2014 (DH)
Exercise 6 – Visualization
At this point we can optimize the This sets the minimum
Vector Parameter Display gradient value for plotting:
Right click in the black space:
This controls the
size of the glyphs:
This plots the glyphs
on the horizon:
226 November, 2014 (DH)
Exercise 6 – Visualization
227 November, 2014 (DH)
Exercise 6 – Visualization
We have now finished the
Azimuthal Analysis exercises.
Close down the Geoview
program by clicking File>Exit
on the upper left corner of the
Geoview window:
End of Exercise 6
228 November, 2014 (DH)
Course Summary
In this course we have seen how
– Fractures and differential horizontal stress introduce velocity anisotropy into
the earth.
– This anisotropy can be observed in P-wave seismic data as azimuthal
amplitude and traveltime variations.
– Have used these azimuthal variations to predict fractures & stress field.
We have examined three different ways to quantify the
azimuthal amplitude variations,
– including the near-offset Rüger equation, Azimuthal Fourier coefficients, and
Azimuthal inversion
Also included in the course notes are an extensive set of
references and Appendices on advanced topics.
229
November, 2014 (DH)
References
230 November, 2014 (DH)
230
References
Aki, K. and Richards, P.G., 1980, Quantitative Seismology: Theory and Methods,
W.H.Freeman and Co., San Francisco, 932 pp.
Angus, D.A., Verdon, J.P., Fisher, Q.J., Kendall, J.M., 2009, "Exploring trends in
microcrack properties of sedimentary rocks: An audit of dry-core velocity-stress
measurements", Geophysics, vol. 74, p. E193-E203.
Al-Dossary, S., and K. J. Marfurt, 2006, Multispectral estimates of reflector curvature and
rotation: Geophysics, 71, P41-P51.
Bakulin, A., V. Grechka, and I. Tsvankin, 2000, Estimation of fracture parameters from
reflection seismic data—Part I: HTI model due to a single fracture set, Geophysics,
65, 1788
Banik, N.C., 1987, An effective anisotropy parameter in transversely isotropic media:
Geophysics, 52, p. 1564.
Bergbauer, .S., T. Mukerji, and P. Hennings, 2003, Improving curvature analyses
of deformed horizons using scale dependent filtering techniques; AAPG bulletin, v.
87, no. 8, 1255-1272
Blangy, J.P., 1994, AVO in transversely isotropic media-An overview: Geophysics, 59, p.
775-781.
231 November, 2014 (DH)
References
Budiansky, B., and R. J. O’Connell, 1976, Elastic moduli of a cracked solid: International
Journal of Solids and Structures, 12, 81–97.
Carcione, J. M., 1992, "Anisotropic Q and velocity dispersion of finely layered media",
Geophysical Prospecting, Vol. 40, p. 761–783.
Cary, P.W., 1999, "Generalized sampling and "beyond Nyquist" imaging", 11th Annual
CREWES Report.
Coulon, J.-P., Lafet, Y., Deschizeaux, B., Doyen, P.M., and Duboz, P., 2006, Stratigraphic
elastic inversion for seismic lithology discrimination in a turbiditic reservoir: 76th
Annual International Meeting, SEG, Expanded Abstracts, 25, no. 1, 2092-2096.
Daley, P.F., and Hron, F., 1977, Reflection and transmission coefficients for transversely
isotropic media: Bull. Seis. Soc. Am., 67, 661-665.
_______, 1979, Reflection and transmission coefficients for seismic waves in
ellipsoidally anisotropic media: Geophysics, 44, p. 27-38.
Delbecq, F., J, Downton, R. Bale, J. Durand, D. Schmidt and B. Roure, 2011, A
comparison of azimuthal seismic techniques: Presented at the CSPG CSEG CWLS
Joint Convention.
232 November, 2014 (DH)
References
Downton, J., and D. Gray, 2006, AVAZ parameter uncertainty estimation; SEG
Expanded Abstracts, 25,234
Downton, J., 2011, Azimuthal Fourier coefficients: A simple method to estimate fracture
parameters, SEG, Expanded Abstracts, 30, 269
Downton, J., D. Holy, D. Trad, L. Hunt, S. Reynolds and S. Hadley, 2010, The effect of
interpolation on imaging and azimuthal AVO: A Nordegg case study, SEG, Expanded
Abstracts, 29, 383
Downton, J. and B. Roure, 2010, Azimuthal simultaneous elastic inversion for fracture
detection: SEG, Expanded Abstracts, 29, 263.
Downton, J., B. Roure, L. Hunt, 2011, Azimuthal Fourier Coefficients: CSEG
recorder, 37, vol. 10, 22-27.
Gassmann, F., 1951, Elastic waves through a packing of spheres: Geophysics, 16,
673-685.
Goodway B., J. Varsek, and C. Abaco , 2006, Practical applications of P-wave AVO for
unconventional gas Resource Plays, Part II: Detection for fracture prone zones with
Azimuthal AVO and coherence discontinuity: RECORDER, 31,4, 52-65
Gray, F. D., and Todorovic-Marinic, D., 2004, Fracture detection using 3D
azimuthal AVO: CSEG Recorder, 29, 5–8.
233 November, 2014 (DH)
References
Grechka V., I Tsvankin, and J.K. Cohen, 1999, Generalized Dix equation and
analytic treatment of normal-moveout velocity for anisotropic media:
Geophysical Prospecting, 47, 117-148
Grechka V., and M. Kachanov, 2006, Effective elasticity of fractured rocks: A
snapshot of the work in progress; Geophysics, 71, W45-W58
Guéguen, Y. amd V. Palciauskas, 1994, Introduction to the Physics of Rocks, Princeton
University Press, New Jersey, 294 pp.
Gurevich, B., 2003, Elastic properties of saturated porous rocks with aligned fractures:
Journal of Applied Geophysics, 54, 203–218.
Gurevich, B., and M. Pervukhina, 2010, An analytical model for stress-induced
anisotropy of a cracked solid: SEG, Expanded Abstracts
Hudson, J. A., 1981, Wave speeds and attenuation of elastic waves in material
containing cracks: Geophys. J. Royal Astronomy. Soc. 64, 133-150.
Hunt, L., S. Reynolds, S. Brown, S. Hadley, J. Downton, and S. Chopra, 2010,
Quantitative estimate of fracture density variations in the Nordegg with azimuthal
AVO and curvature: A case study; TLE, vol. 29, 9, 1122-1137.
234 November, 2014 (DH)
References
Hunt, L., S. Reynolds, S. Hadley, J. Downton, and S. Chopra, 2011, Causal fracture
prediction: Curvature, stress, and geomechanics: TLE, vol. 30, 11, 1274-1286.
Hsu, J.C., and M. Schoenberg, 1993, Elastic waves through a simulated fractured
medium, Geophysics,58,964
Ikelle L.T. 1996. Amplitude variations with azimuths (AVAZ) inversion based on linearized
inversion of common azimuth sections. In: Seismic Anisotropy (eds E. Fjaer, R. Holt
and J.S. Rathore), pp. 601-644. SEG, Tulsa.
Kachanov, M., 1992, Effective elastic properties of cracked solids, Critical review of
some basic concepts: Applied Mechanics Review, 45, 304–335.
——–, 1993, Elastic solids with many cracks and related problems, in J. W.
Hutchinson and T.Wu, eds., Advances in applied mechanics 30: Academic Press Inc,
259–445.
Lancaster, S. and Whitcombe, D., 2000, "Fast-track 'coloured' inversion", SEG 2000
Expanded Abstracts.
Nur, A. and Simmons, G., 1969, Stress-Induced Velocity Anisotropy in Rock: An
Experimental Study, JGR, 74, no. 27, 6667-6674.
235 November, 2014 (DH)
References
Mavko, G., Mukerji, T., and Dvorkin, J., 1998, The Rock Physics Handbook: Cambridge
University Press, Cambridge, 329 pp
Nur, A. and Simmons, G., 1969, Stress-Induced Velocity Anisotropy in Rock: An
Experimental Study, JGR, 74, no. 27, 6667-6674.
O’Connell, R. J., and B. Budiansky, 1974, Seismic velocities in dry and saturated
cracked solids: Journal of Geophysical Research, 79, 5412–5426.
Pšenčik I. and D. Gajewski, 1998, Polarization, phase velocity, and NMO velocity of P-
waves in arbitrary weakly anisotropic media; Geophysics, 63, 1754-1766
Pšenčik I. and J. L. Martins, 2001, Properties of weak contrast PP
reflection/transmission coefficients for weakly anisotropic elastic media: Studia
Geophysica et Geodaetica, 45, 176-197
Reches Z., 1998, "Tensile Fracturing of Stiff Rock Layers Under Traxial Compressive
Stress States", 3rd MARMS, Int. J. of Rock Mech. Y Min. Sci. vol. 35, no. 4-5, Paper
No. 70.
Rich J.P. & M. Ammerman, 2010, Unconventional Geophysics for Unconventional Plays:
SPE 31779
236 November, 2014 (DH)
References
Roberts, A., 2001, Curvature attributes and their application to 3D interpreted horizons.
First Break, 19, 85-99.
Rüger, A., 1997, P-wave reflection coefficients for transversely isotropic models with
vertical and horizontal axis of symmetry: Geophysics, 62, 713-722.
Rüger, A., and I. Tsvankin, 1997, Using AVO for fracture detection: Analytic basis and
practical solutions: The Leading Edge, 10, 1429– 1434.
Rüger, A., 1998, Variation of P-wave reflectivity with offset and azimuth in anisotropic
media: Geophysics, 63, 935-947.
Rüger, A., 2002, Reflection coefficients and azimuthal AVO Analysis in anisotropic
media, SEG geophysical monograph series number 10: Soc. Expl. Geophys.
Sayers, C., and S. Dean, 2001, Azimuth-dependent AVO in reservoirs containing non-
orthogonal fracture sets: Geophysical Prospecting, 2001, 49, 100-106.
Sayers, C. M., and M. Kachanov, 1991, A simple technique for finding effective elastic
constants of cracked solids for arbitrary crack orientation statistics: International
Journal of Solids and Structures, 27, 671–680.
237 November, 2014 (DH)
References
Schoenberg., M., 1980, Elastic behaviour across linear slip interfaces. Journal of the
Acoustical Society of America, 68:1516–1521
Schoenberg, M., and J. Douma, 1988, Elastic wave propagation in media with parallel
fractures and aligned cracks: Geophysical Prospecting, 36, 571–590.
Schoenberg, M. and J. Protázio, 1992, “Zoeppritz” rationalized and generalized to
anisotropy, J. Seismic Expl. 1, 125-144.
Schoenberg, M., and C. Sayers, 1995, Seismic anisotropy of fractured rock, Geophysics,
60, 204–211.
Shaw, R.K., and M. K. Sen, 2006, Use of AVOA data to estimate fluid indicator in a
vertically fractured medium, Geophysics,71,C15
Sheriff, R.E., 2002, Encyclopedic Dictionary of Applied Geophysics, SEG geophysical
reference series number 13: Soc. Expl. Geophys., p13
Thomsen, L., 1986, Weak elastic anisotropy: Geophysics 51, p 1954.
Thomsen, L., 1993, Weak anisotropic reflections, in Castagna, J.P. and
Backus, M, Ed., Offset Dependent Reflectivity, SEG.
238
238 November, 2014 (DH)
References
Thomsen, L., 1986, Weak elastic anisotropy: Geophysics 51, p 1954.
Thomsen, L., 1993, Weak anisotropic reflections, in Castagna, J.P. and
Backus, M, Ed., Offset Dependent Reflectivity, SEG.
Tsvankin, I, 1997, Anisotropic parameters and P-wave parameters for orthorhombic
media: Geophysics, 62, p. 1292-1309.
Tsvankin, I, 2001, Seismic signatures and analysis of reflection data in anisotropic
media: Elsevier Science Publ. Co. Inc.
Vavryčuk V. and Pšenčik I., 1998: PP-wave reflection coefficients in weakly anisotropic
media. Geophysics, 63, 2129-2141.
Vermeer, G.J.O., 2002, "3-D Seismic Survey Design", SEG Geophysical Reference
series, vol. 12, p. 1-205.
Wiggens, R., Kenny, G.S., and McClure, C.D., 1983, A method for determining and
displaying the shear-velocity reflectivities of a geologic formation: European Patent
Application 0113944.
Verdon J.P. and A. Wüstefeld, 2013, Measurement of the normal/tangential fracture
compliance ratio (ZN/ZT) during hydraulic fracture stimulation using S-wave splitting
data: Geophysical Prospecting: Geophysical Prospecting, 61, 461-477.
239
239 November, 2014 (DH)
Appendix 1: Phase velocity
Appendix 2: Bond Transformations and Transverse isotropy
Appendix 3: Orthorhombic, monoclinic, and triclinic anisotropy
Appendix 4: Schoenberg-Protázio extension to the Zoeppritz equation
Appendix 5: Rock Physics of Fractures
Appendix 6: BN/BT ratio
Appendix 7: The Rüger equation
Appendix 8: Azimuthal Geometric Spreading Correction
Appendix 9: Constraints
Appendix10: filtering Azimuths
240 November, 2014 (DH)
Appendix 1: Phase velocity
241 November, 2014 (DH)
Appendix 1: Phase velocity
The stiffness matrix C shown here is form the most general
anisotropic case, triclinic.
c11 c21 c13 c14 c15 c61
c c22 c23 c24 c25 c62
21
c31 c32 c33 c34 c35 c63
C
c41 c42 c43 c44 c45 c64
c51 c52 c53 c54 c55 c65
c61 c62 c63 c64 c65 c66
In this appendix we show how the coefficients cij allow us to
calculate the seismic phase velocities in the general case.
242 November, 2014 (DH)
Appendix 1: Phase velocity
Using the coefficients in C, we set up the Kelvin-Christoffel
equations are as follows:
11 12 13
12 22 23 , where :
13 23 33
11 n12c11 n22c66 n32c55 2n1n2c16 2n1n3c15 2n2n3c56 ,
22 n12c11 n22c66 n32c55 2n1n2c16 2n1n3c15 2n2n3c56 ,
33 n12c11 n22c66 n32c55 2n1n2c16 2n1n3c15 2n2n3c56 ,
12 n1n2 (c12 c66 ) n12c16 n22c26 n32c45 n1n3 ( c14 c56 ) n2n3 (c46 c25 ),
13 n1n3 (c13 c55 ) n12c15 n22 c46 n32 c35 n1n2 ( c14 c56 ) n2n3 ( c36 c45 ),
23 n2 n3 ( c23 c44 ) n12c56 n22c24 n32c34 n1n2 (c46 c25 ) n1n3 ( c36 c45 ),
and n (n1 , n2 , n3 )T propagation direction.
243 November, 2014 (DH)
Appendix 1: Phase velocity
The propagation
direction vector and its
relationship to both
Cartesian and z
q n = [n1,n2,n3]T
spherical coordinates
f y
is shown here, where
the red arrows show x
the coordinate system
and the blue arrow
represents a particular
propagation direction.
244 November, 2014 (DH)
Appendix 1: Phase velocity
The P and S-wave velocities are computed as follows:
First, a propagation direction is chosen and the Kelvin-Christoffel
(KL) matrix is computed using the stiffness matrix values.
Then, a spectral decomposition is done on the KL matrix:
UU 1 , where :
11 12 13 u1 v1 w1 l1 0 0
12 22 23 , U u2 v2 w2 , 0 l2 0 ,
13 23 33 u3 v3 w3 0 0 l3
u1 v1 w1
and u u2 , v v2 , w w2 , and li , i 1,2,3, are the
u3 v3 w3
eigenvecto rs and eigenvalue s, respectively, of .
245 November, 2014 (DH)
Appendix 1: Phase velocity
For orthorhombic and higher symmetries (i.e. isotropic, VTI,
HTI) the stiffness matrix is much simpler:
c11 c21 c13 0 0 0
c c c 0 0 0
21 22 23
c31 c32 c33 0 0 0
C
0 0 0 c44 0 0
0 0 0 0 c55 0
0 0 0 0 0 c66
This allows us to write the Kelvin-Christoffel equations in
much simpler form, as shown next.
246 November, 2014 (DH)
Appendix 1: Phase velocity
For the orthorhombic (or higher symmetry) case we get:
11 12 13
12 22 23 , where :
13 23 33
11 n12c11 n22c66 n32c55, 22 n12c11 n22c66 n32c55,
33 n12c11 n22c66 n32c55, 12 n1n2 ( c12 c66 ),
13 n1n3 ( c13 c55 ), 23 n2n3 (c23 c44 ),
and n ( n1 , n2 , n3 )T propagation direction.
The figure on the next slide shows a VTI example.
247 November, 2014 (DH)
Appendix 1: Phase velocity
l1 l2 l3
The velocities are as follows: VP ,VS1 ,VS 2
r r r
For example, the figure on the right
shows the velocities for the VTI shale
described by Hornby (1998), with Z
4100
parameters given by: 1
4000
VP0= 3442.4, VS0= 1763.6, r=2.540 0.8
= 0.2309, d = 0.1102, g = 0.3544, 0.6
3900
and: 0.4 3800
44.0 17.0 17.4 0 0 0 0.2 3700
17.0 44.0 17.4 0 0 0 0
1 3600
17.4 17.4 30.1 0 0 0 0.5 1
CVTI GPa 0
0
0.5
3500
0 0 0 7.9 0 0 -0.5 -0.5
0 X
0 7.9 0
-1 -1
0 0
Y
0 0 0 0 0 13.5
248 November, 2014 (DH)
Appendix 2: Bond
Transformations and
Transverse isotropy
249 November, 2014 (DH)
Appendix 2: Bond Transformations
It can be shown that the HTI case is identical to the VTI case except
that we rotate from a horizontal symmetry axis to a vertical
symmetry axis using the general Bond rotation given by:
C ' MCM T , where
a112 2
a11 2
a11 2a12a13 2a12a13 2a12a13
2 2 2
a
11 a 11 a 11 2 a a
12 13 2 a a
12 13 2 a a
12 13
a112 2
a11 2
a11 2a12a13 2a12a13 2a12a13
M ,
a21a31 a21a31 a21a31 a22a33 a23a32 a22a33 a23a32 a22a33 a23a32
a21a31 a21a31 a21a31 a22a33 a23a32 a22a33 a23a32 a22a33 a23a32
a21a31 a21a31 a21a31 a22a33 a23a32 a22a33 a23a32 a22a33 a23a32
and aij cosine of angle between i th and j th axes.
250 November, 2014 (DH)
Appendix 2: Bond Transformations
As shown by Carcione (2001), based on the work of Thomsen
(1988), the rotation matrix for a rotation from a vertical
symmetry axis to a horizontal symmetry axis, followed by a
rotation of q around the new z-axis, is given by:
0 sin 2 q cos 2 q sin 2q 0 0
0 cos 2
q sin 2 q sin 2q 0 0
1 0 0 0 0 0
M
0 0 0 0 sin q cos q
0 0 0 0 cos q sin q
1 1
0 sin 2q sin 2q cos 2q 0 0
2 2
251 November, 2014 (DH)
Appendix 2: Bond Transformations
If q = 0, this gives:
0 0 1 0 0 0
0 1 0 0 0 0
1 0 0 0 0 0
M
0 0 0 0 0 1
0 0 0 0 1 0
0 0 0 1 0 0
Application of this rotation matrix transforms a VTI medium
into an HTI medium, as shown in the next slide for the
Mesaverde mudshale.
252 November, 2014 (DH)
Appendix 2: Bond Transformations
55.21 14.99 24.41 0 0 0
14.99 55.21 24.41 0 0 0
24.41 24.41 51.69 0 0 0
CVTI
0 0 0 18.41 0 0
0 0 0 0 18.41 0
0 0 0 0 0 20.11
51.69 24.41 24.41 0 0 0
24.41 55.21 14.99 0 0 0
24.41 14.99 55.21 0 0 0
MCVTI M CHTI
T
0 0 0 20.11 0 0
0 0 0 0 18.41 0
0 0 0 0 0 18.41
253 November, 2014 (DH)
Appendix 3: Orthorhombic,
monoclinic, and triclinic
anisotropy
254 November, 2014 (DH)
Appendix 3: Orthorhombic media
Orthorhombic media are described by 9 independent parameters
c11 c12 c13 0 0 0
c c c 0 0 0
12 22 23
c13 c23 c33 0 0 0
C ,
0 0 0 c44 0 0
0 0 0 0 c55 0
0 0 0 0 0 c66
These media are characterized by 3 mutually orthogonal planes of symmetry
x2 x2
x1 x1
Orthogonal Identical Vertical fractures embedded in a
fractures fractures horizontally-layered shale
255 November, 2014 (DH)
Appendix 3: Orthorhombic media
Tsvankin (1997) extended Thomsen’s parameters to orthorhombic media.
The parameters are similar to Thomsen’s and are defined in specific symmetry
planes.
Each parameter is referenced by a superscript which refers to the normal of the
symmetry plane in which the parameter is defined.
The 2nd normal corresponds to the symmetry plane used by both the VTI and HTI
definitions
C33 C55
reference velocitie s VP 0 , VS 0
r r
(x, z) plane 2 C C33 2 C66 C44
11 ,g
C C55 C33 C55
, d 2 13
2 2
2C33 2C44 2C33 C33 C55
(y,z) plane 1 C C33 2 C66 C55
22 ,g
C C44 C33 C44
, d 2 23
2 2
2C33 2C55 2C33 C33 C44
(x, y) plane
C
d 3 12
C 66 2
C11 C 66 2
2C11C11 C66
256 November, 2014 (DH)
Appendix 3: Orthorhombic media
Tsvankin (1997) derived an approximate expression
for the phase velocity in an orthorhombic medium in
terms of these parameters:
V p q , f VP 0 1 d f sin 2 q cos 2 q f sin 4 q
where
d f d 1 sin 2 f d 2 cos 2 f
f 1 sin 4 f 2 cos 4 f (2 2 d 3 ) sin 2 f cos 2 f
257 November, 2014 (DH)
Appendix 3: Orthorhombic media
P-wave velocity vs azimuth and dip for Dillen sandstone (m/s)
2900
Z
1
18.6 1.60 2.10 0 0 0 0.8
2880
1.60 20.1 2.60 0 0 0
0.6
2.10 2.60 19.8 0 0
2860
0
COrtho , 0.4
0 0 0 8. 6 0 0 2840
0 0 0 0 8.5 0 0.2
0 0 0 0 0 8.9 0 2820
1
0.5 1
0 0.5
2800
X -0.5
0
-0.5 Y
-1 -1
• Phase velocity of orthorhombic medium (Dillen sandstone).
VP0= 2884.3, VS0= 1889.8
1 = 0.0076, 2 = -0.0303
d1 = 0.0000, d2 = -0.0343, d3 = 0.0448,
258 g1 = 0.0235, g2 = 0.0174 November, 2014 (DH)
Appendix 3: Monoclinic media
There are more complex forms of anisotropy such as
monoclinic media. In this case there is only 1 mirror symmetry
plane. Monoclinic media might arise due to 2 sets of unequal
vertical fractures. Monoclinic media are described by 13
parameters.
c11 c12 c13 0 0 c16
c c22 c23 0 0 c26
12
c13 c23 c33 0 0 c36
Cmono
0 0 0 c44 c45 0
0 0 0 c45 c55 0
2 sets of vertical fractures c16 c26 c36 0 0 c66
259 November, 2014 (DH)
Appendix 3: Triclinic media
Multiple fractures at different orientations and dips give rise to
even more complex forms of anisotropy. For general
anisotropy 21 parameters are required to define the triclinic
stiffness matrix.
c11 c21 c13 c14 c15 c61
c c22 c23 c24 c25 c62
21
c31 c32 c33 c34 c35 c63
C .
c41 c42 c43 c44 c45 c64
c51 c52 c53 c54 c55 c65
c61 c62 c63 c64 c65 c66
260 November, 2014 (DH)
Appendix 3: Bond transformation
If HTI anisotropy is caused by a single set of 2900
vertical fractures, the symmetry axis is 1
2890
2880
perpendicular to the fractures. 0.8
2870
Typically the symmetry axis for an HTI 0.6 2860
stiffness matrix is written as the x-axis. If we 0.4 2850
2840
want it to be some other axis or arbitrary 0.2
2830
azimuth then we can perform a Bond 0
1 2820
transformation on it (see Appendix A). 0.5 1
2810
0 0.5
After performing the bond transformation the
0 2800
-0.5 -0.5
-1 -1
symmetry axis (and fractures) can have an
arbitrary azimuth and dip.
2880
The azimuthally rotated stiffness matrix will 1
have the form of a monoclinic stiffness matrix. 0.8
2870
2860
The figures to the right show the phase
0.6
velocity of the Dillen (2001) stiffness matrix 0.4 2850
before and after the application of a 30 0.2
2840
degree counter-clockwise rotation about the 0
1
z-axis. 0.5 1
2830
0 0.5
0 2820
-0.5 -0.5
261 November, 2014 (DH)
-1 -1
Appendix 4: Schoenberg-
Protázio extension to the
Zoeppritz equation
262 November, 2014 (DH)
Schoenberg-Protázio AVO
Schoenberg and Protázio (1992) extended the Zoeppritz
equations to orthorhombic anisotropy by introducing the
following 2X2 block matrices:
αs1 βs3S 2 ραβ 2 s1s3P ρβ Γ
X , and Y
ρα Γ 2 ρβ s1s3S
3
αs 3P βs1
sin q
where : VP , VS , s1 horizontal slowness,
1 cos q
s3P s
2
vertical P - wave slowness,
2 1
1 cos q
s3S s
2
vertical S - wave slowness,
2 1
and : Γ 1 2 β 2 s12.
November, 2014 (DH)
263
Schoenberg-Protázio AVO
Based on the 2x2 impedance matrices just defined, and the
equivalent matrices for the second layer, written as primed
quantities below, the AVO reflection and transmission
matrices are defined:
R X X Y Y X X Y Y
1 1 1 1 1
T 2 X X Y Y
1 1 1
TPP TPS RPP RPS
where : T , and R .
TSP TSS RSP RSS
264 November, 2014 (DH)
Schoenberg-Protázio AVO
Since the Schoenberg/Protázio matrices are difficult to
visualize, it is instructive to look at the zero offset case:
0 1 0 ρβ
X( 0 ) 0
, and Y( 0 )
0
ρα 0 1 0
ρα 1 0
0 ρβ ,
X -1 X ρα , and Y -1Y
0
ρβ
0 1
ρα ρα 2 ρα
ρα ρα 0 ρα ρα 0
R , and T .
ρβ ρβ 2 ρβ
0 0
ρβ ρβ ρβ ρβ
265 November, 2014 (DH)
Schoenberg-Protázio AVO
Schoenberg and Protázio (1992) extended the isotropic
equations to the orthorhombic case in the following way:
eP1 eS 1
X
s e
13 1 P1 33 3 P P 3
c c s e c s e
13 1 S 1 c s e
33 3 S S 3
c55 s1eP 3 s3 P eP1 c55 s1eS 3 s3S eS 1
Y
eP3 eS3
sin q
where s1 , and s3 P and s3S are the eigenvalue solutions,
VP (q )
and eP1, eP 3, eS 1, eS 3 , the associated eigenvecto rs of the
following matrix :
c11s12 c55s32 ρ c55 c13 s1s3 e1 0
c55 c13 s1s3 c55s1 c33s3 ρ e3 0
2 2
266
November, 2014 (DH)
Schoenberg-Protázio AVO
The orthorhombic case is even more difficult to visualize
than isotropic AVO. To simplify the equations, consider the
normal incidence case, where s1 = 0:
1 1
c55 s32 r 0 e1 0 c33 c55 2 2
s3 ,
c33s3 r e3 0 r r
2
0
The two solutions we want are the positive roots, and
solving for the eigenvectors, we get eP=(0,1) and eS=(1,0).
Substituting into X and Y, we simplify to the isotropic case:
0 1 0 1 0 c55 s3 S 0 r
X Y
c33s3 P 0 r 0 1 0 1 0
267 November, 2014 (DH)
Appendix 5: Rock Physics
of Fractures
268 November, 2014 (DH)
Appendix 5: Rock Physics
In this appendix we establish the link between anisotropy and
fluid filled fractures
We first review linear slip deformation theory as a means to
model fractures in an otherwise homogeneous background
media. The case of a single vertical dry fracture is discussed
first. Parameterizations of differing complexities are
considered. The concepts of the previous chapter are used
to show how fractures introduce anisotropy.
The special case of penny-shaped cracks is then considered.
This case is important in that some seismic attributes assume
fractures of this form.
Then multiple fractures are considered.
Lastly, we show how the Gassmann equation can be
generalized to anisotropic media so that the fluid filled
fracture response might be calculated
269 November, 2014 (DH)
Compliances
Hooke’s law may be written in terms of stiffness C
1
– or the compliance matrix : S C
– so that : S
270 November, 2014 (DH)
Fracture rock physics: Linear slip theory
We model the anisotropy introduced by the fractures using linear
slip deformation theory (Schoenberg, 1980)
In this theory the fractures perturb the background rock acting as a
source of extra strain.
Mathematically the total compliance of the fractured rock is the sum
of the background compliance plus the extra compliance due to the
fractures
S Sb S f
where :
Sb is the (isotropic) background compliance matrix
Sf is the excess compliance due to fractures
271 November, 2014 (DH)
LSD Fracture compliance matrix models
For dry parallel fractures and aligned cracks:
Asymmetric fractures (Schoenberg & Douma, 1988)
• three fracture compliance parameter and two orientation parameters
Simplicity of parameterization
• Orthorhombic anisotropy
Rotationally invariant fractures (Schoenberg and Sayers, 1995)
• Two fracture compliance parameter and two orientation parameters
• Vertical fractures HTI anisotropy
Penny-shaped cracks (Hudson et al, 1981; Kachanov, 1993 )
• Two fracture compliance parameters are a function of one parameter for dry
cracks
Scalar fractures (Schoenberg and Sayers, 1995)
• One fracture compliance parameter and two orientation parameters
Multiple (dry) fracture sets: (Sayers & Kachanov, 1995)
scalar fractures lead to orthorhombic anisotropy
Fluid filled fractures: note that the above theory is all for the dry case. The saturated
case can be derived using the anisotropic form of Gassmann’s equations, which is
discussed in Appendix 5.
272 November, 2014 (DH)
Fracture compliance: asymmetric fractures
Dry asymmetric fractures can be modeled as an imperfectly
bonded interface where N BN N
– Traction is continuous but displacement might be
discontinuous V BV V
– Displacement discontinuity is linearly related to the traction H BH H
For a flat dry crack aligned with the reference coordinate system
– fractures described by 3 scalar fracture parameters
– BN normal fracture compliance
– BV vertical fracture compliance
V
– BH horizontal fracture compliance H
– Orthorhombic symmetry
N N
Fracture may be generalized to any orientation by rotating it with a
bond transformation V H
– introduces 2 extra parameters (dip, azimuth)
273 November, 2014 (DH)
Fracture compliance: dry parallel fractures
Dry parallel fractures give rise to orthorhombic symmetry.
Described by 3 fracture parameters and 2 background (isotropic)
parameters
Total of 5 free parameters instead of 9 parameters for the case of
general anisotropy. There is additional parameters when fluids are
considered.
S-wave velocity is different along each of 3 coordinate axes
P-wave velocity
– is the constant in plane perpendicular to symmetry axis.
– is slower in the direction of the symmetry axis.
274 November, 2014 (DH)
Fracture compliance matrix models
For dry parallel fractures and aligned cracks
– Asymmetric fractures (Schoenberg & Douma, 1988)
• 3 fracture compliance parameter & 2 orientation parameters
• Orthorhombic anisotropy
Simplicity of parameterization
– Rotationally invariant fractures (Schoenberg & Douma, 1988)
• 2 fracture compliance parameter & 2 orientation parameters
• Vertical fractures HTI anisotropy
– Penny-shaped cracks (Hudson et al, 1981; Kachanov, 1993 )
• 2 fracture compliance parameters are a function of 1 parameter for dry cracks
– Scalar fractures (Schoenberg & Sayers, 1995)
• 1 fracture compliance parameter & 2 orientation parameters
– Multiple (dry) fracture sets (Sayers & Kachanov, 1995)
– scalar fractures lead to orthorhombic anisotropy
– Fluid filled fractures
275 November, 2014 (DH)
Fracture compliance: Rotationally invariant fracture sets
For rotationally invariant fracture sets the vertical
and horizontal compliances are the same
BT BV BH
2 scalar fracture parameters T
T
– BN normal fracture compliance N BN N
– BT tangential fracture compliance N N
T BT T
T T
Vertical fractures HTI symmetry
Fractures may be generalized to any orientation by
rotating it with a bond transformation
– introduces 2 extra parameters (dip, azimuth)
– TTI symmetry
276 November, 2014 (DH)
Fracture compliance: Rotationally invariant fracture sets
Linear slip deformation (LSD) theory (Schoenberg, 1980)
Fractures perturb the background compliance matrix :
S Sb S f
where :
Sb is the isotropic background compliance matrix :
BN 0 0 0 0 0
0 0
Sb Cb1 0 0 0 0
0 0 0 0 0 0
Sf
Sf is the excess compliance due to fractures : 0 0 0 0 0 0
0 0 0 0 BT 0
0 0 0 0 0 BT
Single set of vertical fractures
embedded in an isotropic background
277
November, 2014 (DH)
Fracture stiffness: Rotationally invariant fracture sets
Rotationally invariant fracture sets added to an isotropic background results in:
M b 1 d N lb 1 d N lb 1 d N 0 0 0
l 1 d M 1 2d
b
lb 1 bd N 0 0 0
N b b N
1
lb 1 d N lb 1 bd N M b 1 b 2d N 0 0 0
CS
0 0 0 mb 0 0
0 0 0 0 mb 1 d T 0
0 0 0 0 0 mb 1 d T
where :
– normal weakness impacts P-wave velocity :
– function of fluid, crack density and aperture
M b BN
dN , d N [0,1)
– tangential weakness impacts S-wave velocity : 1 M b BN
– function of crack density
mb BT
and with: b lb / M b dT , d T [0,1)
1 mb BT
M b lb 2mb
278 November, 2014 (DH)
HTI media
M b 1 d N lb 1 d N lb 1 d N 0 0 0
b
l 1 d M 1 2d lb 1 bd N 0 0 0
N b b N
1 l 1 d N lb 1 bd N M b 1 b d N
2
0 0 0
A b
r 0 0 0 mb 0 0
0 0 0 0 mb 1 d T 0
0 0 0 0 0 mb 1 d T
Slow P-wave
Slow S-wave b lb / M b
M b lb 2mb
After Rüger (1997)
279 November, 2014 (DH)
HTI media
M b 1 d N lb 1 d N lb 1 d N 0 0 0
b
l 1 d M 1 2d lb 1 bd N 0 0 0
N b b N
1 l 1 d N lb 1 bd N M b 1 b d N
2
0 0 0
A b
r 0 0 0 mb 0 0
0 0 0 0 mb 1 d T 0
0 0 0 0 0 mb 1 d T
Fast P-wave Fast S-wave b lb / M b
M b lb 2mb
After Rüger (1997)
280 November, 2014 (DH)
Fracture compliance matrix models
For dry parallel fractures and aligned cracks
– Asymmetric fractures (Schoenberg & Douma, 1988)
• 3 fracture compliance parameter & 2 orientation parameters
• Orthorhombic anisotropy
Simplicity of parameterization
– Rotationally invariant fractures (Schoenberg & Douma, 1988)
• 2 fracture compliance parameter & 2 orientation parameters
• Vertical fractures HTI anisotropy
– Penny-shaped cracks (Hudson et al, 1981; Kachanov, 1993 )
• 2 fracture compliance parameters are a function of 1 parameter for dry cracks
– Scalar fractures (Schoenberg & Sayers, 1995)
• 1 fracture compliance parameter & 2 orientation parameters
– Multiple (dry) fracture sets (Sayers & Kachanov, 1995)
– scalar fractures lead to orthorhombic anisotropy
– Fluid filled fractures
281 November, 2014 (DH)
Penny-shaped cracks (Hudson theory)
For the special case of penny-shaped rotationally invariant cracks Hudson et al., (1981)
showed:
4
d Ndry
3g 1 g
16
dT
33 2 g
3f
where crack density
4a
a: aspect ratio
mb2
g square of S - wave to P - wave velocity ratio
Mb 2
Derived treating fractures as a linear perturbation of stiffness.
Note dN and dT are a function of a single variable crack density.
282 November, 2014 (DH)
Fluid filled penny-shaped cracks (Hudson theory)
For fluid filled penny-shaped cracks Hudson et al., (1981) showed
0.6
16
d Tsat d Tdry
0.5
33 2 g
DELTA T
0.4
0.3
3f
where
0.2
crack density
0.1 4a
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
crack density
gas
4
hudson wet
0.6
Gassmann wet
d Nsat d Ndry
3g 1 g
0.5
DELTA N
0.4
0.3
1
where and 0 1
liq
0.2
1 M
1
ag 1 g M 0
0.1
0
0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1
crack density
The fluid term acts as a fractional scalar to the dry
response
283 November, 2014 (DH)
Penny-shaped cracks (Hudson theory)
Linear Hudson theory:
– For dry cracks dN and dT are a function of one variable, the crack density.
– Fluid filled cracks dN and dT are a function of two variables.
• Crack density
• Fluid term which is a function of the fluid modulus and aspect ratio
– Fractures treated as linear perturbation of stiffness in contrast with rest of the
section which treats fractures as perturbations of compliance
• In this course only looking at linear theory, there is also a nonlinear extension
• In future sections we find that it is easy to parameterize seismic attributes in
terms of Hudson theory
Next, we shall look a penny-shaped cracks derived assuming linear perturbations of
compliance following Kachanov (1992).
284 November, 2014 (DH)
Penny-shaped cracks (Kachanov)
Using linear slip theory Kachanov (1992, 1993) derived relationship linking
BN and BT for penny-shaped cracks:
BN
16a 1 n b2
3Eb
1
BT BN
nb
1
where 2
– nb : background Poisson ratio
– Eb : background Young’s modulus
– a: radius of crack
fracture compliances are a function of only 1 free parameter
• crack radius
• fluid filled penny-shaped cracks are also a function of aspect ratio
Vertical penny shaped cracks HTI symmetry
285 November, 2014 (DH)
Fluid filled penny-shaped cracks (Kachanov)
For flat aligned cracks Gassmann (1951) theory can be used to show that the saturating
fluid does not affect the tangential fracture compliance:
BTSat BTdry
The saturating fluid decreases the normal fracture compliance (O'Connell and Budiansky,
1974).
BNSat BNdry where 0 1
where : Aperture
1
E b : Background Young' s modulus
Eb
1 31 2n b n b : Background Poisson ratio
K
f K f : Fluid modulus
Need at least 2 parameters to describe a fluid filled flat fracture
286
November, 2014 (DH)
Penny-shaped cracks (summary)
Penny shaped cracks parameterization:
– Linear slip theory written in terms of fracture compliances BN and BT.
– Hudson theory written in terms of fracture weaknesses dN and dT .
The two parameterizations are related by the nonlinear relation:
MBN mBT
dN , dT
1 MBN 1 mBT
dN dT
BN , BT
M 1 d N m 1 d T
Seismic attributes can be expressed in terms of dN and dT so it is easier to use
Hudson theory.
Grechka and Kachanov (2006) show through finite element modeling experiments
that treating fractures as a perturbation in compliance is more accurate than using
stiffness.
287 November, 2014 (DH)
Fracture compliance matrix models
For dry parallel fractures and aligned cracks
– Asymmetric fractures (Schoenberg & Douma, 1988)
• 3 fracture compliance parameter & 2 orientation parameters
Simplicity of parameterization
• Orthorhombic anisotropy
– Rotationally invariant fractures (Schoenberg & Douma, 1988)
• 2 fracture compliance parameter & 2 orientation parameters
• Vertical fractures HTI anisotropy
– Penny-shaped cracks (Hudson et al, 1981; Kachanov, 1993 )
• 2 fracture compliance parameters are a function of 1 parameter for dry cracks
– Scalar fractures (Schoenberg & Sayers, 1995)
• 1 fracture compliance parameter & 2 orientation parameters
– Multiple (dry) fracture sets (Sayers & Kachanov, 1995)
– scalar fractures lead to orthorhombic anisotropy
– Fluid filled fractures
288 November, 2014 (DH)
Scalar cracks
For dry penny-shaped cracks:
n
BN 1 b BT
2
thus for rocks with a small background Poisson ratio
BN BT
– Fractures where BN BT are called scalar (Schoenberg & Sayers, 1995).
– Again the fracture compliance matrix is only a function of 1 free
parameter.
– This approximation is useful to help simplify the modeling and analysis of
multiple fractures.
289 November, 2014 (DH)
Fracture compliance matrix models
For dry parallel fractures and aligned cracks:
– Asymmetric fractures (Schoenberg & Douma, 1988)
• 3 fracture compliance parameter & 2 orientation parameters
Simplicity of parameterization
• Orthorhombic anisotropy
– Rotationally invariant fractures (Schoenberg & Douma, 1988)
• 2 fracture compliance parameter & 2 orientation parameters
• Vertical fractures HTI anisotropy
– Penny-shaped cracks (Hudson et al, 1981; Kachanov, 1993 )
• 2 fracture compliance parameters are a function of 1 parameter for dry cracks
– Scalar fractures (Schoenberg & Sayers, 1995)
• 1 fracture compliance parameter & 2 orientation parameters
– Multiple (dry) fracture sets (Sayers & Kachanov, 1995)
– scalar fractures lead to orthorhombic anisotropy
– Fluid filled fractures
290 November, 2014 (DH)
Fracture compliance: Multiple fractures
Using the noninteraction approximation, that the
interactions of different fractures in the stress field
can be ignored, then:
N
S S b S fr
r 1
S fr
where is the kth fracture
Each fracture is potentially described by 3 fracture
compliances & 2 direction parameters (dip,
azimuth).
• Simplification: multiple scalar fractures lead to (Sayers, 2009)
orthorhombic symmetry.
291 November, 2014 (DH)
Fracture compliance matrix models
For dry parallel fractures and aligned cracks:
– Flat cracks (Schoenberg & Douma, 1988)
• 3 fracture compliance parameter & 2 orientation parameters
Simplicity of parameterization
• Orthorhombic anisotropy
– Rotationally invariant fractures (Schoenberg & Douma, 1988)
• 2 fracture compliance parameter & 2 orientation parameters
• Vertical fractures HTI anisotropy
– Penny-shaped cracks (Hudson et al, 1981; Kachanov, 1993 )
• 2 fracture compliance parameters are a function of 1 parameter for dry cracks
– Scalar fractures (Schoenberg & Sayers, 1995)
• 1 fracture compliance parameter & 2 orientation parameters
– Multiple (dry) fracture sets (Sayers & Kachanov, 1995)
– scalar fractures lead to orthorhombic anisotropy
– Fluid filled fractures
292 November, 2014 (DH)
Fluid filled isotropic rocks
Recall for isotropic media Gassmann’s (1951) equation can
be written as:
K sat K dry 2 M
m sat mdry
m sat mdry
Indicates fluid stiffens rock proportional to and M.
– both are functions of porosity
K dry
= Biot coefficient 1 ,
Km
– function of pore compressibility
M = combined modulus of fluid and rock
1 f f
,
M K fl Km
November, 2014 (DH)
293
Fluid filled anisotropic rocks
For anisotropic media Gassmann’s (1951) equation can be written as:
cijsat cijdry i j M
Indicates fluid stiffens rock proportional to and M.
– both are functions of porosity m sat mdry
*
j = anisotropic Biot coefficient 1
*
K dry
,
– function of pore compressibility and fracture compliance Km
M = combined modulus of fluid and rock
1 f * f
,
M K fl Km
This is developed in greater detail in Appendix B.
294 November, 2014 (DH)
Fracture summary
Fractures can be modeled using linear slip theory (Schoenberg, 1980).
Flat thin fractures need 1 to 3 parameters to describe them.
• The fractures can be described in terms of fracture compliances or fracture
weaknesses:
• Fracture compliances have units of compliance
• Fracture weaknesses are fractional values
• If the fracture is fluid-filled, at a minimum one needs 2 parameters to describe
the fracture.
The response due to multiple fractures can be calculated using the NIA summing
the individual fracture compliances.
The Gassmann (1951) equation extends naturally to anisotropic media.
• The presence of fractures amplifies the fluid response.
• In the special case of a single vertical fracture set, fluids influence the response
of epsilon but not gamma.
295 November, 2014 (DH)
Appendix 6: BN/BT ratio
296 November, 2014 (DH)
BN/BT ratio (Verdon and Wüstefeld, 2013)
BN/BT
BN/BT=1
Figure 1 Published values of BN/BT from laboratory and field studies.
297 November, 2014 (DH)
Label numbers correspond the studies listed in Table 1.
BN/BT ratio (Verdon and Wüstefeld, 2013)
298 November, 2014 (DH)
BN/BT ratio - Penny shaped cracks
𝑩𝑵
For Penny shaped cracks = 𝟏 − 𝝂/𝟐
𝐵𝑻
Relationship between Vp/Vs ratio, g and Poisson ratio
𝑽𝒑 𝟐
𝑽𝒔
−𝟐
𝟏−𝟐𝒈
𝝂= 𝑽𝒑 𝟐
= 𝟐(𝟏−𝒈)
𝟐−𝟐 𝑽
𝒔
Vp/Vs Ratio g Poisson ZN/ZT ratio
Ratio
1.414 1/2 0 1
1.7321 1/3 1/4 0.875
2 1/4 1/3 0.8333
299 November, 2014 (DH)
Appendix 7: The Rüger equation
300 November, 2014 (DH)
Appendix 7: The Rüger equation
In this appendix, we discuss the details for inversion of the
AVO and AVAz equations for isotropic and azimuthally
anisotropic media.
For isotropic media, the data is a function of incident angle
q alone, whereas for azimuthally anisotropic media, the
data is a function of both incident angle q and azimuth
angle f, as shown in the figure below:
After Rüger (1998)
301 November, 2014 (DH)
Appendix 7: The Rüger equation
In isotropic AVO, we pick the q
q1 qN
amplitudes of an N-trace angle 600
gather at each time t, as shown
Time t
on the right. These amplitudes (ms)
can be modeled as: 650
R(qi ) A B sin 2 qi , i 1, , N , where :
R(qi ) the reflection coefficien t at incident angle qi ,
A the intercept, and B the gradient.
R(q1 ) 1 sin 2 q1
R(q )
In matrix form we write
2 1 sin 2
q 2
A
B
this as a set of N linear
equations:
302 R (q N
) 1 sin 2
q N
November, 2014 (DH)
Appendix 7: The Rüger equation
To solve for A and B note that we can write the previous
equation more simply as:
d Gm, where :
R(q1 ) 1 sin 2 q1
R(q )
2 1 sin 2
q 2 A
d ,G , and m .
B
R(q N ) 1 sin q N
2
The least-squares solution to this set of equations is:
mˆ G G GT d , where
T 1
mˆ is the best least - squares estimate of m.
303 November, 2014 (DH)
Appendix 7: The Rüger equation
In anisotropic AVO for HTI media, we start with the Rüger
equation:
R(qi , f j ) Aiso [ Biso Bani sin 2 (f j fiso )] sin 2 qi ,
where :
R(qi , f j ) the reflection coeff. at incident angle qi
and azimuth angle f j , where i 1,, N , and j 1,, M ,
fiso azimuth angle in isotropy plane, Aiso isotropic intercept,
Biso isotropic gradient, and Bani anisotropic gradient.
We thus must invert the equation for four unknowns: Aiso, Biso,
Bani and fiso. This is a non-linear equation for fiso, since this
term is contained within a nonlinear equation. A description
of fiso is given in the next figure.
304 November, 2014 (DH)
Appendix 7: The Rüger equation
Rüger (1998) linearized the Zoeppritz equation for the case
of an isotropic half-space over an HTI half-space. As shown
below, the symmetry-axis is at right angles to the fractures
and the isotropy plane is parallel to the fractures:
fsym
Modified from Rüger (1997)
305 November, 2014 (DH)
Appendix 7: The Rüger equation
To linearize the Rüger equation, we note that:
1 cos(2 )
sin 2 ( ) , and :
2 2
cos( ) cos( ) cos( ) sin( ) sin( )
Substitution for the second term in the brackets in the
Rüger equation gives:
Bani cos(2f ) cos(2fiso ) sin(2f ) sin(2fiso )
Bani sin 2 (f fiso ) Bani
2 2
B
ani C cos(2f ) D sin(2f ),
2
B cos(2fiso ) B sin(2fiso )
where : C ani , and D ani .
306
2 November, 2014 (DH)
2
Appendix 7: The Rüger equation
Thus, we can re-write the full Rüger equation as:
R(qi , f j ) A [ B C cos 2f j D sin 2f j ] sin 2 qi , where :
Bani B cos(2fiso ) B sin(2fiso )
A Aiso , B Biso , C ani , and D ani .
2 2 2
This allows us to set up the following set of linear
equations:
d Gm, where :
R(q1 , f1 ) 1 sin 2 q1 cos 2f1 sin 2 q1 sin 2f1 sin 2 q1 A
R(q , f ) B
1 sin 2
q2 cos 2f1 sin 2 q 2 sin 2f1 sin 2 q 2
d 2 1
,G , m .
C
R (q N , fM ) 1 sin 2
qN cos 2fM sin 2 q N sin 2fM sin 2 q N D
307 November, 2014 (DH)
Appendix 7: The Rüger equation
This again has the solution given by:
mˆ GT G GT d
1
Having solved for m, we extract our four constants as follows:
tan 1 ( D / C )
Aiso A, fiso (since : tan 2fiso D / C ),
2
Bani
Bani 2 C D and Biso B
2 2
B C 2 D2 ,
2
cos 2 (2fiso ) sin 2 (2fiso )
2
1 Bani
since C D
2 2 2
Bani .
2 2
Since Bani is the solution of a square root, there is a sign
ambiguity. This is solved using a different technique.
308 November, 2014 (DH)
Appendix 7: The Rüger equation
Often, the least-squares solution given in the previous
slides can be unstable due to data errors. One way to
stabilize the solution is to add pre-whitening, as follows:
mˆ G G lI G T d , where :
T 1
l is a pre - whitening factor, and
I the M x M identity matrix.
Note that in the pre-whitening approach the same weight
is applied to each diagonal element of the GTG matrix.
As the pre-whitening value gets larger, the solution
becomes more stable but less accurate.
However, in most cases, we would prefer a set of values
which apply different weights to different values in GTG.
This more complete solution is shown in the next slide.
309 November, 2014 (DH)
Appendix 7: The Rüger equation
A more complete solution is to compute and apply two
separate weighting matrices, one for the data and one
for the model, as shown here:
mˆ G C G C
T 1
d m
1 1
G T Cd1d , where :
Cd is an NxN data weigh t matrix, and
Cm is an MxM model weight matrix.
Of course, the problem now is to find the best forms of the
Cd and Cm matrices, which is a difficult task.
For Cm we build a set of values that are related to the
statistics of the model data.
For Cd we design an iterative solution which updates the
weights after each iteration based on the error between the
observed and predicted data.
310 November, 2014 (DH)
Appendix 7: The Rüger equation
Since we have four model parameters, Cm is a 4 x 4 matrix
with weights that are related to those parameters (We
define the parameter T later in these notes):
2 2
1 1 rps 0 0
m
2 2 8 2
1 rps 1 ( 1) rps 2
0 0
Cm R2 p ,
m m m
0 5
0 0
2( 256 2 )
5
0 0 0
2( 256 2 )
R2 p variance of P - impedance reflectivi ty, m mudrock slope,
2
rps correlatio n coefficien t between the p and s data, and 2 . T
R p
311 November, 2014 (DH)
Appendix 7: The Rüger equation
In the iterative re-weighted least-squares (IRLS) approach,
we compute the error between the observed and predicted
data after each iteration. For the first iteration we get:
mˆ G G lI G d dˆ Gmˆ
T 1 T
1 1 1
e11 d1 dˆ11
e d ˆ
ˆ
e1 d d1 12
2
d12
,
ˆ
e1N d N d1N
where e1 is the error after the first iteration.
The above equation shows us that after each iteration we
have N error terms, one for each input data value.
Each subsequent iteration gives a new error, where ej is
th iteration.
312the error at the j November, 2014 (DH)
Appendix 7: The Rüger equation
One of the earliest approaches to the IRLS approach was
called the lp norm method, where the weight matrix at the jth
iteration is given by:
e p 2 0 0
j1 p 2
0 e j1 0 0
Cdj1 ,
p 2
0 0 e jN
where means take the absolute value of
the error and p a value between 1 and 2.
When p = 2, the weight matrix reduces to the identity
matrix, and when p = 1, each term is equal to 1/|eji|, where i
= ith term, meaning that larger errors are weighted down.
313 November, 2014 (DH)
Appendix 7: The Rüger equation
The approach used here is called the Cauchy-Gauss
method, were the weights at the jth iteration are given by:
e 2
1
j1 2 1 0 0
2 nj
1
e j2 2
1
Cdj 0 2 1 0 0 ,
nj
2
e jN 2
1
0 0 2 1
2 nj
1 N 2
where nj e ji (subscript n means the noise term).
2
N i 1
314 November, 2014 (DH)
Appendix 7: The Rüger equation
Our complete solution to the Rüger equation is therefore given by
iterating J times through the following set of equations:
1
T 1 2
σ T 1
mˆ j G Cdj G nj
Cm G Cdj d , j 1,, J ,
1
2
σ
Rp
where : Cd11 I N , the NxN identity matrix.
The individual matrices in the above equation have all been defined
earlier, except that you will note the extra scale factor in front of the
model constraint.
Typically, we find that because of the large number of data values (i.e. N
is very large) we only need a couple of iterations to converge to a stable
value of the model parameters.
315 November, 2014 (DH)
Appendix 8: Azimuthal Geometric Spreading
Correction
316 November, 2014 (DH)
Azimuthal Geometric Spreading Correction
Orthorhombic geometric spreading correction (Xu and
Tsvankin, 2005)
– ~Ursin-style correction with V()
X-direction Y-direction
Vp1 Vp1
Vp2x >Vp1 Vp2y>Vp2x >Vp1
Vpx ≠Vpy in one or more layers
x
317 November, 2014 (DH)
Geometric divergence losses
P-velocity S-velocity Vp/Vs density
0 0 0 0
200 200 200 200
400 400 400 400
Depth (m)
600 600 600 600
800 800 800 800 Divergence induced AVO
1000 1000 1000 1000
1200 1200 1200 1200
1400 1400 1400 1400
0 5000 10000 0 2000 4000 1 2 3 1 2 3
m/s m/s g/cc
Note how the geometrical divergence loss
introduces an AVO effect
318 November, 2014 (DH)
Anisotropic Velocities
Hampson-Russell View3D
319 November, 2014 (DH)
Azimuthal Geometrical Divergence Gain Correction
• Travel-times from
• Azimuthal (HTI or orthorhombic)
• Stacking velocity & higher order terms (isotropic or VTI)
• Stacking velocity (isotropic)
-8000
8000
Divergence correction Divergence correction Divergence correction
(0 degree azimuth) (45 degree azimuth) (90 degree azimuth)
320 November, 2014 (DH)
After Xu and Tsvankin (2005)
Gain correction show as COCA cube for one CDP
November, 2014 (DH) 321
Location 3: input
400 ms
700 ms
zero offset
November, 2014 (DH) 322
Location 3: output with azimuthal scaling correction
1400 ms
1700 ms
zero offset
November, 2014 (DH) 323
Appendix 9: Constraints
324 November, 2014 (DH)
Constrained AVAz inversion
In least squares inversion the optimal solution is determined
by minimizing the squared misfit between the model and the
data.
Sometimes the problem is ill-conditioned (poor S/N for the
amount of data) or undetermined (not enough data to find a
unique solution) so that the inversion becomes unstable or
unsolvable. In these cases we can use constraints to find an
optimal solution which is consistent with our a priori
knowledge.
A priori information is our knowledge of what the solution
should look like.
– For example, the mudrock relationship states that the S-wave velocity
should be linearly related to the P-wave velocity.
– Since our knowledge of what the solution should look like is uncertain it is
best to describe this knowledge using probability theory and distributions.
325 November, 2014 (DH)
Constrained Least Squares
Constraints have been implemented for the 3-term AVO and near-
offset Rüger equations. Constraints should help in areas with low
fold and poor azimuthal distributions.
The constraints are calculated from well log statistics using the
process Calculate AVAz Constraints.
The process Azimuthal AVO Volume then optionally reads these
constraints and uses them to run constrained inversion.
The constrained inversion can be run either using Simple Least
Squares (SLS) or Iteratively Reweighted Least Squares (IRLS).
IRLS identifies and down weights outliers to improve the solution.
The following example illustrates the use of constraints and IRLS.
326 November, 2014 (DH)
Anisotropic gradient: Impact of constraints and weighting on solution
Original With With constraints
constraints and IRLS
and SLS
Constraints can be used to help stabilize the solution. In the above example the
original (unconstrained) inversion of the anisotropic gradient becomes unstable
327
at the left edge of the 3D. Applying constraints and weights improve the solution
November, 2014 (DH)
Constraints in 2D
Constraints describe
probabilistic relationships
between the different
parameters
– The P-wave velocity reflectivity is
related to the S-wave velocity
reflectivity through the mudrock
relationship.
– The mudrock slope, the Vs/Vp
ratio and the correlation coefficient
parameters are sufficient to
describe a Bivariate Gaussian
Probability distribution, which can
be used to constrain the AVO
problem.
328 November, 2014 (DH)
Constraints in 3D
Gardner like relationships
can be used to constrain
the density to track the P-
wave velocity reflectivity
and S-wave velocity
reflectivity leading to a
multivariate Gaussian
distribution.
An equiprobability surface
of this looks like an ellipsoid
in 3D parameter space.
329 November, 2014 (DH)
Constraints
The a priori constraints are described by a multivariate Gaussian
distribution:
1 T Cm 1
Pm | I exp m m
2 h 2
where Cm is the covariance matrix:
2 RVp R R R R
Vp Vs Vp den
Cm RVp RVs RVs RVs Rden
2
2
RVp Rden RVs Rden Rden
and h is the global scalar.
330 November, 2014 (DH)
AVO parameter uncertainty
The optimal solution occurs at the
maximum probability of the prior
probability (blue) and the seismic
misfit probability (black).
The potential negative of using
constraints are that they bias the
solution towards the priors. Ideally
the solution should be driven by the
data. To get around this problem,
the relative weighting of the misfit
and prior functions changes
depending on the S/N ratio.
1
T εT ε ~ 1 T
m G WeG C m G Wed
N 1 Rp
2
where Gm-d
331 November, 2014 (DH)
AVAz constraints
The 3-term AVO constraints are extended to the near
offset Rüger equation by assuming
– The anisotropic gradient is proportional to the crack density
– The crack density is uncorrelated with the density, fast P-wave and S-
wave velocities.
Under these assumptions the anisotropic parameter Bani
is uncorrelated with the isotropic parameters.
– One parameter Lambda is sufficient to describe the relative RMS power
of the Bani trace to the RMS power P-wave impedance reflectivity
332 November, 2014 (DH)
Calculate AVAz Constraints
The Calculate AVAz Constraints process
calculates constraints for the azimuthal
AVO plugin. It calculates a covariance
matrix from the well control.
First, the well that the constraints are going
to be calculated on needs to be specified,
in this case well A was chosen.
Then the P-wave velocity curve, S-wave
velocity curve and a Density curves need
to be specified. The time window may be
defaulted from 0 to 99999 ms. The time
interval controls the blocking. Larger block
sizes reduce the variance. Usually the
sampling process rate is a good default.
Lastly the name under which the
constraint file is going to be stored and
accessed needs to be specified. The
default is the name of the well log.
333 November, 2014 (DH)
Calculate AVAz Constraints
The program is run by clicking OK in the bottom right hand corner.
334 November, 2014 (DH)
Calculating AVAz Constraints
• The curves are converted from depth
to time using the default time to depth
mapping.
• The reflectivity for each of these
curves is calculated.
• Then the covariance matrix is
calculated over the time window
specified from Start Time to End Time.
• The covariance matrix is transformed
to the scaled velocity
parameterization.
• The Covariance Matrix describes the
amount of variance for each of the
key variables and the correlations
between these variables.
335 November, 2014 (DH)
Calculate AVAz Constraints
The program actually calculates the
normalized parameter covariance matrix
for the scaled velocity AVO
parameterization.
The covariance matrix is normalized so
that the C11 term is 1.
• This C11 term is the variance of the P-
wave reflectivity and acts as scaling
factor.
• The standard deviation of the P-wave
reflectivity is shown as a parameter
below the covariance matrix and can be
used as a scale factor for the Azimuthal
AVO volume attributes under the
advanced parameters.
• In addition the average Vs/Vp ratio for
the log is shown
336 November, 2014 (DH)
Alternative Parameterization
The program also calculates key rock physics
relationships from the reflectivity including
the:
• Mudrock slope
• Mudrock slope correlation coefficient
• P-wave Gardner coefficient
• P-wave Gardner correlation coefficient
• S-wave Gardner coefficient
The parameter covariance matrix can be
reconstructed from these parameters.
These parameters are written to a file which
can then be read by the Azimuthal AVO volume.
337 November, 2014 (DH)
Azimuthal AVO (using constraints)
These constraints can
then be used in both the
Azimuthal AVO Analysis
and Azimuthal AVO
Volume processes.
In order to do this the
user would first specify
yes for Use Constraints.
338 November, 2014 (DH)
Azimuthal AVO Volume (using constraints)
Then, under the “Advanced Parameter Tab,
Constraint parameters”, They would select
the previously constructed constraints
using previous output, in this example well
A
This reads all the parameters previously
calculated by the Calculate AVAz
Constraints process. The user is free to
change these.
The Azimuthal AVO Volume Process then
calculates the covariance matrix from these
parameters. Different constraint matrices
are calculated depending on the AVO
Parameterization.
The variable Lambda controls the relative
size of the Bani variance compared to the P-
wave impedance reflectivity variance. A
value of 1 means they are the same size.
If well logs are not available, reasonable
default values may be used.
339
November, 2014 (DH)
Appendix10 : filtering Azimuths
340 November, 2014 (DH)
Exercise 4 – Azimuthal analysis
Performing the
azimuthal analysis
results in a
magnitude
attribute
(Anisotropic
gradient) and a
azimuth attribute
(isotropy plane
azimuth).
The azimuth only
makes sense if
the media is
anisotropic. If the
media is isotropic
then the azimuth
is meaningless
and should be
excluded.
341 November, 2014 (DH)
Exercise 4 – Azimuth Display
The next several slides
demonstrate a Trace Math
script to only retain the
azimuths when the media is
anisotropic.
The idea is to make the
azimuths NULL if the media is
isotropic (e.g. if the anisotropic
gradient is less than some user
defined threshold (e.g. 0.02))
If it is smaller than this, then the
media is thought to be isotropic
and the azimuth is defined a
null value (which we define as -
999).
Go to the Processes menu and
select Trace Maths and select
the parameters to the right.
Then click on “Import and Edit”
to create the script following the
text on the next slide.
342 November, 2014 (DH)
Exercise 4 – Azimuth Display
Enter the following script, this can be copied and pasted from
this PowerPoint
/* Input Variables */
/* Bani = Anisotropic Gradient */
/* Azimuth = Azimuth of Isotropy plane */
/* Output */
/* FilteredAzimuth */
/* Initialize the output trace the same size as that of the input */
filteredAzimuth=azimuth;
/* Number of samples */
n = numsamples( azimuth);
i=0;
while (i < n )
{
if (Bani[i] < 0.02)
{
filteredAzimuth[i]=-999; /*output null value*/
}
i=i+1
}
/*Output*/
filteredAzimuth;
After entering the script, test script and if it passes, click OK
and then OK a second time.
343 November, 2014 (DH)
Exercise 4 – Azimuth Display
You should
get a display
that looks
similar to the
right.
Most of the
purple are
null values.
The next
slide shows
how to
modify the
color bar so
as to exclude
the null
values from
the display.
344 November, 2014 (DH)
Exercise 4 – Excluding Null values from the display
Double click on the color below -90
345 November, 2014 (DH)
Exercise 4 – Displaying Azimuth
The null azimuth
values are now
shown in the
background
white color.
Later, we will
show how to
display these
simultaneously,
which is the
subject of the
Visualization
exercise.
346 November, 2014 (DH)
You might also like
- Crooked-Line 2D Seismic Reflection Imaging in Crystalline Terrains: Part 1, Data ProcessingNo ratings yetCrooked-Line 2D Seismic Reflection Imaging in Crystalline Terrains: Part 1, Data Processing12 pages
- Crooked-Line 2D Seismic Reflection Imaging in Crystalline Terrains: Part 2, MigrationNo ratings yetCrooked-Line 2D Seismic Reflection Imaging in Crystalline Terrains: Part 2, Migration11 pages
- Pe20m017 Vikrantyadav Modelling and InversionNo ratings yetPe20m017 Vikrantyadav Modelling and Inversion18 pages
- TIME DEPTH CONVERSION - GROUP 3 - 21jan16No ratings yetTIME DEPTH CONVERSION - GROUP 3 - 21jan1624 pages
- What Is Seismic Interpretation? Alistair BrownNo ratings yetWhat Is Seismic Interpretation? Alistair Brown4 pages
- Use of Seismic Attributes For Sediment Classificat PDFNo ratings yetUse of Seismic Attributes For Sediment Classificat PDF15 pages
- Advance Processing: Institut Teknologi BandungNo ratings yetAdvance Processing: Institut Teknologi Bandung190 pages
- Seismic Data Processing (Part Iii) MigrationNo ratings yetSeismic Data Processing (Part Iii) Migration64 pages
- Free Workshop On Why Seismic Imaging DEC PEANo ratings yetFree Workshop On Why Seismic Imaging DEC PEA2 pages
- Seismic Facies Analysis Applied To P and S Impedances From Pre-Stack InversionNo ratings yetSeismic Facies Analysis Applied To P and S Impedances From Pre-Stack Inversion4 pages
- 1 Fundamentals of Geophysical Data FGDPNo ratings yet1 Fundamentals of Geophysical Data FGDP276 pages
- 2007 - Land Seismic Data Acquisition and ProcessingNo ratings yet2007 - Land Seismic Data Acquisition and Processing214 pages
- Offset Dependent Mis-Tie Analysis at Seismic Line IntersectionNo ratings yetOffset Dependent Mis-Tie Analysis at Seismic Line Intersection11 pages
- MIGRATION - Overview - Fundamentals of Imaging - Brian Schulte - Talisman Energy - Canada - CSEG - 2012-12-RECORDERNo ratings yetMIGRATION - Overview - Fundamentals of Imaging - Brian Schulte - Talisman Energy - Canada - CSEG - 2012-12-RECORDER9 pages
- LH1 - Sources 1 - Vibroseis Theory 1 6341002 02No ratings yetLH1 - Sources 1 - Vibroseis Theory 1 6341002 0229 pages
- 2d 3d 4d Exploracion Sismica para Geologos CH Liner Conex - Mayo 1999No ratings yet2d 3d 4d Exploracion Sismica para Geologos CH Liner Conex - Mayo 1999388 pages
- 2 Petroleum Geology and Geophysics BookNo ratings yet2 Petroleum Geology and Geophysics Book139 pages
- Y014 Advanced 3D Land Internal Mutiple Modeling and Subtraction, A WAZ Oman Case StudyNo ratings yetY014 Advanced 3D Land Internal Mutiple Modeling and Subtraction, A WAZ Oman Case Study5 pages
- List of Contents Subject Page No.: I - Chapter I: Pre - Processing & Sotring Data 1 - Introduction . 1100% (1)List of Contents Subject Page No.: I - Chapter I: Pre - Processing & Sotring Data 1 - Introduction . 112 pages
- Wide-Azimuth Processing For Azimuthal Anisotropy Analysis: Cynthia Gomez, Erika Angerer, CGG Technology LondonNo ratings yetWide-Azimuth Processing For Azimuthal Anisotropy Analysis: Cynthia Gomez, Erika Angerer, CGG Technology London4 pages
- Fractured Reservoir Characterization From Seismic Azimuthal Anisotropy and Well Log AnalysisNo ratings yetFractured Reservoir Characterization From Seismic Azimuthal Anisotropy and Well Log Analysis4 pages
- QVOA Analysis As An Instrument For Fracture CharacterizationNo ratings yetQVOA Analysis As An Instrument For Fracture Characterization7 pages
- '14 Nascimiento Et Al., High-Resolution Acoustic Impedance Inversion To Characterize TurbiditesNo ratings yet'14 Nascimiento Et Al., High-Resolution Acoustic Impedance Inversion To Characterize Turbidites12 pages
- '04 Buenafama, Acoustic Impedance Inversion of Lower Permian Carbonate Buildups in Permian Basin, TexasNo ratings yet'04 Buenafama, Acoustic Impedance Inversion of Lower Permian Carbonate Buildups in Permian Basin, Texas99 pages
- '07 Chopra, Seismic Attributes For Fault Fracture CharacterizationNo ratings yet'07 Chopra, Seismic Attributes For Fault Fracture Characterization5 pages
- '05 Tingdahl Semi-Automatic Detection of Faults in 3D Seismic DataNo ratings yet'05 Tingdahl Semi-Automatic Detection of Faults in 3D Seismic Data10 pages
- '01 Cacas Nested Geological Modelling of Naturally Fractured Reservoir - UnlockedNo ratings yet'01 Cacas Nested Geological Modelling of Naturally Fractured Reservoir - Unlocked10 pages
- IAT-II Question Paper With Solution of 18ME651 Non-Conventional Energy Sources June-2021-Dr. Harish BabuNo ratings yetIAT-II Question Paper With Solution of 18ME651 Non-Conventional Energy Sources June-2021-Dr. Harish Babu6 pages
- Time Allowed: 35 Minutes: Differentiate Each of The Following Functions With Respect To - (2 Marks)No ratings yetTime Allowed: 35 Minutes: Differentiate Each of The Following Functions With Respect To - (2 Marks)12 pages
- Reading 2 Organizing, Visualizing, and Describing Data - AnswersNo ratings yetReading 2 Organizing, Visualizing, and Describing Data - Answers51 pages
- ClassXPreboard1 Schedule and Blue Print 202223No ratings yetClassXPreboard1 Schedule and Blue Print 2022238 pages
- Particle Symbol Charge Relative Mass Location: The AtomNo ratings yetParticle Symbol Charge Relative Mass Location: The Atom4 pages
- Variational Principle For A Particle in A BoxNo ratings yetVariational Principle For A Particle in A Box4 pages
- Influence of Cavity Geometry On The Fracture Strength of Dental Restorations: Finite Element StudyNo ratings yetInfluence of Cavity Geometry On The Fracture Strength of Dental Restorations: Finite Element Study17 pages
- Practical Transformers and Standard Ratings100% (1)Practical Transformers and Standard Ratings30 pages
- Probability On Graphs Random Processes On Graphs and Lattices Second Edition Grimmett All Chapter Instant Download100% (5)Probability On Graphs Random Processes On Graphs and Lattices Second Edition Grimmett All Chapter Instant Download55 pages
- Mcqs For Intermediate Students. Physics Chap 11 (HEAT)No ratings yetMcqs For Intermediate Students. Physics Chap 11 (HEAT)4 pages
- An Experimental Investigation of Cavitation SurgeNo ratings yetAn Experimental Investigation of Cavitation Surge8 pages