Crustal structure of southeast Australia from teleseismic receiver functions

In an effort to improve our understanding of the seismic character of the crust beneath southeast Australia and how it relates to the tectonic evolution of the region, we analyse teleseismic earthquakes recorded by 24 temporary and 8 permanent broadband stations using the receiver function method. Due to the proximity of the temporary stations to Bass Strait, only 13 of these stations yielded usable receiver functions, whereas seven permanent stations produced receiver functions for subsequent analysis. Crustal thickness, bulk seismic velocity properties, and internal crustal structure of the southern Tasmanides – an assemblage of Palaeozoic accretionary orogens that occupy eastern Australia – are constrained by H–κ stacking and receiver function inversion, which point to the following:

1. a ∼ 39.0 km thick crust; an intermediate–high $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratio (∼ 1.70–1.76), relative to ak135; and a broad (> 10 km) crust–mantle transition beneath the Lachlan Fold Belt. These results are interpreted to represent magmatic underplating of mafic materials at the base of the crust.

2. a complex crustal structure beneath VanDieland, a putative Precambrian continental fragment embedded in the southernmost Tasmanides, that features strong variability in the crustal thickness (23–37 km) and $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratio (1.65–193), the latter of which likely represents compositional variability and the presence of melt. The complex origins of VanDieland, which comprises multiple continental ribbons, coupled with recent failed rifting and intraplate volcanism, likely contributes to these observations.

3. stations located in the East Tasmania Terrane and eastern Bass Strait (ETT + EB) collectively indicate a crust of uniform thickness (31–32 km), which clearly distinguishes it from VanDieland to the west.

Moho depths are also compared with the continent-wide AusMoho model in southeast Australia and are shown to be largely consistent, except in regions where AusMoho has few constraints (e.g. Flinders Island). A joint interpretation of the new results with ambient noise, teleseismic tomography, and teleseismic shear wave splitting anisotropy helps provide new insight into the way that the crust has been shaped by recent events, including failed rifting during the break-up of Australia and Antarctica and recent intraplate volcanism.

1 Introduction

The Phanerozoic Tasmanides (Collins and Vernon, 1994; Coney, 1995; Coney et al., 1990) comprise the eastern third of the Australian continent and, through a process of subduction accretion, were juxtaposed against the eastern flank of the Precambrian shield region of Australia beginning in the late Neoproterozoic and early Palaeozoic (Foster and Gray, 2000; Glen, 2005; Glen et al., 2009; Moresi et al., 2014) (Fig. 1). Persistent sources of debate that impede a more complete understanding of the geology of the Tasmanides include (1) the geological link between Tasmania – an island state in southeast Australia – and mainland Australia, which are separated by the waters of Bass Strait; and (2) the presence and locations of continental fragments from Rodinian remnants that are entrained within the accretionary orogens. Furthermore, the lateral boundaries between individual tectonic blocks and their crustal structure are often not well defined. To date, few constraints on crustal thickness and seismic velocity structure have been available for regions such as Bass Strait. Therefore, constraints on the Moho transition, crustal thickness, and velocity structure beneath Bass Strait derived from receiver functions (RFs) can provide fresh insight into the nature and evolution of the Tasmanides.

Figure 1 Regional map of southeastern Australia that shows key geological boundaries and the location of observed or inferred tectonic units (modified from Bello et al., 2019a). Thick black lines delineate structural boundaries, and the thick brown dashed line traces out the boundary of VanDieland. The following abbreviations are used in the figure: HF – Heathcote Fault; GF – Governor Fault; BF – Bootheragandra Fault; KF – Koonenberry Fault; THZ – Torrens Hinge Zone; MA – Macquarie Arc; NVP – Newer Volcanics Province; KI – King Island in Bass Strait; FI – Flinders Island in Bass Strait; WTT – West Tasmania Terrane; ETT – East Tasmania Terrane; AL – Arthur Lineament; TFS – Tamar Fracture System; and RCB – Rocky Cape Block. Outcrop boundaries are sourced from Rawlinson et al. (2016).

Previous estimates of crustal thickness and structure beneath southeastern Australia have been obtained from deep seismic reflection transects, wide-angle seismic data, topography, and gravity anomalies (e.g. Collins, 1991; Collins et al., 2003; Drummond et al., 2006; Kennett et al., 2011). Earlier RF studies in southeast Australia (Shibutani et al., 1996; Clitheroe et al., 2000; Tkalčić et al., 2011; Fontaine et al., 2013a, b) suggested the presence of complex lateral velocity variations in the mid-lower crust that probably reflect the interaction of igneous underplating, associated thinning of the lithosphere, recent hotspot volcanism, and uplift. Furthermore, the intermediate to high crustal $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratio of 1.70–1.78 in this region (Fontaine et al., 2013a), relative to ak135 continental crust where $V_{\mathrm{p}}/V_{\mathrm{s}}$ is ∼1.68, may indicate a mafic composition that includes mafic granulite rocks, granite gneiss, and biotite gneiss. Body- and surface-wave tomography (Fishwick and Rawlinson, 2012; Rawlinson et al., 2015) revealed P- and S-wave velocity anomalies in the uppermost mantle beneath Bass Strait and the Lachlan Fold Belt. Ambient noise surface wave tomography (Bodin et al., 2012b; Young et al., 2012; Pilia et al., 2015b, 2016; Crowder et al., 2019) of the southern Tasmanides revealed significant crustal complexity, but it is unable to constrain crustal thickness or the nature of the Moho transition.

The goal of this paper is to provide fresh insight into the crust and Moho structure beneath the southern Tasmanides using P-wave RFs and to explain the origin of the lateral heterogeneities that are observed. This will allow us to explore the geological relationship between the different tectonic units that constitute the southern Tasmanides and to develop an improved understanding of the region's tectonic history.

2 Geological setting

The Palaeozoic–Mesozoic Tasmanides of eastern Australia form part of one of the most extensive accretionary orogens in existence and evolved from interaction between the East Gondwana margin and the proto-Pacific Ocean. The tectonic evolution of the Tasmanides is complex, and large-scale reconstructions have proven difficult. This is evident from the variety of models that have been suggested to explain how the region formed (Foster and Gray, 2000; Spaggiari et al., 2003; Teasdale et al., 2003; Spaggiari et al., 2004; Boger and Miller, 2004; Glen, 2005; Cawood, 2005; Glen et al., 2009; Cayley, 2011a, b; Gibson et al., 2011; Moresi et al., 2014; Pilia et al., 2015a, b). Particular challenges arise from multiple subduction events, multiple phases of metamorphism, entrainment of exotic continental blocks, the formation of large oroclines, recent intraplate volcanism, and subsequent events, including the separation of Antarctica and Australia and the formation of the Tasman Sea. These challenges are compounded by the presence of widespread sedimentary sequences that hinder direct access to basement rocks (Fig. 1).

The Tasmanides consist of four orogenic belts, namely the Delamerian, Lachlan, Thomson, and New England orogens. The Delamerian Orogen – located in the south – is the oldest part of the Tasmanides and has a southward extension across Bass Strait from Victoria into western Tasmania, where it is commonly referred to as the Tyennan Orogen (Berry et al., 2008). Between about 514 and 490 Ma, the Precambrian and early Cambrian rocks that constitute the Delamerian Orogen were subjected to a contractional orogenic event along the margin of East Gondwana (Foden et al., 2006). Subsequently, the Lachlan Orogen formed in the east, which contains rocks that vary in age from Ordovician to Carboniferous (Glen, 2005). Gray and Foster (2004) argued for a tectonic model of the Lachlan Orogen that involved the interaction of a volcanic arc, oceanic microplates, and three distinct subduction zones. Each subduction zone is linked to the formation of a distinct tectonic terrain: the Stawell–Bendigo zone, the Tabbarebbera zone, and the Narooma accretionary complex. The limited rock exposure in the Tasmanides as a whole has made direct observation of the Lachlan Orogen difficult; this is attributed to a large swath of Mesozoic–Cenozoic sedimentary cover and more recent Quaternary volcanics, which obscure a large portion of the underlying Palaeozoic terrane. However, the Lachlan Orogen contains belts of Cambrian rocks in Victoria and New South Wales that are similar in age to the Delamerian Orogen (Gray and Foster, 2004).

The presence of Precambrian outcrops in Tasmania and the relative lack of rocks that are similar in age in adjacent mainland Australia has led to different models which attempted to explain the existence of Proterozoic Tasmania. For instance, Li et al. (1997) suggested that western Tasmania may be the remnant of a continental fragment set adrift by Rodianian break-up, whereas Calvert and Walter (2000) proposed that King Island, along with western Tasmania, rifted away from the Australian craton around ∼ 600 Ma (Fig. 1). Other researchers have developed scenarios in which the island of Tasmania was present as a separate microcontinental block that was positioned outboard of the eastern margin of Gondwana before reattaching at the commencement of the Palaeozoic (Berry et al., 2008).

A popular model that attempts to reconcile the geology observed in Tasmania and adjacent mainland Australia is that of Cayley (2011a). This model proposes that central Victoria and western Tasmania formed a microcontinental block called “VanDieland” that fused with East Gondwana at the end of the Cambrian, possibly terminating the Delamerian Orogeny. VanDieland became entangled in the subduction–accretion system which built the Palaeozoic orogens that now comprise eastern Australia (Fig. 1). Delineating Precambrian continental fragments within southeast Australia has proven difficult, partly due to more recent sedimentary cover that obscures large tracts of the Tasmanides. However, if present, they likely have distinctive structural and seismic velocity characteristics (Glen, 2013).

3 Previous geophysical studies

To date, a variety of geophysical methods have been employed to study the crustal structure of the Tasmanides. Shibutani et al. (1996) applied a non-linear inversion method to RF waveforms to constrain the shear wave velocity beneath broadband seismic stations in eastern Australia. They found that the Moho is relatively shallow (30–36 km depth) and sharp within the cratonic region, and deeper (38–44 km) and transitional along the axis of the Tasmanides. They suggested that crustal thickening of the fold belt by underplating or intrusion of mantle materials may have contributed to this observation. Clitheroe et al. (2000) built on this earlier work by inverting RFs to map broad-scale crustal thickness and the Moho character across the Australian continent. They found that there was generally good agreement between xenolith-derived estimates of the Moho depth and those determined by RF inversion, except beneath the Lachlan Fold Belt, where a broad Moho transition may be present. Overall, however, the RF results were consistent with those determined by Drummond and Collins (1986) and Collins (1991), who used seismic reflection and refraction transects to determine that the Lachlan Fold Belt includes the thickest crust (∼50 km) in eastern Australia. A more recent study by Fontaine et al. (2013a) employed H–κ stacking and non-linear RF inversion to investigate crustal thickness, shear wave velocity structure, as well as dipping and anisotropy of the crustal layers. Their results also indicated a thick crust (∼ 48 km) and an intermediate (2–9 km) crust–mantle transition beneath the Lachlan Fold Belt, which could be attributed to underplating beneath the crust and/or high concentrations of mafic rocks in the mid-lower crust. Their results also showed a dipping Moho and crustal anisotropy in the vicinity of three seismic stations (YNG, CNB, and CAN).

Over the last decade, ambient noise tomography has become a popular tool for studying the structure of the Australian crust. Saygin and Kennett (2010) produced the first group velocity maps of the Australian continent from Rayleigh wave group velocity dispersion in the period range from 5.0 to 12.5 s. Limited spatial resolution (∼ $\mathrm{2}{}^{\circ }\times \mathrm{2}{}^{\circ }$) in our study region means that this model is only able to represent the structure beneath Bass Strait as a broad, low-velocity anomaly. However, the group velocities exhibit a good correlation with known basins and cratons. Subsequent studies using denser arrays covering southeast mainland Australia (Arroucau et al., 2010), southeastern Australia (Young et al., 2013), and northern Tasmania (Young et al., 2011) show good correlations between both group and phase velocity maps and sedimentary and basement terrane boundaries. In order to account for uneven data distribution, Bodin et al. (2012b) used a Bayesian transdimensional inversion scheme to generate group velocity maps that span the Australian continent from multi-scale ambient noise datasets. However, in our study area, their model is of low resolution due to the limited station coverage; hence, few details on crustal structure can be inferred. Bodin et al. (2012a) subsequently applied Bayesian statistics to reconstruct the Moho geometry of Australia using a variety of seismic datasets, which gave an approximate Moho depth of ∼ 30 km beneath Bass Strait. Pilia et al. (2015a, b) and Crowder et al. (2019) derived 3-D shear wave velocity models of the Bass Strait region using ambient noise data from the same array of temporary stations that we exploit in this study. They were able to constrain the lateral and depth extent of the primary sedimentary basins in the region as well as providing insight into the seismic character of the Precambrian microcontinental block that appears to underpin southern Victoria, northwestern Tasmania, and Bass Strait.

Teleseismic tomography has also been used to image the lithosphere beneath southeast Australia, thanks in part to the prolific deployment of short-period seismometers as part of the WOMBAT transportable array project (Rawlinson and Kennett, 2008; Rawlinson et al., 2015, 2016). While the main focus has been on the upper mantle, in Tasmania, where station spacing was denser, some constraints on crustal velocity structure were possible. Rawlinson and Urvoy (2006) found that the crust beneath the East Tasmania Terrane (ETT) was significantly faster than the crust beneath central Tasmania, which may represent a contrast between crust with oceanic provenance in the east and Precambrian continental provenance in the west. Bello et al. (2019b) built on this work by including teleseismic arrival time data from the same temporary deployment as the current study to generate a detailed upper-mantle model of southeast Australia, which revealed that Bass Strait was underlain by lower velocities, consistent with thinned lithosphere as a result of failed rifting during the break-up of Australia and Antarctica.

Active source seismic profiling has also been widely used in southeast Australia to characterize crustal velocity structure (e.g. Finlayson et al., 1980; Collins, 1991; Finlayson et al., 2002; Drummond et al., 2006; Glen, 2013). This has largely focused on the transition from continental to oceanic crust at passive margins, but it has also been used to image major transition zones or faults between orogens (Glen, 2013) or within orogens (Cayley et al., 2011a, b), the latter of which lead to the VanDieland microcontinental model. Rawlinson and Urvoy (2006) jointly inverted teleseismic arrival times and active source wide-angle travel times in northern Tasmania to constrain crustal velocity, Moho geometry, and upper-mantle velocity structure, and they found that both northeastern and northwestern Tasmania are characterized by thinner (< 28 km) and higher-velocity crust compared with central Tasmania.

Potential field data have also been exploited to study the formation and structure of the Tasmanides. Gunn et al. (1997) integrated potential field data (magnetic and gravity), seismic reflection data, outcrop geology, and well information to study the crustal structure of the Australian continent. Their study found that the occurrence of tensional stress, oriented northeast–southwest along basement structures in the Bass Basin, is able to explain the formation of the three major sedimentary basins that overlie dense mafic material, which in turn was formed by mantle decompression processes associated with crustal stretching. From the interpretation of new aeromagnetic data, Morse et al. (2009) delineated the architecture of the Bass Strait basins and their supporting basement structure. Subsequent studies by Moore et al. (2015, 2016) used gravity, magnetic, seismic reflection, and outcrop data to support the hypothesis of a VanDieland microcontinent. Their study showed that VanDieland comprises seven distinct microcontinental ribbon terranes that appear to have amalgamated by the late Cambrian, with major faults and suture zones bonding these ribbon terranes together.

While the last few decades have seen important advances and insights made into our understanding of the southern Tasmanides, there is still limited data on the deep crustal structure beneath Bass Strait, which is our region of interest. Therefore, it is timely that we can exploit, using the RF technique, teleseismic data recorded by a collection of temporary and permanent seismic stations in the region to study the structure of the crust, Moho, and uppermost mantle beneath mainland Australia, Bass Strait, and Tasmania.

4 Data

A collaboration involving five organizations (the University of Tasmania, the Australian National University, Mineral Resources Tasmania, the Geological Survey of Victoria, and FROGTECH) deployed the temporary BASS seismic array from May 2011 to April 2013. It consisted of 24 broadband, three-component seismic stations that spanned northern Tasmania as well as a selection of islands in Bass Strait and southern Victoria. The instruments used were 23 Güralp 40T and 1 Güralp 3ESP sensors coupled to Earth Data PR6-24 data loggers. The permanent stations consist of eight broadband sensors managed by IRIS, GEOSCOPE, and the Australian National Seismic Network (ANSN). The distribution of all 32 seismic stations used in this study is plotted in Fig. 2. Earthquakes with magnitudes m_b > 5.5 at epicentral distances between 30 and 90^∘ comprise the seismic sources used in this analysis (Fig. 3). This resulted in an acceptable azimuthal coverage of earthquakes between the northwest and east of the array, where active convergence of the Australian and Eurasian plate coupled with westward motion of the Pacific plate has produced extensive subduction zones. To the south and southwest of the array, the absence of subduction zones in the required epicentral distance range means that there are significantly fewer events available for analysis from these regions.

Figure 2 Location of seismic stations used in this study superimposed on a topographic and bathymetric map of southeast Australia (modified from Bello et al., 2019a). The boundary of VanDieland is delineated by a thick black dashed line. The thick red dashed line outlines the boundary of the East Tasmania Terrane and Furneaux Islands. The thick white dashed line highlights the eastern sector of the Lachlan Fold Belt. The topography and bathymetry are based on the ETOPO1 dataset (Amante and Eakins, 2009).

5 Methods

5.1 Receiver functions

The RF technique (Langston, 1979) uses earthquakes at teleseismic distances to enable estimation of the Moho depth and shear wave velocity structure in the vicinity of a seismic recorder. If this technique can be applied to a network of stations with good spatial coverage, it represents an effective way of mapping lateral variations in the Moho depth and crustal structure.

A recorded teleseismic wave field at a broadband station can be described by the convolutional model in which operators that represent the source radiation pattern, path effects, crustal structure below the station, and instrument response are combined to describe the recorded waveform. By using deconvolution to remove the effects of the source, path, and response of the instrument (e.g. Langston, 1979), information on local crustal structure beneath the station can be extracted from P–S wave conversions at discontinuities in seismic velocity (Owens et al., 1987; Ammon, 1991).

P-wave RFs were determined from teleseismic P-wave forms using FuncLab software (Eagar and Fouch, 2012; Porritt and Miller, 2018), following preprocessing using the seismic analysis code (SAC) (Goldstein et al., 2003). RFs were computed by applying an iterative time-domain deconvolution scheme developed by Ligorria and Ammon (1999) with a 2.5 s Gaussian filter width. This is achieved by deconvolution of the vertical component waveform from the radial and transverse waveforms with a central frequency of ∼ 1 Hz. This frequency was selected on account of significant source energy detected in the ∼ 1 Hz range of teleseismic P arrivals, which are sensitive to crustal-scale anomalies. It also provides a favourable lateral sensitivity with respect to Fresnel zone width (∼ 15 km at Moho depth) when the conversions from P to S are mapped as velocity and crustal thickness variations.

The complete set of 1765 events (Fig. 3) and 32 stations produced 21 671 preliminary RFs. These RFs were manually inspected using the FuncLab trace editor, and a subset of 9674 RFs were selected for further analysis using the visual clarity of the direct arrivals as an acceptance criterion. Due to high noise levels and fewer events associated with the temporary BASS array dataset, a modest number of good-quality RFs resulted from the above selection method; therefore, different selection criteria were applied that assessed the P-arrival, Moho conversion, and later amplitudes in conjunction with overall noise levels exhibited by the transverse component RFs. This enabled the temporary BASS stations to yield between 2 and 30 good-quality receiver functions, and increased the number of stations where H–κ stacking and NA inversion could be applied from 13 to 20.

Figure 3 Distribution of distant earthquakes (teleseisms) used in this study. The locations of events that are ultimately used for RF analysis are denoted by yellow dots. Concentric circles are plotted at 30^∘ intervals from the centre of Bass Strait. Topography and bathymetry colours are based on the ETOPO1 dataset (Amante and Eakins, 2009).

5.2 H–κ stacking

Having obtained reliable P-wave RFs, the H–κ stacking technique is used to estimate crustal thickness and bulk $V_{\mathrm{p}}/V_{\mathrm{s}}$ for individual stations. We apply the method of Zhu and Kanamori (2000) to stations where the direct Ps (Moho P to S conversion) phase and its multiples are observed. This technique makes use of a grid search to determine the crustal thickness (H) and V_p/V_s (κ) values that correspond to the peak amplitude of the stacked phases. A clear maximum requires a contribution from both the primary phase (Ps) and the associated multiples (PpPs). In the absence of multiples, the maximum becomes smeared out due to the inherent trade-off between crustal thickness (H) and average crustal velocity properties (κ) (Ammon et al., 1990; Zhu and Kanamori, 2000). The H–κ stacking algorithm reduces the aforementioned ambiguity by summing RF amplitudes for Ps and its multiples – PpPs and PpSs+PsPs – at arrival times corresponding to a range of H and $V_{\mathrm{p}}/V_{\mathrm{s}}$ values. In the H–κ domain, the equation for stacking amplitude is

$$\begin{array}{}\text{(1)}& s\left(H,\mathit{\kappa }\right)=\sum _{j=\mathrm{1}}^N{w}_{\mathrm{1}}{r}_j\left(t_{\mathrm{1}}\right)+w_{\mathrm{2}}{r}_j\left(t_{\mathrm{2}}\right)+w_{\mathrm{3}}{r}_j\left(t_{\mathrm{3}}\right),\end{array}$$

where r_j(t_i); i=1, 2, 3 are the RF amplitude values at the expected arrival times t₁, t₂, and t₃ of the Ps, PpPs, and PpSs+PsPs phases respectively for the jth RF; w₁, w₂, and w₃ are weights based on the signal-to-noise ratio ($w_{\mathrm{1}}+w_{\mathrm{2}}+w_{\mathrm{3}}=\mathrm{1}$); and N is the total number of radial RFs for the station. s(H,κ) achieves its maximum value when all three phases stack constructively, thereby producing estimates for H and $V_{\mathrm{p}}/V_{\mathrm{s}}$ beneath the station (see Figs. 5 and S1–S4). In this study, the weighting factors used are w₁=0.6, w₂=0.3, and w₃=0.1. The H–κ approach requires an estimate of the mean crustal P-wave velocity, which is used as an initial value. Based on the results of a previous seismic refraction study (Drummond and Collins, 1986), we use an average crustal velocity of V_p=6.65 km/s to obtain our estimates of H and κ in the study area, noting that H–κ stacking results are much more dependent on $V_{\mathrm{p}}/V_{\mathrm{s}}$ than V_p (Zhu and Kanamori, 2000). To estimate the uncertainties in the H–κ stacking results, we compute the standard deviation of the H and κ values at each station. When only a small number of RFs are available at a station (e.g. four in the case of MILA), the estimates are unlikely to be particularly robust, and in such instances, they are perhaps best viewed as a lower bounds on uncertainty.

While simple to implement, the Zhu and Kanamori (2000) method can suffer from large uncertainties due to its assumption of a simple flat-laying layer over a half-space with constant crustal and upper-mantle properties. Consequently, there are only two search parameters (H and κ) plus a priori information (V_p, weightings), and this method does not account for variation with back-azimuth. These problems can cause non-unique and inaccurate estimates, which can lead to potentially misleading interpretations; for instance, a low-velocity upper-crustal layer can appear as a very shallow Moho in an H–κ stacking search space diagram. Also, a dipping Moho and/or anisotropic layers within the crust can contribute to uncertainty.

5.3 Non-linear waveform inversion

In an effort to refine the crustal model, we invert a stack of the radial RFs by adopting the workflow described by Shibutani et al. (1996). We divide the waveform data (RFs) into four 90^∘ quadrants based on the back-azimuth of their incoming energy. The first quadrant back-azimuth range is from 0 to 90^∘, and an equivalent range in a clockwise direction defines the consecutive quadrants. The second and third quadrants (southeastern and southwestern back-azimuths) have very small numbers of RFs. Data from the first and fourth quadrants are of better quality, with the first quadrant showing more coherency than the fourth quadrant, which is likely due to the orientation of surrounding tectonic plate boundaries; hence, the pattern of P-wave energy radiated towards Australia. Kennett and Furumura (2008) showed that seismic waves arriving in Australia from the northern azimuths undergo multiple scattering but low intrinsic attenuation due to heterogeneity in the lower crust and mantle; this tends to produce prolonged high-frequency coda. An important assumption in our inversion is that we neglect anisotropy and possible Moho dip, which we assume have a second-order influence on the waveforms that we use to constrain 1-D models of the crust and upper mantle.

Visual examination of coherency in P to S conversions allows us to select a subset of RF waveforms for subsequent stacking. This resulted in groups of mutually coherent waveforms after which a move-out correction is then applied to remove the kinematic effect of different earthquake distances prior to stacking using a cross-correlation matrix approach described in Chen et al. (2010) and Tkalčić et al. (2011). Our visual acceptance criteria yields RFs at only 14 out of the 32 stations used for this study. An example of some stacked RFs is given in Fig. 4.

Figure 4 Stacked receiver functions from the Australian National Seismic Network (ANSN) stations TOO, YNG, and MOO as well as the Global Seismographic Network (GSN) station TAU. Small arrows indicate the arrival of the Ps (black), PpPs (red), and PpPs+PsPs (blue) phases from the Moho.

Figure 5 Results from the H–κ stacking analysis for RFs (Zhu and Kanamori, 2000) at MOO, CAN, and TOO stations. (a) Normalized amplitudes of the stack over all back-azimuths along the travel time curves corresponding to the Ps and PpPs phases for each case. (b) The corresponding stacked receiver function for each station.

We invert RFs for 1-D seismic velocity structure beneath selected seismic stations using the neighbourhood algorithm (NA) (Sambridge, 1999a, b) in order to better understand the internal structure of the crust and the nature of the transition to the upper mantle. NA makes use of Voronoi cells to help construct a searchable parameter space, with the aim of preferentially sampling regions of low data misfit. In the inversion process, a Thomson–Haskell matrix method (Thomson, 1950; Haskell, 1953) was used to calculate a synthetic radial RF for a given 1-D (layered) structure. During the inversion, as in Shibutani et al. (1996) and Clitheroe et al. (2000), each model is described by six layers: a layer of sediment; a basement layer; an upper crust, middle crust and lower crust; and an underlying mantle layer, all of which feature velocity gradients and, potentially, velocity jumps across boundaries. The inversion involves constraining 24 parameters: V_s values at the top and bottom of each layer, layer thickness, and the $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratio in each layer (Table 1). The inclusion of $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratio as an unknown primarily aims to accommodate the effects of a sediment layer with limited prior constraints (Bannister et al., 2003). There are two important controlling parameters required by NA: (1) the number of models produced per iteration (n_s); and (2) the number of neighbourhoods re-sampled per iteration (n_r). After a number of trials, we chose the maximum number of iterations to be 5500, with n_s=13 and n_r=13 for all iterations. We employ a chi-square χ² metric (see Sambridge, 1999a, for more details) to compute the misfit function, which is a measure of the inconsistency between the true (${\mathit{\varphi }}_i^{\mathrm{obs}}$) and predicted (${\mathit{\varphi }}_i^{\mathrm{pre}}\left(m\right)$) waveforms for a given model m:

$$\begin{array}{}\text{(2)}& {\mathit{\chi }}_{\mathit{\nu }}^{\mathrm{2}}\left(m\right)={\displaystyle \frac{\mathrm{1}}{\mathit{\nu }}}\sum _{i=\mathrm{1}}^{N_{\mathrm{d}}}\left({\displaystyle \frac{{\mathit{\varphi }}_i^{\mathrm{obs}}-{\mathit{\varphi }}_i^{\mathrm{pre}}\left(m\right)}{{\mathit{\sigma }}_i}}\right),\end{array}$$

where σ_i represents the noise standard deviation determined from ${\mathit{\varphi }}_i^{\mathrm{obs}}$, following the method described in Gouveia and Scales (1998); and ν represents the number of degrees of freedom (the difference between the number of observations and the number of parameters being inverted for). Using the above-mentioned parameters, the inversion targets the 1-D structure that produces the best fit between the predicted and observed RF. Figures 7–9 and S5–S9 present example results of inversions via density plots of the best 1000 data-fitting S-wave velocity models produced by the NA. The optimum data-fitting model is plotted in red.

Table 1 Model parameter bounds used in the neighbourhood algorithm receiver function inversion. $V_{\mathrm{s}}^{\mathrm{upper}}$ and $V_{\mathrm{s}}^{\mathrm{lower}}$ represent the S- wave velocity at the top and bottom of a layer respectively. $V_{\mathrm{p}}/V_{\mathrm{s}}$ represents the P- and S-wave velocity ratio within a layer.

6 Results

6.1 H–κ stacking results

Maps of crustal thicknesses and average $V_{\mathrm{p}}/V_{\mathrm{s}}$ from H–κ stacking in southeast Australia from 16 stations are shown in Fig. 6. At the remaining stations, we could not detect any clear multiples or Moho conversions in the RFs from any direction. A previous study by Chevrot and van der Hilst (2000) has noted that this region is devoid of clear multiples. The crustal thickness for all analysed stations in the study area varies from 23.2 ± 5.0 km (BA02) beneath northwestern Tasmania to 39.1 ± 0.5 km (CAN) beneath the Lachlan Fold Belt, and the variation is strongly correlated with topography. Crust beneath VanDieland (Fig. 6a) is thin in the north (∼ 37.5 km) and south (∼ 33 km), but it appears to be considerably thinner beneath the Victorian and Tasmanian margin Bass Strait (∼ 25 km). The mountainous region of the Lachlan Fold Belt has the deepest Moho at 39.1 ± 0.5 km (CAN) and a corresponding $V_{\mathrm{p}}/V_{\mathrm{s}}$ value of 1.73 ± 0.02. Crust that is consistently between ∼ 31 and 33 km thick lies beneath the East Tasmania Terrane and eastern Bass Strait (ETT + EB). The $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratio varies between ∼ 1.65 beneath station BA11, which also exhibits the thinnest crust, and ∼ 1.93 beneath stations BA19 and BA20 in southern Victoria. There is no obvious correlation between the number of RFs used in the H–κ stacking and the size of the uncertainty in either the Moho depth or $V_{\mathrm{p}}/V_{\mathrm{s}}$, but as mentioned previously, the uncertainty estimates for stations with a low number of RFs are likely to be less robust. Table 2 shows a summary of the H–κ stacking results for the stations that have been analysed.

Figure 6 (a) Variations in the crustal thickness and (b) $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratio taken from the linear (H–κ) stacking results (Table 2). Crustal thickness varies between ∼ 23 and 39 km. The $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratios vary from ∼ 1.65 to 1.93. The thick black dashed line denotes the boundary of VanDieland. The thick magenta dashed line outlines the boundary of the East Tasmania Terrane and eastern Bass Strait (ETT + EB). The thick cyan dashed line highlights the eastern part of the Lachlan Fold Belt. The illuminated topography and bathymetry are based on the ETOPO1 dataset (Amante and Eakins, 2009).

Table 2 Summary of H–κ stacking and NA inversion results for the current study.

6.2 Non-linear inversion results

Results of the NA inversion were successfully obtained for a selection of permanent and temporary stations, as shown in Table 2 and Fig. 10. If the Moho is defined by a gentle velocity gradient, the base of the velocity gradient is used as a proxy for the Moho depth, as done in previous RF (e.g. Clitheroe et al., 2000; Fontaine et al., 2013a) and seismic refraction (Collins, 1991; Collins et al., 2003) studies. We also adopt an upper-mantle velocity of V_p=7.6 km/s (i.e. V_s=4.3–4.4 km/s for $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratios of 1.73–1.77 at the base of the Moho gradient) following Clitheroe et al. (2000), who used this value for RF studies, and Collins et al. (2003), who used V_p > 7.8 km/s for their summary of both seismic refraction and RF results; these V_p values are consistent with global Earth models (e.g. Kennett et al., 1995). Therefore, we also require the S-wave velocity to be > ∼ 4.4 km/s beneath the Moho. We present the S-wave velocity profiles from the NA inversion for stations CAN, MOO, TOO, YNG, BA13, and BA17 in Figs. 7–9 as well as the observed and predicted RFs. The S-wave velocity inversion results for the remaining stations are included in the Supplement (Figs. S5–S8). In assigning the Moho depth, we consider three criteria to examine the quality of the inversion result: (1) misfit value χ²; (2) the quality of the RF stack (which is based on our ability to pick the direct and multiple phases); and (3) the visual fit between the synthetic and observed RF. Models that fail to fit significant arrivals in the observed RF are rejected. Based on these criteria, the inversion results are classified as follows:

• very good – very low χ² (typically < 0.4), very good visual fit to direct and multiple phases;

• good – low χ² (typically 0.4–0.8), direct phases clearly visible, multiple phases less clear, and a good visual fit to all major identifiable phases;

• poor – medium to high χ² (in the range 0.8–1.2), direct phases visible, multiple phases unclear, and moderate visual fit to some identifiable phases.

Figure 7 (a, c) Seismic velocity models for the CAN and MOO stations obtained from the neighbourhood algorithm (Sambridge, 1999a). The grey area indicates all the models searched by the algorithm. The best 1000 models are indicated by the yellow to green colours; the best model (smallest misfit) is shown using the red line, for both the S-wave velocity (top x axis) and $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratio (bottom x axis), and the white line is the average velocity model. Small black arrows denote the estimated depth of the Moho. (b, d) Waveform matches between the observed stacked receiver functions (black) and predictions (blue) based on the best models. “Misfit” refers to the chi-square estimate as defined by Eq. (2).

In general, the optimum χ² value is normally considered to be 1, as below this value, the tendency is to fit noise rather than signal. However, this is for the ideal case when the number of degrees of freedom and the absolute values of the data uncertainty are well known (e.g. in the case of a synthetic test). In the case of observational data, these values are often poorly constrained, so using the relative χ² values coupled with visual assessment of the data fit appears to be reasonable. With regard to the character of the crust–mantle transition, this study classifies the transition zone as sharp ≤ 2 km, intermediate 2–10 km, or broad ≥ 10 km, as initially proposed by Shibutani et al. (1996) and modified by Clitheroe et al. (2000).

We note that for the seven permanent stations for which we produce receiver function inversion and H–κ stacking results, five have estimates of the Moho depth from previous receiver function studies. Clitheroe et al. (2000) estimated the Moho depth at 49 km beneath CAN based on a non-linear inversion, which is ∼ 10 km greater than the results we obtain for both NA inversion and H–κ stacking (see Sect. 7.1 for further discussion of this discrepancy). Ford et al. (2010) determined the Moho depth beneath the MOO, TOO, TAU, and YNG stations using H–κ stacking and found values (compared with our H–κ stacking results) of 33 ± 3 km (33.0 ± 1.2 km), 34 ± 3 km (37.5 ± 1.2 km), 32 ± 3 km (33.5 ± 1.9 km), and 33 ± 2 km (37.3 ± 0.5 km) respectively. These are all within error, with the slight exception of YNG station, located in Young, on the western flanks of the Great Dividing Range, where we might expect the crust to be slightly thicker than average. Overall, however, these similarities suggest that our results are likely to be robust.

Figure 8 (a, c) Seismic velocity models for TOO and YNG stations obtained from the neighbourhood algorithm. (b, d) Comparison between the observed stacked and the predicted receiver functions from the NA inversion. See the caption of Fig. 7 for more details.

Figure 9 Seismic velocity models for temporary stations BA13 (a) and BA17 (b) obtained from the neighbourhood algorithm. See Figs. S6–S9 in the Supplement for all receiver function inversion results for the temporary BASS network, including waveform fits (Supplement Fig. S7 includes the waveform fit for stations BA13 and BA17).

7 Discussion

For convenience, the seismic stations are separated into three groups (Fig. 2) based on tectonic setting and the results obtained. Stations YNG, CAN, CNB, MILA, and BA13 are located in the Lachlan Fold Belt; stations BA02, BA11, BA19, BA20, TAU, MOO, and TOO sit above the VanDieland microcontinental block; and stations BA07, BA08, BA09, and BA17 lie in the East Tasmania Terrane and eastern Bass Strait (ETT + EB). Stations BA22 and BA24 lie to the west of VanDieland. This discussion focuses on crustal thickness, the nature of the Moho, and crustal velocity and velocity ratio variations from H–κ stacking and the 1-D S-wave velocity models. Overall, the agreement between Moho depths obtained from the H–κ stacking results and NA inversion is generally within error (Table 2), which makes a joint interpretation more straightforward. Comparison is also made to other studies that have examined crustal seismic properties in southeast Australia, and we attempt to integrate our new findings with previous results from teleseismic tomography, SKS splitting, and ambient noise tomography in order to better understand the crust and upper-mantle structure and dynamics beneath this region.

7.1 Lateral variation in crustal thickness and the nature of the Moho

The RF analysis clearly reveals the presence of lateral changes in crustal thickness that span mainland Australia through Bass Strait to Tasmania (Figs. 6 and 10; in the latter case, RF depths from previous studies are also included for reference). The stations located in the Palaeozoic Lachlan Fold Belt reveal a generally thick crust that ranges between ∼ 37 and 40 km. Although the Moho was picked as a velocity jump for stations YNG, CAN, and CNB, the velocity nonetheless tends to continue to increase with depth below the discontinuity. This, coupled with the fact that Clitheroe et al. (2000) estimate the Moho to be almost 10 km deeper beneath CAN, is consistent with the presence of mafic underplating (e.g. Drummond and Collins, 1986; Shibutani et al., 1996; Clitheroe et al., 2000), sourced from the ambient convecting mantle. The top and bottom of such a layer could feature a velocity step with depth, and its internal structure is likely to be layered and/or gradational, resulting in uncertainty in the true Moho depth. Based on deep crustal reflection profiling, Glen et al. (2002) suggested that the deep Moho underlying the Lachlan Orogen results from magmatic underplating that added a thick Ordovician mafic layer at the base of the crust coupled with a thick sequence of Ordovician mafic rocks that can be found in the mid and lower crust. Finlayson et al. (2002) and Glen et al. (2002) also inferred the presence of underplating near CNB and CAN from seismic refraction data. Collins (2002) postulated that the underplating might have occurred in the back-arc region of a subduction zone due to pronounced adiabatic decompression melting in the asthenosphere. The seismic tomography model of Rawlinson et al. (2010, 2011) exhibits an increase in P-wave speed at 50 km depth beneath CAN, CNB, and YNG, and the authors suggest that magmatic underplating may be the cause of the high-velocity anomaly. A recent study by Davies et al. (2015) identified the longest continental hotspot track in the world (over 2000 km total length), which began in north Queensland at ∼ 33 Ma, and propagated southward underneath the present-day Lachlan Fold Belt and Bass Strait. The magmatic underplating could therefore be a consequence of the passage of the continent above a mantle upwelling leading to a more diffuse crust–mantle transition zone. The thickened crust and a transitional Moho observed in the Lachlan Fold Belt are consistent with the proposed delamination models of Collins and Vernon (1994).

Figure 10 Map showing crustal thickness variations based on the S-wave velocity inversion results of this study (circles) and previous studies (triangles) (Clitheroe et al., 2000; Fontaine et al., 2013a, b; Shibutani, 1996; Tkalcic et al., 2012). The topography and bathymetry are based on the ETOPO1 dataset (Amante and Eakins, 2009).

Strong lateral changes in crustal seismic structure (Figs. 6, 10) beneath VanDieland appear to be a reflection of the region's complex tectonic history. The thick crust (∼ 37 km) beneath the Selwyn Block (see Fig. 1 for its location) – within the northern margin of VanDieland in southern Victoria – thins dramatically to ∼ 26 km as it enters Bass Strait, increases to ∼ 30 km beneath King Island (BA11), and then thins to ∼ 23 km beneath northwestern Tasmania, before increasing to ∼ 33 km in southern Tasmania. The results in southern Tasmania agree with those of Korsch et al. (2002) from a seismic reflection profile adjacent to the TAU and MOO seismic stations. The thinner crust beneath Bass Strait and its margins may be a consequence of lithospheric thinning and/or delamination associated with failed rifting that accompanied the break-up of Australia and Antarctica (Gaina et al., 1998). Stations BA07, BA08, BA09, and BA17 (ETT + EB) collectively indicate crust of relative uniform thickness (∼ 31–32 km; Fig. 10a, b). Relative to western Bass Strait, the crust is slightly thicker in this part of the study area, which may suggest underplating associated with a Palaeozoic subduction system (e.g. Drummond and Collins, 1986; Gray and Foster, 2004).

In general, our understanding of crustal thickness variations are limited by station separation, so it is difficult to determine whether smooth variations in thickness or step-like transitions explain the observations.

7.2 $V_{\mathrm{p}}/V_{\mathrm{s}}$ and bulk crustal composition

$V_{\mathrm{p}}/V_{\mathrm{s}}$ can constrain chemical composition and mineralogy more robustly than P- or S-wave velocity in isolation (Christensen and Fountain, 1975). We observe variations in $V_{\mathrm{p}}/V_{\mathrm{s}}$ across the study region, which we can largely equate with variations in composition or melt. Studies in mineral physics and field observations show (1) an increase in $V_{\mathrm{p}}/V_{\mathrm{s}}$ with decreasing SiO₂ content in the continental crust (Christensen, 1996), and (2) partial melt is revealed by elevated $V_{\mathrm{p}}/V_{\mathrm{s}}$, especially if the anomaly is localized to an intra-crustal layer (Owens and Zandt, 1997). A more felsic (SiO₂) composition in the lower crust is represented by a lower $V_{\mathrm{p}}/V_{\mathrm{s}}$, which reflects removal of an intermediate-mafic zone by delamination, whereas a more mafic lower crust is revealed by higher $V_{\mathrm{p}}/V_{\mathrm{s}}$ (> 1.75) which may be due to underplated material (Pan and Niu, 2011). However, lower crustal delamination can also result in decompression melting, which can yield elevated $V_{\mathrm{p}}/V_{\mathrm{s}}$ (He et al., 2015). We interpret the variation of observed $V_{\mathrm{p}}/V_{\mathrm{s}}$ in the southern Tasmanides to be a consequence of compositionally heterogeneous crust and localized partial melt that may likely be sourced from recent intraplate volcanism (Rawlinson et al., 2017).

Figure 6b shows the distribution of bulk $V_{\mathrm{p}}/V_{\mathrm{s}}$ across the study area. The pattern of $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratios appears to delineate three distinct zones of crust. Beneath the Lachlan Orogen, values are ∼ 1.75, which is consistent with the presence of a mafic lower crust, as suggested by a number of other studies (Drummond and Collins, 1986; Shibutani et al., 1996; Clitheroe et al., 2000; Finlayson et al., 2002). Beneath eastern Bass Strait, the $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratios are slightly lower, with BA07, BA08, and BA09 exhibiting values of 1.70, 1.70, and 1.71 respectively. These values are in agreement with constraints from seismic reflection and refraction studies (Finlayson et al., 2002; Collins et al., 2003) and may indicate a felsic to intermediate crustal composition. The geology of Flinders Island, which hosts both BA07 and BA08, is dominated by Devonian granites, which is consistent with this observation. Beneath VanDieland, $V_{\mathrm{p}}/V_{\mathrm{s}}$ is highly variable, with the greatest contrast between BA11 (∼ 1.65) and BA19/20 (∼ 1.93), and BA19/20 and TOO (1.68). BA11 is located on King Island, which is characterized by Precambrian and Devonian granite outcrops, which may help explain the low $V_{\mathrm{p}}/V_{\mathrm{s}}$. The high $V_{\mathrm{p}}/V_{\mathrm{s}}$ beneath BA19/20 is harder to explain but could be caused by melt in the crust associated with the Newer Volcanics Province, which sits along the Cosgrove intraplate volcanic track and last erupted only ∼ 4.6 ka (Rawlinson et al., 2017). The return to lower $V_{\mathrm{p}}/V_{\mathrm{s}}$ beneath TOO over a relatively short distance (∼ 100 km) is also difficult to explain, but we note that this region of Victoria is underlain by granite intrusions.

In summary, the crust beneath VanDieland exhibits the greatest lateral heterogeneity in $V_{\mathrm{p}}/V_{\mathrm{s}}$, which likely reflects considerable variations in composition and the presence of melt. This can partially be explained by the tectonic history of the region, which includes failed rifting in Bass Strait accompanied by widespread magma intrusion and granite emplacement as well as, more recently, the passage of a plume (Rawlinson et al., 2017). Furthermore, Moore et al. (2015) used reflection transects and potential field data to infer that VanDieland is comprised of up to seven continental ribbon terranes that are bounded by major faults and suture zones, which were likely amalgamated by the end of the Proterozoic. Hence, considerable variations in composition and, in turn, $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratio are to be expected.

7.3 Moho depth comparison

Prior to this study, a variety of seismic methods have been used to constrain the Moho depth in southeast Australia, including receiver functions, reflection profiling, and wide-angle refection and refraction experiments. In an effort to combine the results from all of these studies into a single synthesis, Kennett et al. (2011) developed the AusMoho model; this included Moho depth estimates from over 11 000 km of reflection transects across the continent, numerous refraction studies, and 150 portable and temporary stations. Due to irregular sampling, the detail of this model is highly variable; for example, the region beneath Bass Strait is constrained by only five measurements, whereas the central Lachlan Fold Belt around Canberra (see Fig. 1 for location) features relatively dense sampling at ∼ 50 km intervals or less.

AusMoho includes previous receiver function results from Shibutani et al. (1996), Clitheroe et al. (2000), Fontaine et al. (2013a), and Tkalcic et al. (2012), as well as reflection and refraction transects in Tasmania, parts of the Lachlan Orogen, and western Victoria. Figure 11 illustrates AusMoho for our study region, which exhibits large variations in the Moho depth (from ∼ 10 to > 50 km). These extremes are due to the presence of oceanic crust outboard of the passive margin of the Australian continent, and the root beneath the Southern Highlands, which represent the southern extension of the Great Dividing Range in New South Wales. Superimposed on Fig. 11 are Moho depths from the four previous receiver function studies cited above as well as NA inversion and H–κ depth estimates from this study. As expected, the correlation between the previous RF results and AusMoho is generally good, as they were part of the dataset used to build this model. In places where they do not match, this can be attributed to the presence of seismic refraction or reflection lines that were also used to constrain AusMoho.

Figure 11 Comparison between the AusMoho model (background colour map) and the Moho depths determined through RF analysis in this and previous studies. Small coloured circles denote the Moho depths determined from H–κ stacking, whereas large coloured circles correspond to receiver function estimates. When both H–κ- and NA-derived depths are available at a single station, the smaller H–κ circle is superimposed on the larger NA circle so that both depths can be observed on the one plot. Moho depths determined from previous RF studies are denoted by triangles.

In general, the agreement between the results from this study and AusMoho is good, but there are exceptions. For instance, CAN, CNB, YNG, and MILA tend to be somewhat shallower than AusMoho. However, this can be attributed to the likely presence of mafic underplating alluded to earlier, which can effectively yield two options for the Moho transition due to an expected high $V_{\mathrm{p}}/V_{\mathrm{s}}$ (> 1.85) in the underplate layer (e.g. Cornwell et al., 2010). AusMoho Moho depths beneath BA07 and BA08 are considerably shallower than our estimates, which we attribute to a lack of data coverage in this region. Sizable discrepancies also exist beneath BA02, BA19, and BA20: in the former case, the uncertainty in our H–κ stacking estimate is 5 km, which may be a factor here; in the latter case, we also note that there is sparse data coverage southeast of Melbourne to constrain AusMoho, so it would appear that our new Moho depths are more likely to be correct. Overall, while there is good consistency between AusMoho and our new results, any updated version of AusMoho should incorporate the Moho depth estimates from this study.

Figure 12 Composite result of teleseismic tomography (mantle velocity anomalies), ambient noise (crustal velocity anomalies), receiver functions (Moho), and shear wave splitting (inferred mantle flow relative to over-riding plate). Velocity slices are taken at 148^∘ E. Note that the crustal model produced from ambient noise tomography is defined in terms of V_sv, but it was converted to V_p in the study of Bello et al. (2019b) to permit its inclusion in the starting model for the inversion of teleseismic P-wave arrival time residuals. In this figure, the crustal V_p anomalies are shown.

Although AusMoho did make use of results from a 3-D wide-angle reflection and refraction survey of Tasmania (offshore shots and onshore stations), it only used a few sample points for the final Moho model (Kennett et al., 2011); therefore, the resolution of AusMoho is considerably less than the Moho model produced by Rawlinson et al. (2001). Consequently, we plot our three RF results on top of this model in the Supplement (Fig. S10). The agreement between the Moho model and RF depths beneath MOO and TAU is good, but RF estimates beneath BA02 are shallower than the Moho model by about 4 km. However, this is within the margin of error for the H–κ stacking result.

7.4 Synthesis

In this final section, we present a synthesis of results for southeast Australia that are based on (1) our new receiver function results, (2) teleseismic SKS splitting results from Bello et al. (2019a), (3) teleseismic tomography undertaken by Bello et al. (2019b), (4) ambient noise crustal imaging results from Young et al. (2013), and (5) AusMoho (Kennett et al., 2011). This synthesis is encapsulated in the plot shown in Fig. 12, which is a representative transect through the Lachlan Orogen south through Bass Strait and into Tasmania. Moho depths are taken from AusMoho and are refined where additional information is available from our new RF results; crustal P-wave velocity is taken from the ambient noise results (following conversion from S-wave velocity; see Bello et al., 2019b, for more details); and mantle P-wave velocities are taken from Bello et al. (2019b). Arrows are based on interpreted mantle flow patterns undertaken as part of the shear wave splitting study. This previous study used approximately the same temporary and broadband station network as that used in the current study and found that fast axis orientations of anisotropy were aligned with contemporary plate motion (north-northeast) beneath the Lachlan Orogen, but beneath Bass Strait, a radial pattern was observed that is consistent with an upwelling mantle that impinges on the lithosphere and spreads out in all directions. Interestingly, the location of this phenomenon corresponds approximately to the predicted location of the Cosgrove hotspot track source (Davies et al., 2015) and may be caused by an upwelling mantle plume. Thus, the low velocities in the upper mantle beneath Bass Strait may be due to elevated temperatures and melt, although it is not straightforward to explain the higher velocities below 200 km depth in this context.

The thicker Moho boundary beneath the Lachlan Orogen (Fig. 12) reflects the likely presence of underplating, which makes the base of the crust harder to discern seismically. However, the crust is clearly thicker here than beneath Bass Strait or Tasmania. The Moho depth beneath the northern part of the Fig. 12 is not constrained by our RF results, but according to AusMoho, it is relatively flat, which is consistent with Precambrian crust, and there is a faster mantle lithosphere. The strong variations in crustal velocity beneath Bass Strait can be attributed to failed rifting resulting in the formation of thick (> 10 km) sedimentary basins and elevated temperatures (lower velocities) as well as the intrusion of mafic-rich materials into the lower and mid crust (higher velocities).

8 Conclusions

We used H–κ stacking of teleseismic RFs to determine the crustal thickness and $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratio and generate 1-D S-wave velocity profiles of the crust from RF inversion in order to investigate the internal crustal velocity structure beneath the southern Tasmanides in southeast Australia. Our main findings are summarized below.

• The thick crust and broad crust–mantle transition beneath the Lachlan Fold Belt may be caused by magmatic underplating of mafic materials beneath the crust, which is consistent with an elevated $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratio (relative to ak135) of ∼ 1.73. Thicker crust is also to be expected from the elevated topography of the eastern Lachlan Fold Belt.

• The crustal structure is complex beneath VanDieland. It thins considerably from the northern tip of the microcontinent (∼ 37 km) into Bass Strait (∼ 26 km) and northern Tasmania (∼ 23 km), yet the crust is somewhat thicker in southern Tasmania (∼ 33 km) compared with Bass Strait. This may in part be due to the complex origins of the microcontinent, which appears to be comprised of multiple Precambrian continental ribbons, but is also likely due to failed rifting in Bass Strait before and during the separation of Australia and Antarctica. This resulted in lithospheric stretching and delamination, magmatic intrusion, and the deposition of thick sedimentary sequences. Recent intraplate volcanism and the possible progression of a mantle plume beneath VanDieland in the last few thousand years may also have produced compositional heterogeneity and melt in the crust. Such events are likely to contribute significantly to variations in crustal thickness and the pronounced changes in $V_{\mathrm{p}}/V_{\mathrm{s}}$ that we observe.

• Stations within the ETT + EB collectively indicate crust of uniform thickness (∼ 31–32 km) and uniform $V_{\mathrm{p}}/V_{\mathrm{s}}$ (∼ 1.70), which clearly distinguishes it from VanDieland. This region of the crust likely represents a southern continuation of the Lachlan Orogen and is, therefore, underpinned by crust of oceanic origin.

• Comparison of our new Moho depth results with the AusMoho model reveals an overall consistency, although at some of our station locations where AusMoho has few constraints, there are noticeable differences, such as southern Victoria and beneath Flinders Island. The discrepancies beneath the Lachlan Orogen are attributed to the presence of underplated mafic material, which can obfuscate the location of the Moho.

• A synthesis of our new RF results with pre-existing teleseismic tomography, shear wave splitting, and ambient noise studies reveals a complex lithosphere that has clearly been impacted by orogeny (thickened crust), failed rifting beneath Bass Strait (thinned crust and complex crustal velocities), and recent intraplate volcanism (high $V_{\mathrm{p}}/V_{\mathrm{s}}$ ratios and a radial pattern of fast anisotropy patterns above a presumed zone of mantle upwelling).