30

Yang Jinhua
1
, Liu Tao
1
, Tang Genyang
1
, and Hu Tianyue
1
Manuscript received by the Editor July 7, 2008; revised manuscript received January 13, 2009.
*This research is supported by the Natural Science Foundation of China (Grant No. 40574050, 40821062), the National
Basic Research Program of China (Grant No. 2007CB209602), and the Key Research Program of China National Petroleum
Corporation (Grant No. 06A10101).
1. School of Earth and Space Sciences, Peking University, Beijing 100871, China.

Corresponding author (Hu Tianyue, email: tianyue@pku.edu.cn).

APPLIED GEOPHYSICS, Vol.6, No.1 (March 2009), P. 30 - 41 , 18 Figures.
DOI:10.1007/s11770-009-0005-2
Abstract: Seismic modeling is a useful tool for studying the propagation of seismic waves
within complex structures. However, traditional methods of seismic simulation cannot
meet the needs for studying seismic waveelds in the complex geological structures found
in seismic exploration of the mountainous area in Northwestern China. More powerful
techniques of seismic modeling are demanded for this purpose. In this paper, two methods of
nite element-nite difference method (FE-FDM) and arbitrary difference precise integration
(ADPI) for seismic forward modeling have been developed and implemented to understand
the behavior of seismic waves in complex geological subsurface structures and reservoirs.
Two case studies show that the FE-FDM and ADPI techniques are well suited to modeling
seismic wave propagation in complex geology.
Keywords: nite difference, nite element, modeling, arbitrary precise integration
Modeling seismic wave propagation within complex
structures*
Introduction
The complex geology of Northwestern China creates
severe technical difficulties for seismic prospecting,
including: (1) large static corrections, interfering
waveforms from the near-surface, and horizontal velocity
variation due to rapid changes of surface elevation and
lithology (Chang et al., 2002; van Vossen and Trampert,
2007), (2) difficult imaging of complex seismic
wavefields due to velocity inversion in the complex
subsurface area (Zhang and Liu, 2002; Yoon et al.,
2004); (3) velocity modeling problems for prestack depth
migration from low signal-to-noise (S/N) seismic data
(Ekren and Ursin 1999; Li and Peng, 2008); and (4) the
limitations of seismic resolution in reservoir prospecting
for holes and fractures (Coates and Schoenberg
1995; van der Neut et al., 2008). Therefore it is very
important to study wave propagation under complex
structure using seismic simulation. Currently the finite
difference method (FDM) is widely used because of its
computational efciency in seismic modeling, (Carcione
et al., 2002; Wang and Liu, 2007). However, the defects
of FDM reveal themselves in the situations listed above,
such as strong numerical dispersion and the difculty of
dealing with irregular boundaries (Alford et al., 1974;
Lines et al. 1999; Hestholm et al., 2006; Thore 2006; and
Xu et al. 2007). Many improvements have been made in
FDM to solve these problems (Tessmer, 2000; Oliveira,
2003; Ma et al. 2004; Zhang and Chen 2006; and
Moubarak et al. 2007). Two methods of seismic forward
modeling aimed at applying seismic simulation to
complex geological subsurface structures and reservoirs
are developed and implemented in this paper.
The first seismic forward modeling method is the
nite element-nite difference method (FE-FDM), which
31
Yang et al.
combines the advantages of the high computational
efciency of FDM and the excellent adaptation at complex
boundaries of FEM. The FE-FDM technique was initially
introduced by Dong and Yang (2001). Du and Bancroft
(2004) discussed the algorithm for some simple regular
models. The FE-FDM idea is to discretize the seismic
wave equations in the spatial domain with FEM at some
coordinates and with FDM at other coordinates. The aim
of this procedure is to inherit the qualities of both methods
over the entire spatial domain. We developed and applied
it to 2-D acoustic simulation for a complex subsurface
model by introducing a non-linear interpolation function
for FEM and a perfectly matched layer (PML) absorbing
condition to make the algorithm more accurate and stable
(Liu et al., 2008). The second seismic forward modeling
method is the arbitrary difference precise integration
method (ADPI), which employs a FDM scheme in the
spatial domain with an integration scheme in the temporal
domain. This decreases the dispersion error in the time
domain and makes the algorithm more accurate. The
initial theoretical work on ADPI was carried out by Wang
et al. (2004). We developed and applied this method for
simulating seismic wave propagation in a reservoir model
by a velocity-stress scheme with a staggered grid to
improve the efciency and stability (Tang et al. 2007).
Seismic forward modeling methods
Finite element-nite difference method
Consider the 2-D acoustic wave equation:
2 2 2
2 2 2 2
( , , ) ( , , ) 1 ( , , )
0 , in
( , )
( , , )
( , , ) 0 , =0 , when t 0,
u x z t u x z t u x z t
x z c x z t
u x z t
u x z t
t
c c c
+ = O

c c c

= s

c
2 2 2
2 2 2 2
( , , ) ( , , ) 1 ( , , )
0 , in
( , )
( , , )
( , , ) 0 , =0 , when t 0,
u x z t u x z t u x z t
x z c x z t
u x z t
u x z t
t
c c c
+ = O

c c c

= s

c

(1)
where u(x,z,t) denotes the wave function at horizontal
coordinate x, vertical coordinate z (the z axis points
downward), and time t and c(x,z) is the velocity in the
medium.
Discretizing this equation in the z direction by the
Galerkin method, we obtain the weighted equation in
each element e:
2 2 2
2 2 2 2
( , , ) ( , , ) ( , , ) 1
( ) ( , , ) 0,
( , )
e e e
e
e
u x z t u x z t u x z t
w x z t dz
x z c x z t

2 2 2
2 2 2 2
( , , ) ( , , ) ( , , ) 1
( ) ( , , ) 0,
( , )
e e e
e
e
u x z t u x z t u x z t
w x z t dz
x z c x z t

(2)
where u
e
(x,z,t) is the wave function and w
e
(x,z,t) is a
weighting function within an element.
Assume that we have

1 1
( , ) ( ), ( , ) ( )
e e
N N
e ei ei e ei ei
i i
u u x t H z w v x t H z



, (3)
where N
e
is the number of grid points in each element e,
H
e
(z) is an interpolation function, and u
ei
and v
ei
are the
values of u and w at discrete grid points.
Substituting equation (3) into equation (2), we get the
matrix partial differential equations (PDEs) in each element

2 2
2 2
,
t x



e e
e e e e
u u
M K u A (4)
where

2
1
, , .
( , )
e e e
dz dz dz
c x z z z




T
T T e e
e e e e e e e
H H
M H H K A H H

Then, the matrix PDEs can be expanded in the entire
domain in equation (5):

2 2
2 2
,
t x



u u
M K u A
(5)
with
1 1 1
, ,
N N N

e e e
M M K K A A and N is the number
of discrete grid nodes in the z direction.
If we employ 3 nodes in an element, then the
interpolation function turns out to be nonlinear:

2 2
2
, 1- , ,
2 2

( +
=
(

e
H (6)
where
2
3 1
z z
z z

and z
1
, z
2
, and z
3
are the values of grid
nodes in the z direction in each element.
Assuming the velocity in each element is constant, we
can compute the matrices M, K, and A as:

2
4 2 -1 0 0 0 0 0
2 16 2 0 0 0 0 0
-1 2 8 2 -1 0 0 0
0 0 2 16 2 0 0 0
=
0 0 -1 2 8 2 -1 0 15


h
c
M

0 -1 2 4















7 -8 1 0 0 0 0 0
-8 16 -8 0 0 0 0 0
1 -8 14 -8 1 0 0 0
0 0 -8 16 -8 0 0 0 1
0 0 1 -8 14 -8 1 0 6


h
K

0 1 -8 7

4 2 -1 0 0 0 0 0
2 16 2 0 0 0 0 0
-1 2 8 2 -1 0 0 0
0 0 2 16 2 0 0 0
=
0 0 -1 2 8 2 -1 0 15


h
A

.

0 -1 2 4















(7)
32
Modeling seismic wave propagation
We notice that it is unnecessary to constrain the space
step in the z direction to be a constant, which allows
us to deal with irregular grids. Furthermore, a lumped
matrix (Marfurt, 1984) can be employed in place of M to
make the computation convenient

lumped 2
5 0 0 0 0 0 0 0
0 20 0 0 0 0 0 0
0 0 10 0 0 0 0 0
0 0 0 20 0 0 0 0
=
0 0 0 0 10 0 0 0 15


h
c
M

.

0 0 0 5















. (8)
Then, we substitute the matrices M
lumped
, K, and A into
equation (5) to obtain:
2
2
1 1 2
2 2 2
1 1 2
2 2 2
2
2
2 1 2
[ ]
( ) ( [ ] [ ] 2 [ ] )
[ ] [ ] [ ]
( ) (0.1 0.8 0.1 ), if is an even number
[ ]
( ) (0.25 [ ] 2 [ ] 3.5 [ ] 2 [
n
j n n n
j j j
n n n
j j j
n
j n n n
j j j
u i
c
u i u i u i
t h
u i u i u i
c j
x x x
u i
c
u i u i u i u i
t h

1 2
2 2 2 2 2
2 1 1 2 2
2 2 2 2 2
] 0.25 [ ] )
[ ] [ ] [ ] [ ] [ ]
( ) ( 0.1 0.2 0.8 0.2 0.1 ),

n n
j j
n n n n n
j j j j j
u i
u i u i u i u i u i
c
x x x x x

if is an odd number,

j
(9)
if j is an even number

2
2
1 1 2
2 2 2
1 1 2
2 2 2
2
2
2 1 2
[ ]
( ) ( [ ] [ ] 2 [ ] )
[ ] [ ] [ ]
( ) (0.1 0.8 0.1 ), if is an even number
[ ]
( ) (0.25 [ ] 2 [ ] 3.5 [ ] 2 [
n
j n n n
j j j
n n n
j j j
n
j n n n
j j j
u i
c
u i u i u i
t h
u i u i u i
c j
x x x
u i
c
u i u i u i u i
t h

1 2
2 2 2 2 2
2 1 1 2 2
2 2 2 2 2
] 0.25 [ ] )
[ ] [ ] [ ] [ ] [ ]
( ) ( 0.1 0.2 0.8 0.2 0.1 ),

n n
j j
n n n n n
j j j j j
u i
u i u i u i u i u i
c
x x x x x

if is an odd number,

j
2
2
1 1 2
2 2 2
1 1 2
2 2 2
2
2
2 1 2
[ ]
( ) ( [ ] [ ] 2 [ ] )
[ ] [ ] [ ]
( ) (0.1 0.8 0.1 ), if is an even number
[ ]
( ) (0.25 [ ] 2 [ ] 3.5 [ ] 2 [
n
j n n n
j j j
n n n
j j j
n
j n n n
j j j
u i
c
u i u i u i
t h
u i u i u i
c j
x x x
u i
c
u i u i u i u i
t h

1 2
2 2 2 2 2
2 1 1 2 2
2 2 2 2 2
] 0.25 [ ] )
[ ] [ ] [ ] [ ] [ ]
( ) ( 0.1 0.2 0.8 0.2 0.1 ),

n n
j j
n n n n n
j j j j j
u i
u i u i u i u i u i
c
x x x x x

if is an odd number,

j
2
2
1 1 2
2 2 2
1 1 2
2 2 2
2
2
2 1 2
[ ]
( ) ( [ ] [ ] 2 [ ] )
[ ] [ ] [ ]
( ) (0.1 0.8 0.1 ), if is an even number
[ ]
( ) (0.25 [ ] 2 [ ] 3.5 [ ] 2 [
n
j n n n
j j j
n n n
j j j
n
j n n n
j j j
u i
c
u i u i u i
t h
u i u i u i
c j
x x x
u i
c
u i u i u i u i
t h

1 2
2 2 2 2 2
2 1 1 2 2
2 2 2 2 2
] 0.25 [ ] )
[ ] [ ] [ ] [ ] [ ]
( ) ( 0.1 0.2 0.8 0.2 0.1 ),

n n
j j
n n n n n
j j j j j
u i
u i u i u i u i u i
c
x x x x x

if is an odd number,

j
2
2
1 1 2
2 2 2
1 1 2
2 2 2
2
2
2 1 2
[ ]
( ) ( [ ] [ ] 2 [ ] )
[ ] [ ] [ ]
( ) (0.1 0.8 0.1 ), if is an even number
[ ]
( ) (0.25 [ ] 2 [ ] 3.5 [ ] 2 [
n
j n n n
j j j
n n n
j j j
n
j n n n
j j j
u i
c
u i u i u i
t h
u i u i u i
c j
x x x
u i
c
u i u i u i u i
t h

1 2
2 2 2 2 2
2 1 1 2 2
2 2 2 2 2
] 0.25 [ ] )
[ ] [ ] [ ] [ ] [ ]
( ) ( 0.1 0.2 0.8 0.2 0.1 ),

n n
j j
n n n n n
j j j j j
u i
u i u i u i u i u i
c
x x x x x

if is an odd number,

j
if j is an odd number,
where is the time step, i and j denote the position
of the grid node in space, n is the time factor, h
denotes the grid step in the z direction and l for the x
direction.
For the 2-D wave equation, after getting the matrix
PDEs, we employ FDM to solve the equation in the x
and time domains in equations (10) and (11).

2
1 1
2 2
1
( [ ] 2 [ ] [ ] ),
n n n
j j j
u
u i u i u i
t

(10)


2 2
0 2 2
1
1
( ( [ ] [ ] ) [ ] ),
n n n
m j j j
u
c u i m u i m c u i
x l

(11)
where

0 1 2
2.5, =1.333, = - 0.0833 c c c .
Here we use a second-order explicit difference
algorithm in the time domain and a fourth-order
difference algorithm in the x domain. If we use
three nodes in each element for the interpolation
f unc t i on a nd a f our t h- or de r e xpl i c i t f i ni t e
difference scheme in the spatial domain, we obtain
t he confi gurat i on of gri ds shown i n Fi gure 1,
which compares the FE-FDM grid with the FDM
grid.
(a) (b)
Fig.1 Grid congurations in the spatial domain. (a) FDM using nine grid points, (b) FE-FDM using 25 grid points.
x
z
x
z
2. Arbitrary difference precise integration
method (ADPI)
Consider the following system of wave equations (see
Tang et al. 2007):

(0)
(0) ,

0
0
MU CU KU F
U U
U U


(12)
33
Yang et al.
where M is the mass matrix, C is the damping matrix, K
is the stiffness matrix, and F is an equivalent load.
Introducing P=MU+CU/2, we obtain:

, V HV f

(13)
where, ,

U
V
P

0
f
F
, ,

U
V
P

0
f
F
, and
/
.
1
/ 2
4

1 1
1 1
M C 2 M
H
CM C K CM
.
By integration in time domain, the solution of this
equation is:

( )
0
( ) ( ) .
t
t t
t e e d

H H
0
V V f (14)
If
1 n n n
t n t t t t

,
1 n n n
t n t t t t

, the solution at tn+1 is:

1
1
( )
1
( ) ( ) ( ) ( ) ,
n
n
n
t
t
n n
t
t t t e d

H
V T V f (15)
where ( ) .
t
t e


H
T If we consider the value of f is a constant
one at the interval [t
n
,t
n
+1], then we can get the explicit
scheme to solve equation (15): V
n+1
=TV
n
+(T-1)H
-1
f
n
. After
getting the value of V, then we can obtain the value of U.
For the 2-D acoustic equation implemented in the
velocity-stress scheme, we obtain:

1
1
2
2
( )
( )
x y
x
x x
x y
z
z z
x x
x x
z z
z z
u u
v
q v
t x
u u
v
q v
t z
u v
q u v
t x
u v
q u v
t z

c + c
+ =

c c

c +
c
+ =

c c

c c

+ =

c c

c c

+ =

c c
(16)

,
where v
x
,v
z
,u
x
, and u
z
denote the velocity and stress
parameters in the x and z directions, q
x
and q
z
denote the
damped coefficients in the x and z directions, is the
density, and v is the velocity of the medium.
To solve equation (16), we use finite-difference
method to deal with the partial differential equation in
space domain first, and then obtain the value through
integration in time domain. Take vx as an example: after
employing FDM in x direction, we have:


1/ 2, 1/ 2,
1
1
( )
N
n n x
l xi l j xi l j
l
v
a v v
x x

, (17)
where x is the grid step in x direction, i and j denote
the position of the grid node in space, a
l
is the FDM
coefcient in the spatial domain and 2N represents the
order of scheme that we employ for FDM. Here, the
nite difference scheme is based on staggered grids and,
if the fourth-order explicit algorithm is chosen (N=2),
the coefcients should be: a
1
=1.125, a
2
=-0.041667.
Then substitute equation (17) to (16), the equation of
v
x
turns to be:

2
1/ 2, 1/ 2,
1
( )
N
n n x
x x l xi l j xi l j
l
u
q u v a v v
t

. (18)
Compared with equation (13), here

2
1/ 2, 1/ 2,
1
, , f ( )
N
n n
x x l xi l j xi l j
l
u q and v a v v

V H
.
After integration in the time domain with the explicit
scheme, the parameter u
x
can be obtained as:
2
, , 1/ 2 1/ 2
, , 1/ 2, 1/ 2,
1
(1 ) ( ),
xi xi
N
i j i j q t q t n n n n
xi j xi j l xi l j xi l j
l xi
v
u u e e a v v
x q

2
, , 1/ 2 1/ 2
, , 1/ 2, 1/ 2,
1
(1 ) ( ),
xi xi
N
i j i j q t q t n n n n
xi j xi j l xi l j xi l j
l xi
v
u u e e a v v
x q

. (19)
Other parameters such as u
z
, v
x
and v
z
can be easily
obtained in the same way.
Case studies
Two case studies will study the application of FE-
FDM and ADPI to wave propagation studies.
Case 1: A model with complex structure (FE-
FDM)
Figure 2 shows a subsurface model with folds beneath a
weathered surface layer and above a limestone layer and a
sandstone layer based on a currently active seismic survey
in an area typical of northwestern China. The width of the
model is 20 km and the depth is 3.6 km. In the forward
modeling by FE-FDM, 128 shot gathers were calculated,
where each shot gather has 60 receivers to either side in a
split-spread. The receiver interval is 40 m with a minimum
offset of 40 m. A 25 Hz Ricker wavelet is chosen as
the source and the first shot is indicated in the figure.
An absorbing boundary condition (Collino and Tsogka,
2001; Komatictsch and Tromp, 2003) is employed for the
weathering layered at the surface.
Figure 3 depicts the gathers of the rst shot, recorded
at the surface and computed using both the FE-FDM
and FDM methods. From this figure we see that the
amplitude beneath the third layer damps rapidly as
a result of the weathering stratum at the surface.
Reflections, multiples, and diffractions due to the
variation of lithology are clearly visible in the shot
gathers and the results of both methods agree well. The
50th trace of the shot gather is shown in Figure 4, which
also shows good agreement between the FE-FDM and
FDM results.
34
Modeling seismic wave propagation
Fig. 2 Model with an irregular interface, 20 km wide, and 3.6 km depth. It is a typical model from north-western China, where a
weathering layer, folds, limestone layer and sandstone layer can be found.
Fig. 4 Seismogram from the 50th geophone. The red curve is the FE-FDM result and the black curve is the FDM result. The good
match between the wavelets also assures the efciency of FE-FDM.
(a) (b)
Fig. 3 Shot gathers of the rst shot. (a) Shot gather result from FE-FDM; (b) Shot gather result from FDM. The good agreement
between the FE-FDM and FDM results indicate the efciency of FE-FDM.
10
0
-8
0 200 800 1200 1600 1800
Time (ms)
FE-FDM
FDM
Figure 5 is the post-stack depth migrated section from
using the Kirchhoff algorithm for the gathers obtained
from FE-FDM modeling. The velocity was obtained
from the velocity model. It provides a simple test of
2000 6000 10000 14000 18000
100
1000
2000
3000
20000
1500 m/s
1800 m/s
3000 m/s
3500 m/s
3700 m/s
4500 m/s
3600 m/s
4200 m/s
Limestone
Sand stone
1800 m/s
2000 m/s
2500 m/s
3000 m/s
3400 m/s
First shot Weathering layer
Width (m)
D
e
p
t
h

(
m
)
100 125 150162 187 217
100
400
800
1200
1600
1900
T
i
m
e

(
m
s
)
CDP
Reection
Diffraction
Multiple
100 125 150162 187 217
100
400
800
1200
1600
1900
T
i
m
e

(
m
s
)
CDP
Reection
Diffraction
Multiple
35
Yang et al.
the performance of this technique. In comparison with
Figure 2, the migrated section reveals the subsurface
structure of the model clearly, even for the complex
interface indicated on the gure which could identify a
reservoir. In this gure, some small features indicated by
the ellipse can be clearly identied and which conrms
the good match between the FE-FDM profile and the
model. Also, the multiple generated between the second
and the third layers can be identied. So, the accuracy of
FE-FDM is sufcient for the seismic simulation of our
complex area.
Fig. 5 Depth migrated section produced from the synthetic shot gathers. Here we use the Kirchhoff algorithm to implement post-
stack depth migration using the velocities from the velocity model. The good match with the original model demonstrates the
accuracy of FE-FDM.
Now, we study the application of the FE-FDM
technique to an exploration problem. The profile from
the stack of the gathers calculated by FE-FDM is
compared with real processed data to test the original
geological model. Figure 6 is a section of real seismic
data from northwestern China and Figure 7 is the
synthetic stack section obtained from the gathers
computed by FE-FDM from the model in Figure 2. The
events (particularly around 2.0 s) in these two gures do
not match each other very well, which indicates that the
model of Figure 2 cannot describe the structure in this
area correctly.
Fig. 6 A real seismic stacked section obtained from Tarim Oil Field, CNPC.
140 280 420 560 840 700
400
1200
2000
2800
3600
D
e
p
t
h

(
m
)
CDP
160 340 520 700 880
CDP
400
800
1200
1600
2000
2400
T
i
m
e

(
m
s
)
36
Modeling seismic wave propagation
Fig. 7 The raw synthetic stacked section from the model in Figure 2. Compared with the real seismic data, we can see the differences,
particularly around 2000 ms, which indicates that the model we used in Figure 2 could not describe the structure underground well.
We modified the model iteratively to approach a more
reasonable model shown in Figure 8, where the layer
identified in Figure 2 has been changed to be relatively
smooth and concave. The synthetic stack from this revised
model is shown in Figure 9 and approximates the real data
section much better. The application of seismic simulation by
FE-FDM helps us to evaluate the model provided from seismic
interpretation and avoid unnecessary investment in drilling.
Fig. 8 The modied model. We modied the model iteratively to improve the match of the real and synthetic data.
Fig. 9 The synthetic stack section from the modied model. The result is closer to the real data compared with the result in Figure 7. This means the
new model we presented in Figure 8 reveals the real situation better. This procedure helps us to evaluate the model from seismic interpretation.
160 340 520 880 700
400
1200
2000
2400
T
i
m
e

(
m
s
)
CDP
800
1600
0
160 340 520 790 700
400
1200
2000
2400
T
i
m
e

(
m
s
)
CDP
800
1600
0
2000 6000 10000 14000 18000
200
1000
2000
3000
1500 m/s
1800 m/s
3000 m/s
3500 m/s
3700 m/s
4500 m/s
3600 m/s
4200 m/s
1800 m/s
2000 m/s
2500 m/s
3000 m/s
3400 m/s
3500 m/s
D
e
p
t
h

(
m
)
Width (m)
37
Yang et al.
Case 2: A model with small-scale caverns
(ADPI)
We apply the ADPI method to calculate the acoustic
wave propagation in a carbonate reservoir model
featuring small-scale karst caverns and validate the
method through comparison with a physical modeling
result. The model is shown in Figure 10. The larger
caverns size is 80 m in length and 40 m in depth to the
right in the gure and the smaller cavern is 40 m in both
length and depth to the left. We calculate 60 shots to
acquire the seismic data with a shot interval of 50 m and
100 receivers in an off-end spread. The source signature
is a 25 Hz Ricker wavelet. The maximum offset is 1090
m and minimum is 100 m.
Fig. 10 A model characterized by two small-scale caverns. The left cavern is 40 m in both length and depth and the larger cavern
on the right is 80 m long and 40m deep.
To study the effect of the caverns, the computed
results of the model with caverns is compared with a
model without caverns. Figures 11 and 12 shows the
1st, 21st, and 51st shot gathers from the two models
after computation with ADPI. The stacked sections are
compared in Figures 13 and 14. From comparison of the
shot gathers and stacks, the strong diffractions generated
by the caverns can be seen.
Fig. 11 Shot gathers computed from the model of Figure 10. The strong diffractions from the caves can be clearly seen..
3200
51
0 25 50 75 100
21
Shot Number
0 25 50 75 100
1
T
i
m
e

(
m
s
)
400
1200
1600
2400
0 25 50 75 100 Receiver
0 100 200 300 400 500 600 700
0
1000
2000
3000
4000
5000
6000
3000 m/s
2500 kg/m
3
2200 m/s
2000 kg/m
3
3000 m/s
2500 kg/m
3
4400 m/s
3000 kg/m
3
5000 m/s
3500 kg/m
3
2000 m/s
1970 kg/m
3
CDP
D
e
p
t
h

(
m
)
38
Modeling seismic wave propagation
Fig. 13 The stacked section of the model in Figure 10. The velocity used for stacking is derived from the velocity model. There
are strong diffractions generated by the caves.
Fig. 12 Shot gathers computed from the model without caves.
Fig. 14 The stacked section of the model without caves.
51
0 25 50 75 100
21
Shot Number
0 25 50 75 100
1
T
i
m
e

(
m
s
)
400
1200
1600
2400
3200
0 25 50 75 100 Receiver
500
1500
2500
3500
100
CDP
190 280 370 460 550 640
1000
2000
3000
T
i
m
e

(
m
s
)
500
1500
2500
3500
100
CDP
190 280 370 460 550 640
1000
2000
3000
T
i
m
e

(
m
s
)
39
Yang et al.
Figures 15 and 16 are the post-stack depth migration
results (using the Kirchhoff algorithm and velocities
from the original model) for both models which reveals
the original structure well. Figure 17 is the post-stack
time migration result. After time and depth migration,
the caverns still generate strong and short reflections
which are caused by the multiple reflections due to
the caves. Figure 18 is the post-stack time migrated
section from physical modeling (Li and Zhao 2006). The
match between the ADPI and physical modeling results
demonstrates that the ADPI technique can simulate
seismic waveelds for the model of a cavern reservoir.
Through the simulation, the behavior of waves in the
cavern model can be clearly illustrated, such as the
strong diffractions in stacked results and the multiple
reflections in migrated results, which are the typical
signs for the caves.
Fig. 15 Post-stack depth migration of the model in Figure 10. A Kirchhoff algorithm was used for migration and the migration
velocities came from the velocity model. From this result, the multiple reections generated by the larger cave are clear to see.
Fig. 16 The post-stack depth migration of the model without caves.
Conclusions
The finite difference method is widely employed
in seismic forward modeling due to its efficiency
in computation but more powerful techniques are
demanded as the subsurface structure becomes more
complex. In order to minimize dispersion error and adapt
the computations to irregular boundaries, the FE-FDM
method discretizes the spatial domain separately with
500
1500
2500
3500
100
CDP
190 280 370 460 550 640
D
e
p
t
h

(
m
)
4500
5500
500
1500
2500
3500
100
CDP
190 280 370 460 550 640
D
e
p
t
h

(
m
)
4500
5500
40
Modeling seismic wave propagation
Fig. 17 The post-stack time migrated section from the model in Figure 10, also using the Kirchhoff algorithm. The multiple cave
reections can be clearly seen.
FEM and FDM to solve the irregular boundary problem,
while keeping the computational efficiency. For 2-D
models, employing FEM on the z coordinate makes
it possible to discretize freely in that direction, which
allows thin layers to be easily modeled. The method can
be extended to 3-D by discretizing in x and z coordinates
with FEM to simulate the topography and with FDM in
the y coordinate to keep the efciency of computation.
An alternative method called ADPI (arbitrary difference
precise integration) increases the accuracy by replacing
the FDM discretization with integration. It also can
be extended to 3-D. Both methods improve FDM in
different ways and involve extra computing costs but
parallel processing allows us to solve this problem. Two
case studies are given in this paper to test these two
methods. The results confirmed the potential of both
methods. More work is underway to extend the methods
to simulate elastic waves in 3-D in order to simulate real
seismic experiments more closely.
Acknowledgement
We thanks for that the real seismic data and the
velocity model with complex structure for northwestern
China were provided by Tarim Oileld, PetroChina.
References
Alford, R. M., Kelly, K. R., and Boore, D. M., 1974,
Accuracy of finite-difference modeling of the acoustic
wave equation: Geophysics, 39(6), 834 842.
Carcione, J. M., Herman, G. C., and ten Kroode, A. P. E.,
2002, Seismic modeling: Geophysics, 67(4), 1304
1325.
Chang, X., Liu, Y., Wang, H., and Li, F., 2002, 3-D
tomographic static correction: Geophysics, 67(4), 1275
1285.
Coates, R. T., and Schoenberg, M., 1995, Finite-difference
modeling of faults and fractures: Geophysics, 60(5),
1514 1526.
Collino, F., and Tsogka, C., 2001, Application of the
perfectly matched absorbing layer model to the linear
elastodynamic problem in anisotropic heterogeneous
Fig. 18 The post-stack migrated section produced from the
physical model (Li and Zhao, 2006).
1500
2500
3500
CDP
100 200 300 400 500 600
T
i
m
e

(
m
s
)
170
1500
2000
2500
3000
3500
1000
T
i
m
e

(
m
s
)
CDP
280 370 460 550 640
41
Yang et al.
media: Geophysics, 66(1), 294 307.
Dong, Y., and Yang, H. Z., 2001, Solution of 2-d wave
reverse-time propagation problem by the nite element-
nite difference method: Theoretical and Computational
Acoustics:,World Scientic Press, Singapore.
Du, X., and Bancroft, J. C., 2004, 2-D wave equation
modeling and migration by a new finite difference
scheme based on the Galerkin method: 74th Ann.
Int ernat . Mt g. , Soc. Expl . Geophys. , Expanded
Abstracts., 1107 1110.
Ekren, B. O., and Ursin, B., 1999, True-amplitude
frequency-wave number constant-offset migration:
Geophysics, 61(3), 915 924.
Hestholm, S., Moran, M., Ketcham, S., Anderson, T.,
Dillen, M., and McMechan, G., 2006, Effects of free-
surface topography on moving-seismic-source modeling:
Geophysics, 71(6), T159 T166.
Komatictsch, D., and Tromp, J., 2003, A perfectly matched
layer absorbing boundary condition for the second-order
seismic wave equation: Geophys. J. Int., 124, 146-153.
Li, G. F., and Peng, S. P., 2008, Static corrections for
low S/N ratio converted-wave seismic data: Applied
Geophysics, 5(1), 44 49.
Li, J. F., and Zhao, Q., 2006, Figures of seismic physical
modeling for oil and gas exploration: China Petroleum
Industry Press, Beijing.
Lines, L. R., Slawinski, R., and Bording, R. P., 1999, A
recipe for stability of finite-difference wave-equation
computations: Geophysics, 64(3), 967 979.
Liu, T., Hu, T. Y., and Yang, J. H., 2008, Finite element-
nite difference method in 2-D seismic modeling: 70th
EAGE Conference, Extended Abstracts, Rome, P054.
Ma, S., Archuleta, R. J., and Liu, P., 2004, Hybrid modeling
of elastic P-SV wave motion: A combined nite-element
and staggered-grid finite-difference approach: Bulletin
of the Seismological Society of America, 94(4), 1557
1563.
Marfurt, K. J., 1984, Accuracy of finite-difference and
finite-element modeling of the scalar and elastic wave
equations: Geophysics, 49(5), 533 549.
Moubarak, H., Bancroft, J., Lawton, D., Isaac, H.,
Mewhort, L., Emery, D., and Scott, B., 2007, A modeling
study for imaging in structurally complex media: case
history: 77th Ann. Internat. Mtg., Soc. Expl. Geophys.,
Expanded Abstracts, 2002 2005
Oliveira, S. A. M., 2003, A fourth-order finite-difference
method for the acoustic wave equation on irregular grids:
Geophysics, 68(2), 672 676.
Tang, G. Y., Hu, T. Y. and Yang, J. H., 2007, Applications
of a precise integration method in forward seismic
modeling: 77th Ann. Internat. Mtg., Soc. Expl. Geophys.,
Expanded Abstracts, 2130 2133.
Tessmer, E., 2000, Seismic nite-difference modeling with
spatially varying time steps: Geophysics, 65(4), 1290
1293.
Thore, P., 2006, Accuracy and limitation in seismic
modeling of reservoir: 76th Ann. Internat. Mtg., Soc.
Expl. Geophys., Expanded Abstracts, 1674-1677
van der Neut, J., Sen, M. K., and Wapenaar, K., 2008,
Sei smi c refl ect i on coeffi ci ent s of faul t s at l ow
frequencies: a model study: Geophysical Prospecting,
56(3), 287 292.
van Vossen, R., and Trampert, J., 2007, Full-waveform
static correction using blind channel identification:
Geophysics, 72(4), U55 U66.
Wang, X. C., and Liu, X. W., 2007, 3-D acoustic wave
equation forward modeling with topography: Applied
Geophysics, 4(1), 8 15.
Wang, R. Q., Jia, X. F., and Hu, T. Y., 2004, The precise
finite difference method for seismic modeling: Applied
Geophysics, 1(2), 69 74.
Xu, Y. X., Xia, J. H., and Miller, R. D., 2007, Numerical
investigation of implementation of air-earth boundary by
acoustic-elastic boundary approach: Geophysics, 72(5),
SM147 SM153.
Yoon, K, Marfurt, K. J., and Starr, W., 2004. Challenges in
reverse-time migration: 77th Ann. Internat. Mtg., Soc.
Expl. Geophys., Expanded Abstracts, 1057 1060.
Zhang, J., and Liu, T. L., 2002, Elastic wave modeling
in 3D heterogeneous media: Geophysical Journal
International, 150(3), 780 799.
Zhang, W., and Chen, X., 2006, Traction image method
for irregular free surface boundaries in nite difference
sei smi c wave si mul at i on: Geophysi cal Journal
International, 167, 337 353.
Yang Jinhua: See biography and photo in the APPLIED
GEOPHYSICS September 2007 issue, P. 206