TY - JOUR
AU1 - Graham, Kenny M
AU2 - Savage, Martha K
AU3 - Arnold, Richard
AU4 - Zal, Hubert J
AU5 - Okada, Tomomi
AU6 - Iio, Yoshihisa
AU7 - Matsumoto, Satoshi
AB - SUMMARY Large earthquakes can diminish and redistribute stress, which can change the stress field in the Earth’s crust. Seismic anisotropy, measured through shear wave splitting (SWS), is often considered to be an indicator of stress in the crust because the closure of cracks due to differential stress leads to waves polarized parallel to the cracks travelling faster than in the orthogonal direction. We examine spatial and temporal variations in SWS measurements and the Vp/Vs ratio associated with the 2013 Cook Strait (Seddon, Grassmere) and 2016 Kaikōura earthquakes in New Zealand. These earthquake sequences provide a unique data set, where clusters of closely spaced earthquakes occurred. We use an automatic, objective splitting analysis algorithm and automatic local S-phase pickers to expedite the processing and to minimize observer bias. We present SWS and Vp/Vs measurements for over 40 000 crustal earthquakes across 36 stations spanning close to |$5\frac{1}{2}$| yr between 2013 and 2018. We obtain a total of 102 260 (out of 398 169) high-quality measurements. We observe significant spatial variations in the fast polarization orientation, ϕ. The orientation of gravitational stresses are consistent with most of the observed anisotropy. However, multiple mechanisms (such as structural, tectonic stresses and gravitational stresses) may control some of the observed crustal anisotropy in the study area. Systematic analysis of SWS parameters and Vp/Vs ratios revealed that apparent temporal variations are caused by variation in earthquake path through spatially varying media. Spatial analysis, Seismic anisotropy, Continental tectonics: strike-slip and transform, Dynamics: gravity and tectonics, Fractures, faults, and high strain deformation zones 1 INTRODUCTION Stress in the Earth is an important factor in earthquake genesis. Earthquakes are caused by the sudden rupture of rocks along faults exposed to differential stress in the Earth’s crust. Earthquakes occur when the accumulated shear stress exceeds the strength of faults or fractures in rock. Thus, crustal stresses are the immediate driving forces of earthquakes (Zoback & Zoback 2002). A variety of techniques have been devised to measure crustal stress. The most direct technique to determine crustal stress is through strain measurements in boreholes, where both stress magnitudes and orientations can be determined. However, drilling boreholes to seismogenic depths is a difficult and expensive enterprise (Townend & Zoback 2000; Zoback & Zoback 2002). Another well-established method is the inversion of earthquake focal mechanisms to determine the stress orientations in the region where the earthquake occurred (Michael 1984; Hardebeck & Michael 2006; Arnold & Townend 2007). Geodetic techniques such as using GNSS (Global Navigation Satellite System) and InSAR (Interferometric Synthetic Aperture Radar) measurements, which can be used to infer strain rates, can also be interpreted as changes in stress (Hardebeck & Michael 2006; Hardebeck & Okada 2018). In most regions, shear wave splitting (SWS) parameters can also serve to estimate the orientation of the principal horizontal stresses. Thus, stress orientations can be inferred from measuring crustal seismic anisotropy through SWS (Crampin 1984; Cochran & Kroll 2015; Savage et al. 2016) and changes in SWS parameters can also be interpreted as changes in the state of stress. Seismic anisotropy is the dependence of seismic velocity on direction, as caused by the elastic properties of rocks. The estimation of the geometry and strength of seismic anisotropy is often measured through SWS. SWS is a phenomenon that is observed when shear waves propagating through an anisotropic medium split into two nearly perpendicular phases that travel with different velocities. In the crust anisotropy is often considered to be controlled by either stress (stress induced anisotropy) or geological structures (structurally induced anisotropy) (Crampin 1984; Babuska & Cara 1991; Zinke 2000; Boness & Zoback 2006; Crampin et al. 2015; Savage et al. 2016). In the case of stress induced anisotropy, we assume that cracks or microcracks are randomly distributed. Differential stress will close cracks perpendicular to the maximum horizontal stress (SHmax), leaving open cracks only parallel to SHmax. In the case of structurally induced anisotropy, the fast wave is oriented parallel to structural fabric or the orientation of the residual features of palaeostress (Savage 1999; Cochran et al. 2006). SWS can be characterized by two parameters, the fast orientation (ϕ) and the delay time (δt). ϕ is the bearing of the polarization of the faster wave (also referred to as fast orientations) and δt is the time lag between the two split waveforms. δt depends on both the degree of anisotropy and the path length through the anisotropic medium. Although early studies concluded that stress-induced crack alignment is the principal cause of crustal anisotropy (Crampin 1984; Babuska & Cara 1991), several later SWS studies, close to and away from the vicinity of major faults, have concluded that other possible causes of crustal anisotropy exist (Matcham et al. 2000; Zinke 2000; Balfour et al. 2005; Boness & Zoback 2006). A consistent pattern where ϕ measurements in close proximity to faults show fault-parallel orientations and ϕ measurements further away from faults are parallel to the principal horizontal stress (orientations) has been reported around the San Andreas Fault of California, (Zinke 2000; Boness & Zoback 2006), Wellington region and Marlborough Fault System (MFS) of New Zealand (Gledhill 1991; Audoine et al. 2000; Matcham et al. 2000; Balfour et al. 2005; Karalliyadda & Savage 2013). Several researchers have used SWS to measure spatial and temporal variations in seismic anisotropy. Miller & Savage (2001) observed changes in fast polarization orientations before and after an eruption at Mt. Ruapehu volcano, New Zealand. They interpreted this as a variation in principal stress orientation, which they related to an increase in stress related to magmatic intrusion or inflation. Roman et al. (2011) observed rotations of fast orientations that correlated with rotating fault plane solutions for earthquakes associated with volcanic activity at Soufrière Hills volcano in Montserrat. Bianco et al. (2006) observed variations in fast orientation and delay times before the 2001 eruption on Mt. Etna. They observed delay times exhibiting a sudden decrease shortly before the start of the eruption and variations in the fast orientation five days before the onset of the eruption. Zheng et al. (2008) also studied the temporal variations of SWS in the aftershock region of 1999 Chichi Earthquake, Taiwan. They also observed a decrease in delay times shortly before the Chichi main shock and two Chiayi earthquakes. The relation between the ratio of body-wave (P and S waves) velocities, Vp/Vs, fluid content and pore pressure are well established (e.g. Wadati & Oki 1933; Nur 1971; Ito et al. 1979; Dvorkin et al. 1999, and references therein). The Vp/Vs ratio is one of the best parameters to indirectly identify fluid and crack migration within the crust, since pore fluid properties produce variations in seismic velocities (Dvorkin et al. 1999). In fluid-saturated rocks, high Vp/Vs anomalies may indicate high pore pressure (e.g. Nur 1971; Dvorkin et al. 1999). Variations in the Vp/Vs ratio have been interpreted qualitatively by several studies to be related to the movement of fluids in the crust. Li & Vidale (2001) observed an increase in shear velocity of the fault-zone rock around Johnson Valley, after the 1992 Mw7.5 Landers, California, earthquake (from 1994 to 1998). They interpreted the trend as an indication of fault healing by strengthening after the main shock, which is most likely due to the closure of cracks that opened during the 1992 earthquake. They also observed that the fault-zone strength recovery is consistent with a decrease in the apparent crack density within the fault zone. Moreover, they observed a decrease in the ratio of traveltime for P to S waves, which they explained as cracks near the fault zone being partially fluid-filled and becoming more fluid saturated with time (Li & Vidale 2001). Chiarabba et al. (2009), using S- and P-wave arrival times, observed an increase in Vp/Vs before the Mw 5.3 1997 Umbria-Marche earthquake in central Italy and attributed the trend to a pore pressure increase in fluid-filled cracks in the volume around the fault. The combination of variations in seismic anisotropy, Vp/Vs ratio and geodetic data has been reported by recent studies in volcanic areas and a slow-slip region (Savage et al. 2010b; Unglert et al. 2011; Zal et al.2020). Around the Wellington region of New Zealand, Gledhill (1991) and Matcham et al. (2000) used local earthquakes to study the anisotropic structure. Audoine et al. (2000) also characterized the crustal anisotropy around the transition between the subduction and the oblique transform faulting boundaries. In central MFS, they observed that fast polarization orientations were subparallel to the major faults and geological features, and attributed the observation to the presence of metamorphosed schist beneath the MFS. At the eastern edges of the MFS, they observed that fast orientations were parallel to the orientation of the maximum horizontal compressive stress, SHmax which is consistent with crack-induced anisotropy in the crust. In the North Island, their fast polarization orientations were parallel to the strike of the Hikurangi subduction zone as well as to the major geological features (Audoine et al. 2000). Balfour et al. (2005) described the upper plate stress regime above the Southern end of the Hikurangi subduction zone, where the dextral MFS accommodates more than 80 per cent of the relative plate motion (Townend et al. 2012). They examined the relationship between the SHmax orientations and the average fast orientations, and concluded that both stress and structural anisotropy exist around the MFS, but the dominant structural fabric of the MFS controls the seismic anisotropy measurement rather than SHmax (Balfour et al. 2005). Regional stress studies by Townend et al. (2012) and Balfour et al. (2005) around the study area revealed a uniform SHmax of 115° ± 16°, interpreted as the upper plate extending across the northern South Island. Recently, Evanzia et al. (2017) tried to distinguish between stresses around Southern Hikurangi Margin, New Zealand, using SWS and focal mechanism measurements. They suggested that stresses in the overriding plate are likely related to bending stresses, gravitational stresses and tectonic loading. In this study, we perform a systematic analysis of SWS parameters to investigate the physical mechanisms causing crustal anisotropy around the Marlborough and Wellington region, and to explore possible temporal variations in SWS measurements. The 2013 Cook Strait and 2016 Kaikōura earthquake sequences provide a rare data set, where clusters of closely spaced earthquakes with full azimuthal coverage occurred. We use this data set to examine the hypotheses of structure and stress-induced anisotropy that are often used to explain anisotropy in the crust. The stress-induced hypothesis suggests that cracks aligned with the maximum horizontal compressive stress direction, SHmax, remain open with a preferred orientation (Crampin & Booth 1989; Savage 1999) and thus induce anisotropy. The stresses acting on crustal tectonic settings can be characterized into two categories; (i) topographical stresses acting on the shallow crust (often referred to as gravitational stresses) and (ii) stresses driven by the relative plate motions or geodynamic processes (referred to as tectonic stresses). Often, crustal anisotropy studies use SHmax from focal mechanism inversions to explain the stress hypothesis with little focus on SHmax from gravitational stresses. Araragi et al. (2015) and Illsley-Kemp et al. (2017) are two of the few studies that suggested gravitational stress as a possible explanation for crustal seismic anisotropy around Mt. Fuji (Japan) and the Northern Afar region (East Africa) respectively. In theory, the orientation of gravitationally induced horizontal stresses, |$S_{H_{\mathrm{ max}}}^{\mathrm{ Grav}}$|⁠, have different orientations at the peak of a mountain, as opposed to the slopes either side. At the peak |$S_{H_{\mathrm{ max}}}^{\mathrm{ Grav}}$| is parallel to the strike of the mountain ridge, whereas, on the slope it is perpendicular (Flesch et al. 2001; Hirschberg et al. 2018). In basins of two mountains ridges, the expected |$S_{H_{\mathrm{ max}}}^{\mathrm{ Grav}}$| are often perpendicular to the strike of the ridges due to the compressional stresses on the slopes on either side of the basin (Hirschberg et al. 2018). Here, we examine how well both the gravitational and tectonic stresses explain the crustal anisotropy in central New Zealand (southern North Island and northwestern south island; see Fig. 1) by comparing the splitting measurements with SHmax from both gravitational and tectonic stresses (stresses at depth from focal mechanism inversion). Figure 1. Open in new tabDownload slide Map showing the tectonic setting of the study area with the stations used (red and blue triangles) plotted on a basemap of the digital topography and bathymetry (Mitchell et al. 2012). GeoNet CMT solutions of Mw > 6 earthquakes from 2013 to 2018 are shown; blue focal mechanism plots (beach balls) represent reverse mechanisms, with the red representing normal faulting with some dip slip motion and green represents a strike-slip mechanism. Major faults marked include the Wairau (WF), Awatere (AF), Clarence (CF) and Hope (HF) faults. The red lines denote the surface rupture of faults during Kaikōura main shock (Langridge et al. 2016). The yellow star indicates the epicenter of the Mw 7.8 Kaikōura earthquake. Inset: the study region (red frame): the blue line is the plate boundary fault showing the Hikurangi and Puysegur trenches. The structurally controlled hypothesis suggests that the orientation of the anisotropic medium is associated with the orientation of the geological structures (such as fault orientations, shear planes subparallel to faults) (Zinke 2000; Boness & Zoback 2006; Cochran et al. 2006; Okaya et al. 2016). Mineral alignment, associated with shearing in fault zones, or foliation in metamorphosed rocks has been one of the hypotheses used to explain structurally controlled anisotropy in the crust. Around central New Zealand, Audoine et al. (2000) attributed the observed anisotropy to the metamorphosed schist in the crust because recrystallization during metamorphism may form schistosity or a pervasive foliation. Okaya et al. (1995) also observed a high degree (up to 20 per cent) of anisotropy associated with the Haast Schist of the South Island, New Zealand. Okaya et al. (2016) explained the crustal anisotropy in the Taiwan collision zone by deformation-related (due to pervasive shearing in strike-slip fault zones) mineral preferred orientation in the metamorphic rocks (e.g. schist fabric). In this paper, we examine these various hypotheses by comparing our splitting measurements with active fault orientations, orientations of horizontal principal stress (both tectonic and gravitational) and maximum shear plane orientations from the strain rate field (related to crustal deformation). 2 TECTONIC SETTING AND SEISMICITY New Zealand’s tectonic setting is characterized by two subduction systems of opposite polarity, connected by an area of oblique continental convergence (inset in Fig. 1). Westward subduction (and intra-arc rifting) in the North Island at the Hikurangi subduction zone changes to strike-slip at the Marlborough Fault Zone, transpressional collision along the Alpine Fault and finally back to (eastward) subduction at the Puysegur Trench (inset in Fig. 1) offshore (Mortimer 2004). The study area (red box on inset in Fig. 1), central New Zealand, covers the southern end of the Hikurangi subduction zone and the north eastern part of New Zealand’s South Island, which is dominated by the MFS. The MFS marks the transition from a subduction plate boundary in the north to a transpressive plate boundary (the Alpine Fault) in the south. The faults in the MFS are predominantly strike-slip (trending NE-SW parallel to the strike of the Pacific-Australia plate motion), with a relatively small component of reverse motion (Van Dissen & Yeats 1991). The average regional strike of the faults is 55° and their near-surface dips vary from 60° to near-vertical (Van Dissen & Yeats 1991; Anderson et al. 1993). Major geological structures around the southern end of the Hikurangi subduction zone strike in the NE-SW direction with active faulting also following this NE-SW trend (Mortimer 2004; Litchfield et al. 2017). Southwest from Kaikōura lies the North Canterbury fold and thrust belt (Reyners & Cowan 1993; Pettinga et al. 2001). The NE-SW trend of the thrust fault extends through the northeastern part of the Canterbury region. In response to the transition from the subduction related tectonics in the north, thrust faults are evolving (Pettinga et al. 2001). The thrust faults are expressed as topographical ridges with NE-SW strike extending close to the Hope fault (Reyners & Cowan 1993; Pettinga et al. 2001). Nearly two-thirds of the Marlborough region is underlain by greywacke of the Torlesse Supergroup, a laterally uniform wave-speed structure (Reyners & Cowan 1993; Mortimer 2004). The Marlborough Schist found north of the Wairau fault, which is part of the Haast Schist of South Island, is commonly considered to be strongly metamorphosed (Mortimer 1993). The 2013 Cook Strait earthquake sequence started in 2013 July with two foreshocks of Mw 5.7 and 5.8 and climaxed in the 2013 July 21, Mw 6.5 Seddon earthquake and the 2013 August 16, Mw 6.6 Lake Grassmere earthquake (Fig. 1). These large earthquakes, located ∼50 km south of Wellington, New Zealand’s capital, generated significant ground shaking in the Wellington and Marlborough regions (Holden et al. 2013; Hamling et al. 2014). Both earthquakes were strike-slip events with similar magnitudes and characteristics and they were considered to be a ‘doublet’ (Holden et al. 2013). The seismicity in the Cook Strait sequence region is still above the background levels that existed prior to 2013, with more than 16 000 earthquakes (ML ≥ 1) within a 4 yr period from 2013 January to 2017 November (Holden et al. 2013; GeoNet 2019, accessed 2019 February 2). Seismicity in Marlborough (before the 2016 Kaikōura earthquake) was mainly concentrated in the region above the subduction interface and around the north-eastern part of the MFS (Holden et al. 2013). The area was the location of a swarm of M > 4 earthquakes in 2005 and before that, it was the location of the 1977 M 6 Cape Campbell earthquake and the 1966 M 5.8 Seddon earthquake (Holden et al. 2013). Six months after the Mw 6.6 Lake Grassmere earthquake, the 2014 Mw 6.2 Eketahuna earthquake struck 10 km east of Eketahuna. The depth (34 km) and location of the event placed it below the interface between the subducting Pacific plate and overriding Australian plate (Abercrombie et al. 2016). The afterhocks of the 2014 Eketahuna earthquake are characterized by repeating events (Wallace et al. 2017). The Mw 7.8, 2016 November 14, Kaikōura earthquake, is the largest and most complex earthquake recorded on land in New Zealand (Clark et al. 2017; Hamling et al. 2017). The earthquake initiated near the North Canterbury region at a depth of 15 km with an oblique thrust faulting mechanism (Hamling et al. 2017, fig. 1). The rupture propagated from south-west towards the north-east for about 120 s, with an unusual source process, starting with weak radiation and releasing more energy while propagating towards the north-east over a distance of ∼150 km and terminating offshore in the Cook Strait region (Cesca et al. 2017; Duputel & Rivera 2017; Holden et al. 2017; Hamling et al. 2017). A remarkable number (>20) of shallow crustal fault segments (involving a combination of reverse and dextral strike-slip faulting; Fig. 1) ruptured, including vertical motions of more than 10 m and horizontal displacements over 11 m (Clark et al. 2017; Cesca et al. 2017; Hamling et al. 2017). Most of the faults that ruptured have a general NE-SW trend (red lines in Fig. 1), but some have NW-SE orientation, revealing a complex and heterogeneous slip pattern (Clark et al. 2017; Hamling et al. 2017; Kaiser et al. 2017). Surface rupture was mostly associated with known onshore faults. However, some surface traces were produced by faults that had not been previously mapped. Also some faults ruptured offshore, causing a small tsunami (Clark et al. 2017; Wallace et al. 2017). The 2016 Kaikōura earthquake was followed by more than 25,000 local and regional aftershocks (GeoNet catalogue), clustered in three unique spatial patterns (Kaiser et al. 2017). More than 100 000 landslides were triggered by the earthquake and subsequent aftershocks, with 50 of them yielding significant landslide dams (Litchfield et al. 2017; Dellow et al. 2017). Another unusual aspect of the Kaikōura earthquake was the occurrence of large-scale (>15 000 km2) slow-slip events triggered on the central and northern Hikurangi subduction interface (250 and 600 km away from the epicentre, respectively) due to dynamic-stress changes from passing seismic waves (Wallace et al. 2017). 3 DATA AND METHODOLOGY 3.1 Data We determine SWS measurements for over 40 000 local crustal earthquakes that were located around the region of the 2013 Cook Strait earthquake sequence (Fig. 2). Our earthquake catalogue consists of GeoNet detection and locations (Petersen et al. 2011) and spans more than 5 yr: from 2013 January to 2018 June. Although Lanza et al. (2019) relocated some of the events used here, we used the GeoNet catalogue to ensure consistency. The GeoNet catalogue used here has a good azimuthal coverage of events before and after the 2016 Kaikōura main shock (see Fig. S1, Supporting Information). This enables us to search for temporal variations in SWS parameters and Vp/Vs ratios. Figure 2. Open in new tabDownload slide Map of earthquake epicentres used for analysis. Red triangles are the locations of GeoNet stations used. Events are colour-coded by the hypocentral depths. Blue focal mechanism solution (beach ball) represents the Kaikōura earthquake faulting mechanism; red beach ball represents the Eketahuna earthquake faulting mechanism, and green beach balls represent the two Cook Strait sequence earthquake faulting mechanisms (Ristau 2013). Focal mechanisms solutions were obtained from the GeoNet regional moment tensor solution catalog. A cross-section from A to A1 shows the projected depth distributions, coloured by time (timescale on main figure). We use 36 stations deployed around the Wellington and Marlborough region. These stations include 24 permanent seismic stations (combination of broad-band and short-period instruments; red triangles on Fig. 1) operated by GeoNet (the Geological hazard information for New Zealand) and 12 temporary short-period stations installed and operated by DPRI (Disaster Prevention Research Institute, Kyoto University, Japan) (Okada et al. 2019) at different time periods between 2009 and 2018 (blue triangles on Fig. 1). Plots of data continuity for all stations are shown in Supporting Information S2. All the stations were re-sampled with a common sampling frequency of 100 Hz. Since the objective was to study crustal anisotropy, the hypocentral depths were limited to less than 35 km to include only earthquakes whose ray paths travel wholly through the crust. Moreover, we limited events magnitudes to Mw 1.5 and above to remove any poor-quality waveforms. 3.2 Method To estimate SWS parameters for thousands of waveforms, a fully automated and systematic technique was implemented. Here, we used the automatic splitting analysis code. Multiple Filter Automatic Splitting Technique (MFAST, Savage et al. 2010a; Wessel 2010), which is designed to handle large volumes of data. MFAST is based on the eigenvalue-minimization method of Silver & Chan (1991) and the clustering method of Teanby et al. (2004a). MFAST uses an automated workflow to estimate the splitting parameters with an objective grading of measurements (Walsh et al. 2013). A concise description for MFAST is presented in Section 3.2.2 and a detailed description of the method is presented by Savage et al. (2010a), Wessel (2010) and Walsh (2012). We used P-wave arrival times and origin times, to determined by GeoNet (Petersen et al. 2011). We used an automatic algorithm (Section 3.2.1) to estimate local shear wave arrival times for seismographs with good signal-to-noise ratio (SNR > 3). All SWS measurements were made within the ‘shear wave window’ (incidence angles less than a chosen critical angle, defined below). We used 1-D synthetic analysis (determining the effect of angle of incidence on splitting measurements) to resolve the ultimate critical angle for our analysis. This was necessary because, outside the shear wave window, the shear waveforms are susceptible to S to P conversions and scattering at the surface causing nonlinear particle motion (Nuttli 1961; Crampin & Gao 2006; Savage et al. 2016). Neuberg & Pointer (2000) showed that waveforms outside the shear wave window generate elliptical particle motion even without the presence of anisotropy, especially when recording shallow local earthquakes in the vicinity of strong topography. We simulated waveforms in a 1-D anisotropic medium using the Levin & Park (1997) 1-D reflectivity code (which was subsequently modified by Castellazzi et al. 2015; Walsh 2012) and estimated the apparent splitting parameters using MFAST (Wessel 2010; Savage et al. 2010a). Fig. S3 in the Supporting Information shows the 1-D model used for analysis. The effect of angle of incidence on splitting analysis was tested by varying angle of incidence with backazimuth. As shown in Fig. S4 in the Supporting Information, we observed a strong backazimuthal variation of δt and ϕ as angle of incidence is increased (here vertical incidence is 0°). At higher angle of incidence > 30°, we observe significant variations of splitting parameters with backazimuth. Since we observed less variations of splitting parameters with backazimuth at an incidence angle of 35° compared to higher incidence angle, we selected a maximum angle of 35° as the best shear wave window for our analysis. We estimated angles of incidence using the TauP toolkit (Crotwell et al. 1999) and a 1-D velocity model extracted from the Eberhart-Phillips & Fry (2018) 3-D velocity model. 3.2.1 Local phase arrival picking Determination of S-wave arrival times is necessary to estimate SWS measurements. Due to the large volume of our data set, manually picking the S-phase arrival times would have been too time consuming. We therefore used two automatic picking techniques. First, local S-phase arrival times (S picks) for GeoNet stations were automatically picked using the technique of Diehl et al. (2009) modified by Castellazzi et al. (2015), hereafter called Spicker. For stable and reliable S-wave picking, the Spicker technique combines three different detection and picking methods. The STA/LTA (Short Term Averaging / Long Term Averaging) detector (e.g. Allen 1978) and polarization detector (e.g. Flinn 1965) are used to identify the first arriving S phase. The information provided by the detectors is then used to set up the search windows of the autoregressive picker using the Akaike Information Criterion (AR-AIC) as outlined by Takanami & Kitagawa (1988). Finally, all three methods are combined to yield the best arrival time of the first arriving S phase and its corresponding uncertainty. The Spicker makes use of all the three component seismogram to estimate the S-phase arrival. See Diehl et al. (2009) and Castellazzi et al. (2015) for further details on the algorithm. We calibrated the picking algorithm with over 3500 manually picked S phases. These S phases were randomly selected to obtain a uniform station, magnitude and depth distribution. The time difference between manual and Spicker picks revealed a symmetrical unimodal distribution with a mean of 0.030 s and a standard deviation of 0.598 s. In addition, 90 per cent percent of the Spicker picks were within ±0.484 s of the manual Spicks (Fig. 3a). This picking accuracy is adequate for our studies, because MFAST uses multiple windows for its analyses. The best parameters (Table S1, Supporting Information) that yielded more picks and small time difference between manual and Spicker picks were used for the final picking. We compared the common S picks (∼21 000) between GeoNet manual S and Spicker picks across all stations. The time difference between GeoNet manual S and Spicker picks revealed a symmetric unimodal distribution with mean 0.028, standard deviation 0.648 and 90 per cent of the picks were within ±0.678 (Fig. 3b). Over 200 000 S phases were obtained, representing about 23 per cent of the expected (event-station pair) arrivals. To assess the quality of the S picks across our data set, we randomly selected approximately 5 per cent of the S picks across all stations. The selected events were manually inspected and approximately 15 per cent of S picks were unreliable and were either revised or discarded. Based on this analysis we estimate approximately 85 per cent of the S picks are reliable. Figure 3. Open in new tabDownload slide Distributions of local S-phase picking errors relative to analyst pick. (a) Time difference between Spicker picks and manually picked S-phase arrival times. (b) Time difference between GeoNet manual and Spicker picks. Some of components of the short-period stations operated by DPRI were intermittently missing components. We could not apply the Spicker technique to estimate the S-phase arrival, since the technique required all three components for the estimation. Instead we used the Generalized Seismic Phase Detection (GPD), technique developed by Ross et al. (2018) for estimating phase arrivals for the DPRI stations. This technique trained a convolutional neural network (a deep learning framework) to learn the features of seismic waveforms. An already trained model can be applied to waveforms that are not included in the training set, thus a trained model can be applied to different tectonic settings (Kong et al. 2018; Ross et al. 2018). We used the Ross et al. (2018) GPD-framework model which they trained and validated with a total of 4.5 million three-component waveforms across a range of magnitudes. With each component’s time-series, the model predicts the probability that a phase arrival is either a P or Sarrival. We set the minimum probability for phase detection to 95 per cent. An example of a pick example is shown in Fig. S5 in the Supporting Information. To gain confidence in the GDP picks, we estimated P- and S-arrival times for events with GeoNet manual S picks (more than 40 000 picks across the 24 permanent GeoNet stations) and compared the time difference between the GDP and GeoNet manual S picks. The time difference shows a concentrated symmetric uni-modal distribution with mean −0.072 s, standard deviation 0.774 s and 90 per cent percent of the picks were within ±0.261 s (see Fig. S6, Supporting Information). 3.2.2 Measuring splitting parameters SWS parameters were estimated using the MFAST code (Savage et al. 2010a; Wessel 2010). This procedure finds the inverse splitting operator that best removes the splitting. The waveforms are filtered with a series of bandpass filters, and the product of the bandwidth and the SNR of the filtered waveform is used to determine the best three filters. Splitting measurements for all three filtered waveforms are estimated using the Silver & Chan (1991) and Teanby et al. (2004b) method with modification by Walsh et al. (2013). The splitting parameters are estimated using the eigenvalue minimization method repeated over multiple time windows around the S phase. The splitting parameter pair (ϕ, δt) that best removes the splitting as measured by the smallest eigenvalue of the corrected covariance matrix, is chosen as the best measurement for the given time window (Savage et al.2010a). This procedure is then repeated for 75 windows covering slightly different time spans around the S phase. The windows are automatically selected based on the dominant period around the filtered shear wave (Savage et al. 2010a). Cluster analysis over the 75 window measurements is used to select the final splitting parameter for the filter under consideration and to calculate the uncertainty associated with the measurements (Teanby et al. 2004a). Measurements are graded from A to D depending on the consistency between the (ϕ, δt) and measurements in the different windows. MFAST also provides a measure of the incoming polarization orientation (ϕin) by estimating the eigenvalues of the corrected components after splitting is removed (Silver & Chan 1991; Savage et al. 2010a). ϕin corresponds to the polarization orientation of the shear wave before it enters the anisotropic layer responsible for the measured splitting parameters. As a quality control, we removed all null graded measurements (events for which no measurable splitting occurs or event for which ϕ is within 20° of parallel or perpendicular to ϕin) and also kept only A and B graded measurements. We also limited maximum delay times, δt to 0.4 s since local crustal events are often not expected to have larger delay times (Balfour et al. 2005; Cochran & Kroll 2015). Following the quality control, a total of 102 260 (out of 398 169) high-quality measurements were obtained. Wessel (2010), Savage et al. (2010a) and Walsh et al. (2013) give a detailed description of the MFAST method. We estimated Vp/Vs ratios along the ray path of each event that yielded high quality SWS measurements to each station using the approach of Wadati & Oki (1933) and Nur (1971). We assumed a linear ray path and that the Vp/Vs ratio is homogeneous along the seismic wave’s path (Kisslinger & Engdahl 1973). The Vp/Vs ratio of each event-station pair was calculated using Vp/Vs = (ts − to)/(tp − to), where ts and tp are the arrival times of the S and P waves, respectively and to is the event origin time (from the GeoNet catalogue). Only Vp/Vs values between 1.5 and 2.3 were used because values outside of this range are not expected around the study region (Eberhart-Phillips & Fry 2018) and would probably be due to inaccurate arrival time picks. We also estimated the percent anisotropy, κ, following the approach of Babuska & Cara (1991) and Savage (1999). Here, we define |$\kappa = (v_{\mathrm{ max}}-v_{\mathrm{ min}})/ \bar{v}$|⁠. Where vmax is velocity of the shear wave along the fast orientation, vmin is the velocity of the shear wave along the slow orientation and |$\bar{v}$| is an average of vmax and vmin. Assuming both fast and slow waves travel with the same path length, κ can be related to δt as: κ = 200*δt/(2ts + δt). We use κ as well as δt, to characterize the anisotropy structure. 3.2.3 Spatial averaging SWS measurements provide an estimate of the anisotropy along the ray propagation path between the source and station. However, it can be difficult to determine exactly where along the ray propagation path the anisotropy originates. With measurements from several stations and a dense cluster of events, one can probe the spatial variation of ϕ and δt using a spatial averaging technique. We used the TESSA (Tomography Estimate and Shear wave splitting Spatial Average) technique by Johnson et al. (2011) to estimate the spatial averages of ϕ. For the spatial averaging technique, the study area is divided into cells or blocks using the recursive quadtree clustering algorithm described by Townend & Zoback (2001). We set the minimum block size to 5 km2 and the minimum and maximum number of ray paths to be 20 and 80, respectively. The quadtree gridding algorithm works in an iterative process. Cells with less than 20 rays passing through are not used for analysis and those with more than 80 rays passing through are subdivided until a minimum of block size of 5 km2 or until a minimum number of 20 events is attained. These criteria were chosen to ensure that each block contained enough data to give a reliable average measurement. To take into account variations due to heterogeneous anisotropic structure along the ray path (e.g. Rümpker & Silver 1998; Johnson et al. 2011), rays in each grid cell are weighted by the inverse of the distance squared (1/d2) (where d is the distance from the station to the grid cell in question). This weighting scheme is used because we expect splitting to occur later in the ray path and we did not observe a strong correlation between δt and hypocentral depth (Johnson et al. 2011). For comparisons with time, we used the regular gridding to ensure that the data points are always at the same locations. The spatial averages of ϕ are computed using a circular statistics approach (Berens 2009) and are estimated only when the standard deviation of fast orientations in each grid is less than 30° and the standard error of the mean is less than 10°. 3.2.4 Quantitative comparison of ϕ measurements We made quantitative comparisons of our spatially averaged ϕ measurement with GNSS derived principal contraction axes, average fault orientations, gravitational stress and focal mechanism inversion measurements. We also quantitatively compared averaged ϕ measurements between events before and after the 2016 Kaikōura earthquake to observe any spatiotemporal variation associated with the 2016 Kaikōura main shock. For our quantitative comparison, we estimate a test statistic, F, the absolute value of the cosine of the difference between the two angles (F = |cos (ψ1 − ψ2)|) where ψ1 is the mean ϕ measurement and ψ2 is the comparison angle. Since we do not have an identical set of locations at which mean ϕ has been estimated, we use locations that are closest together, and < 10 km to estimate the residual. Values of F range from 0 to 1, where 1 represents parallel orientation and 0 represents perpendicular orientation between ψ1 and ψ2. We contour the residual using the GMT functions grdsample and grdimage (Wessel et al. 2013). We also estimate the mean and median of the F to give an indication of how well the two angles are correlated. 3.2.5 Averaging for time variations To examine variations we employed an averaging technique for values of δt to smooth and minimize the scatter as much as possible. We use a simple moving median as the smoothing function that reduces the noise in the measurement. We used a 20-d moving median filter on the individual results with the aim of removing the scattered values from the data set. We estimate the 95 per cent confidence intervals by using the non-parametric bootstrap approach. In this approach, we randomly sample with replacement from the original sample a data set of the same size as the original sample. We repeat this for 1000 replicates and for each replicate data set we estimate the median value. The 2.5th and 97.5th percentile of the replicate data set is estimated as the respective lower and upper values of the 95 per cent confidence interval for the median. 4 RESULTS We determined SWS measurements for more than 40 000 waveforms recorded on at least one of the 36 stations within the study region (Fig. 1). Tables 1 and 2 give a summary of fast orientations, delay times, percentage of anisotropy and Vp/Vs ratio as well as associated descriptive statistics for each station. Out of ∼102 000 high-quality SWS parameters we estimated, the delay times vary between near zero to 0.372 s with an average of 0.157±0.001 s. The Vp/Vs ratio varies from 1.58 to 2.22 with an average of 1.741±0.001. Percent anisotropy varies from 0 to 5.186 per cent with an average of 0.922±0.004 per cent. Table 1. Summary statistics of fast orientation (ϕ), delay time (δt), Vp/Vs ratio and percent anisotropy (per cent Aniso) with their 95 per cent confidence interval for GeoNet permanent stations Station . Events . Mean ϕ(°) . SD ϕ . Mean δt(s) . SD δt . Mean Vp/Vs . SD Vp/Vs . Mean per cent Aniso . Per cent Aniso SD . BHW 3515 81.14 ±7.40 62.82 0.136±0.002 0.066 1.724±0.002 0.066 0.759±0.015 0.424 BSWZ 4781 52.82 ±1.76 43.38 0.181±0.002 0.076 1.743±0.002 0.084 1.453±0.023 0.801 CAW 6499 17.05 ±3.29 55.91 0.171±0.002 0.081 1.744±0.001 0.051 0.683±0.009 0.342 CMWZ 2162 -82.10±10.43 64.11 0.173±0.003 0.072 1.819±0.008 0.181 2.371±0.065 1.534 CRSZ 1588 67.47 ±4.07 48.34 0.168±0.004 0.085 1.808±0.014 0.281 2.118±0.060 1.167 DUWZ 2639 28.60 ±3.68 50.79 0.203±0.003 0.086 1.738±0.002 0.044 0.623±0.011 0.262 HOWZ 2015 69.93 ±6.58 57.45 0.139±0.003 0.075 1.732±0.002 0.045 0.673±0.019 0.419 KHZ 6343 -32.95±2.89 53.83 0.149±0.002 0.073 1.736±0.001 0.055 0.713±0.011 0.421 KIW 5497 53.85 ±3.11 53.87 0.152±0.002 0.073 1.718±0.001 0.044 0.629±0.008 0.296 MRNZ 2375 52.98 ±3.78 50.37 0.205±0.003 0.081 1.741±0.001 0.030 0.535±0.009 0.214 MSWZ 3481 60.46 ±2.02 43.02 0.142±0.003 0.077 1.784±0.002 0.058 0.560±0.011 0.327 MTW 7300 61.60 ±3.03 55.57 0.150±0.002 0.071 1.771±0.001 0.039 0.625±0.008 0.330 NNZ 1832 50.30 ±10.84 63.55 0.206±0.004 0.082 1.729±0.002 0.040 0.684±0.015 0.298 OGWZ 4232 37.96 ±1.16 34.13 0.123±0.002 0.067 1.736±0.001 0.047 0.478±0.009 0.282 PAWZ 1908 75.94 ±7.12 58.17 0.188±0.004 0.078 1.774±0.002 0.052 0.691±0.015 0.329 PLWZ 1460 -22.50±6.66 55.31 0.166±0.004 0.077 1.780±0.003 0.050 0.700±0.019 0.351 QRZ 1632 46.28 ±4.34 49.59 0.205±0.004 0.077 1.733±0.001 0.025 0.395±0.008 0.156 TCW 9199 14.07 ±2.03 51.22 0.133±0.001 0.069 1.715±0.001 0.054 0.643±0.008 0.362 THZ 2676 26.09 ±4.69 54.62 0.163±0.003 0.078 1.695±0.002 0.044 0.611±0.011 0.281 TKNZ 1798 86.31 ±3.14 44.95 0.226±0.003 0.065 1.735±0.002 0.034 0.543±0.008 0.165 TMWZ 809 22.97 ±19.73 65.95 0.110±0.005 0.064 1.787±0.002 0.030 0.567±0.025 0.334 TRWZ 629 42.11 ±5.33 45.03 0.189±0.006 0.074 1.789±0.006 0.069 0.786±0.033 0.400 TUWZ 4648 51.77 ±1.37 38.43 0.165±0.002 0.074 1.719±0.001 0.045 0.957±0.016 0.528 WEL 4144 -81.07±2.53 48.37 0.144±0.003 0.082 1.747±0.002 0.065 0.711±0.013 0.414 Total 83183 47.59±1.05 57.85 0.158±0.001 0.077 1.744±0.001 0.079 0.847±0.004 0.645 Station . Events . Mean ϕ(°) . SD ϕ . Mean δt(s) . SD δt . Mean Vp/Vs . SD Vp/Vs . Mean per cent Aniso . Per cent Aniso SD . BHW 3515 81.14 ±7.40 62.82 0.136±0.002 0.066 1.724±0.002 0.066 0.759±0.015 0.424 BSWZ 4781 52.82 ±1.76 43.38 0.181±0.002 0.076 1.743±0.002 0.084 1.453±0.023 0.801 CAW 6499 17.05 ±3.29 55.91 0.171±0.002 0.081 1.744±0.001 0.051 0.683±0.009 0.342 CMWZ 2162 -82.10±10.43 64.11 0.173±0.003 0.072 1.819±0.008 0.181 2.371±0.065 1.534 CRSZ 1588 67.47 ±4.07 48.34 0.168±0.004 0.085 1.808±0.014 0.281 2.118±0.060 1.167 DUWZ 2639 28.60 ±3.68 50.79 0.203±0.003 0.086 1.738±0.002 0.044 0.623±0.011 0.262 HOWZ 2015 69.93 ±6.58 57.45 0.139±0.003 0.075 1.732±0.002 0.045 0.673±0.019 0.419 KHZ 6343 -32.95±2.89 53.83 0.149±0.002 0.073 1.736±0.001 0.055 0.713±0.011 0.421 KIW 5497 53.85 ±3.11 53.87 0.152±0.002 0.073 1.718±0.001 0.044 0.629±0.008 0.296 MRNZ 2375 52.98 ±3.78 50.37 0.205±0.003 0.081 1.741±0.001 0.030 0.535±0.009 0.214 MSWZ 3481 60.46 ±2.02 43.02 0.142±0.003 0.077 1.784±0.002 0.058 0.560±0.011 0.327 MTW 7300 61.60 ±3.03 55.57 0.150±0.002 0.071 1.771±0.001 0.039 0.625±0.008 0.330 NNZ 1832 50.30 ±10.84 63.55 0.206±0.004 0.082 1.729±0.002 0.040 0.684±0.015 0.298 OGWZ 4232 37.96 ±1.16 34.13 0.123±0.002 0.067 1.736±0.001 0.047 0.478±0.009 0.282 PAWZ 1908 75.94 ±7.12 58.17 0.188±0.004 0.078 1.774±0.002 0.052 0.691±0.015 0.329 PLWZ 1460 -22.50±6.66 55.31 0.166±0.004 0.077 1.780±0.003 0.050 0.700±0.019 0.351 QRZ 1632 46.28 ±4.34 49.59 0.205±0.004 0.077 1.733±0.001 0.025 0.395±0.008 0.156 TCW 9199 14.07 ±2.03 51.22 0.133±0.001 0.069 1.715±0.001 0.054 0.643±0.008 0.362 THZ 2676 26.09 ±4.69 54.62 0.163±0.003 0.078 1.695±0.002 0.044 0.611±0.011 0.281 TKNZ 1798 86.31 ±3.14 44.95 0.226±0.003 0.065 1.735±0.002 0.034 0.543±0.008 0.165 TMWZ 809 22.97 ±19.73 65.95 0.110±0.005 0.064 1.787±0.002 0.030 0.567±0.025 0.334 TRWZ 629 42.11 ±5.33 45.03 0.189±0.006 0.074 1.789±0.006 0.069 0.786±0.033 0.400 TUWZ 4648 51.77 ±1.37 38.43 0.165±0.002 0.074 1.719±0.001 0.045 0.957±0.016 0.528 WEL 4144 -81.07±2.53 48.37 0.144±0.003 0.082 1.747±0.002 0.065 0.711±0.013 0.414 Total 83183 47.59±1.05 57.85 0.158±0.001 0.077 1.744±0.001 0.079 0.847±0.004 0.645 Notes: Only high-quality measurements (described in Section 3.2.2) are presented in the summary table. ± represents the 95 per cent confidence intervals. SD shows the standard deviation of the measurements. Station averages and standard deviation for fast orientations were performed using circular statistics. Event represents the number of event–station pair used for analysis. Total shows the summary of all stations for GeoNet stations. Open in new tab Table 1. Summary statistics of fast orientation (ϕ), delay time (δt), Vp/Vs ratio and percent anisotropy (per cent Aniso) with their 95 per cent confidence interval for GeoNet permanent stations Station . Events . Mean ϕ(°) . SD ϕ . Mean δt(s) . SD δt . Mean Vp/Vs . SD Vp/Vs . Mean per cent Aniso . Per cent Aniso SD . BHW 3515 81.14 ±7.40 62.82 0.136±0.002 0.066 1.724±0.002 0.066 0.759±0.015 0.424 BSWZ 4781 52.82 ±1.76 43.38 0.181±0.002 0.076 1.743±0.002 0.084 1.453±0.023 0.801 CAW 6499 17.05 ±3.29 55.91 0.171±0.002 0.081 1.744±0.001 0.051 0.683±0.009 0.342 CMWZ 2162 -82.10±10.43 64.11 0.173±0.003 0.072 1.819±0.008 0.181 2.371±0.065 1.534 CRSZ 1588 67.47 ±4.07 48.34 0.168±0.004 0.085 1.808±0.014 0.281 2.118±0.060 1.167 DUWZ 2639 28.60 ±3.68 50.79 0.203±0.003 0.086 1.738±0.002 0.044 0.623±0.011 0.262 HOWZ 2015 69.93 ±6.58 57.45 0.139±0.003 0.075 1.732±0.002 0.045 0.673±0.019 0.419 KHZ 6343 -32.95±2.89 53.83 0.149±0.002 0.073 1.736±0.001 0.055 0.713±0.011 0.421 KIW 5497 53.85 ±3.11 53.87 0.152±0.002 0.073 1.718±0.001 0.044 0.629±0.008 0.296 MRNZ 2375 52.98 ±3.78 50.37 0.205±0.003 0.081 1.741±0.001 0.030 0.535±0.009 0.214 MSWZ 3481 60.46 ±2.02 43.02 0.142±0.003 0.077 1.784±0.002 0.058 0.560±0.011 0.327 MTW 7300 61.60 ±3.03 55.57 0.150±0.002 0.071 1.771±0.001 0.039 0.625±0.008 0.330 NNZ 1832 50.30 ±10.84 63.55 0.206±0.004 0.082 1.729±0.002 0.040 0.684±0.015 0.298 OGWZ 4232 37.96 ±1.16 34.13 0.123±0.002 0.067 1.736±0.001 0.047 0.478±0.009 0.282 PAWZ 1908 75.94 ±7.12 58.17 0.188±0.004 0.078 1.774±0.002 0.052 0.691±0.015 0.329 PLWZ 1460 -22.50±6.66 55.31 0.166±0.004 0.077 1.780±0.003 0.050 0.700±0.019 0.351 QRZ 1632 46.28 ±4.34 49.59 0.205±0.004 0.077 1.733±0.001 0.025 0.395±0.008 0.156 TCW 9199 14.07 ±2.03 51.22 0.133±0.001 0.069 1.715±0.001 0.054 0.643±0.008 0.362 THZ 2676 26.09 ±4.69 54.62 0.163±0.003 0.078 1.695±0.002 0.044 0.611±0.011 0.281 TKNZ 1798 86.31 ±3.14 44.95 0.226±0.003 0.065 1.735±0.002 0.034 0.543±0.008 0.165 TMWZ 809 22.97 ±19.73 65.95 0.110±0.005 0.064 1.787±0.002 0.030 0.567±0.025 0.334 TRWZ 629 42.11 ±5.33 45.03 0.189±0.006 0.074 1.789±0.006 0.069 0.786±0.033 0.400 TUWZ 4648 51.77 ±1.37 38.43 0.165±0.002 0.074 1.719±0.001 0.045 0.957±0.016 0.528 WEL 4144 -81.07±2.53 48.37 0.144±0.003 0.082 1.747±0.002 0.065 0.711±0.013 0.414 Total 83183 47.59±1.05 57.85 0.158±0.001 0.077 1.744±0.001 0.079 0.847±0.004 0.645 Station . Events . Mean ϕ(°) . SD ϕ . Mean δt(s) . SD δt . Mean Vp/Vs . SD Vp/Vs . Mean per cent Aniso . Per cent Aniso SD . BHW 3515 81.14 ±7.40 62.82 0.136±0.002 0.066 1.724±0.002 0.066 0.759±0.015 0.424 BSWZ 4781 52.82 ±1.76 43.38 0.181±0.002 0.076 1.743±0.002 0.084 1.453±0.023 0.801 CAW 6499 17.05 ±3.29 55.91 0.171±0.002 0.081 1.744±0.001 0.051 0.683±0.009 0.342 CMWZ 2162 -82.10±10.43 64.11 0.173±0.003 0.072 1.819±0.008 0.181 2.371±0.065 1.534 CRSZ 1588 67.47 ±4.07 48.34 0.168±0.004 0.085 1.808±0.014 0.281 2.118±0.060 1.167 DUWZ 2639 28.60 ±3.68 50.79 0.203±0.003 0.086 1.738±0.002 0.044 0.623±0.011 0.262 HOWZ 2015 69.93 ±6.58 57.45 0.139±0.003 0.075 1.732±0.002 0.045 0.673±0.019 0.419 KHZ 6343 -32.95±2.89 53.83 0.149±0.002 0.073 1.736±0.001 0.055 0.713±0.011 0.421 KIW 5497 53.85 ±3.11 53.87 0.152±0.002 0.073 1.718±0.001 0.044 0.629±0.008 0.296 MRNZ 2375 52.98 ±3.78 50.37 0.205±0.003 0.081 1.741±0.001 0.030 0.535±0.009 0.214 MSWZ 3481 60.46 ±2.02 43.02 0.142±0.003 0.077 1.784±0.002 0.058 0.560±0.011 0.327 MTW 7300 61.60 ±3.03 55.57 0.150±0.002 0.071 1.771±0.001 0.039 0.625±0.008 0.330 NNZ 1832 50.30 ±10.84 63.55 0.206±0.004 0.082 1.729±0.002 0.040 0.684±0.015 0.298 OGWZ 4232 37.96 ±1.16 34.13 0.123±0.002 0.067 1.736±0.001 0.047 0.478±0.009 0.282 PAWZ 1908 75.94 ±7.12 58.17 0.188±0.004 0.078 1.774±0.002 0.052 0.691±0.015 0.329 PLWZ 1460 -22.50±6.66 55.31 0.166±0.004 0.077 1.780±0.003 0.050 0.700±0.019 0.351 QRZ 1632 46.28 ±4.34 49.59 0.205±0.004 0.077 1.733±0.001 0.025 0.395±0.008 0.156 TCW 9199 14.07 ±2.03 51.22 0.133±0.001 0.069 1.715±0.001 0.054 0.643±0.008 0.362 THZ 2676 26.09 ±4.69 54.62 0.163±0.003 0.078 1.695±0.002 0.044 0.611±0.011 0.281 TKNZ 1798 86.31 ±3.14 44.95 0.226±0.003 0.065 1.735±0.002 0.034 0.543±0.008 0.165 TMWZ 809 22.97 ±19.73 65.95 0.110±0.005 0.064 1.787±0.002 0.030 0.567±0.025 0.334 TRWZ 629 42.11 ±5.33 45.03 0.189±0.006 0.074 1.789±0.006 0.069 0.786±0.033 0.400 TUWZ 4648 51.77 ±1.37 38.43 0.165±0.002 0.074 1.719±0.001 0.045 0.957±0.016 0.528 WEL 4144 -81.07±2.53 48.37 0.144±0.003 0.082 1.747±0.002 0.065 0.711±0.013 0.414 Total 83183 47.59±1.05 57.85 0.158±0.001 0.077 1.744±0.001 0.079 0.847±0.004 0.645 Notes: Only high-quality measurements (described in Section 3.2.2) are presented in the summary table. ± represents the 95 per cent confidence intervals. SD shows the standard deviation of the measurements. Station averages and standard deviation for fast orientations were performed using circular statistics. Event represents the number of event–station pair used for analysis. Total shows the summary of all stations for GeoNet stations. Open in new tab Table 2. Summary statistics of Fast orientation (ϕ), delay time (δt), Vp/Vs ratio and percent anisotropy (per cent Aniso) with their 95 per cent confidence interval for DPRI temporary stations. Station . Events . Mean ϕ(°) . SD ϕ . Mean δt(s) . SD δt . Mean Vp/Vs . SD Vp/Vs . Mean per cent Aniso . Per cent Aniso SD . AGR 277 44.78 ±11.59 51.11 0.167±0.012 0.096 1.648±0.027 0.230 1.222±0.105 0.851 CRF 89 -54.36 ±16.48 47.61 0.033±0.008 0.037 1.650±0.053 0.243 0.456±0.120 0.556 CVR 1593 3.75 ±18.71 69.41 0.183±0.004 0.073 1.789±0.007 0.141 1.508±0.036 0.725 GBR 139 -51.61 ±15.37 50.12 0.137±0.013 0.074 1.739±0.013 0.073 1.572±0.141 0.827 IKR 33 -47.57 ±22.85 44.73 0.146±0.026 0.072 1.767±0.038 0.107 1.027±0.188 0.534 JSP 3726 65.59 ±1.67 40.11 0.155±0.003 0.088 1.641±0.007 0.211 1.407±0.030 0.898 KVR 638 84.34 ±13.42 59.38 0.131±0.006 0.070 1.769±0.012 0.153 1.631±0.071 0.874 MLF 4785 -2.43 ±1.71 42.91 0.117±0.002 0.069 1.697±0.009 0.308 1.053±0.018 0.616 MTV 7281 61.67 ±3.04 55.60 0.150±0.002 0.071 1.771±0.001 0.039 0.623±0.008 0.328 SJQ 2396 57.95 ±2.12 40.41 0.151±0.003 0.066 1.670±0.004 0.091 1.010±0.022 0.536 SVR 1789 49.70 ±5.45 53.85 0.164±0.004 0.084 1.822±0.006 0.126 1.270±0.034 0.712 WJM 1754 8.86 ±2.62 41.47 0.231±0.004 0.077 1.702±0.003 0.066 1.184±0.025 0.524 Total 19076 34.24 ±2.20 57.86 0.157±0.001 0.084 1.686±0.003 0.209 1.199±0.010 0.698 All 102260 45.15 ±0.97 58.09 0.157±0.001 0.079 1.741±0.001 0.088 0.922±0.004 0.688 Station . Events . Mean ϕ(°) . SD ϕ . Mean δt(s) . SD δt . Mean Vp/Vs . SD Vp/Vs . Mean per cent Aniso . Per cent Aniso SD . AGR 277 44.78 ±11.59 51.11 0.167±0.012 0.096 1.648±0.027 0.230 1.222±0.105 0.851 CRF 89 -54.36 ±16.48 47.61 0.033±0.008 0.037 1.650±0.053 0.243 0.456±0.120 0.556 CVR 1593 3.75 ±18.71 69.41 0.183±0.004 0.073 1.789±0.007 0.141 1.508±0.036 0.725 GBR 139 -51.61 ±15.37 50.12 0.137±0.013 0.074 1.739±0.013 0.073 1.572±0.141 0.827 IKR 33 -47.57 ±22.85 44.73 0.146±0.026 0.072 1.767±0.038 0.107 1.027±0.188 0.534 JSP 3726 65.59 ±1.67 40.11 0.155±0.003 0.088 1.641±0.007 0.211 1.407±0.030 0.898 KVR 638 84.34 ±13.42 59.38 0.131±0.006 0.070 1.769±0.012 0.153 1.631±0.071 0.874 MLF 4785 -2.43 ±1.71 42.91 0.117±0.002 0.069 1.697±0.009 0.308 1.053±0.018 0.616 MTV 7281 61.67 ±3.04 55.60 0.150±0.002 0.071 1.771±0.001 0.039 0.623±0.008 0.328 SJQ 2396 57.95 ±2.12 40.41 0.151±0.003 0.066 1.670±0.004 0.091 1.010±0.022 0.536 SVR 1789 49.70 ±5.45 53.85 0.164±0.004 0.084 1.822±0.006 0.126 1.270±0.034 0.712 WJM 1754 8.86 ±2.62 41.47 0.231±0.004 0.077 1.702±0.003 0.066 1.184±0.025 0.524 Total 19076 34.24 ±2.20 57.86 0.157±0.001 0.084 1.686±0.003 0.209 1.199±0.010 0.698 All 102260 45.15 ±0.97 58.09 0.157±0.001 0.079 1.741±0.001 0.088 0.922±0.004 0.688 Notes: Only high-quality measurements (described in Section 3.2.2) are presented in the summary table. ± represents the 95 per cent confidence intervals. SD shows the standard deviation of the measurements. Station averages and standard deviation for fast orientations were performed using circular statistics. Event represents the number of event-station pair used for analysis. Total represents the summary of all stations for DPRI stations. All shows the summary of all stations for both GeoNet and DPRI stations. Open in new tab Table 2. Summary statistics of Fast orientation (ϕ), delay time (δt), Vp/Vs ratio and percent anisotropy (per cent Aniso) with their 95 per cent confidence interval for DPRI temporary stations. Station . Events . Mean ϕ(°) . SD ϕ . Mean δt(s) . SD δt . Mean Vp/Vs . SD Vp/Vs . Mean per cent Aniso . Per cent Aniso SD . AGR 277 44.78 ±11.59 51.11 0.167±0.012 0.096 1.648±0.027 0.230 1.222±0.105 0.851 CRF 89 -54.36 ±16.48 47.61 0.033±0.008 0.037 1.650±0.053 0.243 0.456±0.120 0.556 CVR 1593 3.75 ±18.71 69.41 0.183±0.004 0.073 1.789±0.007 0.141 1.508±0.036 0.725 GBR 139 -51.61 ±15.37 50.12 0.137±0.013 0.074 1.739±0.013 0.073 1.572±0.141 0.827 IKR 33 -47.57 ±22.85 44.73 0.146±0.026 0.072 1.767±0.038 0.107 1.027±0.188 0.534 JSP 3726 65.59 ±1.67 40.11 0.155±0.003 0.088 1.641±0.007 0.211 1.407±0.030 0.898 KVR 638 84.34 ±13.42 59.38 0.131±0.006 0.070 1.769±0.012 0.153 1.631±0.071 0.874 MLF 4785 -2.43 ±1.71 42.91 0.117±0.002 0.069 1.697±0.009 0.308 1.053±0.018 0.616 MTV 7281 61.67 ±3.04 55.60 0.150±0.002 0.071 1.771±0.001 0.039 0.623±0.008 0.328 SJQ 2396 57.95 ±2.12 40.41 0.151±0.003 0.066 1.670±0.004 0.091 1.010±0.022 0.536 SVR 1789 49.70 ±5.45 53.85 0.164±0.004 0.084 1.822±0.006 0.126 1.270±0.034 0.712 WJM 1754 8.86 ±2.62 41.47 0.231±0.004 0.077 1.702±0.003 0.066 1.184±0.025 0.524 Total 19076 34.24 ±2.20 57.86 0.157±0.001 0.084 1.686±0.003 0.209 1.199±0.010 0.698 All 102260 45.15 ±0.97 58.09 0.157±0.001 0.079 1.741±0.001 0.088 0.922±0.004 0.688 Station . Events . Mean ϕ(°) . SD ϕ . Mean δt(s) . SD δt . Mean Vp/Vs . SD Vp/Vs . Mean per cent Aniso . Per cent Aniso SD . AGR 277 44.78 ±11.59 51.11 0.167±0.012 0.096 1.648±0.027 0.230 1.222±0.105 0.851 CRF 89 -54.36 ±16.48 47.61 0.033±0.008 0.037 1.650±0.053 0.243 0.456±0.120 0.556 CVR 1593 3.75 ±18.71 69.41 0.183±0.004 0.073 1.789±0.007 0.141 1.508±0.036 0.725 GBR 139 -51.61 ±15.37 50.12 0.137±0.013 0.074 1.739±0.013 0.073 1.572±0.141 0.827 IKR 33 -47.57 ±22.85 44.73 0.146±0.026 0.072 1.767±0.038 0.107 1.027±0.188 0.534 JSP 3726 65.59 ±1.67 40.11 0.155±0.003 0.088 1.641±0.007 0.211 1.407±0.030 0.898 KVR 638 84.34 ±13.42 59.38 0.131±0.006 0.070 1.769±0.012 0.153 1.631±0.071 0.874 MLF 4785 -2.43 ±1.71 42.91 0.117±0.002 0.069 1.697±0.009 0.308 1.053±0.018 0.616 MTV 7281 61.67 ±3.04 55.60 0.150±0.002 0.071 1.771±0.001 0.039 0.623±0.008 0.328 SJQ 2396 57.95 ±2.12 40.41 0.151±0.003 0.066 1.670±0.004 0.091 1.010±0.022 0.536 SVR 1789 49.70 ±5.45 53.85 0.164±0.004 0.084 1.822±0.006 0.126 1.270±0.034 0.712 WJM 1754 8.86 ±2.62 41.47 0.231±0.004 0.077 1.702±0.003 0.066 1.184±0.025 0.524 Total 19076 34.24 ±2.20 57.86 0.157±0.001 0.084 1.686±0.003 0.209 1.199±0.010 0.698 All 102260 45.15 ±0.97 58.09 0.157±0.001 0.079 1.741±0.001 0.088 0.922±0.004 0.688 Notes: Only high-quality measurements (described in Section 3.2.2) are presented in the summary table. ± represents the 95 per cent confidence intervals. SD shows the standard deviation of the measurements. Station averages and standard deviation for fast orientations were performed using circular statistics. Event represents the number of event-station pair used for analysis. Total represents the summary of all stations for DPRI stations. All shows the summary of all stations for both GeoNet and DPRI stations. Open in new tab Rose diagrams of ϕ orientations recorded at each station, overlain on a basemap of the digital topography, are shown on Fig. 4. SHmax orientations from Townend et al. (2012) and Balfour et al. (2005) are also shown as green bow ties with wedges showing the 95 per cent confidence interval. ϕ orientations for most of the stations show a dominant NE-SW orientation, which is parallel to the strike of the geological structures in the region. A few stations have NW-SE orientations. Figure 4. Open in new tabDownload slide Rose diagrams (circular histograms) showing the ϕ orientation results from local S-phase events, ϕ. Rose diagrams are plotted on the stations at which measurements were made with a basemap of the digital topography and bathymetry. The green bow ties shows the crustal SHmax orientations from Townend et al. (2012) and Balfour et al. (2005) with wedges showing the 95 per cent confidence interval. The black lines are the active faults. The two red rose diagrams are measurements from Audoine et al. (2000) showing bimodal distributions. Around the MFS, ϕ orientations are generally parallel to the strike of the active faults. The consistent trend suggests a relation between ϕ and the structural trend around the central and western MFS. Yet, some stations, KHZ, IKR, CRF and GBR, at the eastern coastal edge of the MFS yield NW-SE orientations of ϕ, which is perpendicular to the structural trend and subparallel to the crustal SHmax orientations. At station CMWZ and KVR, we observe a bi-modal distribution of ϕ orientations, with one subset subparallel to the SHmax orientation and the other subset subparallel to the strike of the faults. At stations BSWZ and TUWZ, ϕ orientations are subparallel to the orientation of the nearby Awatere and Wairau Faults respectively and perpendicular to the nearby SHmax orientation. Around the southern end of the North island, the ϕ orientation at most stations exhibit significant variations from each other. Some stations exhibit a bi-modal distribution while others show a dominant NE-SW direction, with a few stations exhibiting trends that are neither parallel to SHmax nor to the strike the major faults. At stations TMWZ and MTW, a clear bi-modal distribution of fast polarization orientation is observed. One subset trends NE-SW and the other subset trends NW-SE (Fig. 4). To investigate temporal variations resulting from the 2016 Kaikōura main shock, ϕ measurements were split into two subsets: before (from 2013 January 01 to 2016 November 16) and after (from 2016 November 14 to 2018 June 30) the Kaikōura main shock. The 2013 Cook Strait earthquake sequence produced ample seismicity to determine the background ϕ orientation prior to the 2016 Kaikōura main shock. Fig. 5 shows rose diagrams of ϕ measurements before (Fig. 5a) and after the 2016 Kaikōura main shock (Fig. 5b). Figure 5. Open in new tabDownload slide SWS fast orientation results from local S-phase events for before, A (from 01/01/13 to 13/11/16) and after, B (from 14/11/16 to 30/06/18) (bottom) Kaikōura main shock. The ϕ orientation determined from events recorded at each station are represented as rose diagrams and plotted on a basemap of the digital topography. Rose diagrams are plotted on the stations at which measurements were made. The green bow ties shows the SHmax orientations from Townend et al. (2012) and Balfour et al. (2005) with wedges showing the 95 per cent confidence interval. The blue lines are the faults that ruptured during the Kaikōura Earthquake (Langridge et al. 2016). 4.1 Spatial averaging Using the technique described in Section 3.2.3, we performed spatial averaging analysis for more than 102 260 source–receiver ϕ measurements. Based on the criteria outlined in Section 3.2.3, 4579 out of 4720 grid elements created were used for analysis (see Fig. S7, Supporting Information). Circular averages of ϕ measurements in each grid square are represented by bar (coloured by mean ϕ value) and plotted at the centre of each grid square (Fig. 6). The spatial averaging results shows a significant spatial variation. We observe a dominant NE-SW (see yellow bars on Fig. 6) trend, which is subparallel to the strike of the major active faults and also the topographic trends. We also note a clear contrast of perpendicular fast orientations (blue bars on Fig. 6) around the eastern end of the MFS (the region that ruptured during the Kaikōura earthquake). This NW-SE feature existed before the 2016 Kaikōura sequence (see Fig. S8, Supporting Information). Figure 6. Open in new tabDownload slide Spatial averages of fast polarization orientation with weighting inversely proportional to the square of the distance from the station. The bars are coloured by the orientation of the mean ϕ. Black lines are known active faults from the NZ Active Faults Database and the green lines are the faults that ruptured during the 2016 Kaikōura main shock (Litchfield et al. 2017). To investigate the spatiotemporal variation of ϕ measurements associated with the 2016 Kaikōura main shock we performed spatial averaging analysis for events before and after the earthquake. Using the method described in Section 3.2.4, we made quantitative comparisons of our spatially averaged ϕ measurements between events before and after the 2016 Kaikōura main shock to search for potential spatiotemporal variations associated with the 2016 Kaikōura main shock. The orientations of spatially averaged results for before (green bars) and after (blue bars) are shown in Fig. 7. Generally, the orientations of the two angles are in agreement, which indicates no significant spatiotemporal variation (in most cases the blue after bar is invisible under the exactly matching green before bar). The mean and median value of F (Section 3.2.4) is 0.92 ± 0.011 and 0.99 respectively (thus a mean and median difference of 23° and 10° respectively) indicates a strong agreement between the orientations before and after the main shock. Out of the 898 co-located measurements, 88 per cent of them had a test statistic value, F, above 0.8 (37°) with only 1 per cent less than 0.2 (78°) (Fig. 7b). The general agreement observed, also seen from the scatter plot (Fig. 7c) of before and after measurements, suggests that there is no significant temporal variation associated with the 2016 Kaikōura earthquake. Although we observe a few patches of disagreement around the propagation path of the 2016 Kaikōura main shock, we attribute this to spatial variation rather than temporal variation since ray paths for after measurements are slightly different from before measurements (see Fig. S8, Supporting Information). The red patch at the north-eastern edge of the study region may be due to an edge effect from the grid. Due to limited station coverage off-shore, we limit our discussion to onshore results. Figure 7. Open in new tabDownload slide Quantitative comparison between splitting measurements for before and after the 2016 Kaikōura earthquake. The background colour is a contour of the absolute value of the cosine of the difference between the orientations, coloured by the scale at the bottom. Note that 1 represents no difference between the orientations. Before and after Kaikōura earthquake spatial averaging estimates are green and blue bars, respectively. The descriptive statistics of the difference between orientations at the closest gridpoints is displayed on the histogram in (b). A scatter plot of the angles and histogram of each axis are shown in (c). 5 DISCUSSION 5.1 Overall station trends The observed delay times, varying between near 0 and 0.372 s, are in agreement with other SWS studies of crustal earthquakes (Audoine et al. 2000; Balfour et al. 2005; Crampin et al. 2015; Cochran & Kroll 2015; Savage et al. 2016). Vp/Vs ratios varying from 1.58 to 2.22 with an average of 1.746 ± 0.001, are also consistent with tomography results from Eberhart-Phillips & Fry (2018) and close to the global average of 1.76 for continental crust (Christensen 1996). There is however substantial variability in ϕ measurements (the standard deviation of the ϕ measurements is greater than 30° at most stations). The observed large scatter of ϕ measurements suggests that the source of anisotropy may be a combination of more than one mechanism or that the mechanism varies rapidly spatially (Peng & Ben-Zion 2004). The significant spatial variation we observe illustrates the complex regional tectonics and heterogeneous structures around central New Zealand. The average δt we observed (Tables 1 and 2) is about 6–14 per cent of the values reported in SKS studies (∼1.6 s) around New Zealand (Klosko et al. 1999; Zietlow et al. 2014). This small value of δt suggests that there is a deeper source of anisotropy beneath these stations that is measured by SKS and our measurements see only shallow anisotropy. A similar conclusion has been reached by previous studies (Audoine et al. 2000; Karalliyadda & Savage 2013) where they reported delay times for SKS measurement far greater than the observed local S-phase measurements. Similarly most studies of SKS waves assume the mantle is the main source of anisotropy (e.g. Savage 1999). The bi-modal distribution of fast orientation measurements observed at stations TMWZ and MTW on the North Island (Fig. 4) reinforces previous observations. Audoine et al. (2000) reported a similar bi-modal distribution of ϕ measurements at stations LKER and LKIR (red rose diagram on Fig. 4), which were in close proximity to TMWZ and MTW. At station HOWZ, the dominant ϕ orientation is parallel to faulting in the area and also subparallel to SHmax. The correlation between ϕ orientation and the active fault orientations at these stations suggests that the sources of anisotropy beneath them could be structurally controlled. However, the faults are also parallel to the strike of the topography, which is consistent with near surface gravitational stress. Fast orientations at stations WEL and TCW (on either side of the Cook Strait; Fig. 4) are consistent with recent studies (Balfour et al. 2005; Evanzia et al. 2017) around the Marlborough and Wellington region. However, they show neither an agreement with SHmax orientations, nor with tectonic structures in the region. Around the Marlborough region, the faulting is mainly strike slip, thus based on Anderson’s theory of faulting (Anderson 1905; Healy et al. 2012), we expect the faults to generally strike at an angle of ∼30° to the orientation of SHmax. Stations KHZ, IKR, CRF and GBR (and a subset of ϕ measurements at stations KVR and CMWZ), located around the eastern end of the MFS, have ϕ orientations that are subparallel to the orientation of SHmax rather than the fault fabric, suggesting stress-aligned microcracks (stress controlled anisotropy) exist (Figs 4 and 5). However, at these stations, the role of structural anisotropy cannot be completely ruled out since there are subsets of the ϕ measurements that are subparallel to the major structural features. An abrupt change from NE-SW to an NW-SE trend of the averaged fast orientation (blue bars on Fig. 6) is observed at the edge of the high irregular topography boundary around the eastern end of the MFS. This suggests that gravitational stresses could be contributing to the source of anisotropy around this region. In Section 5.5, we quantitively compare averaged ϕ to the orientation of |$S_{H_{\mathrm{ max}}}^{\mathrm{ Grav}}$|⁠. 5.2 Search for temporal variation The use of temporal variation in SWS measurements to monitor precursory changes preceeding major earthquakes is controversial (e.g. Aster et al. 1990; Crampin 1990; Crampin et al. 2015). There have been studies that have used temporal variations to monitor coseismic stress changes (Saiga et al. 2003) and post-seismic fault healing after major earthquakes (Tadokoro et al. 1999; Tadokoro & Ando 2002). Around the central Apennines of Italy, Piccinini et al. (2006) reported temporal variations in anisotropic parameters during the days before and after the occurrence of the 1997 Umbria-Marche seismic sequence. They suggested that the temporal variation is caused by the changes in local stress condition and fluid pressure. However, they concluded that spatial sampling of the selected ray paths may vary with time and the possible contribution of spatial variations of anisotropic parameters can not be excluded. Lucente et al. (2010) also observed a clear temporal variation of seismic anisotropy and velocities before the 2009 Mw 6.3 L’Aquila earthquake, Italy. They inferred that fluids, related to dilatancy-diffusion processes, played a key role in the fault failure process. In southern California, Yang et al. (2011) searched for possible temporal variations of crustal anisotropy associated with big earthquakes, but they did not observe any clear temporal variations in delay times or fast orientation. In addition, studies by Liu (2004) of crustal anisotropy around the San Andreas Fault near Parkfield did not show systematic temporal variations. Moreover, a systematic analysis of crustal anisotropy associated with the 1999 Duzce main shock (Mw 7.1) along the North Anatolian Fault by Peng & Ben-Zion (2004) also showed no temporal variations before, during or after the main shock. One major issue with interpreting variations in observed SWS measurements is distinguishing temporal changes in anisotropy from spatial variations in anisotropy due to differing source locations (e.g. Peng & Ben-Zion 2005). Since the propagation direction through the crust may be different for differing source locations, the wave is likely to sample different rock masses, resulting in different splitting parameters. In a complex setting, like the Marlborough and Wellington regions where the observed anisotropy may be controlled by different mechanisms, the determination of temporal variations is challenging, even if they exist. This is due to contamination associated with the variations of ray paths as a result of changes in source locations. However, temporal variation in stress has been inferred for large earthquakes (Lucente et al. 2010; Nakata & Snieder 2012; Holt et al.2013). 5.2.1 Delay time and Vp/Vs ratios We systematically searched for temporal variations in delay time and Vp/Vs associated with the 2013 Cook Strait and 2016 Kaikōura earthquakes. We applied the moving average technique described in Section 3.2.5 because our δt and Vp/Vs measurements exhibit large scatter without any clear trends. We tried a 10-, 20- and 30-d moving window, and all showed similar results. A time-series plot of δt and Vp/Vs measurements revealed apparent temporal variation at two stations (CAW and KIW, see Figs S9a and b, Supporting Information). We stacked measurements from all 36 station to enhance this variation observed in the individual stations (Fig. 8a). To test whether the apparent variations originate from changes in the medium properties over time rather than changes in the source location or the propagation path, we selected a cluster of measurements within 20 km radius of station CMWZ (which is located in a region with a good number of measurements before and after both the 2013 and 2016 earthquakes sequence (Fig. 8d). The apparent temporal variation of δt and Vp/Vs observed at individual stations with measurements from different locations was not seen when temporal analysis was done using only the small cluster of measurements (Fig. S9a1 and b1). Similarly, the clear temporal variation observed with the stacked measurements from different locations disappears when we stack individual station measurements from earthquakes with similar locations (Fig. 8b). Moreover, analysis of measurements with similar backazimuths revealed no clear temporal variation. This indicates that spatial variations of parameters strongly affect the observed apparent temporal variations. Figure 8. Open in new tabDownload slide Time-series of delay times, δt and Vp/Vs ratio for values from 2013 January to 2018 June are shown in panels (a) and (b). (a) and (b) are plots for all measurements and the small cluster of measurements, respectively, for all stations stacked. (c) and (d) shows the event locations for all measurements and the small cluster of measurements respectively. Blue, red and orange vertical bars represent the main shock of the 2013 Mw 6.5 2013 Seddon, the 2013 Mw 6.6 Lake Grassmere, the 2014 Mw 6.2 Eketahuna and the 2016 Mw 7.8 Kaikōura earthquake. Green horizontal line shows the median value. Green dashed line represents the estimated error using bootstraps. Red and blue lines are the 20-d moving median of Vp/Vs ratio and δt, respectively. 5.2.2 Fast orientations To reveal the state of anisotropy over the |$5\frac{1}{2}$| yr period, we performed the spatial averaging analysis described in Section 3.2.3 with the same parameters for measurements within each year. As shown in Fig. 9, there are significant spatial variations in the fast orientations but no systematic temporal variations. The generally consistent spatial variations over the years demonstrate that neither the 2013 Cook Strait nor the 2016 Kaikōura earthquakes changed the stress directions significantly. Instead, we observe a consistent NE-SW trend (see yellow bars on Fig. 9) that strikes subparallel to the active faulting around the study area. One remarkable feature is the consistent NW-SE orientation of fast orientation (blue bars) around Kaikōura region. This indicates that the NW-SE feature did not originate from the 2016 Kaikōura earthquake sequence. Figure 9. Open in new tabDownload slide State of anisotropy over the |$5\frac{1}{2}$| yr period, The bars show the average fast orientation obtained at the each regular grid. We apply the same criteria as in Section 4.1. The green line are the fault the ruptured during the Kaikōura main shock (Litchfield et al. 2017). The absence of temporal variations in ϕ suggests that any change in the stress field around central New Zealand produced anisotropy changes that are small compared to structural anisotropy. Iidaka & Obara (2013) and Nakata & Snieder (2012) reported a similar observation with the M9.0 Tohoku-Oki megathrust earthquake. Iidaka & Obara (2013)’s analysis of fast polarization orientations before and after the event revealed no changes and they suggested that the orientation of the maximum stress axis in the region of increased seismicity did not change after the 2011 Tohoku earthquake (Iidaka & Obara 2013). In the near-surface (i.e. upper few hundred metres) Nakata & Snieder (2012), using seismic interferometry, reported no significant variations in fast polarization orientation after the M9.0 Tohoku-Oki earthquake. Focal mechanism inversion studies around central New Zealand following the 2016 Kaikōura earthquake (Okada et al. 2019; Okada, personal communication) also revealed no significant variation in the stress field. This corroborates our findings that the 2016 Kaikōura earthquake did not significantly alter the stress field around central New Zealand. 5.3 Depth extent of anisotropy Several studies have attempted to constrain the depth extent of anisotropy. Delay times observed by Audoine et al. (2000) showed no correlation with depth, which they attributed to the higher frequency phases used in their studies, thus sampling mainly the upper crustal anisotropy. In the Marlborough region, Balfour et al. (2005) also did not observe any systematic association between δt and depth and also suggested that the source of anisotropy may reside in the upper crust. Many similar studies conclude that anisotropy is confined to the upper few kilometres of the crust because their analysis often shows no clear correlation between delay time and hypocentral depth (Shih & Meyer 1990; Gledhill 1991; Audoine et al. 2000; Balfour et al. 2005, 2012). We analysed the association between delay time and hypocentral depth and epicentral distance at individual stations (Fig. 10 and Fig. S10, Supporting Information, respectively). The spiky appearance of the depth histogram results from the GeoNet catalogue (they often automatically assign the depth of locations with large uncertainty to the closest depth in the velocity model). That of the delay time histogram is due to rounding δt measurements. The correlation coefficients at the analysed stations were less than 0.2, indicating a lack of clear correlation between delay time and hypocentral depth or epicentral distance. This suggests that we are mainly measuring anisotropy that resides only in the uppermost few kilometres of the crust (Shih & Meyer 1990; Audoine et al. 2000; Balfour et al. 2005). Our results are also consistent with previous studies in this region and other tectonic regions (Gledhill 1991; Audoine et al. 2000; Balfour et al. 2005; Zhang et al. 2007). The shallow sources of anisotropy observed may also suggest that topographic stresses (due to the difference in density between rock and air in the high topography areas) can contribute to the observed anisotropic pattern. Figure 10. Open in new tabDownload slide Delay times versus depth for station KHZ. Distribution of the delay time and depths are plotted on the right and bottom, respectively. The scattered points are coloured by the backazimuth and plotted on a basemap of the kernel density (effectively heat map). High-density regions are shaded white with low regions black. The red line is the least square regression line with the equation and p-value displayed on the lower right corner of the plot. 5.4 Discriminating between stress-induced and structurally controlled seismic anisotropy Discriminating between stress-induced and structurally controlled seismic anisotropy in most tectonic settings is often challenging because SHmax is often subparallel to major structural features (Savage et al. 1989; Peng & Ben-Zion 2004; Johnson et al. 2011). But in places where the SHmax orientation is at a high angle to the orientation of the major structural features (such as around the eastern end of the MFS near the Kaikōura region), it is possible to distinguish between stress-induced and structurally controlled anisotropy (Zinke 2000; Boness & Zoback 2006). Previous crustal anisotropy studies in central New Zealand have suggested that anisotropy in the region is controlled by a combination of both structures and stresses (e.g. Audoine et al. 2000; Balfour et al. 2005; Karalliyadda & Savage 2013; Evanzia et al. 2017). Here, we attempt to distinguish between stress-induced seismic anisotropy and structurally controlled seismic anisotropy. We quantitatively compared our spatially averaged ϕ measurements to the average orientation of the closest active fault segments following the method described in Section 3.2.4. Fig. 11 indicates that our spatially averaged ϕ measurements, ϕavg, are overall oriented similar to the NE-SW strike of the active fault orientation. Out of 2297 measurements we compared, 70 per cent of them had a test statistic value, F, above 0.8 with only 3 per cent less than 0.2. The mean and median value of F was 0.81 ± 0.009 and 0.92, respectively (thus a mean and median difference of 35° and 23°, respectively) indicating a strong agreement. The general agreement of average ϕ measurements with the orientation of active faults in polygon 1 (we exclude measurements in polygon 3 whose borders are contained within polygon 1) suggests that the observed anisotropy in polygon 1 may be controlled by structures. Around the eastern end of the MFS and part of the southern end of the North Island (polygons 2 and 3 on Fig. 11) the azimuth of ϕavg departs from the NE-SW trend to an almost E-W trend. This departure shows a strong disagreement with the orientation of the active faults (red colour on Fig. 11). This suggests that the anisotropy pattern observed within polygons 2 and 3 is not structurally controlled. Agreement of the average ϕ with the orientation of active faults has also been reported elsewhere in areas of active transpression (Savage 1999). However, the structural trends are also parallel to topographic trends. In the following sections we examine how both tectonic stresses (from focal mechanism inversion) and gravitational stresses may explain part of the observed anisotropic pattern. Figure 11. Open in new tabDownload slide Quantitative comparison between spatially averaged ϕ measurements and average orientation of the closest active fault. The background colour is a contour of the absolute value of the cosine of the difference between the orientations, coloured by the scale at the bottom. Note that 1 represents no difference between the orientations. The descriptive statistics of the difference between orientations at the closest gridpoints is displayed on the histogram in B. SHmax are represented by dark green bow ties. A scatter plot of the angles and histogram of each axis are shown in C. 5.5 Comparing spatially averaged ϕ with stress indicators The gravitational forces exerted by dense material placed next to less dense material, such as high mountains next to air, can significantly contribute to the total stress field, which can result in variations in crustal differential stresses (Flesch et al. 2001). This induced stress is referred to as gravitationally induced stress, |$S_{H_{\mathrm{ max}}}^{\mathrm{ Grav}}$| (Flesch et al. 2001; Hirschberg et al. 2018). The orientation and magnitude of gravitationally induced stresses can be modelled with a known crustal rock density, crustal thickness and surface topography. We used the Flesch et al. (2001) method implemented by Hirschberg et al. (2018) using finite-difference methods to calculate the orientation and magnitude of gravitationally induced stresses. The implementation of the Flesch et al. (2001) method involves solving force balance equations, assuming an isotropic viscous medium. Following Hirschberg et al. (2018), we assume that the viscous medium is subjected to stresses and also assumed zero deviatoric stresses for the horizontal boundary condition. For our calculations, we used the 250 m gridded bathymetric data set (Mitchell et al. 2012) around central New Zealand and assumed a crustal thickness and density of 30 km and 2850 kg m−3, respectively. We modelled the orientation and magnitude of gravitationally induced stresses at a 30 km gridpoint spacing. At every gridpoint, we compared averaged ϕ to the orientation of the gravitationally induced stresses (Fig. 12) following the method in Section 3.2.4. Figure 12. Open in new tabDownload slide Quantitative comparison between spatially averaged ϕ measurements and orientation of gravitationally induced stress. The background colour is a contour of the absolute value of the cosine of the difference between the orientations, coloured according to the scale at the bottom. Note that 1 represents no difference between the orientations. The descriptive statistics of the difference between orientations at the closest gridpoints is displayed on the histogram in (b). Part (c) is a 3-D perspective image of the topography and bathymetry map around Kaikōura region. Overall, we observe a general agreement between ϕavg and |$S_{H_{\mathrm{ max}}}^{\mathrm{ Grav}}$|⁠, with the mean and median value of F of 0.77 ± 0.062 and 0.828, respectively (a mean and median difference of 39° and 34°, respectively) (Fig. 12b). Around the Marlborough region (Fig. 12), the estimated |$S_{H_{\mathrm{ max}}}^{\mathrm{ Grav}}$| is parallel to the ϕavg orientation. The NE-SW trend of |$S_{H_{\mathrm{ max}}}^{\mathrm{ Grav}}$| is also in agreement with the trend of the structures. Combined with the results in Section 5.4, this suggests the anisotropy pattern observed in polygon 1 could be controlled by either gravitational stresses or structure or a combination of the two. In polygon 2, where agreement between the ϕavg and the orientation of the active fault is poor (Fig. 11), we observe a general agreement between ϕavg and |$S_{H_{\mathrm{ max}}}^{\mathrm{ Grav}}$| (Fig. 12). We infer that the observed seismic anisotropy in polygon 2 (Figs 11 and 12) is best explained by the gravitationally induced stress that originates from the variations in topography around the region (Fig. 12c). We also examine how tectonic stresses explain the observed anisotropic patterns by comparing ϕ orientation measurements to the orientation of SHmax from focal mechanisms. Fig. S11, in the Supporting Information, illustrates the comparison of ϕ orientations with the nearest SHmax following the quantitative analysis described in Section 3.2.4. Around the northern South Island, we observe a general disagreement between ϕavg and SHmax orientations (mean F of 0.47 ± 0.01, 62° angle difference; see Fig. S11, Supporting Information). The general disagreement observed in the South Island is an indication that a different source of anisotropy, other than tectonic stresses is needed to explain the anisotropic pattern. Since the mean depth distribution of the focal mechanism inversion used for comparison is about 16.12 ± 2.08 km (see Fig. S11c, Supporting Information), we suggest that the tectonic stresses are acting on deeper structures than the measured anisotropic structure. Our results reinforce Balfour et al.’s observations that the tectonic stress field does not control the crustal anisotropy pattern around the Marlborough region. Conversely, in the North Island, the ϕavg orientations generally agree with SHmax orientations (mean F of 0.75 ± 0.02, 42° angle difference). This agreement is an indication that tectonic stresses may contribute to the observed anisotropic pattern, although as shown in Section 5.4 structure may also explain the anisotropy here. Evanzia et al. (2017) observed a similar agreement and suggested that stresses in the overriding plate and gravitational stresses control the anisotropic pattern in the Wellington region. 5.6 Comparing the geodetically determined plane of maximum shear stress with spatially averaged ϕ In the Marlborough and Wellington regions, the tectonics are complicated by the locked part of the Hikurangi subduction zone and a series of parallel strike-slip faults (MFS). In Fig. 13, spatial averaging of ϕ orientations is compared with the plane of maximum shear stress [we assume to be 45° to the principal axes (compressional and extensional plane) of the strain rate field (Ciucu 2010)]. We used the strain rate field result calculated by Lamb et al. (2018). The observed NE-SW trend of the ϕ orientation, which generally shows good agreement with the plane of maximum shear stress (polygon 1 on Fig. 13), probably reflects the effect of the shear deformation along the Australian/Pacific boundary. Figure 13. Open in new tabDownload slide Quantitative comparison between spatially averaged ϕ measurements and the plane of maximum shear stress (from Lamb et al. 2018). The background colour is a contour of the absolute value of the cosine of the difference between the orientations, coloured by the scale at the bottom. Note that 1 represents no difference between the orientations. The descriptive statistics of the difference between orientations at the closest gridpoints is displayed on the histogram in (b). Scatter plot of the angles in shown in (c). SHmax are represented by dark green bow ties. The red and green bars illustrate the plane of maximum shear stress and mean ϕ orientation at each gridpoint, respectively. The blues arrows indicate the principle axes strain rate, with compression pointing inwards and extension outwards. Present deformation along the plate-boundary zone in north eastern South Island, New Zealand is expected to influence anisotropy (Wilson et al. 2004; Balfour et al. 2005; Chen et al. 2013; Karalliyadda & Savage 2013). Okaya et al. (2016) suggested that the subhorizontal foliations with lineations (e.g. schist fabric), developed due to pervasive shearing in strike-slip fault zones, are a possible source of anisotropy. We therefore infer that the margin-parallel orientation of ϕ observed could be caused either by gravitational stress or shear deformation across the Marlborough region as a result of shearing along the active faulting (MFS) or both. The general disagreement of the orientations observed in polygon 2 underpins our earlier interpretation (Section 5.4) that these regions are likely controlled by gravitational stresses. 5.7 What controls seismic anisotropy in central New Zealand? Several studies around the MFS and the Wellington region (e.g. Gledhill 1991; Audoine et al. 2000; Matcham et al. 2000; Balfour et al. 2005; Karalliyadda & Savage 2013; Evanzia et al. 2017) suggest that the anisotropy in the crust is either structurally controlled or stress induced. Here, we try to delineate regions where the proposed source of anisotropy best explains our results by comparing the mean values on the F statistic in each polygon, both North and South Island, and the entire study region. We looked at how the structurally controlled hypothesis (using active fault and maximum shear plane orientations as a proxy) and stress induced hypothesis (relating to gravitational and tectonic stresses) explain the observed anisotropic pattern. In polygon 1 (Fig. 14, and for regions see Fig. 11), both gravitational stresses and structures can explain (with the mean F >0.8) the anisotropic pattern observed. This is because the active faults and foliation, or aligned minerals (due to shearing) in the area are oriented parallel to the strike of the topographic structures, which are also the ridges. Gravitational stress is the hypothesis that best explains the anisotropic pattern in polygon 2 (with the mean F >0.8). This is likely due to the sharp topography contrast (Fig. 12c) observed around the Kaikōura region. In polygon 3, no single hypothesis can explain the observation (with the mean F < 0.7). In the North Island, overall either stresses (tectonic and gravitational) or structures or both explain our results. However, in the South Island gravitational stresses and structures best explain our results. Overall the gravitational stress-induced hypothesis can explain most regions. However, structural controls cannot be ruled out for any region except polygon 2. Figure 14. Open in new tabDownload slide Quantitative comparison of the possible source of crustal anisotropy around the study area. The bar chart is grouped by the delineated regions and in each region the proposed sources of crustal anisotropy (colour coded in the legend) is compared. Green, blue, red and yellow bars represent gravitational stresses, active fault orientation, SHmax orientation from focal mechanism and maximum shear plane orientation, respectively. The thick vertical lines on each bar represent the mean error of the residuals. 6 CONCLUSIONS We have analysed a very large number of measurements of high-quality SWS parameters (∼102 000) around the Marlborough and Wellington region with a detailed systematic analysis of crustal anisotropy using over 40 000 crustal earthquakes recorded on at least one of the 36 stations. The size of our data set has allowed us to show the detailed spatial and temporal variation of anisotropic structure around central New Zealand. We do not observe any spatiotemporal variation of the state of anisotropy over the |$5\frac{1}{2}$| yr period. Apparent temporal variations of δt and VP/VS are mainly due to spatial variations of anisotropy combined with temporally varying paths of seismic waves. We conclude that the crustal anisotropy around the study region is confined to the upper few kilometres of the crust and can be controlled by either one mechanism or a combination of more than one. The high correspondence between SHmax calculated from gravitational potential energy from topography and average fast polarization orientation around the Marlborough and Kaikōura region suggest that gravitationally induced stresses contribute to the crustal anisotropy around the region. In the North Island tectonic stresses inferred from focal mechanisms as well as geological structures could contribute to the observed anisotropy. In the South Island gravitational potential energy gives the best overall match, but inherited geological structures could contribute and control some regions. We suggest that examining the effect of gravitational stresses on crustal seismic anisotropy should not be neglected in future studies. SUPPORTING INFORMATION Figure S1 Map of events (a) before and (b) after the 2016 Kaikōura earthquake: map view with epicenters of earthquakes before (left) and after (right) the Mw 7.8 Kaikōura Earthquake. Red star shows the Mw 7.8 epicentre with the red triangles showing the GeoNet stations in the region. Events are scaled according to their magnitude, and colour coded according to their depth. A cross-section from a to a1 shows the depth distribution and coloured by time (timescale on main figure). Figure S2 Plot of data continuity for all stations used for analysis over the five and half year period. Red bars show the period during which each station was operational. Figure S3 1-D models used in the analysis. The model consists of an anisotropic layer between isotropic layers. p = density, Vp and Vs are P and S velocity, respectively. Theta and phi defines the axis of symmetry. phi is the fast direction. Model A is the baseline model. In model B, we test the effect of incidence angle on splitting results. We used 1-D reflectivity codes Levin & Park (1997) for the waveform simulation in 1-D anisotropic media and MFAST was used to estimate SWS parameters. aninc, i, represents the incidence angle. A delay time of 0.3 s was implemented in the model. Baz is the backazimuth and alpha is the incoming polarization on the shear wave. Figure S4 Comparison of apparent splitting parameters as a function of backazimuth ranging from 0° to 360° with 5° increments. As angle of incidence is increased, there is more variation in δt and ϕ with varying backazimuth. This is because phase conversion, at high angles, may interfere with SWS measurements. Figure S5 Example of generalized phase detection GPD pick. The red and blue colours indicate the probability of P and S waves, respectively. The vertical bars are estimated picks. This is an example of a reliable P and S picks for station CVR. Figure S6 Distributions of time difference between GeoNet manual and GDP S picks. The distribution is normally distributed with mean −0.072 and standard deviation 0.774. Figure S7 Study area is divided into grid cell using quadtree gridding. Red lines show ray paths, and blue triangles are seismic stations. For each grid, a minimum of 20 and maximum of 80 ray paths pass through each grid cell, until a minimum grid size of 5 × 5 km is recorded. Figure S8 Spatial averaging results of ϕ for before (top row) and after (bottom row) Kaikōura earthquake. (a) Study area is divided into grid cells using quadtree gridding. Red lines show ray paths, and blue inverted triangles are seismic stations. A minimum of 20 and maximum of 80 ray paths pass through each block, until a minimum block size of 5 km is recorded. (b) Spatial averages with weighting inversely proportional to the square of the distance from the station. The coloured bars show the mean in each grid. Green line are the fault that ruptured during the Kaikōura main shock. Figure S9 Time-series of delay times, δt and Vp/Vs ratio for values from 2013 January to 2018 June is shown in each panel. a and a1 are plots for all and clustered measurements for station CAW respectively. Similarly, b and b1 are plots for station KIW. Blue, red and orange vertical bars represent the main shock of the 2013 Mw 6.5 2013 Seddon, the 2013 Mw 6.6 Lake Grassmere, the 2014 Mw 6.2 Eketahuna and the 2016 Mw 7.8 Kaikōura earthquakes, respectively. Green horizontal line shows the median value. Green dashed line represents the estimated error using bootstrap. Red and blue lines are the 20-d moving median of Vp/Vs ratio and δt, respectively. Figure S10 Delay times versus path length for station KHZ. Distribution of the delay time and path length are plotted on the right and bottom, respectively. The scattered points are coloured by the depths and plotted on a basemap of the kernel density (effectively heatmap). High-density regions are shaded white with low regions black. The red line is the least-squares regression line with the equation and p-value displayed on the lower right corner of each plot. Figure S11 Quantitative comparison between spatially averaged ϕ measurements and the orientation of the closest (<30 km) SHmax. The background colour is a contour of the absolute value of the cosine of the difference between the orientations, coloured by the scale at the bottom. Note that 1 represents no difference between the orientations. The descriptive statistics of the difference between orientations at the closest gridpoints is displayed on the histogram in (b). SHmax are represented by dark green bow ties. The histogram of the SHmax depths used for analysis is shown in (c). Table S1 Summary of adjusted parameters used for picking S arrivals. Descriptions are after Castellazzi et al. (2015) and Diehl et al. (2009). 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 paper. ACKNOWLEDGEMENTS This work is supported by funding from Victoria University of Wellington (PhD scholarship) and the Earthquake commission (EQC grant number: 18/755). We are grateful to Francesco Pio Lucente and an anonymous reviewer for their suggestions and comments in making this manuscript better. We also are grateful to Simon Lamb, Carolyn Boulton, Katrina Jacobs and Hamish Hirschberg for their discussions and helpful feedback. Waveform data and earthquake catalogues used for this study were provided by New Zealand GeoNet and are available from https://quakesearch.geonet.org.nz/. Most figures were prepared using the Generic Mapping Tools (Wessel et al. 2013). A table of our SWS measurements can be found at https://zenodo.org/record/3738129. Codes used to estimate the orientation and magnitude of gravitationally induced stresses can be found at https://github.com/hamishhirschberg/stress-modelling. Codes used for SWS parameters can be found at https://mfast-package.geo.vuw.ac.nz. REFERENCES Abercrombie R.E. , Bannister S., Ristau J., Doser D., 2016 . Variability of earthquake stress drop in a subduction setting, the Hikurangi Margin, New Zealand , Geophys. J. Int. , 208 ( 1 ), 306 – 320 . 10.1093/gji/ggw393 Google Scholar Crossref Search ADS WorldCat Allen R.V. , 1978 . Automatic earthquake recognition and timing from single traces , Bull. seism. Soc. Am. , 68 ( 5 ), 1521 – 1532 . Google Scholar OpenURL Placeholder Text WorldCat Anderson E.M. , 1905 . The dynamics of faulting , Trans. Edinburgh Geol. Soc. , 8 ( 3 ), 387 – 402 . 10.1144/transed.8.3.387 Google Scholar Crossref Search ADS WorldCat Anderson H. , Webb T., Jackson J., 1993 . Focal mechanisms of large earthquakes in the South Island of New Zealand: implications for the accommodation of Pacific-Australia plate motion , Geophys. J. Int. , 115 ( 3 ), 1032 – 1054 . 10.1111/j.1365-246X.1993.tb01508.x Google Scholar Crossref Search ADS WorldCat Araragi K.R. , Savage M.K., Ohminato T., Aoki Y., 2015 . Seismic anisotropy of the upper crust around Mount Fuji, Japan , J. geophys. Res.: Solid Earth , 120 ( 4 ), 2739 – 2751 . 10.1002/2014JB011554 Google Scholar Crossref Search ADS WorldCat Arnold R. , Townend J., 2007 . A Bayesian approach to estimating tectonic stress from seismological data , Geophys. J. Int. , 170 ( 3 ), 1336 – 1356 . 10.1111/j.1365-246X.2007.03485.x Google Scholar Crossref Search ADS WorldCat Aster R.C. , Shearer P.M., Berger J., 1990 . Quantitative measurements of shear wave polarizations at the Anza Seismic Network, southern California: implications for shear wave splitting and earthquake prediction , J. geophys. Res.: Solid Earth , 95 ( B8 ), 12449 – 12473 . 10.1029/JB095iB08p12449 Google Scholar Crossref Search ADS WorldCat Audoine E. , Savage M.K., Gledhill K., 2000 . 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 , 105 ( B4 ), 8013 – 8033 . Google Scholar Crossref Search ADS WorldCat Babuska V. , Cara M., 1991 . Seismic Anisotropy in the Earth , Vol. 10 , Springer Science & Business Media . Google Scholar Crossref Search ADS Google Preview WorldCat COPAC Balfour N. , Savage M., Townend J., 2005 . Stress and crustal anisotropy in Marlborough, New Zealand: evidence for low fault strength and structure-controlled anisotropy , Geophys. J. Int. , 163 ( 3 ), 1073 – 1086 . Google Scholar Crossref Search ADS WorldCat Balfour N. , Cassidy J.F., Dosso S.E., 2012 . Crustal anisotropy in the forearc of the Northern Cascadia Subduction Zone, British Columbia , Geophys. J. Int. , 188 ( 1 ), 165 – 176 . Google Scholar Crossref Search ADS WorldCat Berens P. , 2009 . CircStat: a MATLAB Toolbox for circular statistics , J. Stat. Softw. , 31 ( 10 ), 10.18637/jss.v031.i10 . Google Scholar OpenURL Placeholder Text WorldCat Bianco F. , Scarfı L., Pezzo E.D., Patanè D., 2006 . Shear wave splitting changes associated with the 2001 volcanic eruption on Mt Etna , Geophys. J. Int. , 167 ( 2 ), 959 – 967 . 10.1111/j.1365-246X.2006.03152.x Google Scholar Crossref Search ADS WorldCat Boness N. , Zoback M., 2006 . Mapping stress and structurally controlled crustal shear velocity anisotropy in California , Geology , 34 ( 10 ), 825 – 828 . 10.1130/G22309.1 Google Scholar Crossref Search ADS WorldCat Castellazzi C. , Savage M., Walsh E., Arnold R., 2015 . Shear wave automatic picking and splitting measurements at Ruapehu volcano, New Zealand , J. geophys. Res.: Solid Earth , 120 , 3363 – 3384 . 10.1002/2014JB011585 Google Scholar Crossref Search ADS WorldCat Cesca S. et al. , 2017 . Complex rupture process of the Mw 7.8, 2016, Kaikōura earthquake, New Zealand, and its aftershock sequence , Earth planet. Sci. Lett. , 478 , 110 – 120 . 10.1016/j.epsl.2017.08.024 Google Scholar Crossref Search ADS WorldCat Chen Y. , Zhang Z., Sun C., Badal J., 2013 . Crustal anisotropy from Moho converted Ps wave splitting analysis and geodynamic implications beneath the eastern margin of Tibet and surrounding regions , Gondwana Res. , 24 ( 3–4 ), 946 – 957 . 10.1016/j.gr.2012.04.003 Google Scholar Crossref Search ADS WorldCat Chiarabba C. , Gori P.D., Boschi E., 2009 . Pore-pressure migration along a normal-fault system resolved by time-repeated seismic tomography , Geology , 37 ( 1 ), 67 – 70 . 10.1130/G25220A.1 Google Scholar Crossref Search ADS WorldCat Christensen N.I. , 1996 . Poisson’s ratio and crustal seismology , J. geophys. Res.: Solid Earth , 101 ( B2 ), 3139 – 3156 . 10.1029/95JB03446 Google Scholar Crossref Search ADS WorldCat Ciucu C. , 2010 . Maximum stress at a tectonic fault plane. Coulomb Law , Rom. Rep. Phys. , 62 ( 2 ), 414 – 421 . Google Scholar OpenURL Placeholder Text WorldCat Clark K. et al. , 2017 . Highly variable coastal deformation in the 2016 Mw 7.8 Kaikōura earthquake reflects rupture complexity along a transpressional plate boundary , Earth planet. Sci. Lett. , 474 , 334 – 344 . 10.1016/j.epsl.2017.06.048 Google Scholar Crossref Search ADS WorldCat Cochran E. , Kroll K., 2015 . Stress- and structure-controlled anisotropy in a region of complex faulting—Yuha Desert, California , Geophys. J. Int. , 202 ( 2 ), 1109 – 1121 . 10.1093/gji/ggv191 Google Scholar Crossref Search ADS WorldCat Cochran E. , Li Y.-G., Vidale J., 2006 . Anisotropy in the shallow crust observed around the San Andreas fault before and after the 2004 M 6.0 Parkfield earthquake , Bull. seism. Soc. Am. , 96 ( 4B ), S364 – S375 . 10.1785/0120050804 Google Scholar Crossref Search ADS WorldCat Crampin S. , 1984 . An introduction to wave propagation in anisotropic media , Geophys. J. Int. , 76 ( 1 ), 17 – 28 . 10.1111/j.1365-246X.1984.tb05018.x Google Scholar Crossref Search ADS WorldCat Crampin S. , 1990 . The scattering of shear-waves in the crust , Pure appl. Geophys. , 132 ( 1–2 ), 67 – 91 . 10.1007/BF00874358 Google Scholar Crossref Search ADS WorldCat Crampin S. , Booth D.C., 1989 . Shear-wave splitting showing hydraulic dilation of pre-existing joints in granite , Sci. Drill. , 1 ( 1 ), 21 – 26 . Google Scholar OpenURL Placeholder Text WorldCat Crampin S. , Gao Y., 2006 . A review of techniques for measuring shear-wave splitting above small earthquakes , Phys. Earth planet. Inter. , 159 ( 1–2 ), 1 – 14 . 10.1016/j.pepi.2006.06.002 Google Scholar Crossref Search ADS WorldCat Crampin S. , Gao Y., Bukits J., 2015 . A review of retrospective stress-forecasts of earthquakes and eruptions , Phys. Earth planet. Inter. , 245 , 76 – 87 . 10.1016/j.pepi.2015.05.008 Google Scholar Crossref Search ADS WorldCat Crotwell H.P. , Owens T.J., Ritsema J., 1999 . The TauP Toolkit: flexible seismic travel-time and ray-path utilities , Seismol. Res. Lett. , 70 ( 2 ), 154 – 160 . 10.1785/gssrl.70.2.154 Google Scholar Crossref Search ADS WorldCat Dellow S. et al. , 2017 . Landslides caused by the 14 November 2016 Mw7.8 Kaikōura earthquake and the immediate response , Bull. New Zealand Soc. Earthq. Eng. , 50 ( 2 ), 106 – 116 . 10.5459/bnzsee.50.2.106-116 Google Scholar Crossref Search ADS WorldCat Diehl T. , Deichmann N., Kissling E., Husen S., 2009 . Automatic S-wave picker for local earthquake tomography , Bull. seism. Soc. Am. , 99 ( 3 ), 1906 – 1920 . 10.1785/0120080019 Google Scholar Crossref Search ADS WorldCat Duputel Z. , Rivera L., 2017 . Long-period analysis of the 2016 Kaikōura earthquake , Phys. Earth planet. Inter. , 265 , 62 – 66 . 10.1016/j.pepi.2017.02.004 Google Scholar Crossref Search ADS WorldCat Dvorkin J. , Mavko G., Nur A., 1999 . Overpressure detection from compressional- and shear-wave data , Geophys. Res. Lett. , 26 ( 22 ), 3417 – 3420 . 10.1029/1999GL008382 Google Scholar Crossref Search ADS WorldCat Eberhart-Phillips D. , Fry B., 2018 . Joint local earthquake and teleseismic inversion for 3-D velocity and Q in New Zealand , Phys. Earth planet. Inter. , 283 , 48 – 66 . 10.1016/j.pepi.2018.08.005 Google Scholar Crossref Search ADS WorldCat Evanzia D. , Wilson T., Savage M.K., Lamb S., Hirschberg H., 2017 . Stress orientations in a locked subduction zone at the Southern Hikurangi Margin, New Zealand , J. geophys. Res.: Solid Earth , 122 ( 10 ), 7895 – 7911 . 10.1002/2017JB013998 Google Scholar Crossref Search ADS WorldCat Flesch L.M. , Haines A.J., Holt W.E., 2001 . Dynamics of the India-Eurasia collision zone , J. geophys. Res.: Solid Earth , 106 ( B8 ), 16435 – 16460 . 10.1029/2001JB000208 Google Scholar Crossref Search ADS WorldCat Flinn E. , 1965 . Signal analysis using rectilinearity and direction of particle motion , Proc. IEEE , 53 ( 12 ), 1874 – 1876 . 10.1109/PROC.1965.4462 Google Scholar Crossref Search ADS WorldCat GeoNet. 2019 . (accessed September 24, 2020) . GeoNet quake search catalog. Available at: https://quakesearch.geonet.org.nz/ . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Gledhill K.R. , 1991 . Evidence for shallow and pervasive seismic anisotropy in the Wellington Region, New Zealand , J. geophys. Res.: Solid Earth , 96 ( B13 ), 21503 – 21516 . 10.1029/91JB02049 Google Scholar Crossref Search ADS WorldCat Hamling I.J. , D’Anastasio E., Wallace L.M., Ellis S., Motagh M., Samsonov S., Palmer N., Hreinsdóttir S., 2014 . Crustal deformation and stress transfer during a propagating earthquake sequence: the 2013 Cook Strait sequence, central New Zealand , J. geophys. Res.: Solid Earth , 119 ( 7 ), 6080 – 6092 . 10.1002/2014JB011084 Google Scholar Crossref Search ADS WorldCat Hamling I.J. et al. , 2017 . Complex multifault rupture during the 2016 Mw 7.8 Kaikōura earthquake, New Zealand , Science , 356 ( 6334 ), 1 – 10 . 10.1126/science.aam7194 Google Scholar Crossref Search ADS WorldCat Hardebeck J.L. , Michael A.J., 2006 . Damped regional-scale stress inversions: methodology and examples for southern California and the Coalinga aftershock sequence , J. geophys. Res.: Solid Earth , 111 ( B11 ), 10.1029/2005jb004144 . 10.1029/2005jb004144 Google Scholar OpenURL Placeholder Text WorldCat Crossref Hardebeck J.L. , Okada T., 2018 . Temporal stress changes caused by earthquakes: A review , J. geophys. Res.: Solid Earth , 123 ( 2 ), 1350 – 1365 . 10.1002/2017JB014617 Google Scholar Crossref Search ADS WorldCat Healy D. , Sibson R.H., Shipton Z., Butler R., 2012 . Stress, faulting, fracturing and seismicity: the legacy of ernest masson anderson , Geol. Soc., Lond., Spec. Publ. , 367 ( 1 ), 1 – 6 . 10.1144/SP367.1 Google Scholar Crossref Search ADS WorldCat Hirschberg H.P. , Lamb S., Savage M.K., 2018 . Strength of an obliquely convergent plate boundary: lithospheric stress magnitudes and viscosity in New Zealand , Geophys. J. Int. , 216 ( 2 ), 1005 – 1024 . Google Scholar Crossref Search ADS WorldCat Holden C. , Kaiser A., Van Dissen R., Jury R., 2013 . Sources, ground motion and structural response characteristics in Wellington of the 2013 Cook Strait earthquakes , Bull. New Zeal. Soc. Earthq. Eng. , 46 ( 4 ), 188 – 195 . Google Scholar Crossref Search ADS WorldCat Holden C. , Kaneko Y., D’Anastasio E., Benites R., Fry B., Hamling I., 2017 . The 2016 Kaikōura Earthquake revealed by kinematic source inversion and seismic wavefield simulations: slow rupture propagation on a geometrically complex crustal fault network , Geophys. Res. Lett. , 44 ( 22 ), 10.1002/2017GL075301 . Google Scholar OpenURL Placeholder Text WorldCat Holt R. , Savage M., Townend J., Syracuse E., Thurber C., 2013 . Crustal stress and fault strength in the Canterbury Plains, New Zealand , Earth planet. Sci. Lett. , 383 , 173 – 181 . Google Scholar Crossref Search ADS WorldCat Iidaka T. , Obara K., 2013 . Shear-wave splitting in a region with newly-activated seismicity after the 2011 tohoku earthquake , Earth Planets Space , 65 ( 9 ), 1059 – 1064 . Google Scholar Crossref Search ADS WorldCat Illsley-Kemp F. et al. , 2017 . Extension and stress during continental breakup: Seismic anisotropy of the crust in Northern Afar , Earth planet. Sci. Lett. , 477 , 41 – 51 . Google Scholar Crossref Search ADS WorldCat Ito H. , DeVilbiss J., Nur A., 1979 . Compressional and shear waves in saturated rock during water-steam transition , J. geophys. Res. , 84 ( B9 ), 4731 – 4735 . 10.1029/jb084ib09p04731 Google Scholar Crossref Search ADS WorldCat Johnson J.H. , Savage M.K., Townend J., 2011 . Distinguishing between stress-induced and structural anisotropy at Mount Ruapehu volcano, New Zealand , J. geophys. Res. , 116 , 10.1029/2011JB008308 . 10.1029/2011JB008308 Google Scholar OpenURL Placeholder Text WorldCat Crossref Kaiser A. et al. , 2017 . The 2016 Kaikōura, New Zealand, Earthquake: preliminary seismological report , Seismol. Res. Lett. , 88 ( 3 ), 727 – 739 . 10.1785/0220170018 Google Scholar Crossref Search ADS WorldCat Karalliyadda S.C. , Savage M.K., 2013 . Seismic anisotropy and lithospheric deformation of the plate-boundary zone in South Island, New Zealand: inferences from local S-wave splitting , Geophys. J. Int. , 193 ( 2 ), 507 – 530 . Google Scholar Crossref Search ADS WorldCat Kisslinger C. , Engdahl E., 1973 . The interpretation of the Wadati diagram with relaxed assumptions , Bull. seism. Soc. Am. , 63 ( 5 ), 1723 – 1736 . Google Scholar OpenURL Placeholder Text WorldCat Klosko E.R. , Wu F.T., Anderson H.J., Eberhart-Phillips D., McEvilly T.V., Audoine E., Savage M.K., Gledhill K.R., 1999 . Upper mantle anisotropy in the New Zealand Region , Geophys. Res. Lett. , 26 ( 10 ), 1497 – 1500 . Google Scholar Crossref Search ADS WorldCat Kong Q. , Trugman D.T., Ross Z.E., Bianco M.J., Meade B.J., Gerstoft P., 2018 . Machine learning in seismology: turning data into insights , Seismol. Res. Lett. , 90 ( 1 ), 3 – 14 . Google Scholar Crossref Search ADS WorldCat Lamb S. , Arnold R., Moore J. D.P., 2018 . Locking on a megathrust as a cause of distributed faulting and fault-jumping earthquakes , Nat. Geosci. , 11 ( 11 ), 871 – 875 . Google Scholar Crossref Search ADS WorldCat Langridge RM. et al. , 2016 . The New Zealand active faults database , New Zealand J. Geol. Geophys. , 59 ( 1 ), 86 – 96 . 10.1080/00288306.2015.1112818 Google Scholar Crossref Search ADS WorldCat Lanza F. et al. , 2019 . Crustal fault connectivity of the mw 7.8 2016 kaikōura earthquake constrained by aftershock relocations , Geophys. Res. Lett. , 46 ( 12 ), 6487 – 6496 . 10.1029/2019GL082780 Google Scholar Crossref Search ADS WorldCat Levin V. , Park J., 1997 . P-SH conversions in a flat-layered medium with anisotropy of arbitrary orientation , Geophys. J. Int. , 253 – 266 . Google Scholar OpenURL Placeholder Text WorldCat Li Y.-G. , Vidale J.E., 2001 . Healing of the shallow fault zone from 1994-1998 After the 1992 M7.5 Landers, California, Earthquake , Geophys. Res. Lett. , 28 ( 15 ), 2999 – 3002 . Google Scholar Crossref Search ADS WorldCat Litchfield N.J. , Benson A ., Bischoff A ., Hatem A., Wandres A., 2017 . Multiple fault ground surface ruptures in the 14 November 2016 Kaikōura Earthquake, New Zealand , Seismol. Res. Lett. , 88 ( 2B ), 624 – 626 . Google Scholar OpenURL Placeholder Text WorldCat Liu Y. , 2004 . Systematic analysis of shear-wave splitting in the aftershock zone of the 1999 Chi-Chi, Taiwan, Earthquake: shallow crustal anisotropy and lack of precursory variations , Bull. seism. Soc. Am. , 94 ( 6 ), 2330 – 2347 . Google Scholar Crossref Search ADS WorldCat Lucente F.P. , Gori P.D., Margheriti L., Piccinini D., Bona M.D., Chiarabba C., Agostinetti N.P., 2010 . Temporal variation of seismic velocity and anisotropy before the 2009 MW6.3 L’Aquila earthquake, Italy , Geology , 38 ( 11 ), 1015 – 1018 . Google Scholar Crossref Search ADS WorldCat Matcham I. , Savage M.K., Gledhill K.R., 2000 . Distribution of seismic anisotropy in the subduction zone beneath the Wellington region, New Zealand , Geophys. J. Int. , 140 ( 1 ), 1 – 10 . Google Scholar Crossref Search ADS WorldCat Michael A.J. , 1984 . Determination of stress from slip data: Faults and folds , J. geophys. Res.: Solid Earth , 89 ( B13 ), 11517 – 11526 . Google Scholar Crossref Search ADS WorldCat Miller V. , Savage M., 2001 . Changes in seismic anisotropy after volcanic eruptions: evidence from Mount Ruapehu , Science , 293 ( 5538 ), 2231 – 2233 . Google Scholar Crossref Search ADS PubMed WorldCat Mitchell J. , Mackay K., Neil H., Mackay E., Pallentin A., Notman P., 2012 . New Zealand Region Bathymetry, 1:4,000000 . Mortimer N. , 1993 . Metamorphic zones, terranes, and Cenozoic faults in the Marlborough Schist, New Zealand , New Zealand J. Geol. Geophys. , 36 ( 3 ), 357 – 368 . Google Scholar Crossref Search ADS WorldCat Mortimer N. , 2004 . New Zealand’s geological foundations , Gondwana Res. , 7 ( 1 ), 261 – 272 . Google Scholar Crossref Search ADS WorldCat Nakata N. , Snieder R., 2012 . Time-lapse change in anisotropy in japan’s near surface after the 2011 tohoku-oki earthquake , Geophys. Res. Lett. , 39 ( 11 ), 10.1029/2012gl051979 . 10.1029/2012gl051979 Google Scholar OpenURL Placeholder Text WorldCat Crossref Neuberg J. , Pointer T., 2000 . Effects of volcano topography on seismic broad-band waveforms , Geophys. J. Int. , 143 ( 1 ), 239 – 248 . Google Scholar Crossref Search ADS WorldCat Nur A. , 1971 . Effects of stress on velocity anisotropy in rocks with cracks , J. geophys. Res. , 76 ( 8 ), 2022 – 2034 . Google Scholar Crossref Search ADS WorldCat Nuttli O. , 1961 . The effect of the earth’s surface on the S wave particle motion , Bull. seism. Soc. Am. , 51 ( 2 ), 237 – 246 . Google Scholar OpenURL Placeholder Text WorldCat Okada T. et al. , 2019 . Comparative tomography of reverse-slip and strike-slip seismotectonic provinces in the northern South Island, New Zealand , Tectonophysics , 765 , 172 – 186 . 10.1016/j.tecto.2019.03.016 Google Scholar Crossref Search ADS WorldCat Okaya D. , Christensen N., Stanley D., Stern T., (term) S.I.G.T., 1995 . Crustal anisotropy in the vicinity of the Alpine Fault Zone, South Island, New Zealand , New Zealand J. Geol. Geophys. , 38 ( 4 ), 579 – 583 . Google Scholar Crossref Search ADS WorldCat Okaya D. , Christensen N.I., Ross Z.E., Wu F.T., 2016 . Terrane-controlled crustal shear wave splitting in Taiwan , Geophys. Res. Lett. , 43 ( 2 ), 556 – 563 . Google Scholar Crossref Search ADS WorldCat Peng Z. , Ben-Zion Y., 2004 . Systematic analysis of crustal anisotropy along the Karadere—Düzce branch of the North Anatolian fault , Geophys. J. Int. , 159 ( 1 ), 253 – 274 . Google Scholar Crossref Search ADS WorldCat Peng Z. , Ben-Zion Y., 2005 . Spatio-temporal variations of crustal anisotropy from similar events in aftershocks of the 1999 M7.4 Izmit and M7.1 Duzce, Turkey, earthquake sequences , Geophys. J. Int. , 160 , 1027 – 1043 . Google Scholar Crossref Search ADS WorldCat Petersen T. , Gledhill K., Chadwick M., Gale N.H., Ristau J., 2011 . The new zealand national seismograph network , Seismol. Res. Lett. , 82 ( 1 ), 9 – 20 . Google Scholar Crossref Search ADS WorldCat Pettinga J. , Yetton M., Van Dissen R., Downes G., 2001 . Earthquake source identification and characterisation for the Canterbury Region, South Island, New Zealand , Bull. New Zealand Soc. Earthq. Eng. , 34 ( 4 ), 282 – 317 . Google Scholar Crossref Search ADS WorldCat Piccinini D. , Margheriti L., Chiaraluce L., Cocco M., 2006 . Space and time variations of crustal anisotropy during the 1997 Umbria-Marche, central Italy, seismic sequence , Geophys. J. Int. , 167 ( 3 ), 1482 – 1490 . Google Scholar Crossref Search ADS WorldCat Reyners M. , Cowan H., 1993 . The transition from subduction to continental collision: crustal structure in the North Canterbury region, New Zealand , Geophys. J. Int. , 115 ( 3 ), 1124 – 1136 . Google Scholar Crossref Search ADS WorldCat Ristau J. , 2013 . Update of regional moment tensor analysis for earthquakes in New Zealand and adjacent offshore regions , Bull. seism. Soc. Am. , 103 ( 4 ), 2520 – 2533 . Google Scholar Crossref Search ADS WorldCat Roman D.C. , Savage M.K., Arnold R., Latchman J.L., De Angelis S., 2011 . Analysis and forward modeling of seismic anisotropy during the ongoing eruption of the Soufrière Hills Volcano, Montserrat, 1996–2007 , J. geophys. Res.: Solid Earth , 116 ( B3 ), 10.1029/2010jb007667 . 10.1029/2010jb007667 Google Scholar OpenURL Placeholder Text WorldCat Crossref Ross Z.E. , Meier M., Hauksson E., Heaton T.H., 2018 . Generalized seismic phase detection with deep learning , Bull. seism. Soc. Am. , 108 ( 5A ), 2894 – 2901 . Google Scholar Crossref Search ADS WorldCat Rümpker G. , Silver P.G., 1998 . Apparent shear-wave splitting parameters in the presence of vertically varying anisotropy , Geophys. J. Int. , 135 ( 3 ), 790 – 800 . Google Scholar Crossref Search ADS WorldCat Saiga A. , Hiramatsu Y., Ooida T., Yamaoka K., 2003 . Spatial variation in the crustal anisotropy and its temporal variation associated with a moderate-sized earthquake in the Tokai region, central Japan , Geophys. J. Int. , 154 ( 3 ), 695 – 705 . Google Scholar Crossref Search ADS WorldCat Savage M. , 1999 . Seismic anisotropy and mantle deformation: what have we learned from shear wave splitting? , Rev. Geophys. , 37 , 65 – 106 . Google Scholar Crossref Search ADS WorldCat Savage M. , Shih X., Meyer R., Aster R., 1989 . Shear-wave anisotropy of active tectonic regions via automated S-wave polarization analysis , Tectonophysics , 165 ( 1–4 ), 279 – 292 . Google Scholar Crossref Search ADS WorldCat Savage M. , Wessell A., Teanby N., Hurst A., 2010a . Automatic measurement of shear wave splitting and applications to time varying anisotropy at Mount Ruapehu volcano, New Zealand , J. geophys. Res. , 115 , 115 , 10.1029/2010jb007722 . 10.1029/2010jb007722 Google Scholar OpenURL Placeholder Text WorldCat Crossref Savage M. , Ohminato T., Aoki Y., Tsuji H., Greve S.M., 2010b . Stress magnitude and its temporal variation at Mt. Asama Volcano, Japan, from seismic anisotropy and GPS , Earth planet. Sci. Lett. , 290 ( 3 ), 403 – 414 . Google Scholar Crossref Search ADS WorldCat Savage M.K. et al. , 2016 . Stress, strain rate and anisotropy in Kyushu, Japan , Earth planet. Sci. Lett. , 439 , 129 – 142 . 10.1016/j.epsl.2016.01.005 Google Scholar Crossref Search ADS WorldCat Shih X.R. , Meyer R.P., 1990 . Observation of shear wave splitting from natural events: South moat of long valley caldera, california, june 29 to august 12, 1982 , J. geophys. Res.: Solid Earth , 95 ( B7 ), 11179 – 11195 . Google Scholar Crossref Search ADS WorldCat Silver P.G. , Chan W.W., 1991 . Shear wave splitting and subcontinental mantle deformation , J. geophys. Res.: Solid Earth , 96 ( B10 ), 16429 – 16454 . Google Scholar Crossref Search ADS WorldCat Tadokoro K. , Ando M., 2002 . Evidence for rapid fault healing derived from temporal changes in S wave splitting , Geophys. Res. Lett. , 29 ( 4 ), 6 – 1-6-4 . Google Scholar Crossref Search ADS WorldCat Tadokoro K. , Ando M., Umeda Y., 1999 . Swave splitting in the aftershock region of the 1995 Hyogo-ken Nanbu earthquake , J. geophys. Res.: Solid Earth , 104 ( B1 ), 981 – 991 . Google Scholar Crossref Search ADS WorldCat Takanami T. , Kitagawa G., 1988 . A new efficient procedure for the estimation of onset times of seismic waves , J. Phys. Earth , 36 ( 6 ), 267 – 290 . Google Scholar Crossref Search ADS WorldCat Teanby N. , Kendall J.-M., Jones R.H., Barkved O., 2004a . Stress-induced temporal variations in seismic anisotropy observed in microseismic data , Geophys. J. Int. , 156 ( 3 ), 459 – 466 . Google Scholar Crossref Search ADS WorldCat Teanby N. , Kendall J.-M., Jones R.H., Barkved O., 2004b . Automation of shear-wave splitting measurements using cluster analysis , Bull. seism. Soc. Am. , 94 ( 2 ), 453 – 463 . Google Scholar Crossref Search ADS WorldCat Townend J. , Zoback M.D., 2000 . How faulting keeps the crust strong , Geology , 28 ( 5 ), 399 . Google Scholar Crossref Search ADS WorldCat Townend J. , Zoback M.D., 2001 . Implications of earthquake focal mechanisms for the frictional strength of the San Andreas fault system , Geol. Soc., Lond., Spec. Publ. , 186 ( 1 ), 13 – 21 . Google Scholar Crossref Search ADS WorldCat Townend J. , Sherburn S., Arnold R., Boese C., Woods L., 2012 . Three-dimensional variations in present-day tectonic stress along the Australia–Pacific plate boundary in New Zealand , Earth planet. Sci. Lett. , 353-354 , 47 – 59 . Google Scholar Crossref Search ADS WorldCat Unglert K. , Savage M.K., Fournier N., Ohkura T., Abe Y., 2011 . Shear wave splitting, Vp/Vs, and GPS during a time of enhanced activity at Aso caldera, Kyushu , J. geophys. Res.: Solid Earth , 116 ( B11 ), 10.1029/2011jb008520 . 10.1029/2011jb008520 Google Scholar OpenURL Placeholder Text WorldCat Crossref Van Dissen R. , Yeats R.S., 1991 . Hope fault, Jordan thrust, and uplift of the Seaward Kaikōura Range, New Zealand , Geology , 19 ( 4 ), 393 , 10.1130/0091-7613(1991)019<0393:HFJTAU>2.3.CO;2 . 10.1130/0091-7613(1991)019&lt;0393:HFJTAU&gt;2.3.CO;2 Google Scholar Crossref Search ADS WorldCat Wadati K. , Oki S., 1933 . On the travel time of earthquake waves. (Part II) , J. Meteorol. Soc. Japan. Ser. II , 11 ( 1 ), 14 – 28 . Google Scholar Crossref Search ADS WorldCat Wallace L.M. , Kaneko Y., Hreinsdóttir S., Hamling I., Peng Z., Bartlow N., D’Anastasio E., Fry B., 2017 . Large-scale dynamic triggering of shallow slow slip enhanced by overlying sedimentary wedge , Nat. Geosci. , 10 ( 10 ), 765 – 770 . 10.1038/NGEO3021 Google Scholar Crossref Search ADS WorldCat Walsh E. , 2012 . Measuring shear wave splitting using the Silver and Chan method , Master’s thesis , Victoria University of Wellington . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Walsh E. , Arnold R., Savage M.K., 2013 . Silver and Chan revisited , J. geophys. Res.: Solid Earth , 118 ( 10 ), 5500 – 5515 . Google Scholar Crossref Search ADS WorldCat Wessel A. , 2010 . Automatic shear wave splitting measurements at Mt. Ruapehu volcano, New Zealand , Master’s thesis , Victoria University of Wellington . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC Wessel P. , Smith W., Scharroo R., Luis J., Wobbe F., 2013 . Generic Mapping Tools: improved version released , EOS, Trans. Am. geophys. Un. , 94 ( 45 ), 409 – 410 . Google Scholar Crossref Search ADS WorldCat Wilson C.K. , Jones C.H., Molnar P., Sheehan A.F., Boyd O.S., 2004 . Distributed deformation in the lower crust and upper mantle beneath a continental strike-slip fault zone: Marlborough fault system, south island, New Zealand , Geology , 32 ( 10 ), 837 – 840 .. 10.1130/g20657.1 . 10.1130/g20657.1 Google Scholar Crossref Search ADS WorldCat Yang Z. , Sheehan A., Shearer P., 2011 . Stress-induced upper crustal anisotropy in southern California , J. geophys. Res.: Solid Earth , 116 ( B2 ), 1 – 11 . Google Scholar Crossref Search ADS WorldCat Zal H.J. et al. , 2020 . Temporal and spatial variations in seismic anisotropy and Vp/Vs ratios in a region of slow slip , Earth planet. Sci. Lett. , 532 , 115970 . https://doi.org/10.1016/j.epsl.2019.115970 . 10.1016/j.epsl.2019.115970 Google Scholar Crossref Search ADS WorldCat Zhang H. , Liu Y., Thurber C., Roecker S., 2007 . Three-dimensional shear-wave splitting tomography in the Parkfield, California, region , Geophys. Res. Lett. , 34 ( 24 ), 10.1029/2007GL031951 . 10.1029/2007GL031951 Google Scholar OpenURL Placeholder Text WorldCat Crossref Zheng X.-F. , Chen C.-H., Zhang C.-H., 2008 . Study on temporal variations of shear-wave splitting in the Chiayi Area, aftershock Zone of 1999 Chichi Earthquake, Taiwan , Chin. J. Geophys. , 51 ( 1 ), 115 – 124 . 10.1002/cjg2.1200 Google Scholar Crossref Search ADS WorldCat Zietlow D.W. , Sheehan A.F., Molnar P.H., Savage M.K., Hirth G., Collins J.A., Hager B.H., 2014 . Upper mantle seismic anisotropy at a strike-slip boundary: South Island, New Zealand , J. geophys. Res.: Solid Earth , 119 ( 2 ), 1020 – 1040 . 10.1002/2013JB010676 Google Scholar Crossref Search ADS WorldCat Zinke J.C. , 2000 . Structure-related and stress-induced shear-wave velocity anisotropy: observations from microearthquakes near the Calaveras Fault in Central California , Bull. seism. Soc. Am. , 90 ( 5 ), 1305 – 1312 . 10.1785/0119990099 Google Scholar Crossref Search ADS WorldCat Zoback M.D. , Zoback M.L., 2002 . Stress in the earth’s lithosphere , in Encyclopedia of Physical Science and Technology , 3rd edn., Vol. 16 , pp. 1221 – 1232 ., Kluwer Academic Publishers . Google Scholar Google Preview OpenURL Placeholder Text WorldCat COPAC © The Author(s) 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. © The Author(s) 2020. Published by Oxford University Press on behalf of The Royal Astronomical Society.
TI - Spatio-temporal analysis of seismic anisotropy associated with the Cook Strait and Kaikōura earthquake sequences in New Zealand
JF - Geophysical Journal International
DO - 10.1093/gji/ggaa433
DA - 2020-11-20
UR - https://www.deepdyve.com/lp/oxford-university-press/spatio-temporal-analysis-of-seismic-anisotropy-associated-with-the-Ih03ndbcIL
SP - 1987
EP - 2008
VL - 223
IS - 3
DP - DeepDyve
ER -