The Influence of Lithospheric Thickness Variations Beneath Australia on Seismic Anisotropy and Mantle Flow

Rapid plate motion, alongside pronounced variations in age and thickness of the Australian continental lithosphere, makes it an excellent location to assess the relationship between seismic anisotropy and lithosphere‐asthenosphere dynamics. In this study, SKS and PKS shear‐wave splitting is conducted for 176 stations covering the transition from the South Australian Craton to eastern Phanerozoic Australia. Comparisons are made with models of lithospheric thickness as well as numerical simulations of mantle flow. Splitting results show uniform ENE‐WSW aligned fast directions over the Gawler Craton and broader South Australian Craton, similar to the orientation of crustal structures generated during an episode of NW‐SE directed compression and volcanism ∼1.6 billion years ago. We propose that heat from volcanism weakened the lithosphere, aiding widespread lithospheric deformation, which has since been preserved in the form of frozen‐in anisotropy. Conversely, over eastern Phanerozoic Australia, fast directions show strong alignment with the NNE absolute plate motion. Overall, our results suggest that when the lithosphere is thin (<125 km), lithospheric contributions are minimal and contributions from asthenospheric anisotropy dominate, reflecting shear of the underlying mantle by Australia's rapid plate motion above. Further insights from geodynamical simulations of the regional mantle flow field, which incorporate Australian and adjacent upper mantle structure, predict that asthenospheric material would be drawn in from the south and east toward the fast‐moving continental keel. Such a mechanism, alongside interactions between the flow field and lithospheric structure, provides a plausible explanation for smaller‐scale anomalous splitting patterns beneath eastern Australia that do not align with plate motion.

Investigations of upper mantle dynamics via seismic anisotropy are, however, not without ambiguity.One additional factor to consider is the potential for frozen-in or fossilized anisotropy within lithospheric mantle (e.g., Fouch & Rondenay, 2006;Silver & Chan, 1991).The lithosphere, Earth's stiff outermost layer, should not be actively generating LPO.However, observations suggest it may preserve an olivine LPO fabric that developed either during its formation (e.g., at the mid-ocean ridge) or during past tectonic/orogenic events that could generate lithospheric deformation (e.g., Debayle & Ricard, 2013;Silver, 1996).The ability to constrain the pattern of fossilized anisotropy therefore holds great potential for revealing the behavior and evolution of Earth's lithosphere.
The Australian plate provides an excellent testing ground for interactions between seismic anisotropy, mantle dynamics, and the lithosphere-asthenosphere system.Australia is the fastest moving continent on Earth with an absolute plate motion (APM) of 7-8 cm per year toward the north-northeast (Kreemer et al., 2014).Such rapid plate motion may exert significant shear on the upper mantle at the base of the tectonic plate (e.g., Debayle et al., 2005).The Australian continental lithosphere also varies substantially both in terms of age and thickness (Figure 1).This likely has implications for the character and amount of lithospheric anisotropy across the continent, considering variations in the thickness of the anisotropic layer, as well as the different tectonic histories experienced.
In addition to varying contributions of lithospheric anisotropy, the existence of lithospheric steps and substantial 3D topography on the lithosphere-asthenosphere boundary (LAB) likely induces deviations of the upper mantle flow field as the uneven basal topography of the plate traverses the underlying asthenosphere (N.Rawlinson et al., 2017).Small-scale convective flow patterns induced by various LAB geometries can be predicted by geodynamic modeling (Duvernay et al., 2021;Farrington et al., 2010), but remain strongly sensitive to uncertain lithospheric structure and uppermost mantle rheology.Nonetheless a number of recent studies have started to link variations in lithospheric thickness/architecture beneath Australia with important and varied surface processes, including Cenozoic volcanism (Davies & Rawlinson, 2014;N. Rawlinson et al., 2017), the localization of critical mineral deposits (Hoggard et al., 2020), dynamic topography (Ball et al., 2021), and intra-plate seismicity (Bezada & Smale, 2019).
Upper mantle anisotropy is typically studied either by surface wave tomography or shear-wave splitting methodologies.Globally the inferences from surface waves and shear-wave splitting tend to agree (Wüstefeld et al., 2009).However, this has not been the case in Australia.Surface wave studies tend to see a strong signal in both azimuthal and radial anisotropy at the base of the continental lithosphere, with fast directions aligned with APM and the expected shear of the underlying asthenosphere (Debayle & Ricard, 2013;Debayle et al., 2005;Fichtner et al., 2010;Fishwick et al., 2008;Simons et al., 2002;Yoshizawa, 2014;Yoshizawa & Kennett, 2015), although not all models agree (e.g., de Laat et al., 2023).Conversely, shear-wave splitting studies using core-refracted phases, such as SKS, have typically not detected plate-motion aligned fast directions (e.g., Bello et al., 2019;Chen et al., 2021;Heintz & Kennett, 2006).
Most previous shear-wave splitting studies in Australia had ≤35 stations available, often sparsely distributed on a continental scale, making it difficult to pin-point the cause and location of spatial variations in anisotropy.The first exception is the study of Heintz and Kennett (2005), which had a large number of stations (>100) but was hindered by unusually short recording times (average < 6 months), and, therefore, a restricted number of events.The second exception is a recent study by Ba et al. (2023), who reported shear-wave splitting results from 116 stations across the continent with spatially complex patterns.Unlike previous studies, Ba et al. (2023) were able 10.1029/2023GC011066 3 of 16 to identify certain locations, such as the peripheral areas of the continent, where the splitting pattern matched well with the direction of APM.Overall, they attributed their results to asthenospheric mantle flow with possible (yet unquantified) lithospheric contributions.They did not however demonstrate a quantitative relationship with lithospheric thickness, nor did they investigate structural/geophysical trends at the surface to see if locations of thick lithosphere could be better explained by frozen-in anisotropy.
Recent results from the BILBY north-south transect (Eakin et al., 2021), however, identified a clear correspondence between the splitting fast axis and prominent structural/gravity trends across the South and North Australia Cratons, as well as the suture zone in-between.This provided strong evidence that the shear-wave splitting results in this particular cratonic region are predominantly reflecting lithospheric frozen-in anisotropy rather than asthenospheric contributions.All BILBY stations, however, were located on thick cratonic lithosphere, and so the question of variable lithospheric contributions could not be assessed.
In this study, we aim to utilize a compilation of over 170 stations that traverse the region extending from thick cratonic lithosphere beneath the South Australia Craton to thinner Phanerozoic lithosphere toward the east.The  study is supported by our two recent deployments in South Australia that provide unprecedented coverage of the eastern Gawler margin and South Australian Craton (Eakin, 2019;O'Donnell, Thiel, Robertson, Gorbatov, & Eakin, 2020).Using new data and shear-wave splitting measurements from across these seismic networks, we can determine, in detail, how variations in lithospheric thickness and architecture exert a first order control on the pattern of seismic anisotropy and upper mantle dynamics beneath Australia.

Tectonic Overview
A large difference in age exists between the western two-thirds of Australia, which is of Precambrian/Proterozoic origin, and the eastern Phanerozoic margin.Precambrian Australia consists of three main crustal components: the North, South, and West Australia Cratons (Figure 1a).These are thought to have existed since ∼1.8 Ga and were assembled by 1.3-1.1 Ga, during the Rodinia supercontinent cycle, to form proto-Australia (Betts & Giles, 2006;Myers et al., 1996).In contrast, the eastern third of Australia can be described as a series of accretionary or orogenic belts, added to the eastern margin of proto-Australia during the Cambrian to Triassic periods (∼550-220 Ma) via subduction (e.g., Glen, 2005).The geological boundary demarcating cratonic Precambrian Australia from eastern Phanerozoic Australia is sometimes referred to as the Tasman Line (Direen & Crawford, 2003) (Figure 1).The precise location of the Tasman Line, inferred from various geophysical and geological observations is, however, often poorly constrained, especially in regions of thick sedimentary cover such as beneath the Lake Eyre region (Agrawal et al., 2022).Nonetheless, it is clear from Figure 1b that western Precambrian Australia is generally underlain by thick lithosphere (>150 km), whereas eastern Australia is characterized by relatively thin continental lithosphere, typically ∼75 km thick (Fishwick et al., 2008;Kennett & Salmon, 2012).
Our study area covers most of the South Australian Craton, the core of which is composed of the Archean-Proterozoic Gawler Craton and the Proterozoic Curnamona Province (Figures 1 and 2).The Gawler Craton is the oldest and largest geological province in South Australia and preserves a complex tectonic history spanning from 3,200 to 1,450 Ma (Myers et al., 1996).The Archean and Paleoproterozoic core of the Gawler Craton forms a folded belt that underwent deformation along discrete shear zones during the Mesoproterozoic (Hand et al., 2007).Throughout the Proterozoic, the region saw multiple major magmatic events, recorded in the geological record by the Donington Suite (ca.1,850), St. Peter Suite (1,620-1,610), and Gawler Range Volcanics-Hiltaba Suite (1,595-1,575 Ma) (Hand et al., 2007).
Compared to the Gawler Craton, the Curnamona Province is largely hidden under cover (Agrawal et al., 2022), but its oval-shaped outline is easily distinguishable from aeromagnetic data (e.g., Williams & Betts, 2007).Sparse outcrops suggest that much of the province is composed of metasedimentary rocks of the late Proterozoic (1,715-1,640 Ma) Willyama Supergroup (Williams & Betts, 2007).Two major orogenic events are associated with the province: the Mesoproterozoic Olarian orogeny (1,600-1,590 Ma, Page et al., 2005) and the early Paleozoic Delamarian orogeny (520-490 Ma, Dutch et al., 2005).Coeval with the Gawler Range Volcanics-Hiltaba Suite recorded in the adjacent Gawler Craton, the Curnamona Province also experienced felsic magmatism with the Willyama Supergroup intruded by granites during the same period (e.g., Page et al., 2005).

Station and Event Availability
In total, 176 stations in Australia were analyzed for SKS and PKS splitting in this study.The distribution of these stations covers a wide area including the South Australia Craton (Gawler and Curnamona Provinces), as well as Phanerozoic regions to the east and south-east (Figure 2a).Two new seismic deployments in South Australia provide increased station coverage over the eastern Gawler Craton (yellow and green symbols in Figure 2a).The first of these deployments is the Lake Eyre Basin Seismic Array (Network: 5G), including 22 broadband stations deployed in the region surrounding Kati Thanda-Lake Eyre (Eakin, 2019).Instruments were installed in several phases between September 2018 and October 2019, and remained in-place until July 2022.Several data gaps exist, including the period from June to October 2020 resulting from COVID-19 related state border closures and travel restrictions that prevented servicing of the network.A second array, AusArray-SA (Network: 6K), was deployed south of the Lake Eyre Basin array by the Geological Survey of South Australia (O'Donnell, Thiel, Robertson, Gorbatov, & Eakin, 2020).This consisted of 38 broadband stations that operated over a similar timeframe from October 2020 to June 2022.
The station coverage is supplemented by data from two permanent national networks.These include 23 stations from the Australia National Seismic Network (red symbols in Figure 2a, Network: AU) and 17 stations from the Australian Seismometers in Schools program (blue symbols Figure 2a, Network: S1).The remaining 76 stations (black symbols) represent a compilation of eight past temporary broadband deployments.Most of these sites operated for 12-24 months.A full list of all temporary (and permanent) networks is provided in the data availability section.The majority of these sites have not been previously analyzed for SKS splitting.Those that have are re-analyzed to ensure consistency in the methodology across the region.The only exceptions are results from the BILBY array (stations BL:15,16,17,19,20,24) and permanent station AU:MULG from our earlier study of Eakin et al. (2021).This study followed the same methodology as is applied here and therefore, these results are directly included in the data set without re-processing.
For event selection, earthquakes of magnitude 6.0 and above were utilized in the epicentral distance range 90°-130° for SKS and 130°-150° for PKS phases (Figure 2b).This equates to around 30-40 events per year for a seismic station in southern Australia.While events originate across a range of backazimuths, they are most plentiful from the South America subduction zone from the south-southeast direction (Figure S4 in Supporting Information S1).

Shear-Wave Splitting Processing and Methodology
Shear-wave splitting analysis was undertaken using the SplitLab software package (Wüstefeld et al., 2008).Multiple methods for estimating the best-fitting shear-wave splitting parameters (the fast direction: Φ and delay time: δt) are available within the SplitLab environment.For quality control purposes, we compare estimates from two independent methods.The first is the rotation correlation (RC) method (Bowman & Ando, 1987), which determines the values of Φ and δt, which generate the maximum cross-correlation between the trial fast and slow components.This method can produce systematic error as a function of initial polarization or as a function of backazimuth for *KS phases (Eakin et al., 2019;Wüstefeld & Bokelmann, 2007).When the initial polarization approaches the fast or slow orientation of the anisotropic medium, the RC method will predict a best fitting value of Φ that deviates 45° from the true value, with a delay time that is close to zero (e.g., Figure S9 in Supporting Information S1).While such systematic measurement error is not ideal, it is well understood and predictable.With this in mind, the true splitting parameters can still be easily retrieved from this method (Eakin et al., 2019).
The second method is the transverse energy minimization (SC) method (Silver & Chan, 1991).Using this method, we seek those values of Φ and δt that best minimize the energy on the transverse component, following a correction for shear-wave splitting.The SC method is thought to produce more stable SKS splitting measurements over a wider range of backazimuths.However, it can also be susceptible to the same systematic error as seen for the RC method when the signal-to-noise ratio is moderately high (Eakin et al., 2019;Wüstefeld & Bokelmann, 2007).Unlike the RC method, the SC method is dependent on the estimated initial polarization.The quality of results from the SC method is therefore particularly sensitive to any misalignment of the station orientation and/or miscalculation of the back azimuth (Eakin et al., 2018).For these reasons, and unless otherwise specified, the results presented in the following sections are from the RC method.
Quality control procedures followed our previous work (Eakin et al., 2015(Eakin et al., , 2019(Eakin et al., , 2021)), including visual inspection of all waveforms and strict quantitative and qualitative criteria.Further details of the criteria used to determine whether a split (or null) measurement was of acceptable quality are provided in Text S1 in Supporting Information S1.As an initial step, waveforms were filtered between 0.04 and 0.125 Hz using a Butterworth bandpass filter.Any event with a signal-to-noise ratio (SNR) of less than 5 for the SKS/PKS phase was discarded.Using this initially curated data set, a station misalignment value was estimated.This estimate follows similar investigations by Eakin et al. (2018) that measured the difference between the initial polarization of SKS/PKS phase (as determined from the orientation of the uncorrected particle motion) and the source-receiver backazimuth.If the station is found to be misaligned with north, then the appropriate orientation correction is applied to all waveforms before further analysis.Station misalignment values applied in this study are available in Table S1 in Supporting Information S1.
To help identify regional variations in shear-wave splitting patterns, a stacked result was calculated for each station that had multiple measurements.This was estimated by stacking of the RC error-matrices from all individual split measurements.The station average Φ and δt values are then found from the global minimum of the stacked error surface (Wolfe & Silver, 1998).A further check is performed between the stacked δt result and the mean δt from the individual measurements with a ratio of <1.75 required.This ensures that the stacked splitting results presented in the final data set are consistent with the average properties of the individual measurements.

Mantle Flow Modeling
We build on the models developed by Davies et al. (2019) to generate a synthetic mantle flow field for comparison with seismic anisotropy observations.Although we focus on the Australian region, our model is global.We solve the equations governing incompressible mantle convection inside a spherical shell using Fluidity (Davies et al., 2011;Kramer et al., 2012).In our simulation, the inner radius corresponds to the Core Mantle Boundary (CMB) and the outer radius to Earth's surface.Free-slip mechanical boundary conditions are specified at the CMB, with present-day plate kinematics from Müller et al. (2016) prescribed at the surface.Our mesh is generated by refining an icosahedron, resulting in a lateral resolution of 50 km at the surface.This mesh is extruded in the radial direction, with radial spacing increasing linearly from 10 km at the surface to 100 km at the CMB.
We determine the present-day density and temperature fields by adopting a robust thermodynamic approach for converting between seismic and physical structure.This approach, described in Ghelichkhan et al. (2021), uses the upper mantle tomography model of de Laat et al. ( 2023) above 300 km depth, transitioning smoothly to the whole mantle shear-wave tomography model of S40RTS (Ritsema et al., 2011) for depths below.The former provides higher resolution on lithospheric structure owing to the sensitivity of surface waves.We use Perple_X (Connolly, 2009) alongside the thermodynamic database of Stixrude and Lithgow-Bertelloni (2011) to determine equilibrium phase assemblages throughout the mantle as a function of temperature, and pressure, and their associated anharmonic V S and density values for pyrolytic mantle.To account for anelastic effects, anharmonic V S is corrected using an updated version of the Q 5 model by Cammarano et al. (2003).We employ a temperature-and depth-dependent viscosity field, with parameters that are identical to those of Davies et al. (2023), constructed to be compatible with observations of Earth's geoid, heat flux, post-glacial rebound, and CMB ellipticity.

SKS and PKS Splitting Results
Overall, our splitting analysis yielded 1,207 measurements across 157 stations, all of which were individually and visually inspected.A full list is provided in Table S1 in Supporting Information S1.While a single measurement by itself can be unreliable, 67 stations had multiple non-null measurements that produced a stacked splitting result (Figure 3).By utilizing this smaller curated data set, we focus our attention on those stations with the most robust results.In general, the pattern of individual measurements across all the stations in the study area agrees well with the pattern retrieved from these stacked station results (Figure S2 in Supporting Information S1).Even with the reduced number of 67 stations, this is still 50% more compared to the most extensive prior study by Ba et al. (2023), who only reported results from 44 stations across the same study area.
Considering the map of stacked split results in Figure 3a, a broad similarity of splitting parameters can be seen between neighboring stations, indicating a spatially coherent source of seismic anisotropy below.Throughout the Gawler and South Australian Craton, there exists a strikingly consistent ENE-WSW splitting pattern (pink/red bars in Figure 3a).The average fast direction (Φ) for the 22 stations within the Gawler Craton is 70°, and 0.80 s for the delay time (δt).This is consistent (within typical errors) with our previous study of permanent station MULG (Φ: 79°, δt: 1.1 s) located within the Gawler Craton (Eakin et al., 2021).Other recent studies have reported a similar NE-SW to ENE-WSW trend over the South Australian Craton (Ba et al., 2023;Birkey & Ford, 2022).The spatial extent of this shear-wave splitting pattern, throughout the Gawler Craton and extending eastwards into the Curnamona Province, can now be seen more clearly.Beyond the craton margins, the pattern changes slightly.To the craton's northwest there are only a handful of stations, but these show a change to WNW-ESE fast directions (yellow/orange bars Figure 3a).To the north of the Gawler Craton, over the Lake Eyre region, the fast direction rotates slightly to become more north-easterly (purple colored bars Figure 3a).A few anomalous stations with a NNE-SSW fast direction are seen to the east of Lake Eyre between the Gawler Craton and the inferred Tasman Line (blue colored bars Figure 3a).These appear to coincide with where there is an indent into the thick cratonic lithosphere in the AuSREM LAB model (Figure 3b).
The most notable change is seen crossing from western cratonic Australia into eastern Phanerozoic Australia.The majority of stations in eastern Australia display larger delay times (mean δt: 1.26 s) and N-S to NNE-SSW fast directions (mean Φ: 20°), represented by blue colored bars in Figure 3a.This orientation is very similar to the APM of the Australian plate, ∼6 cm per year at 21° (clockwise from North) in this location, as indicated by the black arrow in Figure 3. Hints of such a correspondence between the fast direction and the APM for eastern Phanerozoic Australia have been previously noted based on more limited results (Ba et al., 2023;Bello et al., 2019;Birkey & Ford, 2022).The results of our analysis confirm a similarity between the splitting fast direction and the APM for southeastern Australia that is spatially distinct from the pattern over cratonic Australia.
Interestingly, there is a small cluster of five stations in the very south-eastern corner of Australia (∼latitude: 37°S, longitude: 148°E) that show a contrasting but consistent ESE-WNW orientation (ave Φ: −76°, δt: 1.0 s), as indicated by the cluster of orange bars in Figure 3a.Some of these same stations have been previously analyzed in other studies: results from Ba et al. ( 2023) indicated a similar pattern for this sub-group of stations; however, the pattern from Bello et al. (2019) was less spatially consistent.

Comparison With Lithospheric Thickness
Across the study region, lithospheric thickness reduces drastically from west to east (Figure 1b).Beneath the Gawler and South Australian Craton, where the splitting fast directions are consistently ENE-WSW, the lithosphere is ∼200 km thick (Figure 3b).Over eastern Phanerozoic Australia, where stations show a strong similarity between the fast direction and APM (indicated by red bars in Figure 3b), the lithosphere is much thinner with a typical LAB depth ∼70-100 km.
A direct comparison of the splitting parameters against LAB depth is plotted in Figure 4 to better quantify this observation.Across the study region, the magnitude and range of delay times vary relative to the lithospheric thickness (Figure 4a).This varies from median δt values > 1.0 s and maximum values larger than 1.5 s when the LAB is shallower than 125 km, to an average delay time of ∼0.8 s and maximum values generally under 1.5 s when the LAB depth is >175 km.Correspondingly, the misfit between APM and the fast splitting direction is <20° for most stations located where the LAB depth is <125 km (Figure 4b).The cluster of blue points in Figure 4b that defy this relationship corresponds to the anomalous cluster previously identified in the southeastern corner of Australia (i.e., blue bars Figure 3b).For LAB depths greater than 150 km, the misfit between the splitting fast direction and APM tends to grow as the lithosphere gets thicker.
There are a small handful of stations with stacked splitting results that fall along the edge of the thick cratonic lithosphere as it transitions to thinner lithosphere.These stations tend to show fast directions that are approximately parallel to the general geometry of the lithospheric step, such as those that fall within the white-dotted band in Figure 3b.

Comparison With the Predicted Mantle Flow Field
Uppermost mantle structure and lithospheric thickness variations beneath Australia will influence the underlying asthenospheric flow regime.Geodynamical models that incorporate such upper mantle structure can provide insights on what dynamic processes may be occurring, helping to inform shear-wave splitting observations.In Figure 5, our splitting observations are compared to predictions across the region from a mantle flow simulation (as outlined in Section 3.3), which is the first time such a comparison has been made directly for the Australian upper mantle.Where the inferred upper mantle viscosity is high at 200 km depth (indicative of the lithosphere), the predicted flow field (small black arrows in Figure 5) generally aligns with APM, as prescribed at the surface, demonstrating rigid plate motion within the high viscosity lid.The observed ENE-WSW orientated splitting pattern over this region of high viscosity is therefore inconsistent with the present-day mantle flow direction predicted in our simulation.
At the same depth, beneath eastern Australia, the mantle viscosity is predicted to be lower (Figure 5), reflecting the presence of thinner lithosphere and shallower asthenospheric material.In this region the predicted mantle flow field has a greater E-W component, with mantle material from the east being drawn in toward high-viscosity thick lithosphere in the west and rotating to align with plate motion upon interaction with lithospheric structure.In the simulation, as the Australian continent moves rapidly northwards, the motion of the cratonic lithosphere through the upper mantle creates a region of lower pressure, with lower viscosity asthenospheric material from surrounding areas drawn toward this region.This westwardly flow is strongest in the model in the south-east corner of the study region, matching well (particularly at 200 km depth, Figure 5b) with the sub-group of five stations that showed anomalous ESE-WNW orientated splitting (orange bars in Figure 3a).However, other stations located adjacent to this anomalous cluster (minimum station separation: 140 km) still show fast directions parallel to the plate motion.While such lateral changes in upper mantle anisotropy (over length scales of 140 km) are entirely possible to resolve with shear-wave splitting measurements, they are not possible to replicate in our simulation where the mantle flow field, and the viscosity structure, changes more gradually.This is due to the inherent resolution of the tomographic models used in deriving the thermochemical structure.Overall, however, the shear-wave splitting pattern over eastern Australia shows a greater similarity to the APM, rather than the asthenospheric flow field predicted in the mantle flow simulation.

Discussion
Our results demonstrate clear and spatially coherent shear-wave splitting patterns (Figures 3 and 4).This splitting pattern appears to morph from an ENE-WSW alignment over the thicker lithosphere of the South Australian Craton, to a NNE-SSW alignment, that follows the APM for stations located above the thinner Phanerozoic lithosphere in the east.This clear pattern is in contrast to the weak splitting and complex anisotropy often inferred by the earliest studies of SKS splitting in Australia (e.g., Clitheroe and van der Hilst, 1998).Further context regarding the emergence of coherent seismic anisotropy beneath Australia can be found in the Text S2 in Supporting Information S1.As is common practice for *KS splitting studies at the continental scale (e.g., Ba et al., 2023;Eakin et al., 2010), we have presented stacked splitting results at each station.Such stacked results allow for interpretation in terms of a single (and relatively simple) layer of anisotropy.Analysis of the back-azimuthal variations lends support to such a single-layer interpretation (Figures S9 and S10, and we find that complex or multi-layered anisotropy is not required to explain the results (refer to Text S3 in Supporting Information S1 for extended details).

Anisotropic Contributions From the Lithosphere Versus Asthenosphere
Placing the results into the context of simple (i.e., single-layer) anisotropy allows for a straightforward interpretation that is consistent with geological constraints.The relationship between lithospheric thickness and alignment of the fast direction with the APM (Figures 3 and 4) can be easily explained.Where the lithosphere is relatively young and thin (<100 km thick) over eastern Phanerozoic Australia (Figure 3), the contribution from lithospheric anisotropic fabrics is likely small.Instead, the anisotropic signal from the asthenosphere can dominate, resulting in splitting fast directions that match with the APM (Figure 3).This suggests shear of the underlying mantle asthenosphere by the fast plate motion above, as has long been proposed by surface wave studies that imaged strong APM aligned azimuthal anisotropy at the base of the plate (e.g., Debayle et al., 2005).
Conversely, where the lithosphere is relatively old and thick (>175 km thick) over cratonic Precambrian Australia (Figure 3), the consistent ENE-WSW orientations suggest that anisotropic contributions from fossilized lithospheric fabrics dominate the shear-wave splitting signal.This is not to say that strong shear of the asthenosphere below is no longer occurring (as has been imaged by surface waves) but that the shear-wave splitting appears most sensitive to shallower lithospheric anisotropy in the mantle's uppermost 200 km.Such an observation is supported by previous modeling of vertically varying anisotropy and synthetic seismograms by Saltzer et al. (2000).Their work suggested that *KS splitting measurements may be more biased toward the upper portion of the anisotropic medium with the greatest sensitivity at around one-third of the depth of the total extent of the anisotropy.They suggest that this explains why many stable continents display fast splitting directions that mirror surface geology (e.g., Silver, 1996), which is also the case here (Figure 6).
The reduced delay times (<1 s) for stations situated on-top of thick lithosphere (Figure 4a) however would suggest that when the primary source of anisotropy is from the lithosphere, this is not as strong as when the lithosphere is thin (<125 km) and the primary source is more likely from asthenospheric flow.Alternatively, when the lithosphere is thick and lithospheric contributions dominate, a more minor but opposing contribution from the deeper asthenosphere could cause a slight reduction of the delay times.
Considering that the splitting pattern changes drastically from west to east across the study area, the most plausible explanation appears to be a change in the primary source of anisotropy from lithospheric to asthenospheric.In Australia, studies of azimuthal anisotropy from surface waves also tend to show clear east-west trending fast axes at lithospheric depths (∼50-150 km) within the South Australian Craton (e.g., Fishwick et al., 2008;Simons et al., 2002).Moving eastwards as the lithosphere thins, these fast axes tend to change orientation within the same depth slice (e.g., Figure S7 in Supporting Information S1), but it is difficult for surface waves to delineate lateral contrasts in seismic anisotropy.In contrast, Quasi-Love wave scatterers, which are sensitive to lateral gradients in seismic anisotropy at upper mantle depths (∼100-200 km) (Eakin, 2021), can often be traced to the edge of thick lithosphere and to locations where a change in our shear-wave splitting pattern occurs (Figure S6 in Supporting Information S1).

Implications for the South Australian Craton
Over the South Australian Craton, which encompasses the Archean-Proterozoic Gawler Craton and Curnamona Province, the *KS splitting displays a consistent ENE-WSW orientated pattern (Figure 3a).North of our study area, the orientation of the fast direction has previously been shown to flip to match the NW-SE aligned terrane boundaries within the North Australia Craton (Eakin et al., 2021).This likely indicates that the lithospheric anisotropic fabric differs between the two major cratonic domains, reflecting different tectonic histories.
Within the South Australia Craton, upon first impressions, the ENE-WSW orientation of the anisotropy does not appear to match the outline of the major geological provinces or crustal boundaries (Figure S8 in Supporting Information S1), as noted by Birkey and Ford (2022).However, upon closer inspection, the ENE-WSW orientated fast directions agree very well with the structural trends found internally within the Gawler Craton, which are preserved as a series of E-W to NE-SW trending faults (Figure 6).Several of these major crustal shear zones and deformation structures (highlighted by thick black dashed lines in Figure 6) either formed or were reactivated during the mid-Proterozoic around 1.6 to 1.58 billion years ago (Hand et al., 2007).This time period is associated with the last major deformational event that impacted the entire Gawler Craton, due to large-scale NW-SE directed compression (indicated by the thick black arrows in Figure 6) (Hand et al., 2007).Such a tectonic history is consistent with the anisotropic geometry retrieved by this study.
With such strong crustal deformation, anisotropy within the crust could be expected.However, an anisotropic inversion of Rayleigh wave phase velocities between 2.5 and 10 s suggests very weak crustal anisotropy (<0.5%) throughout the Gawler Craton (Pilia et al., 2016).With an average crustal thickness of ∼46 km across the craton (Agrawal et al., 2023), even if the entire crust contributed, this would equate to delay times of less than 0.1 s for 0.5% shear-wave anisotropy.It therefore seems unlikely that crustal anisotropy accounts for the shear-wave splitting pattern seen over the Gawler Craton (Figure 6).Instead, this strongly implies that the splitting pattern is a manifestation of fossilized LPO fabrics frozen into the lithosphere during past tectonic deformation and that this fabric has been preserved for well over a billion years.The consistency of the splitting pattern over such a large area suggests that deformation of the lithosphere during the Proterozoic was uniform over a wide area, likely representing the last deformational event before the region fully cratonized.
Intriguingly, at the time of deformation, the region also experienced significant volcanism and heating, as recorded by the emplacement of the Gawler Range Volcanics (∼1,590) and the Hiltaba Suite (1,595-1,575 Ma).Both magmatic events impacted large swaths of the Gawler Craton, as indicated in Figure 6.The introduction of heat and melts from this volcanism may have substantially weakened the lithosphere and further aided widespread deformation, producing the significant lithospheric LPO fabrics still seen today.After the emplacement of the Hiltaba Suite, volcanism ceased in the region, allowing the lithosphere to cool, strengthen, and preserve the anisotropy over the following 1.6 billion years.

Deviations of the Asthenospheric Flow Associated With LAB Topography
While there appears a clear relationship between the splitting parameters and lithospheric thickness (Figure 4), there are some locations where the splitting fast directions neither follow the APM-aligned mantle flow direction nor the inferred orientation of the lithospheric fabric.One such location is along the eastern edge of the thick cratonic lithosphere (highlighted by the dotted white band in Figure 3b).While there are only a small number of stations in this zone, intriguingly, the orientation of the fast direction appears to follow the general trend of the lithospheric step, as indicated by the AuSREM LAB model.This perhaps hints at a deviation of the mantle flow field around the edge of the cratonic root.
Such a mechanism (of flow diverted around the continental keel) can be seen in our mantle flow simulations, particularly at 120 km depth along the north-eastern boundary of the higher viscosity region (Figure 5a).However, this transition is sensitive to the exact lithospheric structure, the resolution of which remains under refinement by seismic imaging.In the tomographic model Aus22 (de Laat et al., 2023), used as an input for the flow simulations, the lithospheric step is not as sharp as that suggested by the AuSREM LAB model (Figure 3b).Instead in Aus22 there is a more gradual transition, with a region of intermediate lithosphere identified separating the thick cratonic lithosphere in the west with the thin Phanerozoic lithosphere in the east.The location of the eastern boundary of the cratonic lithosphere in the Aus22 model is very similar to the AuSREM LAB model.However, at the depths illustrated, mantle flow in our simulations appears to deflect around the boundary of the intermediate lithosphere (offset to the east), rather than the cratonic lithosphere as our shear-wave splitting measurements suggest (Figure 3b).Regardless, it is clear from the simulations that the existence of lithospheric steps beneath Australia is likely to invoke deviations in the upper mantle flow field, with flow vectors parallel to the boundary (Figure 5a), but that where such deviations occur will strongly depend on the underlying lithospheric structure.
Deviations of flow around the continental keel have previously been invoked to explain departures of the splitting fast direction from the APM of Australia (Ba et al., 2023;Bello et al., 2019;Clitheroe & van der Hilst, 1998;Heintz & Kennett, 2005).However, it has not previously been imaged parallel to, and directly in the vicinity of a proposed lithospheric step, as seen in Figure 3b.Alternatively, some of the stations within this band, particularly in the north, match well with the predicted upper mantle flow from our simulation (Figure 5).This suggests that the pattern could also be reproduced by mantle flow drawn in toward the fast-moving continent, rather than flow diverted around the lithospheric step.
Another interesting location where the splitting fast direction departs from the general trend of APM alignment is the cluster of results in the southeast corner of the study area, highlighted in blue in Figures 3b and 4. In this region, the lithosphere is quite thin (<100 km), so a substantial contribution from lithospheric fabrics would not be expected.Intriguingly, however, the splitting pattern in this location matches well with the westward directed mantle flow predicted by our geodynamic simulation along eastern Australia, especially at ∼200 km depth (Figure 5b).This may suggest that it is the result of a similar process, whereby asthenospheric material from beneath the Tasman Sea is being drawn in toward the higher viscosity lithospheric keel of the Australian continent.Why the process would be localized to just this area remains unclear, but perhaps suggests more complex lithospheric architecture in this region than is presently resolved by tomographic models, such as that by de Laat et al. ( 2023) utilized here.
To the south of this E-W orientated splitting cluster lies a well-imaged low velocity anomaly in the shallow upper mantle of the Bass Strait (e.g., de Laat et al., 2023).This low velocity feature has been variously attributed to localized edge-driven convection (Davies & Rawlinson, 2014;N. Rawlinson et al., 2017) and/or a mantle plume (Davies et al., 2015), both of which would induce small-scale deviations of the mantle flow field.Additionally, strong scattering of Love-to-Rayleigh waves beneath the Bass Strait has been observed (Eakin, 2021), indicative of lateral gradients in upper mantle anisotropy.Together, these observations and models suggest that small-scale convective processes are likely occurring beneath the Bass Strait that may influence the anomalous patterns of seismic anisotropy beneath southeastern Australia.

Conclusion
From analyzing *KS splitting at over 170 stations and focusing on those that could provide high-quality stacked results, coherent but regionally variable anisotropic patterns have emerged beneath south and eastern Australia.In regions of comparatively thin lithosphere (70-125 km), we find fast directions aligned with Australia's rapid plate motion, as predicted from surface wave studies.This demonstrates that lithospheric contributions are minimized beneath eastern Australia, with splitting patterns principally governed by asthenospheric flow.Conversely, where the lithosphere is thicker beneath the South Australian Craton, ancient deformational fabrics appear to dominate 10.1029/2023GC011066 13 of 16 that have likely been preserved within the lithospheric mantle.The trend of this anisotropic fabric is comparable with deformational structures from the Early-Mesoproterozoic found within the Gawler Craton.Such findings imply widespread and uniform deformation of the region, coeval with the addition of heat and the emplacement of the Hiltaba Suite and Gawler Range Volcanics, before the South Australian lithosphere cratonized at ∼1.6-1.5 Ga.
Lithospheric thickness variations therefore appear to exert a first-order control on the anisotropic signal retrieved from shear-wave splitting beneath Australia.While continental scale imaging of the LAB beneath Australia has been achieved by seismic tomography (e.g., Kennett & Salmon, 2012), it is more difficult to resolve sharp lateral contrasts in lithospheric age and thickness that may demarcate important tectonic boundaries such as the enigmatic Tasman Line.With the ongoing expansion of seismic data collection across Australia, further detailed investigations of the anisotropic structure may therefore hold significant potential for unveiling the expansive tectonic history and lithospheric architecture of this ancient continent.

Data Availability Statement
A table of all shear wave splitting measurements made during this study can be found in Table S1 in Supporting Information S1, and is available from the ANU Data Commons repository (https://datacommons.anu.edu.au/;Eakin, 2023).Seismic data from the Lake Eyre Basin (Network: 5G, Eakin, 2018) and AusArray-SA (Network: 6K, O'Donnell, Thiel, Robertson, Gorbatov, & Goleby, 2020) will become publicly available in June 2024 from the Australian Passive Seismic Server (AusPass, http://auspass.edu.au/)hosted at the Research School of Earth Sciences, Australian National University.Waveforms and meta-data from the ensuing list (network name, network code, and corresponding citation) of contributing permanent and temporary seismic networks were accessed through AusPass and/or the Incorporated Research Institutions for Seismology Data Management Center (IRIS DMC, https://ds.iris.edu/ds/nodes/dmc):ANSN, AU (Geoscience Australia, 2021); AUSIS, S1 (Australian National University, 2011); SQSPA, 1E (Tkalčić et al., 2013); BASS, 1P (Reading & Rawlinson, 2011); AQT, 1Q (AusPass, 2016); Marla Line, 3G (Kennett, 2018); ASR, 5J (Salmon & Kennett, 2017); BILBY, 6F (S.Rawlinson & Kennett, 2008); TIGGER BB, 7H (N.Rawlinson & Kennett, 2001); TASMAL, 7I (Kennett, 2003); and SOC, 7K (Fontaine & Kennett, 2007).deployed.SA and CME were supported by Australian Research Council Grant DE190100062.The Lake Eyre Basin seismic array, and many of the previous temporary deployments, were made possible by funding from AuScope (https:// auscope.org.au),instrumentation from the Australian National Seismic Imaging Resource (ANSIR), and contributions from staff at the Research School of Earth Sciences of the Australian National University, most notably Robert Pickle and Michelle Salmon.The AusArray-SA deployment was supported by the Geological Survey of South Australia (GSSA), with instrumentation from ANSIR and Geoscience Australia.The GSSA thanks Isaac Axford, Goran Boren, Ann Goleby, Liz Jagodzinski, Christine Selway, Kate Selway, John Stephenson, Judy and Ed Zajer, Michelle Salmon, Robert Pickle, and Colin Telfer for their invaluable contributions to AusArray-SA.We are incredibly grateful to landholders, traditional owners, and the Department of Defence for granting land access for the seismic arrays.JPOD publishes with the permission of the Director of the Geological Survey of South Australia.Splitting measurements were processed using SplitLab version 1.2 updated by Robert Porritt (https://robporritt. wordpress.com/software/).Figures were made with the aid of Generic Mapping tools (Wessel et al., 2013), scientific color maps from Crameri et al. (2020), and geological maps available via the South Australia Resources Information Gateway (SARIG, https://map.sarig.sa.gov.au/).We thank Janneke de Laat for sharing the Aus22 shear-wave velocity model prior to its subsequent publication (de Laat et al., 2023).Our geodynamical simulations were supported by the Australian Government's National Collaborative Research Infrastructure Strategy (NCRIS), with access to computational resources provided on Gadi through the National Computational Merit Allocation Scheme and the ANU Merit Allocation Scheme.Finally, we thank two anonymous reviewers and the associate editor for their constructive comments, which helped improve the manuscript.Open access publishing facilitated by Australian National University, as part of the Wiley -Australian National University agreement via the Council of Australian University Librarians.

Figure 1 .
Figure 1.Overview of the Australian continent indicating (a) the location of cratons and (b) the variation in lithospheric thickness from AuSREM (Kennett et al., 2012).The oldest Precambrian provinces of Australia from Raymond et al. (2018) are highlighted in beige and pink colors in panel (a).The inferred boundary between Precambrian and Phanerozoic Australia, often referred to as the Tasman Line (Direen & Crawford, 2003), is indicated by the dashed white line.The approximate extent of the West Australian Craton, North Australian Craton, and South Australian Craton is outlined in red.The dotted gray box in panel (b) indicates the study area shown in Figures 2 and 3.The small black arrow represents the absolute plate motion vector of the Australian plate from Kreemer et al. (2014).

Figure 2 .
Figure 2. Map showing the distribution of (a) seismic stations, and (b) earthquakes used in this study.The dashed black line in panel (a) indicates the approximate position of the Tasman Line for this region.The outline of the Gawler Craton and Curnamona (Cu) province are shown in gray.The event map presented in panel (b) shows the typical event distribution at one station, OOD, the location of which is highlighted in panel (a) in cyan.Events are plotted as open circles color-coded by those that produced a splitting measurement (green) and those that did not (black).

Figure 3 .
Figure 3. Stacked *KS splitting results for stations across the study region plotted against (a) surface elevation and (b) an estimate of lithospheric thickness from AuSREM (Kennett et al., 2012).The stacked splitting parameters for each station are represented by colored bars, orientated according to the fast direction, and scaled in length by the delay time.In (a) the bars are colored according to the fast direction, as indicated by the corresponding color bar.In (b) colors represent the misfit between the fast direction and the absolute plate motion.The locations of all stations analyzed in this study are indicated by the small black dots.White dots indicate the location of stations along the BILBY transect previously analyzed by Eakin et al. (2021).

Figure 4 .
Figure 4. Graph illustrating the relationship between splitting parameters and lithospheric thickness, using the data illustrated in Figure3b.Median values for each bin are indicated by pink circles.The five blue dots in panel (b) correspond to the cluster of stations in southeast Australia in Figure3bwhere results were plotted in the same shade of blue.

Figure 5 .
Figure 5. Results of mantle flow simulations at (a) 120 and (b) 200 km depth.Stacked *KS splitting results (same as Figure3) are plotted in yellow and compared to mantle viscosity and flow vectors (tangential velocities indicated by small black arrows).Further details for the mantle flow simulation are provided in Section 3.3, as well as plots of viscosity and radial velocity at various upper mantle depths in FigureS5in Supporting Information S1.

Figure 6 .
Figure 6.Comparison of stacked *KS splitting results (blue bars) with a simplified geological map of the Gawler Craton, modified from Hand et al. (2007).Location of crustal deformation structures are shown (thick dashed black lines) that either formed or reactivated during NW-SE directed compression (as indicated by the large black arrows) during the interval 1,600-1,580 Ma.The similar pattern of Archean-Early Mesoproterozoic faults (thin gray lines), as well as coeval volcanism from the Hiltaba Suite (red polygons), are overlain from Cowley (2006).