Earth Structure Across Many Scales

Over my research career I have worked on many aspects of Earth Structure on scales from the whole globe to the very local, this has provided insights into the nature of heterogeneity and the importance of features with small length scales (up to tens of kilometers) revealed with higher‐frequency seismic waves (>1 Hz). Progress in determining 3‐D structure depends on representing the dominant variation with depth and large‐scale structure as a starting point. Intermediate and fine scales of heterogeneity (less than 100 km), below the typical resolution of seismic tomography, play an important role in shaping the seismic wavefield. The larger scales of variation in seismic structure can be directly linked to geodynamics through the influence of temperature and composition. The smaller scales retain a memory of past processes that can link to geochemical studies.


Plain Language Summary
Heterogeneity is ubiquitous in the Earth's interior on all scales from the planetary to the local to individual minerals. Seismology has made a major contribution to understanding the internal structure of the Earth, particularly through seismic imaging. Nevertheless, it is important to make use of the full character of seismic recordings in understanding the Earth and high frequency content should not be excluded. Although they make seismograms look "messy", the higher frequencies contain important information about smaller scale structures. In understanding the Earth, the full range of heterogeneity needs to be considered. Enforcing apparent simplicity on seismograms by imposing low pass filtering can give misleading results. Multi-scale structures exist throughout the Earth but can usually only be detected indirectly. For the lithosphere dense arrays allow imaging of scales down to a few kilometers superimposed on the larger scales resolved by tomography. I show how combining information from many different aspects of the seismic wavefield can provide a new view on the nature of structures; larger scales relate to geodynamics the smaller to geochemistry. As a result, it is possible to use seismology as a link across a broad range of disciplines.
wave propagation (Kennett, 1984). These approaches have proved very valuable for understanding the interaction of the seismic wavefield with complex structures.
Photographic recording of seismic waves utilized two separate frequency bands, long-period and short-period lying on either side of the microseism peak (∼0.125 Hz). The advent of high-fidelity digital systems with a broad frequency response has meant that the full range of seismic ground motion can be utilized. This improvement has allowed the exploitation of ambient seismic noise using correlation techniques, which has reduced dependence on the specific locations of earthquake sources. Nevertheless, common practice in analyzing broadband seismic records is to apply a low-pass filter. The action of the filter simplifies the appearance of the seismograms, enhancing distinct seismic phases by suppressing high-frequency clutter. Such simplification of seismograms has yielded important results on broad-scale seismic structure but comes at the cost of suppressing information on finer scales.
On arrival at the Australian National University in 1984 I was struck by the unusual character of the records of southern Banda Sea events recorded at the Warramunga Seismic Array in northern Australia. A distinct medium frequency onset was followed a couple of seconds later by the emergence of high frequencies that persisted for minutes. A similar pattern was seen for S superimposed on the P coda. Although the records could be dramatically simplified by low-pass filtering (retaining only frequencies below 0.25 Hz), the high-frequency contribution was too prominent to be ignored. Such characteristics proved also to be typical of propagation across the Australian cratons ( Figure 1) and must be associated with relatively small-scale structures with scales of a few kilometers. Changes in the appearance of the seismograms with source depth hinted at the presence of some class of structural transition in the mid-lithosphere (Kennett, 1987).
In the late 1980s we used relatively dense seismic arrays of portable seismometers (for the time) in northern Australia using events in the Indonesian subduction zone to study upper mantle structure. Evidence then began to Figure 1. Three-component seismograms for an event in the Flores arc recorded at a portable broad-band station in the Northern Territory of Australia at an epicentral distance of 12.60°. All three components show extended high frequency coda for both P and S waves, starting a few seconds after the lower frequency phase onset, indicative of a strong scattering environment. The phases are marked for a regional model, 7 s faster for S than the ak135 reference model (Kennett et al., 1995). In the lower panel a frequency-time analysis is performed for each component to show the slow decay of the high frequency energy for both P and S. emerge for structural variability on horizontal scales of 100-300 km on top of larger scale geographic patterns. The vertical scale needed to be shorter than the horizontal (Kennett & Bowman, 1990). Thus, it was appropriate to think of a complex Earth with multiple scales of variation revealed by different ways of sampling structure though the prevailing picture of the time was that of purely radial variations.

Global to Regional Seismic Tomography
A prerequisite for the definition of 3-D structure is an effective treatment of the main radial variations due to increasing pressure with depth. For low-frequency studies the PREM model (Dziewonski & Anderson, 1981) is still the most used reference because it provides a good match to the frequencies of the normal modes of the Earth and hence long-period seismic waveforms. The effects of 3-D structure can be included via perturbation theory using a modal representation (e.g., Dahlen & Tromp, 1998), with comparisons made with waveforms (e.g., Li & Romanowicz, 1995;Woodhouse & Dziewonski, 1984), phase dispersion of surface waves (e.g., Laske & Masters, 1996) or secondary variables such as structure factors (e.g., He & Tromp, 1996). With the aid of such techniques a relatively consistent model has been established for the largest scale variations in Earth's mantle (spherical harmonic order less than 6), but much greater variability is found for higher harmonics depending on the data sets employed and the style of analysis (e.g., Becker & Boschi, 2002;Ritsema & Lekić, 2020).
With modal analysis the prime emphasis has been on shear-wave structure for which sensitivity is largest (Dahlen & Tromp, 1998). Extensive developments in computational seismology now allow direct calculations of seismic waveform for 3-D structures, with inversion made possible with efficient methods for calculating derivatives with respect to structure (e.g., French & Romanowicz, 2015;Lei et al., 2020). This approach offers the potential for rendering P wave structure and density alongside shear wavespeed. However, much interpretation of 3-D structure is commonly based on just S wavespeed variations (e.g., French & Romanowicz, 2015), where lowered wavespeed are assumed to be linked to elevated temperatures.
For higher-frequency studies the accuracy of source location becomes more important. Although the travel-time tables of Jeffreys and Bullen (1940) had allowed a good definition of global seismicity, there were issues of consistency between different seismic phases and in the assessment of event depth. An international effort to produce improved results employed a velocity model directly rather than creating tables. The first product iasp91 was focused on the main P and S phases (Kennett & Engdahl, 1991) and with its aid travel-time tomography immediately produced striking images of subduction zones (e.g., van der Hist et al., 1991) with much improved variance reduction in the fit to arrival times. A further effort to produce an effective representation of all the major seismic phases produced the ak135 model (Kennett et al., 1995) that is currently in use for event location by all the major international agencies. The model ak135 has been extensively used in studies of 3-D mantle structure using the arrival times of seismic waves (e.g., van der Hilst et al., 1997;Widiyantoro & van der Hilst, 1996;Zhao, 2015).
Nevertheless, careful studies of seismic waves reflected at the underside of the core-mantle boundary (e.g., Kaneshima & Helffrich, 2013) indicated the need for a revision of structure at the top of the core. The sharp reduction in P wavespeed between the mantle and core means that the outer part of the core is only fully sampled by phases such as SKS and so the shear structure in the mantle influences models of P wavespeed in the core. It is not possible to just change the seismic wavespeed in one part of the model without upsetting the fit to many different phases. A revised model ek137 (Kennett, 2020) achieves the difficult task of retaining the fit to mantle phases whilst improving the rendering of the core. Changes from ak135 are small, with an increase in P wavespeed in the top 1,000 km of the core and slight changes to lower mantle structure. An independent study using the cross-correlation of the late coda of large earthquakes (Ma & Tkalčić, 2021) that focusses on seismic waves traveling steeply through the core favors a core structure very close to ek137.
Using the ak135 model Engdahl et al. (1998) set up a procedure to reassess the phase association for all reported phase arrivals followed by relocation to give improved event depth and position. The resulting EHB catalog of selected events provides a high-quality data set for tomographic studies, which has allowed good definition of both medium-scale structure in the main mantle and subduction related features. With the selected events it proves possible to exploit not just P wave picks but also S waves. Although the quality of S picking is lower the larger effects of 3-D structure compensate and good imaging can be achieved, notably of deep subduction. By exploiting stations where both P and S are reported, joint inversion of P and S arrival times can be made in terms 4 of 8 of mapping the distribution of bulk-sound speeds (depending just on bulk modulus) and shear wavespeed (e.g., Gorbatov & Kennett, 2003;Kennett et al., 1998). The resulting 3-D images for the mantle allow the identification of different classes of heterogeneity regimes with depth, based on the relative behavior of the variations in bulk and shear modulus (Kennett, 2004). With only a single tomographic image it is difficult to disentangle the influence of temperature and composition (e.g., Davies et al., 2012), but with multiple images of different physical properties there is more information on which to base a geodynamic interpretation utilizing the character of the bulk and shear modulus variations revealed by mineral physics experiments at high pressure.
Extensive analysis of longer-period arrivals and inversions incorporating allowance for finite frequency effects have confirmed the major features uncovered with the use of bulletin results (e.g., Houser et al., 2008;Sigloch et al., 2008), and offer the prospect of further improvement exploiting the extensive deployments of high-quality temporary stations across the globe. The impact of the various styles of seismic tomography on geodynamic studies has been considerable, with the imaging of subduction related structures extending well into the lower mantle, and complex structures near the core-mantle boundary (e.g., Ritsema & Lekić, 2020). Much effort has been expended on relating large-scale seismic structure to geodynamic processes (e.g., Forte et al., 2010).
At depth the primary influences on seismic wavespeed come from temperature and major element composition and will control the major features such as subducting plates and rising plumes. Variations in the minor elements, which play such an important role in geochemical interpretation, have much less effect, but may well make a significant contribution to the low amplitude smaller-scale features seen in high-resolution tomography of the mantle. It should always be remembered that there will be structural features on too small a scale to be imaged by seismic tomography. Indeed, apparent gaps in subduction attributed to a slab tear may in many cases represent thinning of the slab to a level that the perturbation in a tomographic cell is too small to be imaged (Kennett & Furumura, 2010).
The large amplitude multiple S and surface waves on the seismograms from shallower events provide a major tool for investigating the outer parts of the globe. Although the fundamental modes are restricted in their depth penetration, the addition of higher mode information gives penetration through the upper mantle (e.g., Debayle & Ricard, 2012;Nolet, 1990;Schaeffer & Lebedev, 2013;Visser et al., 2008). On a regional scale, the broad-scale patterns are consistent between different styles of analysis (e.g., Kennett et al., 2013) but significant differences in detail can arise depending, in part, on the specific station distribution used.
Regional structure can also be addressed using the patterns of delay times for teleseismic arrivals across networks of seismometers, with the best results achieved with a fully nonlinear inversion with updating of ray paths at each iteration (e.g., Rawlinson et al., 2015). The steep arrival angle of the incoming waves leads to some degree of vertical smearing but horizontal resolution approaching the interstation separation can be achieved. Starting from a smooth reference derived from surface wave studies the results from body wave inversion show widely distributed intermediate-scale structures irrespective of the geological domains on which the stations sit. Figure 2 illustrates the situation in southeastern Australia (Rawlinson et al., 2015) where 650 stations with an approximately 50 km spacing were employed with over 1,000 teleseismic sources. At 120 km depth significant variability is seen in both the cratonic domains and beneath the Phanerozoic fold belts in the east. The scale of the variations is too small to be able to be imaged with surface wave studies at this depth.

Guided High-Frequency Waves
During a visit to the Earthquake Research Institute, University of Tokyo, in 2003 Takashi Furumura raised with me the question of the origin of anomalously high amplitudes of ground motion on the east coast of Japan from deep earthquakes in the subducting Pacific plate. When we looked at the waveforms, they showed a similar characteristic to those seen in Figure 1 with a lower frequency onset followed by a long high-frequency train. I therefore suggest we look at 2-D models with stochastic heterogeneity described by a few parameters superimposed on the thermally imposed structure of the subduction zone. Almost the first model we tried with a correlation length of 20 km along the slab and 0.5 km transverse produced finite difference seismograms that had a strong resemblance to the observations. We therefore investigated a wide class of different structures but found little further improvement (Furumura & Kennett, 2005). An effort was made to look at 3-D, but the computational limitations of the time meant that the frequency range was rather restricted, and we were not able to reproduce the full effects seen in the observations. It was also clear that the details of the configuration of the subducted plate were important in determining the distribution of ground motion at the surface. The stochastic waveguide model provides a general description of the guiding of high-frequency waves along the subducted plate. The localized patches of lowered wavespeed oriented along the plate serve to keep energy trapped against its natural inclination to escape from the high wavespeed plate. Trapping of waves can also occur in the former oceanic crust (Garth & Rietbrock, 2014), and there may also be complexity of structure near the bend in the plate with faulting and hydration. Yet, the stochastic waveguide model can reproduce the waveforms from deep events recorded at ocean bottom seismometers well away from the subduction zone (Shito et al., 2013). With the recent improvements in computation power, 3-D simulations can be made to high frequencies incorporating 3-D heterogeneity structure (Furumura & Kennett, 2021) that come closer to matching the full character of the guided waves.
For events along the Indonesian subduction zone recorded in northern Australia there are systematic variations in the character of the high-frequency guided waves depending on the position along the arc. These effects could be related to change in structure along the subduction zone as the segment of oceanic lithosphere diminishes to the east as the Australian continental lithosphere impinges on the arc in the southern Banda Sea (Kennett & Furumura, 2008).

Continental Lithosphere
The crustal component of the continental lithosphere shows considerable variability reflecting the complex evolution of the modern crust, but such variation extends through the entire lithosphere. The development of receiver-based methods has helped to reveal complex patterns of internal discontinuities. A favored approach has been to exploit conversions from S to P where the conversions arrive before the main S phase and so avoid contamination from crustal multiples. Such Sp receiver function depends on transmitted P and the frequency range is limited by that of the S arrivals. Analysis of Sp receiver functions provides evidence for a significant reduction in S wavespeed in the middle of the lithosphere (e.g., Selway et al., 2015) though the frequencies employed are low and so vertical resolution is rarely better than 15 km. With dense networks of temporary stations, such as USArray, significant geographic variations in the character of lithospheric variations can be mapped (e.g., Lekić & Fischer, 2014;Liu & Gao, 2018).
Is it possible to use P waves instead at higher frequency to provide greater resolution of the crust and uppermost mantle? Fortunately, a route is provided by correlation techniques (Gorbatov et al., 2013;Ruigrok et al., 2011). The autocorrelation of continuous seismic records at a single station can provide an estimate of the P wave reflectivity beneath a station combining contributions from ambient noise and distant events. With the aid of this approach, I was able to analyze the full suite of available stations across Australia (Kennett, 2015) and demonstrated significant variations in the character of lithospheric reflectivity that correlated with the age of the continental crust on which the stations were placed. Working with relatively high frequencies the reflectivity shows small amplitude rapid vertical variations, and transitions between different domains are marked by a change in the spatial frequency. Thus, there is a significant drop in the amplitude of reflections in the mantle just below the crust and a change to slower variation with depth. Yet, there are also indications of changes with depth in the mantle may be identified with mid-lithospheric discontinuities. Such results are not very convincing on single station records, but by applying autocorrelation to teleseismic P arrivals across a suite of stations, with stacking, it is possible to generate an image of P wave reflectivity through the lithosphere (Liang & Kennett, 2020;Sun & Kennett, 2020) with scales of variation around 20 km horizontally and a few km vertically, superimposed on broader-scale trends. Figure 3 from a profile of 63 stations at 3 km spacing in southern central Australia illustrates the presence of fine-scale features throughout the lithosphere at scales that cannot be resolved by tomography. The nature of the imaged heterogeneity is consistent with that inferred from the interpretation and modeling of the high-frequency waves that have traveled nearly horizontally through the lithosphere, as in Figure 1 (Kennett et al., 2017). An effective representation of the P and S coda requires allowance for the large-scale features imaged by surface wave tomography, but also intermediate-scale structure of the type seen in Figure 2 in addition to the finer scales. The net effect is close to a self-similar wavenumber spectrum.
Thus, we have both direct and indirect evidence for fine scale structure permeating the lithosphere. There are indeed indications for mid-lithospheric features, but their character is rather different from that inferred from the low frequency Sp receiver function studies. Rapid small variations in wavespeed with depth set up complex interference patterns between reverberations within the structure, and the image seen comes from their resultant. Changes in the character of reflectivity are recognized to produce apparent features in crustal reflection profiling without a significant contrast in material properties. A similar effect seems to be operating here. Structures that fit a low-passed version of a seismic record may not actually represent a simplified depth regime because interferences mask the true structure.
The fine scale structure in the lithosphere revealed with higher frequency probes ties well to geochemical profiles with depth where these are available. Kennett et al. (2017) make a comparison of stacked P wave reflectivity near xenolith sites in southeastern Australia analyzed by Gaul et al. (2003) and find a strong correlation between inferred changes in chemistry with depth and seismic reflectivity. The correspondence is the more striking because very different approaches are needed to translate information into depth.
For low frequency waves the net effect of the class of smaller-scale structures with more rapid vertical than horizontal variations will be an effectively anisotropic medium. Frequently it is necessary to use some stochastic representation of smaller scale structure to capture behavior with a few parameters. Such simplified representations are useful, but can give a misleading impression of the true situation.

A Multi-Scale Earth
I have been fortunate to be able to work on many different aspects of seismology over my career and have found considerable benefit in having started from a strong theoretical background and then moving into observational studies. Since suitable analysis tools were not available, they needed to be developed from scratch. I have retained this habit of working things out for myself when looking at new fields and it has often got me into trouble since I am accused of neglecting past work! The conventional tools in seismology have aimed to enhance signal-to-noise and so deliver reliable results. In consequence there has been a tendency to immediately impose a low-pass filter and thereby eliminate the untidy clutter of "messy" high-frequency contributions. As I have tried to demonstrate this has the potential of eliminating important information about the Earth, and in particular some of the links between seismic structure and geochemistry.
Heterogeneity is ubiquitous in the Earth's interior on all scales from the planetary to the local to individual minerals. Different scales tend to be illuminated by different modes of analysis, but all need to be considered to provide a full picture of the nature of the Earth. Smaller scale structures sit on top of the large-scale features commonly interpreted in terms of geodynamic processes, but there can be mutual interaction. For the lithosphere we have a chance to examine a wide range of structural scales directly, but in other boundary zones in the Earth, Figure 3. Image of P wave reflectivity through the full lithosphere derived from the autocorrelation of teleseismic arrivals across a dense array of stations in central Australia (Liang & Kennett, 2020). The position of the Moho is inferred from the change in character of the reflectivity combined with receiver function results. The dashed line indicates a transition zone that tapers out at the position of the edge of the Gawler craton inferred from potential field information (indicated by the orange markers). Organized reflectivity is seen throughout the lithosphere with indications of change near 90 km depth that is probably a mid-lithosphere discontinuity. The band LAT indicates the span of the lithosphere-asthenosphere transition inferred from surface wave tomography (Yoshizawa, 2014).