Quasi-Love wave scattering reveals tectonic history of Australia and its margins reflected by mantle anisotropy

The Australian continental crust preserves a rich geological history, but it is unclear to what extent this history is expressed deeper within the mantle. Here an investigation of Quasi-Love waves is performed to detect scattering of seismic surface waves at mantle depths (between 100–200 km) by lateral gradients in seismic anisotropy. Across Australasia 275 new observations of Quasi-Love waves are presented. The inferred scattering source and lateral anisotropic gradients are preferentially located either near the passive continental margins, or near the boundaries of major geological provinces within Australia. Pervasive fossilized lithospheric anisotropy within the continental interior is implied, on a scale that mirrors the crustal geology at the surface, and a strong lithosphere that has preserved this signal over billions of years. Along the continental margins, lateral anisotropic gradients may indicate either the edge of the thick continental lithosphere, or small-scale dynamic processes in the asthenosphere below. Quasi-Love wave scattering is used to identify gradients in seismic anisotropy in the mantle beneath Australia that may represent major lithospheric boundaries between geological provinces and modified mantle flow along the ocean-continent transition.

T he Australian continent preserves a vast tectonic history, stretching from the Archean to the present day 1,2 . By the Neoproterozoic, supercontinent cycles had assembled the three major Australian cratons (North, South, and West), while later continental accretion added the eastern orogens in the Phanerozoic 3 . Since the break-up of eastern Gondwana, and the formation of its rifted passive margins, the Australian continent has moved northwards colliding with SE Asia 4,5 . Today, Australia is the fastest moving continent on Earth 6 , exerting considerable shear on the underlying mantle asthenosphere. Such deformation is thought to cause a lattice preferred orientation (LPO) in olivine generating strong azimuthal anisotropy aligned with plate motion 7,8 . Intriguingly, this strong azimuthal anisotropy is not detected by shear-wave splitting studies, with weak and complex splitting generally observed across Australia [9][10][11][12] . Instead a contrasting contribution from fossilized anisotropy in the lithosphere has been proposed, frozen-in by past deformational events 10,13 .
Observations of Love-to-Rayleigh surface wave scattering, termed Quasi-Love (QL) waves 14 , can be used to pinpoint the location of lateral gradients in seismic anisotropy in the upper mantle ( Fig. 1), such as those that may exist along present (and/or past) plate boundaries. When a Love wave encounters such a boundary, a portion of its energy is scattered from Love into Rayleigh wave motion, generating a QL wave (Fig. 1). The QL wave retains a similar waveform shape as the fundamental Love wave (G1), hence Quasi-Love, but travels with the slower velocity of the fundamental Rayleigh wave (R1). This principle can be exploited to determine the distance (δx) the QL wave has travelled and thus estimate the location of the anisotropic gradient. An anisotropic medium may also cause a small rotation of surface wave polarization [15][16][17] . These polarization-distortions are sometimes referred to in the literature as quasi-Love waves also 18 , however, their characteristics are distinguishable from QL wave scattering as described here (see 'Methods' for further details). Previously QL wave scattering has been detected across a wide range of tectonic settings from continental collision 19,20 , seafloor spreading 21 , subduction 14,[22][23][24] , hotspot volcanism 25,26 , and former continental rift zones 27 . Most recently QL waves have been recorded along the continental margin of eastern North America 28 , implying localized dynamic processes and a disruption of the underlying mantle flow-field at the margin, potentially induced by the transition from continental to oceanic lithosphere.
As of yet, however, there have been no such studies targeting the rich history of the Australian continent and its extensive passive margins.
Here, in this study, analysis of QL wave scattering is conducted for the Australian plate, creating one of the largest and most geographically widespread such datasets. Across the region QL scatterers are found to be preferentially located along the boundaries of the major geological provinces within Australia, as well as, along the passive continental margins. Such results suggest the presence of lateral gradients in seismic anisotropy at upper mantle depths both within the continental lithosphere along old tectonic boundaries and along the edges of the continent where continental rifting previously took place.

Results and discussion
Quasi-Love wave detections across Australasia. Using seismic stations from the Australian National Seismograph Network, and shallow earthquakes (<50 km depth) at 70°-180°epicentral distances, robust QL wave detections were measured on 22% of the 2248 available event-station pairs (Fig. 2, Fig. S7, Supplementary Data 1) with 275 QL scatterers falling within the study area, as shown in Fig. 3. The QL waves detected are of varying amplitude (represented in Fig. 3 by their symbol size) which is not only a function of the strength of the anisotropic gradient but also the geometry of the anisotropy relative to the angle of the raypath. The back-projected scatterers are widely distributed, comprising one of the largest and most geographically diverse QL datasets produced to date, as well as the first results within Australasia. While the QL scatterers are widespread and numerous (the reliability of any single measurement should be treated with caution), when considered collectively, they display intriguing relationships with certain tectonic features at the Earth's surface (Figs. 3-4 and Figs. S4-S11). Most noteworthy is the continental Fig. 1 Process of QL wave scattering due to an anisotropic gradient. a Schematic illustration of a Love wave with transverse (T) particle motion (indicated in green), which upon encountering a lateral gradient in seismic anisotropy, partially converts into a QL wave with the same phase velocity and vertical (V) particle motion (indicated by blue) as a Rayleigh wave. The fundamental Love (G1) and Rayleigh (R1) waves travel a distance Δ between the source and receiver, whereas the QL wave only travels a distance δx between the scattering location and the receiver. The scattering is caused by a change in strength and/or orientation of the mantle anisotropy, indicated by the grey fabric pattern. b Example recording of a QL wave arrival that is found on the vertical component between G1 and R1. The distance to the scattering point (δx) is related to the time delay (δt) between the QL and G1 via the equation shown. margins (cyan symbols; Fig. 3a, b), that define the boundary between continental and oceanic crust 29 . Almost half of the scatterers (47%) fall within 200 km of this boundary. This is higher than expected considering just under two-fifths (38%) of the study area falls within 200 km (Fig. 4). The pattern is particularly intriguing around the Tasman Sea, along the former rift boundary between Australia and Te Riu-a-Māui/Zealandia (Fig. 3b). The scatterers appear to closely follow curvatures in this boundary and even coincide with micro-continents such as the East Tasman Plateau (ETP) and Gilbert Seamount Complex (GSC), which were isolated from the mainland during the rifting process 30 .
The statistical significance of such a relationship between the location of scatterers and the ocean-continent boundary can be explored further. A two-sample Kolmogorov-Smirnov (K-S) test provides a measure of the difference between two cumulative distribution functions (CDF) [31][32][33] , against which the null hypothesis, that two datasets are drawn from the same continuous distribution, can be tested. The null hypothesis is rejected when the K-S probability is smaller than a critical value based on the significance level, indicating that the two datasets are likely drawn from different distributions. As is standard, a significance level of 0.05 (5%) is applied here 33 . Figure 4a shows the cumulative distribution function (CDF) of the minimum distance between each scatterer and the ocean-continent boundary (cyan line). This is compared with the average CDF (solid black curve) for 100 distributions (thin grey lines) of 275 randomly located points across the study area, the same as the number of QL scatterers in the dataset. The comparison indicates that there are more QL scatterers closer to the ocean-continent boundary than seen in the averaged properties of the 100 random distributions (i.e. the cyan curve is above the upper dashed black curve out to a distance of~600 km). A two-sample K-S test [31][32][33] for the two CDFs returns a K-S probability of 0.04 (4%), indicating that the null hypothesis can be rejected, and thus suggesting that the QL scatterers can be distinguished from a random distribution, although the significance is not especially strong. However, sub-dividing the dataset into oceanic versus continental scatterers (Fig. S10) indicates that the spatial correlation is stronger within the oceans (K-S probability = 0.01), with 38% of the oceanic scatterers located within 100 km of the ocean-continent boundary, which is more than 2.5 standard deviations higher than the mean of the random distributions.
Away from the continental margins, scatterers can also be found within the Australian continental interior coinciding with the boundaries between major geological provinces 34 (yellow symbols; Fig. 3). All 72 scatterers that fall within the Australian mainland are located less than 156 km away from such a boundary, and 94% are located within 100 km (Fig. 4c, d). If the same number of scatterers were distributed randomly across the continental mainland, we would typically expect the maximum distance of scatterers to be twice as far away, up to 310 km from the boundary, and the number within 100 km to be 14% less. A two-sample K-S test [31][32][33] between the CDF of the locations of the observed scatterers and the mean CDF of 100 random distributions returns a K-S probability of 0.0009 (0.09%). This suggests that the null hypothesis; that the randomly distributed points and the continental mainland scatterers are from the same continuous distribution, can be rejected. It can, therefore, be implied that the spatial correlation is significantly better than random.
Looking more closely at where the continental scatterers are found, several appear to be associated with the Proterozoic orogenic blocks and deformation belts of central and northern Australia (Fig. 3b). These tectonic features are associated with prominent gravity anomalies (Fig. 3c), indicating that they represent major crustal structures with potentially deeper expressions in the mantle lithosphere also 35 . Multiple scatterers can be found on the edge of the strongest gravity anomalies, including the east-west trending Musgrave (M) and Arunta (Ar) blocks in central Australia (dashed white circle Fig. 3c), and the north-south trending Darling Fault at the southwest corner of the Yilgarn Craton (Y). Interestingly, there are several large amplitude scatterers associated with the Kimberley (K) region of northwestern Australia, which has been argued to be a small cratonic block underlain by an Archean mantle root 2,36 . Within the submerged continental area of Te Riu-a-Māui/Zealandia a correspondence with topographic variations at the surface can also be seen (Fig. 3a), particularly near New Caledonia (NC) and along the New Caledonia Trough and Norfolk Ridge. This hints at an association with the internal geological architecture of Te Riu-a-Māui/Zealandia similar to that found within Australia.
As well as old tectonic boundaries within the continental interior and along its margins, as would be expected multiple scatterers can be found at present day plate boundaries (white dashed lines and white circles; Fig. 3a), such as the New Hebrides Trench (NHT). Errors and uncertainties on the interpretation of Quasi-Love waves. It is important to first consider the uncertainties associated with the analysis of QL waves before making any physical interpretations. A brief summary is provided here in terms of the mechanism, magnitude, depth, and location (for a more detailed description please see 'Methods'). Numerous studies have demonstrated that QL waves are most efficiently generated by lateral gradients in seismic anisotropy (i.e. changes in either the strength or orientation of seismic anisotropy), and are difficult to reproduce with purely isotropic velocity changes [37][38][39][40] . Higher mode Rayleigh waves with velocities faster than the fundamental mode may potentially interfere within the selected time window, but shallow earthquakes (< 50 km) should preferentially excite the fundamental mode rather than higher modes 41 . In addition all  QL wave detections are visually inspected to ensure similarity in waveform shape with the Love wave (G1). The magnitude of the QL wave (i.e. size of the symbols in Fig. 3), is expressed as the relative amplitude between QL and G1. This magnitude is a function of both the strength of the lateral anisotropic gradient and the propagation azimuth relative to the anisotropic geometry. If the propagation direction is parallel or perpendicular to the symmetry axis then the amplitude will be zero. A particularly large QL amplitude is, therefore, indicative of the presence of a strong gradient but a small amplitude, or even absence of a QL wave entirely, is not conclusive evidence of a weak or non-existent gradient. It is, therefore, not possible to directly interpret the amplitudes of the observations. Depth sensitivity kernels calculated for a QL wave with a dominant frequency of 0.01 Hz (100 s) suggest a broad range of sensitivity throughout the upper mantle, with a peak around showing the fraction of all observed QL scatters as a function of distance from the ocean-continent boundary 29 (cyan line) and major geological provinces (yellow line) 34 . The same CDFs are generated for 100 sets of the equivalent number of points randomly located across the study area (thin grey lines), alongside their mean (solid black line) and one standard deviation (dotted black line). c, d Maps showing spatial distribution of QL scatterers (small white symbols) compared to contours of the distance from either the ocean-continent boundary (c) or the geological province boundaries (d). Fig. 3 Detected locations of QL wave scatterers across Australasia. QL wave scattering points (black circles) compared to a surface topography 59 , b seafloor age 29 , the ocean-continent boundary, and the Precambrian tectonic elements of Australia 60 , c free-air gravity 61 , and d azimuthal anisotropy from 3D2018_08Sv 8 and the lithospheric thickness from AuSREM 32,42 . Circles are scaled in size relative to the amplitude of the QL wave. Those which coincide (within a typical error of 100 km) with either a plate boundary (dashed white lines), the ocean-continent margin (solid cyan lines), or geological provinces within Australia (yellow lines) 34  100-200 km depth (refer to Fig. S1 of Levin et al. 23 , Fig. 9 of Chen and Park 19 ). While potential contributions from shallower or deeper depths cannot be entirely ruled out, the source of the QL waves most likely resides within the upper mantle, and not within the crust. The location of the QL scattering source is found by backprojection along the great-circle path (Fig. 1, S2, S5). Potential location errors are estimated to be on the order of 100 km, considering the typical wavelength of 100 s surface waves, the results of numerical experiments 19 , potential deviations from the great-circle path, and assumptions made the uniformity of the average Rayleigh and Love phase velocities between source and receiver (refer to 'Methods' for details). This is consistent with the spatial distribution of results found relative to certain tectonic and geological features (Figs. 3-4, S9). Nonetheless, the quantification of associated errors and uncertainties suggests they are moderately high, and should be kept in mind when interpreting the results and considering the geological implications.
Implications for the upper mantle beneath the Australian continent and its margins. The results of this study reveal the widespread presence of lateral anisotropic gradients capable of generating QL waves. The location of such gradients appear closely linked to past tectonic boundaries, both within the continental interior and along the rifted margins. Such features visible from the Earth's surface, by inference, must also extend to depth within the upper mantle. Considering the peak sensitivity of Love-to-Rayleigh scattering at 100 s period is expected to be between 100 and 200 km depth 19,23 , this depth range would sample the old thick cratonic lithosphere in central and western Australia (light colours; Fig. 3d), while for eastern Australia the lithosphere is younger (Phanerozoic) and much thinner,~75 km 42 . Except in regions of thick continental lithosphere, the QL scatterers should therefore correspond to lateral anisotropic gradients in the asthenosphere, where olivine latticepreferred orientation is expected to develop through mantle flow 43 .
A large number of scatterers in this dataset are associated with continental margins (Figs. 3a, b, 4a, b), similar to recently published findings 28 in which pervasive scatterers were detected along the eastern North America passive margin. For western and southern Australia the margin is co-located with the edge of the thick continental lithosphere (Fig. 3d, Fig. 5). QL scatterers along this boundary may therefore simply reflect the transition from active flow in the asthenosphere to anisotropy fossilized within the lithosphere, which is supported by the first-order patterns seen in tomographic models of azimuthal anisotropy (Fig. 3d, S11). Alternatively a step in lithospheric thickness may trigger edge-driven convection (EDC) at the margin 44 , which can generate high-strain rates and localize dislocation-creep deformation capable of overprinting latticepreferred orientation anisotropy to beneath the lithospheric step 45 . For the southern Australian margin, geodynamical models have suggested the formation of broad upwellings in the region behind the trailing edge of the continent due to EDC 46 , which may potentially explain the broad zone of QL scatterers in this region. Numerous scatterers along the eastern margin of Australia, and the western margin of Te Riu-a-Māui/Zealandia, are however less easy to explain given that the lithosphere appears consistently thin (<100 km) across this region. This may instead hint at a fundamental difference between how continental versus oceanic lithosphere couples to the underlying asthenosphere, perhaps in terms of differences in viscosity 47 . Alternatively finer-scale variations in lithospheric thickness along the margins may exist than those currently resolved in Fig. 3d. Such small-scale 3D topography on the lithosphereasthenosphere boundary may locally divert or disrupt the mantle flow-field, such as localized EDC processes proposed at the Newer Volcanics Province in SE Australia 45,48 . Recent modeling work 49 suggests that even when the edge of the continent is parallel with the background mantle flow field (as is generally thought to be the case for eastern Australia) dynamical instabilities in mantle flow can still be found at the edges of thick continental lithosphere, and may be enhanced where there are smaller-scale variations in the geometry of the boundary. It remains to be seen however whether the same effect could be generated beneath eastern Australia with only small changes in lithospheric thickness and/or viscosity.
Another intriguing finding of this study is the presence of numerous QL scatterers associated with the major geological provinces of Australia (Fig. 3). Many of these scatterers are located within thick cratonic lithosphere (Fig. 3d), and are associated with Archean-Proterozoic terranes and past orogenies (Fig. 3b). This suggests QL scatterers are sensitive to lithospheric anisotropy frozen-in since the last deformational episode, and therefore that they reflect the evolution of the continent. Furthermore this implies that geological boundaries are preserved over billions of years not only in the crust but also at deeper depths as anisotropic gradients in the lithospheric mantle. If fossilized lithospheric anisotropy is pervasive, yet heterogeneous with strong lateral gradients, including the potential for multiple layers to exist within cratonic lithosphere 50 , this would add interference to the SKS splitting accrued from the strong active flow in the underlying mantle, causing a loss of coherency in the overall signal. This may explain the long-standing conundrum as to why SKS splitting is perplexingly weak and complex across Australia 9-12 , despite fast Australian plate motion, and strong azimuthal anisotropy imaged beneath the continent (Fig. 3d, S11) 7,8,13 . Further studies to constrain lithospheric anisotropy may therefore help to decode the complex SKS splitting, as well as provide further insights on the past deformational history of continent.

Conclusion.
Observations of QL waves have revealed extensive lateral anisotropic gradients beneath the Australian continent and its margins that have otherwise, until now, been impossible to resolve with smooth seismic tomographic methods (e.g. Fig. 3d). Based on the newly acquired results, it is speculated that at the continental margins small-scale convective processes, such as EDC, or localized disruptions to the mantle flow field, may be occurring where the lithosphere transitions from continental to oceanic. Alternatively QL scattering may be explained by the difference between fossil anisotropy in old thick continental lithosphere and active flow in the oceanic asthenosphere. However, this is more difficult to justify where the continental lithosphere is younger and thinner on either side of the Tasman Sea. In the future QL scatterers may be used as a predictive tool, helping to constrain the geophysical signature of major crustal/lithospheric structures at depth, such as those currently hidden under cover 51 , and in particular in identifying those with deep mantle connections that may suggest potential for mineralization 32,52 .

Methods
Overview of data selection and processing. Earthquakes from 1994 to 2020 of magnitude > 6.5, depth < 50 km, and in the epicentral distance range of 70°-180°w ere selected to ensure large amplitude surface waves and enough separation to distinguish QL from R1. All data were low-pass filtered at 100 s (0.01 Hz) to amplify the fundamental Love and Rayleigh wave components and suppress overtones. Following previous studies, the delay time (δt) between the QL and G1 is calculated using the cross-correlation technique 19,21,28 . In the final step the delay time is converted to distance (see Fig. 1b), and back-projected along the great-circle path to determine the location of the scatterer. If the polarization of the surface waves suggests a deviation from the great-circle path of more than 10°then the result is discarded (Fig. S8). In this study, a few modifications have been made to the cross-correlation technique to improve detection and measurement quality. These are outlined in further detail below along with a more detailed description of the methodology.
Strategy for detection and measurement of Quasi-Love waves. Data for each event is requested from the IRIS Data Management Centre (ds.iris.edu/) for 5000 s following the origin time. After low-pass filtering at 0.01 Hz, seismograms are rotated into the vertical, radial, and transverse, and normalized by the maximum amplitude on each component. The maximum on the transverse signifies the fundamental Love wave (G1), and on the vertical the fundamental Rayleigh wave (R1), indicated by the vertical dashed green and dashed blue lines in Figs. S1 and S3. Initial quality control restricts events to those with a signal-to-noise ratio (SNR) of the Love wave greater than 5, where the noise amplitude is determined by the first 300 s of the seismogram. A cut-off of 5 was determined empirically from the data. In general large amplitude Love waves are a necessary pre-requirement to generating observable QL waves, and thus detections are more likely in the Love radiation maxima (and Rayleigh radiation minima).
Most previous studies only measure QL waves on the vertical component even though the particle motion is elliptical and distributed across both the vertical and radial components (Fig. S1b, c). Same as for R1, for QL, the vertical and radial components are 90°out of phase. The first derivative of the radial component (dashed red line Fig. S1d) matches the vertical component. This is a useful diagnostic feature. Any seismograms that do not display this characteristic during the visual inspection process are discarded, and the QL wave must be clearly visible on both components, not just the vertical. Furthermore, stacking of the vertical component with the derivative of the radial component amplifies the QL (and R1) while dampening the noise (magenta line Fig. S1e). This improves both the visibility and confidence of the QL detection.
The next step is to cross-correlate this QL stack with G1 on the transverse component. Fig. S1f shows the absolute amplitude of the cross-correlation (thin black line) as a function of lag time relative to the position of G1. Peaks of this cross-correlation function (i.e. maxima, where the derivative equals zero) are identified by blue crosses and connected by a new blue curve. Peaks in this new curve are then selected and highlighted by blue circles. Typically the largest peak is associated with the correlation between G1 and R1 (as seen in Fig. S1f). The QL peak, if present, is a secondary peak located before the highest R1 peak but after the G1 arrival (indicated by dashed green line). The inferred QL delay time is marked by a green asterisk and dotted green line. The time shifted Love wave is plotted for comparison with the QL stack in Fig. S1e.
All seismograms and cross-correlation functions are visually inspected for quality based on the criteria outlined above (see Fig. S1 and S3). Only those with a clearly visible and distinct QL wave are kept. It is possible however that multiple QL waveforms may appear in one seismogram if multiple scatterers are located along the great-circle path. In such cases, for simplicity, the first QL wave, closest to G1, and therefore closest to Australia is recorded.
The last step in the process is to convert the delay time (δt) into distance and back-project the location of the scatterer along the great-circle path. Considering that the QL wave travels distance, δx, of the total path, Δ, (Fig. 1a) with the slower Rayleigh wave phase velocity (V R ), the delay time can be expressed as: where V L is the Love wave phase velocity. To negate having to use either an average guess for the phase velocities or calculate individually for every seismogram, the following substitution can be made, where Δ is the epicentral distance between source and receiver, and t Rt L is the time delay between the Rayleigh (R1) and Love (G1) waves. Simply rearranging this equation for δx gives, The scatterer distance is therefore independent of the absolute phase velocities or the arrival times, and only depends on the relative time difference between G1 and R1. This is useful given the exact start of the R1 is difficult to determine, especially considering likely overlap with QL. Instead the time difference between the maximum amplitude of G1 and R1 can be utilized, which is much easier to calculate from the seismogram.
Errors and uncertainties on the interpretation of Quasi-Love waves. In this study, the term Quasi-Love waves refers to surface wave polarization anomalies that are characteristic of Love-to-Rayleigh scattering. This is distinct from distortions or rotations of the polarization of surface waves by azimuthal anisotropy that are sometimes also called quasi-Love waves or quasi-Rayleigh waves elsewhere in the literature [15][16][17] . Such polarization distortions are analogous to the terminology of quasi-P and quasi-S when body-wave polarizations are altered by Earth 3D heterogeneity 18 . In such polarization-distortion cases, the quasi-Love wave will exhibit particle motion that is still mostly transverse polarized but with a small added vertical component (see Fig. 2 of ref. 17 ). In contrast, Love-to-Rayleigh scattering generates Quasi-Love waves that are predominantly polarized in the vertical-radial plane, displaying ellipticity characteristic of Rayleigh waves with a 90°phase lag between the vertical and radial components 18 . By comparing the waveform shape of the QL wave on both the vertical component and the firstderivative of the radial component during the measurement process (e.g. Fig. S1d) ensures the presence of characteristic Rayleigh wave motion as expected for QL wave scattering. Incorrect interpretations due to slight polarization-distortions are thus considered unlikely in the present dataset.
There may however be certain scenarios where the region surrounding the seismic station is anisotropic, that such polarization-distortions of the Love wave on the radial and vertical components may still interfere with QL waves generated from a nearby lateral anisotropic gradient. This could potentially induce error in the measurement of the QL delay time, and thus the location of the QL scatterer, when δt is less than the typical period of G1 (i.e.~100 s). However, such interference should also modify the waveform shape of the QL wave, and would disrupt the expected pattern of Rayleigh wave ellipticity between the radial and vertical components. Under such circumstances, the waveforms are unlikely to pass visual inspection during the quality control process (e.g. Fig. S1), when a similar shape between G1 and the QL wave is required, as well as a similarity between the QL wave on the vertical component and first-derivative of the radial component. Such quality control steps should ensure that the effect of interference from polarization-distortions due to anisotropy near the station is minimized for the QL dataset presented here.
Previous studies have shown that Love-to-Rayleigh scattering is most efficiently generated by lateral gradients in seismic anisotropy (i.e. changes in seismic anisotropy either in strength or orientation) rather than changes in isotropic structure [37][38][39][40] . Lateral anisotropic gradients of only 3-5% are enough to generate QL waves with 5-10% amplitude of the fundamental Love wave 20,39 , compared to unrealistic perturbations of 20% or more which would be required to produce the same effect with isotropic velocity changes only 37,38 . Additionally QL waves are not typically seen in synthetic seismograms generated from velocity models without lateral variations in azimuthal anisotropy (Fig. S6) 21,27,28 . It is worth noting that lateral anisotropic gradients will also cause Rayleigh-to-Love scattering, generating Quasi-Rayleigh waves, in addition to Quasi-Love waves 18,53 . Such Quasi-Rayleigh waves will appear on the transverse component trailing G1 and thus are more challenging and less commonly studied compared to QL waves that arrive ahead of R1 on the vertical and radial components 53 . As a relatively large amplitude of G1 is a pre-requirement of this study (SNR > 5), the amplitude of R1 is often comparatively less, and thus Quasi-Rayleigh waves are not expected to be visible on the same seismograms as prominent QL waves (Fig. S1).
The relative amplitude between the QL wave and the fundamental Love wave is represented via the cross-correlation amplitude, which is used to determine the symbol size in Fig. 3. QL amplitude is dependent not only on the strength of the anisotropic gradient but also on the anisotropic geometry relative to the propagation path. Maximum QL conversion occurs when the propagation azimuth is 45°to the symmetry axis, and zero conversion when the propagation azimuth is either parallel or perpendicular 21,54 . This creates a symmetric four-lobed radiation pattern with both positive and negative polarities of the QL wave, which in theory could be used to infer not only the location of an anisotropic gradient but also the geometry of anisotropy. It is safe to say that large QL amplitudes can be interpreted as the existence of a strong gradient in anisotropy but smaller amplitudes do not necessarily equate to a weaker gradient. Likewise, a lack of detection of QL waves is not conclusive evidence that anisotropic gradients do not exist along the great-circle path. The angle between the propagation path and the geometry of the anisotropy may simply be non-optimal for Love-to-Rayleigh scattering to occur.
While lateral variations in seismic anisotropy are the most likely cause of the QL wave scatterers, an important consideration is at what depth they are generated. Depth sensitivity kernels calculated for QL wave scattering with a dominant frequency of 0.01 Hz (100 s) suggest a broad range of sensitivity throughout the upper mantle but with a peak around 100-200 km depth (refer to Fig. S1 of Levin et al. 23 , Fig. 9 of Chen and Park 19 ). While potential contributions from shallower or deeper depths cannot be entirely ruled out, the source of the QL waves most likely resides within the upper mantle, and not within the crust.
As well as the depth, there is also a degree of uncertainty associated with the horizontal location of the QL wave source. Firstly, the width of the lateral anisotropic gradient cannot be exactly determined. For example, a wide gradational zone can produce the same signal as a sharp localized gradient, but both are represented as a point scatterer at the center of the lateral gradient distribution 21 . Secondly, the smallest feature that can be resolved is likely on the order of 1/4th of the wavelength 19 . Surface waves with a period of 100 s have a wavelength of 400 km, suggesting~100 km resolution. This is consistent with numerical experiments based on simple anisotropic models with an equatorial band of 4% azimuthal anisotropy, in which the source of the QL scatterer could be correctly located to within 100 km following a similar methodology to that applied here 19 .
The real Earth however is much more complex than such idealized models and location errors may be influenced by 3D heterogeneity. In particular 3D velocity structure may cause deviations from the great-circle path (e.g. multi-pathing), which could result in a location error when the QL wave source is back-projected. In order to limit this, a comparison is made between the polarization of the surface waves and the source-receiver back-azimuth. Any events with deviations greater than 10°are thrown out from the dataset, but most are in agreement to within a few degrees (Fig. S8). Slight rotations away from the radial/vertical (Rayleigh) and transverse (Love) polarizations are to be expected when in the presence of seismic anisotropy [15][16][17][18] . Nonetheless, a maximum deviation of 10°from the great-circle path would result in a location error of 390 km for a QL wave scatterer at a typical distance of 20°from the receiver (Fig. S8). A typical deviation of 2°would however cause a location error of <100 km.
A location error may also result regarding the assumption that the average Rayleigh and Love phase velocities between the QL scatterer and the seismic station are the same as the average phase velocities along the entire great-circle path. Considering that many of the propagation paths cross an ocean basin before generating a QL wave and then continue through mainly continental regions after scattering at the ocean-continent boundary, this may cause the phase velocities along the respective path segments to substantially differ, either due to changes in isotropic or anisotropic properties or both. If the average phase velocities were 5% higher between the QL scatter and the receiver (i.e. δx), than between the earthquake source and the receiver (i.e. Δ), then this would directly correspond to a 5% error in the scatterer location as a function of distance along the great-circle path. For example, if a typical QL scatterer is estimated to be 20°away from the receiver, then a 5% velocity error would equate to a distance error of 1°(~110 km).
Overall based on the above considerations, typical location errors are expected to be on the order of 100 km but the further away the QL scatterers are the size of the error likely increases. This may be supported by the data as Fig. 4 suggests that QL scatterers are preferentially located near the ocean-continent boundary to within 400-500 km, and with geological provinces to within 100-150 km, based upon the steepness of the CDF curve relative to the random distributions.
Of further note, the methodology applied here is based upon detecting a single scatterer located closer to the receiver than to the source. This ensures enough separation from the Rayleigh wave on the radial and vertical components in order to be able to identify the QL wave. If a QL scatterer is located near the source then the QL arrival will overlap with the Rayleigh wave making it hard to distinguish. Likewise, in the cases where more than one potential QL wave is identifiable on the seismogram, potentially due to multiple lateral gradients in seismic anisotropy along the propagation path, only that which is closest to the receiver is included in the dataset as it is the easiest to distinguish. If the seismogram contains multiple closely arriving QL waves generating a complex overlapping signal, then this will be discarded during the quality control process as the QL wave is required to have a similar waveform shape as G1 (Fig. S1).
Lastly, it is possible for QL waves to resemble higher-mode Rayleigh waves 27,55 . In this study only shallow events less than 50 km depth are considered are these should primarily excite the fundamental modes rather than overtones, and all seismograms are low pass filtered at 100 s which should aid separation of the modes 56 . During the quality control process, the requirement that the QL wave should have a similar waveform shape to G1 (see Fig. S1), should also help distinguish potential interference from higher mode Rayleigh waves, as these would not be expected to possess the same waveform shape as G1. Following these conditions, synthetic seismograms demonstrate that the observed QL waves are unlikely to be the result of higher mode Rayleigh waves (Fig. S6) 19,27,28 .

Data availability
The seismic data used in this study is from the Australian National Seismograph Network (network code AU; www.fdsn.org/networks/detail/AU/), and can be freely obtained via the Data Management Center of Incorporated Research Institutions for Seismology (ds.iris.edu/). A full list of results is provided in Supplementary Data 1. The source data for the cumulative distribution functions displayed in Fig. 4