Dynamic antiplane rupture propagation crossing a material interface: modelling with BIEM

Abstract

Because the Earth's crust is quite heterogeneous and consists of various types of rocks, many faults are located near material interfaces. Recent observational studies have revealed that some faults cross the interfaces, but theoretical treatment of their dynamic behaviour is still undeveloped. In this paper, we develop a boundary integral equation method (BIEM) to analyse the dynamic behaviour of antiplane rupture propagation crossing a material interface. First, we obtain an exact solution of the stress response due to a rectangular slip rate function on one discretized crack element embedded in one of the elastic half-spaces welded along a planar interface. This solution is helpful not only in BIEM simulations of rupture propagation but also in benchmark testing or extension of other numerical schemes. Next, we simulate the dynamic propagation of an antiplane rupture crossing the material interface by using the exact solution for a stress response in a BIEM framework. We find that a shear wave reflected from the material interface generates a significant change in the slip rate on the crack depending on the angle between the interface and the crack, contrast in the elastic properties, and rupture velocity. Moreover, we infer from our numerical results and previous related works that a backward propagating healing front emerges under rate-weakening friction owing to interaction between the reflected wave and friction.

Numerical solutions, Numerical approximations and analysis, Earthquake dynamics, Dynamics and mechanics of faulting

1 INTRODUCTION

Many faults are known to lie in the strongly heterogeneous crust, which is often modelled as a sectionally homogeneous elastic medium. Recent observational studies using dense seismic networks have revealed that some faults cross the interfaces between dissimilar media. Hirano & Yamashita (2011) enumerated some examples of earthquakes along such faults (e.g. the 2004 mid-Niigata Prefecture earthquake (Kato et al.2006), the 2008 Iwate-Miyagi Nairiku earthquake (Okada et al.2010), the 2001 Geiyo earthquake (Wang et al.2004) and an M2 earthquake of 2007 December 27 in a South African deep gold mine (Yabe et al.2009; Naoi et al.2011). Moreover, Ferrari & Bonafede (2012) suggested that the strike-slip fault that caused the M6.6 earthquake of 2000 June 17 in South Iceland cuts not only a free surface but also an interface between a shallow layer and bedrock. In addition, geologically, many mesoscopic outcrops where faults cut strata have been reported (Roche et al.2012).

Although many faults cross material interfaces, as mentioned above, there are very few theoretical studies of the dynamic behaviour of such faults because its mathematical treatment is difficult. Faults are generally considered as cracks embedded in elastic media, whose response can be simulated using numerical schemes. By using an extended finite element method (XFEM, Wada & Goto 2012) analysed spontaneous rupture of a mode-II crack crossing an interface between a surface layer and bedrock. Although their results could explain some patterns of the ground motion spectra, in terms of fault dynamics, their numerical results for the spatio-temporal distribution of the slip velocity on the crack were of very low resolution and had conspicuous numerical oscillations. Hence, the XFEM they used currently has some numerical difficulties, and a more precise method might be needed to treat dynamic rupture in such a layered medium.

On the other hand, another numerical technique has been developed to treat the dynamic rupture behaviour in sectionally homogeneous elastic media. Kame & Kusakabe (2012) and Kusakabe (2013) proposed an extended boundary integral equation method (XBIEM) as an extension of a conventional BIEM (e.g. Cochard & Madariaga 1994; Tada & Yamashita 1997; Tada & Madariaga 2001). In the conventional BIEM, one solves an integral equation corresponding to a representation theorem (Aki & Richards 2002) by discretizing an integration path representing the path of a crack embedded in a homogeneous medium. The BIEM outputs a smoother numerical solution than an FEM because the boundary conditions on the crack surface are strictly satisfied in the BIEM framework. According to Tamura & Ide (2011), an FEM requires tenfold discretized elements along the crack path compared to the BIEM if one attempts to obtain smooth numerical solutions by these two methods at the same level. In the XBIEM framework, Kame & Kusakabe (2012) and Kusakabe (2013) showed a way to treat cracks embedded in sectionally homogeneous elastic media considering a term in the representation theorem that is associated with the traction acting on the material interfaces, so that the integration path represents not only the crack path but also every material interface. Therefore, the XBIEM has both advantages and disadvantages; although multilayered media and non-planar material interfaces can be treated in its framework, such a system requires too many discretized elements and a long computational time because every interface must be discretized.

The XFEM and XBIEM have other problems in addition to the long calculation time. The numerical accuracy of newly suggested computational methods should generally be validated against some accurate solution, and an exact solution is the best benchmark. For a crack embedded in a homogeneous medium, Tada & Madariaga (2001) obtained an exact solution for the stress response to the rectangular slip rate function given on one discretized crack element. In the same way, if an exact solution for the crack element embedded in an inhomogeneous medium is obtained, it can be the best benchmark for any numerical schemes, including the XFEM and XBIEM. Moreover, in the framework of a conventional BIEM, the exact solution is equivalent to the convolution coefficients of the discretized integral equations that represent the stress response of a homogeneous medium due to the slip on a crack (Tada & Madariaga 2001). This implies that even dynamic rupture simulations can be executed in the same framework as the conventional BIEM if the exact solution for an inhomogeneous medium is obtained; that is the numerical integrals along the material interface appearing in the framework of the XBIEM are no longer required.

In this paper, we formulate a general representation of the convolution coefficients when the Green's function with respect to the antiplane displacement field is given for some inhomogeneous media. In particular, we derive an explicit form of the convolution coefficients for a bimaterial system consisting of two welded half-spaces. To do this, we use the Cagniard–de Hoop method, which enables us to obtain the time-domain solution of the Green's function for the 2-D wave equation in a bimaterial system (Aki & Richards 2002). A key characteristic of our mathematical procedures is that the Green's function is treated as a function of a complex ray parameter appearing in the Cagniard–de Hoop method. In addition, we simulate the dynamic behaviour of a rupture crossing a bimaterial interface using the convolution coefficients. We assume a frictionless crack in this paper because various friction laws under a high slip rate have been suggested (e.g. Rice 2006; Ampuero & Ben-Zion 2008; Noda et al.2009), and the most appropriate one is still unclear. Hence, we quantitatively simulate the behaviour of the frictionless crack and qualitatively discuss some effects of the rate-dependent friction laws on the dynamic fault behaviour.

2 CONVOLUTION COEFFICIENTS FOR AN INHOMOGENEOUS MEDIUM

2.1 Representation and physical meaning of the convolution coefficients for a general inhomogeneous medium

Here we derive a general representation of the convolution coefficients for an antiplane crack embedded in an inhomogeneous medium whose Green's function of displacement is given. We consider a 2-D antiplane displacement field u(x, z, t) in a x–y coordinate system at time t. In the antiplane case, the expression for u given by the representation theorem of Aki & Richards (2002, hereinafter AR02) can be reduced to

\begin{eqnarray} u(x,z,t) = \int _\Gamma {\rm d}\gamma \int _{-\infty }^{+\infty } {\rm d}\tau \, \left[ u(\gamma , \tau ) \right] \mu (\gamma ) \frac{\partial }{\partial \zeta } G(x,z,t; \gamma , \tau ), \nonumber\\ \end{eqnarray}

(1)

where γ represents a position on the crack path Γ, and the brackets [·] denote a jump of u between both sides of the crack. We assume here that the rigidity μ depends on the position because we consider an inhomogeneous medium. For the integrand of eq. (1), we set a local Cartesian ξ–ζ coordinate system whose origin is at γ, where the ξ and ζ axes are locally tangent and normal to Γ, respectively [ξ does not appear in eq. (1) explicitly; see Fig. 1a]. Eq. (1) holds even if the crack lies along a sharp material interface because the product of μ(γ) and ∂_ζG, which gives the shear stress acting on the interface, is continuous across the interface, and the displacement field u is determined uniquely (Ampuero & Dahlen 2005).

Figure 1.

(a) ξ-ζ coordinate system defined at the point γ on the crack Γ. (b) schematic illustration of a discretized crack.

Below, the subscript q denotes x or z. By differentiating both sides of eq. (1) with respect to the spatial coordinates, we obtain the qy component of the stress tensor as

\begin{eqnarray} \sigma _{qy}(x,z,t) &=& \mu (x,z) \int _\Gamma {\rm d}\gamma \int _{-\infty }^{+\infty } {\rm d}\tau \,\nonumber\\ &&\times\, \left[ u(\gamma , \tau ) \right] \mu (\gamma ) \frac{\partial ^2}{\partial q \partial \zeta } G(x, z, t; \gamma , \tau ), \end{eqnarray}

(2)

where μ(x, z) denotes the rigidity at (x, z). The integral on the RHS of eq. (2) cannot be calculated numerically because the derivative of the Green's function G is hypersingular (Tada & Yamashita 1997). Many authors have solved this problem by expanding the slip rate with piecewise constant basis functions:

\begin{eqnarray} \frac{\partial }{\partial \tau } \left[u(\gamma , \tau ) \right] &\sim &\sum _{i, n} V^{i,n} \left\lbrace H(\gamma - \gamma _i) - H(\gamma - \gamma _{i-1}) \right\rbrace \nonumber\\ &&\times\, \left\lbrace H(\tau - \tau _n) - H(\tau - \tau _{n-1}) \right\rbrace \end{eqnarray}

(3)

\begin{eqnarray} \Leftrightarrow \frac{\partial }{\partial \gamma }\frac{\partial ^2}{\partial \tau ^2} \left[u(\gamma , \tau ) \right] &\sim &\sum _{i, n} V^{i,n} \left\lbrace \delta (\gamma - \gamma _i) - \delta (\gamma - \gamma _{i-1}) \right\rbrace \nonumber\\ &&\times\,\left\lbrace \delta (\tau - \tau _n) - \delta (\tau - \tau _{n-1}) \right\rbrace \end{eqnarray}

(4)

(e.g. Cochard & Madariaga 1994; Tada & Yamashita 1997; Tada & Madariaga 2001), where H( · ) is the Heaviside function, and δ( · ) is Dirac's delta function. In this approximation, a curvilinear crack path Γ is replaced by a set of straight line segments (see Fig. 1b). The indices i and n denote the ith line segment and nth time step, respectively, and |$V^{i,n} := \left[\dot{u}(\gamma _i,\tau _n) \right]$| without the Einstein summation convention is the piecewise constant slip rate on a spatio-temporal domain (γ_(i−1), γ_i) × (τ_(n−1), τ_n). If we integrate the RHS of eq. (2) by parts once for γ and twice for τ, the approximation (4) yields:

\begin{eqnarray} \sigma _{qy}(x,z,t) \sim \sum _{i,n} V^{i,n} K_{qy}^{i,n}(x,z,t), \end{eqnarray}

(5)

where the convolution coefficients |$K_{qy}^{i,n}$| are represented as follows:

\begin{eqnarray} K_{qy}^{i,n}(x,z,t) & =& - \mu (x,z) \mu \left(\gamma \right) \, \left[ \left[ I_{qy} \right]_{\gamma _{i-1}}^{\gamma _i} \right]_{\tau _{n-1}}^{\tau _n}, \end{eqnarray}

(6)

\begin{eqnarray} I_{qy}(x,z,t) & =& \int {\rm d}\gamma \int {\rm d}\tau \int {\rm d}\tau \frac{\partial ^2}{\partial q \partial \zeta } G(x, z, t; \gamma , \tau ). \end{eqnarray}

(7)

Aside from the explicit representations of the convolution coefficients, eq. (5) has the same form as the equation obtained by discretizing an integral representation of the stress response to the slip rate in the conventional BIEM for a homogeneous medium (e.g. eq. 7 of Tada & Madariaga 2001). Eqs (6) and (7) are the general representations of the convolution coefficients. Thanks to the above integrations by parts, the convolution coefficients become bounded and can be calculated numerically because the integrations reduce the degree of singularities of the Green's function and its derivatives (Cochard & Madariaga 1994; Tada & Madariaga 2001). Hence, dynamic rupture propagation can be numerically simulated following eq. (68) of Tada (2009) for an implicit time-marching scheme or eq. (70) of Tada (2009) for an explicit time-marching scheme.

The physical meaning of the convolution coefficients is the stress response to the rectangular slip rate function given on one discretized crack element, because we can obtain

\begin{eqnarray} \sigma _{qy}(x,z,t) = K_{qy}^{i,n}(x,z,t) \end{eqnarray}

(8)

for every specific set of (i, n) after substituting

\begin{eqnarray} \frac{\partial }{\partial \gamma }\frac{\partial ^2}{\partial \tau ^2} \left[ u(\gamma , \tau ) \right] &=& \left\lbrace \delta (\gamma - \gamma _i) - \delta (\gamma - \gamma _{i-1}) \right\rbrace \nonumber\\ &&\times\,\left\lbrace \delta (\tau - \tau _n) - \delta (\tau - \tau _{n-1}) \right\rbrace \end{eqnarray}

(9)

into eq. (2). This stress response can also be obtained numerically by other schemes that can be adopted for dynamic rupture simulation in an inhomogeneous medium, for example the XBIEM and XFEM. Therefore, the strict solution (8) can be the best benchmark for evaluating the precision of those numerical schemes.

2.2 Convolution coefficients for a bimaterial system

In this subsection, Chapter 6.4 of AR02 will be simply referred to as ‘AR02’. In a bimaterial system, the Green's function is obtained as the imaginary part of a complex function using the Cagniard–de Hoop method [AR02]. Moreover, the method suggests that the Green's function has terms representing the direct, reflected, and transmitted waves. Hence, we denote those terms here as ℑG_D, ℑG_R and ℑG_T, respectively, where ℑ indicates the imaginary part. Below, the material constants, including the rigidity μ, shear wave velocity β and density ρ, have subscripts j or k, where j denotes a half-plane with an impulsive wave source whose support is at (x^′, z^′, τ), and k denotes the other half-plane without the wave source. We give the components of the convolution coefficients due to ℑG_D, ℑG_R and ℑG_T as follows:

\begin{eqnarray} K_{qy}^{i,n} &=& \left\lbrace \begin{array}{rll}K_{qy,D}^{i,n} + K_{qy,R}^{i,n} &\quad (\textrm {if} \, \gamma _i \, \textrm {and} \, (x,z) \, \textrm {are located}\\ &\quad \textrm{in the same half-space}), \\ K_{qy,T}^{i,n} &\quad (\textrm {if} \, \gamma _i \, \textrm {and} \, (x,z) \, \textrm {are located}\\ &\quad \textrm{in the different half-spaces}). \end{array} \right. \end{eqnarray}

(10)

\begin{eqnarray} = \left\lbrace \begin{array}{l}- \mu _j^2 \Im \left[ \left[ I_{qy,D} + I_{qy,R} \right]_{\gamma _{i-1}}^{\gamma _i} \right]_{\tau _{n-1}}^{\tau _n}, \\ - \mu _j \mu _k \Im \left[ \left[ I_{qy,T} \right]_{\gamma _{i-1}}^{\gamma _i} \right]_{\tau _{n-1}}^{\tau _n}. \end{array} \right. \end{eqnarray}

(11)

ℑG_D is the same as the Green's function for a homogeneous medium, and an explicit form of ℑI_qy,D has already been obtained (eqs A2 and A3 of Tada & Madariaga 2001). This could be done only by real analysis because ℑG_D can be represented in the simple form

\begin{eqnarray} \Im G_D = \frac{1}{2 \pi \rho \beta ^2} \frac{H\left( t - \tau - R/\beta \right)}{\sqrt{\left(t - \tau \right)^2 - R^2/\beta ^2}}, \end{eqnarray}

(12)

where |$R=\sqrt{\left(x-x^{\prime }\right)^2+\left(z-z^{\prime }\right)^2}$|⁠. On the other hand, I_qy,R and I_qy,T have not been obtained. As shown in AR02, ℑG_R and ℑG_T are functions of a complex ray parameter p(x, z, t; γ, τ). By considering complex analysis, we successfully executed the multiple integration in eq. (7). In this subsection, we show the entire procedure for calculating I_qy,R and I_qy,T. The details of their derivation are shown in Appendix.

We assume here that the bimaterial interface lies along the plane z = 0. By using the Cagniard–de Hoop method with the one-sided Laplace transform, G_R and G_T are represented in the Laplace domain as

\begin{eqnarray} G_R & =& \int _\tau ^\infty \frac{1}{2 \pi \mu _j s \eta _j} \left\lbrace \frac{\mu _j \eta _j - \mu _k \eta _k}{\mu _j \eta _j + \mu _k \eta _k} \right\rbrace \frac{{\rm d}p}{{\rm d}t} {\rm e}^{-s (t-\tau )} {\rm d}t, \end{eqnarray}

(13)

\begin{eqnarray} t - \tau & = p (x - x^{\prime }) + \eta _j \left| z + z^{\prime } \right| \end{eqnarray}

(14)

and

\begin{eqnarray} G_T & =& \int _\tau ^\infty \frac{1}{2 \pi \mu _j s \eta _j} \left\lbrace \frac{2 \mu _j \eta _j}{\mu _j \eta _j + \mu _k \eta _k} \right\rbrace \frac{{\rm d}p}{{\rm d}t} {\rm e}^{-s (t-\tau )} {\rm d}t, \end{eqnarray}

(15)

\begin{eqnarray} t - \tau & = p (x - x^{\prime }) + \eta _j \left| z^{\prime } \right| + \eta _k \left| z \right|, \end{eqnarray}

(16)

respectively [AR02], where s denotes the Laplace variable; η_j and η_k are defined as

\begin{eqnarray} \eta _j & = \sqrt{\beta _j^{-2} - p^2}, \quad \eta _k = \sqrt{\beta _k^{-2} - p^2}, \end{eqnarray}

(17)

and the real parts of both η_j and η_k are taken to be positive. The complex ray parameters p for G_R and G_T are given as the roots of eqs (14) and (16), respectively. By substituting eqs (13) and (15) into eq. (7), we obtain the following solutions:

\begin{eqnarray} \textrm {sgn}\left(z^{\prime }\right) \pi \mu _j I_{xy,R} & =& - \mu _j \left\lbrace C_j \Phi _{jj}(M) + C_k \Phi _{jk}(M) + B \sin {\theta ^{\prime }} \right\rbrace\nonumber\\ &&+\, \frac{1}{2} \left\lbrace p + A \sin {\theta ^{\prime }} \right\rbrace , \end{eqnarray}

(18)

\begin{eqnarray} \!\pi \mu _j I_{zy,R} & =& + \mu _j \left\lbrace C_j \Psi _j(M) + C_k \Psi _k(M) - B \cos {\theta ^{\prime }} \right\rbrace \nonumber\\&& +\, \frac{1}{2} \left\lbrace \eta _j + A \cos {\theta ^{\prime }} \right\rbrace , \end{eqnarray}

(19)

\begin{eqnarray} \textrm {sgn}\left(z^{\prime }\right) \pi \mu _k I_{xy,T} & = + \mu _k \left\lbrace C_j \Phi _{jj}(M) + C_k \Phi _{jk}(M) + B \sin {\theta ^{\prime }} \right\rbrace , \nonumber\\ \end{eqnarray}

(20)

\begin{eqnarray} \!\pi \mu _k I_{zy,T} & =& + \mu _j \left\lbrace C_j \Psi _j(M) + C_k \Psi _k(M) - B \cos {\theta ^{\prime }} \right\rbrace \nonumber\\&& + \, \left\lbrace \eta _j + A \cos {\theta ^{\prime }} \right\rbrace , \end{eqnarray}

(21)

where

\begin{eqnarray} C_j & =& \frac{\mu _j}{\mu _j^2 - \mu _k^2}, \quad C_k = \frac{\mu _k}{\mu _k^2 - \mu _j^2}, \quad M^2 = \frac{\mu _j \rho _j - \mu _k \rho _k}{\mu _j^2 - \mu _k^2}, \end{eqnarray}

(22)

\begin{eqnarray} \theta ^{\prime } & = \textrm {sgn}\,(z^{\prime }) \theta , \end{eqnarray}

(23)

\begin{eqnarray} A & =& \frac{1}{\beta _j} \log \left\lbrace \frac{\beta _j p - \sin {\theta ^{\prime }}}{\beta _j \eta _j + \cos {\theta ^{\prime }}}\right\rbrace , \end{eqnarray}

(24)

\begin{eqnarray} B & =& \frac{1}{\beta _j} \left\lbrace C_j \Omega _{jj} + C_k \Omega _{jk} \right\rbrace , \end{eqnarray}

(25)

and, for {l, m} = {j, k},

\begin{eqnarray} \Psi _l(X) & =& \sqrt{X^2 - \beta _l^{-2}} \arctan {\frac{\eta _l}{\sqrt{X^2 - \beta ^{-2}}}} - \eta _l, \end{eqnarray}

(26)

\begin{eqnarray} \Phi _{ll}(X) & =& \frac{\beta _l^{-2} - X^2}{i X} \arctan {\frac{i p }{X}} + p, \end{eqnarray}

(27)

\begin{eqnarray} \Phi _{lm}(X) &=& \Phi _{ml}(X) = \frac{1}{\beta _{\rm min}} \bigg\lbrace a R_F\left(c-1, c-b, c\right) \nonumber\\ && \left.-\, \frac{b}{3} R_D\left(c-1, c-b, c\right) + \frac{\left(a - 1\right) \left(a - b\right)}{3} R_J\right. \nonumber\\ &&\times\left(c-1, c-b, c, c-a\right) \bigg\rbrace , \end{eqnarray}

(28)

\begin{eqnarray} \Omega _{ll} & =& \frac{ \left[ \cos {\theta ^{\prime }} \Psi _l(X) - \sin {\theta ^{\prime }} \Phi _{ll}(X) \right]_{X=\sin {\theta ^{\prime }}/\beta _l}^{X=M} }{M^2 - (\sin {\theta ^{\prime }}/\beta _l)^2}, \end{eqnarray}

(29)

\begin{eqnarray} \Omega _{lm} & =& \frac{ \left[ \cos {\theta ^{\prime }} \Psi _l(X) - \sin {\theta ^{\prime }} \Phi _{lm}(X) \right]_{X=\sin {\theta ^{\prime }}/\beta _l}^{X=M} }{M^2 - (\sin {\theta ^{\prime }}/\beta _l)^2}. \end{eqnarray}

(30)

In eq. (28), we employ the following notation:

\begin{eqnarray} \beta _{\rm min} := \min \left\lbrace \beta _j, \beta _k\right\rbrace , \end{eqnarray}

(31)

\begin{eqnarray} a := 1/\left( \beta _{{\rm Max}} X \right)^2, \end{eqnarray}

(32)

\begin{eqnarray} b := \left( \beta _{{\rm min}} / \beta _{{\rm Max}} \right)^2, \end{eqnarray}

(33)

\begin{eqnarray} c & :=& 1/(\beta _{{\rm Max}} p)^2, \end{eqnarray}

(34)

where β_Max ≔ max{β_j, β_k}, and R_F, R_D and R_J are the Carlson symmetric forms of the elliptic integrals defined as

\begin{eqnarray} R_F(x,y,z) & =& \frac{1}{2} \int _0^\infty \frac{1}{\sqrt{(t+x) (t+y) (t+z)}} {\rm d}t, \end{eqnarray}

(35)

\begin{eqnarray} R_D(x,y,z) & =& \frac{3}{2} \int _0^\infty \frac{\sqrt{(t+x) (t+y)}}{(t+z) \sqrt{t+z}} {\rm d}t, \end{eqnarray}

(36)

\begin{eqnarray} R_J(x,y,z,p) & =& \frac{3}{2} \int _0^\infty \frac{1}{(t+p) \sqrt{(t+x) (t+y) (t+z)}} {\rm d}t. \end{eqnarray}

(37)

Even for complex variables, these integrals can be evaluated with arbitrary numerical accuracy following an algorithm proposed by Carlson (1995).

To calculate the wave fields according to the Cagniard–De Hoop method, the Cagniard path p, which is a complex function of (x, z, t; x^′, z^′, τ), has to be calculated. We can obtain the Cagniard path p by solving eqs (14) and (16). The Cagniard path p for the reflected wave is written as

\begin{eqnarray} p & = \left\lbrace \begin{array}{l@\quad l}0, & \left( t_d \le \frac{\left| z_d \right|}{\beta _j} \right) \\ \frac{x_d t_d - \textrm {sgn}(x_d) \left|z_d\right| \sqrt{\left(R/\beta _j\right)^2 - t_d^2}}{R^2}, & \left( \frac{\left| z_d \right|}{\beta _j} < t_d < \frac{\left| R \right|}{\beta _j} \right) \\ \frac{x_d t_d + i \textrm {sgn}(x_d) \left|z_d\right| \sqrt{t_d^2 - \left(R/\beta _j\right)^2}}{R^2}, & \left( t_d > \frac{\left| R \right|}{\beta _j} \right) \end{array} \right., \end{eqnarray}

(38)

where t_d = t − τ, x_d = x − x^′, z_d = z + z^′, and |$R = \sqrt{ x_d^2 + z_d^2}$|⁠. AR02 argues that the second line of eq. (38) holds over the entire range of 0 ≤ t_d < |R|/β_j. The Cagniard path p must, however, depart from the origin because of eq. (6.56) of AR02,

\begin{eqnarray*} G_R(x,z,s) &=& \frac{A}{2 \pi \rho \beta } \Im \int _0^{i\infty } \frac{1}{\eta _j} \left\lbrace \frac{\mu _j \eta _j - \mu _k \eta _k}{\mu _j \eta _j + \mu _k \eta _k} \right\rbrace \nonumber\\ &&\times\, {\rm e}^{-s \left( p x + \eta _j \left| z + z^{\prime } \right| \right)} {\rm d}p, \end{eqnarray*}

where A is the amplitude of the impulse, and x^′ = 0 and τ = 0 are assumed. If p departs not from the origin but from |$p = - (\textrm {sgn}(x_d) |z_d| \sqrt{(R/\beta _j)^2 - t_d^2})/R^2$|⁠, and the inequality |z_d|/(β_jR) ≥ 1/β_k is satisfied, interference with the branch cuts occurs even at t_d = 0, and causality is violated because p passes near one of the branch points ±1/β_k before passing the origin (Fig. 2).

Figure 2.

The Cagniard path of the reflected wave for x_d > 0 provided by AR02 (black arrow). Red and blue dotted lines indicate branch cuts due to η_j and η_k, respectively. AR02 showed only the right half plane of this figure in their fig. 6.15.

The Cagniard path for the transmitted wave obtained from eq. (16) is more complicated and is not given in AR02. Williams & Craster (2000) have pointed out that the Cagniard path for the transmitted wave is represented as a root of a quartic equation. By opening the square roots of eq. (16), we obtain the following quartic equation:

\begin{eqnarray} \left( \begin{array}{c}R^4 - Z^4 \\ - 4 t_d x_d R^2 \\ 2 R^2 T^2 + (2 t_d x_d)^2 + \left( \beta _j^{-2} + \beta _k^{-2} \right) Z^4 \\ - 4 t_d x_d T^2 \\ T^4 - \left( \frac{Z^2}{\beta _j \beta _k} \right)^2 \end{array} \right)^T \;\cdot\; \left( \begin{array}{c}p^4 \\ p^3 \\ p^2 \\ p \\ 1 \end{array} \right) = 0, \end{eqnarray}

(39)

where

\begin{eqnarray} R^2 & = x_d^2 + z^2 + z^{\prime 2}, \end{eqnarray}

(40)

\begin{eqnarray} Z^2 & = 2 z z^{\prime }, \end{eqnarray}

(41)

\begin{eqnarray} T^2 & =& t_d^2 - \left( \frac{z}{\beta _k} \right)^2 - \left( \frac{z^{\prime }}{\beta _j} \right)^2. \end{eqnarray}

(42)

The roots of a quartic equation can be obtained by certain algebraic algorithms (e.g. Ferrari's method). We must, however, choose only one root from among the four to determine a unique solution. Because all the coefficients of eq. (39) are real, the possible sets of roots are four real, two real and a pair of conjugates and two pairs of conjugates. The real roots do not contribute to the transmitted wave because their contribution appears only when wide-angle reflection occurs; thus, the transmitted wave is not observed for the real roots. Hence, it is sufficient to focus on the two pairs of complex conjugates. Even when complex roots exist, one of the conjugates is inappropriate because the imaginary part of p must be non-negative. Moreover, an extra root that does not satisfy eq. (16) exists because eq. (39) can be obtained not only from eq. (16) but also from t_d = px_d − η_j|z^′| + η_j|z|. After rejecting this, only one root remains to represent the wave field.

2.3 Verification of the convolution coefficients

As we show in the previous subsection and Appendix, the convolution coefficients for a bimaterial system are quite complicated, so here we verify the expression for the coefficients by comparing our analytical result with a numerical result obtained using the XBIEM (Kusakabe 2013). Because eq. (6) represents the stress response to a rectangular slip rate function given on one discretized crack element, we consider the crack element and calculate the stress response with these two methods. The configuration for our calculation is as follows:

• β⁻ = 1, β⁺ = 0.5, and ρ⁻ = ρ⁺ = 1, where the superscripts + and − denote z > 0 and z < 0, respectively.

• The source is of length 1 and inclined 45° counterclockwise from the x axis, and its midpoint is at (−5, −10).

• A constant slip velocity is given only at the time interval (τ₀, τ₁) = (0, 0.5).

• The observation points are located at (x, z) = (0, ±2) and (0, ±4).

Fig. 3 shows that the two results are almost identical, so we can conclude that our calculations of the coefficients are verified.

Figure 3.

Comparison of analytical result derived in this study (filled curves) and a numerical result calculated with XBIEM (black lines). (a) K_xy component and (b) K_zy component of stress field in a bimaterial system.

3 NUMERICAL SIMULATIONS OF DYNAMIC RUPTURE PROPAGATION

3.1 Model setup

Using the expression for the convolution coefficients derived in the previous section and an implicit time-marching scheme (Tada 2009), we numerically simulate dynamic rupture propagation along a crack crossing a material interface. We assume the crack has the following parameters and geometry, as illustrated in Fig. 4:

• The contrast of the S-wave velocities, Δβ ≔ (β⁺ − β⁻)/β⁻, is ±30 per cent. Moreover, the density of each medium ρ^± is proportional to β^±. This assumption is approximately valid with an error of less than 10 per cent for the Poisson solid (Birch 1961; Hirano & Yamashita 2011).

• In the half-plane z < 0, where the nucleation point of the rupture is embedded, the angle between the material interface (z = 0) and the orientation of the background antiplane shear stress, written as θ⁻, is assumed to be 90° or 45°.

• In the other half-plane, z > 0, the angle written as θ⁺ and the magnitude of the background stress are automatically determined by the relationships

  \begin{eqnarray} &&{\left. \begin{array}{rl}\sigma _{zy}^+ & = \sigma _{zy}^- \\ \frac{\sigma _{xy}^+}{\mu ^+} & = \frac{\sigma _{xy}^-}{\mu ^-} \end{array} \right\rbrace \Leftrightarrow \left\lbrace \begin{array}{rl}\sigma ^+\left(\theta ^+\right) & = \sigma ^-\left(\theta ^-\right) \sqrt{\cos ^2{\theta ^-} + \left( \frac{\mu ^+}{\mu ^-} \sin {\theta ^-} \right)^2} \\ \frac{\tan {\theta ^+}}{\mu ^+} & = \frac{\tan {\theta ^-}}{\mu ^-} \end{array} \right., }\nonumber\\ \end{eqnarray}

  (43)

  which are equivalent to the welded condition of two dissimilar materials (Rybicki & Yamashita 2008). In eq. (43), σ(θ) denotes the background antiplane shear stress acting on a plane inclined at the angle θ from the x axis.

• The crack consists of two straight segments oriented along the orientations of the background stress, that is the rotation angles θ^± from the interface in each medium.

• The crack is frictionless, so the stress drop is identical to the background stress, and the ratio of the stress drops on the two segments is determined by eq. (43).

• Unilateral rupture propagates at constant velocities |$V_r^\pm$| in each medium. We set |$V_r^+/\beta ^+ = V_r^-/\beta ^- = 90$| per cent for most cases. We set |$V_r^+/\beta ^+ = V_r^-/\beta ^- = 70$| per cent only for the last case to investigate the dependence of the cracking on the rupture velocity.

Figure 4.

Parameters and geometry of the crack crossing a material interface.

3.2 Perpendicular cases: θ^± = 90° and Δβ = ±30 per cent

The results for cracks perpendicular to the material interface, that is θ^± = 90°, show that the effect of the material interface on the cracking can be clearly seen and spatio-temporally localized. For both Δβ = ±30 per cent, snapshots of the amount of slip along the crack (blue lines in Fig. 5(a) for Δβ = −30 per cent and Fig. 6(a) for Δβ = +30 per cent) show approximately elliptical distributions; this result is quite similar to that for a static antiplane cracking perpendicular to a material interface (Bonafede et al.2002).

Figure 5.

Snapshots of behaviour of the crack perpendicular to the material interface. Δβ is assumed to be −30 per cent, that is a nucleation point is embedded in the stiffer medium. Red arrows indicate the material interface. See the text about the red curve in (a).

Figure 6.

Snapshots of behaviour of the crack perpendicular to the material interface. Δβ is assumed to be +30 per cent, that is a nucleation point is embedded in more compliant medium. Red arrows indicate the material interface. See the text about the red curve in (a).

The distributions are distorted, however, after the rupture front reaches the material interface. They can be investigated in detail in terms of the slip velocity, as shown in Fig. 5(b) for Δβ = −30 per cent and Fig. 6(b) for Δβ = +30 per cent. Figs 5(c) and 6(c) are horizontal cross sections of Figs 5(b) and 6(b), respectively. For both Δβ = ±30 per cent, the slip velocity distributions grow self-similarly before the rupture front reaches the material interface. Such self-similar unilateral solutions are quite similar to an analytical solution for unilateral singular cracking in a homogeneous medium (Nielsen & Madariaga 2003). After the rupture front reaches the interface, however, very sharp reflective phases of the slip velocity appear while the front continues to propagate; that is, the slip velocity exhibits a temporary increase for Δβ = −30 per cent and a temporary decrease for Δβ = +30 per cent. For Δβ = +30 per cent, surprisingly, the slip velocity temporarily takes a negative value. We discuss the effect of this phase on faulting in the next section.

We also investigate the change in the stress intensity near the rupture front as the crack extends. The red lines in Figs 5(a) and 6(a) indicate the shear stress on a discretized fault element in front of the rupture front. In other words, when the rupture front is located at x in the abscissa and the length of the discretized fault element is Δx, the red lines show shear stress at |$x+\frac{1}{2} \Delta x$|⁠. In both cases, the stress jumps when the rupture front passes the material interface. This jump is due to the discontinuity of the rigidity and stress drop. Except near an intersection of the crack and the interface, the change in the stress intensity is monotonic, and there is no remarkable transient phase. This might mean that dramatic changes in the stress intensity factor or energy release rate, if any, can occur only near the intersection. This tendency appears for the subsequent cases (Figs 7a, 8a, and 10a).

Figure 7.

Snapshots of behaviour of the crack inclined from the perpendicular case. Δβ is assumed to be −30 per cent, that is a nucleation point is embedded in the stiffer medium. Red arrows indicate the material interface. See the text about the red curve in (a).

Figure 8.

Snapshots of behaviour of the crack inclined from the perpendicular case. Δβ is assumed to be +30 per cent, that is a nucleation point is embedded in more compliant medium. Red arrows indicate the material interface. See the text about the red curve in (a).

3.3 Inclined cases: θ⁻ = 45° and Δβ = ±30 per cent

From eq. (43), we obtain θ⁺ = 19° for Δβ = −30 per cent and θ⁺ = 66° for Δβ = +30 per cent if θ⁻ = 45°. The numerical results for these two cases are shown in Figs 7 and 8, respectively. These results show that the reflective phases of the slip velocity become more diffusive in time than in the perpendicular cases; the arrival of these phases is generally earlier, and their amplitudes are obviously smaller than those observed in the perpendicular cases.

This diffusion in time is reasonable and comprehensible. Consider an arbitrary source–receiver pair on the lower crack segment inclined at an angle θ⁻ from the material interface. Let R be the distance from the intersection of the crack and the material interface to the source, and r be the distance from the intersection to the receiver (Fig. 9). Then the propagation distance of the reflected wave d is |$d = \sqrt{R^2 + r^2 - 2 R r \cos {2 \theta ^-}}$| according to the law of cosines. This is a decreasing function of θ⁻ (≤90°), so as θ⁻ decreases, the traveltime of the reflected wave decreases. Hence, the reflective phases appear earlier in the inclined cases.

Figure 9.

Propagation distance of the reflected wave between an arbitrary source–receiver pair on the lower crack segment.

3.4 Dependence on the rupture velocity: |$V_r^\pm /\beta ^\pm = 70$| per cent

Finally, we investigate the dependence of the cracking on the rupture velocity. In this case, |$V_r^\pm /\beta ^\pm$| is assumed to be 70 per cent, and the other parameters are the same as in the previous case: Δβ = +30 per cent and θ^± = 90° (Fig. 6). Although the amplitude of the reflective phase is smaller than that for a higher rupture velocity |$\left( V_r^\pm /\beta ^\pm = 90\,\mathrm{ per\ cent} \right)$|⁠, a negative slip velocity temporarily appears again. Hence, the material interface affects the cracking even for this slower rupture propagation.

4 DISCUSSION

4.1 Speculation regarding behaviour under rate-weakening friction regimes

Through our numerical simulations, we show that the behaviour of a crack crossing a material interface is obviously different from that of a crack embedded in a homogeneous medium even if a frictionless crack is assumed. Moreover, here we qualitatively discuss the possibility of further complex behaviour under rate-weakening friction. The evolution laws under a rate- and state-dependent friction regime have been widely investigated, mainly by laboratory rock experiments after Dieterich (1972) or Rice & Ruina (1983). Many of these experiments, however, have been done with extremely low slip rates (e.g. less than 1 mm s⁻¹). Experimental studies with high slip rates (e.g. ∼1 m s⁻¹) have shown mainly that the steady-state friction coefficients of rocks are between 0.0 and 0.4 (Di Toro et al.2011); they have not revealed any empirical laws that describe the evolution of friction. Some types of rate-dependent or rate- and state-dependent friction laws under high slip rates have been provided by rather theoretical studies (e.g. Rice 2006; Ampuero & Ben-Zion 2008; Noda et al.2009). Their formulations differ from each other, and the minor differences between them are outside of our scope, so we do not simulate cracking quantitatively by employing one of them. There exists, however, a significant common feature of slip behaviour universally obtained from many numerical simulations with those friction laws: the slip, except near the crack tips, stops growing after the slip rate drops below a cut-off level if the law is assumed to be one of rate weakening (Ampuero & Ben-Zion 2008; Noda et al.2009; Dunham et al.2011; Gabriel et al.2012). This feature should be inevitable because rate-weakening friction laws yield restrengthening of the slipping region when the slip rate decreases, and the restrengthening obviously reduces the slip rate. This positive feedback between the slip rate reduction and restrengthening could vary the slip behaviour, especially when Δβ = +30 per cent and the crack is perpendicular to the material interface, as shown in Fig. 6 or 10. In those cases, a negative slip rate even appears, so rapid restrengthening should occur, and slip growth may die out after the reflective phase passes if a rate-weakening friction law is applied. This means that the reflective phase is no longer a temporary reduction in the slip rate but becomes a backward-healing front, and the slip distribution can be strongly distorted. In contrast, the temporary acceleration of the slip observed in the opposite case, that is Δβ = −30 per cent, may be enhanced under a rate-weakening friction law because of positive feedback between the acceleration of the slip and weakening. The enhancement of the slip rate should change the slip distribution. Although quantitative simulations to verify such speculations are a future problem, we thus find a new possible mechanism that contributes to the complex behaviour of cracking.

Figure 10.

Snapshots of behaviour of the crack perpendicular to the material interface. The rupture velocity is assumed to be |$V_r^\pm /\beta ^\pm = 70$| per cent. Δβ is assumed to be +30 per cent, that is a nucleation point is embedded in more compliant medium. Red arrows indicate the material interface. See the text about the red curve in (a).

4.2 Hybrid methods of our BIEM and FEM or XBIEM

The actual crust is more heterogeneous than a pure bimaterial and is often modelled as a multilayered medium. If a crack is crossing at least two material interfaces, reflective phases like those seen in our simulations might appear several times because of multiple reflections between the interfaces. Hence, the behaviour of cracks embedded in such media might become more complicated and of great interest. Even though the newly developed BIEM in this paper is applicable only to a pure bimaterial system, we will be able to treat more inhomogeneous media, including a multilayered medium, by combining the BIEM with the FEM or XBIEM. We here propose ideas for such hybrid methods.

Goto et al. (2010) proposed a method of combining the FEM and a conventional BIEM for an infinite homogeneous elastic medium to treat a crack embedded in a semi-infinite homogeneous elastic medium. Although wave propagation in inhomogeneous media can be calculated in a FEM framework, the numerical precision of the solutions decreases if the propagation distance of the wave is not sufficiently longer than the mesh size of the FEM. This means that the stress concentration near the rupture front cannot be calculated precisely by using the FEM. On the other hand, a relatively precise solution near the rupture front can be obtained by using the BIEM, as mentioned in the introduction. The procedure proposed by Goto et al. (2010) is as follows. First, they calculated FEM solutions for both infinite homogeneous and semi-infinite homogeneous media. Secondly, they subtracted the former from the latter so that the remainder is the field of only the reflected wave. Finally, they added a solution obtained by using the conventional BIEM to the remainder. In this method, the direct wave propagating a shorter distance is simulated using the BIEM, and the reflected wave propagating a longer distance is simulated using the FEM, so both the difficulty of treating an inhomogeneous medium in the BIEM and the low precision of the FEM solution near the rupture front are overcome. Their method does not work well, however, if a crack is crossing the material interfaces because, in this case, the propagation distance of the reflected or transmitted wave between two arbitrary points of a crack will also become quite short when the points become closer to the material interface. Hence, it again becomes difficult to simulate an interaction between such points by the FEM. Our BIEM based on the exact solution for a bimaterial, however, can be applied to such an interaction with high precision. Hence, a possible hybrid method is as follows. If a crack is embedded in a multilayered medium, we have to calculate the interactions between the crack elements. The direct wave propagating only in each medium can be calculated using the conventional BIEM. The reflected and transmitted waves scattered by a layer boundary only once can be calculated with our BIEM. In addition, multiple reflections or multiple refractions can be calculated using the FEM because their propagation distance is obviously longer than the thickness of each layer. Then the summation of the solutions produced by these three calculations represents the wave field due to cracking in a multilayered medium.

As described in the introduction, the XBIEM developed by Kame & Kusakabe (2012) and Kusakabe (2013) is a promising scheme for considering sectionally homogeneous elastic media. The difficulty of the XBIEM, however, is the huge number of discretized elements aligned along every material interface; the number increases the calculation time, and the discretization reduces the numerical precision. This difficulty can be partially removed by combining the XBIEM and our BIEM. As an example, let us consider a four-layered medium that has different material constants in four regions: −∞ < z < −h, −h < z < 0, 0 < z < +H, and +H < z < +∞. In this case, three integrals of the traction terms are required along z = −h, z = 0, and z = +H in the XBIEM framework. At present, we can consider a bimaterial system in our BIEM framework, so we can regard the four-layered medium as a ‘two-welded bimaterial’ system (Fig. 11). Hence, the interaction between the two regions −∞ < z < −h and −h < z < 0 can be calculated with our BIEM, and, of course, the interaction between 0 < z < +H and +H < z < +∞ can also be handled in the same way. As a result, the integral of the traction terms is required only on z = 0, and the number of discretized elements is dramatically reduced.

Figure 11.

Schematic illustration of a numerical model of crack embedded in a four-layered medium. Number of discretized elements on interfaces is dramatically reduced by using a hybrid scheme of XBIEM and our BIEM.

5 CONCLUSIONS

We developed a BIEM for analysing dynamic propagation of an antiplane rupture crossing a material interface. Because an exact representation of the displacement Green's function associated with the reflected and transmitted waves was provided as a function of the complex ray parameter, we performed a complex analysis to obtain an exact solution of the stress response due to a constant slip rate on a discretized crack element. Our exact solution can be used as a benchmark for other numerical schemes and contributes to the expansion of conventional BIEM analysis.

Using the BIEM developed in this paper, we numerically simulated dynamic propagation of a rupture crossing a material interface. We showed that the wave reflected from the material interface generates complexity in the slip velocity on the crack segment located in a half-plane, where the rupture is nucleated after the rupture front reaches the material interface; the slip velocity is temporarily enhanced when the nucleation point is in the more compliant medium and temporarily reduced, even to negative values, in the opposite case. We found that these phenomena can be seen most clearly when the crack is perpendicular to the interface and become diffusive when the crack crosses the interface obliquely. This tendency can be understood from the relation between the orientation of the crack relative to that of the interface and the propagation distance of the reflected wave.

Although our simulations were done assuming a frictionless crack, the employment of more realistic friction laws suggested by the fault dynamics could vary the behaviour observed in our simulations. Many authors have demonstrated that under rate-weakening friction, positive feedback occurs between the deceleration of the slip and healing on the fault surface. Because of this feedback, the wave reflected from the material interface could freeze the slip evolution when the initiation point of the rupture is embedded in the more compliant medium. This suggests that the slip behaviour, such as the slip duration and final slip distribution, depends on the location of the material interface relative to that of the rupture nucleation point. Note that this qualitative speculation is based not on a specific friction law but on a property universally suggested by various rate-weakening friction models.

Finally, we proposed the possibility of new schemes to treat more complicated heterogeneity with high accuracy. By combining our BIEM and the FEM or XBIEM, we will be able to simulate the behaviour of a crack that crosses multiple material interfaces and interacts with multiple wave reflections. Because the crust is quite heterogeneous, such a simulation will be necessary in the future, and our proposed method could contribute to meeting this need.

We thank Satoshi Ide, Nobuki Kame, Tetsuya Kusakabe and two anonymous reviewers for useful discussions. The XBIEM solution plotted in Fig. 3 was provided by Tetsuya Kusakabe. We acknowledge partial financial support by the MegaQuake project, University of Tsukuba.

APPENDIX: DERIVATION OF I_xy,R, I_zy,R, I_xy,T AND I_zy,T

In this appendix, we employ the same notation as in Section 2. We here describe some algebra for obtaining the representation of I_xy,R written as eq. (18). I_zy,R, I_xy,T and I_zy,T are obtained in a similar manner to the derivation of I_xy,R. In the following calculations, all of the integration constants are neglected because the integrals are ultimately considered as definite integrals, as in eq. (11).

By substituting eq. (13) into eq. (7), the multiple integral appearing in the RHS of (7) can be reduced to single integrals as follows:

\begin{eqnarray} I_{xy,R} & =& \int {\rm d}\gamma \int {\rm d}\tau \int {\rm d}\tau \frac{\partial ^2}{\partial x \partial \zeta } G_R \end{eqnarray}

(A1)

\begin{eqnarray} = \int {\rm d}\gamma \int {\rm d}\tau \int {\rm d}\tau \frac{\partial }{\partial \zeta } \left( \cos {\theta } \frac{\partial }{\partial \xi } - \sin {\theta } \frac{\partial }{\partial \zeta } \right) G_R \end{eqnarray}

(A2)

\begin{eqnarray} =\int {\rm d}\gamma \!\int {\rm d}\tau \!\int {\rm d}\tau \left\lbrace \cos {\theta } \frac{\partial ^2}{\partial \xi \partial \zeta } \!-\! \sin {\theta } \left( - \frac{\partial ^2}{\partial \xi ^2} \!+\! \beta _j^{-2} \frac{\partial ^2}{\partial t^2} \right) \right\rbrace G_R \end{eqnarray}

(A3)

\begin{eqnarray} = \int {\rm d}\gamma \int {\rm d}\tau \int {\rm d}\tau \left\lbrace \frac{\partial ^2}{\partial \xi \partial z} - \frac{\sin {\theta }}{\beta _j^2} \frac{\partial ^2}{\partial \ \tau ^2} \right\rbrace G_R \end{eqnarray}

(A4)

\begin{eqnarray} = - \int {\rm d}\gamma \int {\rm d}\tau \int {\rm d}\tau \left\lbrace \textrm {sgn}\left(z^{\prime }\right) \frac{\partial ^2}{\partial \xi \partial \tau } \eta _j + \frac{\sin {\theta }}{\beta _j^2} \frac{\partial ^2}{\partial \ \tau ^2} \right\rbrace G_R \end{eqnarray}

(A5)

\begin{eqnarray} = - \textrm {sgn}\left(z^{\prime }\right) \left\lbrace \int \eta _j G_R {\rm d}\tau + \frac{\sin {\theta ^{\prime }}}{\beta _j^2} \int G_R {\rm d}\xi \right\rbrace . \end{eqnarray}

(A6)

In eq. (A4), we use the rotational symmetry of the Laplacian, that is |$\beta _j^{-2} \partial _t^2 G_R = ( \partial _x^2 + \partial _z^2 ) G_R = ( \partial _\xi ^2 + \partial _\zeta ^2 ) G_R$|⁠. Moreover, in eq. (A5), we use the relationship

\begin{eqnarray} \frac{\partial G_R}{\partial z} & =& \frac{\partial G_R}{\partial z^{\prime }} = \textrm {sgn}(z^{\prime }) \frac{\partial }{\partial t} (\eta _j G_R) \end{eqnarray}

(A7)

obtained from eq. (14); this is rewritten as |$\partial _z = \textrm {sgn}\left(z+z^{\prime }\right) s \sqrt{\beta ^{-2} - \left( - i k_x / s \right)^2}$| in the Fourier–Laplace domain considered in the Cagniard–De Hoop method, where k_x is the wavenumber along the x axis [AR02]. For I_zy,R, instead, we use the relationship

\begin{eqnarray} \frac{\partial G_R}{\partial x} & =& - \frac{\partial G_R}{\partial x^{\prime }} = \frac{\partial }{\partial t} (p G_R) \end{eqnarray}

(A8)

obtained from eq. (13); this is rewritten as −ik_x = sp in the Fourier–Laplace domain. For I_xy,T and I_zy,T,

\begin{eqnarray} \frac{\partial G_T}{\partial x} & =& - \frac{\partial G_T}{\partial x^{\prime }} = \frac{\partial }{\partial t} \left( p G_T \right), \end{eqnarray}

(A9)

\begin{eqnarray} \frac{\partial G_T}{\partial z} & =& \textrm {sgn}\left(z\right) \frac{\partial }{\partial t} \left( \eta _k G_T \right), \end{eqnarray}

(A10)

\begin{eqnarray} \frac{\partial G_T}{\partial z^{\prime }} & =& \textrm {sgn}\left(z^{\prime }\right) \frac{\partial }{\partial t} \left( \eta _j G_T \right) \end{eqnarray}

(A11)

are used.

Because eq. (13) is a function of the Laplace domain, G_R in the time domain can be written as follows:

\begin{eqnarray} G_R & =& \tilde{G}_R \frac{{\rm d} p}{{\rm d} t}, \end{eqnarray}

(A12)

\begin{eqnarray} \tilde{G}_R & =& \frac{1}{2 \pi \mu _j \eta _j} \left\lbrace \frac{\mu _j \eta _j - \mu _k \eta _k}{\mu _j \eta _j + \mu _k \eta _k} \right\rbrace . \end{eqnarray}

(A13)

By employing this notation and because dt = −dτ, the first integral of eq. (A6) is reduced to

\begin{eqnarray} \int \eta _j G_R\, {\rm d}\tau & =& - \int \eta _j \tilde{G}_R\, {\rm d}p. \end{eqnarray}

(A14)

The second integral of eq. (A6) can also be represented as an integral of p as follows. G_R is rewritten as the following two integrals:

\begin{eqnarray} G_R = \int \frac{\partial }{\partial \xi } G_R\, {\rm d}\xi = \frac{\partial }{\partial t} \int \left( -p \cos {\theta ^{\prime }} + \eta _j \sin {\theta ^{\prime }} \right) G_R \,{\rm d}\xi \end{eqnarray}

(A15)

from eqs (A8) and (A7), and

\begin{eqnarray} G_R = \tilde{G}_R \frac{{\rm d}p}{{\rm d}t} = \frac{\partial }{\partial t} \int \tilde{G}_R {\rm d}p = \frac{\partial }{\partial t} \int \tilde{G}_R \frac{{\rm d}p}{{\rm d}\xi } {\rm d}\xi . \end{eqnarray}

(A16)

By comparing the integrands of these two integrals, we obtain

\begin{eqnarray} G_R &=& \frac{\tilde{G}_R}{- p \cos {\theta ^{\prime }} + \eta _j \sin {\theta ^{\prime }}} \frac{{\rm d}p}{{\rm d}\xi } = \frac{p \cos {\theta ^{\prime }} + \eta _j \sin {\theta ^{\prime }}}{\left( \sin {\theta ^{\prime }} / \beta _j \right)^2 - p^2} \tilde{G}_R \frac{{\rm d}p}{{\rm d}\xi }, \end{eqnarray}

(A17)

that is

\begin{eqnarray} \int G_R \,{\rm d}\xi = \int \frac{p \cos {\theta ^{\prime }} + \eta _j \sin {\theta ^{\prime }}}{\left( \sin {\theta ^{\prime }} / \beta _j \right)^2 - p^2} \tilde{G}_R {\rm d}p. \end{eqnarray}

(A18)

The material constants defined by eq. (22) reduce eq. (A13) to

\begin{eqnarray} \tilde{G}_R &=& \frac{1}{\pi } \left\lbrace \frac{1}{\mu _j \eta _j + \mu _k \eta _k} - \frac{1}{2 \mu _j \eta _j} \right\rbrace = \frac{1}{\pi } \left\lbrace \frac{C_j \eta _j + C_k \eta _k}{M^2 - p^2} - \frac{1}{2 \mu _j \eta _j} \right\rbrace . \end{eqnarray}

(A19)

By substituting eq. (A19) into eqs (A14) and (A18), we find that eq. (A6) can be represented as a linear combination of the following:

\begin{eqnarray} \int \frac{p \cos {\theta ^{\prime }} + \eta _j \sin {\theta ^{\prime }}}{\left(\sin {\theta ^{\prime }}/\beta _j\right)^2 - p^2} \frac{{\rm d}p}{\eta _j} &=& \sin {\theta ^{\prime }} \int \frac{\eta _j}{\left(\sin {\theta ^{\prime }}/\beta _j\right)^2 - p^2} \frac{{\rm d}p}{\eta _j} +\, \cos {\theta ^{\prime }} \int \frac{p}{\eta _j^2 - \left(\cos {\theta ^{\prime }}/\beta _j\right)^2} \frac{{\rm d}p}{\eta _j} \end{eqnarray}

(A20)

\begin{eqnarray} = \beta _j \left\lbrace \textrm {arctanh}{\frac{\beta _j p}{\sin {\theta ^{\prime }}}} + \textrm {arctanh}{\frac{\beta _j \eta _j}{\cos {\theta ^{\prime }}}} \right\rbrace \end{eqnarray}

(A21)

\begin{eqnarray} = \beta _j^2 A, \end{eqnarray}

(A22)

\begin{eqnarray} \int \frac{p \eta _l}{X^2 - p^2} {\rm d}p & =& \int \frac{\eta _l^2}{X^2 - \beta _l^{-2} - \eta _l^2} \frac{p}{\eta _l} {\rm d}p = \int \frac{\eta _l^2}{X^2 - \beta _l^{-2} - \eta _l^2} {\rm d} \eta _l = \Psi _l(X), \end{eqnarray}

(A23)

\begin{eqnarray} \int \frac{\eta _l^2}{X^2 - p^2} {\rm d}p & =& \int \frac{\beta _l^{-2} - X^2 + X^2 - p^2}{X^2 - p^2} {\rm d}p = \Phi _{ll}(X), \end{eqnarray}

(A24)

\begin{eqnarray} \int \frac{\eta _l \eta _m}{X^2 - p^2} {\rm d}p & =& \frac{1}{\beta _{{\rm min}}} \int \frac{\sqrt{1 - q^2} \sqrt{1 - b q^2}}{1/a - q^2} {\rm d}q \end{eqnarray}

(A25)

\begin{eqnarray} = \frac{1}{\beta _{{\rm min}}} \left\lbrace (a-b) \Pi \left(q,\sqrt{b},a\right) + b F\left(q,\sqrt{b}\right) + E\left(q,\sqrt{b}\right) \right\rbrace, \end{eqnarray}

(A26)

where q ≔ β_minp, and

\begin{eqnarray} F\left(q,\sqrt{b} \right) & =& \int \frac{{\rm d}q}{\sqrt{1-q^2} \sqrt{1-bq^2}}, \end{eqnarray}

(A27)

\begin{eqnarray} E\left(q,\sqrt{b} \right) & =& \int \frac{\sqrt{1-bq^2}}{\sqrt{1-q^2}} {\rm d}q, \end{eqnarray}

(A28)

\begin{eqnarray} \Pi \left(q,\sqrt{b},a \right) & =& \int \frac{{\rm d}q}{\left( 1 - a q^2 \right) \sqrt{1-q^2} \sqrt{1-bq^2}} \end{eqnarray}

(A29)

are incomplete elliptic integrals of the first, second, and third kind in Legendre canonical forms, respectively. According to Carlson (1995), there exist the following relationships:

\begin{eqnarray} F\left( q,\sqrt{b} \right) & =& R_F\left( c-1,c-b,c \right), \end{eqnarray}

(A30)

\begin{eqnarray} E\left( q,\sqrt{b} \right) & =& R_F\left( c-1,c-b,c \right) - \frac{b}{3} R_D\left( c-1, c-b, c \right), \end{eqnarray}

(A31)

\begin{eqnarray} \Pi \left( q,\sqrt{b},a \right)\! & =&\! R_F\left( c\!-\!1,c\!-\!b,c \right) \!-\! \frac{a}{3} R_J\left( c\!-\!1, c\!-\!b, c, c\!-\!a \right). \end{eqnarray}

(A32)

Hence,

\begin{eqnarray} \int \frac{\eta _l \eta _m}{X^2 - p^2} {\rm d}p = \Phi _{lm}(X) \end{eqnarray}

(A33)

holds, and, finally, we obtain eq. (18).