Jean-Laurent Mallet. Geomodeling
Подождите немного. Документ загружается.
418
CHAPTER
8.
ELEMENTS
OF
STRUCTURAL GEOLOGY
• The
dilatation
coefficient
OT
depends
only
on the
determinant
of GT and is
defined
by
• The
components
of the
principal
strain
directions
in the
(Ur,
VT)
frame
asso-
ciated
with
the
triangle
T can be
computed
easily
thanks
to
equation
(6.26):
To
analyze
the
deformations
and as
shown
in figures
(8.22)
and
(8.23),
the
dilatation
coefficient
and the
principal strain directions modulated
by the
elongations
defined
by
equation (8.38)
can be
mapped
on
T(H}.
Strains
versus
curvatures:
A
paradox
The
intuition
may
suggest
that
the
principal normal curvatures
of a
surface
should
be
linked with
the
strain tensor
of
this
surface. Unfortunately,
as
demonstrated
below,
this
is not
true
at
all:
•
Consider
a flat
surface
<S[o]
that
is
bent
without
plastic
deformation
to
obtain
a new
shape
S. In
this
case,
the
strain
tensor
of S
vanishes
while
its
principal
normal
curvatures
do not
vanish.
•
Consider
a flat
square
surface
<S[oj
that
is
transformed
(with
plastic
deforma-
tions)
into
a flat
rectangular
surface
S. In
this
case,
the
strain
tensor
of S
does
not
vanish
while
its
principal
normal
curvatures
vanish.
However,
it is
important
to
note
that
this paradox
is due to the
fact
that
we
are
considering
the 2D
strain
tensor
of a
surface
and not the 3D
strain
tensor
of the
volume
bounded
by
this
surface.
For
example,
if the
surface
5[o]
corresponds
to the top
horizon
of a flat
geological layer which
is
transformed
into
a
folded
surface
S
after
some tectonic event, then
it can be
shown
(e.g.,
see
Sokomikoff
[204],
p 198 and
212)
that
the
components
of the 3D
strain
tensor within
the
layer
are
actually
functions
the
principal normal curvatures
of
the top
horizon.
8.6
Unfolding
a
stack
of
layers
8.6.1
Introduction
One of the
most challenging issues
in
numerical structural geology
is the
"balanced"
unfolding/unfaulting
of a
stack
of
sedimentary layers.
In
practice,
this step
is
required
as a
"preprocessing"
by a
number
of
applications such
as,
for
example,
sedimentary process
and fluid-flow
modeling.
As
shown
in figure
(8.24),
the
goal
is to
"restore"
the
stack
of
layers
in
such
a way
that
8.6. UNFOLDING
A
STACK
OF
LAYERS
419
Figure
8.24 Vertical cross section showing
the
unfolding
mechanism
for a
stack
of
layers.
The
constraints consist
of a
given
field of
restoration vectors
r(x[i](a))
associated
with
the
uppermost horizon
H\±\
in
addition
to the
Mass-Preservation
and
the
Minimum-Deformation constraints, which
specify
that
the
mass
of
terrains
must
remain constant while
the
layers
must
be
deformed
as
little
as
possible.
• the
horizon
H^
corresponding
to the
hanging wall
of the top
layer
is
trans-
formed
into
an
unfolded/unfaulted
surface
H[Q]
located
in a
horizontal plane
U
at
altitude
z
u
,
as it was
assumed
to be at the
time
of
deposition before
being deformed
by a
tectonic event;
• the
terrains located below
H^
follow
the
movement
of
H[i]
toward
ff[
0
]
in a
piecewise
continuous way;
• the
faults,
if
any,
are the
only discontinuities
in the
displacement
of the
terrains
located below
H[i].
In
other
words,
for any
particle
a
located
at
point
x(a)
=
X[
1
](o:)
in the
3D
studied
domain,
we
want
to go
backward
through
geological
time
from
the
current
position
X[i](a;),
which
corresponds
to the
current
normalized
geological
time
t = 1, to the
"initial"
position
X[
0
](a)
at
time
t = 0 of the
deposition
of
H[^:
X[i]
(a)
=
current position
of a
X[
t
]
(a)
==
position
of a at
time
t G
[0,1]
X[
0
]
(a)
=
initial position
of a at
deposition time
By
definition,
the
vectors
r(x[i](o;)),
denned
as
follows
at any
point
X[i](a)
of
the
current
3D
geological
domain
to be
restored,
are
called
"restoration
vectors:"
420
CHAPTER
8.
ELEMENTS
OF
STRUCTURAL GEOLOGY
As
soon
as the
restoration vectors
are
known
at any
location
X[
1
](a),
it is
easy
to
restore
all the
terrains
as
follows:
For the
sake
of
simplicity
and
without loss
of
generality,
as in
section
(8.5),
it
is
assumed here
that
geological time
has
been scaled
in
such
a way
that
• the
current
position
of
terrains
corresponds
to
geological
time
t = 1, and
• the
initial position
of
terrains
at the
time
of
deposition
of
H[\]
corresponds
to
geological time
t — 0.
We
propose
to
consider
the
terrains
to be
restored
as a
compressible viscous
fluid
where
each particle
a
moves backward
from
the
current location
X[
X
]
(a)
to
location
X[
0
](a)
during
the
restoration time
lag
Ai
=
—
1 so
that, according
to the
Taylor series expansion
formula
(see page 147),
we can
write
Consequently,
from
equations
(8.70)
and
(8.71),
it can be
concluded
that
the
opposite
of the
restoration vector
r(x^](a))
can be
viewed
as a first-order
approximation
of
the
velocity
of
the
terrain particle
a at
geological time
t = 1:
It
will
be
shown below
how
this
analogy
to
viscous
fluid
mechanics
can be
used
for
determining
the field of
restoration vectors
r(x[!](o:)).
Notion
of
balanced restoration
"Balanced" restoration [53] means
that
the
mass
of
terrains
has to
remain
constant throughout geological time. This
is
equivalent
to
saying
that
the ve-
locity
cbc[
t
](a)/cft
of the
terrain particles must honor
the
well-known "equation
of
continuity"
17
where
p(x[
f
]
(a))
is the
density
of the
terrains
at
location
x^]
(a) and
geological
time
t.
Potentially, this partial
differential
equation
has an
infinite
number
17
For
a
complete presentation
of the
continuity
equation,
the
reader
is
referred
to
[190,
206],
for
example.
If the
components
of any
vector
V of the 3D
space relative
to
an
orthonormal system
of
coordinates
(x,y,z)
are
noted
as
(V
x
,
V
y
,
V
z
),
then
it
will
be
recalled
that
div(V)
can be
expressed
as
follows:
8.6.
UNFOLDING
A
STACK
OF
LAYERS
Figure
8.25
Vertical
cross
section
showing
the
iterative
"back-stripping"
algo-
rithm:
after
flattening
the top
horizon,
the top
layer
is
removed.
of
solutions; therefore,
to
ensure
the
uniqueness
of the field of
restoration
vectors
r(xj
1
]),
we
propose introducing
the
three
following
constraints,
at the
very
least:
1)
the field of
restoration
vectors
is
given
on the
hanging
wall
H[i]
of the top
layer;
2)
the
Minimum-Deformation
principle
defined
by
equations
(8.58)
has to be
honored;
3)
the
mass
of the
terrains
has to be
preserved.
In
practice,
a
given horizontal plane
U
is
assumed
to
contain
the
unfolded
state
H[
0
]
of the
hanging wall
H^
of the top
layer
and any one of the
methods
presented
in
section (8.5)
can be
used
for
constructing
H^.
As
illustrated
in
figure
(8.16),
the
restoration vectors
r(x[!])
at any
location
X[ij
€
H^
on the
top
horizon corresponding
to the
above constraint
(8.74)-l
can
then easily
be
deduced
as
follows
from
the
location
of its
image
X[
0
]
G
H^:
Back-stripping
algorithm
In
practice,
as
suggested
in figure
(8.25),
the
unfolding
of a
stack
of
layers
is
performed
iteratively
in
three steps:
1.
unfold
the top
horizon;
2.
unfold
the
stack
of
layers
to
follow
the
movement
of the top
horizon;
3.
remove
the top
layer
(stripping)
and go
back
to
step
(1)
until
all the top
horizons
have
been
unfolded.
421
422
CHAPTER
8.
ELEMENTS
OF
STRUCTURAL
GEOLOGY
This
iterative algorithm, called "back-stripping," allows
the
unfolding mech-
anism
to be
broken down into
a
series
of
time
steps
{At^\
At^
2
),...,
At^},
all
assumed
to be
normalized
and to
have
a
span equal
to
(—1).
Let
r^
be the field of
restoration vectors determined
at
time step
At^.
It
should
be
noted
that
time step
A^
1
)
allows
the top
horizon
H[^
to be
flattened, but the
underlying layers
are
still
affected
by
residual deformations,
which have
to be
evaluated
at
subsequent
time
steps.
Consequently,
if the
stack
of
horizons
is
such
that
H^
1
'
is the
uppermost horizon
H^
is
immediately below
H^
H
(i}
is
immediately below
H
(i
~V
then
• the
deformation
of the
layer
{H^
l
\H^}
is
characterized
by the field of
restoration vectors
r =
r^
1
';
• the
deformation
of the
layer
{H^
2
\
H^
3
'}
is
characterized
by the field of
restoration vectors
r =
(r^
+
r
1
-
2
-*);
• ...
• the
deformation
of the
layer
{H
(
-
l
~
l
\
H^}
is
characterized
by the field of
restoration vectors
8.6.2
Restoration method based
on
DSI
Nine
steps
are
proposed
for
interpolating
the field of
restoration vectors while
taking
into account
the
three constraints
(8.74):
1.
Cover
the
studied domain with
a
regular
orthogonal
rectilinear
3-grid
(see
section (3.4)
and figures
(8.26)
and
(8.28)).
The
size
of the
grid cells
is
chosen
in
such
a way
that
the
problem
to be
solved
can
reasonably
reduce
to the
determination
of the
restoration vector
r(a)
on the set
£7
of
the
vertices
of the
grid.
In
practice
and for
simplicity's sake
• the
(x,
y,
z)
coordinate system
is
identified with
the
(w,
t>,
w)
axis
of the
grids:
• the z
=
w
axis
is
oriented upward.
2.
Build
a
discrete model
,M
3
(f2,
N,
x,
Cx)
where
is
the
current location
of
node
a
€
0,
of the
grid
at
geological time
t
=
1,
while
(3 €
N(a)
if, and
only
if,
|x(a),x(/5)}
is an
edge
of the
grid
not
cut by a
fault. Note
that
the
neighborhoods
TV
(a) so
defined honor
the
discontinuities introduced
by the
faults,
if any
(see
figure
(8.28)).
8.6.
UNFOLDING
A
STACK
OF
LAYERS
Figure
8.26
Computing
the field
of
restoration
vectors
at the
nodes
of
a
regular
linear
3-grid
covering
the
region
to be
restored.
For the
sake
of
clarity,
only
a
vertical
cross
section
is
shown.
3.
Build
a
second discrete model
Ai
3
($7,
N,r,Cr)
where
r(a)
is the
vector
r(x[!](a))
at
location
x^a)
=
x(o;):
Note
that
Q
and the
neighborhood operator
N are the
same
as
those
of
the
model
A4
3
(O,
TV,
x,Cx)
denned
above.
4.
Sample
the top
horizon
H^
at a
series
of
points
(h[i](i)
: i =
1,...,
m}
where
the
restoration vectors
r
(h[i]
(«))
are
assumed
to be
known. Next,
for
each point
h[
1
](i)
located
in a
3-cell
of the
grid,
add
three
fuzzy
Control-Property constraints (see page 186)
to
CY
thus
specifying
that
the
three components
of the
restoration vector need
to be
equal
to the
respective components
of
r(h[i](i))
at
this location.
5.
Add a
"Minimum-Deformation"
DSI
constraint
to
Cr
to
specify
that
the
components
of the
strain
tensor
have
to be as
close
as
possible
to
zero.
However,
to
allow large displacements along
the
faults,
this type
of
constraint must
not be
installed
for
cells
cut by
faults.
6.
Add a
"Mass-Preservation"
DSI
constraint
to Cr to
specify
that
the
continuity
equation (8.73)
has to be
honored
at any
node
a 6
fJ.
7.
If
need
be,
some additional
DSI
constraints
can be
added
to
C
Y
.
8.
Run DSI on the
model
.M
3
(f2,
N,r,C
r
).
9.
The final
step involves using
the
values
{r(a)
: a 6
£1}
to
interpolate
r(x[j])
at any
point
X[i]
to be
restored
in the
studied domain
K
=
"R^
423
424
CHAPTER
8.
ELEMENTS
OF
STRUCTURAL
GEOLOGY
Figure
8.27
Barycentric
interpolation
of
the
restoration
vector
at a
location
x
within
a
regular
linear
3-grid.
For the
sake
of
clarity,
in
this
figure all the
restoration
vectors
are
assumed
to be
parallel
to a
vertical
plane
and
only
a
vertical
cross
section
is
shown.
and to
apply
the
following backward transformation:
The
only theoretical
difficulty
is in
building
the
DSI
constraints corresponding
to
steps
5 and 6;
this
problem
is
addressed
in
sections
(8.6.4),
and
(8.6.5).
Some
practical
implementation
details
(1)
Let
x
=
X[i]
be a
given point
in the
region
72,
=
T^[i]
to be
restored,
and let
C(x[i])
be the
cell containing
x^]
and
belonging
to the
regular linear 3-grid
used
to
compute
the field of
restoration vectors.
As
suggested
in figure
(8.27),
let
(u,v,w)
be the
local normalized coordi-
nates
(see section
(3.4.3))
of
x^]
relative
to the
local
frame
associated with
the
edges
of the
cell
(^(x^]):
According
to
equation (3.16),
the
barycentric interpolation
r(u,v,w)
of the
restoration vector
at
location
x^]
is
defined
as
follows:
In
this
definition
of
r(w,
v,w),
it is
implicitly assumed
that
the
nodes
of the
3-cell
C(x[i])
are
numbered
in
such
a way
that
8.6.
UNFOLDING
A
STACK
OF
LAYERS
Figure
8.28
Vertical
cross
section
showing
how
geometrical
inconsistencies
can
be
avoided
by
removing
the
grid
cells
and the
triangles
of the
horizons
in the
neighborhood
of
faults.
Using
this barycentric interpolation
(8.76)
of r
allows
any
point
x^]
6
*R,\i\
to
be
easily restored according
to the
following
four-step procedure:
1.
determine
the
3-cell
C(x.[i])
containing
X[ij;
2.
compute
the
local
normalized
coordinates
(u,
v,
w)
of
X[i]
relative
to the
local
frame
associated
with
the
edges
of
C(x[i]);
3.
compute
r(u,v,w)
according
to
equation
(8.76);
4.
restore
X[i]
as
follows:
Some
practical
implementation
details
(2)
As
shown
in
figure
(8.28),
two
types
of
geometrical inconsistencies
may
arise
in
the
neighborhood
of
faults:
• A
horizon
may not
meet
the
fault
surface exactly,
thus
generating some
gaps between
its
boundary
and the
fault.
This kind
of
inconsistency
has
no
effect
on the
interpolation
of the field of
restoration vectors
and can
be
completely ignored.
425
426
CHAPTER
8.
ELEMENTS
OF
STRUCTURAL GEOLOGY
• A
horizon
may go a
short
way
beyond
the
fault,
thus generating over-
laps
with
the
fault.
This
kind
of
inconsistency
may
induce
dramatic
errors
when interpolating
the field of
restoration vectors
on
triangles
corresponding
to
such
an
overlap:
—
should
such
a
problem
arise
for the top
horizon
#[1],
then
sampling
points
h[i]
(i)
located
on
inconsistent boundary triangles
will
induce
erro-
neous
Control-Property
constraints
at
step
4 of the
proposed
restoration
algorithm;
—
should
such
a
problem
arise
for any one of the
horizons,
then
the
vertices
of
the
inconsistent
boundary
triangles
will
not be
correctly
restored
at
step
nine
of the
proposed
restoration
algorithm.
A
practical solution
to
help avoid such inconsistencies
is
suggested
in
figure
(8.28)
and
involves proceeding
as
follows
where
all the
horizons
{H*}
to be
restored
are
assumed
to be
modeled
as
triangulated
surfaces
{T(H
1
}}:
• At
step
4 of the
proposed restoration algorithm, remove
all the
sam-
pling
points located
in one or
several rings
of
boundary triangles
in the
neighborhood
of
faults
from
the top
horizon
H
l
=
H^.
• For
each horizon
T(H
1
}
to be
restored, split
the final
step
of the
pro-
posed restoration algorithm into
the
following
substeps:
—
build
a
discrete model
M
3
(£l
l
,
AP,x,
Cx)
where
x
models
the ge-
ometry
of
T(H
1
}
at
current geological time
t = 1
while
N*
models
its
topology;
—
build
a
twin discrete model
A
/
t
3
(O%
7V"
z
,r,Cr)
where
r
models
the
restoration vector
on
T(H
l
);
—
for
each vertex
x(ct
z
)
ofT(H'
1
)
not
located
in the
neighborhood
of
a
fault,
*
determine
the
cell
C(x(a
t
))
of the
regular
rectilinear
3-grid
that
contains
x(a
l
)
and use the
barycentric
interpolation
(8.76)
for
com-
puting
the
restoration
vector
r(a
l
)
at
this
location;
*
add to Cr a
Control-Node
DSI
constraint
specifying
that
r(a
l
)
is
known;
- run the DSI
algorithm
on
.M
3
(fT,
JV%r,C
r
)
for
computing
r(a
l
)
for
all
a
1
6
fT
located
on the
neighborhood
of a
fault;
- for
each vertex
x(a
z
)
of
T(H
l
),
apply
the
restoration transforma-
tion:
It is
relevant
to
note
that
this
method
for
avoiding side
effects
in the
neigh-
borhood
of
faults
does
not
depend
on the
fact
that
the
objects
to be
unfolded
are
surfaces.
As a
consequence,
this
method
can
easily
be
extended
to the
restoration
of any
object described
by a
discrete model
Ai
3
(0
r
,
-/V%x,C
x
).
8.6.
UNFOLDING
A
STACK
OF
LAYERS
8.6.3
Derivatives
of the
restoration vectors
As
already mentioned,
it can be
observed
that
the
restoration vector
r is the
exact opposite
of the
displacement vector
v
introduced
in
section
(8.4):
Taking
into
account
this
remark,
equations
(8.57)
and
(8.73)
clearly
show
that
derivatives
of the field of
restoration vectors directly control both
the
compo-
nents
of the
strain
tensor
and the
continuity equation.
As a
consequence,
the
Minimum-Deformation
and
Mass-Preservation
DSI
constraints mentioned
in
section
(8.6.2)
should make direct
use of
these derivatives.
For
this reason
and
prior
to
studying these constraints,
it is
important
to
show
how the
derivatives
of
r can be
evaluated numerically: this
is
precisely
the
goal
of
this section.
Notations
For the
sake
of
simplicity,
in the
next sections
the
parameterization u(x)
of
the
folded
region
K
is
assumed
to be
identical
to the
vector
x
itself:
Moreover,
depending
on the
context, components
(x^x^x
3
)
and
(r^r^r
3
)
of
the
vectors
x and r
will
be
noted
(x, y, z) and
(r
x
,
r
y
,
r*),
respectively:
Let
us
consider
a
node
a of the
regular 3-grid used
for
estimating
the field of
restoration vectors.
As
suggested
in figure
(8.29),
if,
starting
from
ct,
a
step
is
performed
in the
positive direction
of the x
axis, then
a new
node
a~£
is
generally
reached except, possibly,
if a is
located
on the
boundary
of the
grid.
Conversely,
if we
perform
a
step
in the
negative direction
of the x
axis, then
we
generally reach
a new
node
a~.
Similarly, moving
one
step
in the
negative
and
positive
directions
of the y and z
axes
allows
new
pairs
of
nodes
(a~,
<x+)
and
(a~,ct^~)
to be
reached.
More
generally,
if we do not
follow
the
directions
of
these steps, then
it
we
can say
that
a
3-cell
C =
C(a)
of the
grid
is
fully
identified
by one of its
vertices
a and
three
of its
direct neighbors
(a£,
a^,
a
c
z
}.
As
suggested
in figure
(8.29),
the
four
remaining vertices
of
such
a
3-cell
are
referred
to as
o^
y
,
a£
2
,
oty
Z
,
and
a
c
xyz
,
respectively.
If the
vertices
(a£,o^,o^)
are
reached
from
a
by
doing (positive
or
negative)
steps
(A£,
A^,
A^)
in the
(x,y,z)
directions,
respectively,
then
the
remaining vertices
of C =
C(a)
can be
reached
as
follows
from
a:
•
of
xy
is
reached
by
doing
steps
A£
and
A£;
•
atxz
is
reached
by
doing steps
A^
and
A^;
•
OL
c
yz
is
reached
by
doing
steps
A^
and
A^;
and
427