Coupled-mode separation of multiply scattered wavefield components in two-dimensional waveguides

For a waveguide that is invariant in one of the horizontal directions, this paper presents a mathematically exact partial-wave decomposition of the wavefield for assessment of multiple scattering by horizontally displaced medium anomalies. The decomposition is based on discrete coupled-mode theory and combination of reflection/transmission matrices. In particular, there is no high-frequency ray approximation. Full details are presented for the scalar case with plane-wave incidence from below. An application of interest concerns ground motion induced by seismic waves, which may be severely amplified by local medium anomalies such as alluvial valleys. Global optimization techniques are used to design an artificial medium termination at depth for a normal-mode representation of the field. A derivation of a horizontal source array to produce an incident plane wave gives, as a by-product, an extension of a previous Fourier-transform relation involving Bessel functions. Like purely numerical methods, such as finite differences and finite elements, the method can handle all kinds of (two-dimensional) anomaly shapes.


Introduction
As a complement to purely numerical methods, analytical and semi-analytical techniques are of interest for giving physical insight and for providing bench-mark solutions to wave-propagation problems.This paper presents a semi-analytical method that can separate multiply scattered wavefield components and aid physical interpretation in connection with laterally displaced medium anomalies in a 2D (two-dimensional) waveguide.A discrete coupled-mode method, based on local modes and modal reflection matrices in a medium that is invariant in one of the horizontal directions, is developed for this purpose.The method has the same flexibility concerning anomaly shapes as purely numerical methods.With an artificial medium termination at depth, the wavefield is expanded in terms of normal modes in each of a number of laterally homogeneous strip regions with vertical boundaries.For simplicity, full details and examples are only given for the scalar case with plane-wave incidence from below.
The coupled-mode method has been applied to wave-propagation problems in various disciplines such as acoustics, seismology, electromagnetics, etc.; see [1] for a brief review.For definiteness, seismic waves are considered in the present paper.The discrete variant of the method, with local modes, i.e. modes adapted to the local structure, has been common in underwater acoustics for a long time [2][3][4][5].It is easy to implement numerically using matrix algebra to solve a two-point boundary-value problem for modal expansion coefficients.For the related continuous variant, without discretization of the medium into laterally homogeneous regions, the corresponding boundary-value problem involves coupled differential equations for the modal expansion coefficients [6].In surface-wave seismology, the continuous variant has been developed for reference modes, i.e. modes for a reference structure [7] as well as local modes [8,9].Maupin [10] gives a good review of both alternatives.With a sloping boundary, the continuous variant involves field expansion with modes that do not fulfil the correct boundary conditions.The implied slow convergence of the mode series can be improved by including artificial boundary modes [1,11].
To avoid issues with numerical stability, the two-point boundary-value problem can be recast as an initial-value problem for generalized modal R/T (reflection/transmission) matrices using invariant embedding.For the continuous coupled-mode variant, this has been done, involving differential equations of Riccati type, for reference modes [7] as well as local modes [12].There are other ways to achieve numerical stability, e.g.[13], but modal R/T matrices are crucial in the present paper to separate wavefield components.Hence, the invariant embedding technique is adapted here for discrete coupled local modes using matrix algebra.Explicit computation of transmission matrices is avoided by stabilized back-propagation of modal expansion-coefficient vectors.Mode matching is applied across vertical stripregion boundaries by a mathematically exact Galerkin approach, while, e.g.[14] applies approximate mode matching invoking Snell's law in each of a number of horizontal sections.
An advantage with the discrete variant of coupled local modes is that the modes need not be tracked carefully as functions of horizontal position to avoid mixing them up.To alleviate this problem for the continuous variant, in the acoustic case with a lossy medium, Pannatoni [15] suggests expansion in terms of local modes for the corresponding lossless medium.
Together with the coupled-mode method, addition rules for R/T matrices, as developed in [16, Sec.6.1] and sometimes called Redheffer star products, are used to achieve the desired separation of multiply scattered wavefield components.In essence, known techniques are combined and adapted to treat multiple scattering among laterally displaced medium anomalies in a 2D waveguide.
The application examples connect to recent seismic hazard research concerning amplification of incident plane SH (shear horizontal) waves, and related dynamic stress concentration, caused by 2D medium anomalies close to the surface.Single and multiple inclusions of various types have been considered [17][18][19], as well as topography irregularities [20] and alluvial valleys [21][22][23].Numerical methods like the boundary-element method (BEM), which is used in several of the mentioned papers, are very flexible concerning variations of anomaly shape and type.Ba & Yin [24] present a multidomain BEM for complex local sites, such as a multilayered half-space with inclusions.They decompose the half-space into suitable regions and set up a linear equation system to solve for densities of fictitious uniformly distributed loads on the region boundaries.Shyo & Teng [21] and Kara [22] apply hybrid methods involving finite elements and finite differences, respectively.The spectral element method is another attractive numerical technique that can be applied in this context [25,26].
An advantage with the coupled-mode approach, compared to the other methods referred to, is the direct applicability of addition rules for R/T matrices to isolate multiply scattered wavefield royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.10: 230352 components.Compared to the other (semi-)analytical methods mentioned, which are specially designed for particular kinds of anomalies, there is no restriction on the anomaly shapes.
The plan of the paper is as follows.Section 2 introduces wavefield decomposition with partial waves in a medium invariant in one of the horizontal directions.The partial waves isolate various types of multiple scattering among laterally displaced medium anomalies.They are defined using discrete coupled-mode theory and combination of R/T matrices.This field decomposition is valid for combined P-SV ( primary, shear vertical) and SH waves.The remaining text only deals with the pure SH case, however.Section 3 provides a representation of a vertically or obliquely incident plane wave by a horizontal source array.An extension of a Fourier-transform relation involving Bessel functions from Watson [45] appears as a by-product.Section 4 presents details of the discrete coupled-mode method.In essence, this is an adaptation to the 2D SH case of the recent coupled-mode method for 2D fluid media with a three-dimensional (3D) point source in [46,Sec. V].The adaptation involves, in particular, handling of an incident plane wave by a horizontal source array, and design of the artificial medium termination.Section 5 continues the discussion from §2 of field decomposition into partial waves, allowing isolation of individual (multiply) scattered waves.Reflection-matrix recursion with successive restarts is here an essential ingredient.Section 6 shows how to obtain the full field in a periodic medium efficiently with computations restricted to a single unit cell.A few computational variants are briefly indicated in §7, before some concluding remarks in §8.Two short appendices provide some additions concerning §4.2 for a different type of upper boundary and [46, Sec.V], respectively.
Examples appear in § §4.4,5.3, 5.4 and 6.2.For convenience, all of them are taken from the recent study of scattering by multiple alluvial valleys in [41].With the exception of the broad-band example in §5.4,handled by frequency synthesis, a harmonic time dependence according to the typically omitted factor exp(−iωt), where ω is the angular frequency and t is the time, is assumed throughout the paper.

Wavefield decomposition in a medium invariant in one of the horizontal directions
Consider a solid medium that is invariant in the y-direction, where x, y, z are Cartesian coordinates increasing to the east, north and downward, respectively.Below an upper free boundary, the medium agrees, for simplicity, with a certain laterally homogeneous reference structure, except in a number of laterally displaced anomaly regions, denoted A, B, … and defined by Except in the anomaly regions, the upper boundary is horizontal at a fixed depth.The medium is homogeneous below a certain depth, from where a plane time-harmonic upwards directed wave with direction vector (k x , k y , k z ) is incident.Figure 1a gives an illustration.The anomalies, i.e. irregularities in the x-direction, can involve variations of surface topography or medium parameters (density and velocities) or a combination of both.For ease of illustration, figure 1 shows the first case.
Between two consecutive anomaly regions, possibly also at the left end and at the right end, there is a connection region agreeing, for simplicity, with the reference structure.Denoting the displacement field by u(x) = (u(x), v(x), w(x)) T , where x = (x, y, z) and u, v, w are the components in the x-, y-, z-directions, respectively, the idea is now to separate different wavefield components or partial waves.A basic partial wave, denoted u 0 , is in each anomaly or connection region as if this medium part was embedded within the reference structure.(A laterally homogeneous medium obviously results for a connection region.)Additionally, each anomaly region A, B, … provides an excitation of its left and right connection regions.In figure 1a these excitations, to be more carefully defined in §2.1, are denoted b X and a X , where X = A, B, … In effect, each anomaly region provides sources for additional partial waves, to be obtained by multiple scattering among the anomaly regions.
Coupled-mode theory provides a useful method to extract these effective sources and additional partial waves.For the discrete variant, consider a discretization of the y-independent solid medium by N + 1 laterally homogeneous strip regions.With N ≥ 0 vertical interfaces 3, …, N, while (when N > 0) strips 1 and N + 1 include x < x 1 and x N < x, respectively.(When N = 0, strip 1 covers the whole x-axis.)The corresponding density, P-wave velocity and shear-wave velocity functions are denoted ρ n (z), α n (z) and β n (z), respectively, for n = 1, 2, …, N + 1.With absorption in the medium, the velocity functions are complex-valued.Since the upper medium boundary is free, free upper horizontal and vertical boundary segments of the strip regions are also appropriate at the medium discretization, cf.[46, the first two paragraphs of Sec.IV].The upper horizontal surface of strip n is at z = z a;n .Figure 1b provides an illustration, with varying z = z a;n to royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.10: 230352 mimic the surface topography in figure 1a.The required fineness of the x-discretization is related to the wavelength and the incidence directions of the waves [47].Of course, only the anomaly regions A, B, … need to be further discretized in this way, each connection region being a single connection strip.For simplicity, the z a;n of the connection strips agree.
To allow a mode representation of the displacement field in each strip region, the medium is artificially truncated with a free horizontal boundary at a finite lower depth z = z b .This is the lockedmode medium approximation, used by, e.g.[48,49].As detailed in §4.1, the medium absorption increases gradually towards the truncation boundary to minimize undesired reflections from it.There is a horizontal plane source array at depth z = z s , below all receivers and below all anomalies, to produce the incident plane wave.As indicated in figure 1b, the source array also produces a downwards directed plane wave.Below z s −, the medium is laterally homogeneous.Except for the artificial medium truncation, it is even homogeneous there.

Modal coefficient vectors
In each strip region, the wavefield is expanded in terms of normal modes, with coefficient column vectors a or a for waves to the right (direction of increasing x) and b or b for waves to the left (direction of decreasing x).The coefficient vectors are a, b and a, b at the left and right ends of the strip, respectively.In terms of modal coefficient vectors, each anomaly region X = A, B, … provides an excitation of its surroundings with a vector b X to its connection strip to the left and a vector a X to its connection strip to the right.These vectors are obtained as the difference between the corresponding coefficient vectors computed for the anomaly region embedded within the reference structure and for the reference structure throughout, in both cases ignoring the other anomaly regions.Figure 1a gives an illustration.
In effect, these b X and a X vectors, one pair for each of the anomaly regions, provide sources for additional partial waves to be computed without further consideration of the sources at z s .Upon (a)  royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.10: 230352 propagation through its connection strip, b X becomes a vector b X and a X becomes a vector a X , for excitation of the neighbouring anomaly regions, respectively.

Anomaly regions as two-ports
Disregarding the sources at z = z s , consider each anomaly region between two surrounding connection strips as a two-port with input field vectors a and b from the left and right connection strip, respectively, and corresponding output field vectors b and a to the left and right connection strip, respectively; cf.figure 2. For the application to a particular partial wave, only one of a and b is nonvanishing.With appropriate reflection matrices R, R and transmission matrices T, T, expressing mode conversions, Still disregarding the sources at z = z s , consider two consecutive anomaly regions, A and B for example, with surrounding and intermediate connection strips.The corresponding R/T matrices are R A , R A , T A , T A and R B , R B , T B , T B , respectively, and Ê denotes the diagonal transmission matrix of the intermediate connection strip.Figure 2 gives an illustration of this composite structure, with input vectors a and b, and output vectors b and a.These vectors are field vectors in the surrounding connection strips, while a c and b c denote field vectors in the intermediate connection strip.Lateral homogeneity of a connection strip allows unambiguous specification of x-direction for the waves there.(With corresponding a c ¼ 0 or b c = 0, a connection strip could be one of the end strips n = 1 or n = N + 1, respectively, of the medium.) There are well-known addition rules for R/T matrices [16,Sec. 6.1].In some texts, they are called Redheffer star products.Concerning R, R, T, T for the composite structure in figure 2, one obtains royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.10: 230352 Note the appearance of the same reverberation operators as in the latter alternatives of equations (2.2)-(2.5).Expansion of the reverberation operators, i.e. expansion of the inverse matrices in geometric series [16], provides the basis for the additional partial waves of the field.With several consecutive anomaly regions, together with surrounding and intermediate laterally homogeneous connection strips, recursive application of equations (2.2)-(2.5)provides R/T matrices for the composite structure in terms of elementary R/T matrices for the individual anomaly regions and diagonal transmission matrices for the connection strips.Reflections back and forth between adjacent composite anomaly region structures are readily incorporated.Again, expansion of the reverberation operators provides the basis for the additional partial waves of the field.

Notation for the additional partial waves
As already mentioned, the basic partial wave is denoted u 0 .The additional partial waves are denoted u j0s j l ,...,j2,j1 where s ¼ À or þ and j 0 , j 1 , …, j l denote anomaly regions (A, B, …).A u j0À j l ,...,j2,j1 wave starts with a b X¼j0 source vector to the left in the connection strip to the left of anomaly region j 0 , while a u j0þ j l ,...,j2,j1 wave starts with an a X¼j0 source vector to the right in the connection strip to the right of anomaly region j 0 ; cf.figure 1a.If l ≥ 1, the partial wave u j0s j l ,...,j2,j1 is subsequently reflected when reaching the anomaly regions j 1 , j 2 , …, j l in this order, and each of these reflections only involves a single anomaly region.Note that non-adjacent anomaly region indices may here coincide to include all multiple reflections (scattering back and forth).
After its last reflection, the partial wave proceeds by transmission through connection strips and anomaly regions towards the left or right end of the medium, thereby contributing to the field.In a connection strip, only waves to the left or right contribute; waves in the opposite direction appear in other partial waves.Within an anomaly region, however, the transmission gives rise to interior reflections which contribute as well.
The total field in a particular anomaly region or connection strip appears by summing the contributions there by u 0 and the additional partial waves.There is no high-frequency ray approximation, and the field decomposition is mathematically exact.

The pure SH case
So far, the discussion of the waves in the y-independent solid medium has been rather general.For simplicity, however, the following development is restricted to the pure SH wave case without conversions between SH (Love) and P-SV (Rayleigh) modes at the vertical x = x n interfaces.Specifically, k y = 0 for the incident plane wave and u(x) = u(x,z) = (0, v(x,z), 0) T .Denoting the medium density and shear-wave velocity by ρ(x,z) and β(x,z), respectively, and introducing the corresponding Lamé parameter μ(x,z) = ρ(x,z) β 2 (x,z), the y-component v(x,z) of the displacement satisfies the SHwave Helmholtz equation: Here, f (x,z) is the y-component of the body force per unit volume, Δ is the Laplacian, and grad is the gradient operator.Moreover, the components τ xy and τ yz of the stress tensor τ appear as τ xy = μ ∂v/∂x and τ yz = μ ∂v/∂z.This follows readily from basic linear elasticity theory (e.g.[50, ch. 2]).The medium velocity β may be complex-valued to accommodate absorption.The coupled-mode SH-wave method of the present paper is related to the method for the fluid case in [46, Sec.V].Actually, a fluid-medium analogy of the SH case appears by setting the density to 1/μ (possibly complex-valued), the sound speed to β, and the pressure to v. A free/rigid solid-medium boundary (with vanishing body force) with normal in the xz-plane then corresponds to a rigid/free fluid-medium boundary.

Generating a plane SH wave by a line array of line sources
The development in [46, Sec.V] concerns a 3D point source (or line source), while SH papers such as [27,41] typically prescribe an incident plane wave.Generation of a plane wave by an royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.10: 230352 array of point or line sources in a homogeneous solid medium with shear-wave velocity β is now considered.
At first, consider equation (2.8) with its right-hand side replaced by r 2 ref dðx À x s Þdðy À y s Þdðz À z s Þ for a 3D point source at x s = (x s , y s , z s ).Here, δ is the one-dimensional Dirac delta function with dimension m −1 , and r ref is a reference length, introduced to achieve consistency of dimensions.The solution for v(x,y,z), now depending on y as well because of the 3D point source, becomes (e.g.[51,) Next, consider a horizontal line source in the y-coordinate direction, according to r ref δ(x − x s ) δ(z − z s ) in the right-hand side of equation (2.8).The solution for the displacement v(x,z) in the y-direction becomes vðx, zÞ q (e.g.[4, Sec.5.2.2]).A verification can be obtained by integrating pointsource solutions according to equation (3.1) over y s and applying an integral representation of the Hankel function [52, 9.1.24].
The solution for a horizontal plane source according to exp(ik where the square root is taken with a non-negative imaginary part.This is readily seen by isolating the dependence on x according to the factor exp(ik x x), applying the radiation conditions and observing the implied step discontinuity of ∂v/∂z across z = z s .Originating from z = z s , cf. figure 1b, there is obviously an upwards directed plane wave and a downwards directed one.The solution according to equation (3.3) can be obtained in two other ways: by integration of line sources over a line or point sources over a plane.Concerning the first alternative, an application of

Coupled-mode SH-wave computations
Consider figure 1b with a discretization of the y-independent solid medium into N + 1 laterally homogeneous strip regions.Let m n ðzÞ ¼ r n ðzÞb 2 n ðzÞ, n = 1, 2, …, N + 1.Before returning to the partial waves in §5, the present §4 deals with coupled-mode computation of the full wavefield.
There is a horizontal plane source array at z = z s corresponding to the right-hand side in equation (2.8).This is a slight generalization of the particular choice w(x) = exp(ik x x) in §3.As in §2, the depth z s is below all receivers and below all anomalies (irregularities in the x-direction), such that the medium is laterally homogeneous below z s −.
For n = 1, 2, …, N + 1, let k m,n and Z m,n (z), m = 1, 2, …, denote the modal horizontal wavenumbers and normalized mode functions, respectively, for SH modes in strip region n.The wavenumbers k m,n appear in the upper complex plane but not on the negative real axis.It is appropriate to order them according to decreasing real parts and (close to the imaginary axis) increasing imaginary parts.The mode normalization means that Z m,n (z) denotes the original mode function Z 0 m,n ðzÞ divided by ½ Ð z b za;n m n ðzÞðZ 0 m,n ðzÞÞ 2 dz 1=2 .According to Sturm-Liouville theory, the modes in strip n are orthogonal with respect to the weight function μ n (z).Well established methods exist to compute the k m,n and the Z m,n (z) (e.g.[4,54]).For robust and automatic computations of the k m,n with winding-number integrals, it can be convenient to use quadruple precision. Introduce for waves to the left (decreasing x) with x-reference at x n , such that, within strip region n, vðx, zÞ where

Artificial medium truncation
The issue here, cf.[46,Sec.IV E], is to design an artificial medium truncation, ending at z = z b , providing negligible reflections.Recalling figure 1b and equation (3.3), the reflections from the downwards directed plane wave must be minimized.Moreover, there are undesired reflections at depth from the downwards directed surface reflections of the upwards directed plane wave.Because of the medium anomalies, these waves may have a slightly different incidence angle onto the artificial medium truncation.Finally, waves are of course scattered downwards in various directions from the anomalies.
The most important incidence angle, or k x , to be handled by the artificial medium truncation is obviously that of the downwards directed plane wave.For rapid truncation-design computations, consider two related variants of the medium, with source array according to exp(ik x x) δ(z − z s ) in the right-hand side of equation (2.8): (i) a homogeneous medium with density ρ 0 , shear-wave velocity β 0 and v(x,z) given by equation (3.3) with β = β 0 , and (ii) its modification with laterally homogeneous artificial medium truncation including absorption between z = z s and z = z b .
For medium (ii), β (and μ) depend on z below z = z s .It follows that v(x,z) = exp(ik x x) Z(z), where Z(z) fulfils ð4:8Þ and x ZðzÞ ¼ 0 ð4:9Þ for z < z s and z s < z < z b , respectively.Here, μ(z) = ρ 0 β 2 (z).In addition, allowing numerical determination of the constant γ for specified complex-valued β(z) and a real k x , there is the source-discontinuity condition Z 0 (z s+ ) − Z 0 (z s− ) = 1 and the free-boundary condition Z 0 (z b ) = 0.For an appropriate medium truncation, γ should apparently be close to 1, making the solutions for media (i) and (ii) similar for z < z s .For a trial z b , evolutionary optimization algorithms, such as differential evolution, are useful to minimize |γ − 1| by varying β(z).It is thereby convenient to specify the complex-valued function β(z) by the following three design parameters: the depth between z s and z b for onset of artificial absorption, β(z b ), and the polynomial degree of the change towards z b (1 for linear change) of β −2 (z).Of course, a large z b allows small artificial reflections, but it necessitates a large number of normal modes for the field representation.Hence, z b (as well as z s ) should be reasonably small, while still providing a negligible |γ − 1| at the minimization.
Numerical experiments with different z b are useful.Small incidence angles (close to normal incidence) may allow a small z b with a large absorption gradient from a shallow onset, while large incidence angles (close to grazing incidence) may require a large z b with a smooth absorption increase.Recalling that waves from the anomalies may be scattered in various directions onto the artificial medium truncation, it is a good idea to include several incidence angles, i.e. several k x and γ, in the minimization.

Reflection-matrix recursion
By physical arguments, considering F n waves to the right and C n waves to the left, there must exist modal reflection matrices R n with x-reference at x n and R n with x-reference at x n−1 , which are royalsocietypublishing.org/journal/rsosR. Soc.Open Sci.10: Accounting for changes of x-reference from x n−1 to x n and vice versa, cf.[46, eqns (37)-( 38)], respectively.Note that equations (4.12) and (4.13) are irrelevant when n = N + 1 and n = 1, respectively, since With the chosen normalization of the normal modes and the basic wave functions, reciprocity arguments of the same type as in [46,Sec. VI] show that all modal reflection matrices R n and R n are symmetric.The integrals in equations (4.6) and (4.12)-(4.13)are easy to compute analytically when w(x) = exp(ik x x).Note that the wavenumbers k m,n have positive imaginary parts because of the artificial absorption above z b .
Corresponding to the Riccati-equation solutions for the related continuous coupled-mode approach in [7], the modal reflection matrices R n and R n may be computed recursively, for decreasing n starting with R N+1 = 0 and for increasing n starting with R 1 ¼ 0, respectively.The recursion equations are derived, by a Galerkin approach, from the continuity of v and τ xy at the vertical interfaces separating the strip regions, cf.

Recursion of R n for decreasing n
By the boundary conditions, τ xy = μ ∂v/∂x vanishes when x = x n and min(z a;n , z a;n+1 ) < z < max(z a;n , z a;n+1 ).Introduce the notation where ð4:21Þ ) subsequently provides the reflection-matrix recursion equation: ð4:23Þ When I n+1 ⊆ I n , on the other hand, integrations over I n+1 and I n provide for continuity of v after multiplication with μ n+1 (z)Z m,n+1 (z) and τ xy after multiplication with Z m,n (z), respectively.This yields where now ð4:26Þ The equation for reflection-matrix recursion follows from equation (4.25) as When I n+1 = I n , a third option is possible by combining equations (4.24) and (4.20).This leads to, cf.[55, Sec.III A], This third option is more efficient computationally, but it requires that the kept finite numbers of modes in the two adjacent strips are equal.

Recursion of R n for increasing n
By the boundary conditions, τ xy vanishes when x = x n−1 and min(z a;n , z a;n−1 ) < z < max(z a;n , z a;n−1 ).

Introduce the notation b
for continuity at x n−1 of v after multiplication with μ n (z)Z m,n (z) and τ xy after multiplication with Z m,n−1 (z), respectively.Insertion of the relation a where now ð4:34Þ royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.10: Rþ n , the equation for reflection-matrix recursion follows as When I n−1 ⊆ I n , on the other hand, integrations over I n−1 and for continuity of v after multiplication with μ n−1 (z)Z m,n−1 (z) and τ xy after multiplication with Z m,n (z), respectively.This yields where now ð4:39Þ The equation for reflection-matrix recursion becomes When I n−1 = I n , and the numbers of modes kept in the two strips n − 1 and n are equal, a third option is possible by combining equations (4.37) and (4.33).This leads to, cf.[55, Sec. ð4:43Þ

Stabilized back-propagation
To compute a n and b n for equation (4.6), substitution of their reflection-matrix expressions from equations (4.10)-(4.11)into equations (4.12)-(4.13)yields an equation system with the solution , respectively, safely control these computations.
royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.10: 230352 .The medium, without absorption and with a flat surface at z = 0 (km), is homogeneous except for three similar semicircular alluvial valleys with radius a, centred at (x,z) = (−4a,0), (0,0) and (4a,0), respectively.With ρ and β denoting alluvium density and shear-wave velocity, and ρ 0 and β 0 denoting corresponding values for the surrounding bedrock, ρ/ρ 0 = 2/3 and β/β 0 = 1/3.A plane SH wave with frequency β 0 /2a is incident from below at three different angles to the horizontal plane: (a) 5°, (b) 45°and (c) 90°(vertical incidence).In (a) and (b), the wave direction is to the right (increasing x).For all numerical results, there are certain control parameters affecting the trade-off between accuracy and efficiency.The main ones for the present coupled-mode computations are the number of strip regions (N + 1), z b and the other parameters for the artificial medium truncation, and the number of included normal modes in each strip region.For figure 3, each alluvial valley is discretized with about 70 strip regions of varying thickness.Compared to the wavelength (2a in the bedrock), a fine xdiscretization is needed, adapted to the local alluvium-bedrock interface slope [47].Furthermore, z b = 50a for figure 3a and z b = 15a for figure 3b,c.In each case, differential-evolution optimization according to §4.1, with z s = 3a, provides suitable parameters for the artificial absorption at depth.The number of included normal modes in each strip region is 503 for figure 3a and 200 for figure 3b,c.

Example
The surface-displacement amplitude results show satisfactory agreement with the corresponding ones (without absorption) in [41, fig. 5], indicated by star symbols in figure 3. Note the significant ground-motion amplification within the soft valleys.Also note the clear wave shielding of the middle valley by the left one in figure 3a with almost horizontal incidence from the left.The shielding of the right valley is less effective because of the slightly non-horizontal incidence, the greater distance from royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.10: 230352 the left valley, and the shielding of the middle one.Using partial-wave decomposition, this example is revisited in §5.3 for quantitative assessment of (multiple) scattering by the alluvial valleys.
After the design of the artificial medium truncation, the computations involve the following four steps: computation of the modal horizontal wavenumbers k m,n , computation of the normalized mode functions Z m,n (z), computation of the mode-coupling matrices F n , G n and F n , G n , and the actual coupled-mode computations with reflection-matrix recursion and stabilized back-propagation according to § §4.2 and 4.3.The first step is typically the most time-consuming.Note that the results from the first three steps can be reused for computations with other directions of plane-wave incidence, for example.
Assume, for simplicity, that the finite number of modes included in the computations, denoted N M , is about the same in all strip regions.In general, the needed N and N M are both roughly proportional to the frequency.For the last of the four computation steps, the computational work is roughly proportional to N N 3 M , since the reflection-matrix recursions involve inversions of matrices of size N M × N M .Some variations may of course occur, depending on the particular recursion option from §4.2 and the number of receivers.For figure 3, with a standard serial home personal computer, the corresponding CPU times are 481 s for figure 3a, and 27 s for each of figure 3b,c.The root-mean-square deviations from results obtained by doubling N and N M are less than 1% of the maximum amplitude, indicating the accuracy.

Field decomposition into partial SH waves
The partial waves from §2 are now simply denoted v 0 (the basic one), and v j0s j l ,...,j2,j1 (the additional ones).As before, s ¼ À or þ and j 0 , j 1 , …, j l denote anomaly regions (A, B, …).A v j0À j l ,...,j2,j1 wave starts with a b X¼j0 source vector to the left involving a C function, while a v j0þ j l ,...,j2,j1 wave starts with an a X¼j0 source vector to the right involving a F function.

Reflection-matrix recursion with successive restarts
It is in fact easy to adapt a computer program for computation of the full field to computation of partial waves for specified anomaly regions and corresponding connection strips.Still apply the reflectionmatrix recursions according to equations (4.23), (4.28), (4.30) and (4.36), (4.41), (4.43) throughout the whole medium, from n = N to n = 1 and from n = 2 to n = N + 1, respectively.Now, however, restart the recursion with a vanishing reflection matrix in the right-hand side upon entry to an anomaly region from one of the connection strips.This procedure automatically provides the elementary reflection matrices R A , R B , … and R A , R B , … for the anomaly regions, cf.figure 2. It is not necessary to compute the transmission matrices T A , T B , … and T A , T B , … explicitly.Typically, transmission is best handled by sequential matrix-vector multiplications according to the technique with stabilized back-propagation from §4.3.

Computation of the additional partial waves
The transmission of a v j0s j l ,...,j2,j1 wave through a connection strip is easily done, using the appropriate Ên matrix.In particular, the pertinent b X¼j0 or a X¼j0 vector is initially transmitted in this way.The transmission through the involved anomaly regions is done with stabilized back-propagation, incorporating the backward-going waves; cf.§4.3.The reflection matrices obtained by recursion with successive restarts according to §5.1 are the appropriate ones for this purpose, as well as for the reflections from the anomaly regions.Since the additional partial waves do not include source terms, except for the initial b X or a X vector, the involved equations (4.18) and (4.31) are simplified.
Concerning the reflections from the anomaly regions, it is instructive to consider a medium with two anomaly regions, A and B, separated by connection strip n.In strip n, vðx, zÞ royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.10: 230352 Note the appearance of the same inverse matrices as in equations (4.44) and (4.45).Expansion of the inverse matrices in geometric series provides the additional partial waves.When A and/or B are split, with more anomaly regions and connection strips, the corresponding refined partial-wave decomposition follows from the two-port algebra according to §2.2 together with the elementary reflection matrices according to §5.1.

Narrow-band example
Returning to the example case from §4.4, introduce three anomaly regions: A for x < −3a, B for |x| < a and C for x > 3a.Note that each anomaly region covers one of the alluvium-valley anomalies.There are two intermediate connection strips: one for −3a < x < −a, and one for a < x < 3a.Restriction is now made to the case from figure 3b, with incident-wave propagation angle 45°.
Figure 4a shows the corresponding basic partial wave v 0 .Note the constant relative amplitude 2 in the two connection strips, resulting from the reflection of the plane SH wave at the free surface.Some modulation of the relative amplitude 2 appears for x < −5a in anomaly region A and for x > 5a in anomaly region C, because of scattering from the alluvial valley in A and C, respectively.There is some asymmetry in the field because of the oblique incidence.As expected, the v 0 fields within the three alluvial valleys agree.
Figure 4b shows the (coherent) sum of all essential additional partial waves v j0s j l ,...,j2,j1 .This field results from remaining effects of single and multiple scattering by the alluvial valleys.(As already noted, some single-scattering effects appear in figure 4a, for |x| > 5a there, because of the extension of anomaly regions A and C beyond their pertinent alluvial valleys.)In each anomaly region, particularly the middle one B, the scattering contributions from the other anomaly regions are significant.The scattering into the two connection strips is also essential.The (coherent) sum of the waves in figure 4a,b agrees very well with figure 3b.
As a complement to figure 4, with a different amplitude scale, figure 5 shows some individual additional partial waves v j0s j l ,...,j2,j1 .Except v A+ and v C− , in figure 5a,e, respectively, all of them involve multiple scattering (reflections) by the alluvial valleys.Note the successively decreasing amplitudes in each of the two upper rows, because of the increasing number of reflections.Figure 5k-n shows, cf.v Aþ A,B,A,B and v CÀ C,B,C,B in figure 5o,p with five reflections, that transmission through anomaly regions may imply larger amplitude losses than reflections.
A plane SH wave is incident from below at the propagation angle 45°relative to the positive x-axis.The spectrum of the source pulse, not exactly as in [41,Sec. 5], is limited to the frequency band (0.07, 1.53) Hz.Coupled-mode computations for time-domain results are performed with Fourier synthesis using about 50 discrete frequencies within this band.Each alluvial valley is thereby discretized with about 35 strip regions, for simplicity the same number for all frequencies.The source-array depth z s is 2 km, while the depth z b varies from 5 km for the highest frequencies to 21 km for the lowest ones.As in the previous example, differential-evolution optimization according to §4.1 furnishes appropriate parameters for the artificial absorption.In this broad-band case, good results are obtained with only 45 normal modes in each strip region, for each frequency.With averaging over frequencies, fewer modes are apparently needed than for narrow-band cases.
Figure 6a shows 100 time traces for the surface displacements between x = −10 km and x = 10 km and a certain time window of length 30 s.The arrivals within |x| < 4 km can be favourably compared to those in [41, fig.10(e)].
To aid the interpretation of the arrivals, using partial waves, introduce three anomaly regions: A for x < −7 km, B for |x| < 1 km and C for x > 7 km.Note that each anomaly region covers one of the alluvium-valley anomalies.There are two intermediate connection strips: one for −7 km < x < −1 km, and one for 1 km < x < 7 km.
Figure 6b shows the corresponding basic partial wave v 0 .Within the connection strips, there is of course only one arrival, the direct one as doubled by the surface reflection.Within an alluvial valley, there is a delayed direct arrival followed by a later reflection from the semicircular valley boundary; cf.royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.10: 230352 Figure 6c shows the (coherent) sum of all essential additional partial waves v j0s j l ,...,j2,j1 .The (coherent) sum of the waves in figure 6b,c agrees very well with figure 6a.Note the rather strong waves at (x, t) ≈ (−6.9 km, −1 s) and (x, t) ≈ (1.1 km, 4.5 s) in figure 6c, caused by scattering by the right ends of the left and middle alluvial valleys, respectively.By destructive interference with the corresponding displacements in figure 6b, some slight wave shielding appears in figure 6a.This wave shielding is of course more significant at more horizontal incidence [41,Sec. 5].
Magnified compared to the previous figure, figure 7 shows some individual additional partial waves v j0s j l ,...,j2,j1 .It is clear that the rather strong waves at (x, t) ≈ (−6.9 km, −1 s) and (x, t) ≈ (1.1 km, 4.5 s) in figure 6c belong to v A+ shown in figure 7a and v B+ shown in figure 7b, respectively.Reflections from the front as well as back sides of an alluvial valley give rise to a clear doublet structure of the partial waves in figure 7c-f.A doublet structure, albeit weak, can be discerned in figure 7a,b too.It is caused by reflections back and forth within anomaly regions A and B, respectively.Forward scattering from the alluvial valley within anomaly region j 0 causes prolongation of the initial arrival for v j0þ j l ,...,j2,j1 .This is clearly seen for v A+ , v B+ and v Bþ C in figure 7. Within a traversed alluvial valley, reflections from its far side appear for each of the additional partial waves in figure 7. Note that the multiply scattered (reflected) waves v BÀ A and v Bþ C , significantly magnified in figure 7e,f, are weak and barely notable in figure 6c.Effects of multiple scattering can be larger when the anomalies are closer together.
For computation of the full field, the aim is now to reduce the computations to a single period or unit cell: the one between x 1 and x N .To that end, regard this part of the medium as a two-port with input field vectors a and b from the strip regions to the left and right, respectively, and corresponding output field vectors b and a to these strip regions, respectively.Figure 8 gives an illustration.

Computations for the unit cell with one medium period
A difference from the two-port discussion in §2.2 is that the sources at z = z s within the two-port now contribute.In the present case, cf.equation ( 2 royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.10: 230352 with b S and a S representing the contributions from the sources.These source vectors are readily computed according to §4 for the related non-periodic medium with strip regions 1 and N + 1 extending to x = −∞ and x = ∞, respectively, and with the source function w(x) set to zero outside the unit cell.Assume, for simplicity, that the periodic medium is laterally continuous across x = x 1 and x = x N , such that z a;1 = z a;2 = z a;N = z a;N+1 , ρ 1

:4Þ
The field in the unit cell follows by summing the solution for the related non-periodic medium and the transmitted fields arising from the vectors a from the left and b from the right, respectively.These transmitted fields are efficiently computed by stabilized back-propagation, this time using the reflection matrices available from the handling of the related non-periodic medium.
Explicit computation of the transmission matrices T and T is actually needed in this case, to compute b and a from equations (6.3) and (6.4).Concerning T, matrices from equations (4.21), (4.26) and (4.29) must be multiplied, while T involves matrices from equations (4.34), (4.39) and (4.42).interrupted by an infinite number of periodically distributed down-going semicircular rigid (!) boundaries with radius a, centred at (x,z) = (8la,0) for l = 0, ±1, ±2, ….As in §4.4, a plane SH wave with frequency β 0 /2a, where β 0 is the shear-wave velocity, is incident from below at three different angles to the horizontal plane: (a) 5°, (b) 45°and (c) 90°(vertical incidence).In (a) and (b), the wave direction is to the right (increasing x).So far, the upper (solid-)medium boundary has been assumed to be free.As detailed in appendix A, the rigid case necessitates some modifications of § §4.2 and 3.

Example
The coupled-mode computations for figure 9 are restricted to the unit cell with |x| < 4a, and the involved semicircular anomaly at |x| < a is discretized with about 70 strip regions of varying thickness; cf. the x-discretization in §4.4.The depths z b and z s , as well as the parameters for artificial absorption, are as in §4.4,for each of the three incidence-angle cases.
The surface-displacement amplitude results in figure 9   royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.10: 230352 displacements at the semicircular boundary.The results in the two previous papers were obtained using high-velocity semicircular inclusions to mimic the desired rigid case, leading to non-vanishing corresponding displacements.Using the technique with partial waves from § §2 and 5, it would be easy to investigate how the response in figure 9 arises from multiple scattering among the down-going semicircular rigidboundary anomalies.Result curves of the same type as in figures 4 and 5 would appear.

Computational variants
The present paper focuses on modal reflection matrices as a convenient tool to transport boundary conditions along the x-axis: R = 0 from the right end, and R ¼ 0 from the left end.These reflection matrices relate coefficient column vectors for an expansion of v(x, z) in each strip n in terms of the row vectors F n ðx, zÞ and C n ðx, zÞ; cf.equations (4.6), (4.10) and (4.11).
By the definition of F n ðx, zÞ and C n ðx, zÞ, expansions of μ n (z) ∂v(x,z)/∂x and −iωv(x,z) in each strip n in terms of the corresponding row vector Z n (z) = {Z m,n (z)} follow from the previous expansion of v(x, z).Linear relations between the corresponding coefficient column vectors appear with so-called impedance and admittance matrices, which are easy to relate to the reflection matrices; cf.[57, eqn (13)].Obviously, related recursions for impedance matrices could be used instead of the reflection-matrix recursions in §4.2, and the boundary conditions could be transported using impedance matrices.Some texts, e.g.[57] for shape optimization of acoustic horns and [58] for acoustic simulation of the vocal tract, apply royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.10: 230352 the impedance-matrix approach.Just as for the reflection matrices [7], differential equations of Riccati type appear in the continuous case without medium discretization.
It is also natural to compare to well-known methods for computation of seismic P-SV waves in multilayered laterally homogeneous media.After Fourier-or Hankel-transformation of a horizontal coordinate to the wavenumber domain, a two-point boundary-value problem appears for a system of ordinary differential equations in the depth variable z.Several numerical methods have been proposed to solve this problem in an unconditionally stable way.For example, Kennett [16] combines 2 × 2 R/T matrices for the different layers recursively, while Wang & Rokhlin [59] compute a 4 × 4 global stiffness matrix recursively from individual layer stiffness matrices.Stiffness and compliance matrices relate stresses to displacements, and vice versa, combining both sides of a (composite) layer.Concerning v(x, z) in the present paper, a corresponding stiffness-matrix method would obviously relate the two expansion column vectors for μ n (z) ∂v(x,z)/∂x at two different x-values, with corresponding n-values, to the two expansion column vectors for v(x, z) at these two x-values.Equation (2.1) would be useful to relate stiffness and compliance matrices to R/T matrices, but these things are not pursued here.
To avoid spurious reflections from down-going waves, an artificial medium truncation involving a classical absorbing layer is carefully designed using global optimization in §4.1.It is a convenient choice, since standard methods for mode expansion with computation of modal wavenumbers and mode functions are directly applicable.There are good alternatives, however, which should be able to remove the spurious reflections with a much thinner artificial layer.Givoli [60] describes some milestones in the development of absorbing boundaries and layers, including Dirichlet-to-Neumann boundary conditions, PML ( perfectly matched layer), and high-order absorbing boundary conditions.The main interest has concerned applications for purely numerical methods, such as finite elements and finite differences.This is true for the PML approach too, with a recent review in [61], but this technique has attracted some interest within a modal framework as well.
With PML, so-called PML modes appear in addition to trapped and leaky modes [62].These PML modes are significant mainly within the PML region.The inclusion of PML modes in the modal basis may need some care, however.For a Pekeris waveguide, Zhu & Lu [63] provide approximate solutions for the PML modes, which replace the pertinent branch-cut integral.According to Zhu & Zhang [64], the eigenfunctions of the modified Helmholtz operator have no orthogonality in a bounded domain with a PML, and the authors derive pertinent conjugate eigenfunctions for the case of a Pekeris waveguide.

Concluding remarks
For a solid medium that is invariant in the horizontal y-coordinate direction, §2 presents a mathematically exact decomposition of the seismic wavefield with partial waves.With § §3 and 4 as additional background, details and examples for the scalar case with pure SH waves follow in §5.The decomposition is defined using discrete coupled-mode theory and combination of elementary reflection matrices, conveniently computed by recursion with successive restarts.It facilitates physical interpretation and allows detailed assessment of multiple scattering among horizontally displaced anomaly regions.Related field decompositions into partial waves have been briefly indicated at the ends of [55, Sec.V B] and [46,Sec. VI].The emphasis there is on reflections (or scattering) from the interiors and exteriors, or sides, of particular source and receiver regions.
Essentially as an adaptation of the 3D point-source case in [46, Sec.V], §4 develops the details of a discrete coupled-mode computation method for 2D SH-wave scattering at plane-wave incidence from below.The medium is discretized into a number of laterally homogeneous strip regions separated by vertical interfaces.A horizontal source array generates the incident plane SH wave according to equation (3.3) in §3.There is an artificial boundary at depth z b , allowing a normal-mode representation of the field in each strip region.Global optimization techniques are applied to design artificial absorption in a layer above this boundary to minimize reflections from it ( §4.1).
Recursion of modal reflection matrices and stabilized back-propagation of modal expansioncoefficient vectors are essential features of the computation method.Compared to the coupled-mode method for a 3D point source in [46, Sec.V], the introduction of a horizontal source array necessitates a double pass of the stabilized back-propagation: a full pass in each direction ( §4.3), to pick up source contributions from the left and from the right, respectively.To add a lot of point-or line-source contributions, with stabilized back-propagation in both directions from each, would not be efficient.

Figure 2 .
Figure 2. Vertical xz-plane with two consecutive anomaly regions A and B with a connection strip in between.R/T matrices, input and output field vectors from and to the surrounding connection strips, as well as field vectors for the intermediate connection strip, are indicated.

Figure 3
Figure 3 concerns an example from Zhang et al.[41, Sec.4.2].The medium, without absorption and with a flat surface at z = 0 (km), is homogeneous except for three similar semicircular alluvial valleys with radius a, centred at (x,z) = (−4a,0), (0,0) and (4a,0), respectively.With ρ and β denoting alluvium density and shear-wave velocity, and ρ 0 and β 0 denoting corresponding values for the surrounding bedrock, ρ/ρ 0 = 2/3 and β/β 0 = 1/3.A plane SH wave with frequency β 0 /2a is incident from below at three different angles to the horizontal plane: (a) 5°, (b) 45°and (c) 90°(vertical incidence).In (a) and (b), the wave direction is to the right (increasing x).For all numerical results, there are certain control parameters affecting the trade-off between accuracy and efficiency.The main ones for the present coupled-mode computations are the number of strip regions (N + 1), z b and the other parameters for the artificial medium truncation, and the number of included normal modes in each strip region.For figure3, each alluvial valley is discretized with about 70 strip regions of varying thickness.Compared to the wavelength (2a in the bedrock), a fine xdiscretization is needed, adapted to the local alluvium-bedrock interface slope[47].Furthermore, z b = 50a for figure3aand z b = 15a for figure3b,c.In each case, differential-evolution optimization according to §4.1, with z s = 3a, provides suitable parameters for the artificial absorption at depth.The number of included normal modes in each strip region is 503 for figure3aand 200 for figure3b,c.The surface-displacement amplitude results show satisfactory agreement with the corresponding ones (without absorption) in[41, fig.5], indicated by star symbols in figure3.Note the significant ground-motion amplification within the soft valleys.Also note the clear wave shielding of the middle valley by the left one in figure3awith almost horizontal incidence from the left.The shielding of the right valley is less effective because of the slightly non-horizontal incidence, the greater distance from

Figure 3 .
Figure 3. Coupled-mode surface-displacement amplitude curves for an example with three alluvium-valley anomalies as described in the text.A plane SH wave is incident from below at three different propagation angles relative to the positive x-axis: (a) 5°, (b) 45°and (c) 90°.The amplitude results are given relative to the incident-wave amplitude, and the star symbols indicate corresponding results from Zhang et al. [41, fig.5].

Figure 4 .
Figure 4. Further surface-displacement amplitude results for the case from figure 3 with incident-wave-propagation angle 45°: (a) the basic partial wave v 0 and (b) the (coherent) sum of the most significant additional partial waves, 16 of which are shown in figure 5.The amplitude results are given relative to the incident-wave amplitude, and the anomaly regions (A, B and C) are indicated.
[41, fig.10(b)] for |x| < 1 km.Of course, the v 0 fields within the three alluvial valleys agree.Because of the extension of anomaly regions A and C beyond their pertinent alluvial valleys, there is a reflected arrival from the valley in A for x < −9 km, and a scattered arrival from the valley in C for x > 9 km.For the middle valley, the corresponding arrivals are clearly seen in the total field of figure6a.In [41, fig.10(b,e)], they are denoted SL01 and SR01, respectively.

Figure 5 .
Figure 5.Some individual additional partial waves included in the sum in figure 4b.The first row shows the partial waves (a) v A+ , (b) v Bþ C , (c) v Aþ A,B and (d ) v Bþ C,B,C .The second row shows (e) v C− , (f ) v BÀ A , (g) v CÀ C,B and (h) v BÀ A,B,A .The third row shows (i) v Aþ C , ( j ) v Aþ B,C , (k) v Aþ C,B,C and (l ) v Aþ A,C .The fourth row, finally, shows (m) v BÀ C,A , (n) v Aþ C,A,B , (o) v Aþ A,B,A,B and ( p) v CÀ C,B,C,B .

Figure 6 .
Figure 6.Coupled-mode surface-displacement time traces for an example with three alluvium-valley anomalies as described in the text.A plane SH wave is incident from below at the propagation angle 45°relative to the positive x-axis.The three panels show: (a) total result, (b) the basic partial wave v 0 and (c) the (coherent) sum of the most significant additional partial waves, six of which are shown in figure 7. The anomaly regions (A, B and C) are indicated in (b,c).

Figure 7 .
Figure 7.Some individual additional partial waves included in the sum in figure 6c.The first and second rows show the partial waves (a) v A+ , (b) v B+ , (c) v B− and (d ) v C− .These waves are magnified three times in relation to those in figure 6.The third row shows (e) v BÀ A and ( f ) v Bþ C , magnified 20 times in relation to the waves in figure 6.
and that w(x) is a regular function at x = x 1 and x = x N .Then it follows by periodicity that a ¼ E Á a and b ¼ E Á b, ð6:2Þ where E = diag m (exp(ik x d)) = exp(ik x d ) I. Solution of equations (6.1) and (6.2) yields b

Figure 9
Figure 9 concerns an example from Zhang et al. [41, Sec.3.2], originally treated in[56, Sec.5.2].The medium, without absorption, is now homogeneous with a flat and free surface at z = 0 (km) interrupted by an infinite number of periodically distributed down-going semicircular rigid (!) boundaries with radius a, centred at (x,z) = (8la,0) for l = 0, ±1, ±2, ….As in §4.4, a plane SH wave with frequency β 0 /2a, where β 0 is the shear-wave velocity, is incident from below at three different angles to the horizontal plane: (a) 5°, (b) 45°and (c) 90°(vertical incidence).In (a) and (b), the wave direction is to the right (increasing x).So far, the upper (solid-)medium boundary has been assumed to be free.As detailed in appendix A, the rigid case necessitates some modifications of § §4.2 and 3.The coupled-mode computations for figure9are restricted to the unit cell with |x| < 4a, and the involved semicircular anomaly at |x| < a is discretized with about 70 strip regions of varying thickness; cf. the x-discretization in §4.4.The depths z b and z s , as well as the parameters for artificial absorption, are as in §4.4,for each of the three incidence-angle cases.The surface-displacement amplitude results in figure9agree well with the corresponding ones in [41, fig.3] and [56, fig.9], which are indicated by star symbols in figure9.Note the vanishing coupled-mode Figure 9 concerns an example from Zhang et al. [41, Sec.3.2], originally treated in[56, Sec.5.2].The medium, without absorption, is now homogeneous with a flat and free surface at z = 0 (km) interrupted by an infinite number of periodically distributed down-going semicircular rigid (!) boundaries with radius a, centred at (x,z) = (8la,0) for l = 0, ±1, ±2, ….As in §4.4, a plane SH wave with frequency β 0 /2a, where β 0 is the shear-wave velocity, is incident from below at three different angles to the horizontal plane: (a) 5°, (b) 45°and (c) 90°(vertical incidence).In (a) and (b), the wave direction is to the right (increasing x).So far, the upper (solid-)medium boundary has been assumed to be free.As detailed in appendix A, the rigid case necessitates some modifications of § §4.2 and 3.The coupled-mode computations for figure9are restricted to the unit cell with |x| < 4a, and the involved semicircular anomaly at |x| < a is discretized with about 70 strip regions of varying thickness; cf. the x-discretization in §4.4.The depths z b and z s , as well as the parameters for artificial absorption, are as in §4.4,for each of the three incidence-angle cases.The surface-displacement amplitude results in figure9agree well with the corresponding ones in [41, fig.3] and [56, fig.9], which are indicated by star symbols in figure9.Note the vanishing coupled-mode

Figure 8 .
Figure 8. Vertical xz-plane as in figure1b, but for the periodic medium with unit cell between x 1 and x N .The vectors a and b provide input from the surrounding strip regions, with corresponding output vectors b and a.

Figure 9 .
Figure 9. Coupled-mode surface-displacement amplitude curves for the example with peridically distributed semicircular rigid boundaries at |x − 8la| < a for l = 0, ±1, ±2, ….A plane SH wave is incident from below at three different propagation angles to the positive x-axis: (a) 5°, (b) 45°and (c) 90°.The amplitude results are given relative to the incident-wave amplitude, and the star symbols indicate corresponding results from Zhang et al. [41] and Ba & Liang [56].The Zhang et al. [41] and Ba & Liang [56] results are often somewhat different for the 90°angle.At each x, the two star symbols in (c) are, therefore, connected by a vertical line segment for visual clarity.
:8Þ By analytic continuation, this relation is valid for complex k and k x such that |Im(k x )| < Im(k).It is also royalsocietypublishing.org/journal/rsos R. Soc.Open Sci.10: 230352 valid (36) the definition of ðF n ðx, zÞÞ m is a kind of normalization, cf. the normalization of the modes.It leads to different dimensions of F n ðx, zÞ and C n ðx, zÞ, but it simplifies some forthcoming equations and it makes the appearing reflection matrices (see §4.2) symmetric.Additionally, introduce the diagonal matrices Ên ¼ diag m ð Êð1Þ m,n ðx n ÞÞ ¼ diag m ð Êð2Þ m,n ðx nÀ1 ÞÞ for transmission, or transfer of x-reference, between the sides of strip region n.Note that Ê1 ¼ ÊNþ1 ¼ I.At each x, expand the field in terms of the local modes there.It follows, cf.[46, eqns (35)-(36)], that there are coefficient column vectors a are boundary conditions.The mode sum in equation (4.6) follows by equation (4.1) and integration over x s of basic modal line-source solutions according to, e.g.[4, Sec.5.2.2].Section 4.1 gives some details on the medium truncation.Subsequently, essentially generalizing the development in [46, Sec.V] to a continuous source array, § §4.2 and 4.3 present a coupled-mode method to determine the column vectors a n , a L n , a R n and b n , b L n , b R n .Of course, only a finite number of modes is kept in each strip at the actual computations.