Seismological Problem, Seismic Waves and the Seismic Mainshock

1. Introduction

A typical seismogram recorded on the Earth’s surface consists of a faible tremor followed by an abrupt motion with a long tail.[1] The precursory tremor is associated to spherical-shell waves, called primary seismic waves, while the abrupt motion is known as the seismic mainshock. Such a seismogram is sketched in Figure 1. The interpretation of the seismograms was recognized since long as the Seismological Problem (or Lamb’s problem).[2] It is known that this seismic motion originates in a very small focal region, where a short, sudden disturbance occurs. In the absence of the knowledge of the force acting in the seismic focus, the primary waves are derived by the so-called double-couple procedure, based on the solution to the Stokes problem.[3] Apart from inconvenient restrictions to particular orientations of the double couple, the result may include unphysical contributions.[4,5] The mainshock, associated to Rayleigh surface waves,[6] is treated as a vibration problem.[7,8]

The force density acting in the seismic focus has been introduced in Ref. [9,10]. It reads

$$f_i=T{M}_{ij}\delta \left(t\right){\partial }_j\delta (\mathit{R}-{\mathit{R}}_0),$$

(1)

where $M_{ij}$ are the cartesian components of a symmetrical tensor, known as the tensor of the seismic moment; this force acts a very short time T in the seismic focus localized at the position ${\mathit{R}}_0$. It corresponds to a shearing fault. The total force and the total angular momentum of this force density are zero, according to the physical requirements. The primary P (longitudinal) and S (transverse) seismic waves produced by this force in a homogeneous and isotropic elastic body have been derived, in agreement with the recorded seismograms.[10] The force given by equation (1) corresponds to a single rupture in the focus; several successive ruptures may appear (for a so-called structured focus), with corresponding oscillations displayed by the seismic waves.

For limited distances the Earth may be approximated by a half-space with a plane surface. Once arrived at the Earth’s surface, the primary waves generates surface sources which, in turn produce secondary waves, according to Huygens’ principle. These secondary waves have the form of an abrupt wall with a long tail (actually two walls, corresponding to the two primary waves), in agreement with the mainshock exhibited by seismograms.[10] A structured focus may generate oscillations in the mainshock.

Also, the static deformations of a (homogeneous and isotropic) elastic half-space generated by the tensorial force density given by equation (1) have been computed,[9] as well as the vibrations of the half-space. Moreover, the seismic-moment tensor was derived from measurements of the primary waves at Earth’s surface (the so-called Inverse Seismological Problem).[11]

2. Elastic wave equation

The seismic waves are governed by the Navier-Cauchy equation[12]

$${\ddot{u}}_i-c_t^2{\partial }_j{\partial }_j{u}_i-(c_l^2-c_t^2){\partial }_i{\partial }_j{u}_j=T{m}_{ij}\delta \left(t\right){\partial }_j\delta \left(\mathit{R}\right),$$

(2)

where $u_i$ are the Cartesian components of the local diplacement vector, $c_{t,l}$ are the transverse and the longitudinal elastic wave velocities (in a homogeneous and isotropic elastic body), $m_{ij}$ is the symmetrical tensor of the seismic moment divided by the density of the medium and T denotes the short duration of the force localized at the initial moment $t=0$ in the seismic focus placed at $\mathit{R}=0$.[10]

Equation (2) was solved for the seismic waves in Ref. [10] by using the decomposition in Helmholtz potentials and the Kirchhoff formula. A certain regularization procedure was needed, in order to remove unphysical contributions arising from the singular nature of the source term in equation (2). We describe here a different method of solving the above equation, which throws more light upon the singular, unphysical behaviour of the solution.

A direct way to solving equation (2) is to perform a Fourier transform, which gives

$$({\omega }^2-c_t^2{k}^2)u_i-(c_l^2-c_t^2)k_i{k}_j{u}_j=-iT{m}_{ij}{k}_j.$$

(3)

From this equation we get

$$k_i{u}_i=-\frac{iT{m}_{ij}{k}_i{k}_j}{{\omega }^2-c_l^2{k}^2},$$

(4)

which inserted in equation (3) gives

$$u_i=-\frac{iT{m}_{ij}{k}_j}{{\omega }^2-c_t^2{k}^2}-(c_l^2-c_t^2)\frac{iT{m}_{jk}{k}_i{k}_j{k}_k}{({\omega }^2-c_t^2{k}^2)({\omega }^2-c_l^2{k}^2)}.$$

(5)

We denote the fractions on the right in this equation by $u_i^{(1,2)}$ and perform the reverse Fourier tranform. We get

$$\begin{array}{c}{u}_i^{\left(1\right)}=-\frac{iT}{{\left(2\pi \right)}^4}\int d\omega d\mathit{k}\frac{m_{ij}{k}_j}{{\omega }^2-c_t^2{k}^2}{e}^{-i\omega t}{e}^{i\mathit{kR}}=\\ \\ =-\frac{T}{{\left(2\pi \right)}^4}{m}_{ij}{\partial }_j\int d\omega d\mathit{k}\frac{e^{-i\omega t}{e}^{i\mathit{kR}}}{(\omega -c_{t}k+i\epsilon )(\omega +c_{t}k+i\epsilon )}=\\ \\ =\frac{T}{4\pi c_t}{m}_{ij}{\partial }_j\frac{\delta (R-c_{t}t)}{R},\end{array}$$

(6)

where we placed the $\omega$-poles in the lower half-plane ($\epsilon \to 0^{+}$), in order to get waves which obey the causality principle (i.e., they are vanishing for $t<0$ ).

A similar procedure for $u_i^{\left(2\right)}$ gives

$$u_i^{\left(2\right)}=\frac{T}{2c_t{\left(2\pi \right)}^2}{m}_{jk}{\partial }_i{\partial }_j{\partial }_k\frac{1}{R}{\int }_{-\infty }^{+\infty }dk\frac{1}{k^2}{e}^{ik(R-c_{t}t)}-(t\to l).$$

(7)

We can see that the k-integral in equation (7) is improper. We need to give a meaning to this integral.

We may use several procedures of regularizing this integral. For instance, we note that its second derivative is a $\delta$-function, so we may integrate twice the $\delta$-function, with two constants which need to be determined. Another procedure might be to use the integral

$${\int }_{-\infty }^{+\infty }dk\frac{1}{k}{e}^{ik(R-c_{t}t)}$$

(8)

as a principal value; it is $i\pi sgn(R-c_{t}t)$, and we can integrate it with a constant to be determined. Also, we may view the integral in equation (8) as giving a $\theta$-function. All these regularization procedures give different results, and they need a justification.

3. Coulomb-potential regularization

Let us introduce the function

$$F=-\frac{T}{2c{\left(2\pi \right)}^2}\frac{1}{R}{\int }_{-\infty }^{+\infty }dk\frac{1}{k^2}{e}^{ik(R-ct)}.$$

(9)

The solution $u_i^{\left(2\right)}$ reads

$$u_i^{\left(2\right)}=m_{jk}{\partial }_i{\partial }_j{\partial }_k\left(F_l-F_t\right),$$

(10)

where $F_{l,t}$ are obtained by replacing c by $c_{l,t}$. The factor $1/k^2$ in equation (9) may be viewed as the Fourier transform of the Coulomb potential. Indeed,

$$\int d\mathit{R}\frac{1}{R}{e}^{-i\mathit{kR}}{e}^{-\mu R}=\frac{4\pi }{k^2+{\mu }^2},$$

(11)

where $\mu \to 0^{+}$. This may suggest that F should be regularized by

$$F=-\frac{T}{2c{\left(2\pi \right)}^2}\frac{1}{R}{\int }_{-\infty }^{+\infty }dk\frac{1}{k^2+{\mu }^2}{e}^{ik(R-ct)}.$$

(12)

On the other hand, by direct calculations we get from equation (9)

$$\begin{array}{c}\ddot{F}=c^2\Delta F=\frac{cT}{4\pi }\frac{1}{R}\delta (R-ct);\end{array}$$

(13)

in the limit $R\to 0$ the term $c^2\Delta F$ dominates, in comparison with $\ddot{F}$, so F satisfies the equation

$$\ddot{F}-c^2\Delta F=-\frac{cT}{4\pi }\frac{1}{R}\delta \left(ct\right)=-\frac{T}{4\pi }\frac{1}{R}\delta \left(t\right).$$

(14)

If we integrate this equation with respect to the time, we get

$$\Delta \int dtF=\frac{T}{4\pi c^2}\frac{1}{R};$$

(15)

indeed, from equation (9) we have

$$\int dtF=-\frac{T}{4\pi c^2}\frac{1}{R}\int dk\frac{\delta \left(k\right)}{k^2+{\mu }^2}{e}^{ikR}$$

(16)

and

$$\begin{array}{c}\Delta \int dtF=\frac{T}{4\pi c^2}\frac{1}{R}\int dk\frac{k^2\delta \left(k\right)}{k^2+{\mu }^2}{e}^{ikR}=\\ \\ =\frac{T}{4\pi c^2}\frac{1}{R}\int dk\delta \left(k\right)e^{ikR}=\frac{T}{4\pi c^2}\frac{1}{R}.\end{array}$$

(17)

It is easy to see that the $\mu$-regularization does not work for the function F given by equation (12), because we have already imposed the retarded-wave condition, while the $\mu$-regularization requires the presence of both retarded and advanced waves (according to the regularization of the static Coulomb potential, equation (11)). Consequently, we must derive the function F from its wave equation (14), and retain only the retarded solutions. By Fourier transforming equation (14), we get

$$F=\frac{T}{(k^2+{\mu }^2)({\omega }^2-c^2{k}^2)}$$

(18)

and

$$\begin{array}{c}F=\frac{T}{{\left(2\pi \right)}^4}\int d\omega d\mathit{k}\frac{e^{-i\omega t}{e}^{i\mathit{kR}}}{(k^2+{\mu }^2)({\omega }^2-c^2{k}^2)}=\\ \\ =-\frac{T}{2c{\left(2\pi \right)}^{2}R}\int \frac{dk}{k^2+{\mu }^2}\left[e^{ik(R-ct)}-e^{ik(R+ct)}\right],\end{array}$$

(19)

where we placed the $\omega$-poles in the lower half-plane, according to the causality principle (such that $F=0$ for $t<0$). We can see that the term $e^{ik(R+ct)}$ gives in fact a damped contribution $e^{-\mu (R+ct)}$, although, formally, it looks like an advanced wave. The result of the integration in equation (19) is

$$F=-\frac{T}{8\pi cR\mu }\left\{\begin{array}{c}{e}^{\mu (R-ct)}-e^{-\mu (R+ct)},0<R<ct,\\ e^{-\mu (R-ct)}-e^{-\mu (R+ct)},0<ct<R.\end{array}\right\}.$$

(20)

Here we may take the limit $\mu \to 0$ and get

$$F=-\frac{T}{4\pi c}\left[\theta (ct-R)+\frac{ct}{R}\theta (R-ct)\right],$$

(21)

for the retarded wave, where $\theta \left(0\right)=1/2$; the value $1/2$ is expected for a series of continuous functions which approximate the $\theta$-function. This result has been previously obtained in Ref. [10], by solving equation (14) with the Kirchhoff retarded potentials (where the function F has been introduced and equation (14) established, by using the Helmholtz potentials for the Navier-Cauchy equation).

4. Seismic-wave regularization

First we note that the function F given by equation (21) for $ct\ne R$ satisfies the free wave equation. Therefore, the $\theta$-contributions should be removed, and only the contributions for $ct=R$ should be retained. This is valid also for the derivatives of the prefactor of the $\theta$ function. The first-order spatial derivative is

$${\partial }_{i}F=\frac{T}{8\pi c}\left[\frac{x_i}{R^2}(R-ct)\delta (R-ct)+\frac{2ct{x}_i}{R^3}\theta (R-ct)\right],$$

(22)

where we introduced the factor $1/2$ for $\theta \left(0\right)=0$. The first term in the above equation is zero, while its derivatives are not; the second term should be disregarded, except for $R=ct$. The second-order derivative is

$$\begin{array}{c}{\partial }_i{\partial }_{j}F=\frac{T}{8\pi c}\left(\frac{{\delta }_{ij}}{R}-\frac{ct{\delta }_{ij}}{R^2}-\frac{x_i{x}_j}{R^3}+\frac{3ct{x}_i{x}_j}{R^4}\right)\delta (R-ct)+\\ \\ +\frac{T}{8\pi c}\left(\frac{x_i{x}_j}{R^2}-\frac{ct{x}_i{x}_j}{R^3}\right){\delta }^{{}^{\prime }}(R-ct)+\\ \\ +\frac{T}{4\pi c}\left(\frac{ct{\delta }_{ij}}{R^3}-\frac{3ct{x}_i{x}_j}{R^5}\right)\theta (R-ct).\end{array}$$

(23)

From equation (23) we get

$$\Delta F=\frac{T}{4\pi c}\frac{\delta (R-ct)}{R}+\frac{T}{8\pi c}\left(1-\frac{ct}{R}\right){\delta }^{{}^{\prime }}(R-ct),$$

(24)

which differs from equation (13) by the ${\delta }^{{}^{\prime }}$-contribution. Consequently, the ${\delta }^{{}^{\prime }}$-contribution must be removed from ${\partial }_i{\partial }_{j}F$. Also, according to the discussion above, the $\theta$-contribution must be removed, so we are left with the regularized expression

$${\partial }_i{\partial }_{j}F=\frac{T}{8\pi c}\left(\frac{{\delta }_{ij}}{R}-\frac{ct{\delta }_{ij}}{R^2}-\frac{x_i{x}_j}{R^3}+\frac{3ct{x}_i{x}_j}{R^4}\right)\delta (R-ct).$$

(25)

The regularization procedure described above amounts to viewing the function $\delta (R-ct)$ as a function peaked on $R=ct$, of the order $1/l$ over a small distance l, and zero otherwise. Similarly, the function ${\delta }^{{}^{\prime }}(R-ct)$ is of the order $1/l^2$ extended over l. Indeed, this way the ${\delta }^{{}^{\prime }}$-function in equation (23) brings a small contribution, which may be neglected.

Now, it is easy to compute $m_{jk}{\partial }_i{\partial }_j{\partial }_k$ in equation (10). In general, $m_{jk}{\partial }_k$ may be replaced by an external force $f_j$, which is applied to ${\partial }_i{\partial }_{j}F$. It is not permissible to set $ct=R$ in equation (25), because $f_j$ may include derivatives, which, for the prefactor in equation (25), should be computed before setting $ct=R$. For the derivative of the $\delta$-function, we should put $ct=R$ in its prefactor, in accordance to the regularization of the quantity ${\partial }_i{\partial }_{j}F$ discussed above. We get

$$\begin{array}{c}{m}_{jk}{\partial }_i{\partial }_j{\partial }_k=\\ \\ =\frac{T}{8\pi c{R}^3}\left(m_{jj}{x}_i+4m_{ij}{x}_j-\frac{9m_{jk}{x}_i{x}_j{x}_k}{R^2}\right)\delta (R-ct)+\\ \\ +\frac{T}{4\pi c}\frac{m_{jk}{x}_i{x}_j{x}_k}{R^4}{\delta }^{{}^{\prime }}(R-ct),\end{array}$$

(26)

where we set $ct=R$ in the ${\delta }^{{}^{\prime }}$-contribution.

By making use of equations (6) and (10), we get the near-field displacement

$$\begin{array}{c}{u}_i^n=-\frac{T{m}_{ij}{x}_j}{4\pi c_t{R}^3}\delta (R-c_{t}t)+\\ \\ +\frac{T}{8\pi R^3}\left(m_{jj}{x}_i+4m_{ij}{x}_j-\frac{9m_{jk}{x}_i{x}_j{x}_k}{R^2}\right)·\\ \\ ·\left[\frac{1}{c_l}\delta (R-c_{l}t)-\frac{1}{c_t}\delta (R-c_{t}t)\right]\end{array}$$

(27)

and the far-field displacement ($R\gg l$)

$$\begin{array}{c}{u}_i^f=\frac{T{m}_{ij}{x}_j}{4\pi c_t{R}^2}{\delta }^{{}^{\prime }}(R-c_{t}t)+\\ \\ +\frac{T{m}_{jk}{x}_i{x}_j{x}_k}{4\pi R^4}\left[\frac{1}{c_l}{\delta }^{{}^{\prime }}(R-c_{l}t)-\frac{1}{c_t}{\delta }^{{}^{\prime }}(R-c_{t}t)\right].\end{array}$$

(28)

We can see that the far-field displacement consists of two spherical-shell waves, with a scissor-like shape, one longitudinal, propagating with velocity $c_l$, the other transverse, propagating with velocity $c_t$. These are the P and S seismic waves.[10] The relevance of the near-field displacement for the derivation of the seismic-moment tensor has been discussed in Ref. [13].

5. Secondary waves

The focus of a typical earthquake is localized both in space and time, in a point beneath the Earth’s surface. During the short time of releasing the seismic energy in an earthquake, the focus produces two primary waves, which look like spherical shells, propagating with the longitudinal and transverse elastic wave velocities $c_{l,t}$. These primary waves are known in seismology as the P and S seismic waves. Once arrived at Earth’s surface, such a primary wave generates a circular wavefront on the surface, which propagates with a velocity v, greater than the velocity c of the primary wave. The difference between the two velocities goes to zero for large epicentral distances. These wavefronts are localized on the surface in an infinitesimal torus. According to the Huygens principle, they generate secondary waves, which give the seismic mainshock. The displacement is given by the derivatives of some potentials, denoted here generically by $\psi$. These potentials satisfy the wave equation

$$\ddot{\psi }-c^2\Delta \psi =\delta (r-vt)\delta \left(z\right),$$

(29)

where $\mathit{r}$ is the position vector parallel to the surface and t denotes the time; the surface is viewed as a plane surface, placed at $z=0$. The velocity v is considered constant. This equation can be called the mainshock equation. The (homogeneous and isotropic) elastic medium occupies the half-space $z<0$. The above equation is valid for a limited range of epicentral distances r, centered on a value of the order of the depth of the focus.[10]

We introduce the notation

$$F(\mathit{R},t)=\delta (r-vt)\delta \left(z\right),$$

(30)

where $\mathit{R}=(\mathit{r},z)$, and compare equation (29) with the same equation with the source $S=\delta \left(\mathit{R}\right)\delta \left(t\right)=\delta \left(\mathit{r}\right)\delta \left(z\right)\delta \left(t\right)$. The solution of this latter equation is the spherical wave $\delta (R-ct)/4\pi cR$. We note that the source S is singular in a point with four coordinates (time included), while the source F is singular in a set of points, each with three coordinates, placed along a line ($r=vt$). Therefore, the source F is more singular than the source S, so we expect a divergent solution of equation (29). The singularity is more effective for a larger size of the length of the line $r=ct$, i.e. for large r, so we need, at least, a small cutoff wavevector.

We show below that equation (29) has a regular solution for $v\simeq c$, when the source may be treated as a boundary condition, in accordance with the standard procedure.

We perform the Fourier transform of equation (29),

$$({\omega }^2-c^2{K}^2)\psi (\omega ,\mathit{K})=-F(\omega ,\mathit{K}),$$

(31)

where $\mathit{K}=(\mathit{k},\kappa )$ and

$$\begin{array}{c}F(\omega ,\mathit{K})=\int dtd\mathit{R}\delta (r-vt)\delta \left(z\right)e^{i\omega t}{e}^{-i\mathit{KR}}=\\ \\ =\frac{1}{v}\int d\mathit{r}{e}^{i\omega r/v}{e}^{-i\mathit{kr}}=\frac{2\pi }{v}\int dr·r{e}^{i\omega r/v}{J}_0\left(kr\right),\end{array}$$

(32)

where $J_0\left(kr\right)$ is the Bessel function. This expression can also be written as

$$\begin{array}{c}F(\omega ,\mathit{K})=\frac{2\pi }{v}\frac{\partial }{\partial (i\omega /v)}\int dr{e}^{i\omega r/v}{J}_0\left(kr\right)=\\ \\ =-\frac{2\pi i}{v{k}^2}\frac{\partial }{\partial \lambda }{\int }_0^{\infty }dx{e}^{i\lambda x}{J}_0\left(x\right),\end{array}$$

(33)

where $\lambda =\omega /vk$. The integral in the second row of equation (33) is the Weyl-Sommerfeld integral,

$$I\left(\lambda \right)={\int }_0^{\infty }dx{e}^{i\lambda x}{J}_0\left(x\right)=\frac{\theta (1-\mid \lambda \mid )}{\sqrt{1-{\lambda }^2}}+isgn\left(\lambda \right)\frac{\theta (\mid \lambda \mid -1)}{\sqrt{{\lambda }^2-1}}.$$

(34)

Therefore, from equation (31) we get

$$\begin{array}{c}\psi (\mathit{R},t)=\frac{i}{{\left(2\pi \right)}^{3}v}\int d\omega d\mathit{K}\frac{\partial I/\partial \lambda }{k^2(\omega -cK+i\epsilon )(\omega +cK+i\epsilon )}{e}^{-i\omega t}{e}^{i\mathit{KR}},\end{array}$$

(35)

where $\epsilon \to 0^{+}$. We place the $\omega$-poles in the lower half-plane in order to have $\psi =0$ for $t<0$, according to the causality principle. Henceforth, we consider $t>0$ only. Equation (35) can also be written as

$$\begin{array}{c}\psi (\mathit{R},t)=\frac{i}{{\left(2\pi \right)}^3{v}^2}\int d\lambda d\mathit{K}\frac{\partial I/\partial \lambda }{k^3(\lambda -{\lambda }_1)(\lambda -{\lambda }_2)}{e}^{-ivkt\lambda }{e}^{i\mathit{KR}},\end{array}$$

(36)

where ${\lambda }_1=cK/vk-i\epsilon$ and ${\lambda }_2=-cK/vk-i\epsilon$. We need to compute $\partial I/\partial \lambda$ for ${\lambda }_{1,2}$, by using equation (34). We get

$${(\partial I/\partial \lambda )}_{{\lambda }_1}=\left\{\begin{array}{c}\frac{{\lambda }_0}{{(1-{\lambda }_0^2)}^{3/2}},{\lambda }_0=cK/vk<1,\\ -\frac{i{\lambda }_0}{{({\lambda }_0^2-1)}^{3/2}},{\lambda }_0=cK/vk>1\end{array}\right.$$

(37)

and

$${(\partial I/\partial \lambda )}_{{\lambda }_2}=\left\{\begin{array}{c}-\frac{{\lambda }_0}{{(1-{\lambda }_0^2)}^{3/2}},{\lambda }_0=cK/vk<1,\\ -\frac{i{\lambda }_0}{{({\lambda }_0^2-1)}^{3/2}},{\lambda }_0=cK/vk>1.\end{array}\right.$$

(38)

It follows

$$\begin{array}{c}\psi (\mathit{R},t)=\frac{1}{{\left(2\pi \right)}^2{v}^2}\int d\mathit{K}\frac{1}{k^3}\left\{\begin{array}{c}\frac{coscKt}{{(1-{\lambda }_0^2)}^{3/2}}\\ -\frac{sincKt}{{({\lambda }_0^2-1)}^{3/2}}\end{array}\right\}{e}^{i\mathit{KR}}=\\ \\ =\frac{v}{{\left(2\pi \right)}^2}\int d\mathit{K}\left\{\begin{array}{c}\frac{coscKt}{{(v^2{k}^2-c^2{K}^2)}^{3/2}}\\ -\frac{sincKt}{{(c^2{K}^2-v^2{k}^2)}^{3/2}}\end{array}\right\}{e}^{i\mathit{KR}},\end{array}$$

(39)

where the expressions under the square root are positive.

The function $\psi (\mathit{R},t)$ given by equation (39) is a superposition of plane waves $e^{icKt+\mathit{iKR}}$ and $e^{-icKt+\mathit{iKR}}$; we must retain only the outgoing wave $e^{-icKt+\mathit{iKR}}$, in accordance with the causality principle, so equation (39) becomes

$$\psi (\mathit{R},t)=\frac{v}{2{\left(2\pi \right)}^2}\int d\mathit{k}\left\{\begin{array}{c}\frac{1}{{(v^2{k}^2-c^2{K}^2)}^{3/2}}\\ -\frac{i}{{(c^2{K}^2-v^2{k}^2)}^{3/2}}\end{array}\right\}{e}^{-icKt}{e}^{i\mathit{kR}}$$

(40)

(where in the second row we must add the complex conjugate and divide by 2). Also, the wave should be progressive, i.e.$ct>R$, a condition which can also be writen as

$$c^2{t}^2>R^2=r^2+z^2;$$

(41)

otherwise, the wave is zero.

6. Mainshock

For v close to c only the second row in equation (40) is valid. By using ${\gamma }^2=v^2/c^2-1$, equation (40) can be written as

$$\begin{array}{c}\psi (\mathit{R},t)=-\frac{iv}{2{\left(2\pi \right)}^2{c}^3}\int d\mathit{k}d\kappa \frac{e^{-ic\sqrt{{\kappa }^2+k^2}t}}{{({\kappa }^2-{\gamma }^2{k}^2)}^{3/2}}{e}^{i\kappa z}{e}^{i\mathit{kr}}=\\ \\ =-\frac{iv}{4\pi c^3}\int dkk{J}_0\left(kr\right)\int d\kappa \frac{1}{{({\kappa }^2-{\gamma }^2{k}^2)}^{3/2}}{e}^{-ic\sqrt{{\kappa }^2+k^2}t}{e}^{i\kappa z}.\end{array}$$

(42)

For $z<0$ the $\kappa$-integration must be carried out in the lower half-plane. The integrand has two branch points at $\kappa =\pm \gamma k$. It is easy to see that the integral along this branch cut is singular, as expected. According to the discussion above, we put $v=c$ and displace the pole $\kappa =0$ slightly below in the lower half-plane. This operation provides the standard procedure of treating the source as a boundary condition. For $v=c$ ($\gamma =0$) the integral in equation (42) becomes

$$\psi (\mathit{R},t)\simeq -\frac{i}{4\pi c^2}\int dkk{J}_0\left(kr\right)\int d\kappa \frac{1}{\kappa {}^3}{e}^{-ic\sqrt{{\kappa }^2+k^2}t}{e}^{i\kappa z}.$$

(43)

The pole placed slightly below $\kappa =0$ plays the role of a lower cutoff wavevector. The calculation of this contribution is performed by writing

$$\begin{array}{c}{e}^{-ic\sqrt{{\kappa }^2+k^2}t}{e}^{i\kappa z}=e^{-ictk}{e}^{-ict{\kappa }^2/2k}{e}^{i\kappa z}\simeq \\ \\ \simeq e^{-ictk}\left(1+i\kappa z-\frac{{\kappa }^2{z}^2}{2}-i\frac{ct{\kappa }^2}{2k}+...\right),\end{array}$$

(44)

which leads to

$$\begin{array}{c}\psi (\mathit{R},t)\simeq \frac{i}{8\pi c^2}\{z^2\int dkk{J}_0\left(kr\right)e^{-ictk}+\\ \\ +ict\int dk{J}_0\left(kr\right)e^{-ictk}\}\int d\kappa \frac{1}{\kappa }.\end{array}$$

(45)

Straightforward calculations give

$$\begin{array}{c}\psi (\mathit{R},t)\simeq \frac{i}{4c^2}\{\frac{z^2}{r^2}\frac{\partial }{\partial \lambda }\int dx{e}^{-i\lambda x}{J}_0\left(x\right)+\\ \\ +\lambda \int dx{e}^{-i\lambda x}{J}_0\left(x\right)\},\end{array}$$

(46)

where $\lambda =ct/r$ (>1 ). By making use of equation (34), we get

$$\psi (\mathit{R},t)\simeq \frac{1}{4c^2}\frac{\left(c^2{t}^2-r^2-z^2\right)ct}{{(c^2{t}^2-r^2)}^{3/2}}$$

(47)

for $ct>r$ (i.e.$c^2{t}^2>r^2$, and $c^2{t}^2>r^2+z^2$). This is precisely the result obtained previously.[10] In order to account for the small difference between the two velocities, from the denominator in equation (42) we may infer that K should be replaced by $Kv/c$ in the exponent of equation (43), which amounts to replacing the time t by the retarded time $\tau =tc/v$.

We can check by direct calculations that

$$\ddot{\psi }-c^2\Delta \psi =\ddot{\psi }-c^2\left(\frac{{\partial }^2\psi }{\partial r^2}+\frac{\partial \psi }{r\partial r}+\frac{{\partial }^2\psi }{\partial z^2}\right)=0,$$

(48)

except for $z\to 0$, $r\to ct$. The singularity at $ct=r$ in equation (47) (for $z=0$), arises from the sharpness of the $\delta$-functions of the source term in equation (29). It can be smoothed out by replacing $ct-r$ by $l,$ where l is an infinitesimal distance.

The solution given by equation (47) has a spherical wavefront $r^2+z^2=c^2{t}^2$; it has a rapid variation with r and a rather slow variation with z, so it may be viewed as a quasi-cylindrical wave, with a wall-like structure. It corresponds to the seismic mainshock. We can see that the potential given by equation (47) and its spatial derivatives (the displacement) look like an abrupt wall with a long tail, propagating with velocity c in the interior of the Earth and on its surface. Actually, we get two such walls, corresponding to the two primary P and S waves, propagating with velocities $c_{l,t}$.

7. Site response and inhomogeneities

The $\delta$-functions occurring in these problems should be viewed as highly-peaked functions over a small region. For instance, $\delta \left(\mathit{R}\right)$ which occurs in the tensorial force acting in the earthquake focus (equation (1)) is approximately $1/l^3$ over a small region with dimension l, where l is of the order of the dimension of the seismic focus. The cutoff length l occurs also in the primary seismic waves derived above (equations (28)). It is related to the elastic energy stored in the seismic focus, and released during an earthquake. A measure of the seismic energy is the earthquake (moment) magnitude $M_w$, such that, for instance, for an earthquake with magnitude $M_w=7$ we get $l=316m$ (for a density $\rho =5g/c{m}^3$ of the Earth and an average velocity $c=5km/s$ of the elastic waves). However, the extension of the spot left by the seismic waves on the Earth’s surface is much larger. This is so, because of the energy loss suffered by the seismic waves during their propagation through the Earth.

The results presented above relate to a homogeneous and isotropic medium, while the Earth is recognized as inhomogeneous and anisotropic. The spatial distribution of a wave is characterized by its Fourier transform. Let us take a far-field seismic wave of the form

$$u=\frac{{\delta }^{{}^{\prime }}(R-ct)}{R};$$

(49)

its Fourier transform is

$$u\left(\mathit{K}\right)=\int d\mathit{R}\frac{{\delta }^{{}^{\prime }}(R-ct)}{R}{e}^{-i\mathit{kR}}=-4\pi coscKt.$$

(50)

When encountering an inhomogeneity, this wave sets in motion its particles, and even the inhomogeneity as a whole. Consequently, the incident wave loses energy, which is transferred to the inhomogeneity, which, in turn, generates secondary waves; part of the energy may be dissipated. Obviously, the effect is larger for small wavelengths, which are numerous, due to the large number of distinct directions for a large K. This is the wave scattering, with a possible energy loss (absorption). It is reasonable to assume that inhomogeneities are distributed relatively uniformly, over their various size and mean separation distance. Consequently, we expect a secondary (scattered) diffuse radiation, with a large content of small wavelengths. It follows that the incident wave content is diminished isotropically, with a larger weigth for small wavelengths. This amounts to modify the Fourier transform given above according to

$$u\left(\mathit{K}\right)=-4\pi e^{-\alpha K}coscKt,$$

(51)

where the parameter $\alpha$ characterizes the inhomogeneities distribution. The resulting reverse Fourier transform is

$$\begin{array}{c}u=-\frac{4\pi }{{\left(2\pi \right)}^3}\int d\mathit{K}{e}^{-\alpha K}coscKt{e}^{i\mathit{KR}}=\\ \\ =\frac{1}{\pi R}\frac{\partial }{\partial R}\frac{\alpha }{{(R-ct)}^2+{\alpha }^2},\end{array}$$

(52)

where we retain only the retarded waves. We can see that we recover the incident wave given by equation (49) in the limit $\alpha \to 0$. The effect of the inhomogeneities, included in the parameter $\alpha$, is a flattening of the ${\delta }^{{}^{\prime }}$-incident wave, which gets a a larger width $l_0=2\alpha /\sqrt{3}>l$ and a smaller heigth $\simeq 1/l_0^2$. The ratio $l_0/l$ may attain values of the order of ten.[14] The scissor-like structure of the primary seismic waves shows that the frequency content of these waves is mainly centered on a single frequency, of the order of the wave velocity to the dimension of the seismic focus ($c/l_0)$. Therefore, the Fourier analysis of the primary waves may give an estimate of the dimension of the focus.[14]

The P and S seismic waves and the mainshock derived above are the seismic motion, generated by a point-like seismic focus acting a short time interval. These results are obtained by assuming a homogeneous and isotropic elastic medium (the Earth). This assumption is valid for an average of the elastic properties of the medium. As long as we are interested in the overall, average behaviour of the elastic motion, this is a satisfactory hypothesis. However, if we are interested in the local elastic motion, the elastic particularities of the site should be taken into account. In a simple model, any site may be viewed as a damped harmonic oscillator, with frequency ${\omega }_g$, connected by elastic forces to its surroundings. The seismic motion acts as an external force upon such an oscillator. The resulting motion is the ground motion. It consists of the original, scissor-like primary seismic waves and the wall-like seismic mainshock, over which the damped ${\omega }_g$-oscillations of the site are superposed; these damped oscillations are the seismic response of the site. In the ground motion the long tail of the seismic mainshock is governed by the damping coefficient of the site. In turn, the ground motion acts as an external force upon the seismographs (or the structures built on the Earth’s surface). As a simple model, we may take a linear damped harmonic oscillator for the seismograph, with frequency ${\omega }_s$, such that the seismograms record the ${\omega }_{g,s}$-oscillations, superposed over the original seismic motion. All these results are included in Ref. [14], where an estimation is also given of the maximum (peak) values of the ground-motion displacement, velocity and acceleration, which may be useful as input parameters for seismic hazard studies.

8. Concluding remarks

The primary seismic waves and the sesimic mainshock are derived for a homogeneous and isotropic half-space with a plane surface (the Earth), by solving the elastic wave equation with the seismic tensorial force acting a very short time (time-impulse) in a localized (point-like) seismic focus. The solution is the seismic motion. It is known as the Seismological, or Lamb’s, Problem. The results reported previously are re-derived by using a new method, which emphasizes the regularization procedure. The time-impulse and point-like tensorial force is a combination of a temporal $\delta$-function and derivatives of a spatial $\delta$-functions (equation (1)). The equation of the elastic motion with such a source term may be called a singular equation. The singular nature of such a source leads to singular, improper solutions and unphysical contributions. Therefore, a regularization procedure is necessary. The regularization procedure employed here includes both a Coulomb-potential type regularization, due to the use of the Fourier expansions, and a seismic-wave regularization, specific to the singular elastic wave equation. This regularization procedure amounts to using a long wavelength cutoff (Coulomb potential), an approximation by continuous-function series for the step $\theta$-function, the removal of superfluous solutions of the free-wave equation and a small cutoff time/length for the peaked temporal/spatial $\delta$-distribution.

In addition, a summary information is provided for the effect of the inhomogeneities and the site response in determining the ground motion and the seismographs’ recordings.