Ambient-noise tomography of the wider Vienna Basin region

We present a new 3-D shear-velocity model for the top 30 km of the crust in the wider Vienna Basin region based on surface waves extracted from ambient-noise cross-correlations. We use continuous seismic records of 63 broad-band stations of the AlpArray project to retrieve interstation Green’s functions from ambient-noise cross-correlations in the period range from 5 to 25 s. From these Green’s functions, we measure Rayleigh group traveltimes, utilizing all four components of the cross-correlation tensor, which are associated with Rayleigh waves (ZZ, RR, RZ and ZR), to exploit multiple measurements per station pair. A set of selection criteria is applied to ensure that we use high-quality recordings of fundamental Rayleigh modes. We regionalize the interstation group velocities in a 5 km × 5 km grid with an average path density of ∼ 20 paths per cell. From the resulting group-velocity maps, we extract local 1-D dispersion curves for each cell and invert all cells independently to retrieve the crustal shear-velocity structure of the study area. The resulting model provides a previously unachieved lateral resolution of seismic velocities in the region of ∼ 15 km. As major features, we image the Vienna Basin and Little Hungarian Plain as low-velocity anomalies, and the Bohemian Massif with high velocities. The edges of these features are marked with prominent velocity contrasts correlated with faults, such as the Alpine Front and Vienna Basin transfer fault system. The observed structures correlate well with surface geology, gravitational anomalies and the few known crystalline basement depths from boreholes. For depths larger than those reached by boreholes, the new model allows new insight into the complex structure of the Vienna Basin and surrounding areas, including deep low-velocity zones, which we image with previously unachieved detail. This model may be used in the future to interpret the deeper structures and tectonic evolution of the wider Vienna Basin region, evaluate natural resources, propagation and improve earthquake locations, among others.


I N T RO D U C T I O N
Accurate seismic velocity models improve our understanding of structures and processes in the solid Earth. At regional scale, such models are useful for seismic hazard assessment, better location of regional and local seismic events, understanding the tectonic evolution of a region, and improving the evaluation of natural resources, among others. They provide insight into greater depth regions, which are not well-sampled by other geophysical methods, and not accessible by near-surface geology. In this study, we present a new seismic velocity model of the wider Vienna Basin region to provide new insight into its complex structure. Because of population density and sensitive infrastructure, understanding the regional * www.alparray.ethz.ch processes and structures in and around the Vienna Basin-one of the seismically most active regions in Austria-is of critical importance.
The Vienna Basin is a thin-skinned pull-apart basin in the Alpine-Carpathian transition zone that spans across eastern Austria, southern Czech Republic and western Slovakia (Fig. 1). Due to its special location in this transition zone it has a complex tectonic history, which has been influenced by the changing tectonic regimes in the last 18 Ma (see Lee & Wagreich (2016) and references therein). The Vienna Basin was formed on top of thrust sheets in the Eastern Alps (Hölzel et al. 2010), which have been caused by the convergence of the African and European plates. Lateral extrusion to the East during the late Oligocene and early Miocene (Ratschbacher et al. 1991;Wölfler et al. 2011) was associated with the formation of strike-slip faults such as the sinistral strike-slip Salzach-Enns-Mariazell-Puchberg (SEMP) and Mur-Mürz Line (MML) faults ( Fig. 1). All these factors have played a role in the formation of the Vienna Basin, leading to a relatively complex tectonic structure.
Previous studies in the region have investigated the Vienna Basin and its surroundings extensively (e.g. Brix & Schultz 1993;Wessely 2006;Behm et al. 2007;Brückl et al. 2010;Ren et al. 2013;Behm et al. 2016). Brix & Schultz (1993) and Wessely (2006) give insight into its geological structure using surface geology classifications, borehole data and data from non-public seismic surveys. Behm et al. (2007) present a crustal 3-D P-wave velocity model based on wideangle reflection and refraction data. They find low velocities associated with the Vienna Basin, but the data did not allow distinguishing it from the Little Hungarian Plain (LHP). Brückl et al. (2010) studied Moho depths using controlled source seismic experiments and elastic plate modelling, and report Moho depths of ∼30-40 km in the study area. Ren et al. (2013) and Behm et al. (2016) both present 3-D shear-velocity models based on ambient-noise tomography. Ren et al. (2013) use data of the regional Carpathian Basins Project (CBP; Dando 2011) and South Carpathian Project (Ren et al. 2012), combined with permanent station data. The authors image low-velocity zones associated with sedimentary basins-including deep low velocities beneath the Vienna Basin-with a lateral resolution of ∼60 km. Behm et al. (2016) use data of the Alpine Lithosphere and Upper Mantle PASsive Seismic Monitoring (AL-PASS; Mitterbauer et al. 2011) and CBP (Dando 2011) projects and measure low Rayleigh and Love group velocities in the Vienna Basin at longer periods up to 20 s. The lateral resolution of these studies in and around the Vienna Basin is limited by relatively poor station coverage (Behm et al. 2016), measurements along only a few profiles (Behm et al. 2007;Brückl et al. 2010) or a broader regional focus (Ren et al. 2013). A high-resolution model of seismic velocities in the region is currently missing.
We compute a new high-resolution crustal 3-D shear-velocity model of the wider Vienna Basin region using ambient-noise tomography. This method is based on the extraction of estimated Green's functions (GFs) from interstation cross-correlations of ambient seismic noise, which allow to create virtual sources at every passive seismic station (see Campillo & Roux (2015) for a review paper). GF retrieval from ambient-noise cross-correlations has revolutionized the use of seismic arrays for imaging and monitoring purposes at various scales and is now an established technique with many proven applications (e.g. Shapiro et al. 2005;Nishida et al. 2009;Poli et al. 2011;Lin et al. 2012;Ren et al. 2013;Brenguier et al. 2014;Molinari et al. 2015;Nakata et al. 2016;Kästle et al. 2018). It allows to gather information about the structure of the Earth between two seismic stations without using an active or earthquake source, as the retrieved GFs contain broad-band information about dispersive surface waves in the microseism period band (Longuet-Higgings 1950;Hasselmann 1963). As the amount of available seismic records is only controlled by the number of stations, noise-based surface wave tomography has improved the apparent resolution of seismic velocity models by capitalizing on the recent expansion of seismological networks (e.g. Lin et al. 2009Lin et al. , 2013Boue et al. 2014;Ben-Zion et al. 2015;Roux et al. 2016). In this study, we will take advantage of the recent deployment of a dense seismic network in the Alpine region as part of the AlpArray project (AlpArray Seismic Network 2015) to image the wider Vienna Basin region with improved resolution.
In the following sections, we present the steps taken to compute the new shear-velocity model in the wider Vienna Basin region: Section 2-data used for this study; Section 3-retrieval of GFs from ambient noise; Section 4-measurement of Rayleighwave group velocities from GFs; Section 5--inversion of group velocities to regionalize the measurements; Section 6-inversion of group-velocity maps for shear-velocity structure. Finally, we will discuss our model with respect to previous seismological and gravitational studies, as well as insights from surface geology and borehole data.

DATA
AlpArray is an international project of 24 institutions across Europe (Hetényi et al. 2018). It aims at advancing our understanding of the Alpine Orogene and surrounding regions with a previously unachieved dense coverage of the entire Alps with broad-band seismometers. This will enable new studies with improved resolution. In total, the network consists of almost 700 seismic stations, comprised of ∼240 newly installed temporary broad-band stations, ∼30 ocean bottom seismometers and ∼400 permanent stations.
The data used in this study consist of continuous seismic records of 16 permanent stations (of the Austrian Seismic Network 1987;Czech Regional Seismic Network 1973;Hungarian National Seismological Network 1992; National Network of Seismic Stations of Slovakia 2004) and 47 temporary broad-band stations of the Al-pArray seismic network (AlpArray Seismic Network 2015). Fig. 1 gives an overview of the study area and the locations of the permanent (purple) and temporary (yellow) stations. The interstation distances range from 20 to 340 km-with an average station spacing of ∼40 km-and an even distribution of interstation azimuths. Faults in the area (black lines in Fig. 1) are compiled from the European Database of Seismogenic Faults (EDSF; Basili et al. 2013) and the International Geological Map of Europe (IGME5000; Aschk 2005). The available seismic records range from 0.5 up to 2 yr in length and have been recorded between February 2015 and April 2017.

A M B I E N T -N O I S E C RO S S -C O R R E L AT I O N S
In this study, we measure interstation surface wave traveltimes on estimated GFs, extracted from interstation cross-correlations of ambient-noise. We estimate group velocities from these traveltimes assuming great circle propagation. The group velocities are then used to image the crustal structure using a tomographic inversion procedure. In this section, we discuss the retrieval of GFs. The continuous seismic records are processed in two major steps to compute the estimated GFs: Pre-processing of the records and GF retrieval.

Pre-processing
Pre-processing aims to render the resulting cross-correlation functions (CCFs) more stable and closer to the true GFs by removing transient sources (e.g. earthquakes) from the continuous seismic records. The wave field produced by a transient source is not diffuse and may introduce spurious arrivals in the estimated GFs (Bensen et al. 2007). To be able to remove these signals, we divide the continuous records into smaller time windows. These time windows need to be long enough, so that the diffuse wave field is sufficiently well-recorded on any two seismic stations. They also should be reasonably short to not remove too much data when removing transient sources (Seats et al. 2011). We tested common pre-processing methods, such as windowing, whitening (Bensen et al. 2007, and references therein), and one-bit-normalization (Cupillard et al. 2011), extensively. To determine the final pre-processing scheme, we compared the signal-to-noise ratio (SNR) of the resulting CCFs for 384 combinations of pre-processing parameters and chose the scheme that yields the highest SNR while providing stable increase in SNR with the number of stacked days (see Supporting Information Section S1). We define SNR as the peak amplitude divided by the standard deviation in a noise-window, where the noise-window is the last 20% of a given CCF. We decided to pre-process the data in the following eight steps: (1) Remove the instrument response for each station.
(2) Remove glitches in the signal by clipping the amplitudes at 15 times the standard deviation of each daily trace.
(4) Remove subtraces that contain (i) large gaps (more than 20% of the subtrace) or (ii) transient sources, detected by an energy threshold. The threshold is defined such that the mean energy of the 30-min subtrace may not exceed 2.5 times the mean energy of the original 24 hr trace of the same day.
(5) Whiten the spectrum of the subtraces, using a water level, to reduce the impact of amplitude variations on the measurements (Bensen et al. 2007).
(6) Dampen remaining transient signals (e.g. small earthquakes that may have passed the energy threshold) by clipping the amplitudes at four times the standard deviation.
(7) Apply a taper to the edges of the subtraces to prevent border artefacts in the cross-correlations.
(8) Downsample all records to 4 Hz to reduce further computational cost.

Green's function retrieval
We extract estimated GFs from the pre-processed subtraces by crosscorrelating and stacking them. For each station pair, we first crosscorrelate the remaining subtraces after pre-processing and then stack all subtraces linearly. We compute daily stacks for quality control and to help identify corrupted data. In Fig. 2 (top views), we show the ZZ-component daily stacks as correlograms, bandpass-filtered between 5-25 s and normalized, for two representative station pairs.
Both examples show arrivals on the causal and acausal parts of the cross-correlation that are stable in time. The station pair KRUC-A014A ( Fig. 2a) is aligned SW-NE and shows similar energy levels in the causal and acausal parts of the cross-correlation. Additionally, we see a slight change in amplitude between days ∼300 and ∼400, concurrent with the summer months. Note that day 0 is the first day of simultaneously available data and not related to days of year. In Fig. 2(b) we show the station pair A014A-A020A, which is aligned N-S and shows strong asymmetry of the causal and acausal parts. Still, this example shows stable phases in both the causal and acausal parts of the cross-correlation. We observe a change in amplitude with time for roughly 100 d (Days ∼200 to ∼300), much more pronounced than in the symmetric example (Fig. 2a). This observation, consistent throughout the whole data set, leads us to assume a strong noise source in the N direction for most of the year and a relative weakening of that source during the summer months (Stehly et al. 2006;Juretzek & Hadziioannou 2016).
One of the main assumptions in GF retrieval is a homogeneous noise-source distribution (e.g. Shapiro & Campillo 2004), which is almost never achieved, except in designed experiments with a controlled noise-source distribution (e.g. Roux et al. 2004). As illustrated above, our study is also affected by a non-uniform noisesource distribution. Therefore, we stack the daily estimated GFs to sample the different noise-source distribution regimes of the seasons and reduce possible effects of seasonality on the traveltime  measurements . For each station pair, we stack all available days to retrieve a final stack (Fig. 2, middle views). In the two shown examples (Fig. 2), we observe a difference in frequency content on the causal and acausal parts of the final CCF-also visible in the daily stacks-that results from the different frequency content of the noise sources in opposite propagation directions (Yang & Ritzwoller 2008). To mitigate this effect and broaden the frequency content, we fold the cross-correlations (Fig. 2, bottom views) to merge low and high frequency information (Verbeke et al. 2012) for the following traveltime measurements.
We compute the full cross-correlation tensor, that is, we crosscorrelate all nine component pairs of Z, N, E for each station pair (ZZ, NN, EE, ZN, NZ, ZE, EZ, NE and EN). Finally, we rotate the cross-correlation tensor to receive radial, transversal and vertical components (ZZ, TT, RR, ZR, RZ, ZT, TZ, RT and TR). We retrieve a total of 17 577 estimated GFs (9 component pairs times 1953 station pairs).
In Fig. 3 we show the combined full cross-correlation tensor for all stations. Each component is represented as a plot of all final folded cross-correlation stacks, sorted and binned by interstation distance (1 km bins). We identify clear wave trains on the components associated with Rayleigh waves (ZZ, RR, RZ and ZR) and Love waves (TT) according to their respective polarization. The observed Love waves with velocities ∼3.7 km s −1 are generally faster than the Rayleigh waves at ∼3.3 km s −1 , as indicated by the steeper slope of the Love wave train. The energy on the Rayleigh wave components is not uniform across components with the ZZ component showing a clearer signal than the ZR, RZ and RR components, where the RR component has the weakest signal. This lack of energy is explained by the sensitivity of the horizontal records in this study to the limitations of temporary installations . Despite such limitations, the Rayleigh waves are identifiable on all four components. On the ZZ component, we also see very fast arrivals near 0 seconds lag time for distances smaller than ∼80 km. These fast arrivals are also visible, although much weaker, on the ZR and RZ components. They may be related to near-vertical incident body phases (Pedersen 2017) and need to be taken into account to not drastically overestimate interstation Rayleigh wave velocities. Our selection criteria presented in Section 4 already remove those measurements, without specifically tuning them for this. The cross terms (TR, RT, ZT and TZ) also carry some energy, but no clear arrivals can be identified. In this study, we focus on the analysis of Rayleigh waves and measure group velocities on all four relevant cross-correlation tensor components (ZZ, RR, RZ and ZR). We do not analyse Love waves because they would require an adapted dispersion curve measurement and selection process (see Section 4), and should ideally be inverted jointly with Rayleigh waves to receive an anisotropic shear-velocity model (Jaxybulatov et al. 2014;Mordret et al. 2015), which is beyond the scope of this paper.

R AY L E I G H -WAV E G RO U P V E L O C I T I E S
We use classic surface wave tomography, which is based on the frequency-dependent wave velocity (dispersion) of surface waves Downloaded from https://academic.oup.com/gji/article-abstract/215/1/102/5046455 by guest on 27 July 2018 (e.g. Stein & Wysession 2003) to image seismic velocity structures in the study area. Because the estimated GFs are dominated by surface waves in the period range of 5 to 25 s (Figs 2 and 3), we can obtain seismic velocity information at crustal depths (e.g. Stein & Wysession 2003). In this section, we describe how we measure and select the Rayleigh group dispersion curves that will be used for the inversion scheme.
We measure group velocities using the Multiple Filter Analysis, first introduced by Dziewonski et al. (1969). We isolate waves of certain periods from the Rayleigh wave train by bandpass-filtering the records with a narrow Gaussian filter around a centre period. The maximum of the envelope of the resulting signal is picked as an estimate for the group arrival time of that centre period. We perform these measurements on all four tensor components for 1953 station pairs, resulting in 7812 Rayleigh group-velocity dispersion curves.
To ensure that we use mainly high-quality fundamental mode measurements, we employ a set of selection criteria. This is necessary, because the picking algorithm itself does not discern between modes and does not check the quality of the measurements (see also Zigone et al. 2015). For each centre period, we employ one station-based criterion and four component-based criteria: (1) There must be at least two wavelengths of the measured wave between the stations (λ = v measured · T, with the period T). This ensures that the wave properly samples the medium.
(2) Measured velocities may not deviate strongly from the mean of all four components (± 10%). We keep only components where the measured velocity does not exceed this threshold. This removes outliers and ensures that the measured velocities on the four components for a single station pair match, thereby avoiding measurements that are biased by noise sources.
(3) The energy of the arriving group must be greater than 1% of the maximum energy measured for that station pair for any period. We estimate the group energy as the peak amplitude of the envelope for a given centre period. We set 1% of the highest measured group-energy for a given component across all centre periods as the threshold. This removes poorly constrained measurements due to very low arrival energy.
(4) The SNR of the filtered cross-correlation used to measure the group velocity must be greater than 4.
(5) Finally, the measurement on the ZZ component must be part of the final set of measurements for a given station pair. This acts as a weighting factor for vertical component measurements. We found that the horizontal components are often less well-resolved, partly due to a large number of stations being temporary installations in sedimentary settings (see Section 3.2).
If at least three of the four components pass all the tests and ZZ is one of them, the mean velocity of these components for that station pair is preliminarily accepted as the interstation group velocity of the given centre period.
We limit the dispersion curves to the range of 5 to 25 s. Below 5 s the measurements are dominated by higher modes and we do not retain enough fundamental mode measurements. Above 25 s the measurements are poorly resolved and are mostly eliminated due to selection criteria.
The original data set consists of 315 208 group-velocity measurements (grey histogram in Fig. 4 a and distribution in Fig. 4b). Of those, the selection criteria remove 119 389 (37.9%) measurements leading to the red distribution in Fig. 4(a). The selection criteria successfully retain the slightly bimodal distribution characteristics, which represent the two distinct dispersion curve trends in the range of 5 s ≤T ≤ 17 s (Fig. 4b), while eliminating most higher-mode measurements and low-velocity artefacts. We remove the remaining outliers, sometimes related to higher modes by keeping only measurements within one standard deviation from the original set of all measurements (before selection criteria were employed) for 5 s ≤T ≤ 7 s, and within two standard deviations for 7.5 s ≤T ≤ 25 s for each centre period respectively. Lower periods are subject to a stricter threshold, because they are less reliable, more likely to be influenced by higher modes, and to stabilize the inversion results. This threshold eliminates an additional 7969 (2.5%) measurements, mostly at the edges of the distribution. In total, we keep 187 850 (59.6%) (black histogram in Fig. 4 a and distribution in Fig. 4c) of the initial 315 208 measurements.
The resulting averaged interstation group velocities (see Supporting Information Section S2) show a spatially coherent trend of faster velocities in the West and slower velocities in the East.

G RO U P -V E L O C I T Y I N V E R S I O N
Combining all measurements of interstation group velocities for a certain period allows to invert for the group velocities associated with regions (cells) instead of paths. We regionalize the measurements by following the inversion routine of Barmin et al. (2001) to obtain isotropic group-velocity maps.
The standard forward problem is posed in matrix notation d = Gm, where the data vector d = t m − t syn consists of the traveltime differences between the measured traveltimes t m and the synthetic traveltimes t syn for a given initial model for each path. The matrix G contains the traveltimes for each path in each cell of the initial model. We choose cells that are 5 km × 5 km in size to balance lateral resolution with the number of measurements per cell and only invert cells with at least three crossing paths. This results in the group-velocity model m = (u − u 0 )/u 0 , with the initial group velocity u 0 and the group velocity after inversion u.
The inversion routine is based on the minimization of a linear combination of data misfit, model smoothness F(m), and convergence speed to the initial model for cells with few measurements H(m) The model smoothness function F(m) is a spatial Gaussian filter with the correlation length σ , given as The third term H(m) describes a weighted norm of the model, which is effective for sparsely sampled cells. It is given as an exponential function with the number of paths crossing the cell ρ and a weighting factor λ.
The inversion is controlled by a total of four regularization parameters. The factors σ and α control the model smoothness. λ and β control the weighted norm. Thanks to the favourable station distribution, path coverage is mostly even and the factors λ and β have only marginal impact on the inversion results (Supporting Information Fig. S3). Therefore, we focus on determining proper model smoothness parameters.

Determination of regularization parameters
Usually, regularization parameters are chosen by a so-called L-curve analysis (e.g. Hansen & O'Leary 1993;Stehly et al. 2009;Mordret et al. 2013). In the L-curve analysis, the inversion is performed for a set of chosen values for a single parameter, while the other parameters are fixed. The variance reductions (or a similar measure) for each inversion result are plotted versus that parameter and a value is picked near the maximum curvature of the resulting L-shaped curve. This analysis aims to give an objective measure of the tradeoff between overfitting and underrepresenting the data. The choice of values for the fixed parameters is arbitrary at first. Additionally, the values of maximum curvature for each parameter are interdependent on the choice of the other parameters. An iterative process of alternating between fixed and varied parameters can help find proper values for the regularization parameters (Hansen & O'Leary 1993).
Here, we propose a 2-D L-curve analysis. Instead of fixing all parameters except one, we fix the parameters that have minimal influence on the inversion in our specific case (λ and β, see Supporting Information Section 3). We vary σ and α simultaneously and retrieve a 2-D surface of variance reduction in parameter space (Fig. 5a). We plot one slice in each direction to illustrate the relationship to a standard L-curve analysis (right and bottom view of Fig. 5a). Fig. 5(b) shows the Gaussian curvature of that surface and slices at the same values for σ and α (right and bottom view). Negative values of curvature are set to 0, as they only appear as artefacts at the edges of the parameter space. This 2-D L-curve analysis, similar to standard L-curve analysis, does not aim to give a final objective answer to the optimization problem at hand (overfitting versus underrepresenting data). It still requires subjective expert judgement for the final choice of regularization parameters, which is not ideal, but still widely used in seismic tomography. We pick the regularization parameters near the maximum of the surface curvature towards lower variance reduction to avoid overfitting the data (σ = 8 km, α =20).

Group-velocity resolution analysis
To interpret and further use the group-velocity maps, estimating their resolution is critical. First, we show path-density maps for selected centre periods (Fig. 6). In Fig. 6(a), we show the path density for 5 s centre period. We achieve an average ∼20 paths per cell in the western part of the study area. In the eastern part, we lose measurements due to dispersion curve selection criteria, because higher mode measurements are more common for paths crossing sedimentary basins and the SNR of the horizontal components is lower for temporary stations in sedimentary settings. At 15 s centre period (Fig. 6b), we observe an even distribution throughout the study region with around 25 paths per cell in most of the centre area. We remove relatively few measurements at this period. In Fig. 6(c), we observe reduced path coverage averaging at ∼15 paths per cell in the centre region for 25 s centre period. Here, we mostly remove measurements due to the group-energy and interstation distance thresholds.
Additionally, we present resolution-length maps for 15 s period (Fig. 7). We define the resolution length as the distance at which the value in the resolution matrix is decreased to half (Barmin et al. 2001;Stehly et al. 2009;Zigone et al. 2015). Because the spatial projection of the individual resolution matrices for each cell are not symmetric, a best and a worst direction exist. We show the mean resolution length (Fig.7a), the resolution length in the best direction (Fig.7b), and the resolution length in the worst direction (Fig.7c) for each cell. The mean resolution length is ∼10 km in the centre of the study area and ∼20 km at the edges (Fig. 7a). In the best direction (Fig. 7b), we see resolution lengths of 6-10 km for most of the study area. In the worst direction (Fig.7c), the resolution length still reaches 12 km in the centre, while dropping to 30 km at the edges. Therefore, we can reliably interpret structures that span at least three cells (15 km length) for most of the study area.

Group-velocity maps
We show our final group-velocity maps for 5, 15 and 25 s centre period (Fig. 8). As major features, we observe two separate low-velocity structures in the eastern part of the map and a more homogeneous high-velocity anomaly in the northwestern part of the map, which qualitatively match the velocity trends expected from topography (Fig. 1). We identify the low-velocity bodies as the Vienna Basin in the centre of the study area (a in Fig. 8) and the LHP in the SE (b in Fig. 8). While the Vienna Basin is clearly visible over the whole period range from 5 to 25 s, the LHP fades away at 15 s. The edge of the Vienna Basin is marked well by the major known faults in the area (black lines in Fig. 8). The western edge of the Vienna Basin seems to move towards east with increasing centre period. Additionally, we observe a smaller low-velocity anomaly at the western edge of the study area (c in Fig. 8), that seems to be bounded by the Alpine Front (AF) to the north and the SEMP fault to the south. The high-velocity anomaly at 5 s centre period in the northwestern part of the study area is identified as the Bohemian Massif (d in Fig. 8). This anomaly is consistently observed at all available centre periods.

S H E A R -V E L O C I T Y I N V E R S I O N
To gain insight into the depth extent of the observed velocity anomalies, we invert for shear-wave structure using the linearized inversion routine of Herrmann (2013). For each cell, we extract all available measurements from the group-velocity maps and combine them to construct new local 1-D group-velocity dispersion curves. Each 1-D curve is then inverted independently and recombined with the other cells to construct the 3-D shear-velocity model.
Because the inversion scheme is linearized, the results may be heavily influenced by the initial model. It is therefore crucial to explore the dependence of the inversion results on the initial model and to choose a proper model. We tested several models, including global models like IASP91, models with constant velocity and published regional models (Ren et al. 2013;Behm et al. 2016). First, we construct a mean model from Behm et al. (2016) by averaging the published 3-D velocity model in the study area to retrieve a representative 1-D model. From Ren et al. (2013) we extract a 1-D velocity model located in the Vienna Basin (at 48.5 • N, 17 • E). To test these initial models (Fig. 9a), we construct a representative dispersion curve by averaging all available 1-D dispersion curves from our group-velocity maps (blue line in Fig. 9b). In Fig. 9 we show the results of inverting a variety of highly different initial models (blue lines in Fig. 9a) in terms of resulting velocity models (black lines in Fig. 9a) and their fitted dispersion curves (black lines in Fig. 9b). We find that the velocities in the depth range of 5-20 km are only marginally influenced by the initial models. Therefore, the results in this depth range seem robust. We choose the mean model from Behm et al. (2016) as the initial model, because all reasonable models (i.e. not constant with depth) show similar results in this depth range and this model was derived near the study area. The model is made up of 42 layers, each 1 km thick, with a half-space beneath. Shear velocities range from 3.1 km s −1 in the top layer to 4.2 km s −1 in the bottom layer, which is extended to the half-space. It does not contain discontinuities. Therefore, we cannot find clear discontinuities at single layers in the final shear-velocity structure (Herrmann 2013).
In Fig. 9(c) we show the selected initial model (mean model from Behm et al. (2016)) and the resulting group-velocity depth sensitivity kernels (Fig. 9d) for the measured period range in this study.
These kernels dictate the depth resolution of the shear-velocity inversion and give insight into the expected resolved depths. The period range in this study (5 s ≤T ≤ 25 s) is sensitive to the top 30 km (Fig. 9d).

Shear-velocity depth resolution
The misfit statistics between the measured local 1-D dispersion curves and synthetic dispersion curves, computed from the final inverted shear-velocity structure, are provided in Supporting Information Section 4. The mean standard deviation of group-velocity misfit for all periods is 0.037 km s −1 with no single measurement deviating more than 0.21 km s −1 . This illustrates a good match between synthetic and observed dispersion curves.
We normalize the resolution matrix of each individual cell and average all of them to compute the average normalized resolution matrix of all cells (Fig. 10). This resolution matrix contains the weights of the linear relationship in which each solution parameter is derived from the weighted averages of nearby true-model parameters (e.g. An 2012). This resolution matrix is useful to measure the solution obtainability for each layer-giving insight into the resolved depths-along with the quality of the inversion based on the degree to which the matrix approximates the identity matrix. The depths in which we achieve good resolution are controlled by the available group-velocity measurements. As previously noted, they are limited mainly by higher modes at shorter periods and poorquality measurements at longer periods. We find good resolution in depths of 4-20 km, indicated by a roughly symmetric resolution matrix in this depth range. At shallower depths our model is likely to overestimate velocities, because they are heavily influenced by the higher velocities at depths around 5 km. At greater depths our model is likely to underestimate velocities somewhat, and low-velocity zones may blur into these greater depths.

Shear-velocity maps
In Fig. 11 we show selected depth slices of our final 3-D shearvelocity model. We provide the model online as Supporting Information. The maps display the same major features as the groupvelocity maps (Fig. 8): Vienna Basin (a in Fig. 11), LHP (b in Fig. 11) and Bohemian Massif (e in Fig. 11). These maps allow to interpret the depth extent of the observed velocity structures.
The Vienna Basin is imaged as a low velocity feature with a lateral extent of ∼80 km across and ∼150 km along the SW-NEstrike of the major faults in the region at 4 km depth (a in Fig. 11). The NW edge of the Vienna Basin and transition to the Bohemian Massif is well-delineated by the AF. The SE edge is marked by the Southern part of the complex Vienna Basin transfer fault system (VBTFS), which delimits the end of the Vienna Basin towards the Leitha Mountains, Little Carpathians and LHP. The lowest shear velocities we observe are located just north of the border triangle of Austria, Czech Republic, and Slovakia with 2.14 km s −1 . At 8 km depth, the SW Vienna Basin is no longer imaged, while the NE part is still clearly mapped. The NE part still shows the lowest velocities in the model at that depth, but the location of the minimum is just East of the border triangle. The NW edge of the basin is no longer delineated as clearly by the surface expression of the AF, the edge shifts ∼8 km towards SE. The SE edge, on the other hand, still seems to be marked quite well by the VBTFS. At 12 km depth, Downloaded from https://academic.oup.com/gji/article-abstract/215/1/102/5046455 by guest on 27 July 2018  the low-velocity body has shifted further towards East with the AF being no longer associated with the NW edge of the Vienna Basin. Still, the VBTFS delimits the SE edge of the Vienna Basin. At 16 km depth, the NW edge of the Vienna Basin is still dipping further towards east, while the SE edge does not move. At 20-24 km, there is a low-velocity anomaly remaining, no longer defined by the lowest velocities for these depths, but instead it shows comparable velocities to other low-velocity features in the study area.
The LHP is imaged as a low-velocity feature with an extent of ∼250 km along SW-NE and ∼120 km across from the Southern edge of this model to the eastern edge at 4 km depth (b in Fig. 11). The lowest velocity is found in the centre of the LHP with 2.36 km s −1 . To the NW, the LHP is limited by the Little Carpathians and the Leitha Mountains. To the SE, the LHP transitions into the Bakony mountain range, which we image with higher velocities (c in Fig. 11). At around 12 km depth, the low velocities beneath the LHP seem to connect to the low velocities beneath the Bakony. At greater depths (16 -24 km), there is no clear low-velocity feature remaining that we would associate with the LHP.
The eastern edge of the northern Calcareous Alps (d in Fig. 11) is imaged as a shallow low-velocity feature, visible at 4-8 km. It is located between the AF to the north and the SEMP fault to the south.
The Bohemian Massif in the NW is represented as widespread high-velocity feature (e in Fig. 11). The SE edge of the Bohemian Massif is marked by the Diendorf Fault (DF) to the Molasse basin in the South at 4 km depth. As the Molasse basin becomes narrower towards East, the edge of the Bohemian Massif lies closer to the AF. At shallow depths (4-8 km), the Bohemian Massif is seen as a relatively homogeneous high-velocity feature, which becomes more complex towards greater depths. At 16-24 km depth, an elongated, ∼40 km wide high-velocity body is visible, which approximately follows the surface expression of the AF (f in Fig. 11). Beneath the Bohemian Massif, slightly reduced velocities are visible at these greater depths.

Shear-velocity profiles
In Fig. 12 we show four depth profiles cutting through the final shearvelocity model. Three profiles are crossing the major structures in our study area at different latitudes (A, B and C in Fig. 12a). Profile D is aligned SW-NE along the strike of the major faults in the region. On the map view (Fig. 12a), we mark the locations of known boreholes (Brix & Schultz 1993;Wessely 2006) that have reached the crystalline basement ( ) and some that have not ( ). Boreholes are plotted at those depths at which they reached the crystalline basement ( ) or the depth at which they were terminated ( ), if they did not reach the basement. We will discuss the boreholes along with other additional observations in Section 7. In the profiles (Fig. 12), main features are labelled with abbreviations at the top, and intersection points with major faults are marked as bold vertical lines.
In profile A (Fig. 12), we see the Bohemian Massif dipping mildly (∼20 • ) towards SE below the Vienna Basin. The low-velocity signature of the Vienna Basin is visible up to depths of 10 km in the NW (at 80 km distance) and up to 20 km in the SE (at 130 km distance). The low-velocity anomaly under the LHP (at 200 km distance), on the other hand, is only visible at shallower depths up to 8 km. The two sedimentary basins are separated in the profile at the Little Carpathians (at 140 km distance), which aligns with the end of the SE extent of very low velocities below 2.3 km s −1 .
In profile B (Fig. 12), the same major structures are visible. The Bohemian Massif is dipping below the Vienna Basin with a dip angle of ∼20 • . The deep low velocities around 15 km depth beneath the Vienna Basin are generally higher than those seen in profile A (v s ∼3.0 km s −1 vs. v s ∼2.7 km s −1 ). The SE edge of this low-velocity anomaly aligns with the SE edge of the Little Carpathians. The LHP shows low velocities (v s ≤3.0 km s −1 ) up to 12 km depth.
In Profile C (Fig. 12), the low velocities of the Vienna Basin can be seen only at shallow depths of less than 5 km. No lowvelocity feature at greater depths beneath the Vienna Basin is visible. Here, we image the expected lower sedimentary thickness of the SW Vienna Basin, compared to the NE Vienna Basin (Wessely 2006). The transition from the Vienna Basin to the LHP in the velocity model aligns with the Leitha Mountains. Between the Bohemian Massif and Vienna Basin, we see very shallow (≤3 km) low velocities, which seem to be associated with the Molasse basin. Note that low velocities do not show up in the near-surface crystalline structure of the Bohemian Massif. The Bohemian Massif itself is harder to follow dipping towards SE and to separate from the overlying structures compared to profiles A and B (Fig. 12).
Finally, in Profile D (Fig. 12), we show a cross-section along the strike of the major faults from the northern Calcareous Alps in the SW to the western Carpathians in the NE. The northern Calcareous Alps show very shallow low velocities, similar in depth extent and magnitude to the Molasse basin (see Profile C). The Vienna Basin, NE of the northern Calcareous Alps, is dipping very mildly (∼10 • ) towards NE until 130-140 km distance. There, we see a steep increase of the low-velocity feature towards greater depths (∼20 km), which is prevalent towards NE into the western Carpathians. From 200 to 230 km distance, we see the velocities of this feature increase from v s ∼2.7 km s −1 to v s ∼3.0 km s −1 at ∼10 km depth.

D I S C U S S I O N
The new 3-D shear-wave velocity model we present in this study correlates with several previously mentioned geological features. For discussion, we will compare our results with previous seismological studies (Tomek & Hall 1993;Behm et al. 2007Behm et al. , 2016Ren et al. 2013;Hrubcová &Środa 2015), as well as ground truth from borehole data (Brix & Schultz 1993;Wessely 2006, and Fig. 12), and with gravity field measurements (Bonvalot et al. 2012). We will not try to interpret the tectonic evolution of the Vienna Basin, surrounding region and underlying structures in detail, and wave the new insights provided by this model to be assessed by more qualified colleagues. In the following, we will interpret the bottom of the observed low-velocity features though, as the interface between sedimentary rocks and crystalline basement.

Vienna Basin
We map the Vienna Basin as a low-velocity feature, 80 km × 150 km in size near the surface at 4 km depth (a in Fig. 11). Its lateral extent is well-delimited by known major faults in the area: towards SE, the southern parts of the complex VBTFS mark the edge towards the Leitha Mountains and Little Carpathians, while the AF coincides with the transition between Vienna Basin and Molasse Basin to the NW (Figs 1a, 12 a and 13a). It matches the lateral extent expected from surface geology (Fig. 13a) and is mapped as a trough with ∼50 mGal in the Bouguer anomaly (a in Fig. 13b), which can be explained by shallow low-density rocks.
Other seismological studies also find low velocities in the area from P-wave wide-angle and refraction seismics (Tomek & Hall 1993;Hrubcová et al. 2005;Behm et al. 2007), ambient-noise tomography using data of the southern Carpathian Project and CBP (Ren et al. 2013), and ambient-noise tomography using data of the ALPASS project (Behm et al. 2016). Behm et al. (2007) map the Vienna Basin as a (P-wave) low-velocity feature, which is not clearly separated from the LHP. They use a grid spacing of 20 km for the regionalization of Pg phase picks. Ren et al. (2013) also image the Vienna Basin as a low-velocity anomaly, the lateral extent roughly matching our model. They use a 0.2 • × 0.2 • grid to regionalize the data and properly resolve structures as small as 60 km in lateral extent. The model of Behm et al. (2016) is limited with respect to the station distribution, especially towards the NE parts of our study area and cannot map the lateral extent of the low velocities in the Vienna Basin.
Our model improves the lateral resolution of seismic velocities in the Vienna Basin and surrounding region greatly. The mean resolution length is ∼15 km for most of the study area in a grid with 5 km × 5 km cells, roughly increasing the resolution by a factor of 4 compared to previous studies (Behm et al. 2007;Ren et al. 2013;Behm et al. 2016). This assessment is further supported by our ability to image the Southern Vienna Basin at 4 km depth (g in  Fig. 11), which we find to be a ∼20 km wide structure at this depth, matching geological interpretations (Wessely 2006).
The improved resolution of our model also allows us to image and interpret the deeper structure (10-20 km) of the Vienna Basin and surrounding region with greater precision than previously possible. The SE Vienna Basin appears as a shallow low-velocity feature with depths of up to ∼5 km, matching ground truth from boreholes Laaer Berg 1 (LB1 in Fig. 12 profiles C, D) and Rauchenwarth 1 (Rau1 in Fig. 12, profile C). Towards NE, the Vienna Basin bottom dips mildly down to ∼10 km depth before a steeper drop is visible in the velocity model NE of borehole Aderklaa Ultratief 1 (AUT1 in Fig. 12, profile D), down to depths of ∼15-20 km. The low velocities at these depths are consistently observed in the NE Vienna Basin and Western Carpathians (Fig. 12, profile D). The reflection profile 8HR (Tomek & Hall 1993) crosses from the SE edge of the Bohemian Massif into the Western Carpathians. It images an SE dipping reflector that is clearly visible down to depths ∼10 km below the Vienna Basin. Further towards SE, where we find a low-velocity anomaly at greater depths, the profile is less conclusive. Here, Ren et al. (2013) find low velocities beneath the Vienna Basin in depths of up to ∼16 km, matching our observations. The refraction profile CEL09 (Hrubcová et al. 2005) of the CELEBRATION2000 project (Guterch et al. 2003) crosses the study area ∼50 km South of the deep low-velocity anomaly in our model and reveals the same SEdipping interface. They find low P-wave velocities up to depths of 8 km just below the central Vienna Basin, consistent with our model.
The available data from most boreholes match well with the shear velocities in our model (Fig. 12). When a borehole reaches the crystalline basement, our model generally transitions from v s ≤ 2.9 km s −1 to higher velocities at similar depths (±∼2 km). This is illustrated by the boreholes Laa 1 (L1 in Fig. 12, profile B), Mauerbach 1 (Mau1 in Fig. 12, profile C) and Aderklaa Ultratief 1 (AUT1 in Fig. 12, profile D), which have reached the crystalline basement at varying depths.
In the NE Vienna Basin, no boreholes have reached the crystalline basement. This matches our model, which shows low shear velocities at the locations and termination depths of boreholes Maustrenk Ubertief 1, ZistersdorfÜbertief 2A (MÜT1 and ZÜT2A in Fig. 12, profile B) and Zavod 93 (Z93 in Fig. 12, profile D) (Wessely 2006). While the crystalline basement is expected ∼1 km below the termination depths of MÜT1 and ZÜT2A, the crystalline basement below Z93 is expected well below 10 km (Wessely 2006). Our model corroborates this and the transition to velocities ≥ 2.9 km s −1 with depth ( Fig. 12, profile B, D) seems to mark the expected crystalline basement depths quite well. Therefore, our model suggests a very deep crystalline basement around 20 km depth at borehole Z93 (profile D in Fig. 12). We are aware of only few other works exploring these depths below the NE Vienna Basin, namely the tomography presented by Ren et al. (2013), which lacks the resolution to give better information on the depth extent of these observed features, and reflection profile 8HR presented in Tomek & Hall (1993), which is not clearly imaging any interfaces at these depths in the specific area. Additionally, Wessely (2006) interpret the sedimentary rocks underlying the Vienna Basin to extend to increasing depths towards NE, matching our observations, but their interpretations do not extend deeper than ∼10 km. Notably, our model does not correlate well with the ground truth from borehole Berndorf 1 (Be1 in Fig. 12, profile D). It is located at the edge of the Southern Vienna Basin and reaches the basement at a depth of 5-6 km (Brix & Schultz 1993). Our model shows shear velocities ≥ 2.9 km s −1 already at 2 km depth. This apparent mismatch may be explained by the location of this borehole in the complex transition zone between the horizontally layered sedimentary basin and the hard rock of the Alpine orogene, illustrating the limitations of our methodology. The inversion routine is based on fitting dispersion curves computed in 1-D layered media. This assumption is less valid in complex, deformed rheology. Still, this approximation seems valid inside the Vienna Basin, supported by the ground truth from boreholes located within the basin.
These observations, combined with the depth resolution analysis (Section 6.1, Fig. 10) give us confidence that we properly resolve the Vienna Basin and its underlying structure in 4-20 km depth. Thus, the deep low-velocity anomaly beneath the NE Vienna Basin (Figs 11 and 12) does not appear to be an artefact introduced during GF retrieval, dispersion curve measurements, or the inversion scheme and thus probably represents a true geological feature. It is unlikely to be caused by smearing of low velocities from shallow to greater depths during the inversion, because the deep low-velocity anomaly is not directly beneath the strongest shallow anomaly, but shifted by ∼25 km towards SE.
While the observed near-surface velocities (in the top 4 km) correlate with expectation from surface geology, the shear velocities at these depths may not necessarily be properly resolved. Our depth resolution analysis (Fig. 10) suggests that these shallow velocities are dominated by the shear velocities at around 5 km depth. Therefore, shallow velocities in our model are likely overestimated, given that seismic velocities generally increase with depth.

Little Hungarian Plain and Bohemian Massif
The LHP is the second major low-velocity feature in our model. At 4 km depth it spans from the Southern to the Western edge of our model (Fig. 11). Towards SE it is limited by the Bakony mountain range geologically (Fig. 13a), which correlates with higher velocities in the shallow crust at the SE edge of our model. Towards NW the LHP ends at the Leitha Mountains and Little Carpathians geologically (Fig. 13a), which our model represents well by velocity contrasts towards higher velocities, not unlike that seen under the Bakony mountain range. The centre of the of the LHP is marked by velocities as low as 2.43 km s −1 . The lateral extent we find in shallow depths coincides well with the model of Ren et al. (2013).
The shape of the strongest low-velocity anomaly in the LHP at 4 km depth fits particularly well with the Bouguer gravity anomaly (b in Fig. 13b). Because decreased Bouguer anomalies can be associated with low density rocks (e.g. sedimentary rocks), this observation corroborates our lateral resolution estimation (Section 5.2, Fig. 7).
We image the low velocities of the LHP as deep as 10 km. It is deeper than the SE Vienna Basin (Fig. 12, profile C), but not as deep as the NE Vienna Basin and its underlying structures, which we image up to 20 km depth (Fig. 12, profile D). Other authors (e.g. Ren et al. 2013) report similar depths of this low-velocity feature.
We map the SE edge of the Bohemian Massif at the NW edge of our model. It is characterized by a relatively homogeneous highvelocity anomaly at shallow depths (Fig. 11). Towards SE we observe a velocity contrast from 3.5 to 3.2 km s −1 near the DF and a very high gradient-from 3.5 to 2.8 km s −1 over ∼20 km-further towards NE, where the Molasse basin is very narrow (Figs 1a, 12 a and 13a). With depth (e.g. at 20 km, Fig. 11), the Bohemian Massif exposes more complex velocity variations in the range of 3.0-3.5 km s −1 . Beneath the high-velocity top (v s ∼ 3.5 km s −1 ) of the Bohemian Massif, lower velocities are visible (v s ∼ 3.2 km s −1 , Fig. 12, profiles A-C). The high-velocity feature along the surface expression of the AF at depths 16-24 km (Fig. 11) seems to align with the top of the BM, dipping below the Vienna Basin (Fig. 12, profiles A-C).
The Bohemian Massif is well-mapped in the Bouguer-gravity anomaly map (c in Fig. 13b) as a positive anomaly relative to the surrounding region. We interpret this as the presence of high-density rocks (e.g. crystalline basement), which we image as high velocities. The lateral extent of the low velocities in our model and the higher gravity anomaly in the Bohemian Massif correlates well (Figs 11  and 13).
Overall, we resolve the crustal structure of the Vienna Basin and surrounding region with previously unachieved resolution. Our model is consistent with and improves upon previous seismological studies (Tomek & Hall 1993;Hrubcová et al. 2005;Behm et al. 2007;Ren et al. 2013;Behm et al. 2016), Bouguer gravity anomaly studies (Bonvalot et al. 2012), represents well-known surface geology (Aschk 2005) and matches ground truth from most boreholes in the region (Brix & Schultz 1993;Wessely 2006).
Our model may be used in the future for several further studies. In seismological applications, the model can be utilized to improve regional wave propagation modelling, as well as the location accuracy of local and regional seismicity by accounting for local and regional heterogeneities affecting wave propagation. Geology and tectonics may profit from our model to gain new insight into the deeper crustal structure beneath the Vienna Basin and surrounding regions. These insights may help to better understand the complex tectonic evolution of the Vienna Basin and Alpine-Carpathian transition zone.
There is potential to further increase the resolution and accuracy of seismic velocities in the region through several means. The station density could be improved further, which would allow even better resolution, either in the whole region or locally with dense seismic arrays (e.g. Ben-Zion et al. 2015;Nakata et al. 2016). The Alpine-Carpathian transition zone has been subject to several seismic studies, which generated plenty of non-simultaneous continuous seismic recordings. These could be utilized using the C3 technique (Stehly et al. 2008), which allows to compute estimated GFs from non-simultaneous records (Spica et al. 2016). Apart from data, there are also opportunities to improve the results by discerning between and incorporating several surface wave modes (especially at short periods), by measuring phase velocities in addition to group velocities, and by incorporating Love waves and jointly inverting Love and Rayleigh waves to derive an anisotropic velocity model.

C O N C L U S I O N S
We computed a detailed 3-D shear-velocity model of the crust in the Vienna Basin and surrounding region using ambient-noise tomography. It complements previously released studies, imaging the wider Vienna Basin region (Behm et al. 2007;Ren et al. 2013;Behm et al. 2016). The model provides new insight into the deep structures of the Vienna Basin and surrounding regions by achieving better resolution thanks to the favourable station distribution of the AlpArray project seismic network (AlpArray Seismic Network 2015). The main outcomes of this study are as follows: (1) We image the main geological structures in the study area clearly, the Vienna Basin, the Bohemian Massif and the LHP (Fig. 11). Their lateral extents match well with known geological features, that is, faults marking the transition between features, such as the DF, the AF and the VBTFS (Figs 1 and 13). We find prominent velocity contrasts near the surface expressions of those faults, which are expected due to the change of lithology between different geological units.
(2) Additional insight into the deep crustal structure (10-20 km) of the wider Vienna Basin region (Figs 11 and 12). We image the northern Vienna Basin and underlying features up to depths of ∼20 km. A shear-velocity contrast with depth (from ∼2.8 to ∼3.2 km s −1 ) is visible near the expected depth of the crystalline basement for most locations, where borehole data are available (±2 km). The Bohemian Massif dips below the Vienna Basin towards southeast at an angle of ∼20 • . Additionally, the Vienna Basin dips towards northeast, mildly at first and then with a steeper slope to the greater depths observed.
(3) The model we provide has a previously unachieved lateral resolution of ∼15 km for most of the study area. This improves on previous work in the area (Behm et al. 2007;Ren et al. 2013;Behm et al. 2016).

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online. Figure S1. Determination of pre-processing parameters using a representative subset of all data.  Figure S2. Interstation group velocities for selected centre periods: 5, 9, 13, 17, 21 and 25 s. For all periods, we see higher velocities in the west and lower velocities in the east of the study area. With increasing period, velocities generally get faster. Figure S3. 2-D L-curve analysis for the regularization parameters λ and β. (a) The main view shows the 2-D surface of variance reduction in parameter space. The right view represents the slice through that surface at λ = 0.4. The bottom view represents the slice through the surface at β = 5. (b) The main view shows the Gaussian curvature of the 2-D surface in (a). The other views are slices through the curvature. The values for the regularization parameters are picked near the maximum curvature (λ = 0.4, β = 5). Figure S4. Group velocity misfit distribution. (a) Group velocity misfit boxplots for each centre period, comparing measured velocities with synthetic group velocities, computed using our final shear velocity model. (b) Stacked misfits over the entire period range 5 s ≤ T ≤ 25 s. Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.