SWRT : A package for semi-analytical solutions of surface wave propagation , including mode conversion , across transversely aligned vertical discontinuities

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 finitedifference 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).


Introduction
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.
Published by Copernicus Publications on behalf of the European Geosciences Union.
A. Datta: Solutions of surface wave propagation across vertical discontinuities 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.

Implementation of methods
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 = 0, σ xx = σ xz = 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 (2) 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 A since it is only required for one of the methods.
2.1 Mode-matching method of Alsop (1966) 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 Geosci.Instrum.Method.Data Syst., 7, 101-112, 2018 www.geosci-instrum-method-data-syst.net/7/101/2018/ 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 2.2 Body wave method of Gregersen and Alsop (1974) 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: 1. 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.
2. 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.
3. 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)   The SWRT package accounts for all the modes that exist in a given medium at a given frequency; note that in this figure the highest modes, where they exist, have negligible amplitudes.Here and in all following figures, the various surface wave modes are labelled according to their mode number, starting from 0 for the fundamental and increasing successively for each higher mode.
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 lefthand 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 2.3 Green's function method of Its and Yanovskaya (1985) 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.
Building the system of equations 1. 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 (1) r = J (2) t = 2i, leading to factors of 1/2 outside the integrals in Eq. ( 15).In this study, working with φ u , P  (c) using the Alsop method.Note that the FD results are the same in all three plots.The last plot (c) is shown to highlight the weaker performance of the Alsop method (at the low-frequency end for the first and second higher modes), compared to the two methods which are the focus of this paper.and S are computed as 2. Rayleigh waves -the expressions have been derived by Keilis-Borok (1989), using a nomenclature similar to that of Aki and Richards (2002) for the stressdisplacement 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 β 1 , β 2 , β 3 , β 4 , and in terms of these quantities the coupling integrals are (Appendix A) with the norm given by  Figure 4. Transmission surface ratio obtained by the GF method (solid lines) and FD calculations (dots), for an incident Rayleigh wave fundamental mode propagating through model L in the forward direction.Note that the conversion to higher modes is much weaker than in the Love wave case, but the SWRT result is in good agreement with the FD result.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.

Solving the system of equations
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: . Schematic representation of the model of Sato (1961).Units of layer parameters are the same as in Fig. 1.
which is a square system of equations of size (T + R + 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 ( = M − 1, 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, t .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 B); in the following section we focus on results from   (stars).Even though higher modes exist in the model at the frequencies considered and SWRT by default computes coupling to all modes that exist at any given frequency, all three methods exactly predict that there is no mode conversion.the latter two methods, which are incomplete in the published literature.

Demonstration of results
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 It-sYan85.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 B).
To validate our results of transmission to multiple higher modes, we compare the SWRT result with results from finitedifference (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.

Discussion and conclusions
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.

Figure 1 .
Figure1.Schematic representation of "Model L".Layer thicknesses are shown in km and units for velocity and density are km s −1 and gm cc −1 respectively.Note that crustal thickness on the right-hand side was originally 7 km but is modified to 5 km for comparisons with finite-difference results in this study.The layer thickness of 5 km facilitates the use of a 2.5 km FD grid, keeping the numerical problem computationally tractable (a 7 km thick layer would require a much finer FD grid spacing of 1 km).

Figure 2 .
Figure2.Transmission surface ratio obtained by the GAl method (solid lines) and GF method (dashed lines with circles), for an incident Love wave fundamental mode (mode 0) propagating through the original model L in the forward direction (a) and backward direction (b).The SWRT package accounts for all the modes that exist in a given medium at a given frequency; note that in this figure the highest modes, where they exist, have negligible amplitudes.Here and in all following figures, the various surface wave modes are labelled according to their mode number, starting from 0 for the fundamental and increasing successively for each higher mode.

Figure 3 .
Figure 3. Transmission surface ratios obtained from the SWRT code (solid lines) and FD calculations (dots), for an incident Love wave fundamental mode propagating through the modified model L in the forward direction.(a) Using the GAl method; (b) using the GF method;(c) using the Alsop method.Note that the FD results are the same in all three plots.The last plot (c) is shown to highlight the weaker performance of the Alsop method (at the low-frequency end for the first and second higher modes), compared to the two methods which are the focus of this paper.

Figure 6 .
Figure6.Transmission surface ratios for an incident Love wave fundamental mode propagating in the forward direction (a) and backward direction (b) in the Sato model.Shown are results from the GAl method (solid lines), the GF method (circles) and the method of Alsop (stars).Even though higher modes exist in the model at the frequencies considered and SWRT by default computes coupling to all modes that exist at any given frequency, all three methods exactly predict that there is no mode conversion.