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 byH–κ stacking and receiver function inversion, which point to the following: 1. a ∼ 39.0 km thick crust; an intermediate–high Vp/Vs 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 Vp/Vs 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.


Introduction 28
The Phanerozoic Tasmanides (Collins and Vernon, 1994;Coney, 1995;Coney et al., 1990)  The goal of this paper is to provide fresh insight into the crust and Moho structure beneath the southern 54 Tasmanides using P-wave RFs, explain the origin of the lateral heterogeneities that are observed and explore the 55 geological relationship between the different tectonic units that constitute the southern Tasmanides, thereby 56 facilitating a better grasp of the region's tectonic history. 57 southern Tasmanides, there still remains limited data on the deep crustal structure beneath Bass Strait, which is 151 our region of interest. It is therefore timely that can exploit, using the RF technique, teleseismic data recorded 152 by a collection of temporary and permanent seismic stations in the region to study the structure of the crust, 153 Moho and uppermost mantle beneath mainland Australia, Bass Strait and Tasmania. Earthquakes with magnitudes ! > 5.5 at epicentral distances between 30 • and 90 • comprise the seismic 163 sources used in this analysis (Fig. 3). This resulted in an acceptable azimuthal coverage of earthquakes between 164 the northwest and east of the array, where active convergence of the Australian and Eurasian plate coupled with 165 westward motion of the Pacific plate has produced extensive subduction zones. To the south and southwest of 166 the array, the absence of subduction zones in the required epicentral distance range means that there are 167 significantly fewer events available for analysis from these regions. 168

Receiver functions 170
The RF technique (Langston, 1979) uses earthquakes at teleseismic distances to enable estimation of Moho 171 depth and shear wave velocity structure in the neighbourhood of a seismic recorder. If this technique can be 172 applied to a network of stations with good spatial coverage, it represents an effective way of mapping lateral 173 variations in Moho depth and crustal structure. The coverage and quality of broadband data available for this 174 study provides a sound basis on which to examine the crustal structure of the southern Tasmanides. 175 A recorded teleseismic wavefield at a broadband station can be described by the convolutional model in which 176 operators that represent the source radiation pattern, path effects, crustal structure below the station and 177 instrument response are combined to describe the recorded waveform. By using deconvolution to remove the 178 effects of the source, path and response of the instrument (e.g. Langston, 1979), information on local crustal 179 structure beneath the station can be extracted from P -to S -wave conversions at discontinuities in seismic 180 velocity (Owens et al., 1987;Ammon, 1991). 181 P-wave RFs were determined from teleseismic P-waveforms using FuncLab software (Eagar and Fouch, 2012; 182 Porritt and Miller, 2018), following preprocessing using the seismic analysis code (SAC) (Goldstein et al., 183 2003). The complete set of 1765 events (Fig. 3) and 32 stations produced 21,671 preliminary RFs. These RFs 184 were manually picked using the FuncLab trace editor, and by using the clarity of the direct arrivals as an 185 acceptance criteria, a total of 9,674 RFs were retained for further analysis. The RFs were computed using an 186 https://doi.org/10.5194/se-2020-74 Preprint. Discussion started: 8 July 2020 c Author(s) 2020. CC BY 4.0 License.
iterative time-domain deconvolution scheme developed by Ligorria and Ammon, 1999 with a 2.5 s Gaussian 187 filter width. This is performed by deconvolution of the vertical component waveform from the radial and 188 transverse waveforms with a central frequency of 1 Hz. This frequency was selected on account of significant 189 source energy detected in the 1 Hz range of teleseismic P arrivals, which are sensitive to crustal-scale 190 anomalies. It also provides a favourable lateral sensitivity with respect to Fresnel zone width (∼15 km at Moho 191 depth) when the conversions from P to S are mapped as velocity and crustal thickness variations. 192

H-κ stacking 193
Having obtained reliable P-wave RFs, the H-κ stacking technique is used to estimate crustal thickness, 194 Poisson's ratio and bulk V p /V s . We apply the method of Zhu and Kanamori (2000) to stations where the direct 195 P s (Moho P-to-S conversion) and its multiples are observed. This technique makes use of a grid search to 196 determine the crustal thickness (H) and V p /V s (κ) values that correspond to the peak amplitude of the stacked 197 phases. A clear maximum requires a contribution from both the primary phase (P s ) and the associated multiples 198 (P p P s and P p S s +P s P s ). In the absence of multiples, the maximum becomes smeared out due to the inherent 199 trade-off between crustal thickness (H) and average crustal velocity properties (κ) (Ammon et al., 1990;Zhu and 200 Kanamori, 2000). The H-κ stacking algorithm reduces the aforementioned ambiguity by summing RF 201 amplitudes for P s and its multiples P p P s and P p S s +P s P s at arrival times corresponding to a range of H and 202 V p /V s values. In the H-κ domain the equation for stacking amplitude 203 where ! ! ; = 1,2,3 are the RF amplitude values at the expected arrival times ! , ! , and ! of the Ps, PpPs, 204 PpSs +PsPs phases respectively for the !! RF, ! , ! , ! are weights based on the signal to noise ratio 205 ! + ! + ! = 1 , and N is the total number of radial RFs for the station. s(H,κ) achieves its maximum value 206 when all three phases stack constructively, thereby producing estimates for H and Vp/Vs beneath the station. In 207 this study, the weighting factors used are ! = 0.6, ! = 0.3 and ! = 0.1 (Zhu and Kanamori, 2000). The H-κ 208 approach requires an estimate of the mean crustal P-wave velocity, which is used as an initial value. Based on 209 the results of a previous seismic refraction study (Drummond and Collins, 1986), we use an average crustal 210 velocity of Vp = 6.65 km/s to obtain our estimates of H and κ in the study area, noting that H-κ stacking results 211 are much more dependent on Vp/Vs than Vp (Zhu and Kanamori, 2000). To estimate the uncertainties in the H-κ 212 stacking results, we compute the standard deviation of the H and κ values at each station. 213 H-κ stacking can also be used to determine Poisson's ratio, which is a useful parameter for inferring the physical 214 and compositional properties of the crust (Christensen, 1996) and providing insight into fractures, fluids, and 215 partial melt (e.g. Mavko, 1980). The Poisson's ratio σ can be determined from κ using the equation 216 where κ = V p /V s . While simple to implement, the Zhu and Kanamori (2000) method can suffer from large 217 https://doi.org/10.5194/se-2020-74 Preprint. Discussion started: 8 July 2020 c Author(s) 2020. CC BY 4.0 License. uncertainties due to its assumption of a simple flat-laying layer over half-space with constant crustal and upper 218 mantle properties. Consequently, there are only two search parameters (H and κ) plus a priori information (V p , 219 weightings) and it does not account for variation with backazimuth. These problems can cause non-unique and 220 inaccurate estimates, which can lead to potentially misleading interpretations; for instance a low velocity upper 221 crustal layer can appear as a very shallow Moho in an H-κ stacking search space diagram. Also, a dipping Moho 222 and/or anisotropic layers within the crust can contribute to uncertainty. 223

Nonlinear waveform inversion 224
In an effort to refine the crustal model, we invert a stack of the radial RFs by adopting the workflow described 225 by Shibutani et al. (1996). We divide the waveform data (RFs) into four 90 • quadrants based on the backazimuth 226 of their incoming energy. The 1 st quadrant backazimuth range is from 0 • and 90 • , and an equivalent range in a 227 clockwise direction defines the consecutive quadrants. The 2 nd and 3 rd quadrants (south-eastern and south-228 western backazimuths) have very small numbers of RFs. Data from the 1 st and 4 th quadrants are of better quality, 229 with the 1 st quadrant showing more coherency than the 4 th quadrant, which is likely due to the orientation of 230 surrounding tectonic plate boundaries and hence the pattern of P-wave energy radiated towards Australia. 231 Kennett and Furumura (2008) showed that seismic waves arriving in Australia from the northern azimuths 232 undergo multiple scattering but low intrinsic attenuation due to heterogeneity in the lower crust and mantle; this 233 tends to produce prolonged high-frequency coda. An important assumption in our inversion is that we neglect 234 anisotropy and possible Moho dip, which we assume have a second order influence on the waveforms we use to 235 constrain 1-D models of the crust and upper mantle. 236 Visual examination of coherency in P to S conversions allows us to select a subset of RF waveforms for 237 subsequent stacking. This resulted in groups of mutually coherent waveforms after which a moveout correction 238 is then applied to remove the kinematic effect of different earthquake distances prior to stacking using a cross-239 correlation matrix approach described in Chen et al. (2010) and . Our strict criteria give 240 reliable RFs at only 6 out of the 32 stations used for this study. An example of some stacked RFs is given in Fig.  241 4. 242

Neighborhood algorithm 243
We invert RFs for 1-D seismic velocity structure beneath selected seismic stations using the Neighbourhood 244 Algorithm or NA (Sambridge, 1999a,b) in order to better understand the internal structure of the crust and the 245 nature of the transition to the upper mantle. NA makes use of Voronoi cells to help construct a searchable 246 parameter space, with the aim of preferentially sampling regions of low data misfit. In the inversion process, a 247 Thomson-Haskell matrix method (Thomson, 1950 andHaskell, 1953) was used to calculate a synthetic radial 248 RF for a given 1-D (layered) structure. During the inversion, as in Shibutani et al. (1996)  iterations. We employ a chi-squared ! metric to compute the misfit function, which is a measure of the 258 inconsistency between the true ∅ ! !"# , and predicted, ∅ ! !"# (m) waveforms for a given model : 259 where ! represents the noise standard deviation determined from ∅ ! !"# , as explained by Gouveia and Scales 260 (1998), and represents the number of degrees of freedom. Using the above stated parameters, the inversion 261 targets the 1-D structure that produces the best fit between the predicted and observed RF.   (Fig. 6a). The crustal thickness is more or less uniform beneath the Lachlan 282 Fold Belt, East Tasmania Terrane and eastern Bass Strait. 283 The majority of our study region has a low-to-intermediate Poisson's ratio. Poisson's ratio is highest (0.262 ± 284 0.014) in the Lachlan Fold Belt (see Table 2 Strong lateral changes in crustal structure and/or composition beneath VanDieland appear to be a reflection of 351 the region's complex tectonic history ( Fig. 6 and 9). The thick crust (37.5 ± 1.2 km) beneath the Selwyn Block -352 within the northern margin of VanDieland in southern Victoria -thins (to 30.5 ± 2.1 km) as it enters Bass Strait, 353 yet in southern Tasmania, at stations TAU and MOO, the crust is thicker (33.5 ± 1.9 km). This is reflected in 354 both the NA inversion and H-κ stacking depth estimates where a sharp Moho is observed beneath this region of 355 the study area ( Fig. 6 and 9) In general, our understanding of crustal thicknesses variations are limited by station separation, so it is difficult 370 to determine whether smooth variations in thickness or step-like transitions explain the observations. 371

Poisson's ratio, V p /V s and average crustal composition 372
Poisson's ratio, which shares an inverse squared relationship to Vp/Vs (Eq. 2) can constrain chemical 373 composition and mineralogy more robustly than P-or S-wave velocity in isolation (Christensen and Fountain, 374 1975

Acknowledgments 466
The work in this paper was performed as part of a PhD study and has been jointly funded by Abubakar Tafawa 467 Balewa University (ATBU), Bauchi, Nigeria and the University of Aberdeen, UK. The authors acknowledge the 468 efforts of staff, students and fieldwork technicians from the Australian National University and University of 469 Tasmania, who deployed the temporary BASS array used in this study. We also thank Qi Li and Armando 470 Arcidiaco for their efforts in BASS data pre-processing and archiving. Australian Research Council Grant 471 LP110100256 supported the BASS deployment. We are grateful to IRIS and Geoscience Australia for providing 472 data from several stations in mainland Australia and Tasmania. Figure 1 was made using Inkscape software 473 (Harrington, et. al., 2005)