Abstract

Insights into plate-boundary deformation are gained from the shear-wave splitting of local S phases that originate within the lithosphere of the South Island, New Zealand. Analysis of the splitting parameters from land stations reveals changes in both delay times and fast azimuths with earthquake depth, earthquake-station back-azimuth and initial polarization azimuth, suggesting both laterally and depth varying anisotropy. When the average results are examined as a whole via tomographic inversion and spatial averaging, consistent patterns in delay times and fast azimuths exist. Spatially averaged fast azimuths reveal a localized high-strain zone in the southern central region of the South Island. Based on fast azimuths observed above 100 km depth, we suggest that the plate-boundary subparallel anisotropy that is produced by pervasive shear is mainly distributed within a zone extending ∼130 km SE of the Alpine fault in southern South Island and is widely distributed in northern South Island. Average station delay times of ∼0.1–0.4 s compared to 1.7 s SKS delay times from previous studies in South Island further suggests a deeper seated anisotropic zone in some areas, but in other regions frequency-dependent anisotropy could explain the observed splitting.

INTRODUCTION

Anisotropy beneath plate-boundary zones provides indirect insights into the tectonic deformation associated with plate interactions. Among geophysical techniques, shear wave splitting (S-wave splitting) has been a well-established and commonly used technique in exploring seismic anisotropy (Silver & Chan 1991; Teanby et al.2004b; Wustefeld et al.2008). Splitting tomography is another aspect that has been utilized to resolve the heterogeneous anisotropy structure from initial splitting parameters and to provide constraints on the distribution of seismic anisotropy (Abt & Fischer 2008; Johnson et al.2011). Usage of such techniques is important in local and regional splitting studies, because it considers the changes in splitting due to various factors (e.g. back azimuth, depth, orientation of crystallographic axes) that could result in complex splitting. However, such methods are highly dependent on the station coverage and the earthquake distribution.

The crustal contributions of anisotropy are usually ignored by splitting analysis that investigates mantle anisotropy using SKS phases (Savage 1999). Yet crustal contributions might have significant effects on SKS measurements and result in complex patterns and small-scale lateral variations (Rumpker & Silver 1998; Mattatall & Fouch 2007; Long & Silver 2009). Thus, S-wave splitting of local and regional S-phases is vital for understanding the crustal and subcrustal contribution to total anisotropy in a given region. S-wave splitting using local events observed along the ray paths in the lithosphere is attributed to different sources of anisotropy that are developed from both shape preferred orientation and lattice preferred orientation. Seismic anisotropy at crustal depths is assumed to result from either stress-induced anisotropy or structural anisotropy (e.g. Boness & Zoback 2006). In the upper crust, fluid-filled microcracks or fractures and structural anisotropy (e.g. crustal faults, rock foliations or beddings) may affect S-wave splitting whereas in the lower crust and the upper mantle, splitting is mainly influenced by intrinsic anisotropy resulting from rock fabrics and preferred orientation of rock-forming minerals such as olivine in the upper mantle and mica and amphibole in the lower crust (Mainprice & Nicolas 1989; Silver 1996; Crampin et al.2004). As suggested by Wilson et al. (2004), subhorizontal foliations with lineations (e.g. schist fabric) developed at lower crustal depths due to pervasive shear in strike-slip fault zones could also act as source of anisotropy.

The interpretation of local splitting is challenging because of the different factors and sources of anisotropy (e.g. tectonic processes) that control splitting. Lattice preferred orientation is mainly a result of finite strain (Nicolas & Christensen 1987). In the lower crust, mineral alignment is predicted to be parallel to the flow or shear direction and in the upper mantle, the [100] axis of olivine is fast and is parallel to the extension (maximum principal stretch) direction or to the flow direction if significant strain has developed in a dry regime (Nicolas & Christensen 1987; Silver 1996; Savage 1999). At shallow crustal depths, differential stress results in preferential closure of microcracks perpendicular to the maximum stress direction creating an anisotropic medium (Crampin et al.2004). The fast azimuth of a split S-wave in such a stress-aligned anisotropic medium is parallel to the cracks that remain open, hence, parallel to maximum principal stress direction (Nur & Simmons 1969; Boness & Zoback 2006).

Plate interactions at the plate-boundary zone in South Island, New Zealand are expected to result in deformation-induced seismic anisotropy. Westward movement of the Pacific Plate with respect to the Australian plate results in an obliquely convergent boundary across New Zealand (Fig. 1). The main surface manifestation of this plate-boundary, the Alpine fault, extends for about 460 km, serving as a transform fault that connects Puysegur and Hikurangi subduction zones with opposite polarities (Fig. 1). The structural evolution of South Island is marked by different tectonic episodes. Some of the key processes include rifting and the breakup of Gondwana in the late Cretaceous, evolution of the two subduction systems and the plate-boundary, the transition of the southern subduction system (Puysegur subduction) into the central Alpine fault zone during the Eocene to Quaternary periods and southward migration of the Hikurangi subduction system increasing the transpressional character of the central collisional (Alpine fault) system (Sutherland et al.2000; Cox & Sutherland 2007; Furlong 2007). Seismic anisotropy could remain from tectonic deformation during these episodes.

Figure 1.

(a) Map of the tectonic setting of the South Island, New Zealand and the main structural features associated with the Australian-Pacific plate-boundary. Approximate location of the schist belt with high-grade metamorphism is marked by the dark-grey shaded region. The black arrow denotes the relative plate motion direction (∼38 mm yr−1) between the Australian and Pacific Plates (Walcott 1998). (b) Schematic diagrams (no vertical exaggeration) showing the proposed models of deformation for central South Island (dashed-line box in the inset of Figs 1a). (i) Pervasive shear or thin viscous sheet model proposed by Molnar et al. (1999). (ii) Localized deformation model.

Figure 1.

(a) Map of the tectonic setting of the South Island, New Zealand and the main structural features associated with the Australian-Pacific plate-boundary. Approximate location of the schist belt with high-grade metamorphism is marked by the dark-grey shaded region. The black arrow denotes the relative plate motion direction (∼38 mm yr−1) between the Australian and Pacific Plates (Walcott 1998). (b) Schematic diagrams (no vertical exaggeration) showing the proposed models of deformation for central South Island (dashed-line box in the inset of Figs 1a). (i) Pervasive shear or thin viscous sheet model proposed by Molnar et al. (1999). (ii) Localized deformation model.

Present deformation along the plate-boundary can also influence anisotropy. In the northern part of South Island, the relative plate motion (∼38 mm yr−1) between the Australian and Pacific Plates is partitioned into a plate-boundary/Alpine fault parallel component (37 mm yr−1 strike-slip movement) and a plate convergent component (11 mm yr−1) orthogonal to the Alpine fault, resulting in significant transpression (Norris et al.1990; Walcott 1998). Northern and north-eastern parts of this region are dominated by crustal faults (Fig. 1) that indicate transition from the Hikurangi subduction zone to the central collisional zone. Continental collision in the central region is marked by structural features such as the Southern Alps mountain range, the Alpine fault and a crustal root (∼37–42 km) and the thickened mantle lithosphere (∼150 km) beneath the Southern Alps (Scherwath et al.2002; Okaya et al.2007). In southern South Island, the Australian plate is being subducted obliquely underneath Fiordland, creating a relatively young subduction zone extending up to maximum depths of ∼120 km (Davey & Smith 1983).

There have been two opposing hypotheses to explain the deformation in the central South Island: the pervasive shear model and the narrow/localized deformation model (Fig. 1b). Structural and geodetic data suggest the partition of total plate displacement to 480 km of right lateral slip on the Alpine fault and ∼200 km wide distributed deformation to the east of the fault (Norris et al.1990). Based on the consistent patterns of fast azimuths observed from SKS studies in South Island, Klosko et al. (1999) suggested the existence of distributed deformation (∼200 km) within 40–400 km depth. Molnar et al. (1999) and Stern et al. (2000) argue that this distributed deformation zone is restricted to the first 100 km of the upper mantle and extends in a ∼100–150-km-wide region either side of the Alpine fault. They proposed a thin viscous sheet model (Molnar et al.1999; Moore et al.2002; Wilson et al.2004) that agrees with the pervasive shear model (Bourne et al.1998) (Fig. 1b-i) to explain such deformation behaviour in the presence of weak lower crust. Conversely, the localized deformation model (Fig. 1b-ii) assumes the presence of a weak upper mantle instead of a weak lower crust (Meade & Hager 2005) and suggests that the strain is localized along the faults within the crustal and upper-mantle depths and is more distributed in the asthenosphere and/or deep mantle due to asthenospheric flow or shear (Avouac & Tapponnier 1993; Devs et al.2011).

Here we focus on splitting of local S waves that originate within the lithosphere of the South Island to determine how anisotropy is distributed laterally and vertically in and around the plate-boundary zone across the South Island. We consider the possible sources of anisotropy and the possible mechanisms/processes associated with the formation of anisotropy. In comparison to previous local splitting studies (Section 4.4) in the northern most region of the South Island (Audoine et al.2000; Balfour et al.2005), we use a larger data set, covering the entire South Island and the complete depth range available, and apply a 2-D delay-time tomography/spatial averaging technique to resolve anisotropy. We compare our results with the available SKS splitting (Klosko et al.1999; Savage et al.2007a), finite strain calculations from GPS data and the structural geology data (Little & Mortimer 2001; Little et al.2002) to examine how the upper-lithospheric deformation varies from the asthenospheric deformation and to test the deformation models proposed for South Island.

DATA AND METHODOLOGY

We use both local (epicentral distance ≤∼100 km) and regional (100 km > epicentral distance >250 km) earthquake data recorded on 26 GeoNet permanent three component seismographs for the period of 2004–2010 and an 11 month deployment of four temporary broad band stations in the northern part of South Island to measure S-wave splitting parameters (Fig. 2). These four temporary PASSCAL seismographs were deployed by Victoria University of Wellington in collaboration with the University of Colorado to enhance the seismic coverage. The complete region has been divided into three geographic units: northern, central and southern South Island. Events in each block are collected from the stations in that particular block (Fig. 2b). Stations that are close to the borders (WKZ, MSZ, ODZ, CASS and LTZ) of a block used the events recorded from either side. This criterion is used to control the epicentral distances. To obtain good signal-to-noise ratio (SNR), we use a threshold magnitude of 3.5 in northern and southern South Island. Due to the lack of seismicity and deep events in central South Island, the threshold magnitude of 3.0 is chosen to obtain more events. We utilize an automatic S-wave splitting (Mfast) technique (Savage et al.2010; Wessel 2010), which was developed on the basis of the eigenvalue minimization method (Silver & Chan 1991) and the cluster analysis method (Teanby et al.2004b).

Figure 2.

(a and b) Maps showing the distribution of earthquakes that yield A and B quality splitting measurements in South Island (SI). The earthquakes are colour coded according to earthquake depth. (a) Events recorded from all possible incidence angles and (b) events that are within the incidence angle (ic) of 35°. Earthquakes within three blocks (northern SI, central SI and southern SI) are analyzed only using the station within that block, but the stations at the block boundaries (MSZ, WKZ, ODZ, CASS and LTZ) used the events recorded on both sides. Squares and triangles display the temporary stations and permanent stations (GeoNet network), respectively.

Figure 2.

(a and b) Maps showing the distribution of earthquakes that yield A and B quality splitting measurements in South Island (SI). The earthquakes are colour coded according to earthquake depth. (a) Events recorded from all possible incidence angles and (b) events that are within the incidence angle (ic) of 35°. Earthquakes within three blocks (northern SI, central SI and southern SI) are analyzed only using the station within that block, but the stations at the block boundaries (MSZ, WKZ, ODZ, CASS and LTZ) used the events recorded on both sides. Squares and triangles display the temporary stations and permanent stations (GeoNet network), respectively.

Splitting of shear waves into two quasi S-phases upon propagating through an anisotropic medium and the directional dependence of velocities of these split phases are the principal phenomena associated with the S-wave splitting method. The delay time (δt) between fast and slow split components and the polarization azimuth of the fast split phase (ϕ) are the two main parameters that are determined during the splitting analysis. These parameters provide important information about the extent and degree of anisotropy and the orientation of the anisotropic medium and the past and/or present deformational processes (Silver & Chan 1991; Crampin 1994). In an isotropic medium, the focal mechanism controls the S-wave particle motion, which is linear (e.g. Aki & Richards 2002). On entering an anisotropic medium, S waves obtain an elliptical particle motion (Fig. 3). Silver & Chan (1991) introduce an inverse splitting operator that removes the splitting from the S-waves to find the linearized particle motion. A grid search over all the possible pairs of fast azimuths (ranging from −90° to 90° with 1° increments) and delay times (ranging from 0 to 1 s at 0.01 s increments) is carried out. The values of ϕ and δt are those that give the most linear incoming particle motion, as measured by minimizing the smaller eigenvalue (λ2) of the corrected 2-D time-domain covariance matrix between the two horizontal components in the reference coordinate system (Fig. 3). Once the smaller eigenvalue is minimized, the corrected particle motion is linear and the eigenvector of the corresponding larger eigenvalue provides an approximate azimuth of the initial polarization or incoming polarization (ϕp) of the unsplit S-wave (Vidale 1986). This method uses the F-test to constrain the uncertainties of the splitting parameters by calculating the 95 per cent confidence limits of the ϕ and δt that correspond to the optimum solutions from the grid search analysis.

Figure 3a.

Outputs from the S-wave splitting analysis using Mfast (Savage et al.2010). Top-left diagram shows the manual S pick (black line) on three components (e, n and z) of the seismograph. Dashed lines indicate the minimum starting time and maximum end time of the different S-wave analysis windows. Top-right figure display the rotated horizontal components before correction for splitting (p and p) and after correction for splitting (corrected p and corrected p). Grey shaded regions in top figures indicate the selected S-wave analysis window. Dashed lines indicate minimum and maximum starting and ending windows. Two middle graphs show the splitting parameters (best is the blue cross) obtained from the Teanby et al. (2004b) cluster analysis. Bottom-left indicates the particle motion plots before and after correction for S-wave splitting. Here the grey region is the region not used in the analysis. Bottom-right plot is the contour diagram indicating the energy surfaces. The solution with the smallest second eigenvalue (λ2) is marked by blue cross and the 95 per cent confidence limit of this solution is denoted by the dark black line.

Figure 3a.

Outputs from the S-wave splitting analysis using Mfast (Savage et al.2010). Top-left diagram shows the manual S pick (black line) on three components (e, n and z) of the seismograph. Dashed lines indicate the minimum starting time and maximum end time of the different S-wave analysis windows. Top-right figure display the rotated horizontal components before correction for splitting (p and p) and after correction for splitting (corrected p and corrected p). Grey shaded regions in top figures indicate the selected S-wave analysis window. Dashed lines indicate minimum and maximum starting and ending windows. Two middle graphs show the splitting parameters (best is the blue cross) obtained from the Teanby et al. (2004b) cluster analysis. Bottom-left indicates the particle motion plots before and after correction for S-wave splitting. Here the grey region is the region not used in the analysis. Bottom-right plot is the contour diagram indicating the energy surfaces. The solution with the smallest second eigenvalue (λ2) is marked by blue cross and the 95 per cent confidence limit of this solution is denoted by the dark black line.

Figure 3b.

Same as figure caption in Fig. 3(a), but with an A grade measurement with higher δt (see Fig. 8).

Figure 3b.

Same as figure caption in Fig. 3(a), but with an A grade measurement with higher δt (see Fig. 8).

Figure 3c.

Same as figure caption in Fig. 3(a), but showing a B grade splitting measurement.

Figure 3c.

Same as figure caption in Fig. 3(a), but showing a B grade splitting measurement.

Before processing data with the automatic method, all the events are visually inspected and S arrivals are hand picked. The automatic method involves searching over a given set of 14 bandpass filters, which are two-pole Butterworth filters with frequency bandwidths within the range of 0.4 to 10 Hz, and calculates SNR to find filters that give the highest SNR-bandwidth product (≥3) for each earthquake-station pair. Then the signals obtained from the selected filters are used to estimate the splitting parameters. The dominant period of each S-phase is determined using the frequency that corresponds to the highest amplitude of the spectrum calculated from the 3 s window directly after the S-pick. Cluster analysis avoids uncertainties in the splitting parameters caused by the selection of S-wave analysis windows (Teanby et al.2004b). This method selects the analysis window that gives the most stable splitting parameters with small errors by performing splitting measurements through a given set of analysis windows (Fig. 3). We chose the splitting measurements that agree with the A and/or B grading criteria (here after AB grade) as good quality splitting measurements (Savage et al.2010). In Mfast, measurement attains AB grading if the measurement's SNR > 3, δt < 0.8*tlagmax (tlagmax is the maximum δt for grid search), maximum error in ϕ < 25°, if it is not a null measurement, and if the measurement graded as a A or B cluster in the cluster grading (see Appendix A1). Note that the A grade criteria [(i) SNR > 4, (ii) δt < 0.8*tlagmax, (iii) not null and dϕ < 10°] are a subset of AB grade criteria, which use more narrow limits for SNR and dϕ compared to AB grade criteria. Null measurements are those for which ϕp is within 20° of parallel or perpendicular to the fast azimuth, which can lead to incorrect δt. The splitting measurements from each event that meet the above quality scheme are again checked considering their uncertainties and errors for different filters. This final quality criterion extracts the measurements that give a similar answer, i.e. if the time difference > tlagmax/8 or if the angular distance > π/8, in each filter (refers to ‘best filter measurements’ in later sections). From the ∼1600 AB quality splitting measurements (Savage et al.2010) we selected the ∼850 measurements that are within the incidence angle of 35° (see Section 2.1) for further analysis (Fig. 2b). The number of good splitting measurements made at each station varies depending on the seismicity distribution (Fig. 2) and the quality criteria used in the analysis.

Limitations and statistical analysis of the measurements

The polarization of an S-wave incident at the free surface is governed by many factors such as the angle of incidence, frequency of the incident wave, the source depth and the layering of the crust (Booth & Crampin 1985). To minimize the effect of incidence angle on S-wave polarization, splitting analysis is usually carried out inside the shear-wave window (i.e. a few degrees less than critical incidence angle) (Nuttli 1961; Evans 1984). Wave fronts with incident angles beyond the critical angle (ic) are likely to display complex particle motion due to phase conversions at the free surface. This will cause particle motion to be non-linear even in the absence of anisotropy (Nuttli 1961) and make the interpretation of S-wave polarization complicated. To accurately calculate the incidence angle, we create a velocity model for South Island by averaging the velocities obtained from velocity-structure studies in South Island (Eberhart-Phillips & Reyners 2001; Eberhart-Phillips & Bannister 2010). This average velocity model has an upper layer with Vs (S-wave velocity) and Vp (P-wave velocity) values of 2.6 and 4.4 km s−1, respectively. The low velocity surface layer allows earthquake waves from a broad epicentral range to approach the surface within the S-wave window (Booth & Crampin 1985). In South Island, most of the seismicity is limited to the plate-boundary zone and this narrow distribution significantly reduces the number of events that can be used to measure the crustal anisotropy far away from the plate-boundary. Therefore, some of the stations in the eastern margin of the central and southern South Island (e.g. OPZ, ODZ and SYZ) are only used for deep events (depth > ∼60 km) (Fig. 4b).

Figure 4.

Straight-line ray paths showing the coverage in the subduction zones in northern South Island (a) and southern South Island (b). Right side shows the cross-sections (depth profiles) along AA and BB lines that are marked on the left side maps. Events and ray paths corresponding to each station are colour coded using the same colour. The approximate slab positions are coloured are in grey.

Figure 4.

Straight-line ray paths showing the coverage in the subduction zones in northern South Island (a) and southern South Island (b). Right side shows the cross-sections (depth profiles) along AA and BB lines that are marked on the left side maps. Events and ray paths corresponding to each station are colour coded using the same colour. The approximate slab positions are coloured are in grey.

In addition to calculating uncertainties and to grading individual splitting measurements, statistical parameters [e.g. mean, standard deviation, standard error (SE)] of the single station measurements are estimated (Table 1). We use directional statistics discussed in Mardia (1972) and Gerst & Savage (2004) for ϕ. Directional statistics are necessary to provide quantitative measures of the ϕ and to estimate errors and accuracy of the mean ϕ (Roman et al.2011). Because of the 180° ambiguity of angular data, directional statistics are different from normal statistics. Mean ϕ, the length of the mean resultant vector of ϕ (R), circular standard deviation (CSD) and SE are the statistical parameters that we estimate. R is a quantitative estimate of circular spread that is calculated from the mean angular azimuth (ϕs). SE of ϕs is calculated using R and the concentration parameter (κ), which is a parameter of a Von Mises distribution ($$VM_{(\phi _s,\kappa )}$$), and SE indicates the accuracy of the ϕs. $$VM_{(\phi _s,\kappa )}$$ is a continuous distribution on the range of 0 − 360° and it is a probability density function that varies with ϕs and κ. Thus, $$VM_{(\phi _s,\kappa )}$$ is considered to be a circular equivalent of the normal distribution that measures the circular spread (Mardia 1972; Philipp 2009). Tables 1(a) and (b) summarize the above statistical parameters that are calculated using splitting measurements made at each station and Section 3.1 explains these results.We have also conducted a detailed analysis of splitting measurements considering the possible relationships between splitting parameters and other variables such as earthquake depth, earthquake-station back azimuth, average maximum frequency (or an average period) in the shear-wave analysis window and ϕp to provide constraints on both lateral distribution and depth extent of anisotropy (see Section 4.1).

Table 1.

Summary statistics of the splitting parameters from: (a) shallow (≤100 km) events and (b) deep (>100 km) events.

(a) 
Station No ϕsa SEb of ϕs CSDc of ϕ Rd Mean δt SE of δt SDe of δt 
DSZ 49 60.526 11.701 48.148 0.244 0.172 0.016 0.115 
MQZ 20 −26.623 14.977 44.727 0.296 0.180 0.019 0.086 
KHZ 68 −49.314 16.878 56.321 0.145 0.196 0.015 0.127 
LTZ 79 −72.066 5.310 38.273 0.410 0.188 0.017 0.148 
NNZ 111 −69.658 8.798 50.158 0.216 0.197 0.013 0.138 
QRZ 66 38.922 8.153 44.532 0.299 0.184 0.011 0.086 
FOZ 15 −84.337 7.642 28.637 0.607 0.103 0.020 0.079 
THZ 58 11.330 17.411 55.619 0.152 0.285 0.019 0.146 
CRLZ 43 −34.399 3.838 25.109 0.681 0.147 0.010 0.066 
WKZ 94 32.646 7.506 46.163 0.273 0.302 0.016 0.158 
JCZ 61 −74.061 14.429 53.199 0.178 0.203 0.015 0.116 
LBZ 79 38.882 9.187 48.098 0.244 0.206 0.015 0.129 
DCZ 139 −9.234 3.246 34.081 0.493 0.159 0.008 0.089 
MSZ 99 12.399 13.200 55.479 0.153 0.218 0.017 0.171 
EAZ 95 8.267 8.333 48.010 0.246 0.222 0.012 0.114 
MLZ 120 41.762 7.044 47.155 0.258 0.132 0.007 0.076 
TUZ 86 −9.634 4.820 37.210 0.430 0.203 0.011 0.099 
WHZ 83 −26.898 6.790 43.319 0.319 0.210 0.013 0.121 
ODZ 38 −25.905 4.077 25.077 0.682 0.269 0.029 0.177 
OXZ 15 −32.126 9.739 33.780 0.499 0.232 0.041 0.160 
APZ 38 −83.567 6.099 33.713 0.500 0.263 0.020 0.124 
OPZ 13 15.006 19.602 45.660 0.281 0.182 0.022 0.080 
PYZ 24 22.666 16.154 47.578 0.252 0.229 0.025 0.121 
SYZ 26 −9.093 11.187 41.845 0.344 0.108 0.010 0.053 
CASS 66.843 10.572 25.730 0.668 0.139 0.027 0.066 
CROE 51 −57.412 8.198 42.323 0.336 0.153 0.012 0.086 
FREW 10 −75.721 22.599 45.851 0.278 0.209 0.046 0.145 
KELY 107 8.748 6.321 44.303 0.302 0.176 0.011 0.114 
RPZ 29 25.981 11.223 42.899 0.326 0.177 0.016 0.088 
DSZ 172 −89.818 4.740 43.407 0.317 0.190 0.008 0.101 
KHZ 131 42.588 17.024 60.993 0.104 0.185 0.009 0.101 
LTZ 91 −58.693 12.466 54.010 0.169 0.190 0.014 0.129 
NNZ 137 22.979 20.417 63.669 0.085 0.181 0.009 0.105 
QRZ 131 −89.716 5.395 43.286 0.319 0.239 0.010 0.118 
THZ 103 7.699 5.271 40.653 0.365 0.185 0.012 0.122 
WKZ 36 −8.316 12.335 46.449 0.269 0.259 0.027 0.161 
JCZ 43 −60.998 2.761 18.326 0.815 0.219 0.009 0.059 
LBZ 124 38.022 1.786 20.145 0.781 0.165 0.006 0.064 
DCZ 91 −12.607 4.413 36.018 0.454 0.162 0.010 0.096 
MSZ 84 −67.490 11.236 51.816 0.195 0.257 0.013 0.115 
EAZ 53 −11.328 22.251 58.472 0.125 0.259 0.024 0.173 
MLZ 84 59.590 9.536 49.211 0.229 0.193 0.013 0.119 
TUZ 74 −41.503 5.457 38.171 0.412 0.206 0.018 0.154 
WHZ 79 −9.397 11.475 51.666 0.197 0.178 0.015 0.129 
ODZ 36 −20.020 8.517 39.788 0.381 0.228 0.026 0.156 
OXZ 20 −34.412 20.972 50.346 0.213 0.308 0.027 0.120 
APZ 23 −86.429 7.795 33.596 0.503 0.191 0.034 0.164 
OPZ 86.165 64.184 60.834 0.105 0.182 0.028 0.084 
PYZ 20 −19.624 11.095 39.224 0.392 0.199 0.030 0.136 
SYZ 52 32.961 11.443 48.271 0.242 0.129 0.010 0.073 
CASS 49.081 35.524 43.222 0.320 0.230 0.084 0.146 
CROE 48 −24.622 44.491 66.860 0.066 0.208 0.020 0.135 
FREW 23 21.890 18.091 49.093 0.230 0.247 0.037 0.178 
KELY 83 13.813 11.706 52.357 0.188 0.209 0.013 0.117 
RPZ 76 −19.418 5.329 37.968 0.416 0.245 0.017 0.147 
(a) 
Station No ϕsa SEb of ϕs CSDc of ϕ Rd Mean δt SE of δt SDe of δt 
DSZ 49 60.526 11.701 48.148 0.244 0.172 0.016 0.115 
MQZ 20 −26.623 14.977 44.727 0.296 0.180 0.019 0.086 
KHZ 68 −49.314 16.878 56.321 0.145 0.196 0.015 0.127 
LTZ 79 −72.066 5.310 38.273 0.410 0.188 0.017 0.148 
NNZ 111 −69.658 8.798 50.158 0.216 0.197 0.013 0.138 
QRZ 66 38.922 8.153 44.532 0.299 0.184 0.011 0.086 
FOZ 15 −84.337 7.642 28.637 0.607 0.103 0.020 0.079 
THZ 58 11.330 17.411 55.619 0.152 0.285 0.019 0.146 
CRLZ 43 −34.399 3.838 25.109 0.681 0.147 0.010 0.066 
WKZ 94 32.646 7.506 46.163 0.273 0.302 0.016 0.158 
JCZ 61 −74.061 14.429 53.199 0.178 0.203 0.015 0.116 
LBZ 79 38.882 9.187 48.098 0.244 0.206 0.015 0.129 
DCZ 139 −9.234 3.246 34.081 0.493 0.159 0.008 0.089 
MSZ 99 12.399 13.200 55.479 0.153 0.218 0.017 0.171 
EAZ 95 8.267 8.333 48.010 0.246 0.222 0.012 0.114 
MLZ 120 41.762 7.044 47.155 0.258 0.132 0.007 0.076 
TUZ 86 −9.634 4.820 37.210 0.430 0.203 0.011 0.099 
WHZ 83 −26.898 6.790 43.319 0.319 0.210 0.013 0.121 
ODZ 38 −25.905 4.077 25.077 0.682 0.269 0.029 0.177 
OXZ 15 −32.126 9.739 33.780 0.499 0.232 0.041 0.160 
APZ 38 −83.567 6.099 33.713 0.500 0.263 0.020 0.124 
OPZ 13 15.006 19.602 45.660 0.281 0.182 0.022 0.080 
PYZ 24 22.666 16.154 47.578 0.252 0.229 0.025 0.121 
SYZ 26 −9.093 11.187 41.845 0.344 0.108 0.010 0.053 
CASS 66.843 10.572 25.730 0.668 0.139 0.027 0.066 
CROE 51 −57.412 8.198 42.323 0.336 0.153 0.012 0.086 
FREW 10 −75.721 22.599 45.851 0.278 0.209 0.046 0.145 
KELY 107 8.748 6.321 44.303 0.302 0.176 0.011 0.114 
RPZ 29 25.981 11.223 42.899 0.326 0.177 0.016 0.088 
DSZ 172 −89.818 4.740 43.407 0.317 0.190 0.008 0.101 
KHZ 131 42.588 17.024 60.993 0.104 0.185 0.009 0.101 
LTZ 91 −58.693 12.466 54.010 0.169 0.190 0.014 0.129 
NNZ 137 22.979 20.417 63.669 0.085 0.181 0.009 0.105 
QRZ 131 −89.716 5.395 43.286 0.319 0.239 0.010 0.118 
THZ 103 7.699 5.271 40.653 0.365 0.185 0.012 0.122 
WKZ 36 −8.316 12.335 46.449 0.269 0.259 0.027 0.161 
JCZ 43 −60.998 2.761 18.326 0.815 0.219 0.009 0.059 
LBZ 124 38.022 1.786 20.145 0.781 0.165 0.006 0.064 
DCZ 91 −12.607 4.413 36.018 0.454 0.162 0.010 0.096 
MSZ 84 −67.490 11.236 51.816 0.195 0.257 0.013 0.115 
EAZ 53 −11.328 22.251 58.472 0.125 0.259 0.024 0.173 
MLZ 84 59.590 9.536 49.211 0.229 0.193 0.013 0.119 
TUZ 74 −41.503 5.457 38.171 0.412 0.206 0.018 0.154 
WHZ 79 −9.397 11.475 51.666 0.197 0.178 0.015 0.129 
ODZ 36 −20.020 8.517 39.788 0.381 0.228 0.026 0.156 
OXZ 20 −34.412 20.972 50.346 0.213 0.308 0.027 0.120 
APZ 23 −86.429 7.795 33.596 0.503 0.191 0.034 0.164 
OPZ 86.165 64.184 60.834 0.105 0.182 0.028 0.084 
PYZ 20 −19.624 11.095 39.224 0.392 0.199 0.030 0.136 
SYZ 52 32.961 11.443 48.271 0.242 0.129 0.010 0.073 
CASS 49.081 35.524 43.222 0.320 0.230 0.084 0.146 
CROE 48 −24.622 44.491 66.860 0.066 0.208 0.020 0.135 
FREW 23 21.890 18.091 49.093 0.230 0.247 0.037 0.178 
KELY 83 13.813 11.706 52.357 0.188 0.209 0.013 0.117 
RPZ 76 −19.418 5.329 37.968 0.416 0.245 0.017 0.147 

a Mean fast azimuth,

b Standard error,

c Circular standard deviation,

d the length of the mean resultant vector ϕ,

e Standard deviation.

Resolving the spatial heterogeneity of splitting measurements

Heterogeneous anisotropic structure observed is further evaluated using 2-D delay-time tomography and a spatial averaging technique (TESSA) introduced by Johnson et al. (2011). This allows us to obtain approximate constraints on regions with low and high anisotropy. A δt from a single measurement is a quantitative estimation of the anisotropic structure along the ray path from source to the station and the δt has been considered to be proportional to path length (L) of the ray through the anisotropic medium (Silver 1996). TESSA is based on a simplified and linear relationship between δt and (L). It assumes that the δtr observed along the ray path is the summation of δt calculated for each grid block (δtb) with variable strength of anisotropy (sb) in 2-D space. First, straight-line ray tracing between station and events determines the lengths of ray segments (Lrb) that lie on each grid block (Fig. A1). Then the total anisotropy measured at the station for a given event is calculated by,  

(1)
\begin{equation} {\delta }t_r = \displaystyle \sum \limits _{b=1}^n ({s_b \times L_{rb}} ). \end{equation}
The inversion part of TESSA determines the sb for resolved blocks and is carried out using the MATLAB medium-scale optimization inversion function that finds the linear least-square solution of the problem. This solution is then used iteratively to find the model solution for the strength of anisotropy of the whole ray (sr) and this process is subjected to boundary conditions. As explained in Johnson et al. (2011), boundary constraints are applied in such a way that the minimum sr will always be greater than 0 s km−1 and the maximum sr will not exceed δtr(max)/Lb, where δtr(max) is the maximum δt observed along a ray path and Lb is the width of the grid block. Tomography results discussed in Section 4.3 are obtained using regular gridding with 25 km square blocks. Minimum and maximum number of rays per block are set up to be 6 and 50, respectively, to obtain reliable results. Tomographic analysis only uses the AB grade measurements that are within the incidence angle of 35° and that satisfy the best filter criteria. This enables us to avoid ambiguous splitting measurements as well as statistical uncertainties that could affect the final solution of the tomography and spatial averages.

We use the tomographic weighting scheme presented in TESSA to determine average ϕ in grid blocks. Tomographic weighting is designed to take into account variations due to heterogeneous anisotropic structure. The tomographic weighting function (ωrb) that is applied in this averaging weights the ϕ considering the anisotropy strength profile determined for each grid block (sb) in the δt tomography, so that more weighting is given to the regions with higher strength of anisotropy. ωrb is defined as the strength profile of the each grid block, normalized by δtr (Johnson et al.2011),  

(2)
\begin{equation} \displaystyle \omega _{rb} = {s_b/\delta {t_r}}. \end{equation}

Then, the weighted mean ϕ of the grid block b with n number of rays is given by azimuth ($$\bar{\phi _b}$$),

 

(3)
\begin{equation} \bar{\phi _b} = \frac{1}{2} \arctan \left( \frac{\displaystyle \sum \limits _{r=1}^n \cos (2\phi _r) \times \omega _{rb}}{\displaystyle \sum \limits _{r=1}^n \sin (2\phi _r) \times \omega _{rb}}\right), \end{equation}
where b and r denote the block and ray number, respectively. Both tomography and spatial averaging results are given in the Section 4.3.

RESULTS

Single-station splitting parameters are mainly discussed with regard to two depth units; shallow (≤100 km) and deep (>100 km). The 100-km depth contour roughly represents the minimum depth limit of the lithosphere-asthenosphere boundary (LAB) in South Island (Stern et al.2000; Scherwath et al.2002). Figs 5(a) and (b)/(ci/ii) show the normalized circular histograms of fast azimuths (ϕ) measured at each station using shallow and deep local earthquakes, respectively. These spatial maps represent ϕ of the high-quality, A grade (Section 2), measurements that were obtained from multiple filters (a maximum of three filters with SNR ≥ 3 for each event). Events that have the same ϕ for different filters are visually more prominent than those with varying azimuths for multiple filters. Tables 1(a) and (b) summarize the average splitting parameters obtained at each station using deep and shallow earthquakes, respectively. To calculate single station averages (Table 1a/b), we use a single A grade result that satisfies the best filter criterion (Section 2) for each station-earthquake pair.

Figure 5.

(a) The spatial distribution of ϕ (A grade) obtained from all events is plotted as circular histograms (rose diagrams) at the stations. (b) and (c i/ii) show the A grade ϕ measured from shallow and deep earthquakes, respectively. Blue line on the Fig. 5(b) outlines the approximate limit of the plate-boundary parallel anisotropy, which we interpreted as the limit of pervasive simple shear. (d) Rose diagrams of ϕ from the stations on the Alpine fault. The consistent WNW–ESE alignment is significant at deeper levels (rose diagrams on the right side of dashed-line). Top-right number denotes the number of measurements at each depth level (bottom-right value).

Figure 5.

(a) The spatial distribution of ϕ (A grade) obtained from all events is plotted as circular histograms (rose diagrams) at the stations. (b) and (c i/ii) show the A grade ϕ measured from shallow and deep earthquakes, respectively. Blue line on the Fig. 5(b) outlines the approximate limit of the plate-boundary parallel anisotropy, which we interpreted as the limit of pervasive simple shear. (d) Rose diagrams of ϕ from the stations on the Alpine fault. The consistent WNW–ESE alignment is significant at deeper levels (rose diagrams on the right side of dashed-line). Top-right number denotes the number of measurements at each depth level (bottom-right value).

Some stations (e.g. MSZ, NNZ, KHZ, CROE) exhibit either large scattering or bimodal distribution of ϕ (or high SE) (Fig. 5a and Table 1a/b). We attribute these complex patterns of ϕ to spatially varying anisotropy and depth dependent ϕ as seen in the detailed analysis of splitting parameters (Fig. 6). In general, ϕ observed from both shallow and deep events indicate significant lateral variations compared to that of SKS azimuths (Section 4.4), but measurements at some closely spaced stations tend to show regional consistencies over short distances suggesting coherent anisotropy in small regions (Fig. 5). Shallow events processed at stations in the northern South Island display two distinct fast polarization azimuths (Fig. 5b): subparallel to the Alpine fault and the Marlborough fault system, mainly in the western and central part of northern South Island (QRZ, KELY, DSZ and THZ); and approximate EW alignment at a high angle to the Alpine fault trace in the eastern side of the northern South Island (LTZ, CRLZ, OXZ and MQZ). ϕ also becomes approximately EW oriented from deeper earthquakes (>100 km) as observed at coastal stations in northern South Island (DSZ,QRZ and KHZ) and the western part of southern South Island (MSZ and JCZ). For shallow earthquakes, NE oriented (plate-boundary subparallel) ϕ are observed at most of the stations in northern South Island and the central region of central South Island and southern South Island (e.g. CASS, RPZ, LBZ, WKZ, MLZ). However, the stations that are very close to the Alpine fault in central South Island (JCZ and FOZ) show either bimodal distributions or WNW–ESE fast azimuths at high angle to the fault (Fig. 5). This nearly fault-perpendicular azimuth also appears to be the prominent azimuth at depths >100 km as seen in JCZ and MSZ (Figs 5d and 6b), indicating a strong WNW–ESE alignment in anisotropy at subcrustal and/or upper-mantle depths in the western side of southern South Island.

Figure 6a.

Changes in ϕ with depth in northern South Island along the profile AA (see Fig. 4a). Station measurements displayed in the upper, middle and bottom diagrams are located in western, central and eastern regions of the northern South Island, respectively. Depth boundary where there is a visible change in ϕ with depth is marked with a black dashed line (i.e DSZ and THZ) and the colours of the bars denote the individual δt measurements. Note that the ϕ are determined in the horizontal plane and are thus plotted relative to N as marked on the figures. The blue cross indicates the reference point of each profile that is used to calculate the relative distance to the station. The approximate position of the slab (grey line), possible depth range for the location of LAB—shaded grey region, and the location of the Alpine/Wairau (AF/Wr) and Murchison-basin (Mu) faults are shown in each profile. The direction of the black bar next to each fault label indicates the approximate strike direction of the fault in map view (not the dip).

Figure 6a.

Changes in ϕ with depth in northern South Island along the profile AA (see Fig. 4a). Station measurements displayed in the upper, middle and bottom diagrams are located in western, central and eastern regions of the northern South Island, respectively. Depth boundary where there is a visible change in ϕ with depth is marked with a black dashed line (i.e DSZ and THZ) and the colours of the bars denote the individual δt measurements. Note that the ϕ are determined in the horizontal plane and are thus plotted relative to N as marked on the figures. The blue cross indicates the reference point of each profile that is used to calculate the relative distance to the station. The approximate position of the slab (grey line), possible depth range for the location of LAB—shaded grey region, and the location of the Alpine/Wairau (AF/Wr) and Murchison-basin (Mu) faults are shown in each profile. The direction of the black bar next to each fault label indicates the approximate strike direction of the fault in map view (not the dip).

Figure 6b.

Changes in splitting parameters with depth in southern South Island along the profile BB (see Fig. 4b). Note the higher δt and complex patterns of ϕ observed at the stations JCZ, WKZ and EAZ that are located in the northern region of the southern South Island close to the Alpine fault. Black dashed line (e.g. JCZ) indicates the depth boundary where there is a clear change in ϕ (colours of the bars denote the individual δt measurements).

Figure 6b.

Changes in splitting parameters with depth in southern South Island along the profile BB (see Fig. 4b). Note the higher δt and complex patterns of ϕ observed at the stations JCZ, WKZ and EAZ that are located in the northern region of the southern South Island close to the Alpine fault. Black dashed line (e.g. JCZ) indicates the depth boundary where there is a clear change in ϕ (colours of the bars denote the individual δt measurements).

INTERPRETATION AND DISCUSSION

Depth extent and backazimuthal variations of splitting parameters

An analysis of the splitting parameters with respect to depth and backazimuth is carried out in order to understand the complex patterns of ϕ observed from the rose diagrams. This analysis is limited to northern and southern South Island where there is a good earthquake distribution with depth due to subducting-slab earthquakes (Fig. 4). Stations that sample different regions of the two subduction zones with considerable numbers of splitting measurements exhibit changes in ϕ with depth (e.g. DSZ, QRZ in Fig. 6a and JCZ, MLZ in Fig. 6b). These depth variations sometimes appear to be combined with lateral variation, indicated by a complex patterns of ϕ (Fig. 6). We also search for variation in delay times (δt) to determine the extent and the degree of the anisotropy. In the presence of homogeneous or slowly varying anisotropy, δt accumulates along the earthquake-station ray path as a result of the difference in travel velocities of two split shear waves (Rumpker & Silver 1998). At most of the stations, δt are within the range of ∼0.15–0.30 s, but some stations recorded high δt (>0.4 s) with ϕ subparallel to the Alpine fault. The west coast close to the Alpine fault in southern South Island is also characterized by varying anisotropic structure as shown in Fig. 6(b) (JCZ, WKZ and EAZ), with inconsistent patterns of ϕ and high δt. Such variations imply a complex tectonic setting in the Fiordland region. However, clusters of closely spaced earthquakes give consistent readings at each station (e.g. Figs 6 and A2).

ϕ perpendicular to the strike of the Alpine fault observed from the intermediate (50–100 km) to deep measurements (>100 km) that are made on some stations (e.g. JCZ, MSZ and FOZ) right above and close to the fault (Fig. 5d) could be explained by two alternative mechanisms. If the upper-most mantle (∼30 km thick upper-mantle layer) contributes most of the anisotropy, it is likely a product of deformation of olivine crystals under different physical and chemical conditions. As proposed by Karato et al. (2008), deformational fabrics of olivine could be influenced by stress, temperature and the water content of the material resulting in changes of olivine slip directions from typical [100] to [010] (B-type olivine). The subduction to collision tectonic setting in the western edge of the central and southern South Island could be a favourable place for high pressure and/or water-rich environment that could change the orientation of ϕ from the expected fault subparallel ϕ and result the ϕ to be orthogonal to the plate-boundary and the Alpine fault. The other possibility is that the deep high-frequency phases are sensitive to crustal anisotropic structures (Audoine et al.2000). The prominent WNW–ESE direction at the above stations is parallel to the principal compressive direction (∼115°) of this region (Boese et al.2012; Townend et al.2012), suggesting the contribution of crustal stress-controlled mechanisms.

Apparent spatial variations of ϕ can also be caused by the geometry of the anisotropic unit. Models of vertically varying structure with two (Silver & Savage 1994) or more layers (Rumpker & Silver 1998) of anisotropy that have different azimuths of horizontal symmetry or have plunging symmetry axes (Babuska et al.1984; Babuska & Plomerova 1993) yield complex but systematic patterns of splitting parameters (see Section 4.2). However, modelling of these effects has been carried out mainly on SKS splitting studies focusing on the anisotropic structure of the upper mantle and have limited applications for the high-frequency local S-wave splitting because of the crustal heterogeneities and scattering that also affect the splitting. Lack of strong correlation between splitting parameters and backazimuth/incidence angle (Fig. A2) in our splitting data set suggests splitting parameters are not affected solely by single dipping layers of anisotropy.

Anisotropy in the two subduction systems

According to the ray geometry, seismic anisotropy sampled by stations LTZ and KHZ in northern South Island and MSZ and DCZ in southern South Island are most likely coming from the slab and subslab regions of the Hikurangi subduction zone (southern most part) and the Fiordland subduction zone, respectively (Figs 4 and 7). LTZ and DCZ show approximately the same azimuths for most of the slab events (Fig. 6) and also indicate an increase in δt with depth (Fig. 7). Assuming 4.5 km s−1 average S-wave velocity for a single anisotropic layer, we have calculated an approximate percent velocity anisotropy (Savage 1999) in the full depth range (∼0–240 km in Hikurangi subduction zone and ∼0–120 km Fiordland subduction zone) of the slab region (Fig. 7) using the highest quality splitting measurements (A grade, best filter and λmax > 8; Appendix A2). Our calculations yield about average 0.20 per cent (±0.10 per cent) shear-wave velocity anisotropy for the Hikurangi subducting slab region (Fig. 7a-ii) and 0.47 per cent (±0.17 per cent) velocity anisotropy for the Fiordland subducting slab (Fig. 7b-ii). Velocity anisotropy calculated for the Hikurangi subduction zone is consistent with the 0.23 per cent (±0.34 per cent) velocity anisotropy found by Audoine et al. (2000) for the northern most region of northern South Island. Yet, they both are small compared to higher shear-wave velocity anisotropy (3.1 ± 2.2 per cent) observed from 70-km-deep zone in the Wellington region that is located further north in the Hikurangi subduction zone (Matcham et al.2000). The ray paths suggest that deep slab events (>70 km) (Fig. 7) in northern South Island are most likely sampling the subslab region and only shallow events recorded on these stations should be sensitive to shallow slab anisotropy. When we only consider shallow slab events (≤70 km), our calculation yields 2.01 ± 0.74 per cent velocity anisotropy (appendix, Fig. A3), within the error range of Matcham et al. (2000). The anisotropy could be under estimated if the re-splitting has occurred (Sections 4.2 and 4.4), in which case the anisotropic region is smaller than the whole path.

Figure 7.

Anisotropy in the two subduction systems. Left- and right-hand plots show ray paths (a and b) and δt measured from events originating in Hikurangi and Fiordland subducting slabs, respectively. Possible depth range at which the Moho and LAB are located is shown by the grey shaded region. Here most of the ray paths sample the subducting slab region. Yellow (LTZ and MSZ) and blue (KHZ and DCZ) dots indicate δt and the black line shows the linear regression fit that is used to calculate the velocity anisotropy of the two subduction zones.

Figure 7.

Anisotropy in the two subduction systems. Left- and right-hand plots show ray paths (a and b) and δt measured from events originating in Hikurangi and Fiordland subducting slabs, respectively. Possible depth range at which the Moho and LAB are located is shown by the grey shaded region. Here most of the ray paths sample the subducting slab region. Yellow (LTZ and MSZ) and blue (KHZ and DCZ) dots indicate δt and the black line shows the linear regression fit that is used to calculate the velocity anisotropy of the two subduction zones.

The steepness of the Fiordland subducting slab results in a narrow spatial extent of the subduction zone. Therefore, complexities due to anisotropy contributions from the mantle wedge and subslab zones are limited in this region. The velocity anisotropy calculated for the Fiordland subducting slab is higher than that calculated for southern most Hikurangi suducting slab (Fig. 7). Shallow slab events (≤100 km) recorded on MSZ yield Alpine fault subparallel ϕ (Fig. 5b), suggesting fault-induced anisotropy at crustal depths in Fiordland. Fairly low percent velocity anisotropy and ϕ that are parallel to structures (e.g. strike of the slab and faults) obtained for both subduction zones suggests that the anisotropy in the subducting slab could be overprinted by crustal anisotropy. The δt measured at MSZ and LTZ are also sensitive to the average period of the shear-wave (Fig. 9). Therefore, the measured anisotropy from the slab events could be governed by more than a single source of anisotropy, such as crustal anisotropy or else the slab itself may consist of more than a single source of anisotropy (e.g. aligned microcracks, bending-induced faulting and mineral alignment).

Evidence for heterogeneous anisotropic structure

Scatter in δt and spatially varying ϕ imply complex anisotropic structure beneath South Island. Complex and/or multiple layered anisotropy can be investigated by studying changes in splitting parameters as a function of backazimuth (BAZ) and ϕp (Silver & Savage 1994; Rumpker & Silver 1998; Long & Silver 2009) or with dominant period (T) of the S-wave. ϕp is the polarization azimuth of the unsplit shear-wave before entering the anisotropic medium. In the case of re-splitting, this polarization azimuth is equivalent to the azimuth of the fast shear-wave that resulted from the preceding anisotropic medium. At DSZ in the western edge of the northern South Island, δt varies with ϕp (Fig. 8b), and varies with depth (Fig. 8c), but other variables failed to yield strong correlation (e.g. backazimuth, magnitude, period of the S-wave).

Figure 8.

Changes in δt with ϕp (b) and changes in ϕ with depth recorded at DSZ in Murcheson basin area (Fig. 1) in the western northern South Island. (a) Map view of the events recorded. (c) Depth profile of the ϕ along the line AA. ϕ and the events are colour coded according to ϕp and the dashed line indicates the change in ϕ with depth. Here the rose diagram shows the histogram of ϕ obtained from events ≤160 km and colour coded according to ϕp. (d) Schematic diagram illustrating the different anisotropic zones (P is a region with E–W aligned anisotropy, Q is a layer with NW–SE anisotropy and R shows anisotropic patches with NE–SW orientation), which result in the observed patterns of ϕ in Fig. 8(c). (e-i) 2-D inversion model of electrical resistivity for northern South Island (Wannamaker et al.2009; doi:10.1038/nature08204). The low resistivity body in Murchison basin area and north of the Marlborough faults are marked by C and B, respectively. (e-ii) Depth profile of the Vp/Vs model for central northern South Island (Eberhart-Phillips & Bannister 2010; doi:10.1111/j.1365-246X.2010.04621.x). Dashed line box shows the rough location of the electrical resistivity model explained in Fig. 8(e-i).

Figure 8.

Changes in δt with ϕp (b) and changes in ϕ with depth recorded at DSZ in Murcheson basin area (Fig. 1) in the western northern South Island. (a) Map view of the events recorded. (c) Depth profile of the ϕ along the line AA. ϕ and the events are colour coded according to ϕp and the dashed line indicates the change in ϕ with depth. Here the rose diagram shows the histogram of ϕ obtained from events ≤160 km and colour coded according to ϕp. (d) Schematic diagram illustrating the different anisotropic zones (P is a region with E–W aligned anisotropy, Q is a layer with NW–SE anisotropy and R shows anisotropic patches with NE–SW orientation), which result in the observed patterns of ϕ in Fig. 8(c). (e-i) 2-D inversion model of electrical resistivity for northern South Island (Wannamaker et al.2009; doi:10.1038/nature08204). The low resistivity body in Murchison basin area and north of the Marlborough faults are marked by C and B, respectively. (e-ii) Depth profile of the Vp/Vs model for central northern South Island (Eberhart-Phillips & Bannister 2010; doi:10.1111/j.1365-246X.2010.04621.x). Dashed line box shows the rough location of the electrical resistivity model explained in Fig. 8(e-i).

Multiple-layered anisotropy with horizontal symmetry axes produce a characteristic pattern between ϕp and apparent splitting parameters with periodicity of π/2 (Silver & Savage 1994). However, the patterns are only detectable with long-period (Tt > 5) waves (Rumpker & Silver 1998) from a wide range of ϕp. Tt values obtained from the local and regional S-phases used in this analysis are within the range of ∼0.3–20. Therefore, seismograms that have long period waveforms (∼11 per cent have Tt > 5) are likely yielding apparent splitting, but those with smaller ratios of Tt ≤ 1 (∼18 per cent) probably measuring splitting parameters corresponding to the upper layer. We do not see a π/2 periodic pattern in splitting parameters with respect to ϕp as observed in long period waveform studies (Silver & Savage 1994), but instead we see a jump in δt at about ϕp 80° (Fig. 8b ). Higher δt tend to have ϕp between 20° and 80° (blue) and the lower δt have ϕp in the range of 80°–170° (red). At depths > ∼160 km both sets of ϕp approximately give the same E–W ϕ. However, at depths ≤160 km, initial polarizations between 20° and 80° yield N–S or NW–SE ϕ, whereas ϕ with 80°–170° initial polarizations are mostly NE (see rose diagram of ϕ in Fig. 8c).

Changes in both δt with ϕp and ϕ with depth may result from a localized subcrustal anisotropic zone. Wannamaker et al. (2009) found that the Murchison basin area fault zone and Marlborough faults are characterized by moderately low resistivity structures from ∼25 to 50 km. They infer that these structures are produced by fluid expelled from the deep (∼100 km) subduction interface and that they are dominated by fluid filled cracks and shear fabrics (Fig. 8e-i, Wannamaker et al.2009), which should yield anisotropy parallel to the plate-boundary. Moderately high Vp/Vs areas identified via tomographic inversion of local earthquake traveltimes by Eberhart-Phillips & Bannister (2010) in the same region also suggest high-fluid pressure (Fig. 8e-ii). Northwest of the major fault systems, low Vp/Vs and high Vp (Eberhart-Phillips & Bannister 2010) in the upper mantle is interpreted to be strong Australian lithosphere, which may not be deformed by present processes (Eberhart-Phillips & Bannister 2010). Thus, the ray paths passing through different zones will measure different ϕ.

We suggest that the waves from earthquakes deeper than 160 km, which travel on the paths mostly to the NW of DSZ and display mainly EW ϕ (Figs 8c and d), are passing through an anisotropic zone with EW fast orientation. This zone could be either in the mantle lithosphere or crust; if it is in the mantle lithosphere it could be left from the Gondwana Cretaceous deformation (Eberhart-Phillips & Bannister 2010) and if it is in the crust it could be caused by the present stress orientation (Audoine et al.2000). A third possibility is that the E–W alignment is a near source effect, because the earthquakes are within a region of high Vp/Vs ratios at the deep slab-edge (Fig. 8e-ii, Eberhart-Phillips & Bannister 2010). We favour the crustal stress interpretation as the ϕ is nearly parallel to the crustal stress direction (yellow bars in Fig. 11, Townend et al.2012), and short period waveforms are often most affected by the upper-most path. The shallower ray paths are in the lithospheric mantle and the crust to the SE of station DSZ, and pass through the Marlborough fault system, which has patches (zone c in Fig 8e-i) of high-fluid content. These rays pass through a region in the lower crust/upper mantle that has roughly NE/SW oriented fast anisotropy induced by plate-boundary parallel shearing (zone S in Fig. 8d), and then pass through another layer with roughly NW/SE fast anisotropy (zone Q in Fig. 8d), which re-splits them. This NW/SE anisotropy could be generated by asthenospheric counterflow perpendicular to the plate-boundary or by some other processes. Some waves continue to DSZ without traveling through any other regions. They have ϕ between N and NW, δt of about 0.3 s and NE/SW incoming polarizations obtained as they were traveling through the lowest layer as fast waves. Other waves pass through a third patch of NE/SW oriented ϕ (patchy area marked by R in Fig. 8d). They have NE/SW ϕ with shorter δt of ∼0.1 s, and NW/SE incoming polarizations (Figs 8b and c). Perhaps these shallow patches are in the fault zone or are regions of schist. Cycle skipping could also produce a dependence between δt and ϕp (Walsh et al.2012). However, the contour plots of the large δt (Fig. 3b) suggest this is not occurring.

Other than at the stations that sample slab anisotropy, values of δt show no clear increase with depth. Instead, there is a scatter in δt that is measured on stations in the central northern South Island (e.g. THZ, KELY) and a few stations (e.g. WHZ, APZ) in southern South Island (appendix, Fig. A4). Both scatter in δt and changes in ϕ at nearby stations imply the presence of heterogeneous anisotropic structure. Regions where changes occur over short distances laterally and vertically are likely characterized by variable sources of anisotropy (e.g. shear-induced structures, microcracks, mineral alignment due to flow and large-scale structural fabric). Some studies suggest that different S-wave frequencies sample different scale-lengths of anisotropy, hence produce variable results depending on heterogeneity of the medium (Marson-Pidgeon & Savage 1997; Wirth & Long 2010). We test this hypothesis by analyzing the highest quality splitting parameters (A grade, best filter and λmax > 8; Appendix A2) at each station (Fig. 9). MSZ, WHZ, APZ and EAZ show moderately high correlation coefficient (CCEF) between period and δt, i.e. frequency-dependent splitting (e.g. CCEF > 0.4) and all of these stations are located in the southern South Island. The dependency of δt on the period of the S-wave reflects heterogeneous anisotropic structure in the Fiordland subduction zone, and SKS splitting can be explained by shallow anisotropy in this area. However the lack of frequency-dependent splitting at the other stations suggests the large SKS δt with consistent ϕ (Klosko et al.1999) can not solely be explained by frequency dependence and therefore require a deep seated source.

Figure 9.

Examples of stations that show changes in δt with period of the S-waves in northern and southern South Island. Black line indicates the least squares fit of the data. Two plots with grey circles show the combined plots of individual station data in the southern and northern South Island. Here the CCEF values for each station is marked on the corresponding station plot. Note the moderately high CCEF (>0.4) calculated at the southern South Island stations.

Figure 9.

Examples of stations that show changes in δt with period of the S-waves in northern and southern South Island. Black line indicates the least squares fit of the data. Two plots with grey circles show the combined plots of individual station data in the southern and northern South Island. Here the CCEF values for each station is marked on the corresponding station plot. Note the moderately high CCEF (>0.4) calculated at the southern South Island stations.

Interpretations from 2-D delay-time tomography and spatial averages

The inconsistency of measurements from closely spaced earthquakes at different stations suggests that lateral variations in anisotropy occur. This in turn suggests that the 2-D delay-time tomographic approach might be useful to determine the average structure. Fig. 10 illustrates the 2-D delay-time tomography results obtained using AB quality splitting parameters calculated from 2615 event-station pairs. Because of the depth dependence of splitting parameters, we produce 2-D distribution maps for two depth sections (depth ≤ 100 km and depth > 100 km) and these maps provide average spatial variations of δt. In Fig. 10 (shallow/deep), warm colours (yellow to red) indicate the regions with high strength of anisotropy (> ∼0.004 s km−1) compared to blue colours with weak

strength of anisotropy (≤ ∼0.002 s km−1). The shading indicates the region with low resolution and/or large residuals. Anisotropy strength profiles that are within the higher resolution region reasonably agree with the possible anisotropic sources in South Island. For depths less than 100 km, crustal fault zones and the areas with schist rock (compare to Fig. 1) appear to be the regions with strong anisotropy. This is clear from the high-anisotropy strength yields around the Alpine/Wairau faults, Marlborough fault zone, Murchison basin area faults in the northern South Island (Fig. 1) and some crustal faults in the southern South Island (Fig. 10-shallow). For the deep earthquakes (>100 km), the region with high strength of anisotropy in the northern part of the southern South Island (SE of the Alpine fault) coincides quite well with the region that has large crustal thickness (∼48 km) (Stern et al.2000; Bourguignon et al.2007) and large negative Bouguer anomaly in and around the Southern Alps (Davey et al.2007). Deformation associated with crustal shortening and mantle lithosphere (upper ∼150 km) thickening as a result of continental collision could be the main mechanisms that governs the anisotropic structure in this region. As suggested by Norris et al. (1990), this region (SE of the Alpine fault) might be dominated by distributed deformation that results from the plate motion mechanisms.

Figure 10.

2-D delay-time tomography from TESSA (Johnson et al.2011). Tomography maps are created for shallow (≤100 km) and deep (>100 km) levels. Warm colours (red to yellowish green) in the maps indicate region with higher strength of anisotropy compared to weak zone (blue colour) and the shaded region shows the regions with low resolution and/or large residuals.

Figure 10.

2-D delay-time tomography from TESSA (Johnson et al.2011). Tomography maps are created for shallow (≤100 km) and deep (>100 km) levels. Warm colours (red to yellowish green) in the maps indicate region with higher strength of anisotropy compared to weak zone (blue colour) and the shaded region shows the regions with low resolution and/or large residuals.

Weighted mean ϕ obtained for each grid block are shown in Fig. 11. Due to lack of spatial resolution of deep measurements (Fig. 4), spatial averaging is limited to depths ≤100 km. Even at shallow depths of eastern-southern South Island, measurements are made in a limited region as a result of sparse distribution of ray paths and therefore, spatial averages are not well constrained at some locations (see Figs 4 and 11). Spatial averages of ϕ obtained from shallow events display consistent regional patterns. ϕ in the red shaded regions are either subparallel to the plate-boundary or aligned parallel to the translational faults (e.g. Marlborough fault zone) in northern South Island. The consistent pattern of ϕ with NS and NNW alignment in the region surrounding the red-shaded zone in southern South Island agree with the SKS measurements that were made in those regions (red bars in Fig. 11, Savage et al.2007a). Thus, it is possible that most of the anisotropy measured in this region is coming from the mantle lithosphere or from coupling between the crustal and the mantle-lithospheric anisotropy. However, the NE-SW alignment in the southern red-shaded zone shows a sharp change in ϕ from the general trend in this region and is aligned parallel to the Pn anisotropy fast azimuth (blue bars in Fig. 11, Scherwath et al.2002), the plate-boundary and the strike of the high topography (e.g. Southern Alps) of the region. This zone suggests an abrupt changes in strain field and is characterized by high strength of anisotropy as seen in Fig. 10 (shallow).

Figure 11.

Spatially (25-km grid blocks) averaged fast azimuths ($$\bar{\phi _b}$$, eq. 3) using TESSA. Map displays the spatial averages from shallow (≤100 km) events. Black colour bars denote the $$\bar{\phi _b}$$ with SE ≤ ±25° and the grey bars indicate the SE of $$\bar{\phi _b} > \pm 25^{\circ }$$. Regions with plate-boundary subparallel $$\bar{\phi _b}$$ are denoted by red shaded regions (note that the sharp change in ϕ from NW–SE trend to NE–SW trend in southern South Island is marked by the shaded red region in the south). Blue dashed line marks the boundary limits of the NE–SW trend that is aligned with strike of plate-boundary and the associated faults. Both individual station measurements (Fig. 5b) and spatial averages have been used to constrained these boundary limits. SKS (Klosko et al.1999; Savage et al.2007a), Pn anisotropy (Bourguignon et al.2007) and SHmax (Townend et al.2012) directions from previous studies are marked by red, blue and yellow bars, respectively.White irregular lines indicate the fault structures in South Island.

Figure 11.

Spatially (25-km grid blocks) averaged fast azimuths ($$\bar{\phi _b}$$, eq. 3) using TESSA. Map displays the spatial averages from shallow (≤100 km) events. Black colour bars denote the $$\bar{\phi _b}$$ with SE ≤ ±25° and the grey bars indicate the SE of $$\bar{\phi _b} > \pm 25^{\circ }$$. Regions with plate-boundary subparallel $$\bar{\phi _b}$$ are denoted by red shaded regions (note that the sharp change in ϕ from NW–SE trend to NE–SW trend in southern South Island is marked by the shaded red region in the south). Blue dashed line marks the boundary limits of the NE–SW trend that is aligned with strike of plate-boundary and the associated faults. Both individual station measurements (Fig. 5b) and spatial averages have been used to constrained these boundary limits. SKS (Klosko et al.1999; Savage et al.2007a), Pn anisotropy (Bourguignon et al.2007) and SHmax (Townend et al.2012) directions from previous studies are marked by red, blue and yellow bars, respectively.White irregular lines indicate the fault structures in South Island.

Spatial averages in the northern part of the central South Island are poorly constrained due to lack of measurements. SE associated with most of the spatial averages in eastern and central part of the central South Island are greater than 25°, indicating more than a single mode of ϕ. Although the prominent trend (WNW–ESE aligned grey bars in Fig. 11) is consistent with the horizontal maximum compressive direction (SHmax) of this region (yellow bars in Fig. 11, Townend et al.2012), spatial averages do not provide robust estimations due to large errors. Single-station measurements at RPZ, LBZ and CASS (Figs 5b) show that there is a plate-boundary subparallel trend in the central region, however this trend is masked by the WNW–ESE ϕ parallel to SHmax when considering spatial averages. Therefore, the contour line connecting the two red-shaded regions is determined comparing both Figs 5(b) and 11 (see Fig. 11).

Comparison with previous anisotropy studies and implications for plate-boundary deformation

Anisotropy in the crust and the lithospheric mantle in the Marlborough fault zone and the upper-most region of the northern South Island was previously investigated by Balfour et al. (2005) and Audoine et al. (2000), but the spatial and the depth resolution of these studies are limited. Balfour et al. (2005) argued that the anisotropy in Marlborough is not controlled by the regional stress and both Balfour et al. (2005) and Audoine et al. (2000) found that ϕ in those regions are roughly parallel to the plate-boundary (NE-SW). In agreement with the above studies, we find a similar average ϕ (spatial) for the same regions in northern South Island (Fig. 11) from shallow events. Most of our deep events Figs 6(a) yield the E–W/WNW–ESE trend presented by Audoine et al. (2000) from deep events at some stations. As suggested in Section 4.1, this E–W/WNW–ESE trend may result from either a deep anisotropic source (e.g. LPO fabrics in the upper mantle or lower crust) or a shallow stress-controlled source or the combination of both. The NE–SW trend is likely caused by pervasive shear associated with the plate-boundary deformation (Fig. 11).

When compared to SKS and Pn anisotropy studies in South Island (Klosko et al.1999; Scherwath et al.2002), ϕ obtained from local/regional data are variable. However, the spatially averaged ϕ from earthquakes with depths ≤100 km (Fig. 11) are mostly consistent with the SKS and Pn fast azimuths observed for the southern central South Island (regions where WKZ and LBZ are located), the eastern part of the southern South Island and northern South Island (compare red and black bars in Fig. 11). Our observations are also consistent with the crustal finite strain estimations from bending related geological markers in South Island (Little et al.2002). As suggested by Little et al. (2002), the consistency indicates coupling of the upper-mantle deformation (as inferred from SKS results) with crustal deformation more towards the Alpine fault and weak coupling in the eastern part of the central and southern South Island (Fig. 5). Based on the ϕ from shallow events (≤100 km) (Figs 5b and 11), we infer that the width of the plate-boundary subparallel ϕ that indicates pervasive simple shear (Molnar et al.1999; Little et al.2002) is limited to about 130 km southeast from the Alpine fault in the southern part of central South Island (regions where WKZ and LBZ are located) and the northern region of southern South Island. A region with high-geodetic shear strain rates in South Island as proposed by Beavan et al. (1999) also lies within the estimated width of the plate-boundary induced deformation in our study. The western boundary limit of the near offshore region of northern and central South Island are not yet known. Also, the lack of deep events in the central South Island does not allow accurate constraints on the depth extent of this trend.

As proposed by Molnar et al. (1999), if the total strain produced by the pervasive shear is confined to the top 100 km of the upper mantle, a large portion of SKS δt anisotropy should be accommodated in the lithosphere. However for local/regional events, even though the ray paths travel through the same region of the upper mantle, the observed station δt average of local/regional splitting measurements is about 9–20 per cent of that observed by SKS studies. Thus, small average station δt (∼0.1–0.4 s) suggest either: (i) the majority of SKS δt are from asthenosphere and/or deeper part of the lithosphere or (ii) SKS and local phases sample the same structure, but the latter have higher frequencies and are seeing smaller δt due to frequency-dependent splitting (Fig. 9 and Section 4.2) or (iii) higher frequency phases are more susceptible to re-splitting from small-scale structures that are present at the last layer of a layered anisotropic media (e.g. Fig. 8). Possibility (ii) appears to be a localized effect (Section 4.2) and therefore either (i) and/or (iii) can be used to explain the δt discrepancies between local (this study) and SKS (previous work) splitting.

CONCLUSIONS

Using both local and regional earthquake data, we have obtained a S-wave splitting data set to provide constraints on the lateral and vertical distribution of seismic anisotropy in the lithosphere of the South Island. The application of quality criteria and statistical analysis allow us to extract well resolved splitting parameters and to provide more confident results from the complex splitting associated with local/regional S-wave splitting analysis. Therefore, the spatial and depth variations observed from the detailed splitting analysis and results from the tomography provide a good approximation for the spatial distribution of anisotropy (Sections 4.1, 4.2 and 4.3). In general, splitting measurements in South Island shows regional consistencies over short distances, which imply localized coherent anisotropy. Spatial and depth variations in splitting parameters suggest the existence of several deformation domains (e.g. Figs 11 and 8d). Anisotropy in these domains appears to be controlled by different source mechanisms: shear-induced fabrics related to present-day plate-boundary deformation (in northern and central-southern South Island), stress-controlled anisotropy linked with prevailing stress (in eastern part of the central South Island) and anisotropy deeper than 50 km due to the mantle anisotropy within the lithosphere and/or asthenosphere that reflects the paleo-anisotropy from previous tectonic episodes or present asthenospheric flow (in southern and eastern parts of the Southern South Island).

2-D delay-time tomography and spatial averaging of ϕ give consistent patterns on a regional scale and suggest that the schist rock may also have a significant contribution to the total splitting (∼0.1–0.4 s) measured. A NE–SW trend of the spatial averages of northern South Island consistent with the strike of the plate-boundary and the strikes of the crustal faults, suggests that the dominant source of anisotropy is plate-boundary induced structures. An abrupt variation in ϕ from NW–SE to NE–SW in southern South Island (Fig. 11) suggests a high-strain structure at crustal or subcrustal depths. This zone also indicates fairly high strength of anisotropy (∼0.006 s km−1) and resides in the region with high-crustal thickness (∼48 km) and inferred thick lithosphere (Stern et al.2000; Bourguignon et al.2007). The origin of this anisotropic structure may be related to pervasive shear from plate-boundary deformation.

Assuming that the plate-boundary subparallel ϕ is a product of distributed shear associated with plate-boundary deformation, we suggest that the zone of shear deformation extends ∼130 km from the Alpine fault in the northern part of the southern South Island (eastern boundary limit of the red shaded region in South), but in northern and central South Island this zone is more distributed than in southern South Island. Unless there is re-splitting due to misaligned anisotropic media, and/or frequency-dependent splitting affecting the local/regional splitting, small δt observed from local S-wave splitting compared to large δt (∼1.7 s) observed from SKS splitting suggest that most of splitting from SKS is accommodated either in the asthenosphere or in the deeper part of thickened mantle lithosphere (depths ∼150–250 km) where we do not have enough depth coverage of splitting measurements. When considering upper-lithospheric depths (∼80–120-km-thick lithosphere) in the central and upper-most southern South Island, our splitting results fit best with the localized deformation model (Fig. 1b).

We acknowledge the support of M. Henderson, A. Seward, C. Boese during the field work and the contribution of S. Greve and K. Jacobs with data processing and analysis. We thank New Zealand GeoNet (http://www.geonet.org.nz/) project for providing seismic data and IRIS, PASSCAL for providing instruments for temporary deployments. Financial support for this project is given by Victoria University of Wellington (PhD scholarship), National Science Foundation and Foundation for Research, Science, and Technology. We are grateful to J. Johnson for providing the TESSA code for public use and for useful suggestion on improving our work. We also thank T. Stern, A. Sheehan, P. Molnar, J. Collins, E. Smith and J. Townend for their useful discussions and comments on this work. J. Plomerova and anonymous reviewer provided helpful comments.

REFERENCES

Abt
D.L.
Fischer
K.M.
Resolving three-dimensional anisotropic structure with shear wave splitting tomography
Geophys. J. Int.
 , 
2008
, vol. 
173
 
3
(pg. 
859
-
886
)
Aki
K.
Richards
P.G.
Quantitative Seismology
 , 
2002
2nd edn
Sausalito, CA
University Science Books
Audoine
E.
Savage
M.K.
Gledhill
K.
Seismic anisotropy from local earthquakes in the transition region from a subduction to a strike-slip plate boundary, New Zealand
J. geophys. Res. (Solid Earth)
 , 
2000
, vol. 
105
 
B4
(pg. 
8013
-
8033
)
Avouac
J.
Tapponnier
P.
Kinematic model of active deformation in central asia
Geophys. Res. Lett.
 , 
1993
, vol. 
20
 
10
(pg. 
895
-
898
)
Babuska
V.
Plomerova
J.
Lithospheric thickness and velocity anisotropy - seismological and geothermal aspects
Tectonophysics
 , 
1993
, vol. 
225
 
1–2
(pg. 
79
-
89
)
Babuska
V.
Plomerova
J.
Sileny
J.
Large-scale oriented structures in the subcrustal lithosphere of central europe
Ann. Geophys.
 , 
1984
, vol. 
2
 (pg. 
649
-
662
)
Balfour
N.J.
Savage
M.K.
Townend
J.
Stress and crustal anisotropy in marlborough, New Zealand: evidence for low fault strength and structure-controlled anisotropy
Geophys. J. Int.
 , 
2005
, vol. 
163
 
3
(pg. 
1073
-
1086
)
Beavan
J.
, et al.  . 
Crustal deformation during 1994-1998 due to oblique continental collision in the central southern Alps, New Zealand, and implications for seismic potential of the alpine fault
Journal of Geophysical Research-Solid Earth
 , 
1999
, vol. 
104
 
B11
(pg. 
25233
-
25255
)
Boese
C.M.
Townend
J.
Smith
E.
Stern
T.
Microseismicity and stress in the vicinity of the Alpine Fault, central Southern Alps, New Zealand
J. Geophys. Res.
 , 
2012
, vol. 
117
 pg. 
B02302
  
doi:10.1029/2011JB008460
Boness
N.L.
Zoback
M.D.
Mapping stress and structurally controlled crustal shear velocity anisotropy in California
Geology
 , 
2006
, vol. 
34
 
10
(pg. 
825
-
828
)
Booth
D.
Crampin
S.
Shear-wave polarizations on a curved wavefront at an isotropic free-surface
Geophys. J. R. astr. Soc.
 , 
1985
, vol. 
83
 
1
(pg. 
31
-
45
)
Bourguignon
S.
Stern
T.A.
Savage
M.K.
Crust and mantle thickening beneath the southern portion of the southern Alps, New Zealand
Geophys. J. Int.
 , 
2007
, vol. 
168
 
2
(pg. 
681
-
690
)
Bourne
S.J.
England
P.C.
Parsons
B.
The motion of crustal blocks driven by flow of the lower lithosphere and implications for slip rates of continental strike-slip faults
Nature
 , 
1998
, vol. 
391
 
0028–0836
(pg. 
655
-
659
)
Cox
S.
Sutherland
R.
Regional geological framework of South Island, New Zealand, and its significance for understanding the active plate boundary
Am. Geophys. Union, Geophys. Monog.
 , 
2007
, vol. 
175
 (pg. 
19
-
46
)
Crampin
S.
The fracture criticality of crustal rocks
Geophys. J. Int.
 , 
1994
, vol. 
118
 
2
(pg. 
428
-
438
)
Crampin
S.
Ga
Y.A.
Volti
T.
Earthquake stress-forecasting: A new understanding of rock deformation
Earthquake: Hazard, Risk, and Strong Ground Motion
 , 
2004
(pg. 
95
-
106
)
Davey
F.J.
Smith
E. G.C.
The tectonic setting of the fiordland region, southwest New-Zealand
Geophys. J. R. astr. Soc.
 , 
1983
, vol. 
72
 
1
(pg. 
23
-
38
)
Davey
F.J.
, et al.  . 
Okaya
D.
Stern
T.
Davey
F.
Geophysical structure of the southern Alps orogen, South Island, New Zealand
A Continental Plate Boundary: Tectonics at South Island, New Zealand
 , 
2007
Washington, D. C.
Geophysical Monograph American Geophysical Union
(pg. 
47
-
73
)
Devs
M.
King
G.C.
Klinger
Y.
Agnon
A.
Localised and distributed deformation in the lithosphere: modelling the dead sea region in 3 dimensions
Earth planet. Sci. Lett.
 , 
2011
, vol. 
308
 
1–2
(pg. 
172
-
184
)
Eberhart-Phillips
D.
Bannister
S.
3-D imaging of Marlborough, New Zealand, subducted plate and strike-slip fault systems
Geophys. J. Int.
 , 
2010
, vol. 
182
 
1
(pg. 
73
-
96
)
Eberhart-Phillips
D.
Reyners
M.
A complex, young subduction zone imaged by three-dimensional seismic velocity, fiordland, New Zealand
Geophys. J. Int.
 , 
2001
, vol. 
146
 
3
(pg. 
731
-
746
)
Evans
R.
Effects of the free-surface on shear wavetrains
Geophys. J. R. astr. Soc.
 , 
1984
, vol. 
76
 (pg. 
165
-
172
)
Furlong
P.K.
Locating the deep extent of the plate boundary along the Alpine fault zone, New Zealnd: implications for patterns of exhumation in the southern Alps
Geol. Soc. Am. Special Papers
 , 
2007
, vol. 
434
 (pg. 
1
-
14
)
Gerst
A.
Savage
M.K.
Seismic anisotropy beneath ruapehu volcano: a possible eruption forecasting tool
Science
 , 
2004
, vol. 
306
 
5701
(pg. 
1543
-
1547
)
Johnson
J.H.
Savage
M.K.
Townend
J.
Distinguishing between stress-induced and structural anisotropy at Mount Ruapehu volcano, New Zealand
J. Geophys. Res.
 , 
2011
, vol. 
116
 pg. 
B12303
  
doi:10.1029/2011JB008308
Karato
S. H., J.
Katayama
I.
Skemer
P.
Geodynamic significance of seismic anisotropy of the upper mantle: new insights from laboratory studies
Ann. Rev. Earth Planet. Sci.
 , 
2008
, vol. 
36
 
1
(pg. 
59
-
95
)
Klosko
E.R.
Wu
F.T.
Anderson
H.J.
Eberhart-Phillips
D.
McEvilly
T.V.
Audoine
E.
Savage
M.K.
Gledhill
K.R.
Upper mantle anisotropy in the New Zealand region
Geophys. Res. Lett.
 , 
1999
, vol. 
26
 
10
(pg. 
1497
-
1500
)
Little
T.A.
Mortimer
N.
Rotation of ductile fabrics across the Alpine fault and cenozoic bending of the new zealand orocline
J. Geol. Soc.
 , 
2001
, vol. 
158
 (pg. 
745
-
756
)
Little
T.A.
Savage
M.K.
Tikoff
B.
Relationship between crustal finite strain and seismic anisotropy in the mantle, pacific-australia plate boundary zone, South Island, New Zealand
Geophys. J. Int.
 , 
2002
, vol. 
151
 
1
(pg. 
106
-
116
)
Long
M.D.
Silver
P.G.
Shear wave splitting and mantle anisotropy: measurements, interpretations, and new directions
Surveys Geophys.
 , 
2009
, vol. 
30
 
4-5
(pg. 
407
-
461
)
Mainprice
D.
Nicolas
A.
Development of shape and lattice preferred orientations: application to the seismic anisotropy of the lower crust
J. Struct. Geol.
 , 
1989
, vol. 
11
 
1/2
(pg. 
175
-
189
)
Mardia
K.V.
Statistics of Directional Data
 , 
1972
London and New York
Academic Press
Marson-Pidgeon
K.
Savage
M.K.
Frequency-dependent anisotropy in wellington, New Zealand
Geophys. Res. Lett.
 , 
1997
, vol. 
24
 
24
(pg. 
3297
-
3300
)
Matcham
I.
Savage
M.K.
Gledhill
K.R.
Distribution of seismic anisotropy in the subduction zone beneath the wellington region, New Zealand
Geophys. J. Int.
 , 
2000
, vol. 
140
 
1
(pg. 
1
-
10
)
Mattatall
L.R.
Fouch
M.J.
Small-scale variations in SKS splitting near Parkfield, California
Eos Trans. AGU
 , 
2007
, vol. 
88
 
52
 
Fall Meet. Suppl., Abstract T53C-02
Meade
B.J.
Hager
B.H.
Block models of crustal motion in southern California constrained by GPS measurements
J. Geophys. Res.
 , 
2005
, vol. 
110
 pg. 
B03403
  
doi:10.1029/2004JB003209
Molnar
P.
, et al.  . 
Continous deformation versus faulting through the continental lithosphere of New Zealand
Science
 , 
1999
, vol. 
286
 
5439
(pg. 
516
-
519
)
Moore
M.
England
P.
Parsons
B.
Relation between surface velocity field and shear wave splitting in the South Island of New Zealand
J. geophys. Res.
 , 
2002
, vol. 
107
 
2198
(pg. 
0148
-
0227
)
Nicolas
A.
Christensen
N.I.
Formation of anisotropy in upper mantle peridotite: a review
AGU Geodyn. Monog. Series
 , 
1987
, vol. 
16
 (pg. 
111
-
123
)
Norris
R.J.
Koons
P.O.
Cooper
A.F.
The obliquely-convergent plate boundary in the South Island of New-Zealand - implications for ancient collision zones
J. Struct. Geol.
 , 
1990
, vol. 
12
 
5–6
(pg. 
715
-
725
)
Nur
A.
Simmons
G.
Stress-induced velocity anisotropy in rock: an experimental study
J. geophys. Res.
 , 
1969
, vol. 
74
 
27
(pg. 
6667
-
6674
)
Nuttli
O.W.
The effects of the earth's surface on the S-wave particle motion
Bull. seism. Soc. Am.
 , 
1961
, vol. 
51
 (pg. 
237
-
246
)
Okaya
D.
Stern
T.
Davey
F.
Henrys
S.
Cox
S.
Continent-continent collision at the pacific/indo-australian plate boundary: background, motivation, and principal results
New Zealand Am. Geophys. Union Geophys. Monog.
 , 
2007
, vol. 
175
 (pg. 
1
-
18
)
Philipp
B.
Circstat: a MATLAB toolbox for circular statistics
J. Stat. Software
 , 
2009
, vol. 
31
 
10
(pg. 
1
-
21
)
Roman
D.C.
Savage
M.K.
Arnold
R.
Latchman
J.L.
De Angelis
S.
Analysis and forward modeling of seismic anisotropy during the ongoing eruption of the Soufrière Hills Volcano, Montserrat, 1996–2007
J. Geophys. Res.
 , 
2011
, vol. 
116
 pg. 
B03201
  
doi:10.1029/2010JB007667
Rumpker
G.
Silver
P.G.
Apparent shear-wave splitting parameters in the presence of vertically varying anisotropy
Geophys. J. Int.
 , 
1998
, vol. 
135
 
3
(pg. 
790
-
800
)
Savage
K.
Duclos
M.
Marson-Pidgon
K.
Seisimic anisotropy in South Island, New Zealand
Geophys. Monog. Series 175
 , 
2007a
Savage
M.K.
Seismic anisotropy and mantle deformation: what have we learned from shear wave splitting?
Rev. Geophys.
 , 
1999
, vol. 
37
 
1
(pg. 
65
-
106
)
Savage
M.K.
Ohminato
T.
Aoki
Y.
Tsuji
H.
Greve
S.M.
Stress magnitude and its temporal variation at Mt. Asama volcano, Japan, from seismic anisotropy and GPS
Earth planet. Sci. Lett.
 , 
2010
, vol. 
290
 
3-4
(pg. 
403
-
414
)
Scherwath
M.
Melhuish
A.
Stern
T.
Molnar
P.
Pn anisotropy and distributed upper mantle deformation associated with a continental transform fault
Geophys. Res. Lett.
 , 
2002
, vol. 
29
 
8
 
doi:10.1029/2001GL014179
Silver
P.G.
Seismic anisotropy beneath the continents: probing the depths of geology
Ann. Rev. Earth Planet. Sci.
 , 
1996
, vol. 
24
 (pg. 
385
-
432
)
Silver
P.G.
Chan
W.W.
Shear-wave splitting and subcontinental mantle deformation
J. geophys. Res. (Solid Earth)
 , 
1991
, vol. 
96
 
B10
(pg. 
16 429
-
16 454
)
Silver
P.G.
Savage
M.K.
The interpretation of shear-wave splitting parameters in the presence of two anisotropic layers
Geophys. J. Int.
 , 
1994
, vol. 
119
 
3
(pg. 
949
-
963
)
Stern
T.
Molnar
P.
Okaya
D.
Eberhart-Phillips
D.
Teleseismic P wave delays and modes of shortening the mantle lithosphere beneath South Island, New Zealand
J. geophys. Res. (Solid Earth)
 , 
2000
, vol. 
105
 
B9
(pg. 
21615
-
21631
)
Sutherland
R.
Davey
F.
Beavan
J.
Plate boundary deformation in South Island, New Zealand, is related to inherited lithospheric structure
Earth planet. Sci. Lett.
 , 
2000
, vol. 
177
 
3-4
(pg. 
141
-
151
)
Teanby
N.A.
Kendall
J.M.
Van der Baan
M.
Automation of shear-wave splitting measurements using cluster analysis
Bull. seism. Soc. Am.
 , 
2004b
, vol. 
94
 
2
(pg. 
453
-
463
)
Townend
J.
Sherburn
S.
Arnold
R.
Boese
C.
Woods
L.
Three-dimensional variations in present-day tectonic stress along the australia pacific plate boundary in New Zealand
Earth planet. Sci. Lett.
 , 
2012
, vol. 
353–354
 (pg. 
47
-
59
)
Vidale
J.E.
Complex polarization analysis of particle motion
Bull. seism. Soc. Am.
 , 
1986
, vol. 
76
 
5
(pg. 
1393
-
1405
)
Walcott
R.I.
Modes of oblique compression: late cenozoic tectonics of the South Island of New Zealand
Rev. Geophys.
 , 
1998
, vol. 
36
 
1
(pg. 
1
-
26
)
Walsh
E.
Savage
M.
Arnold
A.
Brenguier
F.
Rivemale
E.
Monitoring temporal changes with shear wave splitting: testing the methodology
Proceedings of the abstract S43F-2530 presented at 2012 Fall Meeting
 , 
2012
San Francisco, California
AGU
(pg. 
3
-
7
December
Wannamaker
P.E.
, et al.  . 
Fluid and deformation regime of an advancing subduction system at Marlborough, New Zealand
Nature
 , 
2009
, vol. 
460
 
7256
(pg. 
733
-
737
)
Wessel
A.
Automatic shear wave splitting measurements at Mt. Ruapehu volcano, New Zealand
 , 
2010
Awarded research masters thesis, School of Geography, Environment and Earth Sciences, Victoria University of Wellington
Wilson
C.K.
Jones
C.H.
Molnar
P.
Sheehan
A.F.
Boyd
O.S.
Distributed deformation in the lower crust and upper mantle beneath a continental strike-slip fault zone: Marlborough fault system, South Island, New Zealand
Geology
 , 
2004
, vol. 
32
 
10
(pg. 
837
-
840
)
Wirth
E.
Long
M.D.
Frequency-dependent shear wave splitting beneath the Japan and Izu-Bonin subduction zones
Phys. Earth planet. Inter.
 , 
2010
, vol. 
181
 
3-4
(pg. 
141
-
154
)
Wustefeld
A.
Bokelmann
G.
Zaroli
C.
Barruol
G.
Splitlab: a shear-wave splitting environment in MATLAB
Comput. Geosci.
 , 
2008
, vol. 
34
 
5
(pg. 
515
-
528
)

SUPPORTING INFORMATION

Additional Supporting Information may be found in the online version of this article.

Appendix A (http://gji.oxfordjournals.org/lookup/supp1/doi:10.1093/gji/ggt022/-/DC1)

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.