Green’s function representations for Marchenko imaging without up/down decomposition

Wapenaar, Kees; Snieder, Roel; de Ridder, Sjoerd; Slob, Evert

doi:10.1093/gji/ggab220

SUMMARY

Marchenko methods are based on integral representations which express Green’s functions for virtual sources and/or receivers in the subsurface in terms of the reflection response at the surface. An underlying assumption is that inside the medium the wave field can be decomposed into downgoing and upgoing waves and that evanescent waves can be neglected. We present a new derivation of Green’s function representations which circumvents these assumptions, both for the acoustic and the elastodynamic situation. These representations form the basis for research into new Marchenko methods which have the potential to handle refracted and evanescent waves and to more accurately image steeply dipping reflectors.

Controlled source seismology, Seismic interferometry, Wave scattering and diffraction

1 INTRODUCTION

Marchenko redatuming, imaging, monitoring and multiple elimination are all derived from integral representations which express Green’s functions for virtual sources or receivers in the subsurface in terms of the reflection response at the surface (Ravasi et al. 2016; Jia et al. 2018; Staring et al. 2018; Brackenhoff et al. 2019; Lomas & Curtis 2019; Mildner et al. 2019; Elison et al. 2020; Reinicke et al. 2020; Zhang & Slob 2020). These representations, in turn, are derived from reciprocity theorems for one-way wave fields (Slob et al. 2014; Wapenaar et al. 2014), building on ideas presented by Broggini & Snieder (2012). Marchenko methods deal with internal multiples in a data-driven way and have the potential to solve large-scale 3-D imaging and multiple elimination problems (Pereira et al. 2019; Staring & Wapenaar 2020; Ravasi & Vasconcelos 2021). Of course Marchenko methods have also limitations. One of the limitations is caused by the fact that the one-way reciprocity theorems require that the wave field in the subsurface region of interest can be decomposed into downgoing and upgoing fields. Moreover, one of these reciprocity theorems (the correlation-type theorem) is based on the assumption that evanescent waves can be neglected. These assumptions complicate the imaging of steep flanks and exclude a proper treatment of refracted waves and evanescent waves tunnelling through high velocity layers.

To address some of the limitations, Kiraz et al. (2021) propose a Marchenko method without decomposition inside the medium, assuming the input data are acquired on a closed boundary. On the other hand, for reflection data on a single horizontal boundary, a first step has been made towards a Marchenko method that deals with evanescent waves (Wapenaar 2020). This method is restricted to horizontally layered media and uses wave field decomposition inside the medium.

In this paper, we derive more general Green’s function representations which do not rely on wave field decomposition in the subsurface and which hold for an arbitrarily inhomogeneous medium below a single horizontal acquisition boundary. These representations form a starting point for new research on Marchenko methods which circumvent several of the present limitations. Diekmann & Vasconcelos (2021) independently investigate the same problem, but without specifying a focusing condition for their focusing function f and requiring a time-symmetric source function for f. Our derivation follows a different approach, using an explicit focusing condition and requiring no source function for our focusing function f. Moreover, we derive several forms of Green’s function representations, including one for the homogeneous Green’s function between a virtual source and a virtual receiver in the subsurface. We also derive elastodynamic versions of these representations.

This paper is restricted to the derivation of the Green’s function representations; a discussion of their application in new Marchenko methods is beyond the scope of this paper.

2 ACOUSTIC WAVE FIELD REPRESENTATION

We consider a lossless acoustic medium, consisting of a homogeneous isotropic upper half-space and an arbitrary inhomogeneous anisotropic lower half-space, separated by a horizontal surface |${{\partial \mathbb {D}}_R}$|⁠. Coordinates in the medium are denoted by x = (x_H, x₃), with x_H = (x₁, x₂) denoting the horizontal coordinates and x₃ the depth coordinate (the positive x₃-axis is pointing downward). The horizontal surface |${{\partial \mathbb {D}}_R}$| is defined at x₃ = x_{3, R} (in the next section we choose this as the surface at which seismic acquisition takes place). The medium parameters of the lower half-space x₃ > x_{3, R} are the compressibility κ(x) and the mass density tensor ρ_jk(x). At the micro scale (much smaller than the wavelength of the acoustic field) the mass density is isotropic. However, small-scale heterogeneities of the isotropic mass density, for example caused by fine-layering, may manifest themselves as effective anisotropy at the scale of the wavelength (Schoenberg & Sen 1983). The mass density tensor is symmetric, that is, ρ_jk(x) = ρ_kj(x). The parameters of the upper half-space x₃ < x_{3, R} are the constant compressibility κ = κ₀ and the constant isotropic mass density ρ_jk = δ_jkρ₀, where δ_jk is the Kronecker delta function. The propagation velocity of the upper half-space is c₀ = (κ₀ρ₀)^−1/2. At |${{\partial \mathbb {D}}_R}$| we choose the same constant isotropic medium parameters as in the upper half-space.

The basic equations for acoustic wave propagation are the linearized equation of motion

$$\begin{eqnarray*} \rho _{jk}\partial _tv_k+\partial _jp=0 \end{eqnarray*}$$

(1)

and the linearized deformation equation

$$\begin{eqnarray*} \kappa \partial _tp+\partial _iv_i=q, \end{eqnarray*}$$

(2)

respectively. Here p(x, t) is the space (x) and time (t) dependent acoustic pressure, |$v_i({\bf x},t)$| the particle velocity and q(x, t) a source in terms of volume-injection rate density. Operator ∂_i stands for differentiation in the x_i-direction. Lower-case subscripts (except t) take on the values 1, 2 and 3, and the summation convention applies to repeated subscripts. Operator ∂_t stands for differentiation with respect to time. We introduce the specific volume tensor ϑ_ij(x) as the inverse of the mass density tensor, with ϑ_ijρ_jk = δ_ik. Applying the operator ∂_iϑ_ij to eq. (1), operator ∂_t to eq. (2), and subtracting the two equations yields the acoustic wave equation

$$\begin{eqnarray*} \partial _i(\vartheta _{ij}\partial _jp)-\kappa \partial _t^2p=-\partial _t q. \end{eqnarray*}$$

(3)

We introduce a focusing function F(x, x_R, t), in which x_R = (x_{H, R}, x_{3, R}) denotes the position of a focal point at |${{\partial \mathbb {D}}_R}$| (Fig. 1). For fixed x_R and variable x and t, this focusing function is a solution of wave eq. (3) for the source-free situation, hence, for q = 0. We define the focusing condition as

$$\begin{eqnarray*} F({\bf x},{\bf x}_R,t)|_{x_3=x_{3,R}}=\delta ({\bf x}_{\rm H}-{\bf x}_{{\rm H},R})\delta (t), \end{eqnarray*}$$

(4)

and further demand that F(x, x_R, t) is purely upgoing at |${{\partial \mathbb {D}}_R}$| and in the homogeneous isotropic upper half-space. Note that F(x, x_R, t) is similar, but not identical, to the focusing function f₂(x, x_R, t) introduced in Wapenaar et al. (2014). We come back to this in Section 3.2.

$Illustration of the focusing function F(x, xR, t), which focuses at xR. For x at and above ${{\partial \mathbb {D}}_R}$ it is purely upgoing. For x in the half-space below ${{\partial \mathbb {D}}_R}$ it is a complex wave field. The near-horizontal arrows in the thin layer illustrate tunnelling evanescent waves.$

Figure 1.

Illustration of the focusing function F(x, x_R, t), which focuses at x_R. For x at and above |${{\partial \mathbb {D}}_R}$| it is purely upgoing. For x in the half-space below |${{\partial \mathbb {D}}_R}$| it is a complex wave field. The near-horizontal arrows in the thin layer illustrate tunnelling evanescent waves.

Open in new tab Download slide

We define the temporal Fourier transform of a space- and time-dependent function u(x, t) as

$$\begin{eqnarray*} u({\bf x},\omega )=\int _{-\infty }^\infty u({\bf x},t)\exp (i\omega t){\rm d}t, \end{eqnarray*}$$

(5)

where ω is the angular frequency and i the imaginary unit. The integral is taken from t = −∞ to t = ∞ to account for non-causal functions, such as the focusing function F(x, x_R, t) for which no causality condition is implied (hence, in general it can be non-zero for positive and negative time). With this transform, wave eq. (3) transforms to

$$\begin{eqnarray*} \mathcal{L}p=i\omega q, \end{eqnarray*}$$

(6)

with

$$\begin{eqnarray*} {\cal L}=\partial _i\vartheta _{ij}\partial _j+\omega ^2\kappa . \end{eqnarray*}$$

(7)

The focusing function F(x, x_R, ω) obeys in the frequency domain the wave equation

$$\begin{eqnarray*} {\cal L}F=0, \end{eqnarray*}$$

(8)

the focusing condition

$$\begin{eqnarray*} F({\bf x},{\bf x}_R,\omega )|_{x_3=x_{3,R}}=\delta ({\bf x}_{\rm H}-{\bf x}_{{\rm H},R}), \end{eqnarray*}$$

(9)

and it is upgoing at and above |${{\partial \mathbb {D}}_R}$|⁠. We discuss a representation for a wave field p(x, ω), which may have sources in the upper half-space above |${{\partial \mathbb {D}}_R}$|⁠, but which obeys the source-free wave equation |${\cal L}p=0$| for x₃ ≥ x_{3, R}. In the lower half-space we express p(x, ω) as a superposition of mutually independent wave fields that obey the same source-free wave equation as p(x, ω) for x₃ ≥ x_{3, R}. For this purpose, we choose the focusing functions F(x, x_R, ω) and F*(x, x_R, ω) (the asterisk denotes complex conjugation, which corresponds to time-reversal in the time domain). To be more specific, we express p(x, ω) as

$$\begin{eqnarray*} p({\bf x},\omega )&=&\int _{{{\partial \mathbb {D}}_R}} F({\bf x},{\bf x}_R,\omega )a({\bf x}_R,\omega ){\rm d}{\bf x}_R +\int _{{{\partial \mathbb {D}}_R}} F^*({\bf x},{\bf x}_R,\omega )b({\bf x}_R,\omega ){\rm d}{\bf x}_R, \quad\mbox{for}\quad x_3\ge x_{3,R}. \end{eqnarray*}$$

(10)

Here a(x_R, ω) and b(x_R, ω) are as yet undetermined coefficients, which depend on the position x_R at |${{\partial \mathbb {D}}_R}$|⁠. In Appendix A1 we formulate boundary conditions for the acoustic pressure and the vertical component of the particle velocity at |${{\partial \mathbb {D}}_R}$|⁠, from which we solve a(x_R, ω) and b(x_R, ω). We thus obtain

$$\begin{eqnarray*} p({\bf x},\omega )&=&\int _{{{\partial \mathbb {D}}_R}} F({\bf x},{\bf x}_R,\omega )p^-({\bf x}_R,\omega ){\rm d}{\bf x}_R +\int _{{{\partial \mathbb {D}}_R}} F^*({\bf x},{\bf x}_R,\omega )p^+({\bf x}_R,\omega ){\rm d}{\bf x}_R,\quad \mbox{for}\quad x_3\ge x_{3,R}, \end{eqnarray*}$$

(11)

where p⁻(x_R, ω) and p⁺(x_R, ω) represent the upgoing (−) and downgoing (+) parts, respectively, of p(x_R, ω) for x_R at |${{\partial \mathbb {D}}_R}$|⁠. These upgoing and downgoing fields are pressure-normalized, meaning that p⁻ + p⁺ = p at and above |${{\partial \mathbb {D}}_R}$|⁠. Below |${{\partial \mathbb {D}}_R}$| we only consider the total (undecomposed) wave field p.

For an intuitive explanation of the right-hand side of eq. (11) we refer to Fig. 2. First we consider the second integral, which is illustrated in Fig. 2(b). The downgoing field p⁺(x_R, ω) is incident from the homogeneous upper half-space to |${{\partial \mathbb {D}}_R}$|⁠. The complex conjugate focusing function F*(x, x_R, ω) propagates this downgoing field from x_R at |${{\partial \mathbb {D}}_R}$| to x in the lower half-space and the integral superposes the contributions from all x_R at |${{\partial \mathbb {D}}_R}$|⁠. Next we consider the first integral of eq. (11), which is illustrated in Fig. 2(a). Here the upgoing field p⁻(x_R, ω) is backpropagated by the focusing function F(x, x_R, ω) from all x_R at |${{\partial \mathbb {D}}_R}$| to x in the lower half-space. The sum of the two integrals gives the wave field p(x, ω) in the lower half-space. This is a modified form of Huygens’ principle, with the focal points x_R denoting the positions of the secondary sources at |${{\partial \mathbb {D}}_R}$|⁠, and with the forward and backward propagating focusing functions replacing the Green’s functions in the usual form of Huygens’ principle. Note that the two integrals cannot be separately associated with p⁺(x, ω) and p⁻(x, ω) in the lower half-space; only the sum of the two integrals gives the total field p(x, ω), according to eq. (11).

$Explanation of the first (a) and second (b) integral in eq. (11) in terms of Huygens’ principle. The focusing functions F(x, xR, ω) and F*(x, xR, ω) propagate the fields p−(xR, ω) and p+(xR, ω) from all xR at the surface ${{{\partial \mathbb {D}}_R}}$ to x in the lower half-space. The superposition of these propagated fields gives the field p(x, ω) for any x in the lower half-space.$

Figure 2.

Explanation of the first (a) and second (b) integral in eq. (11) in terms of Huygens’ principle. The focusing functions F(x, x_R, ω) and F*(x, x_R, ω) propagate the fields p⁻(x_R, ω) and p⁺(x_R, ω) from all x_R at the surface |${{{\partial \mathbb {D}}_R}}$| to x in the lower half-space. The superposition of these propagated fields gives the field p(x, ω) for any x in the lower half-space.

Open in new tab Download slide

The underlying assumption in the derivation of eq. (11) is that evanescent waves can be neglected at |${{\partial \mathbb {D}}_R}$|⁠. Hence, it only holds for waves that are propagating at |${{\partial \mathbb {D}}_R}$|⁠, having a horizontal slowness s obeying

$$\begin{eqnarray*} |{\bf s}|\le 1/c_0, \quad \mbox{at}\, \, {{\partial \mathbb {D}}_R}. \end{eqnarray*}$$

(12)

This implies that the foci (and hence the secondary sources) at |${{\partial \mathbb {D}}_R}$| are not ideal delta functions [as formulated by eq. (9)] but band-limited versions of delta functions. Note that ignoring evanescent waves at |${{\partial \mathbb {D}}_R}$| does not imply that evanescent waves are not accounted for inside the inhomogeneous medium below |${{\partial \mathbb {D}}_R}$|⁠. For example, in an isotropic horizontally layered medium with depth-dependent velocity c(x₃), in which the horizontal slowness is independent of depth, waves that are propagating at |${{\partial \mathbb {D}}_R}$| become evanescent when they reach a depth at which 1/c(x₃) < |s| ≤ 1/c₀. In Section 3.4, we show with a numerical example that eq. (11) indeed accounts for such evanescent waves. Although for laterally varying media we cannot formulate a similar precise condition for waves becoming evanescent, it is still true that eq. (11) holds for evanescent waves inside the medium, as long as they are related to propagating waves at the surface, as formulated by eq. (12).

Note that we previously derived a representation similar to eq. (11) with heuristic arguments, and used it as the starting point for deriving the Marchenko method (Wapenaar et al. 2013). However, further on in that derivation we applied up/down decomposition to the wave field at an artificial internal boundary in the lower half-space and we neglected evanescent waves throughout space. In the following derivations we avoid up/down decomposition in the lower half-space and evanescent waves are only neglected at |${{\partial \mathbb {D}}_R}$|⁠. From eq. (11) we derive Green’s function representations for the full wave field at any point x in the subsurface, expressed in terms of the reflection response at the surface.

3 ACOUSTIC GREEN’S FUNCTION REPRESENTATIONS

3.1 Representation for the acoustic dipole Green’s function

We introduce the Green’s function G(x, x_S, t) as a solution of eq. (3) for an impulsive monopole source of volume-injection rate density at x_S, hence

$$\begin{eqnarray*} \partial _i(\vartheta _{ij}\partial _jG)-\kappa \partial _t^2G=-\delta ({\bf x}-{\bf x}_S)\partial _t\delta (t). \end{eqnarray*}$$

(13)

We demand that G is the causal solution of this equation, hence G(x, x_S, t) = 0 for t < 0. Note that G obeys source–receiver reciprocity, that is G(x, x_S, t) = G(x_S, x, t). In the frequency domain, G(x, x_S, ω) obeys the following wave equation

$$\begin{eqnarray*} {\cal L} G=i\omega \delta ({\bf x}-{\bf x}_S). \end{eqnarray*}$$

(14)

We choose x_S = (x_{H, S}, x_{3, S}) in the upper half-space, at a vanishing distance ϵ above |${{\partial \mathbb {D}}_R}$|⁠, hence, x_{3, S} = x_{3, R} − ϵ. We define a dipole-source response as

$$\begin{eqnarray*} \Gamma ({\bf x},{\bf x}_S,\omega )=-\frac{2}{i\omega \rho _0}\partial _{3,S} G({\bf x},{\bf x}_S,\omega ), \end{eqnarray*}$$

(15)

where ∂_{3, S} denotes differentiation with respect to the source coordinate x_{3, S}. For x at |${{\partial \mathbb {D}}_R}$| (i.e. just below the source level) we have for the downgoing part

$$\begin{eqnarray*} \Gamma ^+({\bf x},{\bf x}_S,\omega )|_{x_3=x_{3,R}}=\delta ({\bf x}_{\rm H}-{\bf x}_{{\rm H},S}). \end{eqnarray*}$$

(16)

We define the reflection response R(x_R, x_S, ω) of the medium below |${{\partial \mathbb {D}}_R}$| as the upgoing part of the dipole-source response Γ(x_R, x_S, ω), with x_R at |${{\partial \mathbb {D}}_R}$|⁠, hence

$$\begin{eqnarray*} R({\bf x}_R,{\bf x}_S,\omega )&=&\Gamma ^-({\bf x}_R,{\bf x}_S,\omega )\nonumber \\ &=&-\frac{2}{i\omega \rho _0}\partial _{3,S} G^-({\bf x}_R,{\bf x}_S,\omega )\nonumber \\ &=&-\frac{2}{i\omega \rho _0}\partial _{3,S} G^{\rm s}({\bf x}_R,{\bf x}_S,\omega ), \end{eqnarray*}$$

(17)

where superscript s stands for scattered. Substituting p(x, ω) = Γ(x, x_S, ω) and p^±(x_R, ω) = Γ^±(x_R, x_S, ω) into eq. (11), using eqs (16) and (17), gives

$$\begin{eqnarray*} \Gamma ({\bf x},{\bf x}_S,\omega )=\int _{{{\partial \mathbb {D}}_R}} F({\bf x},{\bf x}_R,\omega )R({\bf x}_R,{\bf x}_S,\omega ){\rm d}{\bf x}_R +F^*({\bf x},{\bf x}_S,\omega ), \quad \mbox{for}\quad x_3 \ge x_{3,R}. \end{eqnarray*}$$

(18)

This is a representation for the dipole response Γ(x, x_S, ω) at virtual receiver position x anywhere in the half-space below |${{\partial \mathbb {D}}_R}$|⁠, expressed in terms of the reflection response R(x_R, x_S, ω) at |${{\partial \mathbb {D}}_R}$|⁠. It is similar to our earlier derived Green’s function representations for the Marchenko method, but here it has been derived without applying decomposition in the lower half-space. It only excludes the contribution from waves that are evanescent at |${{\partial \mathbb {D}}_R}$|⁠. Another difference with our earlier representations is that the Green’s function on the left-hand side is a dipole response instead of a monopole response. We address this in Section 3.2.

It may be counter intuitive that in eq. (18) we use focusing functions F(x, x_R, ω) and F*(x, x_S, ω), with their focal points x_R and x_S situated at the surface |${{\partial \mathbb {D}}_R}$|⁠. This is different from the focusing functions in the classical representations for the Marchenko method, which have their focal points in the subsurface. For an intuitive explanation of the right-hand side of eq. (18) we refer again to Fig. 2, this time with p⁺(x_R, ω) and p⁻(x_R, ω) replaced by δ(x_{H, R} − x_{H, S}) and R(x_R, x_S, ω), respectively. In a similar way as in Section 2, the focusing functions propagate the downgoing source field at x_S and the upgoing reflection response at x_R from |${{\partial \mathbb {D}}_R}$| to x in the lower half-space, with the focal points acting again as secondary sources in a modified form of Huygens’ principle. The sum of the two integrals gives Γ(x, x_S, ω).

3.2 Representation for the acoustic monopole Green’s function

In this section we turn eq. (18) into a representation for the monopole Green’s function G(x, x_S, ω). To this end we introduce a modified focusing function |$f({\bf x},{\bf x}_R,\omega )$| via

$$\begin{eqnarray*} F({\bf x},{\bf x}_R,\omega )=\frac{2}{i\omega \rho _0}\partial _{3,R} f({\bf x},{\bf x}_R,\omega ), \end{eqnarray*}$$

(19)

where ∂_{3, R} denotes differentiation with respect to x_{3, R}. According to eqs (8), (9) and (19), |$f({\bf x},{\bf x}_R,\omega )$| obeys the wave equation

$$\begin{eqnarray*} {\cal L}f=0, \end{eqnarray*}$$

(20)

the focusing condition

$$\begin{eqnarray*} \partial _{3,R} f({\bf x},{\bf x}_R,\omega )|_{x_3=x_{3,R}}=\frac{i\omega \rho _0}{2}\delta ({\bf x}_{\rm H}-{\bf x}_{{\rm H},R}), \end{eqnarray*}$$

(21)

and it is upgoing at and above |${{\partial \mathbb {D}}_R}$|⁠. Eq. (19) implies

$$\begin{eqnarray*} F^*({\bf x},{\bf x}_S,\omega )=-\frac{2}{i\omega \rho _0}\partial _{3,S} f^*({\bf x},{\bf x}_S,\omega ). \end{eqnarray*}$$

(22)

Substituting eqs (15), (17), (19) and (22) into eq. (18), applying source–receiver reciprocity to the scattered Green’s function and dropping the operation |$-\frac{2}{i\omega \rho _0}\partial _{3,S}$| from all terms gives

$$\begin{eqnarray*} G({\bf x},{\bf x}_S,\omega )&=& \frac{2}{i\omega \rho _0}\int _{{{\partial \mathbb {D}}_R}} \lbrace \partial _{3,R} f({\bf x},{\bf x}_R,\omega )\rbrace G^{\rm s}({\bf x}_S,{\bf x}_R,\omega ){\rm d}{\bf x}_R+f^*({\bf x},{\bf x}_S,\omega ),\quad\mbox{for}\quad x_3 \ge x_{3,R}. \end{eqnarray*}$$

(23)

We transfer the operator ∂_{3, R} from f to G^s, which is accompanied with a sign change (see Appendix A2). Using the definition of R in eq. (17) (with ∂_{3, S} replaced by ∂_{3, R} and with x_R and x_S interchanged on both sides of the equation) this yields

$$\begin{eqnarray*} G({\bf x},{\bf x}_S,\omega )=\int _{{{\partial \mathbb {D}}_R}} f({\bf x},{\bf x}_R,\omega )R({\bf x}_S,{\bf x}_R,\omega ){\rm d}{\bf x}_R +f^*({\bf x},{\bf x}_S,\omega ), \quad \mbox{for}\quad x_3 \ge x_{3,R}. \end{eqnarray*}$$

(24)

This is the main result of this paper. We discuss a number of aspects of this representation.

Eq. (24) has the same form as eq. (13) in Wapenaar et al. (2014), with f₂ in that paper replaced by f. Using |$\partial_{3,R} f({\bf x},{\bf x}_R,\omega )=-\partial_3 f({\bf x},{\bf x}_R,\omega )$| for x₃ = x_{3, R} (i.e. at the boundary of the homogeneous upper half-space), eq. (21) can be written as
$$\begin{eqnarray*} \partial _3 f({\bf x},{\bf x}_R,\omega )|_{x_3=x_{3,R}}=-\frac{i\omega \rho _0}{2}\delta ({\bf x}_{\rm H}-{\bf x}_{{\rm H},R}). \end{eqnarray*}$$
(25)
This is the same focusing condition as was defined for f₂(x, x_R, ω). An important difference between f and f₂ is the medium in which these focusing functions are defined. Focusing function f₂ is defined in a truncated version of the actual medium, where the medium below some depth level is replaced by a homogeneous medium. It is assumed that up/down decomposition is possible at the truncation level. On the other hand, focusing function f in eq. (24) is defined in the actual medium (similar as F in Fig. 1).
Moreover, the derivation in Wapenaar et al. (2014) of the representation is different: in that paper we start with decomposed focusing functions |$f_1^+({\bf x},{\bf x}_A,\omega )$| and |$f_1^-({\bf x},{\bf x}_A,\omega )$| in the truncated medium, with x_A being a focal point at the truncation depth. Next, we derive representations for decomposed Green’s functions G⁺(x_A, x_S, ω) and G⁻(x_A, x_S, ω) and combine the two into a single representation for G(x_A, x_S, ω) = G⁺(x_A, x_S, ω) + G⁻(x_A, x_S, ω), using the relation |$f_2({\bf x}_A,{\bf x}_R,\omega )=f_1^+({\bf x}_R,{\bf x}_A,\omega )-\lbrace f_1^-({\bf x}_R,{\bf x}_A,\omega )\rbrace ^*$| (note the different order of coordinates in f₁ and f₂). The latter relation is only valid when evanescent waves can be neglected at the truncation level inside the medium. In our current approach we do not make use of decomposition at a truncation level inside the medium and we avoid the approximate relation |$f_2=f_1^+-\lbrace f_1^-\rbrace ^*$|⁠. The only requirement for |$f({\bf x},{\bf x}_R,\omega )$| is that it obeys eqs (20) and (21). Hence, representation (24) gives the full wave field at the virtual receiver position x inside the medium, including multiply reflected, refracted and evanescent waves. It only excludes the contribution from waves that are evanescent at |${{\partial \mathbb {D}}_R}$|⁠, see the condition formulated by eq. (12).
Using another approach, also without applying decomposition inside the medium, Diekmann & Vasconcelos (2021) derive an equation of the same form as eq. (24), but without specifying a focusing condition for f. Their focusing function obeys a wave equation with a non-zero time-symmetric source function which is not explicitly specified. The derivation of eq. (24) in the current paper uses an explicit focusing condition (eq. 21) and does not require a source function for f.
Eq. (24) forms a starting point for deriving the Marchenko method. By applying an inverse Fourier transform we obtain
$$\begin{eqnarray*} G({\bf x},{\bf x}_S,t)-f({\bf x},{\bf x}_S,-t)= \int _{{{\partial \mathbb {D}}_R}}{\rm d}{\bf x}_R\int _{-\infty }^t f({\bf x},{\bf x}_R,t^{\prime })R({\bf x}_S,{\bf x}_R,t-t^{\prime }){\rm d} t^{\prime }, \quad \mbox{for}\quad x_3 \ge x_{3,R}. \end{eqnarray*}$$
(26)
The Marchenko method is based on the separability in time of G(x, x_S, t) and |$f({\bf x},{\bf x}_S,-t )$|⁠. For horizontal plane waves in 1-D media (Burridge 1980; Broggini & Snieder 2012) and for point-source responses at limited horizontal distances |x_H − x_{H, S}| in moderately inhomogeneous 3-D media (Wapenaar et al. 2013), these functions only overlap at t = t_d, which is the time of the direct arrival of the Green’s function. This minimum overlap in time allows the construction of a time-windowed version of eq. (26) with G(x, x_S, t) suppressed and with |$f({\bf x},{\bf x}_S,-t )$| almost completely preserved (this is the 3-D Marchenko equation). From this equation the focusing function |$f({\bf x},{\bf x}_S,t )$| can be resolved, given its direct arrival and the reflection response R(x_S, x_R, t). In essence the separability of the Green’s function and the time-reversed focusing function has been the underlying assumption of all implementations of the Marchenko method. This assumption excludes, among others, the treatment of refracted waves, which may arrive prior to the direct arrival of the Green’s function and interfere with the time-reversed focusing function.
Since we have argued that the representations of eqs (24) and (26) hold for refracted and evanescent waves, it is opportune to start new research on Marchenko methods which exploit the generality of these representations. Care should be taken to account for the overlap in time of the Green’s function and the time-reversed focusing function, particularly when dealing with refracted waves. A further discussion of the development of new Marchenko methods is beyond the scope of this paper.
Eq. (24) is, in principle, suited to retrieve the Green’s function G(x, x_S, ω) for x anywhere in the lower half-space. However, a single type of Green’s function is not a sufficient starting point for imaging. In the classical approach to Marchenko imaging, the downgoing and upgoing parts of the Green’s function are retrieved, from which a reflection image can be obtained, either by a deconvolution (Broggini et al. 2014; Wapenaar et al. 2014) or a correlation method (Behura et al. 2014). In the full-wavefield approach, we need at least one other type of field at x, next to G(x, x_S, ω), which represents the acoustic pressure field at x in response to a volume-injection rate source at x_S. To this end we introduce a Green’s function |$G_i^v({\bf x},{\bf x}_S,\omega )$| which, for i= 1, 2 and 3, stands for the three components of the particle velocity field at x. From the Fourier transform of eq. (1) we derive that the particle velocity |$v_i$| can be expressed in terms of the acoustic pressure as |$v_i=\frac{1}{i\omega }\vartheta _{ij}\partial _jp$|⁠. Similarly, we relate |$G_i^v$| to G via
$$\begin{eqnarray*} G_i^v({\bf x},{\bf x}_S,\omega )=\frac{1}{i\omega }\vartheta _{ij}({\bf x})\partial _jG({\bf x},{\bf x}_S,\omega ). \end{eqnarray*}$$
(27)
Hence, when G(x, x_S, ω) is available on a sufficiently dense grid, |$G_i^v({\bf x},{\bf x}_S,\omega )$| can be obtained via eq. (27). Alternatively, |$G_i^v({\bf x},{\bf x}_S,\omega )$| can be obtained from a modified version of the representation for G(x, x_S, ω). Applying the operation |$\frac{1}{i\omega }\vartheta _{ij}\partial _j$| to both sides of eq. (24) yields
$$\begin{eqnarray*} G_i^v({\bf x},{\bf x}_S,\omega )=\int _{{{\partial \mathbb {D}}_R}} h_i({\bf x},{\bf x}_R,\omega )R({\bf x}_S,{\bf x}_R,\omega ){\rm d}{\bf x}_R -h_i^*({\bf x},{\bf x}_S,\omega ), \quad \mbox{for}\quad x_3 \ge x_{3,R}, \end{eqnarray*}$$
(28)
with
$$\begin{eqnarray*} h_i({\bf x},{\bf x}_R,\omega )=\frac{1}{i\omega }\vartheta _{ij}({\bf x})\partial _jf({\bf x},{\bf x}_R,\omega ). \end{eqnarray*}$$
(29)
The Green’s functions G(x, x_S, ω) and |$G_i^v({\bf x},{\bf x}_S,\omega )$| together provide sufficient information for imaging. For example, one could decompose the field into incident and scattered waves in any desired direction, say in a direction perpendicular to a local interface (Yoon & Marfurt 2006; Liu et al. 2011; Holicki et al. 2019), and use these fields as input for imaging.

3.3 Representation for the homogeneous acoustic Green’s function

The representations in Sections 3.1 and 3.2 give the response to a source at x_S, observed by a virtual receiver at x inside the medium. Here we modify the representation of eq. (24), to create the response at the surface to a virtual source inside the medium. After that, we show how to obtain the response to this virtual source at a virtual receiver inside the medium.

We start by renaming the coordinate vectors in eq. (24) as follows: x_S → x_R, x_R → x_S, x → x_A. This yields, in combination with applying source–receiver reciprocity on the left-hand side of eq. (24),

$$\begin{eqnarray*} G({\bf x}_R,{{\bf x}_A},\omega ) =\int _{{{\partial \mathbb {D}}_R}} R({\bf x}_R,{\bf x}_S,\omega )f({{\bf x}_A},{\bf x}_S,\omega ){\rm d}{\bf x}_S +f^*({{\bf x}_A},{\bf x}_R,\omega ),\, \mbox{for}\quad x_{3,A} \ge x_{3,R}. \end{eqnarray*}$$

(30)

Here R(x_R, x_S, ω) is the reflection response to a dipole source at x_S, observed by a receiver at x_R, both at the surface |${{\partial \mathbb {D}}_R}$|⁠. This is schematically illustrated in Fig. 3(a). The integral in eq. (30) describes redatuming of the sources from all x_S at the surface to virtual-source position x_A in the subsurface, see Fig. 3(b). After adding f^*(x_A, x_R, ω) (according to eq. 30) this gives the Green’s function G(x_R, x_A, ω), which is the response to the virtual monopole source at x_A, observed by the receiver at x_R at the surface.

$Illustration of source and receiver redatuming as a two-step process. Starting with (a) the reflection response R(xR, xS, ω) at the surface, in step one (b) the Green’s function G(xR, xA, ω) is obtained for a virtual source at xA, and step two (c) yields the homogeneous Green’s function Gh(x, xA, ω) for a virtual receiver at x. All functions in this figure are represented by simple rays, but in reality these are wave fields, including primaries, multiples, refracted and evanescent waves.$

Figure 3.

Illustration of source and receiver redatuming as a two-step process. Starting with (a) the reflection response R(x_R, x_S, ω) at the surface, in step one (b) the Green’s function G(x_R, x_A, ω) is obtained for a virtual source at x_A, and step two (c) yields the homogeneous Green’s function G_h(x, x_A, ω) for a virtual receiver at x. All functions in this figure are represented by simple rays, but in reality these are wave fields, including primaries, multiples, refracted and evanescent waves.

Open in new tab Download slide

Our next aim is to derive a representation for the response observed by a virtual receiver at x in the subsurface, given G(x_R, x_A, ω). Eq. (11) cannot be used for this in the same way as before, since G(x, x_A, ω) obeys a wave equation with a singularity at x_A, whereas p(x, ω) in eq. (11) is not allowed to have sources in the lower half-space. To overcome this problem, we define the homogeneous Green’s function (Porter 1970; Oristaglio 1989)

$$\begin{eqnarray*} G_{\rm h}({\bf x},{{\bf x}_A},\omega )=G({\bf x},{{\bf x}_A},\omega )+G^*({\bf x},{{\bf x}_A},\omega ). \end{eqnarray*}$$

(31)

Here G(x, x_A, ω) and G*(x, x_A, ω) obey eq. (14), with source terms iωδ(x − x_A) and −iωδ(x − x_A), respectively, on the right-hand sides. Hence, G_h(x, x_A, ω) obeys the following equation:

$$\begin{eqnarray*} {\cal L} G_{\rm h}=0, \end{eqnarray*}$$

(32)

which confirms that the homogeneous Green’s function is source-free. This time we choose for p(x, ω) in eq. (11)

$$\begin{eqnarray*} p({\bf x},\omega )=G_{\rm h}({\bf x},{{\bf x}_A},\omega ), \end{eqnarray*}$$

(33)

with G_h(x, x_A, ω) defined in eq. (31). For x at |${{\partial \mathbb {D}}_R}$| the Green’s function G(x, x_A, ω) is purely upgoing, since the upper half-space is homogeneous and the virtual source at x_A lies in the lower half-space. Similarly, G*(x, x_A, ω) is downgoing at |${{\partial \mathbb {D}}_R}$|⁠, except for the evanescent field [which we already neglected at |${{\partial \mathbb {D}}_R}$| in the derivation of eq. (11)]. Hence, we may write

$$\begin{eqnarray*} p^-({\bf x},\omega )&=&G_{\rm h}^-({\bf x},{{\bf x}_A},\omega )= G({\bf x},{{\bf x}_A},\omega ),\quad \mbox{for}\quad x_3=x_{3,R}, \end{eqnarray*}$$

(34)

$$\begin{eqnarray*} p^+({\bf x},\omega )&=&G_{\rm h}^+({\bf x},{{\bf x}_A},\omega )= G^*({\bf x},{{\bf x}_A},\omega ), \quad \mbox{for}\quad x_3=x_{3,R}. \end{eqnarray*}$$

(35)

Substitution of eqs (33)–(35) into eq. (11) yields

$$\begin{eqnarray*} G_{\rm h}({\bf x},{{\bf x}_A},\omega )&=&\int _{{{\partial \mathbb {D}}_R}} F({\bf x},{\bf x}_R,\omega )G({\bf x}_R,{{\bf x}_A},\omega ){\rm d}{\bf x}_R +\int _{{{\partial \mathbb {D}}_R}} F^*({\bf x},{\bf x}_R,\omega )G^*({\bf x}_R,{{\bf x}_A},\omega ){\rm d}{\bf x}_R, \quad \mbox{for}\quad x_3\ge x_{3,R} , \end{eqnarray*}$$

(36)

or

$$\begin{eqnarray*} G_{\rm h}({\bf x},{{\bf x}_A},\omega )=2\Re \int _{{{\partial \mathbb {D}}_R}} F({\bf x},{\bf x}_R,\omega )G({\bf x}_R,{{\bf x}_A},\omega ){\rm d}{\bf x}_R, \quad \mbox{for}\quad x_3\ge x_{3,R}, \end{eqnarray*}$$

(37)

where ℜ denotes that the real part is taken. For an intuitive explanation of the right-hand side of eq. (36) we refer again to Fig. 2, this time with p⁺(x_R, ω) and p⁻(x_R, ω) replaced by G*(x_R, x_A, ω) and G(x_R, x_A, ω), respectively. The focusing functions propagate these downgoing and upgoing Green’s functions at x_R from |${{{\partial \mathbb {D}}_R}}$| into the lower half-space, with the focal points acting as secondary sources in a modified form of Huygens’ principle. The two integrals cannot be separately associated with G*(x, x_A, ω) and G(x, x_A, ω) for x in the lower half-space (these functions are singular at x_A); only the sum of the two integrals gives G_h(x, x_A, ω) (which is not singular at x_A).

Hence, eq. (36) describes redatuming of the receivers from all x_R at the surface to virtual-receiver position x in the subsurface, see Fig. 3(c). It gives the homogeneous Green’s function G_h(x, x_A, ω), which is the response to the virtual source at x_A, observed by a virtual receiver at x, plus its complex conjugate. After transforming this to the time domain we obtain

$$\begin{eqnarray*} G_{\rm h}({\bf x},{{\bf x}_A},t)=G({\bf x},{{\bf x}_A},t)+G({\bf x},{{\bf x}_A},-t). \end{eqnarray*}$$

(38)

The two functions at the right-hand side of this equation do not overlap in time (except for x = x_A and only for t = 0), hence, G(x, x_A, t) can be extracted from G_h(x, x_A, t) by selecting its causal part.

Note that there is an asymmetry in the focusing functions used for source redatuming [|$f({\bf x}_A,{\bf x}_S,\omega )$| in eq. (30)] and for receiver redatuming [F(x, x_R, ω) in eq. (37)], see also Fig. 3(c). This is due to the difference in types of responses at the surface [the dipole response R(x_R, x_S, ω)] and in the subsurface [the monopole response G(x, x_A, ω)]. When the response at the surface were also a monopole response, then the focusing function |$f({\bf x}_A,{\bf x}_S,\omega )$| for source redatuming should be replaced by F(x_A, x_S, ω).

Homogeneous Green’s function representations similar to eq. (37) were also derived by Wapenaar et al. (2016a), van der Neut et al. (2017) and Singh & Snieder (2017), but here eq. (37) has been derived without up/down decomposition inside the medium. Hence, it also holds for evanescent waves inside the medium, as long as condition (12) is obeyed. Moreover, the derivation presented here is much simpler than in those references.

The source and receiver redatuming processes can be captured in one equation by substituting eq. (30) into (37). This gives

$$\begin{eqnarray*} G_{\rm h}({\bf x},{{\bf x}_A},\omega )&=& 2\Re \int _{{{\partial \mathbb {D}}_R}}\int _{{{\partial \mathbb {D}}_R}} F({\bf x},{\bf x}_R,\omega )R({\bf x}_R,{\bf x}_S,\omega )f({{\bf x}_A},{\bf x}_S,\omega ){\rm d}{\bf x}_S{\rm d}{\bf x}_R \nonumber \\ &&+2\Re \int _{{{\partial \mathbb {D}}_R}} F({\bf x},{\bf x}_R,\omega )f^*({{\bf x}_A},{\bf x}_R,\omega ){\rm d}{\bf x}_R, \quad \mbox{for}\quad \lbrace x_3,x_{3,A}\rbrace \ge x_{3,R}. \end{eqnarray*}$$

(39)

The double integral on the right-hand side resembles the process of classical source and receiver redatuming (Berkhout 1982; Berryhill 1984), but with the primary focusing functions in those references replaced by full-field focusing functions. It also resembles source–receiver interferometry (Curtis & Halliday 2010), but with the double integration along a closed boundary in that paper replaced by the double integration over the open boundary |${{\partial \mathbb {D}}_R}$|⁠. Hence, via the theories of primary source–receiver redatuming (Berkhout 1982; Berryhill 1984), closed-boundary source–receiver interferometry (Curtis & Halliday 2010) and open-boundary homogeneous Green’s function retrieval using wave field decomposition (Wapenaar et al. 2016a; van der Neut et al. 2017; Singh & Snieder 2017), we have arrived at a representation for open-boundary homogeneous full-field Green’s function retrieval (eq. 39), which accounts for internal multiples, and refracted and evanescent waves in the lower half-space. In Section 5.3 this representation is extended for the elastodynamic situation.

3.4 Numerical examples

We illustrate the representations of Sections 3.2 and 3.3 with numerical examples. Our main aim is to demonstrate that these representations hold for evanescent waves inside the medium. To this end we consider oblique plane waves in a horizontally layered medium, with isotropic depth-dependent medium parameters c(x₃) (propagation velocity) and ρ(x₃) (mass density). We consider a horizontally layered medium because in this case we can unequivocally distinguish between propagating and evanescent waves. However, as discussed in Section 2, the representations also account for evanescent waves in more general inhomogeneous media. We define the spatial Fourier transform of a space- and frequency-dependent function u(x, ω) as

$$\begin{eqnarray*} \tilde{u}({\bf s},x_3,\omega )=\int _{\mathbb {R}^2}\exp \lbrace -i\omega {\bf s}\cdot {\bf x}_{\rm H}\rbrace u({\bf x}_{\rm H},x_3,\omega ){\rm d}{\bf x}_{\rm H}, \end{eqnarray*}$$

(40)

with s = (s₁, s₂), where s₁ and s₂ are horizontal slownesses and |$\mathbb {R}$| is the set of real numbers. This decomposes the function u(x, ω) into monochromatic plane-wave components. Next, we define the inverse temporal Fourier transform per slowness value as

$$\begin{eqnarray*} u({\bf s},x_3,\tau )=\frac{1}{\pi }\Re \int _0^\infty \tilde{u}({\bf s},x_3,\omega )\exp \lbrace -i\omega \tau ){\rm d}\omega , \end{eqnarray*}$$

(41)

where τ is the so-called intercept time (Stoffa 1989).

First we investigate the representation of eq. (24) and take x_S = (0, 0, x_{3, R}). We use the definitions of eqs (40) and (41) to transform this representation to the slowness intercept-time domain. Taking into account that for a horizontally layered, isotropic medium all functions in eq. (24) are cylindrically symmetric, it suffices to consider the transformed representation for one slowness variable only. We thus obtain

$$\begin{eqnarray*} G(s_1,x_3,x_{3,R},\tau ) &=&\int _{-\infty }^\tau f(s_1,x_3,x_{3,R},\tau ^{\prime })R(s_1,x_{3,R},\tau -\tau ^{\prime }){\rm d}\tau ^{\prime }+f(s_1,x_3,x_{3,R},-\tau ),\quad \mbox{for}\quad x_3\ge x_{3,R}. \end{eqnarray*}$$

(42)

For any given value of s₁, the Green’s function G(s₁, x₃, x_{3, R}, τ) is the response to a plane-wave source at x_{3, R} as a function of receiver depth x₃ and intercept time τ. For |s₁| ≤ 1/c(x₃) the plane wave is propagating, whereas for |s₁| > 1/c(x₃) it is evanescent. For propagating waves, the local propagation angle α(x₃) follows from s₁ = sin α(x₃)/c(x₃). The focusing function |$f(s_1,x_3,x_{3,R},\tau )$| obeys the focusing condition formulated by eq. (25), transformed to the slowness intercept-time domain, hence

$$\begin{eqnarray*} f(s_1,x_3,x_{3,R},\tau )|_{x_3=x_{3,R}}=\frac{\rho _0c_0}{2\cos \alpha _0}\delta (\tau ), \end{eqnarray*}$$

(43)

with α₀ = α(x_{3, R}). Consider the horizontally layered medium of Fig. 4(a). Two thin high-velocity layers (c₂ = c₄ = 3000 m s^–1) are embedded in a homogeneous background medium with a velocity of 2000 m s^–1. The mass densities, in kg m⁻³, are assigned the same numerical values as the velocities to get significant contrasts between the different layers. A plane wave is emitted from x_{3, R} into the medium, with slowness s₁ = 1/2800 s m⁻¹, hence, this wave leaves the surface with an angle α₀ = 45.6° and becomes evanescent in the high-velocity layers. For the source function we use a Ricker wavelet with a central frequency of 50 Hz, hence, the wavelength for the central frequency in the high-velocity layers is 60 m. The thickness of the high-velocity layers is 20 m, which is of the same order as the distance over which the evanescent waves decay with a factor 1/e, which is equal to |$1/(\omega _c\sqrt{s_1^2-1/c_2^2})$|=24.8 m. Hence, we may expect that the waves tunnel through these layers. Fig. 4(b) shows the numerically modelled reflection response R(s₁, x_{3, R}, τ) for the chosen slowness. The first two events are composite reflections from the two high-velocity layers (including internal multiples of evanescent waves inside these layers) and the other events are multiple reflections between these layers. Fig. 4(c) shows the numerically modelled focusing function |$f(s_1,x_3,x_{3,R},\tau )$| as a function of x₃ and τ, convolved with the same Ricker wavelet for a clear display. Blue and red arrows indicate upgoing and downgoing waves, respectively, in the homogeneous background medium. The tunnelling of the waves through the high-velocity layers is clearly visible. A single upgoing wave reaches the surface x_{3, R} at τ = 0, conform the focusing condition formulated by eq. (43) [except that in this display δ(τ) is convolved with the Ricker wavelet]. Note that the amplitude increases with increasing depth (to compensate for the evanescent waves in the high-velocity layers), which means that, in practice, the numerically computed focusing function becomes unstable beyond some thickness of the high-velocity layers.

$(a) Horizontally layered medium with two high-velocity layers. (b) Numerically modelled reflection response R(s1, x3, R, τ) at the surface. The horizontal slowness s1 = 1/2800 s m−1 is chosen such that the wave field is evanescent in the high-velocity layers. (c) Numerically modelled focusing function $f(s_1,x_3,x_{3,R},\tau )$. The trace at x3, R = 0 m illustrates the focusing condition of eq. (43).$

Figure 4.

(a) Horizontally layered medium with two high-velocity layers. (b) Numerically modelled reflection response R(s₁, x_{3, R}, τ) at the surface. The horizontal slowness s₁ = 1/2800 s m⁻¹ is chosen such that the wave field is evanescent in the high-velocity layers. (c) Numerically modelled focusing function |$f(s_1,x_3,x_{3,R},\tau )$|⁠. The trace at x_{3, R} = 0 m illustrates the focusing condition of eq. (43).

Open in new tab Download slide

The reflection response of Fig. 4(b) and the focusing function of Fig. 4(c) (the latter without the wavelet) are used as input for the representation of eq. (42). This yields the Green’s function G(s₁, x₃, x_{3, R}, τ) (convolved with the Ricker wavelet) as a function of x₃ and τ, see Fig. 5(a). Blue and red arrows indicate again upgoing and downgoing waves, respectively. This figure shows the expected behaviour of the response to a plane-wave source at x_{3, R} (a downgoing wave leaving the surface, two composite primary upgoing waves and multiple reflections between the high velocity layers). Fig. 5(b) shows G(s₁, x_{3, A}, x_{3, R}, τ) for x_{3, A} = 300 m. The green line is the Green’s function obtained from eq. (42), the red line is the directly modelled Green’s function. Similarly, Fig. 5(c) shows G(s₁, x_{3, B}, x_{3, R}, τ) for x_{3, B} = 210 m, that is inside the first high velocity layer. In both cases the match is perfect, which confirms that the representation of eq. (42) correctly accounts for propagating and evanescent waves inside the medium.

(a) Green’s function G(s1, x3, x3, R, τ) obtained from Figs 4(b) and (c) via the representation of eq. (42). (b) G(s1, x3, A, x3, R, τ), taken from figure (a) for x3, A = 300 m (green), compared with directly modelled Green’s function (red). (c) Similarly, G(s1, x3, B, x3, R, τ), taken from figure (a) for x3, B = 210 m inside the first high velocity layer. (d) Homogeneous Green’s function Gh(s1, x3, x3, A, τ) obtained from Figs 4(c) and 5(b) via the representation of eq. (44).

Figure 5.

(a) Green’s function G(s₁, x₃, x_{3, R}, τ) obtained from Figs 4(b) and (c) via the representation of eq. (42). (b) G(s₁, x_{3, A}, x_{3, R}, τ), taken from figure (a) for x_{3, A} = 300 m (green), compared with directly modelled Green’s function (red). (c) Similarly, G(s₁, x_{3, B}, x_{3, R}, τ), taken from figure (a) for x_{3, B} = 210 m inside the first high velocity layer. (d) Homogeneous Green’s function G_h(s₁, x₃, x_{3, A}, τ) obtained from Figs 4(c) and 5(b) via the representation of eq. (44).

Open in new tab Download slide

Using source–receiver reciprocity we may interpret Fig. 5(b) as G(s₁, x_{3, R}, x_{3, A}, τ), which is the response at the surface x_{3, R} to a virtual plane-wave source at x_{3, A} = 300 m. Hence, G(s₁, x_{3, R}, x_{3, A}, τ) may be seen as the result of redatuming the source from the surface to x_{3, A}. We now discuss receiver redatuming. To this end, we transform the representation of eq. (36) to the slowness intercept-time domain, which yields

$$\begin{eqnarray*} G_{\rm h}(s_1,x_3,x_{3,A},\tau ) &=&\int _{-\infty }^\tau F(s_1,x_3,x_{3,R},\tau ^{\prime })G(s_1,x_{3,R},x_{3,A},\tau -\tau ^{\prime }){\rm d}\tau ^{\prime }\nonumber \\ &&+\int _\tau ^\infty F(s_1,x_3,x_{3,R},-\tau ^{\prime })G(s_1,x_{3,R},x_{3,A},\tau ^{\prime }-\tau ){\rm d}\tau ^{\prime },\quad \mbox{for}\quad x_3\ge x_{3,R}, \end{eqnarray*}$$

(44)

with, analogous to eq. (19),

$$\begin{eqnarray*} F(s_1,x_3,x_{3,R},\tau )=\frac{2\cos \alpha _0}{\rho _0c_0}f(s_1,x_3,x_{3,R},\tau ). \end{eqnarray*}$$

(45)

Note that the right-hand side of eq. (44) contains the Green’s function with the redatumed source at x_{3, A} and the receiver at x_{3, R} at the surface. This representation redatums the receiver from x_{3, R} to any depth x₃ in the subsurface. This yields the homogeneous Green’s function, which consists of G(s₁, x₃, x_{3, A}, τ) plus its time-reversal, see Fig. 5(d). The causal part (right of the green dashed line) is the retrieved Green’s function G(s₁, x₃, x_{3, A}, τ). Conform expectation, we see a virtual source at x_{3, A} emitting downgoing and upgoing plane waves, which reverberate in the wave guide between the two high-velocity layers, but which also emit some energy through tunnelling into the half-spaces above and below the high-velocity layers. This example illustrates the handling of propagating and evanescent waves inside the medium by the homogeneous Green’s function representation of eq. (44).

4 ELASTODYNAMIC WAVE FIELD REPRESENTATION

We derive the elastodynamic equivalent of the representation of eq. (11). We consider the same configuration as in Section 2, except that now the medium parameters of the lower half-space x₃ > x_{3, R} are the stiffness tensor c_ijkl(x) and the mass density tensor ρ_ik(x), with symmetries c_ijkl = c_jikl = c_ijlk = c_klij and ρ_ik = ρ_ki. In the homogeneous isotropic upper half-space x₃ ≤ x_{3, R} the parameters are ρ_ik = δ_ikρ₀ and c_ijkl = λ₀δ_ijδ_kl + μ₀(δ_ikδ_jl + δ_ilδ_jk), with λ₀ and μ₀ the Lamé parameters of the half-space. The P- and S-wave propagation velocities of the upper half-space are c_P = [(λ₀ + 2μ₀)/ρ₀]^1/2 and c_S = (μ₀/ρ₀)^1/2, respectively.

The basic equations in the frequency domain for elastodynamic wave propagation are the linearized equation of motion

$$\begin{eqnarray*} -i\omega \rho _{ik}v_k-\partial _j\tau _{ij}=\,\hat{\!\!f}_i \end{eqnarray*}$$

(46)

and the linearized deformation equation

$$\begin{eqnarray*} i\omega \tau _{ij}+c_{ijkl}\partial _lv_k=0, \end{eqnarray*}$$

(47)

respectively. Here τ_ij(x, ω) is the stress tensor (with symmetry τ_ij = τ_ji), |$v_k({\bf x},\omega)$| the particle velocity and |$\,\hat{\!\!f}_i({\bf x},\omega )$| a source in terms of volume-force density (the circumflex is used to distinguish this source term from the focusing function). Eqs (46) and (47) can be combined into the elastodynamic wave equation

$$\begin{eqnarray*} {\cal L}_{ik}v_k=i\omega \,\,\hat{\!\!f}_i, \end{eqnarray*}$$

(48)

with

$$\begin{eqnarray*} {\cal L}_{ik}=\partial _jc_{ijkl}\partial _l+\omega ^2\rho _{ik}. \end{eqnarray*}$$

(49)

We introduce an elastodynamic focusing function F(x, x_R, ω) as a 3 × 3 matrix, according to

$$\begin{eqnarray*} {\bf F}({\bf x},{\bf x}_R,\omega )= {\begin{pmatrix}F_{1,1} & F_{1,2} & F_{1,3} \\ F_{2,1} & F_{2,2} & F_{2,3} \\ F_{3,1} & F_{3,2} & F_{3,3} \end{pmatrix}}({\bf x},{\bf x}_R,\omega ), \end{eqnarray*}$$

(50)

where x_R denotes again the position of a focal point at |${{\partial \mathbb {D}}_R}$|⁠. Each column of F is a particle velocity vector of which the components, for fixed x_R and variable x, obey the elastodynamic wave eq. (48) for the source-free situation. This is different from the elastodynamic focusing function introduced by Wapenaar & Slob (2014), in which the different elements represent decomposed compressional and shear waves.

We define the focusing condition, analogous to eq. (9), as

$$\begin{eqnarray*} {\bf F}({\bf x},{\bf x}_R,\omega )|_{x_3=x_{3,R}}={\bf I}\delta ({\bf x}_{\rm H}-{\bf x}_{{\rm H},R}), \end{eqnarray*}$$

(51)

(I is the 3 × 3 identity matrix) and demand that F(x, x_R, ω) is purely upgoing at |${{\partial \mathbb {D}}_R}$| and in the homogeneous isotropic upper half-space. Eq. (51) implies that, for the kth column of F, the kth component of the particle velocity vector in that column focuses at x_R and the other two components are zero on |${{\partial \mathbb {D}}_R}$|⁠. Hence, the columns of F are mutually independent.

We discuss a representation for a wave field |$v_k({\bf x},\omega)$|⁠, which may have sources in the upper half-space above |${{\partial \mathbb {D}}_R}$|⁠, but which obeys the source-free wave equation |${\cal L}_{ik}v_k=0$| for x₃ ≥ x_{3, R}. We store the components |$v_k({\bf x},\omega)$| in a 3 × 1 vector v(x, ω) as follows:

$$\begin{eqnarray*} {\bf v}({\bf x},\omega )={\begin{pmatrix}v_1\\ v_2\\ v_3\end{pmatrix}}({\bf x},\omega ). \end{eqnarray*}$$

(52)

In the lower half-space we express v(x, ω) as a superposition of mutually independent wave fields that obey the same source-free wave equation as |${\bf v}({\bf x},\omega)$| for x₃ ≥ x_{3, R}. For this purpose we choose the focusing functions F(x, x_R, ω) and F*(x, x_R,ω), of which the columns are also mutually independent. Hence, analogous to eq. (10) we express v(x, ω) as

$$\begin{eqnarray*} {\bf v}({\bf x},\omega )&=&\int _{{{\partial \mathbb {D}}_R}} {\bf F}({\bf x},{\bf x}_R,\omega ){\bf a}({\bf x}_R,\omega ){\rm d}{\bf x}_R +\int _{{{\partial \mathbb {D}}_R}} {\bf F}^*({\bf x},{\bf x}_R,\omega ){\bf b}({\bf x}_R,\omega ){\rm d}{\bf x}_R, \quad\mbox{for}\quad x_3 \ge x_{3,R}. \end{eqnarray*}$$

(53)

Here a(x_R, ω) and b(x_R, ω) are as yet undetermined 3 × 1 vectors. In Appendix B1, we formulate boundary conditions for the particle velocity and traction vectors at |${{\partial \mathbb {D}}_R}$|⁠, from which we solve a(x_R, ω) and b(x_R, ω). We thus obtain

$$\begin{eqnarray*} {\bf v}({\bf x},\omega )&=&\int _{{{\partial \mathbb {D}}_R}} {\bf F}({\bf x},{\bf x}_R,\omega ){\bf v}^-({\bf x}_R,\omega ){\rm d}{\bf x}_R +\int _{{{\partial \mathbb {D}}_R}} {\bf F}^*({\bf x},{\bf x}_R,\omega ){\bf v}^+({\bf x}_R,\omega ){\rm d}{\bf x}_R, \quad\mbox{for}\quad x_3 \ge x_{3,R}, \end{eqnarray*}$$

(54)

where v⁻(x_R, ω) and v⁺(x_R, ω) represent the upgoing and downgoing parts, respectively, of v(x_R, ω) for x_R at |${{\partial \mathbb {D}}_R}$|⁠. These upgoing and downgoing fields are velocity-normalized, meaning that v⁻ + v⁺ = v at and above |${{\partial \mathbb {D}}_R}$|⁠. Below |${{\partial \mathbb {D}}_R}$| we only consider the total (undecomposed) wave field v.

The explanation of the right-hand side of eq. (54) in terms of Huygens’ principle is similar to that of eq. (11). The main extension is that the matrix–vector products in eq. (54) accomplish a summation over the different components of the secondary sources at |${{{\partial \mathbb {D}}_R}}$| (corresponding to the foci of the different columns of the focusing function F).

As for the acoustic representation of eq. (11), the underlying assumption in the derivation of eq. (54) is that evanescent waves can be neglected at |${{\partial \mathbb {D}}_R}$|⁠. Hence, it only holds for waves which have a horizontal slowness s which obeys

$$\begin{eqnarray*} |{\bf s}|\le 1/c_P, \quad \mbox{at}\, \, {{\partial \mathbb {D}}_R}. \end{eqnarray*}$$

(55)

Using similar arguments as given below eq. (12), it follows that eq. (54) accounts for evanescent waves inside the medium, as long as they are related to propagating waves at the surface, as formulated by eq. (55).

5 ELASTODYNAMIC GREEN’S FUNCTION REPRESENTATIONS

5.1 Representation for a modified elastodynamic Green’s function

We introduce the elastodynamic Green’s function G_{k, n}(x, x_S, ω) as a solution of eq. (48) for a unit point source of volume-force density at x_S in the x_n-direction, hence

$$\begin{eqnarray*} {\cal L}_{ik}G_{k,n}=i\omega \delta _{in}\delta ({\bf x}-{\bf x}_S). \end{eqnarray*}$$

(56)

We demand that the time domain version of G_{k, n}(x, x_S, ω) is causal. Note that G_{k, n} obeys source–receiver reciprocity, that is G_{k, n}(x, x_S, ω) = G_{n, k}(x_S, x, ω). We introduce G(x, x_S, ω) as a 3 × 3 matrix, according to

$$\begin{eqnarray*} {\bf G}({\bf x},{\bf x}_S,\omega )= {\begin{pmatrix}G_{1,1} & G_{1,2} & G_{1,3} \\ G_{2,1} & G_{2,2} & G_{2,3} \\ G_{3,1} & G_{3,2} & G_{3,3} \end{pmatrix}}({\bf x},{\bf x}_S,\omega ). \end{eqnarray*}$$

(57)

Each column is a particle velocity vector of which the components, for fixed x_S and variable x, obey wave eq. (56). The different columns correspond to different directions of the force source at x_S. This is different from the elastodynamic Green’s function used by Wapenaar & Slob (2014), in which the different elements represent decomposed compressional and shear waves. In matrix form, source–receiver reciprocity implies G(x, x_S, ω) = {G(x_S, x, ω)}^t, where superscript t denotes transposition.

We choose x_S = (x_{H, S}, x_{3, S}) again in the upper half-space, at a vanishing distance ϵ above |${{\partial \mathbb {D}}_R}$|⁠, hence, x_{3, S} = x_{3, R} − ϵ. In Appendix B2, we derive a modified version Γ(x, x_S, ω) of G(x, x_S, ω) (eq. B20), of which the downgoing part Γ⁺(x, x_S, ω) for x at |${{\partial \mathbb {D}}_R}$| (i.e. just below the source level) is equal to a spatial delta function. Hence

$$\begin{eqnarray*} {\rm{\boldsymbol \Gamma}}^+({\bf x},{\bf x}_S,\omega )|_{x_3=x_{3,R}}={\bf I}\delta ({\bf x}_{\rm H}-{\bf x}_{{\rm H},S}). \end{eqnarray*}$$

(58)

We define the reflection response R(x_R, x_S, ω) of the medium below |${{\partial \mathbb {D}}_R}$| as the upgoing part of Γ(x_R, x_S, ω), with x_R at |${{\partial \mathbb {D}}_R}$|⁠, hence

$$\begin{eqnarray*} {\bf R}({\bf x}_R,{\bf x}_S,\omega )={\rm{\boldsymbol \Gamma}}^-({\bf x}_R,{\bf x}_S,\omega ). \end{eqnarray*}$$

(59)

Substituting v(x, ω) = Γ(x, x_S, ω) and v^±(x_R, ω) = Γ^±(x_R, x_S, ω) into eq. (54), using eqs (58) and (59), gives

$$\begin{eqnarray*} {\rm{\boldsymbol \Gamma}}({\bf x},{\bf x}_S,\omega )=\int _{{{\partial \mathbb {D}}_R}} {\bf F}({\bf x},{\bf x}_R,\omega ){\bf R}({\bf x}_R,{\bf x}_S,\omega ){\rm d}{\bf x}_R +{\bf F}^*({\bf x},{\bf x}_S,\omega ),\quad \mbox{for}\quad x_3 \ge x_{3,R}. \end{eqnarray*}$$

(60)

This is a representation for the modified version Γ(x, x_S, ω) of the elastodynamic Green’s function. It has been derived without applying decomposition in the lower half-space. It only excludes the contribution from waves that are evanescent at |${{\partial \mathbb {D}}_R}$|⁠.

5.2 Representation for the elastodynamic Green’s function

In Appendix B3 we show that eq. (60) can be reorganized into the following representation for the elastodynamic Green’s function G(x, x_S, ω)

$$\begin{eqnarray*} {\bf G}({\bf x},{\bf x}_S,\omega )=\int _{{{\partial \mathbb {D}}_R}} {\bf f}({\bf x},{\bf x}_R,\omega )\lbrace {\bf R}({\bf x}_S,{\bf x}_R,\omega )\rbrace ^t{\rm d}{\bf x}_R +{\bf f}^*({\bf x},{\bf x}_S,\omega ), \quad \mbox{for}\quad x_3 \ge x_{3,R}. \end{eqnarray*}$$

(61)

Here f(x, x_R, ω) is a modified version of the focusing function F(x, x_R, ω) (eq. B25). This representation gives the full elastodynamic particle velocity field at any virtual receiver position x inside the medium. It is similar to earlier derived elastodynamic representations for the Marchenko method (da Costa Filho et al. 2014; Wapenaar & Slob 2014), but here it has been derived without applying decomposition at a truncation level inside the medium. As a consequence, eq. (61) gives the full wave field at any virtual receiver position x inside the medium, including multiply reflected, converted, refracted and evanescent waves. This representation only excludes the contribution from waves that are evanescent at |${{\partial \mathbb {D}}_R}$|⁠, see the condition formulated by eq. (55).

Applying elastodynamic representations like the one in eq. (61) to derive a Marchenko method is not trivial. The functions G(x, x_S, ω) and f*(x, x_S, ω), transformed back to the time domain, partly overlap and hence they cannot be completely separated by a time window [similar as discussed by Wapenaar & Slob (2014) and Reinicke et al. (2020) for Green’s functions and focusing functions consisting of decomposed compressional and shear waves]. A discussion of elastodynamic Marchenko methods is beyond the scope of this paper.

Similar as in the acoustic situation, the representation of eq. (61) is not a sufficient starting point for imaging. We need at least one other type of field at x, next to G(x, x_S, ω), which represents the particle velocity at x in response to force sources at x_S. To this end, we introduce a Green’s function |${\bf G}_j^{\tau }({\bf x},{\bf x}_S,\omega )$| which, for j = 1, 2, 3, stands for the three traction vectors at x. From eq. (47) we derive that the traction vector |${\rm{\boldsymbol \tau}}_j$| can be expressed in terms of the particle velocity as |${\rm{\boldsymbol \tau}}_j=-\frac{1}{i\omega }{\bf C}_{jl}\partial _l{\bf v}$|⁠, with |$({\rm{\boldsymbol \tau}}_j)_i=\tau_{ij}$| and (C_jl)_ik = c_ijkl. Similarly, we relate |${\bf G}_j^{\tau }$| to G via

$$\begin{eqnarray*} {\bf G}_j^{\tau }({\bf x},{\bf x}_S,\omega )=-\frac{1}{i\omega }{\bf C}_{jl}({\bf x})\partial _l{\bf G}({\bf x},{\bf x}_S,\omega ). \end{eqnarray*}$$

(62)

Hence, when G(x, x_S, ω) is available on a sufficiently dense grid, |${\bf G}_j^{\tau }({\bf x},{\bf x}_S,\omega )$| can be obtained via eq. (62). Alternatively, |${\bf G}_j^{\tau }({\bf x},{\bf x}_S,\omega )$| can be obtained from a modified version of the representation for G(x, x_S, ω). Applying the operation |$-\frac{1}{i\omega }{\bf C}_{jl}\partial _l$| to both sides of eq. (61) yields

$$\begin{eqnarray*} {\bf G}_j^{\tau }({\bf x},{\bf x}_S,\omega )=\int _{{{\partial \mathbb {D}}_R}} {\bf h}_j({\bf x},{\bf x}_R,\omega )\lbrace {\bf R}({\bf x}_S,{\bf x}_R,\omega )\rbrace ^t{\rm d}{\bf x}_R -{\bf h}_j^*({\bf x},{\bf x}_S,\omega ), \quad \mbox{for}\quad x_3 \ge x_{3,R}, \end{eqnarray*}$$

(63)

with

$$\begin{eqnarray*} {\bf h}_j({\bf x},{\bf x}_R,\omega )=-\frac{1}{i\omega }{\bf C}_{jl}({\bf x})\partial _l{\bf f}({\bf x},{\bf x}_R,\omega ). \end{eqnarray*}$$

(64)

The Green’s functions G(x, x_S, ω) and |${\bf G}_j^{\tau }({\bf x},{\bf x}_S,\omega )$| together provide sufficient information for imaging.

5.3 Representation for the homogeneous elastodynamic Green’s function

The representations in Sections 5.1 and 5.2 give the elastodynamic response to a source at x_S, observed by a virtual receiver at x inside the medium. Similar as in Section 3.3, here we modify the representation of eq. (61), to create the response at the surface to a virtual source inside the medium. After that, we show how to obtain the response to this virtual source at a virtual receiver inside the medium.

We start by renaming the coordinate vectors in eq. (61) as follows: x_S → x_R, x_R → x_S, x → x_A. This yields, in combination with transposing all terms and applying source–receiver reciprocity on the left-hand side of eq. (61),

$$\begin{eqnarray*} {\bf G}({\bf x}_R,{{\bf x}_A},\omega )=\int _{{{\partial \mathbb {D}}_R}} {\bf R}({\bf x}_R,{\bf x}_S,\omega ) {\bf f}^t({{\bf x}_A},{\bf x}_S,\omega ){\rm d}{\bf x}_S +{\bf f}^\dagger ({{\bf x}_A},{\bf x}_R,\omega ), \quad \mbox{for}\quad x_{3,A} \ge x_{3,R}. \end{eqnarray*}$$

(65)

Here superscript † denotes transposition and complex conjugation. The integral in eq. (65) describes elastodynamic redatuming of the sources from all x_S at the surface to virtual-source position x_A in the subsurface.

Our next aim is to derive a representation for the response observed by a virtual receiver at x in the subsurface, given G(x_R, x_A, ω). Similar as in Section 3.3, we define the homogeneous elastodynamic Green’s function

$$\begin{eqnarray*} {\bf G}_{\rm h}({\bf x},{{\bf x}_A},\omega )={\bf G}({\bf x},{{\bf x}_A},\omega )+{\bf G}^*({\bf x},{{\bf x}_A},\omega ). \end{eqnarray*}$$

(66)

The components of G(x, x_A, ω) and G*(x, x_A,ω) obey eq. (56), with source terms iωδ_inδ(x − x_A) and −iωδ_inδ(x − x_A), respectively, on the right-hand sides. Hence, the components of G_h(x, x_A, ω) obey this equation without a source on the right-hand side. Following a similar reasoning as in Section 3.3, we substitute

$$\begin{eqnarray*} {\bf v}({\bf x},\omega )&=&{\bf G}_{\rm h}({\bf x},{{\bf x}_A},\omega ), \end{eqnarray*}$$

(67)

$$\begin{eqnarray*} {\bf v}^-({\bf x},\omega )&=&{\bf G}_{\rm h}^-({\bf x},{{\bf x}_A},\omega )= {\bf G}({\bf x},{{\bf x}_A},\omega ),\quad \mbox{for}\quad x_3=x_{3,R}, \end{eqnarray*}$$

(68)

$$\begin{eqnarray*} {\bf v}^+({\bf x},\omega )&=&{\bf G}_{\rm h}^+({\bf x},{{\bf x}_A},\omega )= {\bf G}^*({\bf x},{{\bf x}_A},\omega ),\quad \mbox{for}\quad x_3=x_{3,R}, \end{eqnarray*}$$

(69)

into eq. (54). This gives

$$\begin{eqnarray*} {\bf G}_{\rm h}({\bf x},{{\bf x}_A},\omega )=2\Re \int _{{{\partial \mathbb {D}}_R}} {\bf F}({\bf x},{\bf x}_R,\omega ){\bf G}({\bf x}_R,{{\bf x}_A},\omega ){\rm d}{\bf x}_R, \quad \mbox{for}\quad x_3 \ge x_{3,R}. \end{eqnarray*}$$

(70)

This equation describes elastodynamic redatuming of the receivers from all x_R at the surface to virtual-receiver position x in the subsurface. It gives the homogeneous Green’s function G_h(x, x_A, ω), which is the response to the virtual source at x_A, observed by a virtual receiver at x, plus its complex conjugate. After transforming this to the time domain, G(x, x_A, t) can be extracted from G_h(x, x_A, t) by selecting its causal part.

An elastodynamic homogeneous Green’s function representation similar to eq. (70) was also derived by Wapenaar et al. (2016b) and illustrated with numerical examples by Reinicke & Wapenaar (2019), but here the fields are not decomposed into downgoing and upgoing compressional and shear waves inside the medium. Hence, it also holds for evanescent waves inside the medium, as long as condition (55) is obeyed. Moreover, the derivation presented here is much simpler than in those references.

The source and receiver redatuming processes can be captured in one equation by substituting eq. (65) into (70). This gives

$$\begin{eqnarray*} {\bf G}_{\rm h}({\bf x},{{\bf x}_A},\omega )&=& 2\Re \int _{{{\partial \mathbb {D}}_R}}\int _{{{\partial \mathbb {D}}_R}} {\bf F}({\bf x},{\bf x}_R,\omega ){\bf R}({\bf x}_R,{\bf x}_S,\omega ){\bf f}^t({{\bf x}_A},{\bf x}_S,\omega ){\rm d}{\bf x}_S{\rm d}{\bf x}_R \nonumber \\ &&+2\Re \int _{{{\partial \mathbb {D}}_R}} {\bf F}({\bf x},{\bf x}_R,\omega ){\bf f}^\dagger ({{\bf x}_A},{\bf x}_R,\omega ){\rm d}{\bf x}_R, \quad \mbox{for}\quad \lbrace x_3,x_{3,A}\rbrace \ge x_{3,R}. \end{eqnarray*}$$

(71)

The double integral on the right-hand side resembles the process of classical elastodynamic source and receiver redatuming (Wapenaar & Berkhout 1989; Hokstad 2000), but with the primary focusing functions in those references replaced by full-field focusing functions. It also resembles elastodynamic source–receiver interferometry (Halliday et al. 2012), but with the double integration along a closed boundary in that paper replaced by the double integration over the open boundary |${{\partial \mathbb {D}}_R}$|⁠. Hence, via the theories of elastodynamic primary source–receiver redatuming (Wapenaar & Berkhout 1989; Hokstad 2000), closed-boundary source–receiver interferometry (Halliday et al. 2012) and open-boundary homogeneous Green’s function retrieval using wave field decomposition (Wapenaar et al. 2016b), we have arrived at a representation for elastodynamic open-boundary homogeneous full-field Green’s function retrieval (eq. 71), which accounts for internal multiples, converted, refracted and evanescent waves in the lower half-space.

6 CONCLUSIONS

We have derived acoustic and elastodynamic Green’s function representations in terms of the reflection response at the surface and focusing functions. These representations have the same form as the representations that we derived earlier as the basis for Marchenko redatuming, imaging, monitoring and multiple elimination. However, unlike in our original derivations, we did not assume that the wave field inside the medium can be decomposed into downgoing and upgoing waves and we did not ignore evanescent waves inside the medium. We only neglected the contribution of waves that are evanescent at the acquisition boundary. We have demonstrated with numerical examples that the representations indeed account for evanescent waves inside the medium. The representations form a starting point for new research on Marchenko methods which circumvent the limitations caused by the assumptions underlying the traditional representations. In these new developments, care should be taken to account for the overlap in time of the Green’s function and the time-reversed focusing function, particularly when dealing with refracted waves.

ACKNOWLEDGEMENTS

We thank Marcin Dukalski, Mert Sinan Recep Kiraz and two anonymous reviewers for their comments, which helped us sharpen the message. This work has received funding from the European Union’s research and innovation programme: European Research Council (grant agreement 742703).

DATA AVAILABILITY

This study uses numerical data only. The codes used will be shared on reasonable request to the corresponding author.

REFERENCES

Behura

J.

,

Wapenaar

K.

,

Snieder

R.

,

2014

.

Autofocus imaging: Image reconstruction based on inverse scattering theory

,

Geophysics

,

79

(

3

),

A19

–

A26

.

10.1190/geo2013-0398.1

Berkhout

A.J.

,

1982

.

Seismic Migration. Imaging of Acoustic Energy by Wave Field Extrapolation. A. Theoretical Aspects

,

Elsevier

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Berryhill

J.R.

,

1984

.

Wave-equation datuming before stack

,

Geophysics

,

49

,

2064

–

2066

.

Brackenhoff

J.

,

Thorbecke

J.

,

Wapenaar

K.

,

2019

.

Monitoring of induced distributed double-couple sources using Marchenko-based virtual receivers

,

Solid Earth

,

10

,

1301

–

1319

.

10.5194/se-10-1301-2019

Broggini

F.

,

Snieder

R.

,

2012

.

Connection of scattering principles: a visual and mathematical tour

,

Eur. J. Phys.

,

33

,

593

–

613

.

10.1088/0143-0807/33/3/593

Broggini

F.

,

Snieder

R.

,

Wapenaar

K.

,

2014

.

Data-driven wavefield focusing and imaging with multidimensional deconvolution: numerical examples for reflection data with internal multiples

,

Geophysics

,

79

(

3

),

WA107

–

WA115

.

10.1190/geo2013-0307.1

Burridge

R.

,

1980

.

The Gelfand-Levitan, the Marchenko, and the Gopinath-Sondhi integral equations of inverse scattering theory, regarded in the context of inverse impulse-response problems

,

Wave Motion

,

2

,

305

–

323

.

10.1016/0165-2125(80)90011-6

Corones

J.P.

,

1975

.

Bremmer series that correct parabolic approximations

,

J. Math. Anal. Appl.

,

50

,

361

–

372

.

10.1016/0022-247X(75)90029-3

Curtis

A.

,

Halliday

D.

,

2010

.

Source-receiver wavefield interferometry

,

Phys. Rev. E

,

81

,

046601

.

10.1103/PhysRevE.81.046601

da Costa Filho

C.A.

,

Ravasi

M.

,

Curtis

A.

,

Meles

G.A.

,

2014

.

Elastodynamic Green’s function retrieval through single-sided Marchenko inverse scattering

,

Phys. Rev. E

,

90

,

063201

.

10.1103/PhysRevE.90.063201

Diekmann

L.

,

Vasconcelos

I.

,

2021

.

Focusing and Green’s function retrieval in three-dimensional inverse scattering revisited: a single-sided Marchenko integral for the full wave field

,

Phys. Rev. Research

,

3

,

013206

.

10.1103/PhysRevResearch.3.013206

Elison

P.

,

Dukalski

M.S.

,

de Vos

K.

,

van Manen

D.J.

,

Robertsson

J.O.A.

,

2020

.

Data-driven control over short-period internal multiples in media with a horizontally layered overburden

,

Geophys. J. Int.

,

221

,

769

–

787

.

Fishman

L.

,

McCoy

J.J.

,

1984

.

Derivation and application of extended parabolic wave theories. I. The factorized Helmholtz equation

,

J. Math. Phys.

,

25

(

2

),

285

–

296

.

Halliday

D.

,

Curtis

A.

,

Wapenaar

K.

,

2012

.

Generalized PP + PS = SS from seismic interferometry

,

Geophys. J. Int.

,

189

,

1015

–

1024

.

10.1111/j.1365-246X.2012.05396.x

Hokstad

K.

,

2000

.

Multicomponent Kirchhoff migration

,

Geophysics

,

65

(

3

),

861

–

873

.

Holicki

M.

,

Drijkoningen

G.

,

Wapenaar

K.

,

2019

.

Acoustic directional snapshot wavefield decomposition

,

Geophys. Prosp.

,

67

,

32

–

51

.

10.1111/1365-2478.12700

Jia

X.

,

Guitton

A.

,

Snieder

R.

,

2018

.

A practical implementation of subsalt Marchenko imaging with a Gulf of Mexico data set

,

Geophysics

,

83

(

5

),

S409

–

S419

.

10.1190/geo2017-0646.1

Kennett

B.L.N.

,

Kerry

N.J.

,

Woodhouse

J.H.

,

1978

.

Symmetries in the reflection and transmission of elastic waves

,

Geophys. J. R. astr. Soc.

,

52

,

215

–

230

.

10.1111/j.1365-246X.1978.tb04230.x

Kiraz

M.S.R.

,

Snieder

R.

,

Wapenaar

K.

,

2021

.

Focusing waves in an unknown medium without wavefield decomposition

,

JASA Express Lett.

,

1

(

5

),

055602

.

Liu

F.

,

Zhang

G.

,

Morton

S.A.

,

Leveille

J.P.

,

2011

.

An effective imaging condition for reverse-time migration using wavefield decomposition

,

Geophysics

,

76

,

S29

–

S39

.

Lomas

A.

,

Curtis

A.

,

2019

.

An introduction to Marchenko methods for imaging

,

Geophysics

,

84

(

2

),

F35

–

F45

.

10.1190/geo2018-0068.1

Mildner

C.

,

Broggini

F.

,

de Vos

K.

,

Robertsson

J.O.A.

,

2019

.

Accurate source wavelet estimation using Marchenko focusing functions

,

Geophysics

,

84

(

6

),

Q73

–

Q88

.

10.1190/geo2018-0726.1

Oristaglio

M.L.

,

1989

.

An inverse scattering formula that uses all the data

,

Inverse Probl.

,

5

,

1097

–

1105

.

10.1088/0266-5611/5/6/015

Pereira

R.

,

Ramzy

M.

,

Griscenco

P.

,

Huard

B.

,

Huang

H.

,

Cypriano

L.

,

Khalil

A.

,

2019

.

Internal multiple attenuation for OBN data with overburden/target separation

, in

Proceedings of the 89th Annual Meeting of the Society of Exploration Geophysicists

, pp.

4520

–

4524

.

Porter

R.P.

,

1970

.

Diffraction-limited, scalar image formation with holograms of arbitrary shape

,

J. Opt. Soc. Am.

,

60

,

1051

–

1059

.

10.1364/JOSA.60.001051

Ravasi

M.

,

Vasconcelos

I.

,

2021

.

An open-source framework for the implementation of large-scale integral operators with flexible, modern HPC solutions - enabling 3D Marchenko imaging by least squares inversion

,

Geophysics

,

86

(

early view)

,

doi:10.1190/geo2020–0796.1

.

10.1190/geo2020–0796.1

Google Scholar

OpenURL Placeholder Text

WorldCat

Crossref

Ravasi

M.

,

Vasconcelos

I.

,

Kritski

A.

,

Curtis

A.

,

da Costa Filho

C.A.

,

Meles

G.A.

,

2016

.

Target-oriented Marchenko imaging of a North Sea field

,

Geophys. J. Int.

,

205

,

99

–

104

.

Reinicke

C.

,

Wapenaar

K.

,

2019

.

Elastodynamic single-sided homogeneous Green’s function representation: theory and numerical examples

,

Wave Motion

,

89

,

245

–

264

.

10.1016/j.wavemoti.2019.04.001

Reinicke

C.

,

Dukalski

M.

,

Wapenaar

K.

,

2020

.

Comparison of monotonicity challenges encountered by the inverse scattering series and the Marchenko demultiple method for elastic waves

,

Geophysics

,

85

(

5

),

Q11

–

Q26

.

10.1190/geo2019-0674.1

Schoenberg

M.

,

Sen

P.N.

,

1983

.

Properties of a periodically stratified acoustic half-space and its relation to a Biot fluid

,

J. acoust. Soc. Am.

,

73

,

61

–

67

.

Singh

S.

,

Snieder

R.

,

2017

.

Source-receiver Marchenko redatuming: obtaining virtual receivers and virtual sources in the subsurface

,

Geophysics

,

82

(

3

),

Q13

–

Q21

.

10.1190/geo2016-0074.1

Slob

E.

,

Wapenaar

K.

,

Broggini

F.

,

Snieder

R.

,

2014

.

Seismic reflector imaging using internal multiples with Marchenko-type equations

,

Geophysics

,

79

(

2

),

S63

–

S76

.

10.1190/geo2013-0095.1

Staring

M.

,

Wapenaar

K.

,

2020

.

Three-dimensional Marchenko internal multiple attenuation on narrow azimuth streamer data of the Santos Basin, Brazil

,

Geophys. Prosp.

,

68

,

1864

–

1877

.

10.1111/1365-2478.12964

Staring

M.

,

Pereira

R.

,

Douma

H.

,

van der Neut

J.

,

Wapenaar

K.

,

2018

.

Source-receiver Marchenko redatuming on field data using an adaptive double-focusing method

,

Geophysics

,

83

(

6

),

S579

–

S590

.

10.1190/geo2017-0796.1

Stoffa

P.L.

,

1989

.

Tau-p - A Plane Wave Approach to the Analysis of Seismic Data

,

Kluwer Academic Publishers

.

10.1007/978-94-009-0881-9

Ursin

B.

,

1983

.

Review of elastic and electromagnetic wave propagation in horizontally layered media

,

Geophysics

,

48

,

1063

–

1081

.

Van der Neut

J.

,

Johnson

J.L.

,

van Wijk

K.

,

Singh

S.

,

Slob

E.

,

Wapenaar

K.

,

2017

.

A Marchenko equation for acoustic inverse source problems

,

J. acoust. Soc. Am.

,

141

(

6

),

4332

–

4346

.

Wapenaar

C.P.A.

,

Berkhout

A.J.

,

1989

.

Elastic Wave Field Extrapolation

,

Elsevier

.

Google Scholar

Google Preview

OpenURL Placeholder Text

WorldCat

Wapenaar

K.

,

2020

.

The Marchenko method for evanescent waves

,

Geophys. J. Int.

,

223

,

1412

–

1417

.

Wapenaar

K.

,

Slob

E.

,

2014

.

On the Marchenko equation for multicomponent single-sided reflection data

,

Geophys. J. Int.

,

199

,

1367

–

1371

.

Wapenaar

K.

,

Broggini

F.

,

Slob

E.

,

Snieder

R.

,

2013

.

Three-dimensional single-sided Marchenko inverse scattering, data-driven focusing, Green’s function retrieval, and their mutual relations

,

Phys. Rev. Lett.

,

110

,

084301

.

10.1103/PhysRevLett.110.084301

Wapenaar

K.

,

Thorbecke

J.

,

van der Neut

J.

,

Broggini

F.

,

Slob

E.

,

Snieder

R.

,

2014

.

Marchenko imaging

,

Geophysics

,

79

(

3

),

WA39

–

WA57

.

10.1190/geo2013-0302.1

Wapenaar

K.

,

Thorbecke

J.

,

van der Neut

J.

,

2016a

.

A single-sided homogeneous Green’s function representation for holographic imaging, inverse scattering, time-reversal acoustics and interferometric Green’s function retrieval

,

Geophys. J. Int.

,

205

,

531

–

535

.

Wapenaar

K.

,

van der Neut

J.

,

Slob

E.

,

2016b

.

Unified double- and single-sided homogeneous Green’s function representations

,

Proc. R. Soc. A

,

472

,

20160162

.

10.1098/rspa.2016.0162

Yoon

K.

,

Marfurt

K.J.

,

2006

.

Reverse-time migration using the Poynting vector

,

Explor. Geophys.

,

37

,

102

–

107

.

Zhang

L.

,

Slob

E.

,

2020

.

A fast algorithm for multiple elimination and transmission compensation in primary reflections

,

Geophys. J. Int.

,

221

,

371

–

377

.

APPENDIX A: DERIVATION OF THE ACOUSTIC WAVE FIELD REPRESENTATION

A1 Derivation of the representation of eq. (11)

We derive expressions for the coefficients a(x_R, ω) and b(x_R,ω) in the acoustic wave field representation of eq. (10). We do this by formulating two boundary conditions at |${{\partial \mathbb {D}}_R}$|⁠. First, we consider the acoustic pressure p(x, ω) at |${{\partial \mathbb {D}}_R}$|⁠. To this end, we evaluate eq. (10) for x at |${{\partial \mathbb {D}}_R}$|⁠. Using the focusing condition formulated in eq. (9) we thus obtain

$$\begin{eqnarray*} p({\bf x},\omega )|_{x_3=x_{3,R}}&=&\int _{{{\partial \mathbb {D}}_R}} \delta ({\bf x}_{\rm H}-{\bf x}_{{\rm H},R})a({\bf x}_R,\omega ){\rm d}{\bf x}_R +\int _{{{\partial \mathbb {D}}_R}} \delta ({\bf x}_{\rm H}-{\bf x}_{{\rm H},R})b({\bf x}_R,\omega ){\rm d}{\bf x}_R, \nonumber \\ &=&\lbrace a({\bf x},\omega )+b({\bf x},\omega )\rbrace _{x_3=x_{3,R}}, \end{eqnarray*}$$

(A1)

where we used x_R = (x_H,R, x_3,R). This is our first equation for the coefficients a(x_R, ω) and b(x_R, ω).

Next, we consider the vertical component of the particle velocity |$v_3({\bf x},\omega)$| at |${{\partial \mathbb {D}}_R}$|⁠. From the Fourier transform of eq. (1), using ρ_jk = δ_jkρ₀ at |${{\partial \mathbb {D}}_R}$|⁠, we obtain |$v_3({\bf x},\omega )=\frac{1}{i\omega \rho _0}\partial _3p({\bf x},\omega )$| for x at |${{\partial \mathbb {D}}_R}$|⁠. Substituting eq. (10) gives

$$\begin{eqnarray*} v_3({\bf x},\omega )&=&\frac{1}{i\omega \rho _0}\int _{{{\partial \mathbb {D}}_R}} \partial _3 F({\bf x},{\bf x}_R,\omega )a({\bf x}_R,\omega ){\rm d}{\bf x}_R + \frac{1}{i\omega \rho _0}\int _{{{\partial \mathbb {D}}_R}} \partial _3 F^*({\bf x},{\bf x}_R,\omega )b({\bf x}_R,\omega ){\rm d}{\bf x}_R, \end{eqnarray*}$$

(A2)

for x₃ = x_3,R. Applying the spatial Fourier transformation of eq. (40) to both sides of eq. (A2) gives

$$\begin{eqnarray*} \tilde{v}_3({\bf s},x_3,\omega )&=& \frac{1}{i\omega \rho _0}\int _{{{\partial \mathbb {D}}_R}} \partial _3\tilde{F}({\bf s},x_3,{\bf x}_R,\omega )a({\bf x}_R,\omega ){\rm d}{\bf x}_R + \frac{1}{i\omega \rho _0}\int _{{{\partial \mathbb {D}}_R}} \partial _3\tilde{F}^*(-{\bf s},x_3,{\bf x}_R,\omega )b({\bf x}_R,\omega ){\rm d}{\bf x}_R, \end{eqnarray*}$$

(A3)

for x₃ = x_{3, R}. At this depth level the focusing function is an upgoing field (see Fig. 1), hence it obeys the following one-way wave equation

$$\begin{eqnarray*} \partial _3\tilde{F}({\bf s},x_3,{\bf x}_R,\omega )|_{x_3=x_{3,R}}= -i\omega s_3\tilde{F}({\bf s},x_{3,R},{\bf x}_R,\omega ), \end{eqnarray*}$$

(A4)

with the vertical slowness s₃ defined as

$$\begin{eqnarray*} s_3=\left\lbrace \begin{array}{@{}l@{\quad }l@{}}\sqrt{1/c_0^2-{\bf s}\cdot {\bf s}}, &\mbox{for } {\bf s}\cdot {\bf s}\le 1/c_0^2\\ i\sqrt{{\bf s}\cdot {\bf s}-1/c_0^2}, &\mbox{for } {\bf s}\cdot {\bf s}\gt 1/c_0^2. \end{array}\right. \end{eqnarray*}$$

(A5)

The two expressions in eq. (A5) represent the situation for propagating and evanescent waves, respectively. Applying the spatial Fourier transformation of eqs (40) to eq. (9) we further have

$$\begin{eqnarray*} \tilde{F}({\bf s},x_{3,R},{\bf x}_R,\omega )=\exp \lbrace -i\omega {\bf s}\cdot {\bf x}_{{\rm H},R}\rbrace . \end{eqnarray*}$$

(A6)

Substitution of eqs (A4) and (A6) into eq. (A3) for x₃ = x_{3, R} gives

$$\begin{eqnarray*} \tilde{v}_3({\bf s},x_{3,R},\omega ) &=&-\frac{s_3}{\rho _0}\int _{{{\partial \mathbb {D}}_R}} \exp \lbrace -i\omega {\bf s}\cdot {\bf x}_{{\rm H},R}\rbrace a({\bf x}_R,\omega ){\rm d}{\bf x}_R +\frac{s_3^*}{\rho _0}\int _{{{\partial \mathbb {D}}_R}} \exp \lbrace -i\omega {\bf s}\cdot {\bf x}_{{\rm H},R}\rbrace b({\bf x}_R,\omega ){\rm d}{\bf x}_R \nonumber \\ &=&-\frac{s_3}{\rho _0}\tilde{a}({\bf s},x_{3,R},\omega )+\frac{s_3^*}{\rho _0}\tilde{b}({\bf s},x_{3,R},\omega ). \end{eqnarray*}$$

(A7)

Combining the spatial Fourier transform of eq. (A1) with eq. (A7) gives

$$\begin{eqnarray*} {\begin{pmatrix}\tilde{p}\\ \tilde{v}_3\end{pmatrix}}_{x_3=x_{3,R}}= {\begin{pmatrix}1 & 1 \\ s_3^*/\rho _0 & -s_3/\rho _0\end{pmatrix}} {\begin{pmatrix}\tilde{b}\\ \tilde{a}\end{pmatrix}}_{x_3=x_{3,R}}. \end{eqnarray*}$$

(A8)

For |${\bf s}\cdot {\bf s} \le 1/c_0^2$| at |${{\partial \mathbb {D}}_R}$| we have |$s_3^*=s_3$|⁠, see eq. (A5). Hence, for propagating waves, eq. (A8) is recognised as the well-known system that composes the total wave fields on the left-hand side from downgoing and upgoing fields on the right-hand side (Corones 1975; Ursin 1983; Fishman & McCoy 1984). Hence

$$\begin{eqnarray*} \tilde{b}({\bf s},x_{3,R},\omega )&=&\tilde{p}^+({\bf s},x_{3,R},\omega ), \end{eqnarray*}$$

(A9)

$$\begin{eqnarray*} \tilde{a}({\bf s},x_{3,R},\omega )&=&\tilde{p}^-({\bf s},x_{3,R},\omega ), \end{eqnarray*}$$

(A10)

for |${\bf s}\cdot {\bf s} \le 1/c_0^2$| at |${{\partial \mathbb {D}}_R}$|⁠, where |$\tilde{p}^+({\bf s},x_{3,R},\omega )$| and |$\tilde{p}^-({\bf s},x_{3,R},\omega )$| are downgoing and upgoing fields, respectively, at |${{\partial \mathbb {D}}_R}$|⁠. Transforming these expressions back to the space domain, using

$$\begin{eqnarray*} p^\pm ({\bf x}_R,\omega )=\frac{\omega ^2}{4\pi ^2}\int _{\mathbb {R}^2}\exp \lbrace i\omega {\bf s}\cdot {\bf x}_{{\rm H},R}\rbrace \tilde{p}^\pm ({\bf s},x_{3,R},\omega ){\rm d}{\bf s} \end{eqnarray*}$$

(A11)

and similar expressions for a(x_R, ω) and b(x_R, ω), gives

$$\begin{eqnarray*} b({\bf x}_R,\omega ) &\approx & p^+({\bf x}_R,\omega ), \end{eqnarray*}$$

(A12)

$$\begin{eqnarray*} a({\bf x}_R,\omega ) &\approx & p^-({\bf x}_R,\omega ). \end{eqnarray*}$$

(A13)

The approximation signs signify that evanescent waves are neglected at |${{\partial \mathbb {D}}_R}$| [since eqs (A9) and (A10) hold for propagating waves only, whereas the inverse Fourier transformation involves an integration along all horizontal slownesses]. Substitution of eqs (A12) and (A13) into eq. (10) gives eq. (11).

A2 Analysis of the integral in eq. (23)

We analyze the integral in eq. (23). We show that we can transfer the operator ∂_{3, R} from f to G^s, and that this is accompanied with a sign change. For a function of two space variables, u(x, x_R, ω), we define the spatial Fourier transform along the second space variable as

$$\begin{eqnarray*} \tilde{u}({\bf x},{\bf s},x_{3,R},\omega )= \int _{\mathbb {R}^2}u({\bf x},{\bf x}_{{\rm H},R},x_{3,R},\omega )\exp \lbrace i\omega {\bf s}\cdot {\bf x}_{{\rm H},R}\rbrace {\rm d}{\bf x}_{{\rm H},R}. \end{eqnarray*}$$

(A14)

Note the opposite sign in the exponential, compared with that in eq. (40). Using this Fourier transform and Parseval’s theorem, we obtain for the integral in eq. (23)

$$\begin{eqnarray*} \int _{{{\partial \mathbb {D}}_R}} \lbrace \partial _{3,R} f({\bf x},{\bf x}_R,\omega )\rbrace G^{\rm s}({\bf x}_S,{\bf x}_R,\omega ){\rm d}{\bf x}_R=\frac{\omega ^2}{4\pi ^2}\int _{\mathbb {R}^2} \lbrace \partial _{3,R} \,\,\tilde{\!\!f}({\bf x},-{\bf s},x_{3,R},\omega )\rbrace \tilde{G}^{\rm s}({\bf x}_S,{\bf s},x_{3,R},\omega ){\rm d}{\bf s}. \end{eqnarray*}$$

(A15)

Note that |$\,\,\tilde{\!\!f}({\bf x},-{\bf s},x_{3,R},\omega )$| is differentiated with respect to the focal point depth x_3,R, hence, the one-way wave equation gets a sign opposite to that in eq. (A4), that is |$\partial _{3,R} \,\,\tilde{\!\!f}({\bf x},-{\bf s},x_{3,R},\omega )=i\omega s_3\,\,\tilde{\!\!f}({\bf x},-{\bf s},x_{3,R},\omega )$|⁠, with s₃ defined in eq. (A5). We transfer iωs₃ to the Green’s function and use |$i\omega s_3\tilde{G}^{\rm s}({\bf x}_S,{\bf s},x_{3,R},\omega )=-\partial _{3,R}\tilde{G}^{\rm s}({\bf x}_S,{\bf s},x_{3,R},\omega )$| (which is a differentiation with respect to the source depth x_3,R). Making these substitutions in the right-hand side of eq. (A15) and applying Parseval’s theorem again gives

$$\begin{eqnarray*} -\frac{\omega ^2}{4\pi ^2}\int _{\mathbb {R}^2} \,\,\tilde{\!\!f}({\bf x},-{\bf s},x_{3,R},\omega )\partial _{3,R}\tilde{G}^{\rm s}({\bf x}_S,{\bf s},x_{3,R},\omega ){\rm d}{\bf s}=-\int _{{{\partial \mathbb {D}}_R}} f({\bf x},{\bf x}_R,\omega )\partial _{3,R} G^{\rm s}({\bf x}_S,{\bf x}_R,\omega ){\rm d}{\bf x}_R. \end{eqnarray*}$$

(A16)

Hence, we have transferred the operator ∂_{3, R} under the integral in eq. (23) from f to G^s, which involves a sign change.

APPENDIX B: DERIVATION OF THE ELASTODYNAMIC WAVE FIELD REPRESENTATION

B1 Derivation of the representation of eq. (54)

We derive expressions for the coefficients a(x_R, ω) and b(x_R, ω) in the elastodynamic wave field representation of eq. (53). We do this by formulating two boundary conditions at |${{\partial \mathbb {D}}_R}$|⁠. First, we consider the particle velocity vector v(x, ω) at |${{\partial \mathbb {D}}_R}$|⁠. To this end, we evaluate eq. (53) for x at |${{\partial \mathbb {D}}_R}$|⁠. Using the focusing condition formulated in eq. (51) we thus obtain

$$\begin{eqnarray*} {\bf v}({\bf x},\omega )|_{x_3=x_{3,R}}&=&\int _{{{\partial \mathbb {D}}_R}} \delta ({\bf x}_{\rm H}-{\bf x}_{{\rm H},R}){\bf a}({\bf x}_R,\omega ){\rm d}{\bf x}_R +\int _{{{\partial \mathbb {D}}_R}} \delta ({\bf x}_{\rm H}-{\bf x}_{{\rm H},R}){\bf b}({\bf x}_R,\omega ){\rm d}{\bf x}_R, \nonumber \\ &=& \lbrace {\bf a}({\bf x},\omega )+{\bf b}({\bf x},\omega )\rbrace _{x_3=x_{3,R}}. \end{eqnarray*}$$

(B1)

For the second boundary condition we analyze the traction vector |${\rm{\boldsymbol \tau}}_3({\bf x},\omega)$| at |${{\partial \mathbb {D}}_R}$|⁠. First we establish a relation between |${\rm{\boldsymbol \tau}}_3({\bf x},\omega)$| and v(x, ω) in the homogeneous isotropic upper half-space, including |${{\partial \mathbb {D}}_R}$|⁠. Using eq. (40), we transform v(x, ω) and |${\rm{\boldsymbol \tau}}_3({\bf x},\omega)$| for x₃ ≤ x_3,R to |$\tilde{\bf v}({\bf s},x_3,\omega )$| and |$\tilde{\rm{\boldsymbol \tau}}_3({\bf s},x_3,\omega )$|⁠, respectively. These fields can be related to vectors |$\tilde{\bf p}^+$| and |$\tilde{\bf p}^-$| containing downgoing and upgoing compressional and shear wave fields, according to

$$\begin{eqnarray*} {\begin{pmatrix}\tilde{\bf v} \\ -\tilde{\rm{\boldsymbol \tau}}_3 \end{pmatrix}} = {\begin{pmatrix}\tilde{\bf L}_1^+ & \tilde{\bf L}_1^- \\ \tilde{\bf L}_2^+ & \tilde{\bf L}_2^- \end{pmatrix}} {\begin{pmatrix}\tilde{\bf p}^+\\ \tilde{\bf p}^- \end{pmatrix}},\quad \mbox{for}\quad x_3\le x_{3,R}, \end{eqnarray*}$$

(B2)

(Kennett et al. 1978; Ursin 1983; Wapenaar & Berkhout 1989). Next, we define the downgoing and upgoing parts of |$\tilde{\bf v}$| as |$\tilde{\bf v}^\pm = \tilde{\bf L}_1^\pm \tilde{\bf p}^\pm$| and rewrite eq. (B2) as

$$\begin{eqnarray*} {\begin{pmatrix}\tilde{\bf v}\\ -\tilde{\rm{\boldsymbol \tau}}_3 \end{pmatrix}} ={\begin{pmatrix}\bf I & {\bf I}\\ \tilde{\bf D}^+ & \tilde{\bf D}^- \end{pmatrix}} {\begin{pmatrix}\tilde{\bf v}^+\\ \tilde{\bf v}^- \end{pmatrix}},\quad \mbox{for}\quad x_3\le x_{3,R}, \end{eqnarray*}$$

(B3)

with |$\tilde{\bf D}^\pm = \tilde{\bf L}_2^\pm (\tilde{\bf L}_1^\pm )^{-1}$|⁠. The matrices |$\tilde{\bf L}_1^\pm$| and |$\tilde{\bf L}_2^\pm$| in eq. (B2) are not uniquely defined. They depend on the chosen normalization of the fields contained in |$\tilde{\bf p}^+$| and |$\tilde{\bf p}^-$|⁠. However, independent of the normalization, the matrix |$\tilde{\bf D}^\pm$| in eq. (B3) is uniquely defined. It is given by

$$\begin{eqnarray*} \tilde{\bf D}^\pm ({\bf s})= \frac{\rho_0 c_S^2}{s_3^P s_3^S +{\bf s}\cdot {\bf s}}{\begin{pmatrix}\pm ((c_S^{-2}-s_2^2)s_3^P + s_2^2s_3^S) & \pm s_1s_2(s_3^P -s_3^S) & -s_1(c_S^{-2}-2(s_3^P s_3^S+{\bf s}\cdot {\bf s}))\\ \pm s_1s_2(s_3^P -s_3^S) & \pm ((c_S^{-2}-s_1^2)s_3^P + s_1^2s_3^S)& -s_2(c_S^{-2}-2(s_3^P s_3^S+{\bf s}\cdot {\bf s}))\\ s_1(c_S^{-2}-2(s_3^P s_3^S+{\bf s}\cdot {\bf s})) &s_2(c_S^{-2}-2(s_3^P s_3^S+{\bf s}\cdot {\bf s})) &\pm s_3^S c_S^{-2} \end{pmatrix}} \end{eqnarray*}$$

(B4)

with the vertical slownesses |$s_3^P$| and |$s_3^S$| for P and S waves, respectively, defined as

$$\begin{eqnarray*} s_3^{P,S}=\left\lbrace \begin{array}{@{}l@{\quad }l@{}}\sqrt{1/c_{P,S}^2-{\bf s}\cdot {\bf s}}, &\mbox{for } {\bf s}\cdot {\bf s}\le 1/c_{P,S}^2\\ i\sqrt{{\bf s}\cdot {\bf s}-1/c_{P,S}^2}, &\mbox{for } {\bf s}\cdot {\bf s}\gt 1/c_{P,S}^2, \end{array}\right. \end{eqnarray*}$$

(B5)

where c_P and c_S are the P- and S-wave velocities, respectively, of the upper half-space x₃ ≤ x_{3, R}. Applying the transform of eq. (40) to eq. (53) we obtain for x at |${{\partial \mathbb {D}}_R}$|

$$\begin{eqnarray*} \tilde{\bf v}({\bf s},x_{3,R},\omega )=\int _{{{\partial \mathbb {D}}_R}} \tilde{\bf F}({\bf s},x_{3,R},{\bf x}_R,\omega ){\bf a}({\bf x}_R,\omega ){\rm d}{\bf x}_R +\int _{{{\partial \mathbb {D}}_R}} \tilde{\bf F}^*(-{\bf s},x_{3,R},{\bf x}_R,\omega ){\bf b}({\bf x}_R,\omega ){\rm d}{\bf x}_R. \end{eqnarray*}$$

(B6)

Since |$\tilde{\bf F}({\bf s},x_{3,R},{\bf x}_R,\omega )$| is upgoing, the first term on the right-hand side is the upgoing velocity field |$\tilde{\bf v}^-({\bf s},x_{3,R},\omega )$| and the second term is, for propagating waves (i.e., for |${\bf s}\cdot {\bf s}\le 1/c_P^2$|⁠), the downgoing velocity field |$\tilde{\bf v}^+({\bf s},x_{3,R},\omega )$|⁠. Hence, using eq. (B3) we obtain for the transformed traction vector

$$\begin{eqnarray*} -\tilde{\rm{\boldsymbol \tau}}_3({\bf s},x_{3,R},\omega )=\tilde{\bf D}^-({\bf s})\int _{{{\partial \mathbb {D}}_R}} \tilde{\bf F}({\bf s},x_{3,R},{\bf x}_R,\omega ){\bf a}({\bf x}_R,\omega ){\rm d}{\bf x}_R +\tilde{\bf D}^+({\bf s})\int _{{{\partial \mathbb {D}}_R}} \tilde{\bf F}^*(-{\bf s},x_{3,R},{\bf x}_R,\omega ){\bf b}({\bf x}_R,\omega ){\rm d}{\bf x}_R, \end{eqnarray*}$$

(B7)

for |${\bf s}\cdot {\bf s}\le 1/c_P^2$| at |${{\partial \mathbb {D}}_R}$|⁠. Applying the transform of eq. (40) to the focusing condition of eq. (51) gives

$$\begin{eqnarray*} \tilde{\bf F}({\bf s},x_{3,R},{\bf x}_R,\omega )={\bf I}\exp \lbrace -i\omega {\bf s}\cdot {\bf x}_{{\rm H},R}\rbrace . \end{eqnarray*}$$

(B8)

Substituting this into eq. (B7) we obtain

$$\begin{eqnarray*} -\tilde{\rm{\boldsymbol \tau}}_3({\bf s},x_{3,R},\omega )= \tilde{\bf D}^-({\bf s}) \tilde{\bf a}({\bf s},x_{3,R},\omega )+ \tilde{\bf D}^+({\bf s}) \tilde{\bf b}({\bf s},x_{3,R},\omega ), \end{eqnarray*}$$

(B9)

for |${\bf s}\cdot {\bf s}\le 1/c_P^2$| at |${{\partial \mathbb {D}}_R}$|⁠. Combining this equation with the Fourier transform of eq. (B1) yields

$$\begin{eqnarray*} {\begin{pmatrix}\tilde{\bf v}\\ -\tilde{\rm{\boldsymbol \tau}}_3 \end{pmatrix}}_{x_3=x_{3,R}} ={\begin{pmatrix}\bf I & {\bf I}\\ \tilde{\bf D}^+ & \tilde{\bf D}^- \end{pmatrix}} {\begin{pmatrix}\tilde{\bf b}\\ \tilde{\bf a} \end{pmatrix}}_{x_3=x_{3,R}}, \end{eqnarray*}$$

(B10)

for |${\bf s}\cdot {\bf s}\le 1/c_P^2$| at |${{\partial \mathbb {D}}_R}$|⁠. Comparing this with eq. (B3) we conclude

$$\begin{eqnarray*} \tilde{\bf b}({\bf s},x_{3,R},\omega )&=&\tilde{\bf v}^+({\bf s},x_{3,R},\omega ), \end{eqnarray*}$$

(B11)

$$\begin{eqnarray*} \tilde{\bf a}({\bf s},x_{3,R},\omega )&=&\tilde{\bf v}^-({\bf s},x_{3,R},\omega ), \end{eqnarray*}$$

(B12)

for |${\bf s}\cdot {\bf s}\le 1/c_P^2$| at |${{\partial \mathbb {D}}_R}$|⁠. Transforming these expressions back to the space domain gives

$$\begin{eqnarray*} {\bf b}({\bf x}_R,\omega ) &\approx & {\bf v}^+({\bf x}_R,\omega ), \end{eqnarray*}$$

(B13)

$$\begin{eqnarray*} {\bf a}({\bf x}_R,\omega ) &\approx & {\bf v}^-({\bf x}_R,\omega ). \end{eqnarray*}$$

(B14)

The approximation signs signify that evanescent waves are neglected at |${{\partial \mathbb {D}}_R}$|⁠. Substitution of eqs (B13) and (B14) into eq. (53) gives eq. (54).

B2 Derivation of the modified elastodynamic Green’s function

We derive a modified elastodynamic Green’s function Γ(x, x_S, ω) (with x_3,S = x_3,R − ϵ), such that for x at |${{\partial \mathbb {D}}_R}$|⁠, that is just below the source level, the downgoing part of Γ(x, x_S, ω) obeys eq. (58), that is

$$\begin{eqnarray*} \lim _{x_{3}\downarrow x_{3,S}}{\rm{\boldsymbol \Gamma}}^+ ({\bf x},{\bf x}_S,\omega )={\bf I}\delta ({\bf x}_{\rm H}-{\bf x}_{{\rm H},S}). \end{eqnarray*}$$

(B15)

To this end we first investigate the properties of the downgoing part of G(x, x_S, ω) defined in eqs (56) and (57), just below the source level. Consider the inverse of eq. (B3)

$$\begin{eqnarray*} {\begin{pmatrix}\tilde{\bf v}^+\\ \tilde{\bf v}^- \end{pmatrix}}= {\begin{pmatrix}-\tilde{\rm{\boldsymbol \Delta}}^{-1}\tilde{\bf D}^-& \tilde{\rm{\boldsymbol \Delta}}^{-1}\\ \tilde{\rm{\boldsymbol \Delta}}^{-1}\tilde{\bf D}^+ & -\tilde{\rm{\boldsymbol \Delta}}^{-1} \end{pmatrix}} {\begin{pmatrix}\tilde{\bf v} \\ -\tilde{\rm{\boldsymbol \tau}}_3\end{pmatrix}}, \quad \mbox{for}\quad x_3\le x_{3,R}, \end{eqnarray*}$$

(B16)

with

$$\begin{eqnarray*} \tilde{\rm{\boldsymbol \Delta}}=\tilde{\bf D}^+- \tilde{\bf D}^-. \end{eqnarray*}$$

(B17)

The upper-right matrix in eq. (B16), |$\tilde{\rm{\boldsymbol \Delta}}^{-1}$|⁠, gives the relation between |$-\tilde{\rm{\boldsymbol \tau}}_3$| and the downgoing velocity vector |$\tilde{\bf v}^+$|⁠. The same matrix transforms a unit force source in a homogeneous half-space into the downgoing part of the Green’s function just below this source, hence

$$\begin{eqnarray*} \lim _{x_3\downarrow x_{3,S}}\tilde{\bf G}^+({\bf s},x_3,{\bf 0},x_{3,S},\omega )&=& \tilde{\rm{\boldsymbol \Delta}}^{-1}({\bf s})\nonumber\\ &=&\frac{1}{2\rho_0 }{\begin{pmatrix}\frac{s_1^2}{s_3^P }+\Bigl (\frac{1}{c_S^2}-s_1^2\Bigr )\frac{1}{s_3^S} & \Bigl (\frac{1}{s_3^P }-\frac{1}{s_3^S}\Bigr )s_1s_2 & 0\\ \Bigl (\frac{1}{s_3^P }-\frac{1}{s_3^S}\Bigr )s_1s_2 & \frac{s_2^2}{s_3^P }+\Bigl (\frac{1}{c_S^2}-s_2^2\Bigr )\frac{1}{s_3^S} & 0\\ 0 & 0 & s_3^P +\frac{{\bf s}\cdot {\bf s}}{s_3^S} \end{pmatrix}}. \end{eqnarray*}$$

(B18)

In eq. (B18) the source is located at (0, x_3,S). Next, we consider G(x, x_S, ω) for a laterally shifted source position (x_H,S, x_3,S). Applying a spatial Fourier transform along the horizontal source coordinate x_{H, S}, using eq. (A14) with x_H,R replaced by x_H,S, yields |$\tilde{\bf G}({\bf x},{\bf s},x_{3,S},\omega )$|⁠. For the downgoing part just below the source we obtain a phase-shifted version of the Green’s function of eq. (B18), according to

$$\begin{eqnarray*} \lim _{x_3\downarrow x_{3,S}}\tilde{\bf G}^+({\bf x},{\bf s},x_{3,S},\omega ) =\tilde{\rm{\boldsymbol \Delta}}^{-1}({\bf s})\exp \lbrace i\omega {\bf s}\cdot {\bf x}_{\rm H}\rbrace . \end{eqnarray*}$$

(B19)

Comparing this with the desired condition of eq. (B15) suggests to define the modified Green’s function (for arbitrary x) as

$$\begin{eqnarray*} \tilde{\rm{\boldsymbol \Gamma}}({\bf x},{\bf s},x_{3,S},\omega )=\tilde{\bf G}({\bf x},{\bf s},x_{3,S},\omega )\tilde{\rm{\boldsymbol \Delta}}({\bf s}), \end{eqnarray*}$$

(B20)

such that

$$\begin{eqnarray*} \lim _{x_{3}\downarrow x_{3,S}}\tilde{\rm{\boldsymbol \Gamma}}^+ ({\bf x},{\bf s},x_{3,S},\omega ) ={\bf I}\exp \lbrace i\omega {\bf s}\cdot {\bf x}_{\rm H}\rbrace . \end{eqnarray*}$$

(B21)

The inverse Fourier transform from s to x_{H, S} gives indeed eq. (B15).

We define the reflection response |$\tilde{\bf R}({\bf x}_R,{\bf s},x_{3,S},\omega )$| of the medium below |${{\partial \mathbb {D}}_R}$| as the upgoing part of the modified Green’s function |$\tilde{\rm{\boldsymbol \Gamma}}({\bf x}_R,{\bf s},x_{3,S},\omega )$|⁠, with x_R at |${{\partial \mathbb {D}}_R}$|⁠, hence

$$\begin{eqnarray*} \tilde{\bf R}({\bf x}_R,{\bf s},x_{3,S},\omega )=\tilde{\rm{\boldsymbol \Gamma}}^-({\bf x}_R,{\bf s},x_{3,S},\omega ), \end{eqnarray*}$$

(B22)

or, using eq. (B20),

$$\begin{eqnarray*} \tilde{\bf R}({\bf x}_R,{\bf s},x_{3,S},\omega )&=&\tilde{\bf G}^-({\bf x}_R,{\bf s},x_{3,S},\omega )\tilde{\rm{\boldsymbol \Delta}}({\bf s})\nonumber \\ &=&\tilde{\bf G}^{\rm s}({\bf x}_R,{\bf s},x_{3,S},\omega )\tilde{\rm{\boldsymbol \Delta}}({\bf s}), \end{eqnarray*}$$

(B23)

where superscript s stands for scattered. The inverse Fourier transform of eq. (B22) from s to x_H,S yields eq. (59).

B3 Derivation of the representation of eq. (61)

To obtain a representation for G(x, x_S, ω) we start by transforming all terms in eq. (60) along x_H,S, using eq. (A14), with x_H,R replaced by x_H,S, hence

$$\begin{eqnarray*} \tilde{\rm{\boldsymbol \Gamma}}({\bf x},{\bf s},x_{3,S},\omega )&=&\int _{{{\partial \mathbb {D}}_R}} {\bf F}({\bf x},{\bf x}_R,\omega )\tilde{\bf R}({\bf x}_R,{\bf s},x_{3,S},\omega ){\rm d}{\bf x}_R +\tilde{\bf F}^*({\bf x},-{\bf s},x_{3,S},\omega ),\quad\mbox{for}\quad x_3 \ge x_{3,R}. \end{eqnarray*}$$

(B24)

We introduce a modified focusing function |$\tilde{\bf f}({\bf x},{\bf s},x_{3,S},\omega )$| via

$$\begin{eqnarray*} \tilde{\bf F}({\bf x},{\bf s},x_{3,S},\omega )=\tilde{\bf f}({\bf x},{\bf s},x_{3,S},\omega )\tilde{\rm{\boldsymbol \Delta}}({\bf s}). \end{eqnarray*}$$

(B25)

According to eq. (B18) we have for propagating waves (i.e. for |${\bf s}\cdot {\bf s}\le 1/c_P^2$|⁠)

$$\begin{eqnarray*} \tilde{\rm{\boldsymbol \Delta}}({\bf s})=\tilde{\rm{\boldsymbol \Delta}}(-{\bf s})=\tilde{\rm{\boldsymbol \Delta}}^*({\bf s})=\tilde{\rm{\boldsymbol \Delta}}^t({\bf s}). \end{eqnarray*}$$

(B26)

Hence, for |$\tilde{\bf F}^*({\bf x},-{\bf s},x_{3,S},\omega )$| we obtain

$$\begin{eqnarray*} \tilde{\bf F}^*({\bf x},-{\bf s},x_{3,S},\omega )=\tilde{\bf f}^*({\bf x},-{\bf s},x_{3,S},\omega )\tilde{\rm{\boldsymbol \Delta}}({\bf s}). \end{eqnarray*}$$

(B27)

Multiplying all terms in eq. (B24) from the right by |$\tilde{\rm{\boldsymbol \Delta}}^{-1}({\bf s})$|⁠, using eqs (B20), (B23) and (B27), and transforming the resulting expression back from s to x_H,S gives

$$\begin{eqnarray*} {\bf G}({\bf x},{\bf x}_S,\omega )=\int _{{{\partial \mathbb {D}}_R}} {\bf F}({\bf x},{\bf x}_R,\omega ){\bf G}^{\rm s}({\bf x}_R,{\bf x}_S,\omega ){\rm d}{\bf x}_R +{\bf f}^*({\bf x},{\bf x}_S,\omega ), \quad \mbox{for}\quad x_3 \ge x_{3,R}. \end{eqnarray*}$$

(B28)

We modify the integral step by step. First we use source–receiver reciprocity for the scattered Green’s function G^s(x_R, x_S, ω) and we apply Parseval’s theorem. We thus obtain for the integral in eq. (B28)

$$\begin{eqnarray*} \frac{\omega ^2}{4\pi ^2}\int _{\mathbb {R}^2} \tilde{\bf F}({\bf x},-{\bf s},x_{3,R},\omega )\lbrace \tilde{\bf G}^{\rm s}({\bf x}_S,{\bf s},x_{3,R},\omega )\rbrace ^t{\rm d}{\bf s}. \end{eqnarray*}$$

(B29)

Substituting eq. (B25), using eq. (B26), gives

$$\begin{eqnarray*} \frac{\omega ^2}{4\pi ^2}\int _{\mathbb {R}^2} \tilde{\bf f}({\bf x},-{\bf s},x_{3,R},\omega )\lbrace \tilde{\bf G}^{\rm s}({\bf x}_S,{\bf s},x_{3,R},\omega )\tilde{\rm{\boldsymbol \Delta}}({\bf s})\rbrace ^t{\rm d}{\bf s}. \end{eqnarray*}$$

(B30)

Using eq. (B23) this gives

$$\begin{eqnarray*} \frac{\omega ^2}{4\pi ^2}\int _{\mathbb {R}^2} \tilde{\bf f}({\bf x},-{\bf s},x_{3,R},\omega )\lbrace \tilde{\bf R}({\bf x}_S,{\bf s},x_{3,R},\omega )\rbrace ^t{\rm d}{\bf s}. \end{eqnarray*}$$

(B31)

Applying Parseval’s theorem again and inserting the resulting integral in eq. (B28) yields eq. (61). It has been derived without applying decomposition in the lower half-space, but it excludes the contribution from waves that are evanescent at |${{\partial \mathbb {D}}_R}$|⁠.

This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://dbpia.nl.go.kr/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Download all slides

Month:	Total Views:
June 2021	57
July 2021	38
August 2021	22
September 2021	46
October 2021	34
November 2021	13
December 2021	10
January 2022	3
February 2022	12
March 2022	12
April 2022	12
May 2022	13
June 2022	13
July 2022	4
August 2022	11
September 2022	10
October 2022	11
November 2022	1
December 2022	1
February 2023	2
March 2023	7
April 2023	2
May 2023	1
July 2023	8
August 2023	3
September 2023	6
October 2023	1
November 2023	14
December 2023	4
January 2024	25
February 2024	15
March 2024	20
April 2024	4
May 2024	22
June 2024	8
July 2024	13
August 2024	16
September 2024	16
October 2024	13
November 2024	8
December 2024	6
January 2025	8
February 2025	5
March 2025	11
April 2025	13
May 2025	2

Article Contents

Green’s function representations for Marchenko imaging without up/down decomposition

SUMMARY

1 INTRODUCTION

2 ACOUSTIC WAVE FIELD REPRESENTATION

3 ACOUSTIC GREEN’S FUNCTION REPRESENTATIONS

3.1 Representation for the acoustic dipole Green’s function

3.2 Representation for the acoustic monopole Green’s function

3.3 Representation for the homogeneous acoustic Green’s function

3.4 Numerical examples

4 ELASTODYNAMIC WAVE FIELD REPRESENTATION

5 ELASTODYNAMIC GREEN’S FUNCTION REPRESENTATIONS

5.1 Representation for a modified elastodynamic Green’s function

5.2 Representation for the elastodynamic Green’s function

5.3 Representation for the homogeneous elastodynamic Green’s function

6 CONCLUSIONS

ACKNOWLEDGEMENTS

DATA AVAILABILITY

REFERENCES

APPENDIX A: DERIVATION OF THE ACOUSTIC WAVE FIELD REPRESENTATION

A1 Derivation of the representation of eq. (11)

A2 Analysis of the integral in eq. (23)

APPENDIX B: DERIVATION OF THE ELASTODYNAMIC WAVE FIELD REPRESENTATION

B1 Derivation of the representation of eq. (54)

B2 Derivation of the modified elastodynamic Green’s function

B3 Derivation of the representation of eq. (61)

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

Article Contents

Green’s function representations for Marchenko imaging without up/down decomposition

SUMMARY

1 INTRODUCTION

2 ACOUSTIC WAVE FIELD REPRESENTATION

3 ACOUSTIC GREEN’S FUNCTION REPRESENTATIONS

3.1 Representation for the acoustic dipole Green’s function

3.2 Representation for the acoustic monopole Green’s function

3.3 Representation for the homogeneous acoustic Green’s function

3.4 Numerical examples

4 ELASTODYNAMIC WAVE FIELD REPRESENTATION

5 ELASTODYNAMIC GREEN’S FUNCTION REPRESENTATIONS

5.1 Representation for a modified elastodynamic Green’s function

5.2 Representation for the elastodynamic Green’s function

5.3 Representation for the homogeneous elastodynamic Green’s function

6 CONCLUSIONS

ACKNOWLEDGEMENTS

DATA AVAILABILITY

REFERENCES

APPENDIX A: DERIVATION OF THE ACOUSTIC WAVE FIELD REPRESENTATION

A1 Derivation of the representation of eq. (11)

A2 Analysis of the integral in eq. (23)

APPENDIX B: DERIVATION OF THE ELASTODYNAMIC WAVE FIELD REPRESENTATION

B1 Derivation of the representation of eq. (54)

B2 Derivation of the modified elastodynamic Green’s function

B3 Derivation of the representation of eq. (61)

Citations

Views

Altmetric

Email alerts

Astrophysics Data System

Citing articles via

Latest

Most Read

Most Cited

This Feature Is Available To Subscribers Only