# A note on the dynamic and static displacements from a point source in multilayered media

Lupei Zhu, Luis Rivera

► To cite this version:

Lupei Zhu, Luis Rivera. A note on the dynamic and static displacements from a point source in multilayered media. Geophysical Journal International, 2002, 148 (3), pp.619-627. 10.1046/j.1365-246X.2002.01610.x. hal-03200897

HAL Id: hal-03200897

https://hal.science/hal-03200897

Submitted on 27 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A note on the dynamic and static displacements from a point source in multilayered media

Lupei Zhu<fnref n="1" />,\* and Luis A. Rivera<fnref n="2" />,<fnref n="3" />

Accepted 2001 September 16. Received 2001 May 7; in original form 2000 July 13

## SUMMARY

A simple and unified approach is presented to solve both the elasto-dynamic and elasto-static problems of point sources in a multi-layered half-space by using the Thompson-Haskell propagator matrix technique. It is shown that the apparent incompatibility between the two is associated with the degeneracy of the dynamic problem when $\omega = 0$ and both can be handled uniformly using the Jordan canonical forms of matrices. We re-derive the propagator matrices for both the dynamic and static cases. We then show that the dynamic propagator matrix and the solution converge to their static counterparts as $\omega \to 0$. Satisfactory static deformation can be obtained numerically using the dynamic solution at near-zero frequency.

**Key words:** deformation, elastic wave theory, seismic wave propagation.

## 1 INTRODUCTION

Historically, the problems of dynamic displacement and static deformation generated by various sources in a layered half-space were treated separately. Thomson (1950) and Haskell (1953) developed a matrix formalism for studying elastic wave propagation in a multi-layered medium. The solution of surface displacement due to a point source in the medium was given by Haskell (1963, 1964) and Harkrider (1964) using the propagator matrix technique. Kennett (1974) partitioned the Thomson-Haskell propagator matrix into submatrices and introduced them as the generalized reflection-transmission coefficient matrices. The elasto-static problem of a dislocation source in a uniform half-space was first solved by Steketee (1958) using the image source technique. He found that a total of 6 sets of the Green's function are needed to describe an arbitrary dislocation and gave a compact solution for the vertical strike-slip fault case. Maruyama (1964) presented solutions for other cases. Ben-Menahem & Singh (1968) proposed a general approach to handle both the dynamic and static problems using the Hansen eigenvector expansion. It was used by Singh (1970) to formalize the static problem in the form of the propagator matrix and obtain the solution for multi-layered media in the same fashion as in the dynamic problem. Similar solution was also given by Sato (1971).

In theory, the dynamic solution should converge to the static solution as the frequency, $\omega$, goes to zero. However, this intrinsic connection between the two were often obscured by different conventions used by different authors and, sometimes, lengthy algebra. For example, in Haskell's definition of the displacement-stress vector (Haskell 1964), the stress component was normalized by $\omega^2$, which introduced an apparent singularity at $\omega = 0$ in the source terms and the propagator matrix. Some of the mathematical tricks in dealing with one case were given in an *ad hoc* way that they are difficult to apply to other cases. In this paper, we present a simple and unified approach to solve both the dynamic and static problems and show that the static solution is indeed the $\omega = 0$ component of the dynamic solution. We first re-derive the propagator matrices for both dynamic and static cases separately using the Jordan canonical forms of matrices. The surface displacement solution of a point source in multi-layered half-space is obtained. We then show that the dynamic propagator matrix and the solution converge to their static counterparts as $\omega \to 0$.

## 2 THEORY

### 2.1 Displacement-stress in source-free homogeneous medium

We set up a cylindrical coordinate system ($e_r, e_\theta, e_z$), with $e_z$ pointing upward. The displacement in a vertically heterogeneous medium can be expanded in terms of three orthogonal vectors (e.g. Takeuchi & Saito 1972):

<footnote n="1">¹Department of Earth and Atmospheric Sciences, Saint Louis University, St Louis, MO 63103, USA. E-mail: lupei@eas.slu.edu</footnote>

<footnote n="2">²Ecole et Observatoire des Sciences de la Terre; 5, rue Rene Descartes; 67084 Strasbourg, France</footnote>

<footnote n="3">³Seismological Laboratory, Caltech, Pasadena, CA 91125, USA</footnote>

<footnote>\*Contact: Tel +1 314 977 3118; Fax +1 314 977 3117.</footnote>
$$ \mathbf{u}(r, \theta, z, t) = \frac{1}{2\pi} \sum_{m=0, \pm 1, \dots} \int e^{-i\omega t} d\omega \int_0^\infty k \, dk (U_z \mathbf{R}_m^k + U_r \mathbf{S}_m^k + U_\theta \mathbf{T}_m^k), \quad (1) $$

$\mathbf{R}_m^k, \mathbf{S}_m^k, \mathbf{T}_m^k$ are called the surface vector harmonics. Similar expansion can be done to the traction on the horizontal plane

$$ \sigma(r, \theta, z, t) = \frac{1}{2\pi} \sum_{m=0, \pm 1, \dots} \int e^{-i\omega t} d\omega \int_0^\infty k^2 dk (T_z \mathbf{R}_m^k + T_r \mathbf{S}_m^k + T_\theta \mathbf{T}_m^k). \quad (2) $$

Note that a factor $k$ is drawn from the $T_z$, $T_r$ and $T_\theta$ to simplify the later derived matrices. Eqs (1) and (2) separate the $z$-variation of the displacement and stress from the $(r, \theta)$ variations. Under this expansion, the second-order differential equation of motion

$$ (\lambda + 2\mu) \nabla (\nabla \cdot \mathbf{u}) - \mu \nabla \times \nabla \times \mathbf{u} + \rho \omega^2 \mathbf{u} = 0 \quad (3) $$

is reduced to a set of first-order ordinary differential equations

$$ \frac{d}{dz} \begin{bmatrix} U_r \\ U_z \\ T_z \\ T_r \\ U_\theta \\ T_\theta \end{bmatrix} = k \begin{bmatrix} 0 & -1 & 0 & \frac{1}{\mu} & 0 & 0 \\ 1 - 2\xi & 0 & \frac{\xi}{\mu} & 0 & 0 & 0 \\ 0 & -\rho\left(\frac{\omega}{k}\right)^2 & 0 & 1 & 0 & 0 \\ 4\mu\xi_1 - \rho\left(\frac{\omega}{k}\right)^2 & 0 & 2\xi - 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{\mu} \\ 0 & 0 & 0 & 0 & \mu - \rho\left(\frac{\omega}{k}\right)^2 & 0 \end{bmatrix} \begin{bmatrix} U_r \\ U_z \\ T_z \\ T_r \\ U_\theta \\ T_\theta \end{bmatrix}, \quad (4) $$

or, in the vector-matrix form:

$$ \frac{d\mathbf{b}(z)}{dz} = \mathbf{M}\mathbf{b}(z). \quad (5) $$

Here, $\rho$ is the density; $\xi = \mu/(\lambda + 2\mu)$; $\xi_1 = 1 - \xi$; $\lambda, \mu$ are the Lame constants. The vector $\mathbf{b}$ is often called the displacement-stress vector. Note that $\mathbf{M}$ can be partitioned into a $4 \times 4$ submatrix describing the motion in the $(z, r)$ plane and a $2 \times 2$ submatrix for the motion in the $\theta$-direction. They are often referred to as the $P - SV$ system and the $SH$ system. We will concentrate on the $P - SV$ problem. The corresponding solution for the $SH$ problem can be found in Appendix A.

For a homogeneous medium, $\mathbf{M}$ is constant. In this case, the general solution of (5) is

$$ \mathbf{b}(z) = e^{\mathbf{z}\mathbf{M}}\mathbf{b}_0. \quad (6) $$

To calculate the matrix exponential, we use the Jordan decomposition of $\mathbf{M}$ (e.g. Turnbull & Aitken 1952; Ganter 1960):

$$ \mathbf{M} = \mathbf{E}\mathbf{J}\mathbf{E}^{-1}, \quad (7) $$

where $\mathbf{E}$ is a similarity matrix and $\mathbf{J}$ is the Jordan canonical form of $\mathbf{M}$. Using (7) and the definition of matrix exponential, we have

$$ e^{\mathbf{z}\mathbf{M}} = \mathbf{E}e^{-\mathbf{z}\mathbf{J}}\mathbf{E}^{-1}. \quad (8) $$

From (6) and (8), the displacement-stress at any $z$ can be expressed as

$$ \mathbf{b}(z) = \mathbf{E}\Lambda(z)\mathbf{w}, \quad (9) $$

where

$$ \Lambda(z) = e^{\mathbf{z}\mathbf{J}}, \quad (10) $$

$\mathbf{w}$ is a constant vector that is to be determined by boundary conditions.

If the matrix $\mathbf{M}$ has a complete set of independent eigenvectors, the Jordan decomposition is reduced to the eigenvalue decomposition. $\mathbf{J}$ is simply a diagonal matrix with the eigenvalues as the diagonal elements and the columns of $\mathbf{E}$ are the eigenvectors. This is the case for the elasto-dynamic problem with $\omega \neq 0$. For the $P - SV$ system, the $4 \times 4 \mathbf{M}$ has 4 eigenvalues and

$$ \mathbf{J} = \begin{bmatrix} -\nu_{\alpha} & 0 & 0 & 0 \\ 0 & -\nu_{\beta} & 0 & 0 \\ 0 & 0 & \nu_{\alpha} & 0 \\ 0 & 0 & 0 & \nu_{\beta} \end{bmatrix}, \quad (11) $$

where $\nu_\alpha = \sqrt{k^2 - (\frac{\omega}{\alpha})^2}$ and $\nu_\beta = \sqrt{k^2 - (\frac{\omega}{\beta})^2}$. Here, $\alpha, \beta$ are the compressional and shear velocities of the medium. Correspondingly,
$$ \mathbf{E} = \begin{bmatrix} -1 & -\frac{v_{\beta}}{k} & 1 & \frac{v_{\beta}}{k} \\ \frac{v_{\alpha}}{k} & 1 & \frac{v_{\alpha}}{k} & 1 \\ -2\mu\gamma_1 & -2\mu\frac{v_{\beta}}{k} & 2\mu\gamma_1 & 2\mu\frac{v_{\beta}}{k} \\ 2\mu\frac{v_{\alpha}}{k} & 2\mu\gamma_1 & 2\mu\frac{v_{\alpha}}{k} & 2\mu\gamma_1 \end{bmatrix}, \qquad (12) $$

$$ \Lambda(z) = \begin{bmatrix} e^{-v_{\alpha}z} & 0 & 0 & 0 \\ 0 & e^{-v_{\beta}z} & 0 & 0 \\ 0 & 0 & e^{v_{\alpha}z} & 0 \\ 0 & 0 & 0 & e^{v_{\beta}z} \end{bmatrix}, \qquad (13) $$

$$ \mathbf{E}^{-1} = \frac{\gamma}{2} \begin{bmatrix} -1 & -\gamma_1 \frac{k}{v_{\alpha}} & \frac{1}{2\mu} & \frac{k}{2\mu v_{\alpha}} \\ \gamma_1 \frac{k}{v_{\beta}} & 1 & -\frac{k}{2\mu v_{\beta}} & -\frac{1}{2\mu} \\ 1 & -\gamma_1 \frac{k}{v_{\alpha}} & -\frac{1}{2\mu} & \frac{k}{2\mu v_{\alpha}} \\ -\gamma_1 \frac{k}{v_{\beta}} & 1 & \frac{k}{2\mu v_{\beta}} & -\frac{1}{2\mu} \end{bmatrix}, \qquad (14) $$

where $\gamma = 2k^2\beta^2/\omega^2$, $\gamma_1 = 1 - 1/\gamma$.

## 2.2 Surface displacement of layered half-space from an embedded point source

In case of a multi-layered half-space (see Fig. 1), the displacement-stress vectors at the top ($z = z_{n-1}$) and the bottom ($z = z_n$) of a layer can be related as

$$ \mathbf{b}(z_n) = \mathbf{a}_n \mathbf{b}(z_{n-1}), \qquad (15) $$

where

$$ \mathbf{a}_n = \mathbf{E}_n \Lambda_n (z_n - z_{n-1}) \mathbf{E}_n^{-1}, \qquad (16) $$

is often called the Thomson-Haskell propagator matrix.

For the P-SV system

$$ \mathbf{a}_n = \gamma \begin{bmatrix} C_\alpha - \gamma_1 C_\beta & \gamma_1 Y_\alpha - X_\beta & \frac{C_\beta - C_\alpha}{2\mu} & \frac{X_\beta - Y_\alpha}{2\mu} \\ \gamma_1 Y_\beta - X_\alpha & C_\beta - \gamma_1 C_\alpha & \frac{X_\alpha - Y_\beta}{2\mu} & \frac{C_\alpha - C_\beta}{2\mu} \\ 2\mu\gamma_1(C_\alpha - C_\beta) & 2\mu(\gamma_1^2 Y_\alpha - X_\beta) & C_\beta - \gamma_1 C_\alpha & X_\beta - \gamma_1 Y_\alpha \\ 2\mu(\gamma_1^2 Y_\beta - X_\alpha) & 2\mu\gamma_1(C_\beta - C_\alpha) & X_\alpha - \gamma_1 Y_\beta & C_\alpha - \gamma_1 C_\beta \end{bmatrix} \qquad (17) $$

<figure><img src="image_2.png" /><figcaption>**Figure 1.** The layered half-space consists of N layers over a half-space at the bottom. The source is located at a depth of h between layer m and m + 1 with identical elastic properties.</figcaption></figure>
where $C_\alpha = \cosh(v_\alpha d)$, $X_\alpha = v_\alpha \sinh(v_\alpha d)/k$, and $Y_\alpha = k \sinh(v_\alpha d)/v_\alpha$, and similarly for $C_\beta$, $X_\beta$, and $Y_\beta$. $d = z_{n-1} - z_n$ is the thickness of the layer. Our $a_n$ is different from the original one given by Haskell (1964) for the traction related terms. This stems from the difference of our definition of the displacement-stress vector from the Haskell's in which he multiplied the traction by $\omega^2$. As shown later, the Haskell's definition introduces the apparent $\omega$-dependence of source terms and causes difficulty to unify the elasto-dynamic solution and the elasto-static solution.

Using the traction free condition at the surface and the radiation condition in the half-space, the displacement at the surface is obtained as (e.g. Haskell 1964)

$$ \begin{bmatrix} U_r \\ U_z \end{bmatrix} = \frac{1}{\mathbf{R}|_{12}^{12}} \begin{bmatrix} R_{22} & -R_{12} \\ -R_{21} & R_{11} \end{bmatrix} \begin{bmatrix} L_{1i} s_i \\ L_{2i} s_i \end{bmatrix}, \qquad (18) $$

where $\mathbf{s}$ is the displacement-stress jump produced by the source (see Appendix B),

$$ \mathbf{L} = \mathbf{E}_{N+1}^{-1} \mathbf{a}_N \cdots \mathbf{a}_{N+1}, \qquad (19) $$

$$ \mathbf{R} = \mathbf{E}_{N+1}^{-1} \mathbf{a}_N \cdots \mathbf{a}_1. \qquad (20) $$

Here $\mathbf{R}|_{12}^{12}$ is a subdeterminant of $\mathbf{R}$ defined by:

$$ \mathbf{R}_{ij}^{kl} = R_{ik} R_{jl} - R_{il} R_{jk}. \qquad (21) $$

## 2.3 Static deformation of layered half-space

In the static case, the same first-order differential equation (4) still holds except that the $\rho(\omega_k)^2$ term vanishes. The $4 \times 4 \mathbf{M}$ of the $P-SV$ system is degenerate with only two distinct eigenvalues $\pm k$ and two independent eigenvectors, see (12). The canonical form of $\mathbf{M}$ is

$$ \mathbf{J} = \begin{bmatrix} -k & 1 & 0 & 0 \\ 0 & -k & 0 & 0 \\ 0 & 0 & k & 1 \\ 0 & 0 & 0 & k \end{bmatrix}, \qquad (22) $$

and the exponential matrix of $z\mathbf{J}$ is:

$$ e^{z\mathbf{J}} = \begin{bmatrix} e^{-kz} & z e^{-kz} & 0 & 0 \\ 0 & e^{-kz} & 0 & 0 \\ 0 & 0 & e^{kz} & z e^{kz} \\ 0 & 0 & 0 & e^{kz} \end{bmatrix}. \qquad (23) $$

So, the corresponding matrices to construct the displacement-stress as in (9) for the static case are:

$$ \mathbf{E} = \begin{bmatrix} -1 & \xi & 1 & \xi \\ 1 & 1 & 1 & -1 \\ -2\mu & -2\mu\xi_1 & 2\mu & -2\mu\xi_1 \\ 2\mu & 0 & 2\mu & 0 \end{bmatrix}, \qquad (24) $$

$$ \Lambda(z) = \begin{bmatrix} e^{-kz} & \xi_1 k z e^{-kz} & 0 & 0 \\ 0 & e^{-kz} & 0 & 0 \\ 0 & 0 & e^{kz} & \xi_1 k z e^{kz} \\ 0 & 0 & 0 & e^{kz} \end{bmatrix}, \qquad (25) $$

$$ \mathbf{E}^{-1} = \frac{1}{2} \begin{bmatrix} -\xi_1 & 0 & -\frac{\xi}{2\mu} & \frac{1}{2\mu} \\ 1 & 1 & -\frac{1}{2\mu} & -\frac{1}{2\mu} \\ \xi_1 & 0 & \frac{\xi}{2\mu} & \frac{1}{2\mu} \\ 1 & -1 & -\frac{1}{2\mu} & \frac{1}{2\mu} \end{bmatrix}. \qquad (26) $$

Here, we deliberately multiply the off-diagonal elements of $\Lambda(z)$ by $\xi_1 k$ to simplify $\mathbf{E}$. The propagator matrix for the static case is
$$ \mathbf{a}_n = \begin{bmatrix} C + Sx & Cx + S\xi & -\frac{Sx}{2\mu} & -\frac{Cx + S(1+\xi)}{2\mu} \\ S\xi - Cx & C - Sx & \frac{Cx - S(1+\xi)}{2\mu} & \frac{Sx}{2\mu} \\ 2\mu Sx & 2\mu(Cx - S\xi_1) & C - Sx & -Cx - S\xi \\ -2\mu(Cx + S\xi_1) & -2\mu Sx & Cx - S\xi & C + Sx \end{bmatrix} \quad (27) $$

where $x = dk\xi_1$, $C = \cosh(kd)$, and $S = \sinh(kd)$. Eq. (27) is the same as obtained by Singh (1970) with the only exception that the third and fourth rows (and also the third and fourth columns) are swapped due to the difference in **b**.

By following the same procedure as in the dynamic case, one can obtain the solution (18) of surface static deformation of a layered half-space. The only difference is that the $\mathbf{E}^{-1}$ and $\mathbf{a}_n$ matrices are replaced by their counterparts for the static case.

## 2.4 Dynamic solution at zero frequency

In the previous section, we derived the propagator matrix and displacement for both the dynamic and static cases. It is interesting to recognize the similarity between (17) and (27). In fact, it can be shown that at zero frequency (17) becomes (27). Note that all elements of the dynamic propagator matrix are in the form of $\gamma(F_\alpha - F_\beta)$, where $\gamma \sim \frac{1}{\omega^2}$ and $F_{\alpha,\beta}$ are functions of $\omega$ through $v_\alpha = \sqrt{k^2 - (\frac{\omega}{\alpha})^2}$ and $v_\beta = \sqrt{k^2 - (\frac{\omega}{\beta})^2}$. They become 0/0 as $\omega \to 0$. Using the L'Hôpital rule, one can get

$$ \lim_{\omega \to 0} \gamma(C_\alpha - C_\beta) = xS, \quad (28) $$

$$ \lim_{\omega \to 0} \gamma(Y_\alpha - X_\beta) = xC - \xi_1 S, \quad (29) $$

$$ \lim_{\omega \to 0} \gamma(X_\alpha - Y_\beta) = xC + \xi_1 S. \quad (30) $$

From these relations it is clear that the dynamic propagator matrix becomes the static propagator matrix at zero frequency. This is expected because that both of them are defined as the connection between the displacement-stresses at two depths. What is less obvious is that the displacement solution (18) for the dynamic case also becomes identical to the static solution at $\omega = 0$, even though the matrices $\mathbf{E}_N^{-1}$, and consequentially the **R** and **L**, are different for the two cases.

To prove this, we use the radial displacement solution as an example:

$$ U_r = \frac{(R_{22}L_{1i} - R_{12}L_{2i})s_i}{R_{11}R_{22} - R_{21}R_{12}}. \quad (31) $$

From (17)–(20), we can write $R_{ij} = e_{ik}r_{kj}$ and $L_{ij} = e_{ik}l_{kj}$, where $e_{ik}$ is the matrix elements of $\mathbf{E}^{-1}$. Matrices ($r_{kj}$) and ($l_{kj}$) are multiplications of propagator matrices, so they are the same for both the dynamic and static cases at $\omega = 0$. By rearranging order, eq. (31) becomes

$$ U_r = \frac{(e_{11}e_{2k} - e_{1k}e_{2l})r_{l2}l_{ki}s_i}{(e_{1l}e_{2k} - e_{1k}e_{2l})r_{k1}r_{l2}} = \frac{\mathbf{E}^{-1}|_{lk}^{12}r_{l2}l_{ki}s_i}{\mathbf{E}^{-1}|_{lk}^{12}r_{k1}r_{l2}}. \quad (32) $$

The non-zero subdeterminants of $\mathbf{E}^{-1}$ for the dynamic case are

$$ \mathbf{E}^{-1}|_{lk}^{12} = \frac{\gamma k^2}{8\mu v_\beta v_\alpha} \left( 2\mu(\delta - \gamma_1), \frac{v_\alpha}{k}, -\delta, \delta, -\frac{v_\beta}{k}, \frac{\delta + 1}{2\mu} \right), \quad (33) $$

where $lk = \{12, 13, 14, 23, 24, 34\}$, $\delta = \gamma(1 - v_\alpha v_\beta / k^2) - 1$. For the static case,

$$ \mathbf{E}^{-1}|_{lk}^{12} = \frac{1}{8\mu} \left( 2\mu(\xi - 1), 1, -\xi, \xi, -1, \frac{\xi + 1}{2\mu} \right). \quad (34) $$

One can show that

$$ \lim_{\omega \to 0} \delta = \xi, \quad (35) $$

so eqs (33) and (35) become identical at $\omega = 0$, except for a constant which is cancelled among the numerator and denominator in (32). Therefore, the static solution is nothing but a special case in the dynamic solution at $\omega = 0$.

Since the dynamic Thomson-Haskell matrix elements become 0/0 at zero frequency, numerical computation of the static deformation using the dynamic solution needs a small frequency to offset the evaluation from exact $\omega = 0$. An alternative is to design the code that switches between static and dynamic propagator matrices based on the frequency. However, numerical tests show it is not necessary for general purposes and digital computers can handle near-zero frequency accurately enough. Fig. 2 shows the static displacement from an explosion source in a homogeneous half-space. We use the dynamic solution with a small imaginary frequency 0.01i Hz to approximate the zero frequency and compare the result with the exact static solution. The relative errors are less than 1 per cent up to a distance of 50 km.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Comparison of numerically computed static deformation (dots) using the dynamic solution at near-zero frequency with the analytical solutions (solid lines) for an 5 km deep explosion source of $10^{25}$ dyn·cm in a homogeneous half-space ($Vp = 6.3$ km s$^{-1}$, $Vs = 3.5$ km s$^{-1}$, $\rho = 2.8$ g cm$^{-3}$).</figcaption></figure>

### 3 DISCUSSION AND CONCLUSIONS

When dealing with wave propagation in a layered medium, it is useful to write the wave field as the superposition of three terms: $\nabla\Psi_1$, $\nabla \times \nabla \times (\Psi_2 e_z)$, and $\nabla \times (\Psi_3 e_z)$ standing for the *P*, *SV* and *SH* components of the field respectively. The three scalar functions (potentials) obey the scalar wave equations (with propagation velocities of $\alpha, \beta$ and $\beta$, respectively). It was early realized however that this decomposition is inadequate for the $\omega \to 0$ limit (Hansen 1936; Morse & Feshbach 1953). This problem is also present in other fields such as electromagnetism. It is particularly clear in elasticity where the ‘P’ and the ‘SV’ fields become linearly dependent and are no longer sufficient to represent a general solution (along with the *SH* field) for the static problem (Takeuchi 1959). To solve the elasto-static problem, different authors modified the basic fields by replacing $\nabla \times \nabla \times (\Psi_2 e_z)$ with $\Psi_2 e_z + \frac{1-\xi}{1+\xi} z \nabla \Psi_2$ (Takeuchi 1959; Ben-Menahem & Singh 1968; Singh 1970), or changed the potentials (Sato 1971). Either way, an extraneous term linearly dependent on $z$ was introduced.

It is shown in this paper that this apparent incompatibility is associated with the deficiency of **M** when $\omega=0$ and can be handled uniformly using the Jordan decomposition. The non-diagonal terms appearing in the static $\Lambda(z)$ (see 23) correspond to the *ad hoc* linear $z$-dependent field previously introduced to solve the elasto-static problem. With this clarification, the static solution finds its place as a natural component of the dynamic solution. The operation of reducing a matrix into its Jordan canonical form is unstable (Arnold 1981) in the sense that even if the entries of the original matrix vary continuously, **J** and the associated similarity matrix **E** can vary discontinuously, as demonstrated by eqs (11) and (22) in the transition from the elasto-dynamic case to the elasto-static case.

It should be pointed out that the degeneracy of **M** discussed above occurs for the isotropic medium and will disappear for general anisotropic medium. For a special case of the transversely isotropic layered medium, one can derive the same first-order differential equation (5) (Takeuchi & Saito 1972) where **M** has a complete set of independent eigenvectors at zero frequency. Pan (1989) gave the static propagator matrix and the surface response for a transversely isotropic layered medium.

In conclusion, we found that the elasto-dynamic and elasto-static problems for a vertically heterogeneous medium can be treated in a unified way using the Jordan decomposition. We have derived separately the propagator matrices and the displacement responses of the medium for the both cases and showed that the dynamic solution converges to the static solution as $\omega \to 0$. Satisfactory static deformation can be obtained numerically using the dynamic solution at near-zero frequency.

### ACKNOWLEDGMENTS

This work was inspired by related discussions during the daily coffee-break at the seismological laboratory, Caltech. The authors like to thank H. Kanamori, D.L. Anderson, and D.V. Helmberger for their encouragement. Comments by R.B. Herrmann, S.N. Ward, E. Brodsky, M. Simons and J. Ritsema help us to improve the manuscript greatly. All the mathematical derivations were carried out using the Wolfram
Research's Mathematica package. LZ was supported partially by a USGS NEHRP grant 00HQGR0007 and a NSF grant EAR-0001016. This is Division of Geological and Planetary Sciences, Caltech contribution 8817.

## REFERENCES

Aki, K. & Richards, P.G., 1980. *Quantitative Seismology, Theory and Methods*, Freeman and Company.

Arnold, V.I., 1981. On matrices depending on parameters, in *Singularity Theory*, vol. 53 of Lecture Note Series, pp. 46–60, Cambridge University Press.

Ben-Menahem, A. & Singh, S.J., 1968. Multipolar elastic fields in a layered half-space, *Bull. seism. Soc. Am.*, **58**, 1519–1572.

Gant matcher, F.R., 1960. *Matrix Theory*, Chelsea Publishing Company, New York, NY.

Hansen, W.W., 1936. A new type of expansion in radiation problems, *Phys. Rev.*, **47**, 139–143.

Harkrider, D.G., 1964. Surface waves in multilayered elastic media. i. Rayleigh and Love from buried sources in a multilayered elastic half-space, *Bull. seism. Soc. Am.*, **54**, 627–679.

Haskell, N.A., 1953. The dispersion of surface waves on multilayered media, *Bull. seism. Soc. Am.*, **43**, 17–34.

Haskell, N.A., 1963. Radiation pattern of Rayleigh waves from a fault of arbitrary dip and direction of motion in a homogeneous medium, *Bull. seism. Soc. Am.*, **51**, 495–502.

Haskell, N.A., 1964. Radiation pattern of surface waves from point sources in a multi-layered medium, *Bull. seism. Soc. Am.*, **54**, 377–393.

Helmberger, D.V., 1983. Theory and application of synthetic seismograms, in *Earthquakes: Observation, Theory and Interpretation*, pp. 174–222, Soc. Italiana di Fisica, Bolgna, Italy.

Kennett, B.L.N., 1974. Reflections, rays, and reverberations, *Bull. seism. Soc. Am.*, **64**, 1685–1696.

Maruyama, T., 1964. Statical elastic dislocations in a infinite and semi-infinite medium, *Bull. Earthq. Res. Inst., Tokyo Univ.*, **42**, 289–368.

Morse, M.M. & Feshbach, H., 1953. *Methods of Theoretical Physics*, McGraw-Hill Book Company, New York.

Pan, E., 1989. Static response of a transversely isotropic and layered half-space to general surface loads, *Phys. Earth planet. Inter.*, **54**, 353–363.

Sato, R., 1971. Crustal deformation due to dislocation in a multi-layered medium, *J. Phys. Earth*, **19**, 32–46.

Singh, S.J., 1970. Static deformation of a multilayered half-space by internal sources, *J. geophys. Res.*, **75**, 3257–3263.

Steketee, J.A., 1958. Volterra's dislocations in a semi-infinite elastic medium, *Can. J. Phys.*, **36**, 192–205.

Takeuchi, H., 1959. General solutions of equations of some geophysical importance, *Bull. seism. Soc. Am.*, **49**, 273–283.

Takeuchi, H. & Saito, M., 1972. Seismic surface waves, in *Methods in Computational Physics*, Vol. 11, pp. 217–295, ed. Bolt, B.A., Academic Press, New York.

Thomson, W.T., 1950. Transmission of elastic waves through a stratified solid medium, *J. Appl. Phys.*, **21**, 89–93.

Turnbull, H.W. & Aitken, A.C., 1952. *An introduction to the theory of canonical matrices*, Blackie and Son, Ltd., London and Glasgow.

Wang, C.Y. & Herrmann, R.B., 1980. A numerical study of P, SV, and SH-wave generation in a plane layered medium, *Bull. seism. Soc. Am.*, **70**, 1015–1036.

## APPENDIX A: SH PROBLEM

For the dynamic SH problem, the $2 \times 2$ M is

$$ M = k \begin{bmatrix} 0 & \frac{1}{\mu} \\ \mu - \rho \left(\frac{\omega}{k}\right)^2 & 0 \end{bmatrix}. \qquad (A1) $$

It has two eigenvalues $\pm\nu_\beta$, which yields

$$ E = \begin{bmatrix} -1 & 1 \\ \mu \frac{\nu_\beta}{k} & \mu \frac{\nu_\beta}{k} \end{bmatrix}, \qquad (A2) $$

$$ \Lambda(z) = \begin{bmatrix} e^{-\nu_\beta z} & 0 \\ 0 & e^{\nu_\beta z} \end{bmatrix}, \qquad (A3) $$

$$ E^{-1} = \frac{1}{2} \begin{bmatrix} -1 & \frac{k}{\mu\nu_\beta} \\ 1 & \frac{k}{\mu\nu_\beta} \end{bmatrix}. \qquad (A4) $$

The Thomson-Haskell propagator matrix is

$$ a_n = \begin{bmatrix} C_\beta & -\frac{\nu_\beta}{\mu} \\ -\mu X_\beta & C_\beta \end{bmatrix}. \qquad (A5) $$

The solution of the surface displacement kernel for multi-layered half-space

$$ U_\theta = \frac{L_{1i}s_i}{R_{11}}, \qquad (A6) $$

R and L are $2 \times 2$ matrices in the forms as given in (19)–(20).

In the static case, there is no degeneration in M. All corresponding matrices can be obtained from above by setting $\omega = 0$:

$$ E = \begin{bmatrix} -1 & 1 \\ \mu & \mu \end{bmatrix}, \qquad (A7) $$
$$ \Lambda(z) = \begin{bmatrix} e^{-kz} & 0 \\ 0 & e^{kz} \end{bmatrix}, \qquad (A8) $$

$$ \mathbf{E}^{-1} = \frac{1}{2} \begin{bmatrix} -1 & \frac{1}{\mu} \\ 1 & \frac{1}{\mu} \end{bmatrix}, \qquad (A9) $$

and

$$ \mathbf{a}_n = \begin{bmatrix} C & -\frac{S}{\mu} \\ -\mu S & C \end{bmatrix}. \qquad (A10) $$

# APPENDIX B: SOURCE TERMS AND HORIZONTAL RADIATION PATTERN

The displacement-stress discontinuities to represent the source can be found by expanding the solution of the whole-space problem with the cylindrical spherical harmonics (1) and (2) (Haskell 1953; Takeuchi & Saito 1972). Takeuchi & Saito (1972) have listed the displacement-stress discontinuities produced by various sources. Since our definition of s is slightly different from theirs, we will give below the non-zero terms for several types of sources often encountered in seismology.

## B1 Explosion source

Only $m = 0$ term exists for the isotropic source

$$ \mathbf{s}^0 = (0, \xi/\mu, 0, 2\xi, 0, 0)^T. \qquad (B1) $$

## B2 Single force

None-zero terms exist for $m = 0, \pm 1$. We use the symmetry between $m = -1$ and $m = 1$ terms and factor out the common source geometry independent term. This reduces the number of source vectors s from 3 to 2:

$$ \mathbf{s}^0 = \frac{1}{k}(0, 0, -1, 0, 0, 0)^T, \qquad (B2) $$

$$ \mathbf{s}^1 = \frac{1}{k}(0, 0, 0, -1, 0, 1)^T. \qquad (B3) $$

They will produce five components of ground displacement, $u_z^0, u_r^0, u_z^1, u_r^1, u_\theta^1$ ($u_\theta^0$ is always zero). The actual displacement is obtained by adding the force orientation when summing over azimuthal modes $m$. By re-arranging terms, the summation can be expressed as

$$ u_z = \cos \phi \cos \delta u_z^1 - \sin \delta u_z^0, \qquad (B4) $$

$$ u_r = \cos \phi \cos \delta u_r^1 - \sin \delta u_r^0, \qquad (B5) $$

$$ u_\theta = -\sin \phi \cos \delta u_\theta^1, \qquad (B6) $$

where $\delta$ is the dip angle of the force, measured from the horizontal plane; $\phi$ is the azimuth of the station, measured clockwise from the direction of the force.

## B3 Double-couple without torque

Similar to the single force, the five none-zero source vectors ($m = 0, \pm 1, \pm 2$) can be reduced to three:

$$ \mathbf{s}^0 = (0, 2\xi/\mu, 0, 4\xi - 3, 0, 0)^T, \qquad (B7) $$

$$ \mathbf{s}^1 = (1/\mu, 0, 0, 0, -1/\mu, 0)^T, \qquad (B8) $$

$$ \mathbf{s}^2 = (0, 0, 0, 1, 0, -1)^T, \qquad (B9) $$

The displacement for arbitrary double-couple is

$$ u_z = \frac{1}{2} \sin 2\delta \sin \lambda u_z^0 - (\sin \phi \cos 2\delta \sin \lambda - \cos \phi \cos \delta \cos \lambda) u_z^1 - \left( \sin 2\phi \sin \delta \cos \lambda + \frac{1}{2} \cos 2\phi \sin 2\delta \sin \lambda \right) u_z^2, \qquad (B10) $$

$$ u_r = \frac{1}{2} \sin 2\delta \sin \lambda u_r^0 - (\sin \phi \cos 2\delta \sin \lambda - \cos \phi \cos \delta \cos \lambda) u_r^1 - \left( \sin 2\phi \sin \delta \cos \lambda + \frac{1}{2} \cos 2\phi \sin 2\delta \sin \lambda \right) u_r^2, \qquad (B11) $$
$$u_{\theta} = -(\sin \phi \cos \delta \cos \lambda + \cos \phi \cos 2\delta \sin \lambda)u_{\theta}^{1} + \left(\frac{1}{2} \sin 2\phi \sin 2\delta \sin \lambda - \cos 2\phi \sin \delta \cos \lambda\right)u_{\theta}^{2}, \quad (B12)$$

where $\delta$ is the dip angle of the fault plane; $\lambda$ is the slip direction measured counterclockwise from the strike of the fault; $\phi$ is the azimuth of the station, measured clockwise from the strike of the fault. The above coefficients are often called the horizontal radiation patterns and were given by different authors (e.g. Aki & Richards 1980; Wang & Herrmann 1980; Helmberger 1983).