- Articles & preprints
- Submission
- Policies
- Peer review
- Editorial board
- About
- EGU publications
- Manuscript tracking
Research article 12 Mar 2018
Research article | 12 Mar 2018
Correspondence: Arjun Datta (arjun.datta@tifr.res.in)
HideCorrespondence: Arjun Datta (arjun.datta@tifr.res.in)
We present a suite of programs that implement decades-old algorithms for computation of seismic surface wave reflection and transmission coefficients at a welded contact between two laterally homogeneous quarter-spaces. For Love as well as Rayleigh waves, the algorithms are shown to be capable of modelling multiple mode conversions at a lateral discontinuity, which was not shown in the original publications or in the subsequent literature. Only normal incidence at a lateral boundary is considered so there is no Love–Rayleigh coupling, but incidence of any mode and coupling to any (other) mode can be handled. The code is written in Python and makes use of SciPy's Simpson's rule integrator and NumPy's linear algebra solver for its core functionality. Transmission-side results from this code are found to be in good agreement with those from finite-difference simulations. In today's research environment of extensive computing power, the coded algorithms are arguably redundant but SWRT can be used as a valuable testing tool for the ever evolving numerical solvers of seismic wave propagation. SWRT is available via GitHub (https://github.com/arjundatta23/SWRT.git).
There is a vast body of seismological literature on analytical methods for describing surface wave propagation in laterally inhomogeneous media where the lateral heterogeneity takes simple forms such as inclusions of regular shape or plane interfaces. In particular, models with plane vertical interfaces where the elastic parameters have a first-order discontinuity received considerable attention in the 1960s and 1970s as idealizations of sharp, large-scale heterogeneities that actually exist in the Earth such as continental margins and grabens. The convenience of a vertical boundary is that the wavefield can be expressed in terms of incident, reflected and transmitted waves – this not only yields physical insight into the problem but also reduces it to the evaluation of reflection and transmission coefficients at the lateral boundary. In general the reflection and transmission of surface waves at a plane vertical interface cannot be solved for exactly because the boundary conditions for stresses and displacements cannot be satisfied exactly by a field of surface waves alone, which in turn is because the surface wave eigenfunctions for a (layered) half-space do not form a complete set of basis functions for the elastic wave equation. The complete basis set consists of a discrete number of surface wave modes plus a continuum of homogeneous body waves (see e.g. Levshin, 1989; Alsop et al., 1974). A range of approximate methods were therefore developed by theoretical seismologists to study the interaction of surface waves with geological discontinuities and these have been comprehensively reviewed (Malischewsky, 1987; Keilis-Borok, 1989). In this paper we revisit some of the original algorithms with the aim of demonstrating their hitherto undemonstrated ability to model mode conversions. In the process we also present them under the unified framework of the Herrera (1964) scalar product.
We focus on the body wave or ray-theoretical method of Gregersen and Alsop (1974) for Love waves and the Green's function method of Its and Yanovskaya (1985) for both Love and Rayleigh waves. For comparison we also show results obtained with the seminal mode-matching method of Alsop (1966). Although the methods highlighted in this study have since been advanced and improved upon (Vaccari et al., 1989; Its, 1991; Its and Lee, 1993; Romanelli et al., 1996), their ability to model mode conversions – a critical aspect of surface wave propagation through laterally homogeneous media – even in the simplest of scenarios has never been shown.
Since SWRT is based on previously published techniques, we only discuss the implementation of those techniques in the SWRT package. Throughout this section, we use a_{r} and b_{t} to denote reflection and transmission coefficients respectively at any given frequency. Mode indices r and t span all the modes that exist, at that frequency, on the incidence and transmission sides of the vertical discontinuity. A_{r} and B_{t} denote the corresponding amplitudes of modes (with respect to normalized eigenfunctions) in the two media. In this paper we focus on transmission results (rather than reflection) because the transmitted wavefield is easily measured from the frequency-domain finite-difference solver which we use for testing. Accordingly, the end result of each algorithm is presented in terms of a “transmission surface ratio” (Alsop, 1966, hereafter Alsop66). This quantity – shorthand for surface displacement ratio of transmitted to incident waves – provides a physically meaningful representation, independent of normalization conventions, of surface wave propagation through a discontinuity.
On the issue of normalization, the main difference between the SWRT implementation and the original literature is in the normalization of eigenfunctions used. The original publications used eigenfunctions ϕ normalized by the Herrera (1964) scalar product, whereas in this study we have worked with eigenfunctions ϕ^{u}, which are normalized to unit displacement at the surface; this leads to the appearance of normalization factors in the equations defining the algorithms.
As the Herrera normalization and the consequent relationship between ϕ and ϕ_{u} underscore the entire SWRT implementation, these are first discussed before getting into the details of individual algorithms. Herrera (1964) provided an orthogonality relation for surface wave eigenfunctions which, for plane waves in 2-D and suitably chosen coordinate axes, can be written as
where m and n are mode indices, and ^{*} denotes the complex conjugate. u_{q} and σ_{xq} are components of the displacement and traction vectors respectively on a surface normal to the propagation direction, x.
For the types of problems considered in this study, with the lateral discontinuity in the yz plane and Love waves propagating along the x direction, ${u}_{x}={u}_{z}=\mathrm{0}$, ${\mathit{\sigma}}_{xx}={\mathit{\sigma}}_{xz}=\mathrm{0}$ and σ_{xy}=ikμu_{y}. The displacement u_{y} is proportional to the eigenfunction Φ, so in this case Eq. (1) yields the simple norm
With the exception of possible scalar factors, the original authors of all the algorithms considered in this study used Love wave eigenfunctions normalized as
which implies H_{m}=i. SWRT on the other hand uses Φ=ϕ^{u}, and
From Eqs. (3) and (4) it is clear that ϕ and ϕ_{u} are related as
With this introduction we go on to discussing the individual methods. The Rayleigh equivalent of the development following Eq. (1) is deferred to Appendix Appendix A since it is only required for one of the methods.
Of the three methods included in this study, this is the one where mode conversions have actually been demonstrated in the published literature. Nonetheless, the method, hereafter referred to simply as the Alsop method, is included here for its remarkable simplicity.
a_{r} and b_{t} are obtained as solutions to the matrix equations
where
and [X]_{i} denotes the ith column of any matrix X. The matrices P, S, T and V in the above equations are composed of the following integrals:
To compute the transmission surface ratio (Y), we note that, for incident mode s and transmitted mode t, the ratio of transmitted to incident displacement is
In this study ϕ^{u(1)} and ϕ^{u(2)} have unit displacement at the surface; therefore
This method is founded on representing the field of surface waves as a superposition of plane homogeneous and inhomogeneous body waves. Body wave reflection and transmission coefficients are computed over the vertical interface to obtain outgoing stresses and displacements on the interface. This outgoing stress-displacement system is regarded as the source of outgoing surface waves, whose amplitudes are computed via its projection onto the outgoing surface wave eigenfunctions. In this study the method, hereafter referred to as the GAl method, is implemented, as follows for Love waves:
SH wave reflection and transmission coefficients over the vertical contact are calculated using the phase velocity of medium 1.
Elementary formulae are used for the SH reflection and transmission coefficients (Γ_{R} and Γ_{T}). The textbook expressions (e.g. Lay and Wallace, 1995, p. 102)
for plane SH rays at an interface between two solid media (the symbols have their usual meanings) are rewritten in terms of phase velocity because cos θ_{i} and cos θ_{t} can be rewritten (by Snell's law) in terms of the phase velocity c_{1} of the incidence-side medium (Keilis-Borok, 1989):
These equations are applied individually to each section of the vertical interface, i.e. each layered contact (there is really a layer index in addition to the medium index in Eqs. 10 and 11 but it has been left out for notational convenience) to obtain a vector of body wave coefficients over the entire vertical interface.
Outgoing stresses and displacements as a function of depth over the vertical interface are calculated using the eigenfunction of the incident Love wave mode and the SH reflection/transmission coefficients.
Coupling between the interface stress-displacement system and the eigenfunctions of outgoing Love waves (eigenfunctions of medium 1 for reflection, those of medium 2 for transmission) is calculated by means of the Herrera (1964) scalar product.
Since Herrera's scalar product may be computed for any two stress-displacement vectors S^{m} and S^{n}, this method uses Love–Love and SH–Love systems. In the former case Eq. (1) defines the orthonormalization of Love wave eigenfunctions; in the latter it yields the coupling coefficient for the GAl method. This means that when the surface wave modes are normalized by Eq. (1), the coupling computed (by the left-hand side of the same) between any stress-displacement system and a surface wave mode is the amplitude generated of that normalized surface wave mode.
So if S_{R} and S_{T} are the SH stress-displacement vectors on the interface computed in step 2 of the algorithm, Love wave reflection and transmission coefficients are given by (S_{R},ϕ^{(1)}) and (S_{T},ϕ^{(2)}). Therefore in this study they are computed as
In the above expressions the norms N have the usual definition given by Eq. (4), whilst the numerators are similarly evaluated using μ and k values corresponding to the incidence-side medium for a_{r} and to the transmission-side medium for b_{t}. Finally, as with the Alsop method, transmission surface ratio is given by
This method, hereafter referred to as the GF method, has been implemented in SWRT for both Love and Rayleigh waves. We follow the algorithm of Its and Yanovskaya (1985, hereafter ItsYan85), who obtained the coupled system of equations:
where P_{rt} and S_{st} are the coupling integrals
where u and σ are the displacement and traction vectors on the vertical plane of discontinuity, and J_{m}=2H_{m}. Using the appropriate u and σ in Eq. (15) leads to explicit expressions for the Love and Rayleigh cases.
Love waves — these expressions are easily derived as there is only one component of displacement and traction relevant for Eq. (15). The norm J_{m} reduces to 2i×N_{m} and for eigenfunctions normalized by Eq. (3), ${J}_{\mathrm{r}}^{\left(\mathrm{1}\right)}={J}_{\mathrm{t}}^{\left(\mathrm{2}\right)}=\mathrm{2}i$, leading to factors of 1∕2 outside the integrals in Eq. (15). In this study, working with ϕ^{u}, P and S are computed as
Rayleigh waves — the expressions have been derived by Keilis-Borok (1989), using a nomenclature similar to that of Aki and Richards (2002) for the stress-displacement vector, although with a sign error in one of the terms. The expressions used in this study are different, because they are based on the stress-displacement vector as defined by Gomberg and Masters (1988). The stress-displacement vector is ${\mathit{\beta}}_{\mathrm{1}},{\mathit{\beta}}_{\mathrm{2}},{\mathit{\beta}}_{\mathrm{3}},{\mathit{\beta}}_{\mathrm{4}}$, and in terms of these quantities the coupling integrals are (Appendix Appendix A)
with the norm given by
and ψ_{xx} defined by Eq. (A2). Note that in the Rayleigh case, the eigenfunctions themselves have subscripts 1–4; these should not be confused with the medium indices 1–2 which are shown as superscripts and in parentheses.
In our implementation P_{rt} and S_{st} are both real-valued, so the system of Eq. (14) is solved as
where R and T are the highest mode numbers that exist in the incidence-side and transmission-side media respectively at the given frequency. Equation (14) is easily recognized as a matrix equation of the form Cx=d:
which is a square system of equations of size $(T+R+\mathrm{2})$. In general R≠T and this can lead to cumbersome code because of the variability in the declaration of array sizes. However SWRT achieves efficiency with a single declaration regardless of the relationship between R and T. In the special case when $R=T\text{(}=M-\mathrm{1}\text{, say})$, the coefficient matrix of this system is itself composed of square sub-matrices:
Finally, we note that in computing the transmission surface ratio (Y), there is a subtle difference between the Love and Rayleigh cases. In the Love wave case, the appearance of the normalization factors for both media outside the brackets in Eq. (16) ensures that the a_{r} and b_{t} obtained from Eq. (19) are defined with respect to eigenfunctions normalized according to Eq. (3); therefore, as before, ${Y}_{\mathrm{st}}={b}_{\mathrm{t}}\times \sqrt{{N}_{\mathrm{s}}^{\left(\mathrm{1}\right)}/{N}_{\mathrm{t}}^{\left(\mathrm{2}\right)}}$. However the same is not true of the Rayleigh case – the implementation Eq. (17) ensures that the reflection and transmission coefficients obtained from Eq. (19) are in fact defined with respect to the eigenfunctions β, which have unit surface displacement. Hence, in the Rayleigh case, Y_{st}=b_{t}.
The three methods described in this section are each implemented as Python programs in the SWRT code. The only inputs these programs require are the eigenfunctions for the two media on either side of the lateral discontinuity and the depths of horizontal interfaces in these media. The programs have been tested using published results in the literature (Appendix Appendix B); in the following section we focus on results from the latter two methods, which are incomplete in the published literature.
We use a single model to demonstrate the capabilities of the GAl and GF methods. This is a model of an ocean–continent boundary known as “model L” (Fig. 1), and was used by both Gregersen and Alsop (1974, hereafter GregAl74) and ItsYan85. Both of these studies considered the propagation of a Love wave fundamental mode through this model, in both directions and over a range of frequencies in which a number of higher modes exist on the transmission side. However both studies have shown only the transmission of the fundamental mode, with no mention of any conversion to higher modes. Using the SWRT code, we find that considering periods as short as ∼ 7 s, there is significant conversion of the fundamental mode to up to the fourth overtone, for propagation in either direction (Fig. 2). Concurrently, the fundamental mode result obtained from SWRT is in good agreement with the published results (see Appendix Appendix B).
To validate our results of transmission to multiple higher modes, we compare the SWRT result with results from finite-difference (FD) simulations. The FD simulations are performed using the method of Roecker et al. (2010) duly adapted for teleseismic surface wave modelling. However for the FD comparison, model L is modified slightly in that the crustal thickness on the oceanic side is taken as 5 km rather than 7 km (see Fig. 1). In the rest of this section, all “model L” references imply this modified version.
Figure 3 shows excellent agreement between the SWRT results (GAl and GF methods) and FD results for Love waves in model L. Only forward propagation is considered. The corresponding Rayleigh case is shown in Fig. 4.
Figures 2–4 represent the core contribution of this study. The reader is reminded that model L has only been used as a demo example in this study. The published literature does not contain any examples of mode conversions modelled by the GAl method or the GF method, whilst the SWRT code provides complete solutions (with coupling to all possible modes) in any model containing a single vertical contact between two laterally homogeneous media. To emphasize this point as well as the general accuracy of the SWRT code, we provide a second demonstration using the model of Sato (1961). This is a unique model (Fig. 5) in which the Love wave reflection and transmission problem has an exact solution – the Love wave eigenfunctions on both sides of the discontinuity have the same shape; therefore, boundary conditions on either side of the discontinuity can be exactly satisfied by these eigenfunctions, therefore there is no mode conversion in this model, in either direction. Figure 6 demonstrates this with SWRT results.
This paper has described the author's SWRT code and used it to demonstrate cross-branch mode coupling, of both Love and Rayleigh waves, in simple 2-D models for which previously published results showed only self coupling. SWRT provides a consolidated implementation of sophisticated algorithms by previous authors, which is fully capable of modelling mode conversions arising from normal incidence at a single lateral discontinuity. These algorithms are more accurate than that of Alsop66, whose ability to model mode conversion has been known since the original author's work.
We conclude with remarks on the transferability and flexibility of the SWRT code. Whilst the code has been written to work with a specific format of input file (containing the surface wave eigenfunctions), it is modular with respect to file input: the eigenfunction-reading module of SWRT, which has nothing to do with the algorithms discussed in this paper, can simply be replaced in order to adapt the code to differently formatted input files. Possible additional work may be to modify the expressions for the coupling coefficients in the main programs, depending on the theoretical definitions of the eigenfunctions. As a further comment on flexibility, there is no restriction on the type of normalization that the input eigenfunctions must obey; only the expressions for surface ratios will need modification if the eigenfunctions are not normalized to unit surface displacement. Finally, whilst results have been presented in this paper only for the case of an incident fundamental mode, application of SWRT to the case of arbitrary mode number being incident can be trivially achieved.
SWRT is available at https://doi.org/10.5281/zenodo.1061008 (Datta, 2017).
The Rayleigh wave stress-displacement vector defined by Gomberg and Masters (1988) is $\mathit{\beta}=\mathit{\{}{\mathit{\beta}}_{\mathrm{1}},{\mathit{\beta}}_{\mathrm{2}},{\mathit{\beta}}_{\mathrm{3}},{\mathit{\beta}}_{\mathrm{4}}\mathit{\}}$, the four components being proportional to ${u}_{z},{u}_{x},{\mathit{\sigma}}_{zz}$ and σ_{xz} respectively:
For the type of problem under consideration, with the plane interface being normal to the x direction, the required components of traction are σ_{xz} and σ_{xx}. The latter can be written in terms of β_{1} and β_{2} (with the phase factor omitted for convenience):
The first requirement is to obtain an expression for the scalar product J_{mn}=2H_{mn}. Using Eq. (A2) and the first three parts of Eq. (A1),
So the norm required in Eq. (15) is
Similarly the expressions for the coupling integrals P_{rt} and S_{st} can be evaluated using Eq. (15) and the result is Eq. (17) of the main text.
The author declares that he has no conflict of interest.
I would like to acknowledge the Dr. Manmohan Singh Scholarship provided by
St. John's College in support of a PhD studentship at the University of
Cambridge as well as a Society of Exploration Geophysicists (SEG) Foundation
scholarship (SEG/Chevron Scholarship and SEG/John Bookout Scholarship)
received during the last year of my study. I would also like to thank Keith
Priestley for introducing me to the literature on which this work is
based.
Edited by: Luis Vazquez
Reviewed by: three anonymous referees
Aki, K. and Richards, P. G.: Quantitative Seismology, 2nd Edn., University Science Books, 2002. a
Alsop, L.: Transmission and reflection of Love waves at a vertical discontinuity, J. Geophys. Res., 71, 3969–3984, 1966. a, b
Alsop, L., Goodman, A., and Gregersen, S.: Reflection and transmission of inhomogeneous waves with particular application to Rayleigh waves, B. Seismol. Soc. Am., 64, 1635–1652, 1974. a
Datta, A.: SWRT, https://doi.org/10.5281/zenodo.1061008, 2017. a
Gomberg, J. S. and Masters, T. G.: Waveform modelling using locked-mode synthetic and differential seismograms: application to determination of the structure of Mexico, Geophys. J. Int., 94, 193–218, 1988. a, b
Gregersen, S. and Alsop, L.: Amplitudes of horizontally refracted Love waves, B. Seismol. Soc. Am., 64, 535–553, 1974. a, b
Herrera, I.: On a method to obtain a Green's function for a multi-layered half space, B. Seismol. Soc. Am., 54, 1087–1096, 1964. a, b, c, d
Its, E.: Surface wave propagation across a narrow vertical layer between horizontally homogeneous quarter-spaces, Wave Motion, 14, 11–23, 1991. a
Its, E. and Lee, J.: Reflection and transmission of surface waves at a vertical interface in anisotropic elastic media, B. Seismol. Soc. Am., 83, 1355–1372, 1993. a
Its, E. and Yanovskaya, T.: Propagation of surface waves in a half-space with vertical, inclined or curved interfaces, Wave Motion, 7, 79–94, 1985. a, b
Keilis-Borok, V.: Surface Waves in Media Involving Vertical Contacts, in: Seismic Surface Waves in a Laterally Inhomogeneous Earth, edited by: Keilis-Borok, V., vol. 9 of Modern Approaches in Geophysics, Springer Netherlands, 71–98, 1989. a, b, c
Lay, T. and Wallace, T. C.: Modern Global Seismology, vol. 58, Academic Press, 1995. a
Levshin, A.: Surface Waves in Vertically Inhomogeneous Media, in: Seismic Surface Waves in a Laterally Inhomogeneous Earth, edited by: Keilis-Borok, V., vol. 9 of Modern Approaches in Geophysics, Springer Netherlands, 3–33, 1989. a
Malischewsky, P.: Surface waves and discontinuities, Elsevier Science Pub. Co. Inc., New York, NY, 1987. a
Roecker, S., Baker, B., and McLaughlin, J.: A finite-difference algorithm for full waveform teleseismic tomography, Geophys. J. Int., 181, 1017–1040, 2010. a
Romanelli, F., Bing, Z., Vaccari, F., and Panza, G.: Analytical computation of reflection and transmission coupling coefficients for Love waves, Geophys. J. Int., 125, 132–138, 1996. a
Sato, R.: Love waves propagated across transitional zone, Japanese Journal of Geophysics, 2, 117–134, 1961. a, b
Vaccari, F., Gregersen, S., Furlan, M., and Panza, G.: Synthetic seismograms in laterally heterogeneous, anelastic media by modal summation of P–SV waves, Geophys. J. Int., 99, 285–295, 1989. a
An interactive open-access journal of the European Geosciences Union