Teleseismic S-wave tomography of South Island, New Zealand upper mantle

The South Island of New Zealand, where oblique convergence between the Australian and Pacific plates occurs, is an excellent location to investigate subsurface deformation. The tectonic history is relatively simple, and a number of previous studies of the region have occurred. In this investigation, we present an S-wave tomogram of the mantle underlying the South Island that incorporates data collected from ocean bottom seismometers deployed on either side of the island. This allowed for broader and deeper resolution of mantle structure under the South Island, as well as the chance to discuss the challenges of tomographic studies in this region.

In the South Island, the Alpine fault (Fig. 1), a predominantly dextral strike-slip fault running the length of the island, is the major tectonic feature of the area and accommodates the majority of motion between the Australian and Pacific plates. Based on plate reconstructions, more than 800 km of displacement occurred between the plates during the past 45 m.y., with slip on the Alpine fault accommodating at least 450 km of this displacement (Sutherland, 1995, 1999; Sutherland et al., 2000) and possibly more than 700 km (Lamb et al., 2016). During this time, the relative motion between the Australian and Pacific plates appears to have changed little, though a component of convergence that initiated ca. 11 Ma (Cande and Stock, 2004) resulted in ∼70 km of shortening across the South Island and the uplift of the Southern Alps. Geodetic measurements show ∼40 mm/yr of right-lateral strike slip along the Alpine fault and ∼10 mm/yr of convergence perpendicular to the fault (Beavan et al., 2002).

Subduction zones to the north and south of the South Island bookend the Alpine fault. To the north, oblique subduction of the Pacific plate beneath the Australian plate is of opposite polarity to subduction of the Australian plate beneath the Pacific plate in the south (Sutherland, 1995). Convergence between the Australian and Pacific plates initiated ca. 45 Ma (Sutherland, 1995, 1999; Cande and Stock, 2004). Seismicity under the North Island to depths of ∼300 km and under the northern end of the South Island to depths of ∼250 km suggests the subducting slab turns subvertical at ∼100 km depth (Anderson and Webb, 1994). Subduction in the south is highly oblique and younger than in the north, initiating during the Miocene with ∼170 km of plate subducted since this time (Sutherland, 1995; Eberhart-Phillips and Reyners, 2001).

Previous studies found evidence for a nearly vertical, high seismic-wave–speed anomaly under the central portion of the Southern Alps that may extend to at least 200 km deep (Stern et al., 2000; Kohler and Eberhart-Phillips, 2002). These studies indicate that the high-speed anomaly is relatively symmetric and is either subvertical or dipping steeply to the west. Such studies attribute high speeds directly under the Southern Alps to thickening of the lithosphere during the convergence that initiated ca. 11 Ma. Alternatively, other investigations (e.g., Reyners et al., 2017) suggest that the high-speed structure in the mantle in the central portion of the island is the southern extent of the Hikurangi plateau, an oceanic plateau that is being subducted beneath the Australian plate.

Recently, Zietlow et al. (2016) undertook a teleseismic P-wave study of the South Island using data from an onshore-offshore seismic experiment with station spacing ∼100 km (larger aperture but also larger station spacing than in the experiments described above). This study found a lack of high P-wave speeds in the central portion of the South Island under the Southern Alps, suggesting that high-speed structure in this region does not extend deeper than 150–200 km. Zietlow et al. (2016) also inferred high P-wave speeds to depths of ∼450 km under the northwestern South Island, consistent with oblique westward subduction of Pacific lithosphere since 45 Ma. The depths that high speeds reach, combined with the known convergence across New Zealand, suggest that ∼200–300 km of slab is required below intermediate depth seismicity in order to produce sufficient weight to induce the intermediate-depth earthquakes (Zietlow et al., 2016).

The lack of resolvable high-speed structure in the uppermost mantle in the central region of the South Island in the P-wave study of Zietlow et al. (2016) is surprising. Furthermore, a surface-wave study using data recorded on the same instruments used in Zietlow et al. (2016) and the study presented here shows a high-speed mantle feature (Vs ∼4.7 km/s) in the central South Island. While resolution from S waves is typically less than that from P waves, this did prompt a question. Given the same temporary deployment of ocean bottom seismometers and on-land stations used by Zietlow et al. (2016), can an S-wave teleseismic tomographic study resolve any high-speed structure under the central South Island in the upper mantle? Additionally, comparing Vp and Vs tomograms allows for investigation of the cause of seismic heterogeneity in the region.

Consisting of a one-year deployment from late January 2009 to early February 2010, the Marine Observations of Anisotropy Near Aotearoa (MOANA) experiment comprised 30 ocean bottom seismometers (OBSs) installed off both coasts of the South Island of New Zealand (Fig. 1) and complemented the permanent New Zealand National Seismograph Network (Petersen et al., 2011). The MOANA array used ∼100 km station spacing with a full-array aperture of over 900 km NW–SE, increasing the extent of the seismic network to over four times the width of the South Island. Seventeen stations (NZ05–NZ21) were deployed to the west of the South Island on the Challenger Plateau, and nine stations were deployed to the east on the Campbell Plateau. Both plateaus are submerged continental crust of the Australian and Pacific plates, respectively. Four stations (NZ01–NZ04) lay atop oceanic crust to the west of the South Island and south of the Challenger Plateau. Amongst the 30 OBSs deployed, station NZ01 did not yield usable seismic data, and we did not recover station NZ17. To orient the horizontal components of the OBSs, we used Rayleigh wave polarizations (Stachnik et al., 2012). In addition to the OBSs, the MOANA experiment also included four temporary broadband land stations (CASS, CROE, KELY, and FREW; Fig. 1).

We examined seismograms (Fig. 2) recorded across MOANA from all events of m_b ≥ 5.5 with epicentral distances between 25° and 100° occurring from February 2009 to February 2010. Similar to other OBS teleseismic body-wave studies (e.g., Wolfe et al., 2011), we examined these seismograms in multiple frequency bands, although we found few usable events for frequencies greater than 0.12 Hz (<∼8 s). Of the seismograms examined, we observed clear S-wave arrivals on more than 50% of stations from 35 earthquakes in the filter band 0.05–0.1 Hz (20–10 s) and 12 earthquakes in the filter band 0.08–0.12 Hz (12.5–8.3 s). To help fill undersampled back azimuths and epicentral distances, we utilized SKS phases where available, finding four events in the band 0.05–0.1 Hz and one in the band 0.08–0.12 Hz with clear arrivals.

Many of the observed earthquakes that yielded clear S-wave arrivals occurred in the subduction zones offshore Fiji-Tonga, Japan, and Sumatra. So as not to preferentially weight a particular back azimuth, we separated the earthquakes into 10° back azimuth and 10° distance bins and selected one earthquake from each bin. This yielded 25 earthquakes with usable picks in the band 0.05–0.1 Hz and 11 earthquakes in the band 0.08–0.12 Hz (Fig. 3). In total, we measured 1321 S and 213 SKS arrival times (RMS ∼2.16 s) using dbxcor, a waveform cross-correlation, beam-forming program (Fig. S1 in the Supplemental Files¹; Pavlis and Vernon, 2010).

To calculate crustal thickness beneath each station, we used recent Pn travel-time measurements (Collins and Molnar, 2014) from earthquakes recorded across the MOANA array and the permanent land stations. A Pn time term is the difference in travel time between a horizontally traveling P wave in the mantle and one refracted to the surface. So, to calculate crustal thickness, we solved Equation 1 for h_c with k = 90°, ∆T equal to the Pn time term determined by Collins and Molnar (2014), h_s equal to the sediment thickness given by Whittaker et al. (2013), v_m = 8.1 km/s, v_c = 6.2 km/s, v_s = 2.5 km/s, and i and j from Snell’s Law. At stations where Pn measurements were not published (NZ02, NZ03, NZ04, APZ, DCZ, KHZ, NNZ, PYZ, and THZ), we assumed a crustal thickness given by Salmon et al. (2013).

With the calculated crustal thickness and assumed sediment thickness (Table S1 [^{footnote 1}]), we calculated travel-time corrections for our data on a by-event basis assuming a shear-wave speed in the mantle of 4.5 km/s and in the crust of 3.6 km/s for continental crust and 3.8 km/s for oceanic crust, and appropriate angles of incidence from the TauP calculator (Crotwell et al., 1999). Although shear-wave speed in sediment varies greatly with depth (e.g., Hamilton, 1976; Ruan et al., 2014), we found that varying the shear-wave speed in sediment between 50 m/s and 2 km/s did not greatly alter the final tomographic image of shear-wave speeds. Thus, we assumed a shear-wave speed in sediment of 1 km/s for the sediment correction.

After application of corrections for crust and sediment, patterns of teleseismic S-wave travel-time residuals reveal lateral variations in S-wave speeds (Fig. 4 and Fig. S2 [^{footnote 1}]) and vary in magnitude between approximately +6 and –6 s. Recordings of earthquakes that occurred to the northeast of the South Island show early arrivals at stations on land, at OBSs just offshore the east coast, and at NZ14 and NZ15 compared to all other OBSs that record late arrivals (Figs. 4A and 4B). In the northwestern part of the South Island, residuals are upwards of ∼5 s fast (e.g., NNZ). This pattern suggests that relatively high speeds underlie the northwest part of the island, though rays to these stations passing directly through the Pacific slab beneath the North Island may also contribute to these observed early arrivals (e.g., Williams et al., 2013).

For earthquakes from the south (Figs. 4C and 4D and Fig. S2a [^{footnote 1}]), residuals show late arrivals, ranging between +2 and +3 s, on most land stations, as well as on the OBSs offshore the east coast, relative to other stations. Most OBSs offshore the west coast show early arrivals. These early arrivals, with residuals <–2 s, at northwestern land stations and OBSs are again consistent with a high-speed zone beneath the northwestern part of the South Island and its offshore region. Residuals from earthquakes from the southwest (Fig. S2b [^{footnote 1}]) show a similar pattern to residuals from earthquakes to the south, though residuals in the south-central region (e.g., EAZ and MLZ) record early arrivals instead of late.

Ocean bottom seismometers deployed off the west coast record late arrivals compared to most land stations that record early arrivals for earthquakes from the northwest (Figs. 4E and 4F and Figs. S2c and S2d [^{footnote 1}]). The residuals recorded on some eastern land stations (e.g., OPZ, SYZ, and TUZ) and the OBSs off the east coast differ depending on angle of incidence. S waves from earthquakes occurring within ∼60° of the island (Fig. 4E and Figs. S2c and S2d [^{footnote 1}]) arrive early at the most eastern land stations and OBSs off the east coast, whereas those from earthquakes greater than 60° away arrive late (Fig. 4F). The different distributions of residuals for different azimuths indicate a depth limit to the high speeds under the South Island. Residuals range between approximately +3 and –3 s from earthquakes less than 60° away and between approximately +5 and –5 s from earthquakes greater than 60° away.

The region for which we performed our tomographic inversion ranges in latitude from ∼35°S to 51°S and in longitude from ∼158°E to 177°W. The model structure extends in depth from 60–600 km. We did not invert for structure shallower than 60 km because we corrected the travel-time measurements for crustal thickness and sediment prior to the inversion. A trapezoidal mesh (that accounts for the curvature of the Earth) of 70 × 71 horizontal nodes and 13 vertical nodes (representing the number of depth slices) parameterizes the structure. This resulted in a node spacing of ∼20 km in the center of the structure and ∼30–35 km toward its edges. All node spacings are less than the approximate width of the first Fresnel zone sensitivity kernel (∼100 km at the top of the model and ∼300 km at the bottom) at the frequency band of the measured S waves.

Given our residual measurements and damping and smoothing parameters, we tested model resolution using synthetic hypothesis and checkerboard tests (Figs. 5 and 6). To demonstrate the spatial coverage given our array geometry and travel-time data, we also show plots of the “hit quality” parameter (Fig. 5) at depth slices ranging from 60–465 km. “Hit quality” is a measure of both hit count and back-azimuth coverage ranging from 0–1, where a node rated “1” is sampled by several rays from each back-azimuth sextant (Schmandt and Humphreys, 2010). Measurements of “hit quality” indicate that the on-land portion of our model, as well as ∼300 km offshore either side of the central part of the South Island, is well sampled (i.e., hit quality >0.8) to 465 km deep.

The synthetic checkerboard test (Figs. 5 and 6) consisted of cubes ∼200 km on a side with alternating positive and negative amplitudes of +10% and –10% in a neutral background of the reference speeds at each depth, AK135. We chose the dimensions of the cubes to be approximately the width of the first Fresnel zone at 300 km depth. In general, with the same source locations and stations as in the inversion of real data and no added synthetic noise, the locations of the checkers are recovered, but the amplitudes of anomalies are underestimated. Maximum amplitude recovery (meaning the highest recovered value in a given depth slice) ranges from ∼20% to 88%, with recovered amplitudes greater than 50% limited to the top 165 km of the model. Amplitude recovery is ∼30% or less below depths of 330 km.

Checkerboard tests are common to test model resolution but are poor for identifying possible artifacts in the results of the inversion (e.g., Eilon et al., 2015) and typically represent the extreme limit of model resolution. Thus, we also determined expected model resolution of a 200-km-wide hypothetical structure representative of a high-speed slab spanning nearly the length of the western South Island and slightly offshore and extending either to 330 km (the depth below which the checkerboard test indicated resolution became less than 30%) or 565 km depth (Figs. 5 and 6). We set the amplitude of the high-speed anomaly at +10% of the reference speed at each depth, imbedded in a background average of –0.02%. This test also included synthetic noise, approximated as a normal distribution with standard error equal to the root-mean-square (RMS) misfit of travel times for the best-fitting model (RMS = 1.35 s) and assigned randomly to the arrival times. Typically, the recovered anomaly was larger in magnitude with the addition of synthetic noise compared to the no-noise case. In general, though, these hypothesis tests show that amplitudes are still underestimated. Furthermore, the shape of structures is only roughly recovered with some degree of vertical smearing, inherent in many body-wave tomography studies (e.g., Foulger et al., 2013).

For a synthetic slab structure reaching to 330 km depth (Figs. 5 and 6), the average recovery of speed magnitude for depths of 60–330 km reached a maximum of ∼86% of the input wave speed anomaly. In general, though, the anomalous structure is recovered in the 50%–60% range. Recovery is 50% or less toward the top and bottom of the synthetic structure. The best recovery of the structure occurs from ∼100–200 km depth. Considering a synthetic structure extending to 565 km, recovery over the depth range 60–330 km reached a maximum of ∼91% of the input wave speed anomaly. This test suggests that travel times through deep structure could leak into shallower structure and enhance inferred speeds within the shallower structure. Similar to the shallower structure, the best amplitude recovery occurs from ∼100–250 km depth, with much of the speed of the structure in this depth range ∼70%–80% resolved; however, the shape is still only roughly recovered. Below ∼350 km depth, the wave speed anomaly is recovered at 50% or less.

The hypothesis tests of a slab extending to either 330 km or 565 km depth also show that the width of the inferred high-speed zone increases with depth and appears to dip slightly to the west. This is likely the result of decreasing resolution with depth, changes in the radius of the first Fresnel zone with depth (Equation 2), and more rays arriving from the west compared to other back azimuths. As we found above, the radius of the first Fresnel zone is ∼70 km at 100 km depth and over 100 km at 300 km depth. Thus, we can resolve lateral variations in speed below ∼350 km, though not as well as those at shallower depths, and due to the increasing width of Fresnel zones, lateral dimensions of resolvable structure will also broaden with depth.

While the width of the first Fresnel zone in the upper 300 km of the model is on the order of 200 km wide, we investigated our model’s ability to resolve structure that is only 100 km wide (Figs. 5 and 6). Similar to our tests with a 200-km-wide synthetic structure, the amplitudes of a 100-km-wide structure are underestimated and the shape only roughly recovered with evidence of vertical smearing. The recovered 100-km-wide structures also, probably unsurprisingly, become wider. Thus, while structure smaller than the width of the first Fresnel zone is definitely not as well resolved, our model is still sensitive to structures with widths on the order of 100–200 km.

When considering which tomography code to use to invert the travel-time measurements, we also explored the code of Roecker et al. (2006). This is a ray-based, iterative, spherical, finite-difference tomography code. Above 200 km, we found the resolution of a checkerboard test from the code of Roecker et al. (2006) (Fig. S4 [^{footnote 1}]) to be similar to that from the code of Schmandt and Humphreys (2010) (Fig. 5) in terms of magnitude and location recovery. Below 200 km, recovery of the location of the checkerboard from Roecker et al. (2006) is clearly less than that from Schmandt and Humphreys (2010). Northeast-southwest streaking of the recovered features in the Roecker et al. (2006) resolution tests is also evident below 200 km. Considering the low frequencies at which we made travel-time measurements (0.05–0.12 Hz), it is not surprising that recovery of the magnitude and location of features deep in the model (i.e., greater than 200 km depth) is greater for the Schmandt and Humphreys (2010) finite-frequency code. The model volume sampled by long-period kernels in finite-frequency tomography would be larger than that sampled by crossing rays in a ray-based code at these depths. Despite this, we did find the resulting S-wave model from the ray-based code of Roecker et al. (2006) (Fig. S5) broadly consistent with the models produced from the finite-frequency code of Schmandt and Humphreys (2010) presented later in this paper. This suggests that ray-based and finite-frequency tomography codes yield models consistent with one another.

To summarize the limitations of our model based on these resolution tests, we can best resolve structure that is at least 200 km wide. In general, the wave speed anomaly is ∼60%–80% resolved in the depth range of 100–250 km, but it is >90% in the best resolved regions. In the most shallow portions of the model (i.e., <100 km) and below ∼350 km, resolution of the wave speed anomaly is roughly 50% or less. Structure 100–200 km wide is also recoverable, though appears wider due to the greater width of the first Fresnel zone. The shape of structure is only roughly recovered and in some cases, a vertical structure may appear slightly dipping to the west. Despite this, structure is still distinguishable in the better resolved areas of the model; however, we do not emphasize absolute anomaly amplitudes and only discuss the most robust features.

We show the S-wave tomogram of the mantle under the South Island of New Zealand (Fig. 7) in 11 depth slices ranging from 60 to 465 km (Fig. 8) and five cross sections, four traversing the central portion of the South Island perpendicular to the surface expression of the Alpine fault and one traversing the South Island NE-SW parallel to the surface expression of the Alpine fault (Fig. 9). The full 3D grid is provided in the Supplemental Files (^{footnote 1}). This model tomogram accounts for ∼68% of the ∼4.7 s² variance of the measured residuals. With an average error in the S-wave arrival time picks of ∼0.9 s and an RMS travel-time residual of ∼2.2 s, the best expected variance reduction is ∼83%. Thus, the inferred structure explains much of the variance in the observed travel-time residuals. Furthermore, the main features of the tomogram discussed below are evident in a model produced using the finite-frequency code of Schmandt and Humphreys (2010) and the ray-based code of Roecker et al. (2006), suggesting they are robust features of the travel-time data. Again, we only discuss in detail the model from the finite-frequency code of Schmandt and Humphreys (2010).

In the mantle beneath the northwestern South Island, a nearly vertical high-speed body reaches depths of ∼400–450 km (Figs. 7 and 8 and cross sections A–A′ through C–C′ in Fig. 9; labeled “1” in Fig. 7A). Speeds associated with this feature reach a maximum of at least 9%–10% higher than the reference model in the top 200 km and decrease in magnitude below this depth. While some degree of vertical smearing is evident in the hypothesis test simulating a slab of high-speed reaching depths of 565 km (Fig. 6), the test also showed that speeds at least 5% faster for a 10% fast input anomaly are recoverable to a depth of ∼400 km. Thus, while vertical smearing certainly contributes, our tomogram is also consistent with high-speed structure existing at great depths.

High speeds are also imaged in the upper mantle extending off the central west coast to a latitude of ∼43°S, though the high speeds in this offshore region are weaker and shallower, reaching to only ∼200 km depth. Such speeds are consistent with the observed travel-time data. For example, many arrival times measured from earthquakes with incidence angles at the Moho between 15°–30° at station NNZ are upwards of 6 s early (Fig. 10A). Assuming that a 6 s residual developed from passing through a 300-km-thick layer at an incidence angle of 15° (corresponding to a typical S-wave travel time of ∼69 s), the average speed must be ∼9% faster than reference speeds. Given that the resolution tests suggest a magnitude recovery of ∼60%–80%, high speeds from this feature could likely be 11%–15% faster than reference speeds. Underneath the southwestern South Island, we image another high-speed body in the upper mantle of ∼6% higher than average speeds (Figs. 7 and 8; labeled “2” in Fig. 7A), though only to depths of ∼250 km.

A region with speeds ∼5%–8% slower than the reference model lies ∼150–200 km deep to the west and south of the city of Christchurch (Figs. 7 and 8; labeled “3” in Fig. 7A). Farther south near the city of Dunedin, a similar low-speed region also reaches depths of ∼200 km in the upper mantle. Measured travel-time residuals at station ODZ (Fig. 10D), located between Christchurch and Dunedin, show mostly late arrivals to the station, consistent with a low wave-speed region beneath the east coast of the South Island. Off the west coast under the Challenger plateau, similar low speeds also exist, though they are less reliable since they are on the edge of the array.

The tomogram shows a high-speed feature ∼7%–9% faster than reference speeds that is stronger at shallow depths compared to deeper in the mantle under the Southern Alps centered near 45°S, 169°E (Figs. 7 and 8 and cross sections D–D′ and E–E′ in Fig. 9; labeled “4” in Fig. 7A). The feature is ∼150 km in width and observed to no deeper than ∼150 km, and well defined by the observed travel times at station WKZ (Fig. 10G). Assuming an incidence angle of 30° and a mantle speed of 4.5 km/s, an S wave passing through a 150-km-thick layer of 9% fast material yields a delay time of ∼3.5 s, consistent with delays observed at WKZ. To ensure that this feature is not an artifact or the result of a poor crustal correction at station WKZ, we performed an inversion with travel-time measurements made at WKZ removed (Fig. 11B). Although the magnitudes of the high speeds under the Southern Alps decreased with the removal of arrival times at WKZ, the travel times at other stations still require high speeds in this region to depths of 100–150 km. Including travel-time measurements from this station in the inversion, but varying the crustal correction to account for error in the Pn time term (4.44 ± 0.13 s [Collins and Molnar, 2014], which maps to a crustal thickness of 42.8 ± 1.3 km beneath WKZ) also had little effect. To force high speeds to not exist under the Southern Alps, we had to vary the crustal thickness at WKZ, as well as at the neighboring stations EAZ, LBZ, JCZ, MLZ, and MSZ, by ∼10 km (Fig. 11C). A 10 km change in crustal thickness requires a change in the measured Pn time term of over 1 s, an error ∼10 times that reported in Collins and Molnar (2014) at those stations.

Cross section D–D′ in Figure 9 shows the high-speed structure under the Southern Alps dipping to the northwest. To investigate whether this apparent dip is the result of ray coverage, we inserted a synthetic vertical structure of 10% faster than reference speeds to 200 km depth under the Southern Alps (Fig. 12). Resolution of such a synthetic structure shows a dip toward the northwest as well as smearing of high speeds along this direction, similar to the dip and vertically smeared high speeds seen in the tomogram (cross section D–D′ in Fig. 9). Because of the vertical smearing, we investigated the minimum depth the model requires high speeds to reach under the Southern Alps. Plotting variance reduction as a function of model depth (Fig. 13A) suggests that variance reduction does not significantly change (i.e., change by >10%) when the model depth exceeds 285 km. Increases in variance reduction for model depths greater than ∼300 km are likely due to better fitting the high speeds under the northwestern South Island. We then performed a squeeze test (Figs. 13C–13E). First, we forced all the variations in the observed travel times into the top 285 km of the model. We then differenced the predicted travel times produced by the squeezed inversion from the observed travel-time measurements and inverted the differenced travel times for a model with a relaxed depth constraint. Squeezing the speed anomalies into the top 285 km (Fig. 13C) yields a structure similar to the inferred model (cross section D–D′ in Fig. 9), with high speeds (i.e., >5% of the reference speed) under the Southern Alps reaching 200–250 km. Relaxing the squeezed condition (Fig. 13E) reveals that these data do not require high-speed structure much deeper than ∼200 km under the Southern Alps.

Observed travel-time residuals at stations in the offshore region of the central west coast (e.g., NZ14 and NZ15) suggest that high speeds in the upper mantle do extend offshore. Almost all measured residuals at NZ14 (Fig. 10J) and NZ15 are more than ∼2.5 s early. Zietlow et al. (2014) found large SKS split times (>2 s) in this region, with orientations of fast quasi-S waves trending approximately parallel to N60°E, suggesting anisotropy may enhance the high speeds observed extending off the central west coast of the South Island. To test if anisotropy possibly enhances the southwestward extent of this high-speed body, we removed stations NZ14 and NZ15, and therefore rays traveling directly through this anisotropic region, from the inversion. With (tomograms labeled “Vs” in Fig. 8) or without (tomograms labeled “No NZ14 or NZ15” in Fig. 8) these stations, the nearly identical inferred structures show high speeds in the mantle under the northwestern South Island to depths below 350 km, as well as in the mantle offshore the central west coast.

We also investigated the effect that structure outside of the model region, such as subducted Pacific plate under the North Island, may have on the high-speed region in the northwest of the South Island. Rays passing through lithospheric slabs should arrive earlier at NZ14, NZ15, QRZ, NNZ, and other stations nearby. To do so, we removed travel-time residuals measured from earthquakes occurring near Tonga (back azimuths from 10°–30° in Fig. 3). Even after eliminating travel-time residuals from the Tonga events, high speeds remained in the upper mantle underneath the northwestern part of the South Island as well as offshore the central west coast, but the magnitude of the high-speed anomalies, especially below depths of ∼250 km, decreased (tomograms labeled “No Tonga” in Fig. 8). These two tests suggest that, although anisotropy and the passage of waves through high-speed regions outside the model space may contribute to the high-speed anomalies observed in the northwest, they can account for only a small fraction of the observed travel-time residuals.

Others have investigated the sediment thicknesses offshore the South Island of New Zealand (e.g., Ball et al., 2014) or produced regional or global data sets with relevant data (e.g., Divins, 2003; Whittaker et al., 2013). Such studies report some consistency in sediment thickness but also some variance. For instance, Ball et al. (2014) estimated a sediment thickness of 800 m at NZ16 using seafloor compliance, which compares well to the acoustic basement depth of 711 m reported by Divins (2003) and the 710 m reported by Whittaker et al. (2013) used in this study. At NZ06, however, Ball et al. (2014) predicted a sediment layer 1.7 km thick compared to 466 m and 450 m determined by Divins (2003) and Whittaker et al. (2013), respectively. Above we investigated the effect that varying S-wave speeds through sediment had on our model. To investigate the effect that varying sediment thickness had on our model with a fixed S-wave sediment speed, we performed inversions where the sediment layer was assumed to be zero in the travel-time measurement corrections, as well as one in which the sediment layer was assumed to be twice as thick as the values used from Whittaker et al. (2013) (“No Sediment” and “2x Sediment” in Fig. 8).

Varying the sediment thickness affects the top ∼300 km of the model. Assuming no sediment layer enhances and increases the width of high-speed features while diminishing slow-speed features. The opposite is true when assuming sediment thicknesses twice those given in Whittaker et al. (2013). Adjusting the sediment layer thickness also affects predicted onshore wave speeds. This is not surprising since the mean is removed from all the travel-time measurements after correcting for crust and sediment. While the magnitude and size of the recovered structures shift slightly with varying sediment thickness, none of the large-scale features disappear altogether, suggesting that errors in sediment thickness do not impact our model drastically.

All of the features discussed in the Vp tomogram of Zietlow et al. (2016) are evident in this new Vs tomogram, namely the high-speed structures in the northwestern and southwestern parts of the South Island, as well as the low-speed features along the eastern coast. The key difference between the Vp and Vs tomograms is the resolution of high-speed structure in the mantle under the Southern Alps in the Vs tomogram. A number of reasons could explain the ability to resolve these high speeds in the S-wave model and not the P-wave model, given the same data set. In a finite-frequency inversion, the radius of the first Fresnel zone generated by a P wave is wider than that generated by an S wave. For instance, from Equation 2 for Vp = 8.1 km/s versus Vs = 4.5 km/s (typical upper mantle speeds), the difference in the radius of the first Fresnel zone is ∼25 km at 100 km depth. This means that speeds in a P-wave tomogram generated from a finite-frequency inversion are distributed over a larger volume than in an S-wave tomogram and may appear more diffuse. Also, while slower traveling S waves could yield more diffractions off a structure, slower wave speed means more time to generate larger delays and thus make a structure more resolvable.

Regarding the main features of the Vs tomogram, the 150–200-km-deep ∼5%–8% low-speed features along the east coast near Christchurch and Dunedin (Figs. 7 and 8; labeled “3” in Fig. 7A) coincide with regions that experienced Miocene and younger surface volcanism (Godfrey et al., 2001, and references therein). The Akaroa and Lyttelton volcanoes near Christchurch last erupted ca. 6 Ma, and the Otago volcano near Dunedin last erupted ca. 10 Ma. Both regions also correspond to elevated heat flux, ∼70–74 mW m^–2 around Christchurch and ∼90–92 mW m^–2 around Dunedin (Godfrey et al., 2001). The elevated heat flux and recent surface volcanism imply a thinned lithosphere beneath the east coast of the South Island, suggesting that the observed low-speed regions in this area are anomalously warm. Other studies also propose a thinned lithosphere under the east coast, namely those based on geochemical investigations that invoke lithospheric removal to explain the composition and origin of the Cenozoic intraplate volcanism (Hoernle et al., 2006; Timm et al., 2009).

In the southwestern South Island, observed high speeds result from subduction of Australian lithosphere beneath the Pacific plate at the Puysegur Trench. The nearly vertical high-speed zone extends to depths of ∼200 km (Figs. 7 and 8; labeled “2” in Fig. 7A), consistent with the tomogram of Eberhart-Phillips and Reyners (2001), as well as with studies of seismicity that show a nearly vertical seismic zone (e.g., Anderson and Webb, 1994). This is also consistent with the P-wave tomogram of Zietlow et al. (2016) (tomograms labeled “Vp and Vs” in Fig. 8 and Fig. 14), though the region is not as well defined in width and depth in this study possibly because of differing ray coverage in this southwestern area compared to the P-wave study.

A significant finding from Zietlow et al. (2016) was the depth that high speeds, interpreted as marking subducted Pacific lithosphere, reached beneath the northwestern South Island to ∼400–450 km. In this S-wave study, we again see evidence for near-vertical high speeds to depths of at least 400 km in the mantle under the northwestern South Island (Figs. 7 and 8 and cross sections A–A′ and B–B′ in Fig. 9), even when ray paths passing through Pacific lithosphere beneath the North Island are removed. New Zealand’s permanent network, GeoNet (Petersen et al., 2011), has recorded some earthquakes as deep as 450 km under the northwestern region of the South Island, with earthquakes occurring ∼100 km deep to the southwest. This again supports the inference of subducted Pacific lithosphere to depths of ∼400–450 km beneath the northwestern portion of the South Island. Furthermore, the lateral extent of the high shear-wave speeds is similar to that found in a surface-wave study also using data from the MOANA experiment (Ball et al., 2016). In Ball et al. (2016), they imaged a zone of high shear-wave speeds (Vs > 4.6 km/s) from 50 to 150 km depth offshore the northwestern South Island and extending southwards.

The enhanced high speeds under stations CROE, NZ14, NZ15, and WVZ in the S-wave tomogram compared to the P-wave tomogram of Zietlow et al. (2016) (tomograms labeled “Vp and Vs” in Fig. 8; and Fig. 14) may be due to the inclusion of additional travel-time measurements on earthquakes not present in the P-wave study. For instance, SKS phases (incidence angles <15° in Fig. 10J) show travel-time residuals more than ∼2.5 s fast. Ray paths from an additional earthquake occurring near Tonga (back azimuth of ∼30° in Fig. 10J) with travel-time residuals more than 3 s fast may also contribute.

Pacific lithosphere to depths of ∼400 km is consistent with estimates of finite rotations that show 850–1000 km of Pacific lithosphere has been obliquely subducted beneath the northwestern South Island since 45 Ma (Sutherland, 1995). We estimate ∼500–600 km of gently dipping subducted slab lies between the Hikurangi Trench and the observed high speeds, and reach a depth of ∼100 km. For oblique convergence of 850–1000 km, the vertical extent of the subducted slab should reach 400–600 km depth, as seen here. Furthermore, as proposed by Zietlow et al. (2016), the depth that high speeds reach in the S-wave tomogram below the location of intermediate depth (subcrustal to 200 km) seismicity suggests that ∼200–300 km of slab is needed below the deepest earthquakes to generate enough weight to produce the intermediate depth seismicity. The depths that subducted slab reach under the northern South Island may also play a contributing factor in large crustal earthquakes in the region, such as the 2016 Kaikoura event. This crustal earthquake was complex, rupturing across multiple faults (e.g., Duputel and Rivera, 2017; Hollingsworth et al., 2017); however, our model presented here does not accurately resolve the shallow subsurface (i.e., <60 km), so we can only speculate.

Other studies (e.g., Stern et al. 2000; Kohler and Eberhart-Phillips, 2002; Reyners et al., 2017) report high-speed structure under the Southern Alps similar to the structure presented here (Figs. 7 and 8; cross sections D–D′ and E–E′ in Fig. 9; feature labeled “4” in Fig. 7A). The 7%–9% faster than reference speeds structure observed in this study appears no wider than ∼150 km to ∼150 km deep in the mantle. Again, while this high-speed region is less than the width of the first Fresnel zone, resolution tests do show that such a feature is resolvable with our data but will appear wider (Figs. 5 and 6).

The exact shape of this structure is difficult to resolve, but we can put some limitations on it. The squeeze test (Fig. 13) indicates that high-speed structure under the Southern Alps deeper than ∼200 km is not required by these travel-time data. Furthermore, travel-time residuals at station WKZ (Fig. 10G) are consistent with a high-speed structure under the Southern Alps to no deeper than ∼200 km. Measurements from earthquakes occurring to the northwest with incidence angles at the Moho of ∼30° show residuals of ∼3 s faster than reference speeds. Assuming an S-wave speed of 4.5 km/s, a 140-km-thick 9% fast structure produces about a 3 s fast S-wave travel-time residual. Likewise, a 180-km-thick 7% fast structure also produces a 3 s fast residual.

From synthetic tests, we found that a vertical, high-speed structure will appear dipping to the northwest due to our ray geometry (Fig. 12). Considering only the shallow central portion (i.e., <100 km deep) of at least 7% faster than reference speeds, the high-speed feature under the Southern Alps appears near vertical. Furthermore, if a shallowly northwestward dipping slab (i.e., dipping <20°) existed beneath the South Island, stations along and off the west coast would record early arrivals from earthquakes to the northwest. From such earthquakes (Fig. 4E and Figs. S2c and S2d [^{footnote 1}]), OBS stations deployed off the west coast record late arrivals. On land, stations such as JCZ and MSZ also record late arrivals for earthquakes to the northwest less than 60° away. Since stations such as JCZ and MSZ record late arrivals, this suggests high speeds in the central portion of the South Island do not reach farther west than these stations and further indicates that this high-speed region is near vertical. Thus, this high-speed feature is near vertical, less than 150 km wide, and no deeper than 200 km.

Reyners et al. (2017) image high speeds in this region to 105 km depth. Likewise, Kohler and Eberhart-Phillips (2002) show a high-speed structure 60–100 km wide to no deeper than 200 km, close to where we observe a similar structure (cross section D–D′ in Fig. 9). Stern et al. (2000) suggested an ∼100-km-wide structure reaching to depths of ∼170 km beneath the thickest crust, but not highest topography, of the Southern Alps from measurements of teleseismic P-wave delays. This high-speed feature discussed in Stern et al. (2000) is farther northeast than the one discussed in this study (closer in location to cross section C–C′ in Fig. 9 where the dominance of the high speeds associated with the subducting slab offshore may be swamping any signal from high speeds under the Alps). It is possible, though, that the high-speed structure seen here is a southward continuation of the one inferred by Stern et al. (2000). The speeds inferred from Stern et al. (2000) and this study are certainly consistent. Assuming a Vp/Vs value of 1.7 (Reyners et al., 2006; Eberhart-Phillips et al., 2014) and an anomalous speed recovery of 60%, the speed advance in Vs seen in this study maps to an ∼7%–9% advance in Vp, similar to that reported by Stern et al. (2000).

Stern et al. (2000) and Kohler and Eberhart-Phillips (2002) associated these high speeds with Cenozoic shortening across the South Island. With ∼70 km of shortening across the South Island (Cande and Stock, 2004) and with a lithospheric mantle 100 km thick, a high-speed body formed due to lithospheric thickening in a convergence zone should reach a depth of 150–200 km, as seen here. Alternatively, Reyners et al. (2017) suggest that high speeds in this region at depths greater than 65 km are associated with a subducted Hikurangi plateau that buttresses a lower-speed root associated with the Southern Alps. From this, it would follow that a strong and rigid lithosphere underlies much of the central South Island. In fact, they attribute the 2010 Mw 7.1 Darfield earthquake aftershock sequence to the existence of the Hikurangi plateau beneath the city of Christchurch, the top of which, they argue, is ∼10 km below the surface (Reyners et al., 2013).

Based on the Vs model presented here, we cannot distinguish whether the imaged high-speed structure in the central portion of the island is associated with the oblique convergence across the island or a subducted piece of a strong, rigid Hikurangi plateau. Our Vs model begins at 60 km depth, below where the bulk of the putative Hikurangi plateau hypothetically lies; however, both the Vp tomogram of Zietlow et al. (2016) and the Vs tomogram presented here show a clear region in the upper mantle with lower than reference speeds between the cities of Christchurch and Dunedin. These low speeds are evident at 60 km depth and possibly shallower. Given this, a region of low-speed material along the eastern coast of the South Island and extending ∼100 km inland seems inconsistent with a cold, rigid plateau underlying much of the central portion of the island.

We inverted teleseismic S-wave travel times at seismographs on and offshore the South Island of New Zealand to yield a 3D S-wave speed structure. Similar to the teleseismic P-wave tomogram of Zietlow et al. (2016), a low S-wave speed feature along the east coast in the mantle to depths of 150–200 km indicates thinned lithosphere underlying Miocene and younger volcanism. Subduction at the Puysegur trench is imaged to ∼200 km beneath the southwestern South Island. In the mantle beneath the northwestern South Island, high S-wave speeds >10% higher than average continue to depths of ∼400 km, and apparently mark Pacific lithosphere subducted obliquely westward since 45 Ma. Unlike in the Vp tomogram, the Vs tomogram shows a high-speed region of ∼7%–9% faster than average speeds no wider than ∼150 km and reaching to ∼150 km into the mantle below the central portion of the island. This structure is consistent with either Cenozoic shortening across the island or the subducted remains of the Hikurangi plateau. The low-speed features along the eastern coast of the island, however, are not consistent with a cold, rigid plateau lying just below the surface throughout the central South Island. A linear regression of Vp and Vs travel-time residuals yields a ratio of relative speed perturbations of ∼1.68, suggesting wave speed heterogeneity in the mantle under the South Island of New Zealand is due mostly to variations in temperature.

We thank Justin Ball, Carolin Boese, Cailey Condit, Danny Feucht, Jenny Nakai, Colin O’Rourke, Martha Savage, Vera Schulte-Pelkum, Steven Plescia, Joshua Stachnik, Laura Wallace, Stuart Wier, and especially Craig Jones and Peter Molnar for helpful discussions, Brandon Schmandt and Steven Roecker for sharing their tomography codes and advising us on how to use them, and John Collins for leading the MOANA OBS experiment. We thank Tim Stern for a helpful and thorough review of the manuscript and Associate Editor Raymond Russo. The ocean bottom seismic data used in this research were provided by instruments from the Ocean Bottom Seismograph Instrument Pool (www.obsip.org), which is funded by the National Science Foundation. The Incorporated Research Institutions for Seismology (IRIS) PASSCAL program provided three of the land broadband stations. We thank Martha Savage, Mark Henderson, Carolin Boese, and Sapthala Karaliyadda for deploying and maintaining these stations. The OBSIP and PASSCAL data are available for download from the IRIS Data Management Center. The New Zealand GeoNet project and its sponsors, the Earthquake Commission (EQC), GNS Science, and Land Information New Zealand (LINZ), provided additional data, which are also available for download. The figures were made with the Generic Mapping Tools (GMT) software by Wessel and Smith (1998). This study was supported by the National Science Foundation under grant EAR-0409835, and in part by a Cooperative Institute for Research in Environmental Sciences (CIRES) Graduate Research Fellowship Award and by a gift from the Crafoord Foundation.