**4. Free boundary conditions in the Cartesian and curvilinear coordinate systems**

Three-Dimensional Wavefield Simulation in

Heterogeneous Transversely Isotropic Medium with Irregular Free Surface 111

; ;

q qq 0 i,j,k + i,j,k - i,j,k

<sup>1</sup> D u = (D u + D u ) <sup>2</sup>

<sup>1</sup> D u = (D u + D u ) <sup>2</sup>

<sup>1</sup> D u = (D u + D u ) <sup>2</sup>

r rr 0 i,j,k + i,j,k - i,j,k

s ss 0 i,j,k + i,j,k - i,j,k

(14)

(15)

**5. A discretization scheme on curvilinear grid** 

**Figure 2.** Grids distributions in curvilinear coordinate.

components of the wave field are given by

q i+1,j,k i,j,k

u -u

r i,j+1,k i,j,k

s i,j,k+1 i,j,k

q

h u -u

r

q q - i,j,k + i -1,j,k r r - i,j,k + i,j-1,k s s - i,j,k + i,j,k-1

D u =D u D u =D u D u =D u

h u -u

s

h

discretized according to the following equations

and the derivation operators as

+ i,j,k

Du =

+ i,j,k

Du =

+ i,j,k

Du =

i q j r k s q r s

where l, w, h are the length of the rectangular solid in q-, r- and s-directions, respectively; qrs h ,h ,h >0 define the grid size in q-, r- and s-directions, respectively. The three

i,j,k i,j,k i,j,k ijk ijk ijk [u (t), v (t),w (t)] = [u(q ,r ,s , t), v(q ,r ,s , t),w(q ,r ,s , t)]

The right hand sides of eqs. (7)-(9) contain spatial derivatives of nine basic types, which are

q q r r s s

h = l (N - 1) h = w (N - 1) h = h (N - 1)

i = 1,...,N j = 1,...,N k = 1,...,N

q = (i - 1)h r =(j - 1)h s = (k - 1)h

To approximate (7)-(9) we discretize the rectangular solid (Figure 2)

At the free surface, the boundary conditions in the Cartesian coordinates are given by:

$$\begin{bmatrix} \mathbf{c}\_{11} \frac{\partial \mathbf{u}}{\partial \mathbf{x}} + \mathbf{c}\_{12} \frac{\partial \mathbf{v}}{\partial \mathbf{y}} + \mathbf{c}\_{13} \frac{\partial \mathbf{w}}{\partial \mathbf{z}} & \mathbf{c}\_{66} \frac{\partial \mathbf{u}}{\partial \mathbf{y}} + \mathbf{c}\_{66} \frac{\partial \mathbf{v}}{\partial \mathbf{x}} & \mathbf{c}\_{44} \frac{\partial \mathbf{u}}{\partial \mathbf{z}} + \mathbf{c}\_{44} \frac{\partial \mathbf{w}}{\partial \mathbf{x}}\\\ \mathbf{c}\_{66} \frac{\partial \mathbf{v}}{\partial \mathbf{x}} + \mathbf{c}\_{66} \frac{\partial \mathbf{u}}{\partial \mathbf{y}} & \mathbf{c}\_{11} \frac{\partial \mathbf{v}}{\partial \mathbf{y}} + \mathbf{c}\_{12} \frac{\partial \mathbf{u}}{\partial \mathbf{x}} + \mathbf{c}\_{13} \frac{\partial \mathbf{w}}{\partial \mathbf{z}} & \mathbf{c}\_{44} \frac{\partial \mathbf{v}}{\partial \mathbf{z}} + \mathbf{c}\_{44} \frac{\partial \mathbf{w}}{\partial \mathbf{y}}\\\ \mathbf{c}\_{44} \frac{\partial \mathbf{w}}{\partial \mathbf{x}} + \mathbf{c}\_{44} \frac{\partial \mathbf{u}}{\partial \mathbf{z}} & \mathbf{c}\_{44} \frac{\partial \mathbf{w}}{\partial \mathbf{y}} + \mathbf{c}\_{44} \frac{\partial \mathbf{v}}{\partial \mathbf{z}} & \mathbf{c}\_{33} \frac{\partial \mathbf{w}}{\partial \mathbf{x}} + \mathbf{c}\_{13} \frac{\partial \mathbf{u}}{\partial \mathbf{y}} + \mathbf{c}\_{13} \frac{\partial \mathbf{v}}{\partial \mathbf{y}}\end{bmatrix$$

Here T xyz n ,n ,n is the inward normal of the free surface. Using relationships (2), the above boundary conditions in the curvilinear coordinates can be re-written as

 x 11 x q x r x s 12 y q y r y s 13 z q z r z s y 66 y q y r y s 66 x q x r x s z 44 z q z r z s 44 x q x r x s s c (q u + r u + s u ) + c (q v + r v + s v ) + c (q w + r w + s w ) +s c (q u + r u + s u ) + c (q v + r v + s v ) +s c (q u + r u + s u ) + c (q w + r w + s w ) = 0 (11)

 x 66 x q x r x s 66 y q y r y s y 11 y q y r y s 12 x q x r x s 13 z q z r z s z 44 z q z r z s 44 y q y r y s s c (q v + r v + s v ) + c (q u + r u + s u ) +s c (q v + r v + s v ) + c (q u + r u + s u ) + c (q w + r w + s w ) +s c (q v + r v + s v ) + c (q w + r w + s w ) = 0, (12)

 x 44 x q x r x s 44 z q z r z s y 44 y q y r y s 44 z q z r z s z 33 z q z r z s 13 x q x r x s 13 y q y r y s s c (q w + r w + s w ) + c (q u + r u + s u ) +s c (q w + r w + s w ) + c (q v + r v + s v ) +s c (q w + r w + s w ) + c (q u + r u + s u ) + c (q v + r v + s v ) = 0 (13)

Note that here the normal is represented by the normalized metric (evaluated along the free surface)

$$\mathbf{\dot{s}}\_{\mathbf{x}} = \frac{\mathbf{s}\_{\mathbf{x}}}{\sqrt{\mathbf{s}\_{\mathbf{x}}^2 + \mathbf{s}\_{\mathbf{y}}^2 + \mathbf{s}\_{\mathbf{z}}^2}}, \quad \mathbf{\dot{s}}\_{\mathbf{y}} = \frac{\mathbf{s}\_{\mathbf{y}}}{\sqrt{\mathbf{s}\_{\mathbf{x}}^2 + \mathbf{s}\_{\mathbf{y}}^2 + \mathbf{s}\_{\mathbf{z}}^2}}, \quad \mathbf{\dot{s}}\_{\mathbf{z}} = \frac{\mathbf{s}\_{\mathbf{z}}}{\sqrt{\mathbf{s}\_{\mathbf{x}}^2 + \mathbf{s}\_{\mathbf{y}}^2 + \mathbf{s}\_{\mathbf{z}}^2}}.$$

#### **5. A discretization scheme on curvilinear grid**

110 Earthquake Engineering

Here 










surface)

T

**systems** 

**4. Free boundary conditions in the Cartesian and curvilinear coordinate** 

At the free surface, the boundary conditions in the Cartesian coordinates are given by:

uvw uv uw

xy yx z z y

11 12 13 66 66 44 44

c +c +c c +c c +c

above boundary conditions in the curvilinear coordinates can be re-written as

 

y 66 y q y r y s 66 x q x r x s

+s c (q u + r u + s u ) + c (q v + r v + s v )

z 44 z q z r z s 44 x q x r x s

 

x 66 x q x r x s 66 y q y r y s

s c (q v + r v + s v ) + c (q u + r u + s u )

 

+s c (q v + r v + s v ) + c (q w + r w + s w ) = 0,

,

s =


s =

s

x y z

s +s +s

z 44 z q z r z s 44 y q y r y s

 

y 44 y q y r y s 44 z q z r z s

+s c (q w + r w + s w ) + c (q v + r v + s v )

x 44 x q x r x s 44 z q z r z s

s c (q w + r w + s w ) + c (q u + r u + s u )

+s c (q u + r u + s u ) + c (q w + r w + s w ) = 0

44 44 44 44 33 13 13

 

xyz n ,n ,n is the inward normal of the free surface. Using relationships (2), the

66 66 11 12 13 44 44 y

wu wv wuv c +c c +c c +c +c

x z y z z xy

 

 

y 11 y q y r y s 12 x q x r x s 13 z q z r z s

 

Note that here the normal is represented by the normalized metric (evaluated along the free

s

,

s +s +s .

s =


s

x y z

s +s +s

+s c (q w + r w + s w ) + c (q u + r u + s u ) + c (q v + r v + s v ) = 0

z 33 z q z r z s 13 x q x r x s 13 y q y r y s


+s c (q v + r v + s v ) + c (q u + r u + s u ) + c (q w + r w + s w )

x 11 x q x r x s 12 y q y r y s 13 z q z r z s

s c (q u + r u + s u ) + c (q v + r v + s v ) + c (q w + r w + s w )

xy z yx z x <sup>n</sup> vu vuw vw c +c c +c +c c +c n =0

x

(10)

(11)

(12)

(13)

z

n

To approximate (7)-(9) we discretize the rectangular solid (Figure 2)

**Figure 2.** Grids distributions in curvilinear coordinate.

$$\begin{aligned} \mathbf{q}\_{\mathbf{i}} &= (\mathbf{i} - \mathbf{1}) \mathbf{h}\_{\mathbf{q}} \quad \mathbf{i} = \mathbf{1}, \dots, \mathbf{N}\_{\mathbf{q}} \quad \mathbf{h}\_{\mathbf{q}} = \mathbf{1} \Big\rangle \langle \mathbf{N}\_{\mathbf{q}} \ - \mathbf{1} \rangle; \\ \mathbf{r}\_{\mathbf{j}} &= (\mathbf{j} - \mathbf{1}) \mathbf{h}\_{\mathbf{r}} \quad \mathbf{j} = \mathbf{1}, \dots, \mathbf{N}\_{\mathbf{r}} \quad \mathbf{h}\_{\mathbf{r}} = \mathbf{w} \Big\rangle \langle \mathbf{N}\_{\mathbf{r}} \ - \mathbf{1} \rangle; \\ \mathbf{s}\_{\mathbf{k}} &= (\mathbf{k} - \mathbf{1}) \mathbf{h}\_{\mathbf{s}} \quad \mathbf{k} = \mathbf{1}, \dots, \mathbf{N}\_{\mathbf{s}} \quad \mathbf{h}\_{\mathbf{s}} = \mathbf{h} \Big\rangle \langle \mathbf{N}\_{\mathbf{s}} \ - \mathbf{1} \rangle \end{aligned} \tag{14}$$

where l, w, h are the length of the rectangular solid in q-, r- and s-directions, respectively; qrs h ,h ,h >0 define the grid size in q-, r- and s-directions, respectively. The three components of the wave field are given by

$$[\mathbf{u}\_{\mathrm{i},\mathrm{j},\mathrm{k}}(\mathbf{t}), \mathbf{v}\_{\mathrm{i},\mathrm{j},\mathrm{k}}(\mathbf{t}), \mathbf{w}\_{\mathrm{i},\mathrm{j},\mathrm{k}}(\mathbf{t})] = [\mathbf{u}(\mathbf{q}\_{\mathrm{i}}, \mathbf{r}\_{\mathrm{j}}, \mathbf{s}\_{\mathrm{k}}, \mathbf{t}), \mathbf{v}(\mathbf{q}\_{\mathrm{i}}, \mathbf{r}\_{\mathrm{j}}, \mathbf{s}\_{\mathrm{k}}, \mathbf{t}), \mathbf{w}(\mathbf{q}\_{\mathrm{i}}, \mathbf{r}\_{\mathrm{j}}, \mathbf{s}\_{\mathrm{k}}, \mathbf{t})]$$

and the derivation operators as

$$\begin{split} \mathbf{D}^{q}\_{+} \mathbf{u}\_{i,j,\mathbf{k}} &= \frac{\mathbf{u}\_{i+1,j,\mathbf{k}} \cdot \mathbf{u}\_{i,j,\mathbf{k}}}{\mathbf{h}\_{q}} \\ \mathbf{D}^{r}\_{+} \mathbf{u}\_{i,j,\mathbf{k}} &= \frac{\mathbf{u}\_{i,j+1,\mathbf{k}} \cdot \mathbf{u}\_{i,j,\mathbf{k}}}{\mathbf{h}\_{r}} \cdot \mathbf{D}^{q}\_{+} \mathbf{u}\_{i,j,\mathbf{k}} = \mathbf{D}^{q}\_{+} \mathbf{u}\_{i,j,\mathbf{k}} \cdot \mathbf{D}^{q}\_{-} \mathbf{u}\_{i,j,\mathbf{k}} = \frac{1}{2} (\mathbf{D}^{q}\_{+} \mathbf{u}\_{i,j,\mathbf{k}} + \mathbf{D}^{q}\_{-} \mathbf{u}\_{i,j,\mathbf{k}}) \\ \mathbf{D}^{r}\_{+} \mathbf{u}\_{i,j,\mathbf{k}} &= \frac{\mathbf{u}\_{i,j,\mathbf{k}+1} \cdot \mathbf{u}\_{i,j,\mathbf{k}}}{\mathbf{h}\_{s}} \cdot \mathbf{D}^{s}\_{-} \mathbf{u}\_{i,j,\mathbf{k}} = \mathbf{D}^{s}\_{+} \mathbf{u}\_{i,j,\mathbf{k}-1} \\ \mathbf{D}^{s}\_{+} \mathbf{u}\_{i,j,\mathbf{k}} &= \frac{\mathbf{u}\_{i,j,\mathbf{k}+1} \cdot \mathbf{u}\_{i,j,\mathbf{k}}}{\mathbf{h}\_{s}} \end{split}$$

The right hand sides of eqs. (7)-(9) contain spatial derivatives of nine basic types, which are discretized according to the following equations

 q q q q- + 1 2 r q q 00 s q q 00 (aω ) D (E (a)D ω) q (dω ) D (dD ω) r (gω ) D (gD ω) s q r r 00 rr r r - 12 + s r r 00 (bω ) D (bD ω) q (eω ) D (E (e)D ω) r (mω ) D (mD ω) s q s s 00 r s s 00 ss s s - 12 + (cω ) D (cD ω) q (fω ) D (fD ω) r (pω ) D (E (p)D ω) s (16)

Here *ω* represents u, v or w; a, b, c, d, e, f, g, m and p are combinations of metric and material coefficients. We introduce the following averaging operators:

$$\begin{aligned} \mathbf{E}^{\rm q}\_{1/2}(\mathbf{y}\_{i,\mathbf{j},\mathbf{k}}) &= \mathbf{y}\_{i+1/2,\mathbf{j},\mathbf{k}} = \frac{\mathbf{Y}\_{i,\mathbf{j},\mathbf{k}} + \mathbf{Y}\_{i+1,\mathbf{j},\mathbf{k}}}{2} \\ \mathbf{E}^{\rm r}\_{1/2}(\mathbf{y}\_{i,\mathbf{j},\mathbf{k}}) &= \mathbf{y}\_{i,\mathbf{j}+1/2,\mathbf{k}} = \frac{\mathbf{Y}\_{i,\mathbf{j},\mathbf{k}} + \mathbf{y}\_{i,\mathbf{j}+1,\mathbf{k}}}{2} \\ \mathbf{E}^{\rm s}\_{1/2}(\mathbf{y}\_{i,\mathbf{j},\mathbf{k}}) &= \mathbf{y}\_{i,\mathbf{j},\mathbf{k}+1/2} = \frac{\mathbf{Y}\_{i,\mathbf{j},\mathbf{k}} + \mathbf{y}\_{i,\mathbf{j},\mathbf{k}+1}}{2} \end{aligned} \tag{17}$$

Three-Dimensional Wavefield Simulation in

(20)

(21)

(22)

Heterogeneous Transversely Isotropic Medium with Irregular Free Surface 113

 

<sup>v</sup> Jρ = D [E (M )D v + E (M )D u + E (M )D w]

q qs sss sq qs 0 50 20 40 r rs s sr s rs s 0 50 20 40 s sq q qs q sq q 0 50 20 40 s sr r rs r sr r 0 50 20 40

+D [M D v + M D u + M D w] +D [M D v + M D u + M D w] +D [M D v + M D u + M D w] +D [M D v + M D u + M D w]

r rq q qr q rq q 0 50 20 40

+D [M D v + M D u + M D w]

q qq q qqq qq q qq q 2 - 5+ 2+ 4+ 1 2 1 2 1 2

> rq r r qr 20 40

M D u + M D w]

r r rr r r rr r r rr r - 12 5 + 12 2 + 12 4 + s s ss s s ss s s ss s - 12 5 + 12 2 + 12 4 +

+D [E (M )D v + E (M )D u + E (M )D w] +D [E (M )D v + E (M )D u + E (M )D w]

q qq q qqq qq q qq q 2 - 3+ 4+ 6+ 1 2 1 2 1 2

> rq r r qr 40 60

M D v + M D w]

r r rr r r rr r r rr r - 12 3 + 12 4 + 12 6 + s s ss s s ss s s ss s - 12 3 + 12 4 + 12 6 +

+D [E (M )D u + E (M )D v + E (M )D w] +D [E (M )D u + E (M )D v + E (M )D w]

in the grid points <sup>i</sup> <sup>j</sup> <sup>k</sup> (q ,r ,s ) , qrs (i, j,k) [1,N ] [1,N ] [1,N ]. We have introduced the following notations for the material and metric terms in order to express the discretized

5 x x 66 y y 11 z z 44 6 x x 44 y y 44 z z 33

We discretize in time using second-order accurate centered differences. The full set of

M = Jk l c + Jk l c + Jk l c M = Jk l c + Jk l c + Jk l c

1 x x 11 y y 66 z z 44 2 x y 12 y x 66

M = Jk l c + Jk l c + Jk l c M = Jk l c + Jk l c

3 x z 13 z x 44 4 y z 13 z y 44

M = Jk l c + Jk l c M = Jk l c + Jk l c

 

<sup>w</sup> Jρ = D [E (M )D u + E (M )D v + E (M )D w]

q sq sss sq qs 0 30 40 60 r sr s sr s rs s 0 30 40 60 s qs q qs q sq q 0 30 40 60 s rs r rs r sr r 0 30 40 60

+D [M D u + M D v + M D w] +D [M D u + M D v + M D w] +D [M D u + M D v + M D w] +D [M D u + M D v + M D w]

r qr q qr q rq q 0 30 40 60

+D [M D u + M D v + M D w]

 

 2

t

(v)

 

equations in a more compact form:

discretized equations is

(w)

kl kl

where k and l represent the metric coefficients q, r or s.

L (u,v,w)

kl kl

kl kl

q rq r 0 30

+D [M D u +

L (u,v,w)

q qr r 0 50

+D [M D v +

 2

t

The cross terms which contain a normal derivative on the boundary are discretized onesided in the direction normal to the boundary:

$$
\widehat{\mathbf{D}\_{0}^{s}\mathbf{u}\_{i,\mathbf{j},\mathbf{k}}} = \begin{cases}
\mathbf{D}\_{+}^{s}\mathbf{u}\_{i,\mathbf{j},\mathbf{k}'} & \mathbf{k} = \mathbf{1}, \\
\mathbf{D}\_{0}^{s}\mathbf{u}\_{i,\mathbf{j},\mathbf{k}'} & \mathbf{k} \ge \mathbf{2}. \end{cases} \tag{18}
$$

#### **5.1. A discretization on curvilinear grid: Elastic wave equations**

We approximate the spatial operators in eqs. (7)-(9) by (16). After suppressing grid indexes, this leads to

 2 q qq q qqq qq q qq q 2 - 1+ 2+ 3+ 1 2 1 2 1 2 q qs sss qs qs 0 10 20 30 r rs s rs s rs s 0 10 20 30 s sq q sq q sq q 0 10 20 30 s sr r sr r sr r 0 10 20 30 q qr r 0 10 <sup>u</sup> Jρ = D [E (M )D u + E (M )D v + E (M )D w] t +D [M D u + M D v + M D w] +D [M D u + M D v + M D w] +D [M D u + M D v + M D w] +D [M D u + M D v + M D w] +D [M D u + qr r r qr 20 30 r rq q rq q rq q 0 10 20 30 r r rr r r rr r r rr r - 12 1 + 12 2 + 12 3 + s s ss s s ss s s ss s - 12 1 + 12 2 + 12 3 + (u) M D v +M D w] +D [M D u + M D v + M D w] +D [E (M )D u + E (M )D v + E (M )D w] +D [E (M )D u + E (M )D v + E (M )D w] L (u, v,w) (19)

 2 q qq q qqq qq q qq q 2 - 5+ 2+ 4+ 1 2 1 2 1 2 q qs sss sq qs 0 50 20 40 r rs s sr s rs s 0 50 20 40 s sq q qs q sq q 0 50 20 40 s sr r rs r sr r 0 50 20 40 q qr r 0 50 <sup>v</sup> Jρ = D [E (M )D v + E (M )D u + E (M )D w] t +D [M D v + M D u + M D w] +D [M D v + M D u + M D w] +D [M D v + M D u + M D w] +D [M D v + M D u + M D w] +D [M D v + rq r r qr 20 40 r rq q qr q rq q 0 50 20 40 r r rr r r rr r r rr r - 12 5 + 12 2 + 12 4 + s s ss s s ss s s ss s - 12 5 + 12 2 + 12 4 + (v) M D u + M D w] +D [M D v + M D u + M D w] +D [E (M )D v + E (M )D u + E (M )D w] +D [E (M )D v + E (M )D u + E (M )D w] L (u,v,w) (20) 2 q qq q qqq qq q qq q 2 - 3+ 4+ 6+ 1 2 1 2 1 2 q sq sss sq qs 0 30 40 60 r sr s sr s rs s 0 30 40 60 s qs q qs q sq q 0 30 40 60 s rs r rs r sr r 0 30 40 60 q rq r 0 30 <sup>w</sup> Jρ = D [E (M )D u + E (M )D v + E (M )D w] t +D [M D u + M D v + M D w] +D [M D u + M D v + M D w] +D [M D u + M D v + M D w] +D [M D u + M D v + M D w] +D [M D u + rq r r qr 40 60 r qr q qr q rq q 0 30 40 60 r r rr r r rr r r rr r M D v + M D w] +D [M D u + M D v + M D w] +D [E (M )D u + E (M )D v + E (M )D w] (21)

$$\begin{aligned} & \mathbf{^0D\_{-1}}\mathbf{^1D\_{1/2}}(\mathbf{^0W\_+}\mathbf{^1D\_{-1/2}}\mathbf{^0D\_+}\mathbf{^1D\_+}\mathbf{^0D\_+}\mathbf{^1D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^0D\_+}\mathbf{^$$

in the grid points <sup>i</sup> <sup>j</sup> <sup>k</sup> (q ,r ,s ) , qrs (i, j,k) [1,N ] [1,N ] [1,N ]. We have introduced the following notations for the material and metric terms in order to express the discretized equations in a more compact form:

kl kl 1 x x 11 y y 66 z z 44 2 x y 12 y x 66 kl kl 3 x z 13 z x 44 4 y z 13 z y 44 kl kl 5 x x 66 y y 11 z z 44 6 x x 44 y y 44 z z 33 M = Jk l c + Jk l c + Jk l c M = Jk l c + Jk l c M = Jk l c + Jk l c M = Jk l c + Jk l c M = Jk l c + Jk l c + Jk l c M = Jk l c + Jk l c + Jk l c (22)

where k and l represent the metric coefficients q, r or s.

112 Earthquake Engineering

q

r

s

this leads to

(gω ) D (gD ω)

(dω ) D (dD ω)

q q q q- + 1 2

q

r

s

material coefficients. We introduce the following averaging operators:

1 2 i,j,k i+1 2,j,k

1 2 i,j,k i,j+1 2,k

1 2 i,j,k i,j,k+1 2

 

Du =

**5.1. A discretization on curvilinear grid: Elastic wave equations** 

Here *ω* represents u, v or w; a, b, c, d, e, f, g, m and p are combinations of metric and

q i,j,k i+1,j,k

<sup>γ</sup> <sup>+</sup> <sup>γ</sup> E (γ ) = <sup>γ</sup> <sup>=</sup> <sup>2</sup>

<sup>γ</sup> <sup>+</sup> <sup>γ</sup> E (γ ) = <sup>γ</sup> <sup>=</sup> <sup>2</sup>

<sup>γ</sup> <sup>+</sup> <sup>γ</sup> E (γ ) = <sup>γ</sup> <sup>=</sup> <sup>2</sup>

r i,j,k i,j+1,k

s i,j,k i,j,k+1

The cross terms which contain a normal derivative on the boundary are discretized one-

s s + i,j,k 0 i,j,k s

q qq q qqq qq q qq q 2 - 1+ 2+ 3+ 1 2 1 2 1 2

> qr r r qr 20 30

M D v +M D w]

r r rr r r rr r r rr r - 12 1 + 12 2 + 12 3 + s s ss s s ss s s ss s - 12 1 + 12 2 + 12 3 +

+D [E (M )D u + E (M )D v + E (M )D w] +D [E (M )D u + E (M )D v + E (M )D w]

D u , k = 1,

0 i,j,k

We approximate the spatial operators in eqs. (7)-(9) by (16). After suppressing grid indexes,

<sup>u</sup> Jρ = D [E (M )D u + E (M )D v + E (M )D w]

 

q qs sss qs qs 0 10 20 30 r rs s rs s rs s 0 10 20 30 s sq q sq q sq q 0 10 20 30 s sr r sr r sr r 0 10 20 30

+D [M D u + M D v + M D w] +D [M D u + M D v + M D w] +D [M D u + M D v + M D w] +D [M D u + M D v + M D w]

r rq q rq q rq q 0 10 20 30

+D [M D u + M D v + M D w]

(mω ) D (mD ω)

(eω ) D (E (e)D ω)

(bω ) D (bD ω)

q r r 00

rr r r - 12 + s r r 00

q

r

s

D u , k 2. (18)

q s s 00

(cω ) D (cD ω)

(fω ) D (fD ω)

ss s s - 12 +

(17)

(19)

(16)

r s s 00

(pω ) D (E (p)D ω)

r q q 00

(aω ) D (E (a)D ω)

s q q 00

sided in the direction normal to the boundary:

 

(u)

L (u, v,w)

q qr r 0 10

+D [M D u +

 2

t

We discretize in time using second-order accurate centered differences. The full set of discretized equations is

$$\rho(\frac{\mathbf{u}^{\text{n}^{n+1}} \cdot 2\mathbf{u}^{\text{n}} + \mathbf{u}^{\text{n}^{n-1}}}{\delta\_{\text{t}}^{2}}) = \mathbf{L}^{(\text{u})}(\mathbf{u}^{\text{n}}, \mathbf{v}^{\text{n}}, \mathbf{w}^{\text{n}})$$

$$\rho(\frac{\mathbf{v}^{\text{n}^{n+1}} \cdot 2\mathbf{v}^{\text{n}} + \mathbf{v}^{\text{n} \cdot 1}}{\delta\_{\text{t}}^{2}}) = \mathbf{L}^{(\text{v})}(\mathbf{u}^{\text{n}}, \mathbf{v}^{\text{n}}, \mathbf{w}^{\text{n}}) \tag{23}$$

$$\rho(\frac{\mathbf{w}^{\text{n}^{n+1}} \cdot 2\mathbf{w}^{\text{n}} + \mathbf{w}^{\text{n} \cdot 1}}{\delta\_{\text{t}}^{2}}) = \mathbf{L}^{(\text{w})}(\mathbf{u}^{\text{n}}, \mathbf{v}^{\text{n}}, \mathbf{w}^{\text{n}})$$

Three-Dimensional Wavefield Simulation in

Heterogeneous Transversely Isotropic Medium with Irregular Free Surface 115

order error on the boundary in the differential equations (19)-(21) can be absorbed as a

The accuracy of the proposed method is examined by comparing numerical results with the analytical solution of the Lamb's problem, for a transversely isotropic medium with a vertical symmetry axis (VTI medium). The elastic parameters describing the VTI medium are given in Table 1. The analytical solution is obtained by convolving the free-surface Green-function with the source function (Payton, 1983). A vertical point source of the type

25.5 2.0 14.0 18.4 5.6 2.4

2 2 0 0 -0.5f (t -t )

2000 m) at the surface, which is marked as an asterisk in Figure 3. It should be mentioned that Carcione et al. (1992) and Carcione (2000) presented an analytical comparison of the point-source response in a 3-D VTI medium in the absence of the free surface. The comparisons are performed by first transforming the 3-D numerical results into a line-source response by carrying out an integration along the receiver line (Wapenaar et al., 1992) and then comparing the emerging results with the 2-D Lamb's analytical solutions. The numerical model contains 401 x 401 x 191 grid nodes in the x-, y- and z-direction, respectively. The grid spacings are 10 m in all directions. The solution is advanced using a

33 c (GPa)

0 0 f(t) = e cosπf (t - t ) (27)

f = 10 Hz, is assumed to be located at (300 m,

44 c (GPa)

*ρ* <sup>3</sup> (g cm )

13 c (GPa)

second-order perturbation of the boundary conditions (24)-(26).

**Figure 3.** Model of a half-space with a planar free surface.

12 c (GPa)

with 0t = 0.5 s and a high cut-off frequency <sup>0</sup>

time step of 1.25 ms for 3.5 s.

**Table 1.** Medium parameters in the homogeneous half-space

**6. Accuracy and efficiency tests** 

**6.1. Accuracy** 

11 c (GPa)

where <sup>t</sup> δ represents the time step.

#### **5.2. A discretization on curvilinear grid: Free boundary conditions**

The boundary conditions (11)-(13) are discretized by

$$\begin{aligned} &\frac{1}{2} [\mathbf{M}\_{1}^{\rm ss})\_{i,j,j/2} \mathbf{D}\_{+}^{\rm s} \mathbf{u}\_{i,j,1} + (\mathbf{M}\_{1}^{\rm ss})\_{i,j,1/2} \mathbf{D}\_{+}^{\rm s} \mathbf{u}\_{i,j,0}] + (\mathbf{M}\_{1}^{\rm sq})\_{i,j,1} \mathbf{D}\_{0}^{\rm s} \mathbf{u}\_{i,j,1} + (\mathbf{M}\_{2}^{\rm sq})\_{i,j,1} \mathbf{D}\_{0}^{\rm s} \mathbf{v}\_{i,j,1} \\ &+ (\mathbf{M}\_{3}^{\rm ss})\_{i,j,1} \mathbf{D}\_{0}^{\rm s} \mathbf{w}\_{i,j,1} + \frac{1}{2} [(\mathbf{M}\_{2}^{\rm ss})\_{i,j,2} \mathbf{D}\_{+}^{\rm s} \mathbf{v}\_{i,j,1} + (\mathbf{M}\_{2}^{\rm ss})\_{i,j,1} \mathbf{D}\_{+}^{\rm s} \mathbf{v}\_{i,j,0}] + (\mathbf{M}\_{1}^{\rm ss})\_{i,j,1} \mathbf{D}\_{0}^{\rm s} \mathbf{u}\_{i,j,1} \\ &+ (\mathbf{M}\_{2}^{\rm ss})\_{i,j,1} \mathbf{D}\_{0}^{\rm s} \mathbf{v}\_{i,j,1} + (\mathbf{M}\_{3}^{\rm ss})\_{i,j,1} \mathbf{D}\_{0}^{\rm s} \mathbf{w}\_{i,j,1} + \frac{1}{2} [(\mathbf{M}\_{3}^{\rm ss})\_{i,j,2} \mathbf{D}\_{+}^{\rm s} \mathbf{w}\_{i,j,1} + (\mathbf{M}\_{3}^{\rm ss})\_{i,j,1} \mathbf{D}\_{+}^{\rm s} \mathbf{w}\_{i,j,0}] = 0 \end{aligned} \tag{24}$$

$$\begin{aligned} &\frac{1}{2}[\mathbf{M}\_{5}^{ss}]\_{i,j,j/2}\mathbf{D}\_{+}^{s}\mathbf{v}\_{i,j,1}+(\mathbf{M}\_{5}^{ss})\_{i,j,j/2}\mathbf{D}\_{+}^{s}\mathbf{v}\_{i,j,0}]+(\mathbf{M}\_{5}^{sq})\_{i,j,1}\mathbf{D}\_{0}^{q}\mathbf{v}\_{i,j,1}+(\mathbf{M}\_{2}^{qs})\_{i,j,1}\mathbf{D}\_{0}^{q}\mathbf{u}\_{i,j,1} \\ &+(\mathbf{M}\_{4}^{sq})\_{i,j,1}\mathbf{D}\_{0}^{q}\mathbf{w}\_{i,j,1}+\frac{1}{2}[(\mathbf{M}\_{2}^{ss})\_{i,j,j/2}\mathbf{D}\_{+}^{s}\mathbf{u}\_{i,j,1}+(\mathbf{M}\_{2}^{ss})\_{i,j,j/2}\mathbf{D}\_{+}^{s}\mathbf{u}\_{i,j,0}]+(\mathbf{M}\_{5}^{sr})\_{i,j,1}\mathbf{D}\_{0}^{r}\mathbf{v}\_{i,j,1} \\ &+(\mathbf{M}\_{2}^{sr})\_{i,j,1}\mathbf{D}\_{0}^{r}\mathbf{u}\_{i,j,1}+(\mathbf{M}\_{4}^{sr})\_{i,j,1}\mathbf{D}\_{0}^{r}\mathbf{w}\_{i,j,1}+\frac{1}{2}[(\mathbf{M}\_{4}^{ss})\_{i,j,3/2}\mathbf{D}\_{+}^{s}\mathbf{w}\_{i,j,1}+(\mathbf{M}\_{4}^{sr})\_{i,j,1}\mathbf{D}\_{2}^{s}\mathbf{w}\_{i,j,0}]=0 \end{aligned} \tag{25}$$

$$\frac{1}{2} [\mathbf{M}\_3^{\rm ss})\_{i,j,j/2} \mathbf{D}\_+^s \mathbf{u}\_{i,j,1} + (\mathbf{M}\_3^{\rm ss})\_{i,j,j/2} \mathbf{D}\_+^s \mathbf{u}\_{i,j,0}] + (\mathbf{M}\_3^{\rm ss})\_{i,j,1} \mathbf{D}\_0^s \mathbf{u}\_{i,j,1} + (\mathbf{M}\_4^{\rm ss})\_{i,j,1} \mathbf{D}\_0^s \mathbf{v}\_{i,j,1}$$

$$+ (\mathbf{M}\_6^{\rm ss})\_{i,j,1} \mathbf{D}\_0^4 \mathbf{w}\_{i,j,1} + \frac{1}{2} [(\mathbf{M}\_4^{\rm ss})\_{i,j,j/2} \mathbf{D}\_+^s \mathbf{v}\_{i,j,1} + (\mathbf{M}\_4^{\rm ss})\_{i,j,1} \mathbf{D}\_+^s \mathbf{v}\_{i,j,0}] + (\mathbf{M}\_3^{\rm ss})\_{i,j,1} \mathbf{D}\_0^s \mathbf{u}\_{i,j,1} \tag{26}$$

$$+ (\mathbf{M}\_4^{\rm ss})\_{i,j,1} \mathbf{D}\_0^s \mathbf{v}\_{i,j,1} + (\mathbf{M}\_6^{\rm ss})\_{i,j,1} \mathbf{D}\_0^s \mathbf{w}\_{i,j,1} + \frac{1}{2} [(\mathbf{M}\_6^{\rm ss})\_{i,j,2} \mathbf{D}\_\*^s \mathbf{w}\_{i,j,1} + (\mathbf{M}\_6^{\rm ss})\_{i,j,1} \mathbf{D}\_\*^s \mathbf{w}\_{i,j,0}] = 0$$

$$\mathbf{i} = \mathbf{1}, \dots, \mathbf{N}\_{\rm q}; j = 1, \dots, \mathbf{N}\_{\rm r}.$$

The key step in obtaining a stable explicit discretization is to use the operator <sup>s</sup> D0 ( which is one-sided on the boundary) for the approximation of the normal derivative in qs rs sq , , and *s r* cross derivatives. At first glance, it may appear that using a onesided operator the accuracy of the method would be reduced to the first-order. However, as it was theoretically shown by Nilsson et al. (2007) (for a Cartesian discretization), a firstorder error on the boundary in the differential equations (19)-(21) can be absorbed as a second-order perturbation of the boundary conditions (24)-(26).

**Figure 3.** Model of a half-space with a planar free surface.

#### **6. Accuracy and efficiency tests**

#### **6.1. Accuracy**

114 Earthquake Engineering

where <sup>t</sup> δ represents the time step.

sr r 2 i,j,1 0 i,j,1

rs r 2 i,j,1 0 i,j,1

rs r 4 i,j,1 0 i,j,1

1

2

1

2

1

2

The boundary conditions (11)-(13) are discretized by

n+1 n n-1

u - 2u + u

v - 2v + v

2 t n+1 n n-1

ρ( ) = L (u , v ,w ) <sup>δ</sup>

ρ( ) = L (u , v ,w ) <sup>δ</sup>

ρ( ) <sup>=</sup> <sup>L</sup> (u ,v ,w ) <sup>δ</sup>

ss s ss s sq q sq q 1 i,j,3 2 + i,j,1 1 i,j,1 2 + i,j,0 1 i,j,1 0 i,j,1 2 i,j,1 0 i,j,1

[(M ) D u + (M ) D u ]+ (M ) D u + (M ) D v

<sup>1</sup> +(M ) D w + [(M ) D v + (M ) D v ]+ (M ) D u <sup>2</sup>

+(M ) D v sr r ss s ss s

sq q ss s ss s sr r 3 i,j,1 0 i,j,1 2 i,j,3 2 + i,j,1 2 i,j,1 2 + i,j,0 1 i,j,1 0 i,j,1

ss s ss s sq q qs q 5 i,j,3 2 + i,j,1 5 i,j,1 2 + i,j,0 5 i,j,1 0 i,j,1 2 i,j,1 0 i,j,1 sq q ss s ss s sr r 4 i,j,1 0 i,j,1 2 i,j,3 2 + i,j,1 2 i,j,1 2 + i,j,0 5 i,j,1 0 i,j,1

[(M ) D v + (M ) D v ]+ (M ) D v + (M ) D u

<sup>1</sup> +(M ) D w + [(M ) D u + (M ) D u ]+ (M ) D v <sup>2</sup>

+(M ) D u sr r ss s ss s

ss s ss s qs q qs q 3 i,j,3 2 + i,j,1 3 i,j,1 2 + i,j,0 3 i,j,1 0 i,j,1 4 i,j,1 0 i,j,1

[(M ) D u + (M ) D u ]+ (M ) D u + (M ) D v

<sup>1</sup> +(M ) D w + [(M ) D v + (M ) D v ]+ (M ) D u <sup>2</sup>

+(M ) D v sr r ss s ss s

sq q ss s ss s rs r 6 i,j,1 0 i,j,1 4 i,j,3 2 + i,j,1 4 i,j,1 2 + i,j,0 3 i,j,1 0 i,j,1

q r i = 1,...,N ; j = 1,...,N .

The key step in obtaining a stable explicit discretization is to use the operator <sup>s</sup> D0 ( which is one-sided on the boundary) for the approximation of the normal derivative in qs rs sq , , and *s r* cross derivatives. At first glance, it may appear that using a onesided operator the accuracy of the method would be reduced to the first-order. However, as it was theoretically shown by Nilsson et al. (2007) (for a Cartesian discretization), a first-

2 t n+1 n n-1

w - 2w + w

2 t

**5.2. A discretization on curvilinear grid: Free boundary conditions** 

(u) nn n

(v) nn n

3 i,j,1 0 i,j,1 3 i,j,3 2 + i,j,1 3 i,j,1 2 + i,j,0 <sup>1</sup> + (M ) D w + [(M ) D w + (M ) D w ] = 0 <sup>2</sup>

4 i,j,1 0 i,j,1 4 i,j,3 2 + i,j,1 4 i,j,1 2 + i,j,0 <sup>1</sup> + (M ) D w + [(M ) D w + (M ) D w ] = 0 <sup>2</sup>

6 i,j,1 0 i,j,1 6 i,j,3 2 + i,j,1 6 i,j,1 2 + i,j,0 <sup>1</sup> + (M ) D w + [(M ) D w + (M ) D w ] = 0 <sup>2</sup>

(w) nn n

(23)

(24)

(25)

(26)

The accuracy of the proposed method is examined by comparing numerical results with the analytical solution of the Lamb's problem, for a transversely isotropic medium with a vertical symmetry axis (VTI medium). The elastic parameters describing the VTI medium are given in Table 1. The analytical solution is obtained by convolving the free-surface Green-function with the source function (Payton, 1983). A vertical point source of the type


**Table 1.** Medium parameters in the homogeneous half-space

$$\mathbf{f(t) = e^{-0.5t\_0^2 \left(t - t\_0\right)^2} \cos \mathbf{n} f\_0 \left(t - t\_0\right)}\tag{27}$$

with 0t = 0.5 s and a high cut-off frequency <sup>0</sup> f = 10 Hz, is assumed to be located at (300 m, 2000 m) at the surface, which is marked as an asterisk in Figure 3. It should be mentioned that Carcione et al. (1992) and Carcione (2000) presented an analytical comparison of the point-source response in a 3-D VTI medium in the absence of the free surface. The comparisons are performed by first transforming the 3-D numerical results into a line-source response by carrying out an integration along the receiver line (Wapenaar et al., 1992) and then comparing the emerging results with the 2-D Lamb's analytical solutions. The numerical model contains 401 x 401 x 191 grid nodes in the x-, y- and z-direction, respectively. The grid spacings are 10 m in all directions. The solution is advanced using a time step of 1.25 ms for 3.5 s.

Three-Dimensional Wavefield Simulation in

Heterogeneous Transversely Isotropic Medium with Irregular Free Surface 117

**Figure 6.** Snapshots of the vertical component of the wavefield at the surface (xy-plane) of the planar

Synthetic seismograms are computed at Line 3. The seismograms in Figure 5 show the direct quasi-P wave (qP) and a high-amplitude Rayleigh wave (R). Snapshots of the vertical component of the wavefield in horizontal (xy-) plane at the propagation time of 1.4 s are displayed in Figure 6. We define the incidence plane by the propagation direction and the zaxis, quasi-P wave and quasi-SV wave (qSV) motions lie in this plane, while SH motion is normal to the plane. Hence, the z-component does not contain SH motion. The xy-plane of a transversely isotropic medium is a plane of isotropy, where the velocity of the qP wave is about 3260 -1 ms and the velocity of the qSV wave is about 1528 -1 ms . The amplitude of the qP wave is so weak compared with that of the Rayleigh and qSV wave that one can hardly identify it in the snapshot (Figure 6a). In order to observe the qP wave, a gain has been given to the amplitudes of the wavefield. Owing to this, side reflections also appear in the photo, as shown in Figure 6b. As the velocity of the Rayleigh wave is very close to that of the qSV wave, the two waves are almost superimposed and it is difficult to distinguish

**Figure 7.** Snapshots of the x-component of the wavefield in the vertical (xz-) plane which contains the receiver line and the source at 1.4 (a) and 2.3 s (b) propagation times for the planar surface model.

between the two in synthetic seismograms and snapshots.

surface model.

**Figure 4.** Comparison between numerical and analytical vertical components of the displacement for the VTI medium.

Three receiver lines are positioned on the free surface, two of which are parallel to the ydirection with respective normal distances of 130 (Line 1) and 1000 m (Line 2) away from the point source, the other crosses the source location and parallels to the x-direction (Line 3). The integrations are performed along the first two receiver lines, these represent 2-D results of 130 and 1000 m away from the source, respectively. Figure 4 shows the comparisons between the resulting numerical and 2-D analytical z-components of the displacement for the VTI medium. In spite of the errors resulting from the transformation of the point-source response into the line-source one, numerical and analytical results agree well in Figure 4. These comparisons demonstrate the accuracy of our corresponding algorithm.

**Figure 5.** Seismogram sections at Line 3 for the planar surface model. Symbol qP indicates the qP wave and R indicates Rayleigh wave.

Three-Dimensional Wavefield Simulation in Heterogeneous Transversely Isotropic Medium with Irregular Free Surface 117

116 Earthquake Engineering

the VTI medium.

and R indicates Rayleigh wave.

**Figure 4.** Comparison between numerical and analytical vertical components of the displacement for

Three receiver lines are positioned on the free surface, two of which are parallel to the ydirection with respective normal distances of 130 (Line 1) and 1000 m (Line 2) away from the point source, the other crosses the source location and parallels to the x-direction (Line 3). The integrations are performed along the first two receiver lines, these represent 2-D results of 130 and 1000 m away from the source, respectively. Figure 4 shows the comparisons between the resulting numerical and 2-D analytical z-components of the displacement for the VTI medium. In spite of the errors resulting from the transformation of the point-source response into the line-source one, numerical and analytical results agree well in Figure 4.

**Figure 5.** Seismogram sections at Line 3 for the planar surface model. Symbol qP indicates the qP wave

These comparisons demonstrate the accuracy of our corresponding algorithm.

**Figure 6.** Snapshots of the vertical component of the wavefield at the surface (xy-plane) of the planar surface model.

Synthetic seismograms are computed at Line 3. The seismograms in Figure 5 show the direct quasi-P wave (qP) and a high-amplitude Rayleigh wave (R). Snapshots of the vertical component of the wavefield in horizontal (xy-) plane at the propagation time of 1.4 s are displayed in Figure 6. We define the incidence plane by the propagation direction and the zaxis, quasi-P wave and quasi-SV wave (qSV) motions lie in this plane, while SH motion is normal to the plane. Hence, the z-component does not contain SH motion. The xy-plane of a transversely isotropic medium is a plane of isotropy, where the velocity of the qP wave is about 3260 -1 ms and the velocity of the qSV wave is about 1528 -1 ms . The amplitude of the qP wave is so weak compared with that of the Rayleigh and qSV wave that one can hardly identify it in the snapshot (Figure 6a). In order to observe the qP wave, a gain has been given to the amplitudes of the wavefield. Owing to this, side reflections also appear in the photo, as shown in Figure 6b. As the velocity of the Rayleigh wave is very close to that of the qSV wave, the two waves are almost superimposed and it is difficult to distinguish between the two in synthetic seismograms and snapshots.

**Figure 7.** Snapshots of the x-component of the wavefield in the vertical (xz-) plane which contains the receiver line and the source at 1.4 (a) and 2.3 s (b) propagation times for the planar surface model.

Figure 7 shows the x-component of the wavefield in the vertical (xz-) plane at 1.4 and 2.3 s propagation times. The xz-plane contains the receiver line (Line 3) and the source location. Both snapshots show the wave front of the qP-wave and the qSV-wave. The former snapshot (1.4 s) shows the qSV-wave with the cusps. A headwave (H) can also be found in the photos, the headwave is a quasi-shear wave and is guided along the surface by the qP-wave.

Three-Dimensional Wavefield Simulation in

Heterogeneous Transversely Isotropic Medium with Irregular Free Surface 119

hill-shaped free surface (and compare with the synthetic seismograms in Figure 5), the amplitudes of the quasi-P wave and Rayleigh wave are reduced in the right-hand part of the sections. In addition, after the ordinary quasi-P wave a secondary quasi-P wave (RqPf) induced by the scattering of the direct Rayleigh wave can be observed. Similarly, a secondary Rayleigh wave (qPRf) which travels in front of the ordinary Rayleigh wave induced by the scattering of the direct quasi-P wave can also be distinguished. Some energy is scattered back

**Figure 9.** The gridding scheme which shows the detailed cross-section of the grids along Line3 in the Gaussian shape hill topography model. For clarity, the grids are displayed with a reducing density

**Figure 10.** Seismograms along the receiver line for the Gaussian shape hill topography model: (a) xcomponent (horizontal) of the displacement; (b) z-component (vertical). Symbols mean the following: (qPd) qP wave diffracts to qP wave; (Rd) Rayleigh wave diffracts to Rayleigh wave; (qPRf) qP wave scatters to Rayleigh wave and propagates forward; (qPRb) qP wave scatters to Rayleigh wave and propagates backward; (RqPf) Rayleigh wave scatters to qP wave and propagates forward; (RqPb) Rayleigh wave scatters to qP wave and propagates backward; (RR) Rayleigh wave reflectes to Rayleigh

factor of 3.

wave.

to the left-hand side as a Rayleigh wave (qPRb, RR) and a quasi-P wave (RqPb).

#### **6.2. Numerical simulations on an irregular (non-flat) free surface**

Three numerical experiments with irregular free surfaces are now investigated. The first example is a test on smooth boundaries, while the second example consists of a hemispherical depression to test the ability of the method for non-smooth topography. For sake of simplicity both these examples are based on homogeneous half-spaces, i.e., the medium parameters are the same as in the case of flat surface (Table 1). The same source is located at the same place as in the planar surface model, the time step is 0.8 ms. The total propagation time is 3.5 s for the two models. Finally, we consider a two-layered model with a realistic topography.

#### *6.2.1. Topography simulating a shaped Gaussian hill*

The first model considered here is a half-space whose free surface is a hill-like feature (Figure 8). The shape of the hill resembles a Gaussian curved surface given by the function

**Figure 8.** Model of a half-space with Gaussian shape hill topography.

The computational domain extends to depth z(x, y)=2000 m. The volume is discretized with equal grid nodes in each direction as in the planar surface one. The grid spacings are 10 m in the x, y-directions and about 10.5 m in the z-direction for average. The vertical spacing varies with depth, it is smaller toward the free surface and larger toward the bottom of the model. The minimum and maximum of the vertical spacings are 6 and 12 m, respectively.

The gridding scheme which shows the detailed cross-section of the grids along Line3 is shown in Figure 9. Synthetic seismograms are also computed at Line3 (Figure 10). As a result of the hill-shaped free surface (and compare with the synthetic seismograms in Figure 5), the amplitudes of the quasi-P wave and Rayleigh wave are reduced in the right-hand part of the sections. In addition, after the ordinary quasi-P wave a secondary quasi-P wave (RqPf) induced by the scattering of the direct Rayleigh wave can be observed. Similarly, a secondary Rayleigh wave (qPRf) which travels in front of the ordinary Rayleigh wave induced by the scattering of the direct quasi-P wave can also be distinguished. Some energy is scattered back to the left-hand side as a Rayleigh wave (qPRb, RR) and a quasi-P wave (RqPb).

118 Earthquake Engineering

a realistic topography.

Figure 7 shows the x-component of the wavefield in the vertical (xz-) plane at 1.4 and 2.3 s propagation times. The xz-plane contains the receiver line (Line 3) and the source location. Both snapshots show the wave front of the qP-wave and the qSV-wave. The former snapshot (1.4 s) shows the qSV-wave with the cusps. A headwave (H) can also be found in the photos,

Three numerical experiments with irregular free surfaces are now investigated. The first example is a test on smooth boundaries, while the second example consists of a hemispherical depression to test the ability of the method for non-smooth topography. For sake of simplicity both these examples are based on homogeneous half-spaces, i.e., the medium parameters are the same as in the case of flat surface (Table 1). The same source is located at the same place as in the planar surface model, the time step is 0.8 ms. The total propagation time is 3.5 s for the two models. Finally, we consider a two-layered model with

The first model considered here is a half-space whose free surface is a hill-like feature (Figure 8). The shape of the hill resembles a Gaussian curved surface given by the function

The computational domain extends to depth z(x, y)=2000 m. The volume is discretized with equal grid nodes in each direction as in the planar surface one. The grid spacings are 10 m in the x, y-directions and about 10.5 m in the z-direction for average. The vertical spacing varies with depth, it is smaller toward the free surface and larger toward the bottom of the model. The minimum and maximum of the vertical spacings are 6 and 12 m, respectively.

The gridding scheme which shows the detailed cross-section of the grids along Line3 is shown in Figure 9. Synthetic seismograms are also computed at Line3 (Figure 10). As a result of the

x - 2000 2 2 y - 2000 z(x,y) = -150exp(-( ) - ( ) )m 150 150 <sup>2</sup> (x, y) [0m, 4000m] (28)

the headwave is a quasi-shear wave and is guided along the surface by the qP-wave.

**6.2. Numerical simulations on an irregular (non-flat) free surface** 

*6.2.1. Topography simulating a shaped Gaussian hill* 

**Figure 8.** Model of a half-space with Gaussian shape hill topography.

**Figure 9.** The gridding scheme which shows the detailed cross-section of the grids along Line3 in the Gaussian shape hill topography model. For clarity, the grids are displayed with a reducing density factor of 3.

**Figure 10.** Seismograms along the receiver line for the Gaussian shape hill topography model: (a) xcomponent (horizontal) of the displacement; (b) z-component (vertical). Symbols mean the following: (qPd) qP wave diffracts to qP wave; (Rd) Rayleigh wave diffracts to Rayleigh wave; (qPRf) qP wave scatters to Rayleigh wave and propagates forward; (qPRb) qP wave scatters to Rayleigh wave and propagates backward; (RqPf) Rayleigh wave scatters to qP wave and propagates forward; (RqPb) Rayleigh wave scatters to qP wave and propagates backward; (RR) Rayleigh wave reflectes to Rayleigh wave.

Three-Dimensional Wavefield Simulation in

Heterogeneous Transversely Isotropic Medium with Irregular Free Surface 121

two parts, one travels forward (qPRf), and the other travels backward (qPRb). These can be seen clearly in the later snapshots (1.4 - 2.3 s). In addition, a reflected Rayleigh wave (RR) can be observed. The direct quasi-P wave (qP) and Rayleigh wave (R) are also marked in the figure. By the way, side reflections from the boundaries can also be noted in the plane. Figure 12 shows the x-component of the wavefield in the vertical (xz-) plane. The xz-plane contains the receiver line and source location. The snapshots show the diffracted quasi-P

**Figure 12.** Snapshots of the x-component of the wavefield in the vertical (xz-) plane which contains the receiver line and the source at different propagation times for the Gaussian shape hill topography model.

and quasi-SV waves clearly in the vertical plane.

**Figure 11.** Snapshots of the vertical component of the wavefield at the surface (xy-plane) at different propagation times for the Gaussian shape hill topography model.

Snapshots of the wavefield in the horizontal (xy-) plane at different propagation times are displayed in Figures 11. The amplitudes are also gained. In the beginning the wavefield propagates undisturbed along the free surface. At 1.1 s the direct quasi-P wave hits the hill and generates a circular diffracted wave. This wave is a Rayleigh wave, which is marked as two parts, one travels forward (qPRf), and the other travels backward (qPRb). These can be seen clearly in the later snapshots (1.4 - 2.3 s). In addition, a reflected Rayleigh wave (RR) can be observed. The direct quasi-P wave (qP) and Rayleigh wave (R) are also marked in the figure. By the way, side reflections from the boundaries can also be noted in the plane. Figure 12 shows the x-component of the wavefield in the vertical (xz-) plane. The xz-plane contains the receiver line and source location. The snapshots show the diffracted quasi-P and quasi-SV waves clearly in the vertical plane.

120 Earthquake Engineering

**Figure 11.** Snapshots of the vertical component of the wavefield at the surface (xy-plane) at different

Snapshots of the wavefield in the horizontal (xy-) plane at different propagation times are displayed in Figures 11. The amplitudes are also gained. In the beginning the wavefield propagates undisturbed along the free surface. At 1.1 s the direct quasi-P wave hits the hill and generates a circular diffracted wave. This wave is a Rayleigh wave, which is marked as

propagation times for the Gaussian shape hill topography model.

**Figure 12.** Snapshots of the x-component of the wavefield in the vertical (xz-) plane which contains the receiver line and the source at different propagation times for the Gaussian shape hill topography model.

#### *6.2.2. Topography simulating a shaped hemispherical depression*

In the second model, we consider a hemispherical depression model as illustrated in Figure 13. The first model that we have considered is of smooth topography, that is, with continuous and finite slopes everywhere. However, the shaped hemispherical depression here taken as reference is a case of extreme topography, such that the vertical-to-horizontal ratio of the depression is very large (1:2) and the slopes of the edges tend to infinity. The hemispherical depression is at the middle of the free surface and the radius is 150 m.

Three-Dimensional Wavefield Simulation in

Heterogeneous Transversely Isotropic Medium with Irregular Free Surface 123

of the body wave when propagating through the hemispherical depression, this indicating

**Figure 15.** Seismograms at the receiver line for the hemispherical shape depression topography model: (a) x-component (horizontal) of the displacement; (b) z-component (vertical). Symbols mean the same as

The photos in Figure 16 show the vertical component of the wavefield in the horizontal (xy-) plane. Compared with the photos of the hill topography model, we can see the Rayleigh wave scattering at the edges of the hemispherical depression; it seems as if the reflected Rayleigh wave propagating faster in the hemispherical depression model than that in the hill topography model. What's more, the back scattered waves of Rayleigh wave in the hemispherical depression model are much stronger, this may also indicating that such sharp

It is also interesting to study a realistic example. We consider a model in Tibet (Figure 17). The length and width of the model are 21.6 km, and the "average" height of the topography is roughly -3560 m (3560 m in the geodetic coordinate system). The computational domain is extended to depth z(x, y)= 7200 m. For simplicity we use a two-layered model with parameters given in the model sketch (Figure 17) instead of the "real" velocity structure under the realistic topography. It consists of 241×241×121 grid nodes in the x-, y-, and zdirection, respectively, with equal vertical grid nodes in each layer. A vertical point source like the used in previous models is loaded in the middle of the free surface (indicated by the asterisk in Fig. 17), where the high cut-off frequency has been changed to 2.7 Hz and the time-shift is 1.5 s. Two lines of receivers crossing the source location and paralleling to the xand y-direction respectively are placed at the free surface. The time step is 5 ms, and the

depression blocks the propagation of Rayleigh wave more significantly.

in Fig. 10.

*6.2.3. Real topography simulating* 

total propagation time is 8 s.

that such sharp depression can affect the propagation of Rayleigh wave significantly.

**Figure 13.** Model of a half-space with hemispherical shape depression topography.

**Figure 14.** The gridding scheme which shows the detailed cross-section of the grids along Line3 in the hemispherical shape depression topography model. For clarity, the grids are displayed with a reducing density factor of 3.

The numerical model is discretized in the same way as in the hill topography model. The gridding scheme which shows the detailed cross-section of the grids along Line3 is shown in Figure 14. Owing to the existence of model edges with strong slopes at x=1850 and x=2150 m along the receiver line, both body and Rayleigh waves scattered by sharp changes in the topography can be clearly observed on the synthetic seismograms shown in Figure 15. Owing to its shorter wavelength, the scattering of Rayleigh wave is much stronger than that

of the body wave when propagating through the hemispherical depression, this indicating that such sharp depression can affect the propagation of Rayleigh wave significantly.

**Figure 15.** Seismograms at the receiver line for the hemispherical shape depression topography model: (a) x-component (horizontal) of the displacement; (b) z-component (vertical). Symbols mean the same as in Fig. 10.

The photos in Figure 16 show the vertical component of the wavefield in the horizontal (xy-) plane. Compared with the photos of the hill topography model, we can see the Rayleigh wave scattering at the edges of the hemispherical depression; it seems as if the reflected Rayleigh wave propagating faster in the hemispherical depression model than that in the hill topography model. What's more, the back scattered waves of Rayleigh wave in the hemispherical depression model are much stronger, this may also indicating that such sharp depression blocks the propagation of Rayleigh wave more significantly.

#### *6.2.3. Real topography simulating*

122 Earthquake Engineering

density factor of 3.

*6.2.2. Topography simulating a shaped hemispherical depression* 

In the second model, we consider a hemispherical depression model as illustrated in Figure 13. The first model that we have considered is of smooth topography, that is, with continuous and finite slopes everywhere. However, the shaped hemispherical depression here taken as reference is a case of extreme topography, such that the vertical-to-horizontal ratio of the depression is very large (1:2) and the slopes of the edges tend to infinity. The

hemispherical depression is at the middle of the free surface and the radius is 150 m.

**Figure 13.** Model of a half-space with hemispherical shape depression topography.

**Figure 14.** The gridding scheme which shows the detailed cross-section of the grids along Line3 in the hemispherical shape depression topography model. For clarity, the grids are displayed with a reducing

The numerical model is discretized in the same way as in the hill topography model. The gridding scheme which shows the detailed cross-section of the grids along Line3 is shown in Figure 14. Owing to the existence of model edges with strong slopes at x=1850 and x=2150 m along the receiver line, both body and Rayleigh waves scattered by sharp changes in the topography can be clearly observed on the synthetic seismograms shown in Figure 15. Owing to its shorter wavelength, the scattering of Rayleigh wave is much stronger than that It is also interesting to study a realistic example. We consider a model in Tibet (Figure 17). The length and width of the model are 21.6 km, and the "average" height of the topography is roughly -3560 m (3560 m in the geodetic coordinate system). The computational domain is extended to depth z(x, y)= 7200 m. For simplicity we use a two-layered model with parameters given in the model sketch (Figure 17) instead of the "real" velocity structure under the realistic topography. It consists of 241×241×121 grid nodes in the x-, y-, and zdirection, respectively, with equal vertical grid nodes in each layer. A vertical point source like the used in previous models is loaded in the middle of the free surface (indicated by the asterisk in Fig. 17), where the high cut-off frequency has been changed to 2.7 Hz and the time-shift is 1.5 s. Two lines of receivers crossing the source location and paralleling to the xand y-direction respectively are placed at the free surface. The time step is 5 ms, and the total propagation time is 8 s.

Three-Dimensional Wavefield Simulation in

Heterogeneous Transversely Isotropic Medium with Irregular Free Surface 125

**Figure 17.** A two-layered model with a realistic topography. The medium parameters of each layer are also given in the figure. The units for the elasticity and density are GPa and <sup>3</sup> g cm , respectively.

Snapshots of the z-component of the wavefield in the vertical plane which contains receiver line Line1 and the source location are presented on Figure 18, and the seismograms of the zcomponent are also computed at the two receiver lines (Figure 19). We can see that the effect of the topography is very important, with strong scattered phases that are superimposed to the direct and reflected waves, which make it very difficult to identify effective reflections from subsurface interface. The scattering in the seismograms also reflect different features of the surface. The scattering in the seismograms at Line 1 (Figure 19a) is much stronger than that in the seismogram at Line 2 (Figure 19b), indicating that the surface along Line 1 is much rougher than that along Line 2, which also can be observed in Figure 17. What's more, the scattering in Figure 19a is approximately uniformly distributed while in Figure 19b it is mostly distributed in the vicinity of the shot. These may due to different distributions of the

surface topography along these two lines.

**Figure 16.** Snapshots of the vertical component of the wavefield at the surface (xy-plane) at different propagation times for the hemispherical shape depression topography model.

**Figure 16.** Snapshots of the vertical component of the wavefield at the surface (xy-plane) at different

propagation times for the hemispherical shape depression topography model.

**Figure 17.** A two-layered model with a realistic topography. The medium parameters of each layer are also given in the figure. The units for the elasticity and density are GPa and <sup>3</sup> g cm , respectively.

Snapshots of the z-component of the wavefield in the vertical plane which contains receiver line Line1 and the source location are presented on Figure 18, and the seismograms of the zcomponent are also computed at the two receiver lines (Figure 19). We can see that the effect of the topography is very important, with strong scattered phases that are superimposed to the direct and reflected waves, which make it very difficult to identify effective reflections from subsurface interface. The scattering in the seismograms also reflect different features of the surface. The scattering in the seismograms at Line 1 (Figure 19a) is much stronger than that in the seismogram at Line 2 (Figure 19b), indicating that the surface along Line 1 is much rougher than that along Line 2, which also can be observed in Figure 17. What's more, the scattering in Figure 19a is approximately uniformly distributed while in Figure 19b it is mostly distributed in the vicinity of the shot. These may due to different distributions of the surface topography along these two lines.

Three-Dimensional Wavefield Simulation in

Heterogeneous Transversely Isotropic Medium with Irregular Free Surface 127

We propose a stable and explicit finite difference method to simulate with second-order accuracy the propagation of seismic waves in a 3D heterogeneous transversely isotropic medium with non-flat free surface. The method is an extension of the 2D method proposed by Appelo and Petersson (2009) to the 3D anisotropic case. The surface topography is introduced via mapping rectangular grids to curved grids. The accurate application of the free surface boundary conditions is done by using boundary-modified difference operators to discretize the mixed derivatives in the governing equations of the problem. Several numerical examples under different assumptions of free surface are given to highlight the complications of realistic seismic wave propagation in the vicinity of the earth surface. Synthetic seismograms and snapshots explain diffractions, scattering, multiple reflections, and converted waves provoked by the features of the free surface topography. The typical cuspidal triangles of the quasi-transverse (qS) mode also appear in the snapshots of the

The future directions of our research will include an extension of the schemes to the viscoelastic case. This will allow a realistic attenuation of the seismic waves due to the

We are grateful to Jinlai Hao and Jinhai Zhang for their assistance and the facilities given in the course of this work. Fruitful discussions with Tao Xu, Kai Liang, Wei Zhang are also greatly appreciated. The Ministry of Science and Technology of China (SINOPROBE-02-02)

Al-Shukri, H. J.; Pavlis, G. L. & Vernon F. L. (1995). Site effect observations from broadband

Appelo, D. & Petersson, N. A. (2009). A stable finite difference method for the elastic wave equation on complex geometries with free surfaces. *Commun. Comput. Phys.,* 5, 84–107. Ashford, S. A.; Sitar, N.; Lysmer, J. & Deng N. (1997). Topographic effects on the seismic

Backus, G. E. (1962). Long-wave elastic anisotropy produced by horizontal layering. *J.* 

*State Key Laboratory of Lithosphere Evolution, Institute of Geology and Geophysics,* 

**7. Conclusion** 

anisotropic medium.

**Author details** 

**Acknowledgement** 

supported this research.

**8. References** 

Haiqiang Lan and Zhongjie Zhang

*Chinese Academy of Sciences, China P.R.* 

presence of a weathered layer to be included (2000).

arrays. *Bull. Seismol. Soc. Am.* 85, 1758-1769.

*Geophys. Res.,* 67, 4427-4440.

response of steep slopes. *Bull. Seismol. Soc. Am.,* 87, 701-709.

**Figure 18.** Snapshots of the vertical component of the wavefield in the vertical (xz-) plane along Line 1 at different propagation times for the two-layered model with a realistic topography.

**Figure 19.** Vertical-component synthetic seismograms coming from the two-layered model with real topography represented in Figure 19: (a) at the receiver line (Line 1) that crosses the source location and is parallel to x-direction (Fig. 19 ); (b) at the receiver line (Line 2) that crosses the source location and is parallel to y-direction.

## **7. Conclusion**

126 Earthquake Engineering

parallel to y-direction.

**Figure 18.** Snapshots of the vertical component of the wavefield in the vertical (xz-) plane along Line 1

**Figure 19.** Vertical-component synthetic seismograms coming from the two-layered model with real topography represented in Figure 19: (a) at the receiver line (Line 1) that crosses the source location and is parallel to x-direction (Fig. 19 ); (b) at the receiver line (Line 2) that crosses the source location and is

at different propagation times for the two-layered model with a realistic topography.

We propose a stable and explicit finite difference method to simulate with second-order accuracy the propagation of seismic waves in a 3D heterogeneous transversely isotropic medium with non-flat free surface. The method is an extension of the 2D method proposed by Appelo and Petersson (2009) to the 3D anisotropic case. The surface topography is introduced via mapping rectangular grids to curved grids. The accurate application of the free surface boundary conditions is done by using boundary-modified difference operators to discretize the mixed derivatives in the governing equations of the problem. Several numerical examples under different assumptions of free surface are given to highlight the complications of realistic seismic wave propagation in the vicinity of the earth surface. Synthetic seismograms and snapshots explain diffractions, scattering, multiple reflections, and converted waves provoked by the features of the free surface topography. The typical cuspidal triangles of the quasi-transverse (qS) mode also appear in the snapshots of the anisotropic medium.

The future directions of our research will include an extension of the schemes to the viscoelastic case. This will allow a realistic attenuation of the seismic waves due to the presence of a weathered layer to be included (2000).

#### **Author details**

Haiqiang Lan and Zhongjie Zhang *State Key Laboratory of Lithosphere Evolution, Institute of Geology and Geophysics, Chinese Academy of Sciences, China P.R.* 

### **Acknowledgement**

We are grateful to Jinlai Hao and Jinhai Zhang for their assistance and the facilities given in the course of this work. Fruitful discussions with Tao Xu, Kai Liang, Wei Zhang are also greatly appreciated. The Ministry of Science and Technology of China (SINOPROBE-02-02) supported this research.

#### **8. References**


Boore, D. M. (1972). A note on the effect of simple topography on seismic SH waves. *Bull. Seismol. Soc. Am.,* 62, 275-284.

Three-Dimensional Wavefield Simulation in

Heterogeneous Transversely Isotropic Medium with Irregular Free Surface 129

Komatitsch, D. & Tromp, J. (1999). Introduction to the spectral element method for three-

Komatitsch, D. & Tromp, J. (2002). Spectral-element simulations of global seismic wave

Lan, H. & Zhang, Z. (2011a). Comparative study of the free-surface boundary condition in two-dimensional finite-difference elastic wave field simulation. *Journal of Geophysics and* 

Lan, H. & Zhang, Z. (2011b). Three-Dimensional Wave-Field Simulation in Heterogeneous Transversely Isotropic Medium with Irregular Free Surface. *Bull. Seismol. Soc. Am.*,

Levander, A. R. (1990). Seismic scattering near the earth's surface. *Pure Appl. Geophys.,* 132,

Liu, E.; Crampin, S.; Queen, J. H. & Rizer, W. D. (1993). Velocity and attenuation anisotropy caused by microcracks and microfractures in a multiazimuth reverse VSP. *Can. J. Explor.* 

Lombard, B.; Piraux, J. ; Gélis, C. & Virieux, J. (2008). Free and smooth boundaries in 2-D finite-difference schemes for transient elastic waves. *Geophys. J. Int.,* 172, 252-261. Moczo, P.; Bystricky, E.; Kristek, J.; Carcione, J. & Bouchon, M. (1997). Hybrid modelling of P-SV seismic motion at inhomogeneous viscoelastic topographic structures. *Bull. Seism.* 

Nielsen, P.; If, F.; Berg, P. & Skovgaard, O. (1994). Using the pseudospectral technique on curved grids for 2D acoustic forward modelling. *Geophys. Prospect.,* 42, 321-342. Nilsson, S., Petersson, N. A.; Sjogreen, B. & Kreiss, H. O. (2007). Stable difference approximations for the elastic wave equation in second order formulation. *SIAM J.* 

Payton, R. G. (1983). *Elastic wave propagation in transversely isotropic medium*. Martinus Nijhoff

Rial, J. A.; Saltzman, N. G. & Ling, H. (1992). Earthquake-induced resonance in sedimentary

Robertsson, J. O. A. (1996). A numerical free-surface condition for elastic/viscoelastic finite-

Sanchez-Sesma, F. J. & Campillo, M. (1993). Topographic effects for incident P, SV and

Sanchez-Sesma, F. J., Ramos-Martinez, J. & Campillo, M. (2006). An indirect boundary element method applied to simulate the seismic response of alluvial valleys for incident

Schoenberg, M. & Muir, F. (1989). A calculus for finely layered anisotropic medium.

Tessmer, E.; Kosloff, D. & Behle, A. (1992). Elastic wave propagation simulation in the

Tessmer, E. & Kosloff, D. (1994). 3-D elastic modeling with surface topography by a

difference modeling in the presence of topography. *Geophysics,* 61, 1921.

P, S and Rayleigh waves. *Earthquake Eng. Struct. Dynam.,* 22, 279-295.

presence of surface topography. *Geophys. J. Int.*, 108, 621-632.

Chebychev spectral method. *Geophysics,* 59, 464-473.

dimensional seismic wave propagation. *Geophys. J. Int.,* 139, 806-822.

propagation-I. Validation. *Geophys. J. Int.,* 149, 390-412.

*Engineering*, 8, 275-286.

101(3), 1354–1370.

*Geophys.,* 29, 177-188.

*Soc. Am.,* 87, 1305-1323.

*Numer. Anal.,* 45, 1902-1936.

*Geophysics,* 54, 581-589.

basins. *American Scientist,* 80, 566-578.

Rayleigh waves. *Tectonophysics,* 218, 113-125.

21-47.

Publ.


Komatitsch, D. & Tromp, J. (1999). Introduction to the spectral element method for threedimensional seismic wave propagation. *Geophys. J. Int.,* 139, 806-822.

128 Earthquake Engineering

*Seismol. Soc. Am.,* 62, 275-284.

*porous medium*, Pergamon Amsterdam.

*Wave motion,* 3, 343-391.

2074.

888.

364-373.

*Earth Planet. Inter.,* 25, 297-356.

calculations. *Geophysics,* 53, 625-637.

mantle. *Geophys. J. Int.,* 43, 103-162.

surface topography. *Geophys. Prosp.,* 42, 371-390.

surface topography. *Geophysics,* 63, 613-622.

cracks. *Geophys. J. Int.,* 64, 133-150.

University of Denmark.

*Geophysics,* 53, 1045.

Boore, D. M. (1972). A note on the effect of simple topography on seismic SH waves. *Bull.* 

Bouchon, M.; Campillo, M. & Gaffet, S. (1989). A boundary integral equation-discrete wavenumber representation method to study wave propagation in multilayered

Bouchon, M.; Schultz, C. A. & Toksoz, M. N. (1995). A fast implementation of boundary integral equation methods to calculate the propagation of seismic waves in laterally

Campillo, M. & Bouchon, M. (1985). Synthetic SH seismograms in a laterally varying medium by the discrete wavenumber method. *Geophys. J. R. Astr. Soc.,* 83, 307-317. Carcione, J. M.; Kosloff, D.; Behle, A. & Seriani G. (1992). A spectral scheme for wave propagation simulation in 3-D elastic-anisotropic medium. *Geophysics,* 57, 1593-1607. Carcione, J. M. (2000). *Wave fields in real medium: Wave propagation in anisotropic, anelastic and* 

Crampin, S. (1981). A review of wave motion in anisotropic and cracked elastic-medium.

Dziewonski, A. M., and Anderson, D. L. (1981). Preliminary reference Earth model\* 1, *Phys.* 

Fornberg, B. (1988). The pseudo-spectral method: Accurate representation in elastic wave

Forsyth, D. W. (1975). The early structural evolution and anisotropy of the oceanic upper

Frankel, A. & Vidale, J. (1992). A three-dimensional simulation of seismic waves in the Santa Clara Valley, California, from a Loma Prieta aftershock. *Bull. Seism. Soc. Am.,* 82, 2045-

Gao, H. & Zhang, J. (2006). Parallel 3-D simulation of seismic wave propagation in heterogeneous anisotropic medium: a grid method approach. *Geophys. J. Int.,* 165, 875-

Galis, M.; Moczo, P. & Kristek J. (2008). A 3-D hybrid finite-difference—finite-element

Helbig, K. (1984). Anisotropy and dispersion in periodically layered medium. *Geophysics,* 49,

Hestholm, S. & Ruud, B. (1994). 2-D finite-difference elastic wave modeling including

Hestholm, S. & Ruud, B. (1998). 3-D finite-difference elastic wave modeling including

Hudson, J. A. (1981). Wave speeds and attenuation of elastic waves in material containing

Hvid, S. L. (1994). Three dimensional algebraic grid generation. *Ph.D. thesis,* Technical

Jih, R. S.; McLaughlin, K. L. & Der, Z. A. (1988). Free-boundary conditions of arbitrary polygonal topography in a two-dimensional explicit elastic finite-difference scheme.

viscoelastic modelling of seismic wave motion. *Geophys. J. Int.,* 175, 153-184.

medium having irregular interfaces. *Geophysics,* 54, 1134-1140.

varying layered medium. *Bull. Seismol. Soc. Am.,* 85, 1679-1687.


Thompson, J. F.; Warsi, Z. U. A. & Mastin, C. W. (1985). *Numerical grid generation: foundations and applications*. North-holland Amsterdam.

**Chapter 5** 

© 2012 Lee and Chen, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 Lee and Chen, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

**Full-Wave Ground Motion** 

Additional information is available at the end of the chapter

En-Jui Lee and Po Chen

http://dx.doi.org/10.5772/50114

**1. Introduction** 

Los Angles region[5].

**Forecast for Southern California** 

The damages and loses caused by earthquakes are increased as urbanization increased in past decades. However, scientists currently cannot predict the time, location and magnitude of an earthquake accurately. Currently, the earthquake predictions are not yet reliable, so long-term probabilistic earthquake hazard analysis and rapid post-earthquake early

For long-term ground motion forecasts, seismic hazard maps are widely used for estimating probabilities of ground motion exceeding certain amount in different areas in 50 years. In addition, the ground motion estimations in seismic hazard maps are useful for different applications, including, for instance, building codes, insurance rates, land-use policies and education of earthquake response. In United States, the U.S. Geological Survey (USGS) periodically updates the National Seismic Hazard maps that provide a 50 years ground motion estimations for United States[4, 3]. To consider effects form different aspects, the National Seismic Hazard maps incorporate both geological and geophysical information in ground motion estimations[3]. Recent advances in computational seismology allow us to simulate wave propagation in complex velocity structure models and then open the probability of physics-based long-term ground motion estimations. The Southern California Earthquake Center (SCEC) has developed a methodology which considers both source and structure effects in ground motion simulations for long-term ground motion estimations in

The earthquake early warning systems are designed for short-term ground motion forecasts. The idea of earthquake early warning systems is based on the transmission speed of electromagnetic signal is much faster than the propagation speed of seismic shear-waves and surface waves that usually generate strong ground motion [2]. The Earthquake Alarms Systems, ElarmS, is the earthquake early warning systems designed for California region [1,

warning are two alternative solutions to reduce potential earthquake damages [1-3].


**Chapter 5** 
