Introduction

The Australian continent is a significant part of the (Indo-) Australian plate (Figure 1a), composed of three Archean to Paleo-Mesoproterozoic tectonic blocks (west, east, and north Australian Cratons) that were amalgamated during several Proterozoic orogeneses (Pirajno & Bagas, 2008) (Figure 1b). The eastern continental margin, known as the Tasmanides, is composed of several younger units accreted to the plate during the Phanerozoic (Glen, 2005). Seismic structures beneath the Australian continent have been a subject of investigation for several decades, revealing insights into both crustal and mantle structures. Early studies, such as those by Bolt et al. (1958) and Cleary et al. (1972), suggested that the seismic wavespeed under the eastern Australian continent is slower than in the western and central parts of the plate, providing initial glimpses into the regional variations. However, these early observations were limited by the sparse distribution of seismic stations and low imaging resolution. Significant progress in unraveling the finer details of the seismic landscape occurred with the increased deployment of permanent seismic stations and temporary arrays across the continent, starting from the operation of the SKIPPY portable array in 1993 (van der Hilst et al., 1994).

Regional surface-wave studies successfully image the high wavespeed cratonic keels underneath western and central Australia, with slow anomalies along the eastern Phanerozoic margin (Debayle & Kennett, 2000a;Fishwick et al., 2005Fishwick et al., , 2008;;Fishwick & Rawlinson, 2012;Simons et al., 1999;Yoshizawa & Kennett, 2004;Zielhuis & Hilst, 1996). In addition to seismic wavespeed heterogeneity, several regional studies have imaged of radial and azimuthal anisotropy beneath the Australian continent (Debayle & Kennett, 2000b;Fichtner et al., 2010;Simons et al., 2002;Yoshizawa, 2014;Yoshizawa & Kennett, 2015). These studies suggest complex lateral variations in seismic anisotropy at shallow depths (above >150 km) and reveal a strong correlation between the absolute plate motion direction and the fast-axis orientation around 200-300 km depths (Debayle & Kennett, 2000b;Simons et al., 2002). The simplified tectonic structures of the (Indo-) Australian plate. Red lines denoting plate boundaries according to Bird (2003). Subduction zones are labeled as follows: Ba (Banda), Br (New Britain), Cr (San Cristobal), Hb (New Hebrides), Hk (Hikurangi), Jv (Java), Ke (Kermadec), To (Tonga), Pu (Puysegur), and Su (Sumatra). (b) The main tectonic structures of the Australian continent, modified from Y. Chen et al. (2023). Dashed black lines outline the boundaries of the west, north, and south Australian cratons, while red lines represent the Tasman Line. Tectonic provinces include: AF (Albany-Fraser belt), Ar (Arunta Block), Am (Amadeus basin), Ca (Canning basin), Cp (Capricorn Orogen), Cu (Curnamona Craton), Er (Eromanga basin), Eu (Eucla basin), Ga (Gawler craton), Ge (Georgetown inlier), Ha (Hamersley basin), Ki (Kimberley block), La (Lachlan Orogen), Mc (MacArthur basin), MI (Mt. Isa block), Mu (Musgrave block), NE (New England Orogen), Of (Officer basin), PC (Pine Creek Inlier), Pi (Pilbara craton), Pj (Pinjarra Orogen), T (Tennant Creek block), and Yi (Yilgarn craton).

In addition to continental-scale studies, global surface wave tomography studies (Chang et al., 2016;Y. Chen et al., 2023;Debayle et al., 2005;Gung et al., 2003;Panning & Romanowicz, 2006;Ritsema et al., 2011) have provided comprehensive maps of long-wavelength and deep structures beneath the Australian plate. This region encompasses not only the Australian continent but also parts of the submerged continent Zealandia, as well as New Guinea, and several island arcs situated along convergent plate boundaries. For instance, Gung et al. (2003) observe significantly faster SH radial anisotropy under the Australian continent at 250-400 km depth. Subsequent studies, such as Debayle et al. (2005), note a unique increase in azimuthal anisotropy at 175-300 km depths in the Australian continent relative to other continents, suggesting horizontal shear in the asthenosphere. Chang et al. (2016) contribute to this understanding by finding that the Tonga slab exhibits faster SV radial anisotropy above 660 km, underlain by faster SH radial anisotropy in the lowermost mantle, possibly linked to subductioninduced shear deformation. While investigations employing body wave tomography techniques (Amaru, 2007;Li et al., 2008;Lu et al., 2019;Obayashi et al., 2013) have revealed fine-scale details of subducting slabs around the Australian plate. These studies unveil the stagnation of the Java slab above the 1,000-km discontinuity and the juxtaposition of the Tonga and Kermadec slabs in the lowermost mantle (Fukao & Obayashi, 2013).

Enhanced data coverage has significantly improved seismic imaging quality, yet most of the studies cited above predominantly rely on classical ray theory or finite-frequency theory (F. A. Dahlen et al., 2000;Marquering et al., 1999;Yoshizawa & Kennett, 2004). Advancements in numerical modeling techniques, notably the spectral element method (Komatitsch & Tromp, 1999;Komatitsch & Vilotte, 1998;Peter et al., 2011), coupled with increased computational capabilities, now enable more cost-effective simulation of wave propagation using realistic anisotropic and anelastic Earth models. Full waveform inversion (FWI) is an advanced technique that incorporates the physics of wave propagation to address seismic inverse problems, building upon the foundational principles used in other methods (F. A. Dahlen et al., 2000;Hung et al., 2000;Lailly, 1983;Liu & Gu, 2012;Tarantola, 1984;Tromp et al., 2005). Moreover, adjoint state methods have demonstrated efficacy in integrating the 3-D seismic wave propagation within iterative inversion procedures (Akcelik et al., 2003;Komatitsch et al., 2016;Tromp et al., 2005Tromp et al., , 2008;;Virieux & Operto, 2009). This approach has been successfully applied, including exploration (Gauthier et al., 1986;Mora, 1987;Pratt et al., 1998;Virieux & Operto, 2009) and addressing largescale continental problems (Afanasiev et al., 2016;Bozda et al., 2016;M. Chen et al., 2015;Fichtner et al., 2013;Tape et al., 2009;Zhu et al., 2013;Zhu & Tromp, 2013). Efficient computation of misfit gradients using the adjoint state method (Liu & Tromp, 2006, 2008) allows the utilization of gradient-based optimization methods, including preconditioned conjugate gradient (Fletcher, 1964) and Limited-memory BFGS approach (Matthies & Strang, 1979;Nocedal, 1980). The FWI technique has been applied to construct global Earth models (Bozda et al., 2016;French & Romanowicz, 2014;Lei et al., 2020;Thrastarson et al., 2024) and image regional and continental-scale seismic structures (M. Chen et al., 2015;Fichtner et al., 2009Fichtner et al., , 2010;;Lloyd et al., 2020;Tao et al., 2018;Tape et al., 2009;Tape et al., 2010;Zhu et al., 2012Zhu et al., , 2015Zhu et al., , 2020)). In comparison with previous continental-scale FWI studies for the Australasian region (Fichtner et al., 2009(Fichtner et al., , 2010)), a significant improvement in this study is the extensive use of three-component body and surface waves to jointly constrain shallow and deep upper mantle structures (Moulik & Ekström, 2014;Ritsema et al., 1999;Zhu et al., 2012). This approach enables the imaging of finer-scale structures down to approximately 800 km depths. In addition, we have benefited from significantly improved data coverage over the past decades, providing approximately a tenfold increase in wavepaths (Figure 2c) compared to the data set used in Fichtner et al. (2010). The Australian Passive Seismic Array project (AusArray) has gathered collaboration between the government and academia, seeking to achieve a comprehensive national half-degree data coverage and an updatable 3D national velocity model. Geoscience Australia, in partnership with Australian State and Territory Geological Surveys, academia, and AuScope, has extended the AusArray project to a national scale through the Exploring For The Future program (EFTF). The investment from EFTF has facilitated the integration of available data sets across Australia, resulting in a doubling of the national rate of new data acquisition (Gorbatov et al., 2019;Kennett et al., 2023).

The primary objective of this study is to construct a comprehensive tomography model for radially-anisotropic shear-wavespeed in the crust and upper mantle beneath the Australian continent and its surrounding regions, refined through the incorporation of new data from the newly deployed AusArray. The intent is to gain deeper insights into the evolutionary processes of the area through the application of FWI. Our approach commences with a detailed description of the data set employed and a thorough explanation of the methodology applied in constructing the AU21 model. Subsequently, we present a series of cross-sectional views of AU21, accompanied by a resolution analysis utilizing point-spread functions (PSFs) (Fichtner &Trampert, 2011) tests to assess the reliability and precision of the model. We will examine the details of the shear-wavespeed, and radial anisotropy structures embedded in model AU21. To conclude, we conduct a comparative analysis, aligning our results with previously published models to highlight distinctive features and contribute to the evolving understanding of the subsurface structure beneath the study region.

Data Set and Method

Distribution of Earthquakes and Stations

In this study, we utilize data from 248 earthquakes from 1993 to 2019, with moment magnitudes (Mw) ranging from 4.7 to 7.0 (Figure 2a). The lower bound in the chosen magnitude range is established to ensure a sufficient signal-to-noise ratio in recorded data, while the upper bound is set to exclude seismic sources with complex rupture processes that violate the point source approximation utilized in the modeling. Notably, earthquakes with relatively smaller magnitudes are concentrated along the Southeast Indian Ridge (SEIR) and 90°E Ridge. The considered selection of these seismic events ensures a reasonably uniform source distribution to the south and west of our targeted study region. The data set is recorded by 1,102 seismographic stations covering the Australian plate (Figure 2b), mostly consisting of permanent networks (e.g., AU, NZ, II, IU, G, and GE), as well as published arrays (e.g., networks 6F, 7B, 7G, 7J, 7K, and S1) from the Australian Passive Seismic Server (AUSPASS). In addition, we incorporate data from island stations and ocean bottom seismometers to enhance our understanding of subduction zones surrounding the Australian plate. The corresponding ray coverage map is shown in Figure 2c.

Inversion Strategy

The cycle-skipping is one of the limitations of FWI since its origin (Gauthier et al., 1986). To address this problem (Virieux & Operto, 2009), we start from relatively longer period bands for both body and surface waves, and gradually incorporate higher frequency signals into the inversion. During the initial iterations (1 st to 6 th ), our primary focus is on constructing long-wavelength structures using relatively longer period surface waves (50-150 s) and body waves (20-50 s). In subsequent iterations (7 th to 12 th ), we gradually reduce the short-period corner of surface waves from 50 to 40 s, simultaneously adjusting the short-period corner for body waves from 20 to 18 s. As we advance to iterations 13 th to 16 th , the frequency bands for body and surface waves are further refined to 16-45 s and 32-125 s, respectively. In the final iterations (17 th to 21 st ), observing that misfit evolution curves become relatively flat, we decide to apply 25-100 s surface waves and 15-40 s body waves to update fine-scale features.

The selection of an initial model is critical in FWI. Several studies have addressed this challenge by choosing starting models and utilizing appropriate measurements to mitigate the risk of convergence to local minima (Brossier et al., 2009;Pratt & Shipp, 1999;Prieux et al., 2013;Y. O. Yuan & Simons, 2014;Y. O. Yuan et al., 2015). In addition, nonlinearities can be mitigated by commencing with smooth models and low-frequency signals, systematically integrating higher frequency content in subsequent iterations (Nolet et al., 1986;Zhu et al., 2012). The one-dimensional radial structure of the Earth, particularly the extensive long-wavelength shearwave-speed structure in the mantle, is widely acknowledged (Becker & Boschi, 2002;Ritzwoller & Lavely, 1995). Current iterative inversions which start from radially symmetric models support this consensus (Lekić & Romanowicz, 2011). Moreover, widely recognized global 3-D crustal models such as the 2 °× 2°Crust2.0 (Bassin, 2000) and its subsequent version Crust1.0 (Laske et al., 2013), featuring a resolution of 1 °× 1°, are readily accessible. In this study, we employed a coupled 3-D mantle model, S362ANI (Kustowski et al., 2008), along with the Crust2.0 (Bassin, 2000). S362ANI is a global radially anisotropic model which was constructed using surface-wave phase speeds, body-wave traveltimes, and long-period body and mantle waveforms. The integrated model, S362ANI + Crust2.0, achieves an adequate fit to body and surface waves with low frequencies and represents compelling advancements over a 1-D model (Bozdağ et al., 2016;Tromp et al., 2010;Zhu et al., 2015).

We discretize the simulation domain using 284 × 192 elements, with lateral dimensions of 102 °× 67°. The average size of one element on the Earth's surface is approximately 40 km, and the minimum period resolved by this mesh is about 15 s. Employing two computing nodes with 96 cores on the Lonestar5 cluster at the Texas Advanced Computing Center (TACC), it takes approximately 1 hr for a 30-min forward simulation and 3.5 hr for an adjoint simulation to construct the misfit gradients. Iterative updates to the transversely isotropic model involve parametrization in terms of compressional and shear-wave speeds, along with a dimensionless parameter, while density perturbations are scaled to isotropic shear-wave-speed perturbations (Additional information on the misfit function and model parametrizations can be found in Supporting Information S1 file). The utilization of this transversely isotropic parametrization plays a crucial role in addressing the Rayleigh/Love discrepancy (Dziewonski & Anderson, 1981).

Improvements of Travel Time Histograms and Misfit Functions

The improvements in the six-component travel time histogram for both short-period (15-40 s) body waves and long-period (25-100 s) surface waves are shown in Figure 3. These travel times are computed using crosscorrelations between observed and predicted seismograms with windows selected by FLEXWIN (Maggi et al., 2009). All six-component histograms for the final model AU21 show significant improvements compared to the initial S362ANI model. For example, the mean value for Rayleigh waves on the radial component has been reduced from 3.96 to 0.05, and the standard deviation is decreased from 4.01 to 3.13. These enhancements demonstrate that synthetic seismograms generated from the final AU21 model simultaneously match well with both observed body and surface waves. Figure 4 illustrates the evolution curves for the six sub-categories, along with the curve for the total misfit, defined as the average of all sub-categories. It is observed that misfit values decrease significantly after the first three iterations and gradually flatten. We ensure that misfits are decreasing within each stage, and it is normal to experience abrupt jumps in misfit values when the frequency band is changed at different stages. The total number of body-wave measurements has increased from 17,549 to 32,655 after 21 iterations, while the number of surfacewave measurements has increased from 30,431 to 35,897.

Comparisons of Three-Component Waveforms

Figures 5 and 6 compare three-component observed (black) and synthetic (red) seismograms for two earthquakes. The first event is an Mw 5.8 earthquake that occurred at a depth of 20.7 km beneath western Australia (CMT201812161426A). We apply a bandpass filter (25-100 s) to both observed and synthetic data to compare long-period surface waves. The second event is an Mw 6.2 earthquake that occurred at 603 km depth beneath Fiji (CMT201901261956A). We bandpass filter both observed and predicted seismograms from 15 to 40 s for the second event. Theoretical arrival times of P, S, and ScS phases from the PREM model (Dziewonski & Anderson, 1981) are indicated in the seismograms. In addition, we present comparisons for both long-period (25-100 s) and short-period (15-40 s) waveforms from another deep earthquake (CMT201608310311A) occurred beneath the Solomon Islands Figures S1-S2 in Supporting Information S1. All these results show that long-period Rayleigh waves in vertical and radial components, and Love wave in transverse component are well-matched in phases, and the waveform fittings around P and S arrivals in short-period data are also reasonably good. These data comparisons suggest that β v and β h models in AU21 are well-constrained. However, significant deviations of the horizontal component from the N-S direction have been reported in some stations in Australia (e.g., AU. STKA) (Eakin et al., 2023). While we did not apply corrections for such deviations, our evaluation indicates that these uncorrected deviations might introduce localized inaccuracies in waveform fits but do not significantly compromise the overall model accuracy. Journal of Geophysical Research: Solid Earth 10.1029/2024JB029260

Approximate Hessian and Resolution Analysis

Evaluating the reliability of tomographic models is important for seismic tomography. We initiate this process by analyzing seismic illumination, computing the approximated diagonal Hessian using the formulation in Luo et al. ( 2013): (Dziewonski & Anderson, 1981). Z and R component fits worsen at far epicentral distances and late times, likely influenced by coda waves.

where s(x, t and s (x, Tt) denote the forward and adjoint wavefields, respectively.

In Figure 7, horizontal depth slices of the approximated diagonal Hessian at depths ranging from 100 to 800 km are shown. The results indicate that illumination is notably good at depths shallower than 450 km beneath major continental areas, however, it starts to diminish at 550 km depth beneath western Australia. Depths greater than Journal of Geophysical Research: Solid Earth The PSFs is a widely employed tool for evaluating resolution and inter-parameter trade-offs (Fichtner & Trampert, 2011;Zhu et al., 2015), which can be approximated as the difference between the gradients of the current model m and a perturbed model m + δm.

Hδm ≈ g(m + δm) -g(m).

(

where H denotes the Hessian and g denotes the misfit gradients of any model parameter.

However, implementing this method for testing numerous locations in the inverted model is not feasible due to the computational cost, which is equivalent to one FWI iteration per PSF. To address this challenge and efficiently evaluate PSFs at multiple locations, we adopt a strategy inspired by Rawlinson and Spakman (2016), as well as Tao et al. (2018). Here, we input sparsely distributed Gaussian perturbations and compute all the PSFs simultaneously. Generally, anisotropic materials are characterized by a fourth-order elastic tensor (c ijkl ) involving 21 independent elements. For materials with radial symmetry axis, this complexity is reduced to five primary parameters: A, C, L, N, and F, known as the Love parameters (Love, 2013). Since S362ANI is based on six key parameters: mass density (ρ), two compressional wave speeds (α v and α h ), two shear wave speeds (β v and β h ), and a dimensionless parameter (η), their relationship with the Love parameters are used (F. Dahlen & Tromp, 2020) (Equations S2-S6 in Supporting Information S1).

Figure 8 shows the resolution test results beneath the Australian continent. Each column represents input Gaussian perturbations on model parameters N (shear modulus with σ horizontal = 130 km and σ vertical = 40 km), and the PSFs for parameters N, L (bulk modulus), A (P-wave speed), and C (density). Notably, there is a robust recovery of input perturbations beneath the Australian continent down to 600 km, and inter-parameter leakages from N, L, A, and C are quite weak. However, at 800 km depth, the PSFs for parameter N become less distinct, indicating reduced resolution at this depth.

We conducted similar resolution tests at various locations, including New Zealand, the Tonga subduction zone, the Coral and Tasman Sea, and the SEIR (Figures S3-S6 in Supporting Information S1). While the wave path coverage in these regions is not as extensive as that beneath the Australian continent, a significant portion of input Gaussian perturbations is reasonably well recovered, but the PSFs beneath the Samoan plume are missing. It is essential to highlight that the PSFs for parameter N beneath Fiji are smearing along the NE-SW direction, leading to recovered perturbations larger than the inputs, especially at 200 km depth. These results suggest a potential overestimation of the size of the Tonga slab in model AU21. Journal of Geophysical Research: Solid Earth 10.1029/2024JB029260 BODUR ET AL.

Results

We present model AU21 along several horizontal and vertical cross-sections, where Figures 9 and 10 It is important to note that radial anisotropy can also result from the alignment of partially molten layers or finescale stratifications within the mantle (Aki, 1968;Babuska & Cara, 1991). However, in our interpretation of ξ, we primarily consider mantle flow directions due to their significant influence on seismic anisotropy and our interest to investigate the dynamic processes. Journal of Geophysical Research: Solid Earth 10.1029/2024JB029260

The Australasian Region

The dominant feature in the shear-wavespeed model is the presence of relatively large fast anomalies (up to > +7%) beneath the Precambrian regions of western and central Australia for depths exceeding 200 km (Figures 9a-9c). Contrastingly, the easternmost parts of Australia and the adjacent Tasman and Coral Seas are dominated by band-shaped slow anomalies, except several weak fast anomalies are embedded at 100 km depth beneath the Lord Howe Rise region, which might be related to the submerged plateau in these regions (Figure 9a). The location and extent of these anomalies align with previous surface wave studies (Debayle & Kennett, 2000a, 2000b;Fichtner et al., 2009;Fishwick et al., 2005Fishwick et al., , 2008;;Simons et al., 1999;Yoshizawa, 2014;Yoshizawa & Kennett, 2004). The northern extent of the fast anomalies reaches New Guinea and Timor Island, marking the forefront of the continental collision between Australia and Southeast Asia. Around these depths, we do not observe a prominent dichotomy in radial anisotropy beneath the Australian continent and surrounding oceans, furthermore, long-wavelength faster SH radial anisotropy (ξ > 1.05) dominates in most study regions (Figures 10a-10c). Fast anomalies beneath western and central Australia diminish at depths of 250-300 km, suggesting the transition from the lithosphere to the asthenosphere (Figures 9d and 9e), while slow anomalies persist beneath the Tasman and Coral Seas. Around these depths, the Western Australian Craton (WAC) region changes to ξ > 1, while southeast Australia and east of New Guinea are characterized as ξ < 1 and collocated with slow wavespeed perturbations (Figures 10d and 10e). Journal of Geophysical Research: Solid Earth Journal of Geophysical Research: Solid Earth 10.1029/2024JB029260

At greater depths, large fast anomalies (+2>3%) emerge at 450 km underneath the Australian continent and New Guinea, extending down to 700 km. This phenomenon may be related to the fossil Welford slab (van der Meer et al., 2018) from past subduction during 10-20 Ma (Schellart & Spakman, 2015). Slow anomalies are observed beneath the east of New Guinea and the Coral Sea down to 800 km depth, which is also indicated by French and Romanowicz (2014). Another group of slow anomalies shows up beneath Tasmania and the Tasman Sea down to greater depths, which might be associated with the thermal origin of the Tasmanid hotspot (McDougall & Duncan, 1988;Sutherland, 1983;Wellman, 1983). Beneath 350 km depth, the Australian continent is observed with laterally varying ξ > 1 and ξ < 1 areas, but faster SV radial anisotropy seems to be dominant (Figures 10f-10l).

Adjacent Subduction Zones

At the Tonga-Kermadec subduction zone (TKSZ), the Tonga slab is seen as a strong, fast anomaly (>+4%) at a depth of 200 km, continuing down to 700 km (Figure 9c). The Kermadec slab is less noticeable above 400 km, but a band-shaped fast anomaly shows up below 450 km depth (Figure 9h). At 200 km depth, the Tonga slab appears as a "dot" due to its steep subduction angle and narrow profile, reflecting the intense tectonic processes in the region. The Tonga slab is stagnant within the transition zone, while the Kermadec slab penetrates through the 660km discontinuity without pooling in the transition zone. Slab-shaped fast anomalies generally correlate with ξ < 1 above 400 km depth. However, at some places in the mantle transition zone, fast anomalies are collocated with ξ > 1. As for the Samoan plume, it is illuminated down to 600 km depth, with remarkably slow shear-wavespeed perturbations (> -2%), and generally ξ > 1 radial anisotropy.

Stripe-shaped fast anomalies emerge at the Java subduction zone above 400 km depth, showing ξ < 1 and ξ > 1 in the mantle wedge and sub-slab regions, respectively. Beyond the 400 km, there is a significant incensement in the extent of fast anomalies, hinting at potential slab stagnation within the mantle transition zone. However, it is worth noting to consider the potential smearing artifacts due to inherent limitations in resolution beneath the Java subduction zone at greater depths (Figure S4 in Supporting Information S1).

Another significant slab-shaped fast anomaly is observed reaching the 660 km beneath New Guinea (Figure 12g), possibly corresponding to the subducted New Britain slab or the Solomon Sea plate lithosphere (Bijwaard et al., 1998;Hall & Spakman, 2002). The curvature of fast perturbations and seismicity suggests the New Britain slab bends backward as it goes through the mantle transition zone. Potential slab stagnation may occur, indicated by the enlarged volume of fast anomalies above the 660.

The Southeast Indian Ridge

Slow anomalies are observed along the SEIR at depths below 150 km (Figures 9a and 9b) due to the thermal conditions around the mid-ocean ridges. In addition, a slightly fast anomaly is found near the Australia-Antarctica Discordance (AAD), consistent with the findings of Kuo et al. (1996). Large fast anomalies at depths exceeding 450 km are feasible, possibly originating from an ancient subduction zone that retreated across the Tethys Ocean (Gurnis et al., 1998;Hall & Spakman, 2002). The suggested average subduction speed of the Mesozoic Tethyan oceanic subduction is approximately 2 cm/yr, a rate that would have reached the core-mantle boundary (CMB), making it impossible to detect in current mantle tomography (van der Meer et al., 2018). However, long-term slab stagnation in the transition zone has been imaged in multiple regions (e.g. (Fukao et al., 2009;van der Hilst & Seno, 1993),), furthermore, this stagnation is visible in vertical cross-sections (Figure 12a), consistent with Simmons et al. (2015).

Radial anisotropy suggests a change in flow direction from horizontal to vertical across depths of 150-200 km, assuming that ξ > 1 corresponds to horizontal flow, and ξ < 1 corresponds to vertical flow (H. Yuan et al., 2011;Zhu et al., 2015). Considering the expectation of horizontal-dominated flow near mid-ocean ridges and verticaldominated flow near subducting slabs, these observations align with the transition between these two features at the displayed depth in vertical cross-section (Figure 12a).

Journal of Geophysical Research: Solid Earth

10.1029/2024JB029260

Discussion

In this section, our goal is to identify important features in shear-wavespeed perturbations and radial anisotropy by comparing our model AU21 with previously published models. Subsequently, we focus on interpreting isotropic shear-wavespeed perturbations and radial anisotropy beneath the Australian continent.

Comparisons With Previous Models

To compare seismic structures at depths shallower than 300 km, three regional tomographic models were selected: AuSREM (Kennett et al., 2013), CSEM (Fichtner et al., 2018), andSimons et al. (2002) (referred to as Simons02). These models compute average shear-wavespeeds as well as wavespeed perturbations with respect to the average values (see Figure 13). At depths of 100-170 km, all three models reveal consistent long-wavelength structures for shear-wavespeed perturbations. Fast anomalies are observed for the cratonic lithosphere, while slow anomalies are delineated under the eastern continent margin, Coral Sea, and Tasman Sea. However, slight discrepancies exist in the location of the boundary between Precambrian and Phanerozoic Australia. Furthermore, significant differences arise in the location and shape of the craton keel at depths of 250 and 300 km, with notable variations in the magnitude of fast anomalies across the models. Regarding radial anisotropy (see Figure 14), shared characteristics among the models include a transition from ξ > 1 to ξ < 1 beneath southeast Australia from (Fichtner et al., 2018), AuSREM (Kennett et al., 2013), and Simons02 (Simons et al., 2002) at depths ranging from 80 to 300 km (top to bottom).

80 to 300 km depth and ξ ≈ 1 beneath the WAC at approximately 170 km depth. However, inconsistencies arise, such as AU21 showing ξ > 1 beneath the Java subduction zone and the west of the South Island, New Zealand at certain depths, which are not observed in other models.

At greater depths, three radially anisotropic shear-wave models were compared: SPani (Tesoniero et al., 2015), SAVANI (Auer et al., 2014), and SGLOBE_rani (Chang et al., 2015) (Figure 15). Long-wavelength fast wave speed perturbations are observed along convergent plate boundaries, correlating with subducting plates. Discrepancies in radial anisotropy among these models become more apparent, with variations in the polarity of ξ at different depths and regions (Figure 16). Specifically, AU21 shows certain features, such as ξ > 1 beneath the Tonga slab at depths greater than 400 km, which differ from other models that display varying anisotropy orientations. Such differences highlight the unique aspects of our model but also underscore the complexity and variability in interpreting subsurface structures. Inconsistencies may be related to methodological differences and resolution or geological complexity with prior studies are observed. To assess short-wavelength structures beneath the TKSZ in model AU21, five published global P-wave travel-time tomography models were collected: GAP_P4 (Fukao & Obayashi, 2013), MITP08 (Li et al., 2008), UU-P07 (Amaru, 2007), and TX2019P (Lu et al., 2019). Although P-wave and S-wave velocities reflect different elastic properties-bulk and/or shear moduli-and can be influenced by factors such as water content and state of stress, this comparison provides a comprehensive understanding of the subsurface structures. Despite these differences, valuable insights into geological and tectonic features can still be obtained. In addition, a P-wave model was constructed (Figure S7 in Supporting Information S1) to consider the differences and similarities in velocity perturbations between P and S waves, ensuring a robust interpretation. Figure 17 illustrates wavespeed perturbations for three profiles across the TKSZ, showing significant differences in results beneath the 660 among these models. Notably, while all models suggest the Kermadec slab penetrates through the 660 and changes to shallow dipping below it, differences arise in the illustration of fast anomalies pooling above the 1,000 km depth across profile L-l, possibly attributed to degraded resolution at greater depths in AU21.

Given these observations, it is important to acknowledge the inherent uncertainties and limitations in our resolution assessment. Point-spread functions and similar methods offer valuable insights but are not capable of fully capturing the complexities and uncertainties of the Earth's subsurface structures. Therefore, careful interpretation of the differences observed between various models is warranted. While our model AU21 provides significant advancements and detailed insights, it is essential to consider these results within the broader context of existing models and their respective limitations. (Tesoniero et al., 2015), SAVANI (Auer et al., 2014), and SGLOBE_rani (Chang et al., 2015) at depths ranging from 400 to 700 km (top to bottom).

Interpretation of Seismic Structures

In our analysis of seismic wavespeed variations, we discern a noteworthy discrepancy in lithospheric thickness between major Archean cratonic regions and the Proterozoic central Australia. Specifically, our findings indicate a thinner lithosphere beneath Archean cratons, contrasting with the more substantial lithospheric thickness observed in Proterozoic central Australia. This discrepancy may be attributed to complex geological processes occurring over extended geological time scales. Subcontinental downwelling, for instance, could contribute to the thinning of the lithosphere beneath Archean cratons, while erosions of cratonic roots over time may also play a significant role in shaping the observed lithospheric variations (Fichtner et al., 2010;Gung et al., 2003;Simons et al., 1999). In addition, the viscosity contrast between the stable cratons and the underlying asthenosphere, as well as the influence of basal drag, are critical factors that further clarify the dynamic interactions between these layers (Yoshida & Yoshizawa, 2021).

Within the Archean cratonic regions, our analysis reveals a noteworthy pattern of radial anisotropy, particularly between depths of 80-150 km. We observe a substantial decrease in radial anisotropy within this depth range, suggesting the presence of compositional variations and deformation processes during continental mantle formation. The observed decrease in radial anisotropy hints at potential complexities in the underlying mantle structure, possibly influenced by the tectonic history and geological evolution of the region. These findings provide valuable insights into the dynamic processes shaping the continental lithosphere (Fichtner et al., 2010;Yoshizawa, 2014). The weak radial anisotropy observed at these depths beneath the WAC can be attributed to the alignment of dipping layers or the tilted symmetry axes of anisotropic minerals. This interpretation is supported by similar findings in other cratonic regions, where seismic anisotropy is often influenced by the structural alignment of minerals such as olivine (Karato et al., 2008). The dipping layers could represent ancient tectonic Journal of Geophysical Research: Solid Earth 10.1029/2024JB029260

features, such as strike-slip faults or rifting processes, consistent with the tectonic history of Western Australia (Cassidy et al., 2006;Czarnota et al., 2010). Alternatively, the tilted symmetry axes of anisotropic minerals could result from deformation processes related to past tectonic events. These geological features lead to variations in seismic wave propagation velocities, contributing to the observed anisotropic patterns.

At depths ranging from 150 to 200 km, our analysis uncovers moderate faster SH anisotropy, coinciding with the lithosphere-asthenosphere transition. The Lithosphere-Asthenosphere boundary depth was estimated using shear Journal of Geophysical Research: Solid Earth 10.1029/2024JB029260

wavespeed perturbation data. The base of faster anomalies in the velocity perturbation profiles, typically positioned at depths of 250-300 km, was identified to define the LAB, indicating the transition from the more rigid lithosphere to the more ductile asthenosphere. This observation suggests significant horizontal shear deformation occurring at the base of the lithosphere, potentially influenced by the basal drag of the Australian plate. The presence of moderate faster SH anisotropy within this depth range highlights the dynamic interplay between tectonic forces and mantle deformation processes, underscoring the complex nature of lithospheric dynamics in the region (Debayle & Kennett, 2000a;Simons et al., 2002) This interpretation supports the detailed discussions by Yoshizawa and Kennett (2015) regarding the role of horizontal shear deformation in defining the lithosphereasthenosphere transition.

In comparing our FWI model with Yoshizawa (2014)'s radially anisotropic shear wave model, several distinctions in the characterization of the LAB become apparent. Yoshizawa (2014), which utilizes multi-mode surface wave tomography, provides a detailed picture of LAB variations across the Australian Plate by examining the vertical gradient of shear wave speed. This approach highlights a gradual transition from the lithosphere to the asthenosphere, which is particularly thick beneath the Archean cratons and comparatively thinner under Proterozoic regions. In contrast, our FWI model incorporates finite-frequency effects, allowing for higher-resolution imaging that reveals sharper velocity gradients across the LAB. This difference suggests that our model captures finer structural details, which may provide new insights into the precise depth and characteristics of the LAB beneath Australia. By contrasting the results of our model with those from Yoshizawa (2014), we can further assess the regional variability of the LAB, potentially indicating differing lithospheric compositions and tectonic histories within Australia. To illustrate, in central Australia, our FWI results estimate the LAB at approximately 250 km depth, where a sharp transition in wavespeed is observed. Yoshizawa's surface wave inversion results, however, display a smoother transition in this region, making it challenging to pinpoint the LAB depth as precisely. This example highlights the advantage of FWI in resolving finer-scale features of the lithospheric structure that are less discernible in traditional surface wave inversion models.

As we explore deeper into the mantle, our analysis reveals intriguing patterns of radial anisotropy below 250 km depth. Here, radial anisotropy diminishes, possibly indicative of the transition from dislocation to diffusion creep or olivine lattice-preferred orientation (LPO) transition. These findings underscore the intricate nature of mantle deformation mechanisms operating at greater depths and highlight the dynamic interplay between various rheological processes shaping the underlying mantle structure (Fichtner et al., 2010;Yoshizawa, 2014).

The integration of observations from Eastern Australia further enriches our understanding of regional seismic characteristics. For instance, the correlation between high heat flux regions and observed seismic anomalies suggests the ascent of mantle materials from the asthenosphere, potentially related to cryptic mantle plumes.

Comparison with the Australian mantle temperature model by (Goes et al., 2005), estimated from earlier tomography models, further supports these observations, indicating that regions with high heat flux correspond to elevated mantle temperatures, highlighting the concepts of significant regional thermal activity beneath Eastern Australia. Recent studies on Cenozoic volcanism in eastern Australia (Davies et al., 2015;Rawlinson et al., 2017) have shown that variations in lithospheric thickness can lead to vigorous mantle upwellings, promoting decompression melting. This process likely contributes to the seismic anomalies we observe, particularly beneath Tasmania and the Tasman Sea, where the lithosphere is thinner, facilitating the ascent of mantle materials and subsequent volcanic activity. The presence of slow anomalies beneath specific regions, such as Tasmania and the Tasman Sea, may be attributed to putative plume sites, indicating complex mantle dynamics beneath Eastern Australia (Kennett & Davies, 2020;Morgan & Morgan, 2007).

In the interpretation of radial anisotropy presented above, we implicitly rely on the relation derived from studies on olivine under simple shear deformation (Zhang & Karato, 1995). However, it is crucial to note that this assumption is not universally valid for all olivine LPO types. Jung and Karato (2001) demonstrate that the transition from A-type to B-or C-type LPO can result in complex radial anisotropy patterns, highlighting the need for careful consideration of mantle conditions in interpreting seismic observations.

Conclusions

We implemented a continental-scale full-waveform inversion to investigate radially anisotropic shear-wavespeed structures in the upper mantle beneath the Australian plate and surrounding regions. This involved analyzing three-component data generated by 248 regional earthquakes and recorded by 1,102 seismic stations, both

Journal of Geophysical Research: Solid Earth 10.1029/2024JB029260 permanent and temporary installations. After 21 preconditioned conjugate gradient iterations, we observed a significant reduction in misfits and improvements in the histograms of travel-time differences. Simultaneously, the waveform fitting between observed data and synthetics using the AU21 model demonstrated exceptional accuracy. Resolution tests based on PSFs indicated that a majority of features in AU21 are well-resolved. Our findings lead to the following conclusions:

1. We identified sharp contrasts in shear-wavespeed perturbations between Phanerozoic and Precambrian Australia above >200 km depth. This boundary is located to the east of the Tasman Line, estimated based on outcrop age. The estimated LAB from shear-wavespeed perturbations beneath western and central Australia is positioned at 250-300 km. 2. A layer with weak radial anisotropy at 80-150 km depth is observed beneath the WAC, potentially resulting from alignments of dipping layers or tilted axes of symmetry of anisotropic minerals. This observation may have potential implications for understanding the formation processes of the Archean craton. Moreover, between depths of 150-250 km, the presence of faster SH radial anisotropy suggests the occurrence of relatively significant horizontal strain during the lithosphere-to-asthenosphere transition. 3. Slow anomalies, accompanied by faster SV radial anisotropy and a high Vp/Vs ratio, extend into the uppermost lower mantle beneath the eastern regions of New Guinea, the Tasman Island, and the Tasman Sea. This observation suggests that deep thermal activities contribute to the formation of the low wavespeed band at the eastern continent margin. 4. Several slabs are imaged beneath the surrounding subduction zones. The Tonga slab remains stagnant within the mantle transition zone, while the Kermadec slab penetrates through the 660-km discontinuity. The Java slab appears to bend backward; however, the resolution is degraded beneath the Java subduction zone at depths exceeding 600 km.