Systematic variation in anisotropy beneath the mantle wedge in the Java–Sumatra subduction system from shear-wave splitting

Physics of the Earth and Planetary Interiors 178 (2010) 189–201

Contents lists available at ScienceDirect

Physics of the Earth and Planetary Interiors journal homepage: www.elsevier.com/locate/pepi

Systematic variation in anisotropy beneath the mantle wedge in the Java–Sumatra subduction system from shear-wave splitting J.O.S. Hammond a,b,∗ , J. Wookey a,c , S. Kaneshima b , H. Inoue d , T. Yamashina d , P. Harjadi e a

Department of Earth Sciences, University of Bristol, Bristol, UK Department of Earth and Planetary Sciences, Kyushu University, Hakozaki 6-10-1, Fukuoka 812-8581, Japan Department of Earth Sciences, University College London, Gower Street, London, UK d Badan Meteorologi, Klimatologi dan Geofisika, Jakarta, Indonesia e National Research Institute for Earth Science and Disaster Prevention, Tsukuba, Japan b c

a r t i c l e

i n f o

Article history: Received 6 July 2009 Received in revised form 3 October 2009 Accepted 13 October 2009 Edited by: D. Rubie Keywords: Seismic anisotropy Subduction zones Indonesia S-wave splitting

a b s t r a c t The tectonic context of south-east Asia is dominated by subduction. One such major convergent boundary is the Java-Sunda trench, where the Australian–Indian plates are being subducted beneath the Eurasian plate. We measure shear-wave splitting in local and teleseismic data from 12 broadband stations across Sumatra and Java to study the anisotropic characteristics of this subduction system, which can provide important constraints on dynamical processes involved. Splitting in S-waves from local earthquakes between 75 and 300 km deep show roughly trench parallel fast directions, and with time-lags 0.1–1.3 s (92% ≤0.6 s). Splitting from deeper local events and SKS, however, shows larger time-lags (0.8–2.0 s) and significant variation in fast direction. In order to infer patterns of deformation in the slab we apply a hybrid modelling scheme. We raytrace through an isotropic subduction zone velocity model, obtaining event to station raypaths in the upper mantle. We then apply appropriately rotated olivine elastic constants to various parts of the subduction zone, and predict the shear-wave splitting accrued along the raypath. Finally, we perform grid searches for orientation of deformation, and attempt to minimise the misfit between predicted and observed shear-wave splitting. Splitting from the shallow local events is best explained by anisotropy confined to a 40 km over-riding plate with horizontal, trench parallel deformation. However, in order to explain the larger lag times from SKS and deeper events, we must consider an additional region of seismic anisotropy in or around the slab. The slab geometry in the model is constrained by seismicity and regional tomography models, and many SKS raypaths travel large distances within the slab. Models placing anisotropy in the slab produce smaller misfits than those with anisotropy outside for most stations. There is a strong indication that inferred flow directions are different for sub-Sumatran stations than for sub-Javanese, with >60◦ change over ∼375 km. The former appear aligned with the subduction plate motion, whereas the latter are closer to perpendicular, parallel to the trench direction. There are significant differences between the slab being subducted beneath Sumatra, and that beneath Java: age of seafloor, maximum depth of seismicity, relative strength of the bulk sound and shear-wave velocity anomaly and location of volcanic front all vary along the trench. We speculate, therefore, that the anisotropy may be a fossilised signature rather than due to contemporary dynamics. Crown Copyright © 2009 Published by Elsevier B.V. All rights reserved.

1. Introduction At the Java-Sunda trench the Indian–Australian plates are being subducted beneath the Eurasian plate. This is one of the most seismically active regions in the world, with earthquakes occurring from the surface to the base of the transition zone. This makes

∗ Corresponding author at: University of Bristol, Department of Earth Science, Wills Memorial Building, Queens Road, Bristol, BS8 1RJ, UK. Tel.: +44 117 954 5245; fax: +44 117 925 3385. E-mail address: [email protected] (J.O.S. Hammond).

the region an ideal natural laboratory to study subduction dynamics using seismic imaging techniques. Here, we present a study of seismic anisotropy (the variation of seismic wavespeed with direction) beneath the Java-Sunda subduction system, providing some insights into the deformation history and active processes in a subduction setting. As a seismic shear-wave propagates through an anisotropic medium it is split into two quasi shear-waves, one with particle motion in the direction of fast seismic wavespeed, and another perpendicular to this. This results in the two shear-waves travelling at different speeds and arriving at a seismic station separately; a phenomenon called shear-wave splitting (see Crampin and Booth,

0031-9201/$ – see front matter. Crown Copyright © 2009 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.pepi.2009.10.003

190

J.O.S. Hammond et al. / Physics of the Earth and Planetary Interiors 178 (2010) 189–201

1985; Silver and Chan, 1991). At a 3-component seismic station it is possible to measure these separate shear-waves and determine the fast direction of wave polarisation and the time between the shear-waves (see Silver, 1996; Savage, 1999; Kendall, 2000, for reviews). In many cases anisotropy in the upper mantle is assumed to be caused by the lattice preferred orientation (LPO) of olivine. Typically, the olivine fast axis (a-axis) aligns in the direction of upper-mantle flow (e.g., Babuska and Cara, 1991; Mainprice et al., 2000; Mainprice, 2007). As a result shear-wave splitting can provide direct information about dynamic processes, such as mantle flow, and accumulated strain due to previous deformation events which have ‘frozen’ a source of anisotropy into the lithosphere beneath a seismic station. This technique has been applied to many subduction zone settings around the world, both using local slab related events, and teleseismic events (e.g., Fouch and Fischer, 1996; Long and Silver, 2008). The results are far from simple with fast directions appearing trench parallel, perpendicular and oblique, often at the same subduction zone [e.g., New Zealand, (Morley et al., 2006), Japan (Nakajima and Hasegawa, 2004), South America (Helffrich et al., 2002) and the Carribean (Pinero-Felicangeli and Kendall, 2008)]. This suggests a complex region of anisotropy and many different mechanisms have been invoked to explain these results. Deformation in the overriding plate and return flow in the mantle wedge are discussed as causing trench parallel and perpendicular fast directions for New Zealand (Morley et al., 2006). The presence of water has been suggested to change the alignment mechanism of olivine, thus producing fast directions perpendicular to the direction of maximum stress (Jung and Karato, 2001) and this has been invoked to explain trench parallel splitting beneath Japan (Nakajima and Hasegawa, 2004), although the authors also acknowledge that trench parallel flow, caused by along strike dip variations might also explain these results. Evidence for deeper sources of anisotropy has also been described to explain splitting seen in teleseismic shear-wave arrivals, for example slab anisotropy in New Zealand (Brisbourne et al., 1999), sub-slab anisotropy due to toroidal flow beneath the slab in the Carribean (Pinero-Felicangeli and Kendall, 2008) and basal drag beneath the slab in South America (Helffrich et al., 2002). Some attempts have been made to explain subduction zone shear-wave splitting observations with uniform global models. Long and Silver (2008) compiled global observations of shearwave splitting at subduction zones. They postulated that the largely trench parallel splitting observations which arise from sub-mantle wedge regions (i.e. slab or sub-slab) are dominated by lateral flow in the sub-slab region. This is induced by return flow due to trench migration. Long and Silver (2008) explain mantle wedge anisotropy (above the slab) as being due to the interaction of lateral flow due to trench migration and corner flow driven by viscous coupling between the slab and the wedge. Such corner flow would induce a trench perpendicular fast direction in the splitting. The authors argue that the observed supra-slab trench parallel or trench perpendicular splitting results observed globally are due to the dominance of one or the other of these mechanisms. Another model invoked to explain global observations of subduction zone anisotropy suggests preferentially-oriented hydrated faults in the subducting slab (Faccenda et al., 2008). During the subduction process, faulting in the downgoing plate allows water to migrate into the crust and mantle. This results in the formation of sheets of silicates in the top few kms of the slab, inducing a shape preferred orientation, resulting in high degrees of anisotropy. Superimposed on this, an alignment of serpentinized rocks will cause an additional crystallographic preferred orientation, thus giving the potential for large degrees of shear-wave splitting from a thin layer at the top of the slab (Faccenda et al., 2008).

The variability in shear-wave splitting results, and the wide variety of mechanisms described to explain them shows that subduction zone settings are very complicated. Possible regions of anisotropy exist in the crust, lithosphere and mantle wedge, and in or below the slab. In this study we try to understand the problem in more detail by forward modelling our shear-wave splitting results with realistic subduction zone models. 2. Data and methodology We use data from 9 stations of the JISNET (Japan Indonesia Seismic Network) array (Ohtaki et al., 2000) to determine the anisotropic characteristics beneath Sumatra and Java. The stations were deployed in 1998 and some stations are still recording. We use data up until the end of 2003. We also use 2 stations from the GEOFON array (UGM, 9 years, GSI, 3 years) and 1 from the Ocean hemisphere network project (PSI, 13 years), a station used by Long and Silver (2008) to investigate subduction zone shear-wave splitting. Station locations are shown in Fig. 1. We employ earthquakes with distances between 85◦ and 140◦ , events ideal for SKS/SKKS-wave splitting analysis. A total of 47 earthquake-station paths provided SKS/SKKS-phases of sufficient quality to enable splitting to be estimated (2% of data used). All SKS/SKKS data were filtered between 0.05 and 0.3 Hz. We also measure splitting in S-waves from local earthquakes where the incidence angle is less than 45◦ [i.e., within the shear-wave window (Evans, 1984)]. Earthquake locations from the catalogue of Engdahl et al. (1998) are used where possible, others are taken from the USGS catalogue, which could introduce some significant mislocation errors into our results. This results in 127 measurements of splitting from local earthquakes (49% of data used). Individual results are shown in Tables S1–S4. All local earthquake data were filtered between 0.1 and 1 Hz. The overlap between SKS/SKKS- and local S-wave filter bands helps to minimise frequency dependent effects between the two data sets. We estimate shear-wave splitting in both the core phases and S-phases from local earthquakes using the method of Teanby et al. (2004), which is based on the methodology of Silver and Chan (1991). It uses a grid search over fast direction () and time-lag (ıt), minimising the second eigenvalue of the covariance matrix for the particle motion around a given time window. The advantage of the Teanby et al. (2004) method is that it automates window selection, making 100 measurements around the relevant phase, and applying cluster analysis to identify the most stable result. Fig. 2 displays a typical result of a local S- and teleseismic SKS-wave splitting result at PSI. 3. Results From our initial analysis we reject results with 1 errors greater than 0.4 s in ıt and 30◦ in . Furthermore, the polarisation of the shear-wave after the anisotropy correction is applied (source polarisation, hereafter spol) is compared with the back azimuth for SKS arrivals. For radially polarised phases (i.e. a P–S conversion at the core/mantle boundary) these should be similar, thus we reject results with differences >30◦ . Local earthquakes are not radially polarised as they have not been converted to a P-wave at any point along their ray path. However some earthquakes are large enough to have a previously determined focal mechanism (we use the Global CMT and USGS catalogue), from which the initial polarisation can be calculated. For these earthquakes the spol measurement and initial polarisation are compared to provide an extra quality assessment (10% of local splitting measurements included this check). Again spol and the initial polarisation must agree within 30◦ to be included. Those results coming from earth-

J.O.S. Hammond et al. / Physics of the Earth and Planetary Interiors 178 (2010) 189–201

191

Fig. 1. Map showing stations and events used in this study. White inverted triangles mark stations, blue circles show earthquakes (local events on main map, teleseismic events on top right map). Also shown are quaternary volcanos (red triangle), absolute plate motions (Gripp and Gordon (1990), black arrows), slab contours at 100 km depth intervals (Gudmundsson and Sambridge, 1998), the Great Sumatra Fault, and Isochrons (Müller et al., 1997). (For interpretation of references to colour, the reader is referred to the web version of this article.)

quakes without an available focal mechanism are assumed to be correct. Results are shown in Fig. 3 and Tables S1–S4. We discuss the results for Sumatra and Java separately in the following sections.

deeper source of anisotropy, within or beneath the slab, must be present.

3.1. Sumatra

3.2.1. Local S-wave splitting Local earthquakes, within the shear-wave window, have been recorded down to depths of 700 km at most stations in Java. Splitting results from these events display considerable variation in both fast direction and time lag on individual stations as well as from different stations (Fig. 3). Most stations show both trench parallel and trench perpendicular fast directions. The trench perpendicular results are in general from intermediate depth earthquakes (250–350 km), with both shallower and deeper earthquakes having more trench parallel fast directions. Some deep events produce large delay times and fast directions oblique to the strike of the trench (UGM, BRB). These seem isolated to one region of the subduction zone, with adjacent splitting measurements from nearby earthquakes showing much smaller delay times, albeit with similar fast directions. Splitting at BMI is very complicated, possibly because it lies above a part of the slab which bends beneath Sulawesi.

3.1.1. Local S-wave splitting Seismic stations BSI, PSI, PPI, KSI and KOTA are located approximately above the 150 km contour of the subducting Indian plate (based on the slab contours of Gudmundsson and Sambridge, 1998). As a result splitting is estimated in shear-waves generated from earthquakes in the ∼100–200 km depth range. One station, GSI, is located closer to the trench, and thus provides estimates of splitting from earthquakes in the depth range ∼30–80 km. All these splitting results show little variation with depth (Fig. 4). Time-lags for all stations, including GSI, are generally between 0.2 and 0.4 s, and fast directions roughly trench, parallel, with some scatter, and parallel to the strike of the Sumatran fault (Figs. 3 and 4). This indicates that the source of this anisotropy is in the over-riding plate and that the mantle wedge above the slab is largely isotropic. 3.1.2. SKS/SKKS-wave splitting The stations located above the 150 km slab contour (BSI, PSI, PPI, KSI and KOTA) show ∼1 s of splitting, with a general northsouth trend. However some variation in the fast direction is evident, particularly at PSI (Fig. 3). The station close to the trench, GSI, is markedly different with a large time-lag of 1.8 s and a trench perpendicular fast direction of 52.5◦ . The obvious discrepancy between the time-lags for local and teleseismic splitting suggests that a

3.2. Java

3.2.2. SKS/SKKS-wave splitting In general SKS-wave results have an east-west fast direction, with ∼1.6 s of splitting (Fig. 3). There does seem to be some subtle variation in the splitting results with the estimated fast directions occurring in two groups, ∼75◦ and ∼105◦ . The magnitude of splitting is again larger than the splitting seen in local events, except for a few deep earthquakes at UGM, and BRB. Again, similarly to

192

J.O.S. Hammond et al. / Physics of the Earth and Planetary Interiors 178 (2010) 189–201

Fig. 2. Example of the Teanby et al. (2004) method of shear-wave splitting. (Top) SKS-wave splitting result at PSI and, bottom) local S-wave splitting result for PSI. (a and g) Original traces (E, N, and Z). (b and h) Traces rotated into R and T directions before and after the anisotropy correction. R component is the initial shear-wave polarization before entering the anisotropic region. T component is perpendicular to the R component. Energy on the corrected transverse trace should be minimised in the analysis window. (c and i) The top traces show the fast/slow shear waveforms for uncorrected (left) and corrected (middle(normalised)/right(real amplitudes)) seismograms. The bottom panels show the particle motion for uncorrected (left) and corrected (right) seismograms. A good result will show similar fast/slow shear waveforms and any elliptical

J.O.S. Hammond et al. / Physics of the Earth and Planetary Interiors 178 (2010) 189–201

193

Fig. 4. Variation in local splitting delay time as a function of depth. Coloured crosses show local S-wave splitting delay times as a function of depth. SKS-wave splitting delay times are shown by coloured squares in the grey box at the bottom of the figure. ˙ deep, and the difference Note the lack of variation with depth for events <300km between delay time for events <300 km deep and events >300 km deep and SKSevents. (For interpretation of references to colour, the reader is referred to the web version of this article.)

Fig. 3. Local and SKS-wave splitting results beneath Sumatra (top), and Java (bottom). Yellow and blue bars show splitting result for SKS-wave arrivals and local S-waves respectively. Bars are oriented in the direction of the fast shear-wave polarisation, and the length of the bar is proportional to the amount of splitting. SKS-wave splitting results are plotted at the station location, local S-wave splitting results are plotted at the event location. Other symbols are as in Fig. 1. (For interpretation of references to colour, the reader is referred to the web version of this article.)

the results at Sumatra, this suggests a slab or sub-slab source of anisotropy. 4. Modelling It is evident that more than one region of anisotropy exists beneath Sumatra and Java. Shallow earthquakes (<300 km deep) beneath Sumatra and Java show little variation in time-lag with depth, and show consistent trench parallel fast directions. This suggests that the mantle wedge is isotropic, and that the splitting observed from these earthquakes is accrued in the over-riding plate. However larger time-lags observed in SKS/SKKS splitting at Sumatra and Java, and splitting from deeper (>300 km deep) earth-

quakes beneath Java suggest that anisotropy must exist in the slab, or sub-slab as well. To understand the anisotropy in more detail we have developed a technique based on ray-tracing through an inhomogeneous subduction zone model (Guest and Kendall, 1993). We build a model consisting of four regions: layer 1, the over-riding plate and mantle wedge; layer 2, a thin 10 km layer above the slab; layer 3, the subducting slab; and layer 4, the sub-slab mantle (Fig. 5). The surface of the slab is calculated using the regionalized upper mantle (RUM) seismic model slab contours (Gudmundsson and Sambridge, 1998) which is consistent with a relocated local seismicity study (Fauzi et al., 1996) and deep events from the USGS earthquake catalogue. The model is also consistent with regional tomography models for Indonesia (Gorbatov and Kennett, 2003). Uncertainty in these slab models could introduce error into these results, however recently published slab models of Syracuse and Abers (2006) show only minor changes (∼20 km shift east) in Sumatra–Java subduction zone to those calculated by Gudmundsson and Sambridge (1998). We use ATRAK, a ray-tracer capable of tracking seismic rays in multilayered 3D media with general anisotropy (Guest and Kendall, 1993). This technique has been utilised to investigate anisotropy in other subduction zones [e.g., Kurils (Kendall and Thomson, 1992) and New Zealand (Brisbourne et al., 1999)]. Our model differs from previous studies by the fact that we trace rays through an isotropic model, assuming that anisotropy only minimally distorts the raypath (a reasonable assumption for relatively weakly anisotropic media). We trace rays for each station separately (except UGM and SWH which are very close together, and so are incorporated into one model). We appeal to reciprocity and shoot rays from the station, saturating the model (>12,000 rays per station, Fig. 5). We then match each event to its nearest ray, thus defining the raypath

particle motion will have been linearised. (d and j) The results of the grid search over ıt and . The method used minimises the second eigenvalue of the particle motion (i.e. the best result occurs where the particle motion is linear after removing the splitting). The optimum splitting parameters are represented by the cross, and the 1st surrounding contour denotes the 95% confidence region. (e and k) Measurements of ıt and  obtained from 100 different analysis windows plotted against window number. A stable result will have a similar solution independent of analysis window. (f and l) Cluster analysis of ıt and  obtained from 100 different analysis windows. A good result will have solutions clustered in one area.

194

J.O.S. Hammond et al. / Physics of the Earth and Planetary Interiors 178 (2010) 189–201

Fig. 5. (a) Ray tracing through an isotropic slab model to find the path travelled from local events to the surface. Events are matched to the nearest point on a ray (local) or to their ak135 predicted piercing point at a depth of 800 km (SKS/SKKS). The colour scale relates to the velocity structre in the model. (b) The raypaths for stations UGM/SWH, showing the proportion of the ray spent in each part of the model (represented by the various colors shown in the key). Splitting is predicted using the Christoffel equation ˇ ` 2001) for an appropriately rotated olivine tensor. (For interpretation of references to colour, the reader is referred to the web version of this article.) (Cerven y,

taken from the earthquake to the station. For local earthquakes we use the latitude, longitude and depth to match the ray, and for SKS/SKKS phases we calculate the piercing point of the phase at a depth of 800 km beneath the station [through the AK135 1-D velocity model (Kennett et al., 1995)], and use this information as the

event location (Fig. 5). This ray-path for each station-event pair is then broken down into a length, time and angle spent in each part of the model. With this information we can then impose an anisotropy on a specified region of the model (layers 1–4), represented by an appropriately rotated elastic tensor (we use olivine elastic con-

Fig. 6. An example of the shallow and slab anisotropy modelling result for PPI. For the shallow model anisotropy is confined to the uppermost 40 km of layer 1, representing the over-riding plate, and events from depths <300 km only are modelled. For the slab model the modelled anisotropy from the shallow anisotropy is fixed in layer 1, and a second layer of anisotropy exists in the slab. Events from depths >300 km and SKS/SKKS-phases are modelled. (a and d) The calculated misfit surfaces for the grid search of anisotropy strength (i.e. degree of alignment—VFA) and a-axis orientation () in a horizontal plane. Nodes in the grid are evaluated as the summed misfit across all event-station pairs, calculated for each ıt and  separately and also an equally weighted combined misfit. (b and e) The comparison of modelled splitting (blue) and observed (black) for the optimum model, plotted at event location. (c and f) A density rotorgram of the misfit in a-axis direction for shallow and slab anisotropy respectively (; black is lowest misfit, white is highest), at the best fitting anisotropy strength VFA. This shows graphically the confidence in orientation of our best fitting model. Note the good fit to the data, and the well constrained direction. (For interpretation of references to colour, the reader is referred to the web version of this article.)

J.O.S. Hammond et al. / Physics of the Earth and Planetary Interiors 178 (2010) 189–201

195

196

J.O.S. Hammond et al. / Physics of the Earth and Planetary Interiors 178 (2010) 189–201

stants from Abramson et al. (1997)). Anisotropy is confined to the olivine stability field (i.e. <410 km), and can be restricted to a particular depth range within a layer. For each contiguous section of ˇ ` 2001) the raypath we can use the Christoffel equation (Cerven y, to calculate the fast and slow shear-wave velocities and particle motion, thus, along with the distance in the layer, we can calculate a time-lag and fast direction through the layer. We then perform a grid search over orientation (,  = 5◦ ), and anisotropic strength (by Voigt–Reuss–Hill dilution of the elastic tensor; denoted volume fraction aligned—VFA, VFA = 0.01–0.05). The angle  is the orientation of the olivine a-axis in a plane which is either horizontal (shallow models) or inclined at the average dip of the subducting slab (deep models). The b-axis is held normal to this plane. For each node in the grid we calculate ıt and  for each event-station raypath in the data. The model is then evaluated as the summed misfit between these and the real measurements (for ıt and  separately, and for a combination). Typically, VFA is poorly constrained, so is not interpreted in great detail, instead we primarily analyse the best fitting orientation , illustrated by the density rotorgram in Fig. 6. This can be performed for each layer in the model. Four models are run in total, the first constraining the layer 1 anisotropy, and the last three incorporate this result with layers 2, 3 and 4 respectively. Multiple layers of anisotropy can be included, though only one may vary, as a result we can not formerly asses trade-off between these layers. Where two layers of anisotropy are included we combine the splitting operators using the n-layer splitting equations of Silver and Savage (1994), and the initial source polarisation from the splitting measurements. An example result is shown in Fig. 6 for station PPI with anisotropy in layer 1 (Fig. 6a–c) and anisotropy in both layers 1 and 3 (Fig. 6d–f). The locations of stations and earthquakes used in this study are ideal to image the various parts of the subduction system. We image paths in the sub-slab, slab, supra-slab, mantle wedge and overriding plate. This enables us to go one step further than typical shear-wave splitting studies and place constraints on the location of the anisotropy, and thus better understand the dynamics of JavaSunda subduction. The modelling scheme outlined here allows us to extract information on which part of the subduction system each event is sampling. Therefore, we are able to investigate further the observation made in Section 3, that the mantle wedge appears largely isotropic. Fig. 7 shows the variation in delay time as a function of path length through various parts of the model. It is evident that there is no clear relationship between splitting delay times and path length in the mantle wedge. In fact, those data with paths almost exclusively in the wedge display the smallest delay times (<0.4 s, except two SKS events at PSI which, in our model, never sample the region below the wedge). An estimate can be placed on the amount of anisotropy in a 40 km thick over-riding plate and the mantle wedge below. Fig. 7a shows a linear fit to the data with paths exclusively in the wedge (removing the two SKS-wave arrivals at PSI). The line has equation ıt = 0.0007l + 0.1306, where l is the path length in the mantle wedge and over-riding plate. Assuming an average S-wave velocity of 4.5 km/s, this equates to an anisotropy of 1.8% in a 40 km thick plate, and 0.3% anisotropy in the mantle wedge below this plate. This shows that the anisotropy in the over-riding plate is dominating the shear-wave splitting. In comparison, some correlation is evident between splitting delay time and path length spent in the region beneath the mantle wedge, but above 410 km (layers 2–4 in our model) (Fig. 7b). The splitting results can be explained by a model which has 20% aligned olivine crystals along this path length. This observation supports the idea that the mantle wedge is largely isotropic, and that anisotropy is largely confined to the over-riding plate, slab and sub-slab region. To further constrain the

Fig. 7. Variation in delay time as a function of (a) path length in the mantle wedge and over-riding plate (layer 1 in Fig. 5) and b) path length above 410 km but beneath the mantle wedge (layers 2–4 in Fig. 5) for all local and SKS splitting results. Red crosses in (a) show splitting estimates for events which spend most of their path in the mantle wedge and over-riding plate (i.e. less than 10 km outside of this region). The red line shows the best fitting slope to the red crosses estimating an anisotropy of 1.8% in a 40 km thick over-riding plate, and 0.3% anisotropy in the mantle wedge below. The grey shaded region shows the range of potential splitting for an olivine fabric with 20% alignment (based on the minimum and maximum shear-wave anisotropy of the olivine elastic tensor). Note the poor correlation between the shaded area and the splitting results in (a), and the consistently low delay times for events which spend most of their time in the mantle wedge and over-riding plate. The splitting results are explained much better by putting anisotropy in the slab or subslab region (panel b). This shows that the mantle wedge is largely isotropic. (For interpretation of references to colour, the reader is referred to the web version of this article.)

location of the anisotropy we split the modelling into two steps, shallow and deep anisotropy. To constrain the characteristics of the shallow anisotropy we use events with depths <300 km, and to model deeper anisotropy we use events with depths >300 km and SKS-events.

J.O.S. Hammond et al. / Physics of the Earth and Planetary Interiors 178 (2010) 189–201

197

Table 1 Best fitting models for shallow earthquakes (<300 km). Anisotropy is confined to a 40 km thick layer at the surface, and fits the data well with trench parallel anisotropy orientations. Station

Number of events (n)

Over-riding plate (layer 1) 

VFA

2

2 /n

BSI GSI PSI PPI KSI KOTA Sumatra misfits

7 8 35 5 8 4 67

-20 -15 -55 -20 -55 -20

0.20 0.30 0.40 0.20 0.25 0.50

46.92 38.38 1486.22 26.98 126.90 268.79 1994.19

6.70 4.80 42.46 5.40 15.86 67.20 29.76

LEM UGM + SWH BMI KHK Java misfits

6 15 18 2 41

-90 -35 -80 -75

0.45 0.25 0.30 0.25

147.14 1426.24 1107.10 95.14 2775.62

24.52 95.08 61.50 47.57 67.70

4769.81

44.16

Overall misfits

108

4.1. Shallow anisotropy We impose anisotropy in the uppermost 40 km of layer 1 of the model, simulating the overriding plate lithosphere. The choice of this thickness is fairly arbitrary (i.e. a layer half as thick with twice the anisotropy strength will fit the data just as well), but we base our model on the limit of our splitting results (the shallowest events at GSI are ∼30 km). In these cases the olivine a-axis is horizontal and the misfit is calculated only using events shallower than 300 km. Table 1 and Fig. 8 show that the best fitting anisotropy orientations at all stations are well constrained and show trench parallel orientations. BRB is not modelled as we have too few splitting measurements. Errors in depth locations could have an effect on our estimated VFA, however we expect this effect to be minimal due to the fact that we impose a 40 km thick layer of anisotropy at the top of the model. Events which are mislocated from below 40 km to above 40 km, and events which are mislocated to below 40 km from above 40 km will cause the VFA to be over and underestimated respectively. This is only a problem at GSI, where event depths range from 30 to 80 km. At other stations only one event occurs shallower than 70 km (SWH, 33 km). 4.2. Deep anisotropy To fit local slab events with depths >300 km, and SKS/SKKS phases we test 3 classes of model where anisotropy is confined to a thin layer above the slab (layer 2), the slab (layer 3) or the sub-slab (layer 4). We impose an olivine anisotropy with its a-axis oriented () in a plane inclined at the average dip of the subduction zone (baxis is normal to the plane). In all deep models we also impose the best fitting shallow anisotropy for each station (determined above). Thus, each model has two anisotropic regions, one fixed in orientation and strength (shallow; layer 1) and one varying (deep; layers 2–4). Results are shown in Figs. 8 and 9 and Table 2, and it is evident that it is hard to confidently discriminate between the slab and sub-slab models (layers 3 and 4), with both giving compatible anisotropy orientations. Anisotropy in a thin layer at the top of the slab (layer 2) consistently has the largest misfit, and requires a very high degree of alignment, suggesting this is a less likely cause of our observations (Fig. 9 and Table 2). For both Java and Sumatra the misfit is considerably lower for a model where anisotropy is confined to the slab (assuming a thickness of 100 km). A robust feature in all models (with the exception of GSI) is that the deeper anisotropy orientation beneath Sumatra is more north-south, and rapidly changes at the Sunda strait becoming more east-west

Fig. 8. Density rotorgrams (see Fig. 6) for models of subduction zone anisotropy showing the best fitting orientations obtained from forward modelling. (a) Olivine orientations best fitting the splitting results (red bars) obtained from local events <300 km deep; these are predominantly trench parallel. (b) Olivine orientations best fitting splitting results from local events >300 km deep and SKS/SKKS-wave splitting results (red bars). Note the rotation from north-south orientations beneath Sumatra (GSI is an exception), and the east west anisotropy orientations beneath Java. The complicated pattern at BMI reflects the scatter observed in the data. We speculate that this is due to complications arising from the slab bending beneath this station. (For interpretation of references to colour, the reader is referred to the web version of this article.)

beneath Java (Figs. 8 and 9). GSI shows a trench perpendicular anisotropy orientation, suggesting that a different mechanism is causing anisotropy at the trench compared to that beneath the volcanic/back-arc region (Fig. 8). GSI has the longest path in the sub-slab region. The anomalous trench perpendicular anisotropy orientation, and large time delays at this station may suggest that sub-slab anisotropy is dominant in the forearc. BRB and KHK do not have enough observations to model this. Mislocated events for deep and SKS/SKKS-wave events will have minimal effects on their raypaths. 5. Discussion From the modelling it is evident that anisotropy in the overriding plate is required to explain the splitting observations from events <300 km deep at Java and Sumatra. The modelled anisotropy orientation is near parallel to the strike of the trench. The simplest explanation for this anisotropy is vertically aligned cracks associated with deformation in the over-riding Eurasian plate, with a largely isotropic mantle wedge below. This observation is supported by a recent local earthquake tomography study for central Java, where high P-wave anisotropy (7–10%) is observed in the crust

198

J.O.S. Hammond et al. / Physics of the Earth and Planetary Interiors 178 (2010) 189–201

Fig. 9. A map showing the misfits (bars) for models with anisotropy in a thin layer above the slab (red), in the slab (green), and below the slab (blue), along with the best fitting anisotropy orientations (arrows). Note, the best fitting model places anisotropy in the slab with north-south orientations beneath Sumatra, and east-west orientations beneath Java. (For interpretation of references to colour, the reader is referred to the web version of this article.)

and lithospheric mantle (Koulakov et al., 2009). Similar observations of over-riding plate anisotropy and isotropic mantle wedges have been seen at New Zealand (Morley et al., 2006), the Carribean (Pinero-Felicangeli and Kendall, 2008) and South America (Polet et al., 2000). The fact that the wedge is isotropic suggests that either 2D corner flow caused by viscous coupling of the slab and wedge is of similar magnitude as lateral flow from trench migration, or that they are both weak in this region, thus meaning no coherent LPO can be developed (Long and Silver, 2008). This observation is in drastic contrast to other subduction zones such as Japan, where Hiramatsu et al. (1998) observe two 75–175 km deep anisotropic zones that are caused by melt filled cracks, Tonga where Smith et al. (2001) suggest that along strike flow in the mantle related to slab rollback and the Samoan plume is inducing an anisotropy and the Aleutians where Long and Silver (2008) suggest that 2D corner flow is causing a large anisotropy. In all these regions they observe large ıt ∼1–2 s, much larger than we observe beneath Java and Sumatra. Models of deeper anisotropy, constrained by deep subduction events, and SKS/SKKS show that a slab source of anisotropy is most

likely, although a sub-slab source can not be ruled out and gives similar best fitting anisotropy orientations. However, in all these models a sharp change in anisotropic characteristics is observed between Sumatra and Java (Fig. 8). This change is also evident in the raw data (Fig. 10). SKS fast direction estimates for the same event recorded at KOTA, the southernmost Sumatran station, and LEM the westernmost Javan station vary by ∼60◦ over ∼375 km (Fig. 10). This change in anisotropic characteristics is mirrored by changes in other geophysical parameters close to the Sunda Strait. The most obvious change near this location is in the age of the downgoing slab (Fig. 1). Magnetic anomalies show a change in age from ∼ 60 to ∼100 Ma from southern Sumatra to Western Java (Sdrolias and Muller, 2006). The sharpness, and exact location of this boundary is unclear, especially at depth in the subducted slab. Syracuse and Abers (2006) observed that this age boundary co-incides with a change in volcanic characteristics. They found that the volcanic front in Sumatra is located above the 90 km slab contour, with an abrupt change close to the Sunda Strait where volcanism is present

Table 2 Best fitting models for deep earthquakes (>300 km) and SKS phases. A model with anisotropy in the slab fits splitting results best. In all models the best result from the shallow models (Fig. 8 and Table 1) are included using the n-layer splitting equations of Silver and Savage (1994). Station

Number of events (n)

Supra-slab (layer 2) 

Slab (layer 3)

VFA

2

2 /n



VFA

Sub-slab (layer 4) 2

2 /n



VFA

2

2 /n

BSI GSI PSI PPI KSI KOTA Sumatra misfits

11 2 11 6 4 2 36

0 30 15 85 0 5

1.00 0.65 0.55 1.00 1.00 0.45

43.79 248.01 531.06 269.93 57.74 315.19 1465.72

3.98 124.00 48.29 44.99 14.44 157.60 40.71

0 55 -30 10 5 10

0.30 0.80 0.10 0.95 0.70 0.60

44.88 3.92 568.75 98.37 21.55 33.25 770.72

4.08 1.96 51.70 16.40 5.39 16.63 21.41

0 55 -35 0 5 0

0.1 0.3 0.10 0.85 0.45 0.40

48.01 5.20 575.47 195.40 10.91 221.97 1056.96

4.36 2.60 52.32 32.57 2.72 110.99 29.36

LEM UGM + SWH BMI KHK Java misfits

6 7 10

-90 -35 -15

1.00 0.95 0.85

295.74 147.33 296.16

49.29 21.05 29.62

-65 -75 55

0.45 0.5 0.55

81.89 114.37 152.79

13.65 16.34 15.28

75 -85 -75

0.5 0.25 0.25

125.68 244.74 146.57

20.95 34.96 14.66

23

739.23

32.14

349.05

15.18

516.99

22.48

Overall misfits

59

2204.95

37.37

1119.77

18.98

1573.95

26.68

J.O.S. Hammond et al. / Physics of the Earth and Planetary Interiors 178 (2010) 189–201

199

Fig. 10. SKS-wave splitting results from one event (1999/01/28 08:10:05, 52.89◦ N, 169.12◦ W) recorded at KOTA and LEM. Note the completely different splitting results even though the stations are ∼375 km apart.

above the 150 km slab contour. This change occurs over <150 km, similarly to the abrupt change in anisotropy, and also correlates with the abrupt reduction in deep (>500 km) seismicity beneath Sumatra (e.g., Engdahl et al., 1998), although no change in the shape of the slab is observed. Syracuse and Abers (2006) also correlate this change in the location of the volcanic front with variation in the ratio of potassium to silica (Wheller et al., 1987), but the reasons for this are unclear (Syracuse and Abers, 2006). Regional tomography profiles of Gorbatov and Kennett (2003) show an abrupt change in the bulk sound speed in the slab close to the Sunda Strait further suggesting a sharp change in the slab characteristics in this region. The possible causes of the observed rapid change in anisotropy, with reference to the other changes in geophysical characteristics are discussed further below.

5.2. Water in the mantle

5.1. Change in geometry

Fig. 5 shows that the deep and teleseismic events travel some distance in the subducting slab, and thus anisotropy in the slab will manifest itself in the observed splitting results. Interestingly, the sharp change in anisotropic characteristics between Java and Sumatra co-incides with major changes in slab characteristics. The age of formation of the subducting plate beneath Java and Sumatra varies between 100 Ma in the north to 43 Ma at 2◦ S and up to 160 Ma in the east (Sdrolias and Muller, 2006). This transition in ages occurs near the region where we see an abrupt change in anisotropy. The plate subducting beneath Sumatra formed 40–80 Ma and the plate subducting beneath Java formed 80–140 Ma. Other features of the slab are also different between these regions. Seismicity beneath Sumatra only extends to a depth of ∼300 km and as the slab bends beneath Java it extends to the base of the transition zone and this change in seismicity is located approximately at the same location as an abrupt change in the bulk sound speed (Gorbatov and Kennett, 2003) and change in the location of the volcanic front (Syracuse and Abers, 2006). The abrupt changes in age and seismic properties suggests that differences between the Javan and Sumatran slabs may be the cause for the abrupt change in anisotropic characteristics. Faccenda et al. (2008) show that preferentially oriented hydrated faults in the uppermost part of the slab can produce a high degree of shear-wave splitting. A change in age between the slab beneath western Java and southern Sumatra may effect the efficacy of fracture generation in the slab, and thus hydrated faults may change in style or quantity across this region. However, our models indicate that a thin highly anisotropic layer at the top of the slab poorly fits our results, mainly due to the fact that rays spend such a small amount of time in this region. If this layer behaves as a

The Indian and Australian plates are subducting northwards beneath the Java-Sunda arc. As a result subduction moves from normal to the trench beneath Java to oblique subduction beneath Sumatra (most obviously characterised by the Great Sumatran transform fault). This change in the geometry of subduction across the region could potentially explain the change in anisotropy observed in the splitting results. The alignment of the deep anisotropy orientations at all stations beneath Sumatra, except GSI, roughly matches with the APM of the subducting slab and oblique to the orientation of the trench. This would suggest that in this case shear associated with the subducting plate moving through the mantle is aligning sub-slab olivine, and causing APM aligned anisotropy. However, beneath Java, the opposite is true. The deep anisotropy orientation is perpendicular to the APM of the subducting slab. The slab between Sumatra and Java appears continuous (Syracuse and Abers, 2006), so it seems unlikely that basal drag would create an anisotropy beneath southern Sumatra, and not western Java <375 km away. Thus, basal drag does not seem a likely cause of anisotropy. Long and Silver (2008) suggest that lateral flow beneath the slab, produced by trench migration, can induce a trench parallel anisotropy orientation, and this is observed beneath Java. However, this would predict trench parallel anisotropy orientations beneath Sumatra, and we observe trench oblique orientations. Again, this change from lateral flow beneath western Java, to an oblique flow beneath southern Sumatra seems unlikely over such a short distance.

Jung and Karato (2001) have suggested that the presence of water can make the b-axis of olivine align with the flow direction, manifesting in a anisotropy orientation perpendicular to the direction of flow. Thus a change in the amount of water beneath Java and Sumatra could induce a change from APM parallel anisotropy orientations to APM normal anisotropy orientations. The most likely region of water rich mantle is above the slab, in the mantle wedge. We observe that the mantle wedge in this region is largely isotropic, with little variation in time-lags with depth, suggesting that this is an unlikely scenario to cause the sharp change in anisotropy. 5.3. Fossil anisotropy

200

J.O.S. Hammond et al. / Physics of the Earth and Planetary Interiors 178 (2010) 189–201

strong waveguide, a feature not modelled by our simple subduction models, then this may explain the observed shear-wave splitting. Another possible location for a fossil anisotropy is in the slab lithosphere. As oceanic lithosphere is formed and moves over the asthenosphere it creates a large anisotropic signature (e.g. at the East Pacific Rise Wolfe and Solomon, 1998; Harmon et al., 2004), in the direction of plate motion. These anisotropic characteristics can be frozen into the lithosphere, and may remain present in the subducted slab beneath Java and Sumatra. This mechanism has been suggested to be responsible for at least part of the observed shear-wave splitting seen at the Marianas (Fouch and Fischer, 1998). Also, Brisbourne et al. (1999) show that significant anisotropy, with trench parallel anisotropy orientations, exists in the Hikurangi subduction zone, New Zealand. The Pacific plate in this region has a trench parallel age progression [i.e. the age of the Pacific oceanic lithosphere changes in a direction parallel to the trench (Sdrolias and Muller, 2006)]. This is similar to that seen beneath Java, suggesting that a fossilised anisotropic signature, frozen in as the plate formed, could explain anisotropy in this region. The oceanic lithosphere subducting beneath Java formed 80–140 Ma, and Lithgow-Bertelloni and Richards (1998) show that plate motions for the Australian plate between 74 and 119 Ma are generally east-west. The oceanic lithosphere subducting beneath Sumatra formed 40–80 Ma, and Lithgow-Bertelloni and Richards (1998) show that the plate motions for the Indian plate at this time were north-south. The shear-wave splitting change observed between Java and Sumatra could, therefore, be an effect of sampling these differing frozen anisotropic signatures.

6. Conclusion Shear-wave splitting has been measured in SKS-wave arrivals and S-waves from local subduction zone earthquakes beneath the Java-Sunda subduction zone. Results show that the mantle wedge is mainly isotropic, with anisotropy largely confined to the overriding plate. This has been further emphasised with modelling, which shows a best-fitting trench-parallel anisotropy orientation, suggesting anisotropy seen in splitting from local earthquakes is derived from anisotropy caused by vertically aligned cracks associated with deformation in the over-riding Eurasian plate. Splitting from deeper earthquakes, and SKS-wave arrivals have much larger delay times, suggesting that a deeper source of anisotropy is present. Modelling constrains the location of this anisotropy, with the most likely anisotropic region being the slab. The best fitting anisotropy orientations in the slab were predominantly northsouth beneath Sumatra (except GSI which lies above the trench, and shows trench perpendicular anisotropy orientations, suggesting a different mechanism to that seen at the volcanic arc/back arc). This change in orientation occurs over less than 375 km and at the same location as changes in other geophysical characteristics such as slab age (Sdrolias and Muller, 2006), bulk sound speed (Gorbatov and Kennett, 2003) and location of the volcanic front (Syracuse and Abers, 2006), suggesting that a change in the slab characteristics is the reason for the change in anisotropy. We speculate that the anisotropy has been frozen into the subducting lithosphere as it was formed, or moved over the asthenosphere. The subducting lithosphere beneath Java is older, and would have formed with eastward plate motions (LithgowBertelloni and Richards, 1998), resulting in east-west anisotropy orientations, compared with the younger lithosphere subducting beneath Sumatra which formed when the Indian–Australian plate was moving north (Lithgow-Bertelloni and Richards, 1998), resulting in north-south anisotropy orientations.

Acknowledgments We would like to thank Mike Kendall and George Helffrich for valuable discussion which along with thoughtful discussion and reviews by Frederik Tilmann and one anonymous reviewer helped improve the paper. Also, we would like to thank Alexei Gorbatov for the use of his tomographic models. All the data used in this paper was provided by the Japan Indonesian Seismic Network (JISNET), IRIS, GEOFON and the Ocean Hemisphere Network Project, courtesy of the Earthquake Research Institute, University of Tokyo. This research was funded by a Japan Society for the Promotion of Science (JSPS) postdoctoral fellowship (short term), JSPS/FF1/367. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.pepi.2009.10.003. References Abramson, E., Brown, J., Slutsky, L., Zaug, J., 1997. The elastic constants of San Carlos olivine to 17 GPa. J. Geophys. Res. 102 (6), 12253–12263. Babuska, V., Cara, M., 1991. Seismic Anisotropy in the Earth. Kluwer Academy, Norwell, MA. Brisbourne, A., Stuart, G., Kendall, J.-M., 1999. Anisotropic structure of the Hikurangi subduction zone, New Zealand-integrated interpretation of surface-wave and body-wave observations. Geophys. J. Int. 137 (1), 214–230. ˇ ` V., 2001. Seismic Ray Theory. Cambridge University Press. Cerven y, Crampin, S., Booth, D.C., 1985. Shear-wave polarizations near the North Anatolian Fault, II, Interpretation in terms of crack-induced anisotropy. Geophys. J. R. Astr. Soc. 83, 75–92. Engdahl, E., van der Hilst, R., Buland, R., 1998. Global teleseismic earthquake relocation with improved travel times and procedures for depth determination. Bull. Seism. Soc. Am. 88 (3), 722–743. Evans, R., 1984. Effects of the free surface on shear wavetrains. Geophys. J. R. Astr. Soc. 76 (1), 165–172. Faccenda, M., Burlini, L., Gerya, T., Mainprice, D., 2008. Fault-induced seismic anisotropy by hydration in subducting oceanic plates. Nature 455, 1097–1100. Fauzi, McCaffery, R., Wark, D., Sunaryo, Haryadi, P.Y.P., 1996. Lateral variation in slab orientation beneath Toba Caldera, northern Sumatra. Geophys. Res. Lett. 23 (5), 443–446. Fouch, M., Fischer, K., 1996. Mantle anisotropy beneath northwest Pacific subduction zones. J. Geophys. Res. 101 (B7), 15987–16002. Fouch, M., Fischer, K., 1998. Shear wave anisotropy in the Mariana subduction zone. Geophys. Res. Lett. 25 (8), 1221–1224. Gorbatov, A., Kennett, B.L.N., 2003. Joint bulk-sound and shear tomography for Western Pacific subduction zones. Earth Planet. Sci. Lett. 210 (3), 527–543. Gripp, A.E., Gordon, R.G., 1990. Current plate velocities relative to the hotspots incorporating the Nuvel-1 global plate motion model. Geophys. Res. Lett. 17, 1109–1112. Gudmundsson, O., Sambridge, M., 1998. A regionalized upper mantle (RUM) seismic model. J. Geophys. Res. 103 (B4), 7121–7136. Guest, W.S., Kendall, J.-M., 1993. Modelling seismic waveforms in anisotropic inhomogeneous media using ray and Maslov asymptotic theory: applications to exploration seismology. J. Can. Soc. Expl. Geophys., 78–92. Harmon, N., Forsyth, D., Fischer, K., Webb, S., 2004. Variations in shear-wave splitting in young Pacific seafloor. Geophys. Res. Lett., 31 (15). Helffrich, G., Wiens, D.A., Vera, E., Barrientos, S., Shore, P., Robertson, S., Adaros, R., 2002. A teleseismic shear-wave splitting study to investigate mantle flow around South America and implications for plate-driving forces. Geophys. J. Int. 149, F1–F7. Hiramatsu, Y., Ando, M., Tsukuda, T., Ooida, T., 1998. Three-dimensional image of the anisotropic bodies beneath central Honshu, Japan. Geophys. J. Int. 135 (3), 801–816. Jung, H., Karato, S., 2001. Water-induced fabric transitions in olivine. Science 293, 1460–1463. Kendall, J.-M., 2000. Seismic anisotropy in the boundary layers of the mantle. In: Karato, S., Forte, A.M., Liebermann, R.C., Masters, G., Stixrude, L., (Eds.), Earth’s Deep Interior: Mineral Physics and Tomography From the Atomic to the Global Scale. No. 117. Geophysical Monograph, 133–159. Kendall, J.-M., Thomson, C.J., 1992. Seismic modelling of subduction zones with inhomogeneity and anisotropy-I. Teleseismic P-wavefront tracking. Geophys. J. Int. 112, 39–66. Kennett, B.L.N., Engdahl, E.R., Buland, R., 1995. Constraints on seismic velocities in the Earth from traveltimes. Geophys. J. Int. 122 (1), 108–124. Koulakov, I., Jakovlev, A., Luehr, B.G., 2009. Anisotropic structure beneath central Java from local earthquake tomography. Geochem. Geophys. Geosyst., 10 (2). Lithgow-Bertelloni, C., Richards, M., 1998. The dynamics of Cenozoic and Mesozoic plate motions. Rev. Geophys. 36 (1), 27–78.

J.O.S. Hammond et al. / Physics of the Earth and Planetary Interiors 178 (2010) 189–201 Long, M.D., Silver, P.G., 2008. The subduction zone flow field from seismic anisotropy: a global view. Science 319, 315–318. Mainprice, D., 2007. Seismic anisotropy of the deep Earth from a mineral and rock physics perspective. Treat. Geophys. 2, 437–492. Mainprice, D., Barruol, G., Ben Ismail, W., 2000. The seismic anisotropy of the Earth’s mantle: from single crystal to polycrystal. In: Karato, S., Forte, A.M., Liebermann, R.C., Masters, G., Stixrude, L., (Eds.), Earth’s Deep Interior: Mineral Physics and Tomography From the Atomic to the Global Scale. No. 117. Geophysical Monograph, AGU, Washington, DC, pp. 237–264. Morley, A.M., Stuart, G.W., Kendall, J.-M., Reyners, M., 2006. Mantle wedge anisotropy in the Hikurangi subduction zone, central North Island, New Zealand. Geophys. Res. Lett., 33. Müller, R.D., Roest, W.R., Royer, J.Y., Gahagan, L.M., Sclater, J.G., 1997. Digital isochrons of the world’s ocean floor. J. Geophys. Res. 102 (B2), 3211–3214. Nakajima, J., Hasegawa, A., 2004. Shear-wave polarization anisotropy and subduction-induced flow in the mantle wedge of northeastern Japan. Earth Planet. Sci. Lett. 225, 365–377. Ohtaki, T., Kanjo, K., Kaneshima, S., Nishimura, T., Ishihara, Y., Yoshida, Y., Harada, S., Kamiya, S., 2000. Sunarjo, broadband seismic network in Indonesia-JISNET. Bull. Geol. Surv. Jpn. 51, 189–203. Pinero-Felicangeli, L., Kendall, J.-M., 2008. Sub-slab mantle flow parallel to the Caribbean plate boundaries: inferences from SKS splitting. Tectono-physics 462, 22–34. Polet, J., Silver, P.G., Beck, S., Wallace, T., Zandt, G., Ruppert, S., Kind, R., Rudloff, A., 2000. Shear wave anisotropy beneath the Andes from the BANJO, SEDA and PISCO experiments. J. Geophys. Res. 105 (B3), 6287–6304.

201

Savage, M.K., 1999. Seismic anisotropy and mantle deformation: what have we learned from shear wave splitting? Rev. Geophys. 37 (1), 65–106. Sdrolias, M., Muller, R.D., 2006. Controls on back-arc basin formation. Geochem. Geophys. Geosyst., 7 (4). Silver, P.G., 1996. Seismic anisotropy beneath the continents: probing the depths of geology. Annu. Rev. Earth Planet. Sci. 24, 385–432. Silver, P.G., Chan, W.W.J., 1991. Shear-wave splitting and subcontinental mantle deformation. J. Geophys. Res. 96, 16429–16454. Silver, P.G., Savage, M.K., 1994. The interpretation of shear wave splitting parameters in the presence of two anisotropic layers. Geophys. J. Int. 119, 949–963. Smith, G., Wiens, D., Fischer, K., Dorman, L., Webb, S., Hildebrand, J., 2001. A complex pattern of mantle flow in the Lau backarc. Science 292 (5517), 713–716. Syracuse, E., Abers, G., 2006. Global compilation of variations in slab depth beneath arc volcanoes and implications. Geochem. Geophys. Geosyst., 7. Teanby, N.A., Kendall, J.-M., Van der Baan, M., 2004. Automation of shear-wave splitting measurements using cluster analysis. Bull. Seis. Soc. Am. 94 (April 2), 453–463. Wheller, G., Varne, R., Foden, J., Abbott, M., 1987. Geochemistry of quaternary volcanism in the Sunda-Banda arc, Indonesia, and three-component genesis of island-arc basaltic Magmas 32 (1–3), 137–160. Wolfe, C.J., Solomon, S.C., 1998. Shear-wave splitting and implications for mantle flow beneath the MELT region of the East Pacific rise. Science 280, 1230–1232.