SUMMARY

As seismic imaging moves towards the imaging of more complex media, properly modelling elastic effects in the subsurface is becoming of increasing interest. In this context, elastic wave conversion, where acoustic, pressure (P-) waves are converted into elastic, shear (S-) waves, is of great importance. Accounting for these wave conversions, in the framework of forward and inverse modelling of elastic waves, is crucial to creating accurate images of the subsurface in complex media. The underlying mechanism of wave conversion is well understood and described by the Zoeppritz equations. However, as these equations are highly nonlinear, approximations are commonly used. The most well-known of these approximations is Shuey’s approximation. However, this approximation only holds for small angles and small contrasts, making it insufficient for realistic forward and inverse modelling scenarios, where angles and contrasts may be large. In this paper we present a novel set of approximations, based on Taylor expansions of the Zoeppritz equations, which we name the extended Shuey approximations. We examine the quality of these approximations to the Zoeppritz equations and compare them to existing approximations described in literature. We then apply these extended Shuey approximations to the elastic full-wavefield modelling algorithm for a simple, synthetic, 1.5-D example, where we show that we can accurately model the P- and S-wavefields in a forward modelling case. Finally, we apply our approximations to the elastic full-wavefield migration algorithm for a simple, synthetic, 1.5-D example, where we show that we can recover an accurate image in an inverse modelling case.

1 INTRODUCTION

Properly modelling the propagation and scattering of elastic waves throughout the subsurface of the Earth is an important aspect in many seismic applications. In this context, the process of wave conversion, which occurs due to the elastic nature of the Earth’s subsurface, has received increased interest in recent years. Wave conversion is a process where acoustic, pressure waves, known as P waves, convert into elastic, shear waves, known as S waves. This conversion takes place when a wavefield strikes an interface between two different media. This process is highly relevant in multiple areas, such as subsurface imaging in complex media, anisotropy analysis, and reservoir characterization (Stewart et al. 2003).

The underlying mechanism by which wave conversion occurs is well understood in literature, and exact reflection and transmission coefficients can be derived for flat interfaces between two media. These exact reflection and transmission coefficients are known as the Zoeppritz equations (see Zoeppritz 1919 and Aki & Richards 2002). However, while these equations are exact, they are also notoriously difficult to work with. This is due to the nonlinearity of the Zoeppritz equations, which presents itself in two different ways. First of all, the Zoeppritz equations depend on the medium properties of the two media in a very nonlinear fashion. This makes it difficult to use these equations in inversion settings, where one wishes to recover the medium properties from the measurement data. Secondly, as an additional complication, the Zoeppritz equations need to account for critical angles, where the reflection and transmission coefficients diverge, which can cause a multitude of issues.

Due to these challenges, when working with the Zoeppritz equations in practice, approximations are commonly applied. By far the most well-known of these approximations is Shuey’s approximation (Shuey 1985). With Shuey’s approximation the Zoeppritz equations are linearized, which makes them simple to invert, and removes the critical points from the Zoeppritz equations. Due to this ease of inversion, it has become the standard in amplitude variation with offset (AVO) analysis. However, Shuey’s approximation is only accurate for small angles and weak contrasts, which limits its applicability. Also, it only describes the |$PP$|-reflection coefficient, which means that wave conversions are not taken into account.

To address these issues, alternative approximations have been developed. The most influential of these alternative approximations are the linearized approximations of Aki and Richards (Aki & Richards 2002). These approximations do address wave conversions, extending their applicability. However, these approximations are also only applicable in situations where the contrasts are weak. Other approximations, such as by Wang (Wang 1999), have also been developed. However, most of these approximations focus on the |$PP$|-reflection coefficient, as this coefficient has been of the most interest historically.

In recent years, the forward and inverse modelling of wave paths including wave conversions have become of greater interest. As the focus of seismic imaging shifts towards more challenging scenarios, specifically those with strong reflectors (Jones & Davison 2014), or towards cases where the goal is the recovery of the elastic parameters (Prieux et al. 2013), it is important to accurately describe these wave conversions. To do this, high quality approximations of the full set of Zoeppritz equations are required.

In this paper, we present a novel set of accurate approximations to the full set of Zoeppritz equations, which we name the extended Shuey’s approximations, based on applying Taylor expansions to the Zoeppritz equations. We show that the result of this approximation reduces to the well-known Shuey’s approximation in the linear case. In addition, we show how it can be extended to retrieve more accurate approximations for both the |$PP$|-reflection coefficient as well as the other reflection and transmission coefficients. We then apply our approximation to a simple forward and inverse modelling scenario to give a proof-of-concept of how this approximation can improve the accuracy of modelling wave conversions in practice.

2 THEORY

This section is split into three parts. First, in Section 2.1, we derive the extended Shuey’s approximations, where we focus on the |$R_{PP}$| term, and show that they are a natural extension of the traditional Shuey approximation. Next, in Section 2.2, we apply this approximation to a simple, elastic, forward modelling algorithm called elastic full-wavefield modelling (FWMod). Finally, in Section 2.3 we discuss the inverse problem associated with the forward model of Section 2.2, referred to as elastic full-wavefield migration.

2.1 Derivation of the extended Shuey’s approximations

We begin by deriving the extended Shuey’s approximations of the Zoeppritz equations. Consider a flat, laterally invariant interface between two media located at a depth level |$z = z_n$|⁠. Following the notation used in Aki & Richards (2002), each medium is characterized by its P-wave velocity |$\alpha$|⁠, S-wave velocity |$\beta$|⁠, and mass density |$\rho$|⁠. The subscript 1 is used to denote properties within the medium above the interface (at a depth |$z < z_n$|⁠), while properties within the medium below the interface (at a depth |$z > z_n$|⁠) are denoted by the subscript 2.

In this situation, the reflection and transmission coefficients are given analytically by the Zoeppritz equations. As these equations are unwieldy to write down directly, we once again follow the notation used in Aki & Richards (2002) and introduce a set of intermediate variables. We begin by introducing the dimensionless contrast parameters

(1)

Rewriting the medium parameters above and below the interface in terms of these contrast parameters gives

(2)

where

(3)

Next, we introduce the intermediate variables a, b, c, and d, which we write in terms of the contrast parameters |$c_\alpha$|⁠, |$c_\beta$|⁠, and |$c_\rho$|⁠, viz.

(4)

with generalized angle of incidence |$\theta$|⁠, and where we use the shorthand notation |$\hat{V} = \bar{\beta }/\bar{\alpha }$|⁠. Also note the use of the |$\, \tilde{}\,$| symbol in eq. (4) to indicate dimensionless variables.

In eq. (4) we have introduced the generalized angle of incidence |$\theta$|⁠, which is defined implicitly using Snell’s law, viz.

(5)

where |$i_1$| and |$i_2$| are the angles of incidence and refraction of the P waves in the medium above and below the interface, respectively, while |$j_1$| and |$j_2$| are the angles of incidence and refraction of the S waves in the medium above and below the interface. Using eq. (5), along with the trigonometric identity |${\sin ^2}\left( \theta \right) + {\cos ^2}\left( \theta \right) = 1$|⁠, we write the associated cosine terms as

(6)

Next, we introduce the additional intermediate variables

(7)

where we have once again indicated dimensionless variables with the |$\displaystyle \, \tilde{}\,$| symbol. Finally, using eqs (2), (4), (6) and (7), the Zoeppritz |$\displaystyle PP$|-reflection coefficient for waves coming from above equals

(8)

As expected, |$\displaystyle R^\cup _{PP}\left( {\sin \left( \theta \right),{c_\alpha },{c_\beta },{c_\rho }} \right)$| is dimensionless, and only depends on the dimensionless variables |$\displaystyle c_\alpha$|⁠, |$\displaystyle c_\beta$|⁠, |$\displaystyle c_\rho$|⁠, |$\displaystyle \hat{V}$|⁠, as well as the angle |$\displaystyle \sin \left( \theta \right)$|⁠. In a similar way, the expressions for the remaining reflection and transmission coefficients can be found.

The goal of our work is to find accurate approximations of the Zoeppritz reflection and transmission coefficients. To do this, we begin by noting that |$\displaystyle c_\alpha$|⁠, |$\displaystyle c_\beta$|⁠, |$\displaystyle c_\rho$| and |$\displaystyle \sin \left( \theta \right)$| all have a magnitude smaller than 1 for normal seismic applications. Therefore, a natural approximation is given by taking the Taylor expansion of eq. (8) with respect to |$\displaystyle c_\alpha$|⁠, |$\displaystyle c_\beta$|⁠, |$\displaystyle c_\rho$| and |$\displaystyle \sin \left( \theta \right)$|⁠. For a function of multiple variables, the Taylor expansion is given by

(9)

where the argument |$\displaystyle {\left. {} \right|_{\left( {\mathbf {0}} \right)}}$| is used as shorthand for |$\displaystyle \sin \left( \theta \right) = 0$|⁠, |$\displaystyle {c_\alpha } = 0$|⁠, |$\displaystyle {c_\beta } = 0$| and |$\displaystyle {c_\rho } = 0$|⁠. As the Zoeppritz equations are well-behaved up to the critical angle |$\displaystyle \sin \left( {{\theta _c}} \right) = \bar{\alpha }/\max \left( {{\alpha _1},{\alpha _2}} \right)$|⁠, we know that the Taylor series of eq. (9) will converge to |$\displaystyle R_{PP}$| for angles |$\displaystyle \theta < \theta _c$|⁠, given enough terms. For ease of interpretation, we order the terms of eq. (9) in terms of the power |$\displaystyle n$| of the |$\displaystyle \sin ^n \left( \theta \right)$| term and the total power |$\displaystyle \lambda = m + k + l$| of the contrast terms |$\displaystyle c_\alpha ^m$|⁠, |$\displaystyle c_\beta ^k$| and |$\displaystyle c_\rho ^l$|⁠, hence

(10)

with

(11)

While eq. (11) is difficult to evaluate by hand, by using mathematical software such as Maple it is straightforward to calculate the necessary terms. The first few terms of eq. (11) are given by

(12)

Also note that |$\displaystyle {\left( {\tilde{R}_{PP}^ \cup } \right)_\lambda ^n} = 0$| for odd values of |$\displaystyle n$|⁠, as we know that |$\displaystyle R_{PP}^ \cup$| is an even function with respect to |$\displaystyle \theta$| [meaning that it contains no |$\displaystyle {{\sin }^n}\left( \theta \right)$| terms for odd values of |$\displaystyle n$|].

Let us now examine the terms of eq. (11) in more detail. We begin by writing the two-term Shuey approximation (Shuey 1985), which, using the notation used in this paper, is denoted as

(13)

However, this is exactly equal to the approximation given by taking |$\displaystyle N = 2$| and |$\displaystyle \Lambda = 1$| in eq. (10). Continuing in this manner, we examine the three-term Shuey approximation (Shuey 1985), which can be written as

(14)

We note that the Taylor expansion of |$\displaystyle {{{\tan }^2}\left( \theta \right) - {{\sin }^2}\left( \theta \right)}$| is given by

(15)

which, once again, matches the approximation given in eq. (10) for |$\displaystyle N = 8$| and |$\displaystyle \Lambda = 1$|⁠. Therefore, we conclude that the approximation presented here can be viewed as an extension to Shuey’s approximation.

In a similar way, we can construct approximations to the rest of the Zoeppritz equations. For example, the |$\displaystyle SP$|-reflection coefficient |$\displaystyle {R^\cup _{SP}}$|⁠, which describes the conversion of |$\displaystyle P$| waves to |$\displaystyle S$| waves for waves coming from above, is given by

(16)

Once again, we take the Taylor expansion

(17)

with

(18)

The first few terms of eq. (18) are given by

(19)

In a similar way as before for the |$\displaystyle R_{PP}^ \cup$| terms, we note that |$\displaystyle {\left( {\tilde{R}_{SP}^ \cup } \right)_\lambda ^n} = 0$| for even values of |$\displaystyle n$|⁠, as we know that |$\displaystyle R_{SP}^ \cup$| is an odd function with respect to |$\displaystyle \theta$|⁠.

A thorough analysis of the quality of these extended Shuey approximations are given in Sections 3.1 and 4. In Section 3.1 we examine how many terms of the Taylor expansion are required to produce an accurate approximation. In Section 4 we compare the approximation introduced here with well-known approximations from literature, such as the standard Shuey’s approximation and the linearized approximation of Aki & Richards (2002).

2.2 Elastic full wavefield modelling

We now apply these extended Shuey approximations, described in Section 2.1, to the elastic full-wavefield modelling (FWMod) algorithm (Berkhout 2014a). For simplicity’s sake we consider a 1.5-D medium: a laterally homogeneous, but potentially vertically changing medium, characterized by a |$\displaystyle P$|- and |$\displaystyle S$|-wave velocity profile |$\displaystyle \alpha \left( z \right)$| and |$\displaystyle \beta \left( z \right)$|⁠, respectively, and a mass density profile |$\displaystyle \rho \left( z \right)$|⁠. As the medium is laterally homogeneous, we will work in the |$\displaystyle \left( k_x , \omega \right)$|-domain, where |$\displaystyle k_x$| is the lateral wave number and |$\displaystyle \omega$| is the angular frequency. Our goal is to describe the propagation and scattering of the |$\displaystyle P$| and |$\displaystyle S$| wavefields |$\displaystyle p_P \left( k_x, \omega , z \right)$| and |$\displaystyle p_S \left( k_x, \omega , z \right)$|⁠, respectively.

Assume the existence of an interface at a depth level |$\displaystyle z_n$|⁠, illustrated in Fig. 1. At a depth level |$\displaystyle z_n^-$| right above the interface, the |$\displaystyle P$| and |$\displaystyle S$| wavefields are given by

(20)

where we have split the |$\displaystyle P$| and |$\displaystyle S$| wavefields |$\displaystyle {p_{P/S}}\left( {{k_x},\omega ,z_n^ - } \right)$| into downgoing components, |$\displaystyle p_{P/S}^ + \left( {{k_x},\omega ,z_n^ - } \right)$|⁠, travelling towards the interface, and upgoing components, |$\displaystyle q_{P/S}^ - \left( {{k_x},\omega ,z_n^ - } \right)$|⁠, travelling away from the interface. Similarly, at a depth level |$\displaystyle z_n^+$| right below the interface, the P and S wavefields are given by

(21)

where we have once again split the |$\displaystyle P$| and |$\displaystyle S$| wavefields into upgoing components, |$\displaystyle p_{P/S}^ - \left( {{k_x},\omega ,z_n^ + } \right)$|⁠, travelling towards the interface, and downgoing components, |$\displaystyle q^ +_{P/S} \left( {{k_x},\omega ,z_n^ + } \right)$|⁠, travelling away from the interface. Note the notation used here. Superscripts above the field quantities |$\displaystyle p$| and |$\displaystyle q$| are used to indicate the direction of propagation, where |$\displaystyle +$| and |$\displaystyle -$| denote downwards and upwards propagation, respectively. The characters |$\displaystyle p$| and |$\displaystyle q$| are used to denote wavefields propagating towards or away from the interface, respectively. Finally, the subscripts |$\displaystyle P$| and |$\displaystyle S$| are used to differentiate between |$\displaystyle P$| and |$\displaystyle S$| wavefields, respectively.

Schematic representation of the wavefields and operators at depth level $\displaystyle z_n$ and its neighbouring levels.
Figure 1.

Schematic representation of the wavefields and operators at depth level |$\displaystyle z_n$| and its neighbouring levels.

At the interface at |$\displaystyle z = z_n$|⁠, the relationships between these wavefields are given by

(22)
(23)

where, for ease of legibility, we have omitted the dependence on |$\displaystyle \left( k_x , \omega \right)$|⁠, as well as the dependence on the depth level |$\displaystyle z_n$|⁠. For similar reasons, we also introduce the following shorthand notation for eqs (22) and (23), viz

(24)
(25)

which we will use throughout the rest of this paper.

In eqs (22) and (23) the terms |$\displaystyle R_{\ldots } ^ \cup$| and |$\displaystyle T_{\ldots } ^ +$| represent the reflection and transmission coefficients for waves coming from above, respectively, while the terms |$\displaystyle R_{\ldots } ^ \cap$| and |$\displaystyle T_{\ldots }^ -$| represent the reflection and transmission coefficients for waves coming from below. In principle, it is possible to use the true Zoeppritz reflection and transmission coefficients for these terms. We will, however, use the extended Shuey approximations detailed in Section 2.1 instead, to simplify the inverse problem described in Section 2.3. Applying the extended Shuey approximations to |$\displaystyle {R_{PP}^ \cup }$|⁠, for example, we write

(26)

where we have used the relationship

(27)

In a similar way, we can write the approximations to the other reflection and transmission coefficients.

Next, we write the relationship between wavefields at different depth levels, viz.

(28)
(29)

where we have introduced the propagation operators |$\displaystyle W_{P}$| and |$\displaystyle W_{S}$| for |$\displaystyle P$| and |$\displaystyle S$| wavefields, respectively. These operators are given by

(30)

where |$\displaystyle \Delta z$| is the distance between neighbouring depth levels. Note that the quantity |$\displaystyle \sqrt{{{\left( {\omega /\alpha \left( {z_n^ - } \right)} \right)}^2} - k_x^2}$| may become imaginary for certain values of |$\displaystyle \omega$|⁠, |$\displaystyle \alpha \left( {z_n^ - } \right)$| and |$\displaystyle k_x$|⁠. To avoid problems during modelling and inversion, we make the following substitution

(31)

in practice. Finally, we once again introduce shorthand notations for eqs (28) and (29), viz.

(32)
(33)

With the building blocks of eqs (24), (25), (32) and (33) in place, we now examine the forward modelling algorithm. To initialize the algorithm, we set all upgoing and downgoing wavefields to zero, that is, |$\displaystyle {{\mathbf {p}}^{ - ,0}}\left( {{z_n}} \right) = {{\mathbf {p}}^{ + ,0}}\left( {{z_n}} \right) = {\mathbf {0}}$|⁠. Note the notation of |$\displaystyle {{\mathbf {p}}^{ - ,0}}$| and |$\displaystyle {{\mathbf {p}}^{ + ,0}}$|⁠, where we have introduced an extra number in the superscript, which denotes how many so-called ‘round-trips’ have been modelled. Each round-trip increases the maximum order of multiples which are taken into account by one, up to the chosen number of round-trips to be modelled. Furthermore, for simplicity’s sake, we assume that there are no sources within the subsurface, only at the surface. In that case, we set |$\displaystyle {{\mathbf {p}}^{ + ,m}}\left( {{z_0}} \right) = {{\mathbf {s}}_0}$| for all |$\displaystyle m$| round-trips, where |$\displaystyle {{\mathbf {s}}_0}$| is a vector containing the source wavefield for one shot at |$\displaystyle z = z_0$|⁠. Finally, we assume that there are no upgoing waves coming from below the deepest depth level |$\displaystyle z_n = z_{N_z}$|⁠. Therefore, we write |$\displaystyle {{\mathbf {p}}^{ - ,m}}\left( {{z_{{N_z}}}} \right) = {\mathbf {0}}$| for all |$\displaystyle m$| round-trips.

We then begin by computing the downgoing wavefields for the first round-trip. Starting at |$\displaystyle z_n = z_0$|⁠, we use eq. (24) to write |$\displaystyle {{\mathbf {q}}^{ + ,1}}\left( {{z_0}} \right) = {{\mathbf {T}}^ + }\left( {{z_0}} \right){{\mathbf {p}}^{ + ,1}}\left( {{z_0}} \right) + {{\mathbf {R}}^ \cap }\left( {{z_0}} \right){{\mathbf {p}}^{ - ,0}}\left( {{z_0}} \right)$|⁠. Next, we apply the propagation operators, using eq. (32) to write |$\displaystyle {{\mathbf {p}}^{ + ,1}}\left( {{z_1}} \right) = {\mathbf {W}}\left( {{z_1},{z_0}} \right){{\mathbf {q}}^{ + ,1}}\left( {{z_0}} \right)$|⁠. At the depth level |$\displaystyle z_n = z_1$| we simply repeat this process with the appropriate transmission, reflection and propagation operators. In this way, we model the downgoing wavefield |$\displaystyle {{\mathbf {p}}^{ + ,1}}\left( {{z_n}} \right)$| at all depth levels.

Next, we compute the upgoing wavefield |$\displaystyle {{\mathbf {p}}^{ - ,1}}\left( {{z_n}} \right)$| at all depth levels. We start at the deepest depth level |$\displaystyle z_n = z_{N_z}$|⁠. We then use eq. (25) to write |$\displaystyle {{\mathbf {q}}^{ - ,1}}\left( {{z_{{N_z}}}} \right) = {{\mathbf {T}}^ - }\left( {{z_{{N_z}}}} \right){{\mathbf {p}}^{ - ,1}}\left( {{z_{{N_z}}}} \right) + {{\mathbf {R}}^ \cup }\left( {{z_{{N_z}}}} \right){{\mathbf {p}}^{ + ,1}}\left( {{z_{{N_z}}}} \right)$| and apply the propagation operators of eq. (33) to write |$\displaystyle {{\mathbf {p}}^{ - ,1}}\left( {{z_{{N_z} - 1}}} \right) = {\mathbf {W}}\left( {{z_{{N_z} - 1}},{z_{{N_z}}}} \right){{\mathbf {q}}^{ - ,1}}\left( {{z_{{N_z}}}} \right)$|⁠. Once again, we continue to apply these operators to find the upgoing wavefield |$\displaystyle {{\mathbf {p}}^{ - ,1}}\left( {{z_n}} \right)$| at each depth level.

We now repeat this process as many times as we wish to account for higher order scattering, where the scattering order is increased by one after each round-trip. The full process for finding the wavefields |$\displaystyle {{\mathbf {p}}^{ + ,M}}\left( {{z_n}} \right)$| and |$\displaystyle {{\mathbf {p}}^{ - ,M}}\left( {{z_n}} \right)$| for up to |$\displaystyle M$| scattering orders is illustrated in algorithm 1.

2.3 Elastic full wavefield migration

Next, we examine the elastic full wavefield migration (FWM) algorithm (Berkhout 2014b), which is the inversion process associated with the elastic FWMod algorithm described in Section 2.2. As in Section 2.2, we limit our analysis to the 1.5-D case. Note the acronyms used here, FWMod refers to the forward modelling algorithm, while FWM refers to the associated inversion process.

Consider a situation with known, |$\displaystyle P$|- and |$\displaystyle S$|-wave measurement data |$\displaystyle {{\bf d}} \left( k_x, \omega \right)$| at the surface. We also assume the source wavefield |$\displaystyle {{{\bf s}}_0}$| and the propagation operators |$\displaystyle {{\bf W}}\left( {{z_{n + 1}},{z_n}} \right)$| and |$\displaystyle {{\bf W}}\left( {{z_{n - 1}},{z_n}} \right)$| to be known. Our goal is to recover the contrast parameters |$\displaystyle {c_\alpha }\left( {{z_n}} \right)$|⁠, |$\displaystyle {c_\beta }\left( {{z_n}} \right)$| and |$\displaystyle {c_\rho }\left( {{z_n}} \right)$| at all depth levels |$\displaystyle z_n$|⁠.

To achieve this, we first define a cost function |$\displaystyle J$|⁠, which measures the mismatch between the measured data |$\displaystyle {{\bf d}} \left( k_x, \omega \right)$| and the forward modelled data at the surface after |$\displaystyle M$| roundtrips, |$\displaystyle {{{\bf p}}^{ - ,M}}\left( {{k_x , \omega , z_0}} \right)$|⁠, viz.

(34)

To retrieve the contrast parameters, we now apply a gradient-descent algorithm with respect to these parameters to minimize |$\displaystyle J$|⁠. Taking the gradient of the cost function with respect to the contrast parameter |$\displaystyle {c_\alpha }$| at depth level |$\displaystyle z_n$|⁠, for example, yields

(35)

where |$\displaystyle {{\bf e}}\left(k_x , \omega \right) = {{{\bf d}}\left( {{k_x},\omega } \right) - {{{\bf p}}^{ - ,M}}\left( {{k_x},\omega ,{z_0}} \right)}$|⁠.

We now examine the derivative of |$\displaystyle {{{\bf p}}^{ - ,M}}\left( {{k_x , \omega , z_0}} \right)$| with respect to |$\displaystyle {c_\alpha }\left( {{z_n}} \right)$|⁠. To evaluate this, we first examine the contribution of a scatterer located at the depth level |$\displaystyle z = z_n$| to the forward modelled wavefield |$\displaystyle {{{\bf p}}^{ - ,M}}\left( {{k_x , \omega , z_0}} \right)$|⁠, viz

(36)

where we have assumed |$\displaystyle {{{\bf p}}^{ - ,M-1}}\left( {{z_n}} \right) \approx {{{\bf p}}^{ - ,M}}\left( {{z_n}} \right)$|⁠. In eq. (36) we have also introduced the generalized propagation operators |$\displaystyle {{{{\bf \bar{W}}}}^ - }\left( {z_0 , z_n} \right)$| and |$\displaystyle {{{{\bf \bar{W}}}}^ \cup }\left( {z_0 , z_n} \right)$|⁠. These operators are constructed by applying sequences of propagation, reflection and transmission operators and are defined as

(37)
(38)
(39)

If we assume that the contrast parameters are independent at each depth level, we can write

(40)

with

(41)

where |$\displaystyle {\partial {{\bf R}}_{{c_\alpha }\left( {{z_n}} \right)}^ \cup }$|⁠, |$\displaystyle {\partial {{\bf R}}_{{c_\alpha }\left( {{z_n}} \right)}^ \cap }$|⁠, |$\displaystyle {\partial {{\bf T}}_{{c_\alpha }\left( {{z_n}} \right)}^ + }$| and |$\displaystyle {\partial {{\bf T}}_{{c_\alpha }\left( {{z_n}} \right)}^ - }$| are 2|$\displaystyle \times$|2 matrices. These matrices are constructed by calculating the derivatives of the approximate reflection and transmission coefficients introduced in Section 2.1. For example, the |$\displaystyle PP$| component of the matrix |$\displaystyle {\partial {{\bf R}}_{{c_\alpha }\left( {{z_n}} \right)}^ \cup }$| is given by

(42)

with similar definitions for the other components. Using eqs (41) and (42) we can now evaluate eq. (35) for all depth levels |$\displaystyle z_n$|⁠.

In a similar fashion, we can calculate the gradients with respect to |$\displaystyle {c_\beta }\left( {{z_n}} \right)$| and |$\displaystyle {c_\rho }\left( {{z_n}} \right)$|⁠. Using the full gradient, we calculate the linearized wavefield perturbation at the surface, viz.

(43)

Finally, using eq. (43), we calculate the update to the contrast parameters using

(44)

for |$\displaystyle {c \in {c_\alpha },{c_\beta },{c_\rho }}$|⁠, where the (real) constant |$\displaystyle \gamma$| is given by

(45)

Repeating the process described above for as many iterations as desired, one can retrieve the contrast parameters |$\displaystyle {c_\alpha }\left( {{z_n}} \right)$|⁠, |$\displaystyle {c_\beta }\left( {{z_n}} \right)$| and |$\displaystyle {c_\rho }\left( {{z_n}} \right)$| at all depth levels |$\displaystyle z_n$|⁠. This process is summarized in algorithm 2

3 NUMERICAL RESULTS

This section is, once again, split into three parts. First, in Section 3.1, we examine how many terms of the extended Shuey’s approximation introduced in Section 2.1 are required to accurately approximate the Zoeppritz equations. Next, in Section 3.2, we compare the forward modelled wavefields using the FWMod algorithm of 2.2 to true, synthetic, data generated using an elastic Kennett algorithm. Finally, in Section 3.3, we examine the inversion results of the elastic FWM algorithm of Section 2.3.

3.1 Extended Shuey’s approximations

We begin by examining the quality of the approximation of eq. (26) of Section 2.1. Specifically, we are interested in how many terms |$\displaystyle n \le N$| and |$\displaystyle \lambda \le \Lambda$| of the Taylor expansion are required for a good approximation. The Taylor expansions of |$\displaystyle R^\cup _{PP}$| and |$\displaystyle R^\cup _{SP}$| for different values of |$\displaystyle N$| and |$\displaystyle \Lambda$| for both a low- and a high-contrast interface are shown in Fig. 2. The associated medium parameters used for the low- and high-contrast interface are displayed in Table 1.

A comparison of the extended Shuey approximation to the exact Zoeppritz equations. On the left-hand side, results for a low-contrast interface are shown, while the right-hand side shows the results for a high-contrast interface. The medium parameters used for both interfaces are displayed in Table 1. Within each figure, the number of contrast terms $\displaystyle \Lambda$ taken into account is constant, while the different coloured lines represent the results for different values of $\displaystyle N$. These results are compared to the true, Zoeppritz, equations, which are represented by the solid black lines. Also note the dashed vertical lines, which indicate the maximum angle $\displaystyle \theta _{max} = 0.8 \cdot \theta _c$ up to which the approximation holds.
Figure 2.

A comparison of the extended Shuey approximation to the exact Zoeppritz equations. On the left-hand side, results for a low-contrast interface are shown, while the right-hand side shows the results for a high-contrast interface. The medium parameters used for both interfaces are displayed in Table 1. Within each figure, the number of contrast terms |$\displaystyle \Lambda$| taken into account is constant, while the different coloured lines represent the results for different values of |$\displaystyle N$|⁠. These results are compared to the true, Zoeppritz, equations, which are represented by the solid black lines. Also note the dashed vertical lines, which indicate the maximum angle |$\displaystyle \theta _{max} = 0.8 \cdot \theta _c$| up to which the approximation holds.

Table 1.

Parameters Fig. 2

Low contrast High contrast 
|$\displaystyle \alpha _1 =$| 1500 m s−1|$\displaystyle \alpha _2 =$| 1800 m s−1|$\displaystyle \alpha _1 =$| 1500 m s−1|$\displaystyle \alpha _2 =$| 3000 m s−1
|$\displaystyle \beta _1 =$| 750 m s−1|$\displaystyle \beta _2 =$| 1000 m s−1|$\displaystyle \beta _1 =$| 750 m s−1|$\displaystyle \beta _2 =$| 1600 m s−1
|$\displaystyle \rho _1 =$| 1000 m s−1|$\displaystyle \rho _2 =$| 1400 m s−1|$\displaystyle \rho _1 =$| 1000 m s−1|$\displaystyle \rho _2 =$| 1700 m s−1
Low contrast High contrast 
|$\displaystyle \alpha _1 =$| 1500 m s−1|$\displaystyle \alpha _2 =$| 1800 m s−1|$\displaystyle \alpha _1 =$| 1500 m s−1|$\displaystyle \alpha _2 =$| 3000 m s−1
|$\displaystyle \beta _1 =$| 750 m s−1|$\displaystyle \beta _2 =$| 1000 m s−1|$\displaystyle \beta _1 =$| 750 m s−1|$\displaystyle \beta _2 =$| 1600 m s−1
|$\displaystyle \rho _1 =$| 1000 m s−1|$\displaystyle \rho _2 =$| 1400 m s−1|$\displaystyle \rho _1 =$| 1000 m s−1|$\displaystyle \rho _2 =$| 1700 m s−1
Table 1.

Parameters Fig. 2

Low contrast High contrast 
|$\displaystyle \alpha _1 =$| 1500 m s−1|$\displaystyle \alpha _2 =$| 1800 m s−1|$\displaystyle \alpha _1 =$| 1500 m s−1|$\displaystyle \alpha _2 =$| 3000 m s−1
|$\displaystyle \beta _1 =$| 750 m s−1|$\displaystyle \beta _2 =$| 1000 m s−1|$\displaystyle \beta _1 =$| 750 m s−1|$\displaystyle \beta _2 =$| 1600 m s−1
|$\displaystyle \rho _1 =$| 1000 m s−1|$\displaystyle \rho _2 =$| 1400 m s−1|$\displaystyle \rho _1 =$| 1000 m s−1|$\displaystyle \rho _2 =$| 1700 m s−1
Low contrast High contrast 
|$\displaystyle \alpha _1 =$| 1500 m s−1|$\displaystyle \alpha _2 =$| 1800 m s−1|$\displaystyle \alpha _1 =$| 1500 m s−1|$\displaystyle \alpha _2 =$| 3000 m s−1
|$\displaystyle \beta _1 =$| 750 m s−1|$\displaystyle \beta _2 =$| 1000 m s−1|$\displaystyle \beta _1 =$| 750 m s−1|$\displaystyle \beta _2 =$| 1600 m s−1
|$\displaystyle \rho _1 =$| 1000 m s−1|$\displaystyle \rho _2 =$| 1400 m s−1|$\displaystyle \rho _1 =$| 1000 m s−1|$\displaystyle \rho _2 =$| 1700 m s−1

From the results of Fig. 2 we can draw a number of conclusions. First of all, we see that increasing |$\displaystyle N$| and |$\displaystyle \Lambda$| improves the approximation up to the critical angle, as we would expect. If the contrasts at the interface are low, we see that the quality of the approximation is mostly determined by the number of |$\displaystyle \sin ^n \left( \theta \right)$| terms |$\displaystyle N$| taken into account. By contrast, if the contrasts are high, we see that the quality is mostly determined by the number of contrast terms |$\displaystyle \Lambda$| taken into account.

In general, we conclude that, in order to obtain a high-quality approximation, one should choose a value of |$\displaystyle N > 2$|⁠, as well as a value of |$\displaystyle \Lambda > 1$|⁠, at the very least, especially for interfaces with high contrasts. In this case, the extended Shuey’s approximation yields reasonable results for both |$\displaystyle R^\cup _{PP}$| and |$\displaystyle R^\cup _{SP}$|⁠. Note that, while Fig. 2 only shows the results for |$\displaystyle R^\cup _{PP}$| and |$\displaystyle R^\cup _{SP}$|⁠, a similar analysis has been performed for the other reflection and transmission coefficients. These results can be found in the appendix.

3.2 Elastic FWMod

Next, we examine the results of the elastic FWMod algorithm, as described in Section 2.2. To benchmark the method, we compare the results to synthetic data generated using an elastic Kennett algorithm (Kennett 1984), which takes the full, Zoeppritz, reflection and transmission coefficients into account. The medium parameters used are equivalent to the high-contrast scenario described in Table 1. Note that, as we know that the approximation presented in Section 2.1 only holds up to the critical angle |$\displaystyle \theta _c$|⁠, data corresponding to angles |$\displaystyle \theta > 0.8 \cdot \theta _c$| have been removed from the data.

In Figs 3 and 4 we have displayed the difference between the forward modelled data and the true, synthetic, data at the surface for |$\displaystyle P$| and |$\displaystyle S$| waves, respectively. From these figures, we see similar behaviour as in Fig. 2. The difference between the true data at the surface and the forward modelled data decreases as |$\displaystyle N$| and |$\displaystyle \Lambda$| increase. Once again, we see that high-quality results are only achieved for |$\displaystyle N > 2$| and |$\displaystyle \Lambda > 1$|⁠.

Comparison of the FWMod results to true, synthetic, data for $\displaystyle P$ waves. In the top-left corner the synthetic data, generated using an elastic Kennet algorithm, is displayed. The top-centre image shows the underlying model parameters used. In the top-right corner the cost function $\displaystyle J$ for the $\displaystyle P$ waves is displayed for different values of $\displaystyle N$ and $\displaystyle \Lambda$ after five iterations. The remaining figures show the data misfit between the synthetic data and the forward modelled data at the surface for different values of $\displaystyle N$ and $\displaystyle \Lambda$ after five iterations. Note that these figures have been clipped to a value of 10 per cent of the maximum value of the data.
Figure 3.

Comparison of the FWMod results to true, synthetic, data for |$\displaystyle P$| waves. In the top-left corner the synthetic data, generated using an elastic Kennet algorithm, is displayed. The top-centre image shows the underlying model parameters used. In the top-right corner the cost function |$\displaystyle J$| for the |$\displaystyle P$| waves is displayed for different values of |$\displaystyle N$| and |$\displaystyle \Lambda$| after five iterations. The remaining figures show the data misfit between the synthetic data and the forward modelled data at the surface for different values of |$\displaystyle N$| and |$\displaystyle \Lambda$| after five iterations. Note that these figures have been clipped to a value of 10 per cent of the maximum value of the data.

Comparison of the FWMod results to true, synthetic, data for $\displaystyle S$ waves. In the top-left corner the synthetic data, generated using an elastic Kennet algorithm, is displayed. The top-centre image shows the underlying model parameters used. In the top-right corner the cost function $\displaystyle J$ for the $\displaystyle S$ waves is displayed for different values of $\displaystyle N$ and $\displaystyle \Lambda$ after five iterations. The remaining figures show the data misfit between the synthetic data and the forward modelled data at the surface for different values of $\displaystyle N$ and $\displaystyle \Lambda$ after five iterations. Note that these figures have been clipped to a value of 10 per cent of the maximum value of the data.
Figure 4.

Comparison of the FWMod results to true, synthetic, data for |$\displaystyle S$| waves. In the top-left corner the synthetic data, generated using an elastic Kennet algorithm, is displayed. The top-centre image shows the underlying model parameters used. In the top-right corner the cost function |$\displaystyle J$| for the |$\displaystyle S$| waves is displayed for different values of |$\displaystyle N$| and |$\displaystyle \Lambda$| after five iterations. The remaining figures show the data misfit between the synthetic data and the forward modelled data at the surface for different values of |$\displaystyle N$| and |$\displaystyle \Lambda$| after five iterations. Note that these figures have been clipped to a value of 10 per cent of the maximum value of the data.

This conclusion is supported by the plots of the cost function |$\displaystyle J$|⁠, seen in the top-right corner of Figs 3 and 4. In these plots we see that the cost function sharply decreases up to |$\displaystyle N = 4$|⁠, after which it remains more or less constant. Also, we see a considerable improvement when comparing the curves for |$\displaystyle \Lambda = 1$| (blue lines) to the curves for |$\displaystyle \Lambda = 2$| (red lines).

3.3 Elastic FWM

Finally, we examine the results of the FWM algorithm using the extended Shuey approximations, as described in Section 2.3. We apply the FWM algorithm to synthetic data, generated using the same elastic Kennett algorithm as we used in Section 3.2 for the high-contrast case of Table 1. The results of this process are shown in Fig. 5.

Results of the FWM method after 30 iterations, applied to synthetic data generated using an elastic Kennett algorithm. The underlying medium parameters are displayed in Table 1. On the left-hand side, the results for $\displaystyle c_\alpha$, $\displaystyle c_\beta$ and $\displaystyle c_\rho$ are displayed. The band-limited, ground truth contrasts are displayed in black, with the results for different values of $\displaystyle N$ and $\displaystyle \Lambda$ displayed in different colours. In the top-right corner, the cost function $\displaystyle J$ after 30 iterations is displayed for different values of $\displaystyle N$ and $\displaystyle \Lambda$. Finally, the data residual for both P and S waves after 30 iterations is displayed in the centre-right and bottom-right figures, respectively, for $\displaystyle N = 4$ and $\displaystyle \Lambda = 1$.
Figure 5.

Results of the FWM method after 30 iterations, applied to synthetic data generated using an elastic Kennett algorithm. The underlying medium parameters are displayed in Table 1. On the left-hand side, the results for |$\displaystyle c_\alpha$|⁠, |$\displaystyle c_\beta$| and |$\displaystyle c_\rho$| are displayed. The band-limited, ground truth contrasts are displayed in black, with the results for different values of |$\displaystyle N$| and |$\displaystyle \Lambda$| displayed in different colours. In the top-right corner, the cost function |$\displaystyle J$| after 30 iterations is displayed for different values of |$\displaystyle N$| and |$\displaystyle \Lambda$|⁠. Finally, the data residual for both P and S waves after 30 iterations is displayed in the centre-right and bottom-right figures, respectively, for |$\displaystyle N = 4$| and |$\displaystyle \Lambda = 1$|⁠.

From this figure, we initially conclude that the FWM algorithm recovers the contrasts to a reasonable accuracy, irrespective of the value of |$\displaystyle N$| and |$\displaystyle \Lambda$| chosen. Looking closer, however, we note that increasing the value of |$\displaystyle N$| to at least |$\displaystyle N = 4$| improves the result somewhat. From the figures of the contrasts, we see a slightly improved recovery of the ground truth for |$\displaystyle N = 4$| compared to |$\displaystyle N = 2$|⁠. Also, examining the cost function (displayed in the top-right corner of Fig. 5), we see that the cost function still decreases between |$\displaystyle N = 2$| and |$\displaystyle N = 4$|⁠.

Also note the case |$\displaystyle N = 0, \Lambda = 1$|⁠, displayed as blue lines in the plots on the left-hand side of Fig. 5. This case is of special interest, as it represents the result of applying conventional, angle-independent, acoustic FWM to the data. Note that, in this case, the contrasts |$\displaystyle c_\alpha$| and |$\displaystyle c_\rho$| are impossible to separate. This can be seen in Fig. 5, where the results for |$\displaystyle c_\alpha$| and |$\displaystyle c_\rho$| are identical for the acoustic result. Also, we see that |$\displaystyle c_\beta$| cannot be recovered in this case. These effects may lead to artefacts, which can be seen in the area between the two reflectors, as well as the area below the lowest reflector. While some artefacts are also visible for the other cases, they are reduced compared to the acoustic result.

Finally, we note that, in this case, taking |$\displaystyle \Lambda > 1$| does not appear to improve the results significantly. This is most likely due to the fact that, due to the band-limitation of the source wavefield, it is impossible to recover the sharp interfaces present in the medium. Instead, a spread-out, band-limited approximation of the contrasts is recovered. As this reduces the value of the contrasts, higher order powers of the contrasts do not notably impact the result. In cases where this band-limitation is compensated for, however, taking |$\displaystyle \Lambda > 1$| will improve the results.

4 DISCUSSION

In this section, we compare the extended Shuey’s approximations derived in this paper to two existing sets of approximations in literature. The first set of approximations we consider are the two- and three-term Shuey approximations, given by

(46)
(47)

Secondly, we also consider the linearized approximations of Aki and Richards (Aki & Richards 2002), which are given by

(48)
(49)

where we have introduced the average |$\displaystyle P$|-wave angle |$\displaystyle i = \left( i_1 + i_2 \right) / 2$| and the average |$\displaystyle S$|-wave angle |$\displaystyle j = \left( j_1 + j_2 \right) / 2$|⁠, following the notation of Section 2.1.

The reflection coefficients for the approximations under consideration at different angles are shown in Fig. 6. From this figure, we can draw a number of conclusions. First of all, we note that for high contrasts, both the 2- and 3-term Shuey’s approximations, shown in blue, as well as the linearized Aki-Richards approximations, shown in red, fail for angles up to the critical angle. In this case, one requires nonlinear terms with respect to the contrasts to properly approximate the reflection coefficients, which are absent in these approximations. This can be seen most clearly in the bottom right of Fig. 6, where we see that one requires third-order terms with respect to the contrasts to accurately approximate the reflection coefficient up to the critical angle.

Comparison of the extended Shuey approximations with the standard two- and three-term Shuey approximations and the linearized Aki-Richards approximations. On the left-hand side the $\displaystyle PP$-reflection coefficients are shown for a low and high contrast case. The parameters used are given in Table 1. On the right-hand side, the $\displaystyle SP$-reflection coefficients are shown.
Figure 6.

Comparison of the extended Shuey approximations with the standard two- and three-term Shuey approximations and the linearized Aki-Richards approximations. On the left-hand side the |$\displaystyle PP$|-reflection coefficients are shown for a low and high contrast case. The parameters used are given in Table 1. On the right-hand side, the |$\displaystyle SP$|-reflection coefficients are shown.

Secondly, we see that the standard, 2-term Shuey approximation (solid blue line), as well as the equivalent approximation for |$\displaystyle R_{SP}^ \cup$| using |$\displaystyle N = 1, \Lambda = 1$|⁠, only holds for small contrasts and angles far from the critical angle. While this behaviour is expected, one should note the very small range of angles for which this approximation holds, especially for the |$\displaystyle R_{SP}^ \cup$| term. This limits the applicability of this approximation for forward and inverse modelling, where good performance of the approximation over a range of angles is required.

Finally, we note that the linearized Aki-Richards approximation (solid red line) performs very well for low-contrast scenarios, showing a very good approximation for all angles. This is due to the |$\displaystyle 1 / \cos ^2 \left( i \right)$| and |$\displaystyle 1 / \cos \left( j \right)$| terms in eqs (48) and (49). While these terms lead to a high-quality approximation, they also introduce critical points in the approximations, making them difficult to invert for directly. The extended Shuey’s approximations of this paper avoid this problem, while showing similar behaviour up to the cutoff angle |$\displaystyle \theta _c$| for |$\displaystyle N \ge 4$| and |$\displaystyle \Lambda \ge 2$|⁠. For high-contrast scenarios, the linearized Aki-Richards approximation deteriorates, due to the lack of higher order terms with respect to the contrasts. In these situations, the extended Shuey’s approximations yield a more accurate approximation when taking |$\displaystyle \Lambda \ge 2$|⁠.

5 CONCLUSION

In this paper, we present an alternative set of approximations to the Zoeppritz equations for the elastic reflection and transmission coefficients, which we call the extended Shuey approximations. We show that these approximations can be applied to the elastic FWMod forward modelling algorithm as well as the elastic FWM migration algorithm, where we have shown that using the extended Shuey approximations improves the forward and inverse modelling results compared to using the standard Shuey approximations. Finally, we have compared the extended Shuey approximations to the conventional approximations in literature, where we have shown that the extended Shuey approximations achieve comparable performance in cases where the contrasts are low, up to the critical angle, and achieve better results in cases with large contrasts. Based on these results, we conclude that the extended Shuey approximations presented in this paper are a useful addition to the existing approximations to the Zoeppritz equations, specifically in their application in forward and inverse modelling.

ACKNOWLEDGMENTS

The authors thank the sponsoring companies from the Delphi Consortium at Delft University of Technology for their support to this research.

DATA AVAILABILITY

The synthetic data used in this paper is available upon request. Due to agreements made with the sponsoring companies of the Delphi Consortium, we are unable to share the code used to generate the results of this paper.

REFERENCES

Aki
 
K.
,
Richards
 
P.G.
,
2002
.
Quantitative Seismology
, 2nd edn,
University Science Books
.

Berkhout
 
A.J.
,
2014a
.
Review paper: an outlook on the future of seismic imaging, Part I: forward and reverse modelling
,
Geophys. Prospect.
,
62
,
911
930
.

Berkhout
 
A.J.
,
2014b
.
Review paper: an outlook on the future of seismic imaging, Part II: full-wavefield migration
,
Geophys. Prospect.
,
62
,
931
949
.

Jones
 
I.F.
,
Davison
 
I.
,
2014
.
Seismic imaging in and around salt bodies
,
Interpretation
,
2
,
1N
T284
.

Kennett
 
B.L.N.
,
1984
.
An operator approach to forward modeling, data processing and migration
,
Geophys. Prospect.
,
32
,
1074
1090
.

Prieux
 
V.
,
Brossier
 
R.
,
Operto
 
S.
,
Virieux
 
J.
,
2013
.
Multiparameter full waveform inversion of multicomponent ocean-bottom-cable data from the Valhall field. Part 2: imaging compressive-wave and shear-wave velocities
,
Geophys. J. Int.
,
194
,
1665
1681
.

Stewart
 
R.R.
,
Gaiser
 
J.E.
,
Brown
 
R.J.
,
Lawton
 
D.C.
,
2003
.
Converted–wave seismic exploration: applications
,
Geophysics
,
68
,
40
57
.

Shuey
 
R.T.
,
1985
.
A simplification of the Zoeppritz equations
,
Geophysics
,
4
,
530
743
.

Wang
 
Y.
,
1999
.
Approximations to the zoeppritz equations and their use in AVO analysis
,
Geophysics
,
64
,
1673
1955
.

Zoeppritz
 
K.
,
1919
.
VIIb. Über reflexion und durchgang seismischer wellen durch unstetigkeitsflächen
,
Nachr Königl Gesell. Wissensch. Göttingen Math.-phys. Klasse
,
1919
,
66
84
.

APPENDIX A: ADDITIONAL REFLECTION AND TRANSMISSION COEFFICIENTS

In this section we show the approximations to the full set of reflection and transmission coefficients. Fig. A1 shows the four terms of |${\bf R}^\cup$|⁠, Fig. A2 shows the four terms of |${\bf R}^\cap$|⁠, Fig. A3 shows the four terms of |${\bf T}^+$| and Fig. A4 shows the four terms of |${\bf T}^-$|⁠.

A comparison of the extended Shuey approximation to the exact Zoeppritz equations for ${\bf R}^\cup$. The left-most column shows the 1,1-element, which corresponds to the $PP$ term. The centre-left column shows the 1,2-element, which corresponds to the $PS$ term. The centre-right column shows the 2,1-element, which corresponds to the $SP$ term. Finally, the right-most column shows the 2,2-element, which corresponds to the $SS$ term. The medium parameters used for the interface correspond to the low-contrast parameters displayed in Table 1. Within each figure, the number of contrast terms $\Lambda$ taken into account is constant, while the different coloured lines represent the results for different values of N. These results are compared to the true, Zoeppritz equations, which are represented by the solid black lines. Also note the dashed vertical lines, which indicate the maximum angle $\theta _{max} = 0.8 \cdot \theta _c$ up to which the approximation holds.
Figure A1.

A comparison of the extended Shuey approximation to the exact Zoeppritz equations for |${\bf R}^\cup$|⁠. The left-most column shows the 1,1-element, which corresponds to the |$PP$| term. The centre-left column shows the 1,2-element, which corresponds to the |$PS$| term. The centre-right column shows the 2,1-element, which corresponds to the |$SP$| term. Finally, the right-most column shows the 2,2-element, which corresponds to the |$SS$| term. The medium parameters used for the interface correspond to the low-contrast parameters displayed in Table 1. Within each figure, the number of contrast terms |$\Lambda$| taken into account is constant, while the different coloured lines represent the results for different values of N. These results are compared to the true, Zoeppritz equations, which are represented by the solid black lines. Also note the dashed vertical lines, which indicate the maximum angle |$\theta _{max} = 0.8 \cdot \theta _c$| up to which the approximation holds.

Same as Fig. A1, but for the ${\bf R}^\cap$ terms.
Figure A2.

Same as Fig. A1, but for the |${\bf R}^\cap$| terms.

Same as Fig. A1, but for the ${\bf T}^+$ terms.
Figure A3.

Same as Fig. A1, but for the |${\bf T}^+$| terms.

Same as Fig. A1, but for the ${\bf T}^-$ terms.
Figure A4.

Same as Fig. A1, but for the |${\bf T}^-$| terms.

This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.