Interacting Seismic Waveguides: Multimode Surface Waves and Leaking Modes

Recentdevelopmentsinseismicrecordingprovidedensesamplingoftheseismicwavefieldthat allows the extraction of higher surface wave modes as well as the fundamental, and in some circumstances also P -dominatedmodes. Thecharacterofthemodalbranchesandtheirdispersionwithfrequencyandphase speed depends on interactions between different aspects of the seismic structure and between wavetypes. The influence of P waves becomes quite strong in the very near surface when there are strong wavespeed gradients. Aneffectivetoolforunderstandingthenatureofthemodalinteractionsduetostructureisprovidedby the seismic response in frequency–phase speed space, the kernel for seismogram calculation. Such displays for all three-componentsof the surfaceresponseextractthe full modal responseforRayleighand Lovewaves, including leaking mode effects, and an indication of the way that the modes are excited. Non-technical summary Seismic energy trapped between the free surface and the rapid increase in seismic wavespeeds with depth can propagate to substantial distances. With increasing density of seismic observations the full character of such guided waves can be revealed, with strong dependencies on interaction between different parts of the structure and wavetypes. In favourable circumstances, P -wave dominated features can be tracked in addition to the more commonly studied S dominated modes. Understanding the nature of such modes aids improved rendering of structure at depth.


Introduction
Surface waves have received increased attention in recent years, in part because of the success in exploiting ambient noise (see, e.g., Nakata et al., 2019, and references therein). With increasing station density and careful processing it has proved possible to extract the dispersion of both the fundamental mode and two or three higher modes from broad scale deployments of seismometers, at frequencies up to 0.5 Hz (e.g., Chen et al., 2022). Such S-dominated modes are fully trapped between the free surface and the increasing wavespeeds at depth. With very high density observations, as in the LASSO deployment of 1800 seismic nodes in Oklahoma, USA, the dispersion characteristics of partially trapped P-dominated modes have also been imaged in frequency-phase speed space (Li et al., 2022). Other classes of high-density sampling of the seismic wavefield such as the use of distributed acoustic sensing (DAS) also can allow the extraction of multi-mode dispersion including P modes in favourable circumstances (e.g., Fichtner et al., submitted 2023).
Although the concepts of surface wave dispersion are well established, the complexities of modal interaction at high frequency are rarely considered. Sediments, the crust, and the upper mantle each have the capacity to trap seismic energy and interactions between these different waveguides produce complexity in the charac- * Corresponding author: brian.kennett@anu.edu.au ter of surface wave modes. The dispersion curves for modes with different character apparently cross. Likewise, although Rayleigh waves are dominantly linked to the behaviour of the SV wavespeed with depth, where low P wavespeeds occur near surface there can be significant interaction between SV -and P-dominated modes. The nature of the waveguide interactions can be understood in terms of the controls on modal dispersion.
In a full three-dimensional body, such as the Earth, all deformations can be completely described by a superposition of the set of normal modes of the structure. However, when we consider a truncated portion of the structure in cartesian geometry, the properties at depth set a limit on the range of phase speeds for which surface wave modes are fully trapped. Waves with higher phase speed lose energy as they propagate horizontally by radiation into the lower halfspace, but can sometimes be recognised as 'leaking' modes. For a fully elastic medium, the trapped modes correspond to poles of the response lying in the real frequency-slowness plane. When leakage occurs, the singularities for the leaking modes occur on lower sheets in complex slowness at each frequency (e.g., Haddon, 1984). In the presence of attenuation all poles are shifted into the complex plane. However, for typical levels of seismic attenuation the trajectories of the poles in frequency-slowness space can be retrieved with a purely elastic approximation.
For a realistic seismic structure, including attenua-tion, the character of the modal field can be extracted by direct calculation in the frequency-slowness domain with modal branches defined by localised increases in amplitude reflecting the presence of poles off the top slowness sheet. In this way both trapped modes and leaking modes can be visualised, and the effects associated with the interactions of multiple seismic waveguides analysed.
Commonly the near-surface has sediment cover with lower wavespeeds than the basement. Wavespeeds increase steadily through the crust, and then jump again on entry to the upper mantle at the Moho. The wavespeed increases again through the mantle punctuated by the upper-mantle discontinuities near 410 and 600 km depth. Wavespeed gradients are generally not strong, but are sufficient to turn back energy to the surface. When a cartesian geometry is used, earth-flattening transformations increase the gradients to compensate for the decrease in horizontal distance scale with depth due to sphericity (e.g., Chapman, 1973).
The sediment, crustal, and mantle structures are each linked to a set of reverberations whose interference controls the associated surface wave dispersion. But, such multiple reflections with the same slowness (phase speed) can also span multiple zones so that dispersion behaviour is linked. An individual mode will then partake mostly of the character of the multiple guides in different segments with link points where mode branches nearly meet. The presence of low velocity channels within the structures produces further opportunities for guided waves with further classes of modal interaction (Kerry, 1981).
I here bring together a range of results based on the reflection properties of the free surface and the structure at depth to provide physical insight into the nature of surface wave dispersion and the way in which different modes of propagation interact to produce the nature of the different modes. Following a presentation of the nature of coupling between different types of wave guiding, I show how the properties of the full modal field can be readily examined by constructing and displaying the response of the full structure as a function of frequency and phase speed (slowness).

The character of the modal field
Consider an stratified elastic medium with surface P and S wavespeeds α 0 , β 0 , underlain by a uniform halfspace with properties α L , β L . The response of the medium to a surface source in the frequency-slowness (ω − p) domain is (Kennett, 1983): where R 0L D is the response of the entire structure below the surface, R F the matrix of free-surface reflection coefficients, W F the displacement response including free-surface effects and I the identity matrix. Here the source has been left arbitrary. Different aspects of the seismic wave field can be isolated by selecting the components of R 0L D to be modulated by the reverberations in the structure (the inverse term) and the free surface response W F . For isotropic, and transversely isotropic, media the response (Eq. 1) breaks into two separate parts. SH waves propagate independently and the characteristics of the medium are expressed through the SH reflection coefficient [R 0L D ] HH . The P and SV components are coupled so that we need to treat 2×2 matrices of reflection coefficients and free-surface terms where we have used a compressed notation so that (3) Here Z is the vertical and R the radial component. The P wave response of the structure can be emphasised by selecting just the P P and P S elements of the matrix R 0L D in Eq. 1 and the SV response by selecting the SS and SP elements.
The response of the medium specified by Eq. 1 will be singular when the determinant of the reverberation term vanishes, i.e., det I − R 0L D R F = 0. The locus of such values in frequency-slowness space determines the trajectory of the dispersion branches for seismic modes. We will concentrate on the response directly in the frequency-slowness domain. Seismograms may be constructed by introducing the source spectrum and relevant horizontal phase terms followed by integration over both slowness and frequency.
The dispersion relation for the entire structure takes a relatively simple form for Love waves, involving only the SH (tangential) component: since SH waves do not couple to P and the free-surface reflection coefficient for SH waves is unity. This means that the dispersion condition requires the SH reflection coefficient for the entire structure has to be equal to unity and so its phase χ D (p, ω) must be a multiple of 2π: This property can be used to set up an efficient recursive system for models without localised low velocity zones (Kennett and Clarke, 1983).
For the P-SV case the dispersion relation for Rayleigh waves takes the form where R F (p) is the free-surface reflection matrix that depends on slowness p alone and R 0L D (p, ω) is the reflection matrix from the full stratification. Other equivalent representations exist where the structure is broken at different levels (Kerry, 1981;Chen, 1993).
In terms of the components of the free surface reflection matrix R F and the reflection matrix R 0L D from the structure beneath the surface, the dispersion relation (Eq. 6) can be written as With the aid of Eq. 7 we can identify different propagation regimes in terms of slowness p or its reciprocal phase speed c = 1/p.

Both P and S evanescent
In this regime there is only a single possible Rayleigh mode, the fundamental. At moderate to high frequencies only the coefficient R SS D has significant amplitude and R P P D , R P S D , Rd SP are negligible. The dispersion relation for the fundamental mode (Eq. 7) then reduces to so that the small amplitude of R SS D has to be compensated by the growth of R SS F (p) in this evanescent regime. With increasing frequency R SS D tends to zero, and the dispersion relation for the fundamental mode tends to that for a Rayleigh wave on a half space with the surface properties (see, e.g., Kennett, 2001).

S propagating, P evanescent
Now, in principle there is an infinite sequence of Rayleigh modes, and the higher modes have the asymptotic limit p = β −1 0 , so their phase speed is always greater than the surface S wavespeed β 0 . For this slowness range S waves have travelling wave character at the surface and a turning level at depth. P waves are still evanescent throughout the structure. The spacing of the modal branches at fixed slownesss p (phase speed c) is primarily controlled by the inverse of the delay time τ (p), so that the interval between successive modes decreases as the slowness decreases and phase speed increases (see, e.g., Kennett, 1983).
For high frequency propagation, the decay of evanescent P means that R P P D , R SP D , R P S D are small compared to R SS D . We neglect all terms involving R P P D , and all processes involving more than a single conversion of wavetype in the structure at depth. Then, at moderate frequencies we can reduce the dispersion relation (Eq. 7) to the approximate form The single conversions from depth will be weak because of the evanescence of the P wave legs. At high frequencies these conversion terms in braces in Eq. 9 can be neglected.
Equation 9 applies to the full set of Rayleigh modes, so that for fixed slowness p there will be a a sequence of modes with increasing frequency and mode number, with similar spacing in frequency dictated by the phase of the R SS D R SS F combination.

P and S both trapped
In this regime both P and S waves are trapped in the structure, and the full form of the dispersion relation (Eq. 7) needs to be used. The principal term comes from the product of two elements that correspond to dominant S or P propagation. The condition for the first bracket to vanish is a direct continuation of the surface wave branches from Eq. 9 corresponding to trapped S waves. The vanishing of the second bracket represents a suite of modes associated with multiple P-wave reverberations between the surface and structure at depth, which will be most evident when there is a rapid increase in P wavespeed in the near surface. Thus, the dispersion relation (Eq. 7) can be read as the interaction of two separate dispersion systems for SV and P waves linked by the conversion terms with a double surface interaction. The conversions R P S D , R SP D will have a strong dependence on the V p /V s ratio through the structure.
The modes with dominantly S character continue from the previous regime. For phase speeds c just greater than the surface P wave speed α 0 the freesurface reflection coefficient R SS F is quite small whilst R SP F is close to unity. This means that the near-surface conversions between P and S play an important role in controlling the modal dispersion.
Through this regime the main Rayleigh mode branches are controlled by the influence of S-wave structure. Even when P dominated branches are not directly evident their existence modulates the structure of the other mode branches by influencing their spacing (see, e.g., Figure 11.3 of Kennett, 1983;Sun et al., 2021). As noted above the spacing of branches at a given slowness is approximately proportional to the inverse of the delay time. This means that the spacing for the P-dominated modes with small delay is much larger than that for the S-dominated modes.
The full dispersion relation (Eq. 7) can be thought of as comprising the product of terms that would represent specifically S or P dispersion coupled by interconversions between wavetypes. A mode branch will switch character according to which dispersion relation is most closely satisfied. There will be a rapid switch in slope near where both dispersion relations approach zero. At such zones mode branches nearly touch ('osculation' points).

P propagating, S not trapped
Once the phase speed is larger than the shear wavespeed at the base of the structure, S waves leak into the underlying halfspace by radiation, and the S-dominated branches head into the leaking mode regime on lower sheets in the complex slowness plane at each frequency. Just beyond the cut off at c = β L such modes can still have a muted effect on the top sheet from the proximity of the poles. Even so partial trapping is possible for P waves and the the response is singular when but this only occurs for complex p at real frequency. The dominant conversions will occur at the free surface. For a smooth wavespeed profile with gentle gradients with depth, internal conversion from P to S will be negligible, so the dispersion relation reduces to a direct analogue of the condition for acoustic modes in a layered fluid but still |R P P F |< 0 and so there are no real p roots for real frequency. With strong wavespeed gradients, such as in the near-surface zone in the presence of sediments, inter-conversions between P and S become significant and such P-dominated modes still have a dependence on the S wavespeed structure.

P and S not trapped
Now the phase speed exceeds the P wavespeed at the base of the layering so that P waves radiate into the underlying half space and a modal treatment is no longer very useful. Only leaking modes are left well onto the underlying complex slowness sheets that have little influence on the real frequency-slowness response.

Interacting waveguides
The dependence of the seismic response for a surface source on structure at depth comes through the combination in Equation 1. In the case of a structure containing two distinct parts separated by a level J we can expand the reflection term R 0L D into contributions associated with shallow (0J )and deep structure (JL): We can think of these two parts of the structure as representing separate waveguides linked to surface reflections, with waves always having to pass through the shallower zone to get to the deeper. In consequence, there will be interaction between the surface wave propagation in the two waveguides.
With the split (Eq. 15) of the reflection response, the matrix inverse in Eq. 14 can be expanded as so that each reflection from the deeper part of the structure is accompanied by multiple reverberations in the shallow structure. Although Eq. 16 provides some insight into the wave processes at work it is not directly useful for assessing dispersion.
Using the decomposition (Eq. 15), the dispersion relation (Eq. 6) takes the form We can rearrange this singular term as (18) where we extract the product of two components that have the form of dispersion elements for the shallower and deeper parts of the structure, together with a coupling term that involves both elements of the reflection response and the free surface reflections. Thus if we identify A = R 0J D R F and B =R JL D R F our dispersion relation can be cast into the form where the coupling term C J involves a minimum of two surface reflection terms and propagation in both parts of the structure. This is a comparable scenario to that encountered above for the coupling of SV -and Pdominated modes. Consider then the general coupled dispersion system Unless det(I − B) itself approaches zero, the frequency shift δω will be small. Hence the dispersion characteristics will be close to that for the A system alone. A comparable effect arises near frequencies that satisfy det(I − B) = 0, but now involving the frequency derivative for the B relation. For such a coupling scenario we therefore have two sets of largely independent mode segments satisfying the separate dispersion relations, with strong interaction only where both dispersion relation are close to being satisfied.

Figure 1
Schematic representation of multiple modes associated with a shallow low speed waveguide (red), crossing the modal curves for a deeper guide with higher wavespeeds (blue). The spacing of modes is inversely proportional to delay time and so larger for the lower wavespeed guide.
For our situation with a separation into shallower and deeper parts of the structure, we therefore have a mode branch suite with relatively wide spacing associated with the shallow component where the delay time is small, these will appear to cross the modes linked to the deeper structure with longer delay times and closer spacing. This situation is schematically rendered in Figure 1. There are no actual crossing points, though the various branches come close to touching. In the neighbourhood of these osculation points there is strong curvature as an individual mode branch switches character between shallow and steeper slopes. Even with just a separation into two parts we see that a formal description of the behaviour becomes quite complex, but we are able to recognise the main characteristics. With the addition of additional component of structure with depth, we will get a similar situation with suites of apparently distinct modes associated with each part of the structure that interact strongly locally.
For Rayleigh waves, the decomposition of the dispersion relation for each part of the structure will have a strong variation with slowness, as discussed above, because of the presence of both Sand Pdominated modes. Thus, superimposed on the basic mode structure controlled primarily by the S wavespeed distribution, we have further effects linked to the (partial) trap-ping of P waves.
For Love waves in an isotropic model, the situation is simpler because we no longer have wavetype coupling and there is no slowness variation of the free-surface reflection coefficient. Nevertheless, there is still the potential for interaction between waves dominantly propagating in different parts of the stratification.
Each discontinuity in wavespeed structure or zone of strong wavespeed gradients will have a modest influence on the modal structure, but the largest effects are associated with significant jumps in physical properties, as in models comprised of a few uniform layers. In frequency-phase speed space the presence of a discontinuity induces an inflection in the slope of the modal trajectory around the wavespeed associated with the faster layer, even when modal branches do not approach closely. In effect, at each discontinuity we are seeing the interaction of the waveguides above and below the interface.
Strong effects on modal dispersion arise in the presence of pronounced zones of lowered seismic wavespeeds within the structure that act as strong waveguides, particularly at high frequencies. Channelled waves can be trapped in such waveguides, with only weak linkage to the external structure via evanescent waves (Kerry, 1981). The transition between the regular surface wave modes and the channel waves is very abrupt as a function of slowness (or phase speed) and can cause substantial difficulties in tracking individual mode branches.

Illustration of multi-mode surface waves and leaking modes
In the discussion above I have focussed on the dispersion of surface waves and the way that this is affected by coupling between different components of the waveguide and between wavetypes. The tracking of dispersion for complex models for even a few mode branches can become a tricky exercise in numerical root finding, particularly when it involves leaking modes where the singular points are no longer on the top sheet. However, when seeking to understand the structure of the wavefield it is more effective to look directly at the response in the frequency-slowness or frequencyphase speed domain, since this will capture all the different mode branches and their relative excitation (e.g., Kennett, 2001;Dal Moro, 2020). All modes, including leaking ones, that make an effective contribution to the wavefield can thereby be captured. A number of authors have constructed such diagrams, principally in context of analysing models with a few uniform layers. For example, Li et al. (2022) have constructed 'theoretical dispersion spectra' as an aid to the interpretation of dispersion diagrams constructed from very dense LASSO array in Oklahoma, with attention to the difference in characteristics seen on the vertical and radial components.
A similar approach is adopted here to examine the separate behaviour of the SH and P-SV systems and the way in which the different components of the wavefield separate modes with dominant SV or P character. I con- struct the three components of the seismic response from Eq. 1 for a surface source evaluated on a preassigned grid with even sampling in frequency and slowness or phase speed. The use of phase speed provides greater separation of the different aspects of the modal field and a more direct relation to the wavespeed structure and so is adopted in the figures below.
To illustrate the interaction of different waveguides, I have built a composite model with horizontal stratification combining realistic sedimentary, crustal and upper mantle structures (Figure 2). The shallow structure was derived from an ambient noise study on a dry lake bed to the north of Canberra in southeastern Australia (C. Jiang -personal communication, 2022). The crustal structure is derived from receiver function studies of stations in eastern Australia on Phanerozic fold belts. The mantle structure below 45 km depth comes from ak135 (Kennett et al., 1995) and is continued to 1000 km depth. Earth flattening is applied to compensate for the sphericity of the Earth. The S wavespeed at the base of the flattened model is 7.513 km/s, and only modes with phase speed higher than this will lie in the true leaking mode domain. At smaller phase speeds the extension of the structure to depth allows S waves lost from the crust to be returned to the surface by refraction and so ensure full trapping of energy, described by modes on the top slowness sheet.
This model has strong gradients in seismic wavespeed at the surface with a high ratio of P to S wavespeeds. The rapid changes in seismic wavespeed with depth allow strong trapping of both S and P energy above the basement structure. The wavespeed gradients are more muted in the crust down to 35 km depth, and then a modest jump in wavespeed at the Moho ties into the upper mantle structure. Fine discretisation with 5 m thick uniform layers was used to 0.5 km depth, and then for the crust layers 0.5 km thick were employed. In the mantle below 45 km, layer thickness increases from 10 to 14 km with increasing depth. Weak attenuation has been applied throughout the model with Q −1 α = 0.001 and Q −1 β = 0.002. The presence of attenuation means that no surface wave poles actually reside on the real phase speed axis.
In Figures 3 and 4 I show the frequency-phase speed response for the composite model, with 257 layers, for SV and SH excitation via a grey-tone display. The responses were calculated using recursive construction of the reflection matrix R 0L D for the entire structure for each phase speed and frequency. I have used 200 frequencies and 400 phase speeds for these displays and have set a common threshold for the transition to black. The various modes appear directly through the elevated amplitude associated with the trajectories of the modal dispersion. The first few modal dispersion curves are labelled for Rayleigh waves in Figure 3 and Love waves in Figure 4. In each case the width of the band around each mode provides a measure of the relative excitation.
A comparison of the response of the full model and that for just the sediment zone is made in Appendix A. The dispersion curves for both Rayleigh waves and Love waves out to 5.3 km/s, for the full model are displayed in Appendix B.
There is a notable difference between the SV and SH scenarios. The gradient zones in the sediments have a stronger effect on the Rayleigh modes than the Love modes, in part because of the rather low P wavespeeds at the surface. In each case the multiple mode branches associated with crustal trapping of S waves appear once the phase speed exceeds 3.5 km/s. For the Rayleigh modes the transition between the continuation of the sediment mode branches and the crustal modes is so sharp that the first higher mode branch appears to be dissected into several pieces. For the second higher mode it is more difficult to track the continuation of the Figure 3 The radial component of the surface response in frequency-phase speed for the composite model, corresponding to SV excitation. The radial component emphasises Rayleigh modes for lower phase speeds, but still displays noticeable contributions from P-dominated modes for the interval of crustal P wavespeeds. The markers indicate the wave speeds associated with the different segments of the model (Figure 2) with S structure in green and P in red. The solid symbols mark the wave speeds at the surface and the base of the full model. mode branch but its presence is apparent from the enhanced visibility of the crustal modes. A similar effect is seen for the Love modes, where the crustal modes appear more distinctly in the zone where the extension of the sediment branches occurs, even though the continuation of such branches cannot be tracked directly. These features associated with interaction between waveguides are marked with a circled C in both Figures 3 and 4. Such complications will not be readily picked up by conventional dispersion calculations unless they are accompanied by detailed assessments of relative amplitudes (Sun et al., 2021), but this involves much more work than the direct approach employed here.
With the threshold employed for the plots in Figures 3 and 4, the modal branches associated with mantle propagation are very difficult to discern for either Rayleigh or Love waves. Such modes are much more effectively excited for sources at depth, and can be handled with a similar frequency-phase speed approach with a source set at some depth. The trajectories of the modal branches are unaffected by such a change of reference level, but the amplitudes will be modified.
To the right of Figure 3, for wavespeeds greater than 5.0 km/s, we see a group of mode branches associated with P wave propagation in the crust that extend into the leaking mode zone and even begin to pick up mantle effects, though finer sampling in frequency would be needed to see the details. It is intriguing to find that these P-dominated modes are most visible along the general location of the extension of the SV sedimentary branches. The delay times for P waves in the crust are around 1.8 times smaller than those for S waves so the spacing between the mode branches in frequency is greater.
In the frequency range displayed in Figure 3 there is no direct interaction between Pand SV -dominated modes. To see such effects we need to look at somewhat higher frequencies. In Figure 5 I have confined attention to phase speeds between 2.0 and 6.0 km/s for frequencies up to 10 Hz. This figure compares the vertical component for P excitation with the radial component for SV excitation. The two images of the frequency-phase speed response strikingly separate the P-dominated modes on the vertical component from their SV -dominated counterparts on the radial component.
The complementary character of the P-dominated and SV -dominated modes is well illustrated in Figure  5. The P-dominated trajectories are comprised of segments of many different mode branches. The net effect is to cut across the SV -dominated modes with just short intervals in phase speed where character is indeterminate, as expected from Eq. 22. In the neighbourhood of the close approach between mode branches with different physical character, the modes are visible in both displays. However, away from such locations the modes tend to be more distinct in one or other of Figure 4 The tangential component of the surface response in frequency-phase speed corresponding to SH excitation displaying the Love modes in the sediment and crust. The markers indicate the S wave speeds associated with the different segments of the model (Figure 2). The solid symbols mark the wave speeds at the surface and the base of the full model. the panels. We have noted before the way in which the crustal modes apparently cut the extension of the sedimentary SV branches; a similar effect can be seen for the P mode branches. The blurring of the P-dominated mode branches beyond 4.5 km/s arises from the influence of the upper mantle on the higher branches of the SV -dominated modes.
In this complex near-surface structure the first higher mode has a higher P wave content than the fundamental mode. Such behaviour has been found for other models with strong near-surface gradients and would not be expected from studies that employ relatively thick uniform layers.

Discussion and Conclusions
Modern developments in seismic recording open up the possibility of a full rendering of the wavefield, and this needs to be complemented by holistic methods that can present the full character of the seismic response of a medium. As the higher modes of surface waves become more readily accessible, the full character of seismic structures can be revealed rather than the muted versions derived from just the dispersion of the fundamental mode.
The direct use of frequency-phase speed as presented here enables effective analysis of complex structures, without any of the issues associated with the presence of, e.g., low velocity channels, in dispersion calculations. All modes are present and can be examined in detail by appropriate choice of frequency and phase speed ranges. On even a modest computer a 200×200 frequency-phase speed array can be calculated for a 400 layer model in a few seconds to visualise the full suite of modes and their interactions.
The effect of sediment compaction means that the near-surface frequently contains quite strong wavespeed gradients, and such structures contribute strongly to the trapping of both S and P waves that can be recognised as surface wave branches. For shallow structure it is tempting to make an approximation with a few uniform layers to simplify inversion, but this approach runs the risk of an inadequate representation of the physical processes.
The modulation of SV -dominated dispersion curves by the influence of the P wavespeed structure so clearly seen in Figure 5 has been suggested as a diagnostic for estimating shallow P wavespeeds (Li et al., 2021). For such studies it is important that a full elastic treatment be made, because the acoustic approximation will tend to under-estimate the modal frequencies, though the discrepancy is reduced for very low shear wavespeeds compared to those for P waves (Roth et al., 1998). With considerable effort it is possible to track leaking modes on the other sheets (Shi et al., 2021), but it is preferable to work with an extension of structure to depth so that only propagating modes need be considered. As noted in Appendix A, such an extension also means that the low-frequency parts of the dispersion characteristics can be correctly captured.
The use of the frequency-phase speed diagrams enables an improved understanding of the characteristics of the interactions of surface wave modes due to the presence of waveguides formed by different parts of the structure or the effects of the structure associated with different wavetypes. Where fine sampling of the seismic wavefield has been achieved, observed seismograms can be transformed into the frequency-phase speed domain using, e.g., the F-J transform (Li et al., 2021). It is then possible to contemplate direct compar-ison of the theoretical and observed response in this transform domain, as the basis for a nonlinear inversion using an exploration of a suitable parameter space. Where specific dispersion results are needed, the techniques developed for handling observed spectra (e.g., Dong et al., 2021) can be used to extract the trajectories of mode branches in frequency-phase speed space.
The results we have derived are for a stratified model with no horizontal variation in seismic wavespeeds within each layer, for which the various surface wave modes propagate independently. The influence of lateral heterogeneity can be described by introducing coupling between mode branches, which will have most effect when modal branches are tightly spaced. Thus, although variations in the near surface are stronger than at depth, the broad spacing of the modal branches mean that their general character will be maintained. In general, the strongest features seen in the frequency -phase speed response diagrams in Figure 3 and 4 can be expected to be robust under the influence of moderate horizontal changes in wavespeed.

Appendix A: Comparison of coupled and single waveguide results
In Section 4 I have shown the dispersion characteristics for the composite model, where there are three major divisions of seismic structure that each act as a waveguide. How then does the response compare with the situation where the sediment zone is considered on its own? The sedimentary structure is taken down to 1 km depth (as indicated in the right hand panel in Figure 1) and then extended to a half space with no contrast in physical properties, rather than the further jump as in the case where the full crust and mantle are included.
In Figure 6 I show a comparison of the frequencyphase speed response for the two cases for frequencies up to 5 Hz and phase speeds out to 6 km/s. The modes for the sediment-only case are superimposed in magenta on the grey tone for the full model. The underlying half-space for the sediment-only model now has S wavespeed of 3.0 km/s and P wavespeed of 5.2 km/s, so that everything with phase speeds greater than 3 km/s represents leaking modes for this case. However, in the full model this regime will still contain trapped modes.
For phase speeds below 2.2 km/s there is a complete correspondence between the Rayleigh modes for the sediment-only and full cases. This regime represents trapping by the strong near-surface S wavespeed gradients. The change in structure at depth between the two models manifests itself in a change in the character of the fundamental mode and the first higher mode with frequency. The first higher mode shows a distinct sensitivity to P wave structure and extends well into the leaking mode domain. There is slight P dependence for the second higher mode that only just nudges into the leaking mode field. The other higher modes for the sediment-only case effectively truncate at the 3.0 km/s cut off associated with the underlying half space. It is interesting to note that the extension of these branches in the full model to the entry into the main crustal waveguide at 3.5 km/s shows weak P wavespeed dependence.
In the frequency range up to 5 Hz, only the lowest frequency mode for P-dominated propagation is visible. For the phase speed range less than 4.5 km/s where the P waves interact with the wavespeed gradient in the sediment, the leaking mode for the sediment-only case and the trapped mode for the full model coincide apart from the interruptions by the crustal modes. The behaviour near the P cut off for the sediment only case, as near the S cut off, is modified by the absence of the deeper structure.
The changes in the dispersion behaviour imposed by the truncation of the model in the sediment-only case present a warning for sensitivity to the assumptions made about the extent of the model. As noted earlier it is preferable for models to be allowed to extend to substantial depth, so that the evanescent behaviour of the modes is well represented, rather than to force a specific, shallow, depth for the range of the structure. Inversion of the true dispersion curves with a truncated model is liable to produce a distortion of structure, particularly if only the first couple of modes are available. The dispersion curves for slower phase speeds should be reliable, and discrepancies will only occur as the assumed S wavespeed for the uniform underlying half space is approached.
In many studies the choice of the depth at which a transition is made into a uniform half space is made based on a priori assumptions about the nature of structure. Thus it would be common to make a truncation at the expected base of sediments and often to use a model composed of a few uniform layers. Strong gradient zones, such as in compacting sediments, require fine layering for an accurate description of the wave behaviour and there is the potential for an over-simplified model to be very misleading. Complications will ensue if there is a significant interface in the actual scenario below the level of truncation. This can be seen in Figure 6 in the behaviour of the fundamental and first higher Rayleigh modes.

Appendix B: Dispersion Curves
In Section 4 we have demonstrated the interaction of multiple wave guide effects by examining the re-  sponse in frequency -phase speed space, rather than the more conventional approach of working with dispersion curves. In Figure 7 we show the dispersion curves for both Rayleigh and Love wave modes for phase speeds out to 5.3 km/s, for the same frequency band as in Figures 3 and 4. A slightly simplified model has been used with 100 layers extending down to 410 km depth. For higher phase speeds more than 200 modes occur in the frequency interval up to 2.25 Hz and it is no longer at all easy to track features that depend on subtle modulation of the spacing of the modal branches.
This portion of frequency -phase speed space was computed for Figures 3 and 4 with 200×200 evaluations of the response, i.e., 40,000. Whereas, for the dispersion behaviour the roots of the secular equation have to be found for each mode. The dispersion results in Figure 7 were computed with the scheme of Kerry (1981) which uses fixed phase speed and then searches for successive modes in frequency. The choice of phase speed step is adaptive so that the changes in slope of the mode branches can be followed. More than 50,000 roots need to be found for the Love wave diagram and over 200,000 for the Rayleigh case. For each root the frequency has to be bracketed and refined. This means that many response evaluations are required and so computation times are typically several orders of magnitude greater than for the direct approach used in Figures 3 and 4. Very high precision is needed to disentangle the modal branches and their changes of slope.
We can see a direct correspondence between the visible features in Figures 3 and 4 and the structure of the dispersion branches for the span of phase speeds covered in Figure 7. Where apparent extensions of modal branches from the sediments to higher wavespeeds occur we can see that they are related to the fine structure of the dispersion curves and, notably, tight osculation points. This is very distinctive for the projection of the mode R1 through the crustal and mantle branches. The switches between branches occur on such a fine scale that some roots were missed or misidentified, producing a small blank zone for phase speeds between 4.6 and 4.7 km/s associated with the depth range between 120 and 210 km in the mantle.