Complex mantle ﬂow in the Mariana subduction system: evidence from shear wave splitting

) for events < 250 km depth and APM-parallel for deeper events. Ray paths sampling west of the spreading centre show fast orientations ranging from arc-perpendicular to APM-parallel. Arc-parallel fast orientations characterize the southern part of the arc with variable orientations surrounding Guam. These results suggest that the typical interpretation of mantle wedge ﬂow strongly coupled to the downgoing slab is valid only at depths greater than ∼ 250 km and at large distances from the trench. We conclude that the arc-parallel fast orientations are likely the result of physical arc-parallel mantle ﬂow and are not due to recently proposed alternative lattice preferred orientation mechanisms and fabrics. This ﬂow pattern may result from along-strike pressure gradients in the mantle wedge, possibly due to changes in slab dip and/or convergence angles.


I N T RO D U C T I O N
The flow pattern of subduction zone mantle wedges is a matter of great importance for understanding the dynamics of subduction and backarc spreading processes. Models of subduction systems suggest that mantle wedge flow may be dominated by viscous coupling to the downgoing slab, producing flow directions parallel to the presentday absolute plate motion (APM) of the downgoing plate (McKenzie 1979;Ribe 1989;van Keken 2003). However, more complex flow patterns, including arc-parallel flow above the slab, may be produced by significant amounts of slab rollback (Buttles & Olson 1998), shearing or extension in the arc-parallel direction (Hall et al. 2000), changes in slab dip (Buttles & Olson 1998;Hall et al. 2000), or changes in convergence angle (Blackman & Kendall 2002;Honda & Yoshida 2005). Mantle flow may also align parallel to the orientation of maximum extension beneath the backarc spreading centre .
It is commonly assumed that observations of seismic anisotropy can provide strong constraints on mantle flow patterns. Both petrophysical studies of mantle xenoliths (Nicolas & Christensen 1987;Mainprice & Silver 1993) and laboratory studies of deforming rocks (Zhang & Karato 1995;Zhang et al. 2000) suggest that, in most cases, flow of mantle materials should produce a 'fast orientation' of upper-mantle anisotropy aligned close to or along the flow orientation, although some recent studies suggest that under a restricted range of conditions such as high stress and high water content, the 'fast orientation' may be perpendicular to the flow direction (Jung & Karato 2001;Karato 2003). The Mariana subduction system is a highly complex tectonic environment. Active volcanism throughout the forearc, arc and backarc  [2003][2004] OBSs that returned data (see Table 1), yellow triangles are [2001][2002] OBSs used in the current study. Large black vector is Pacific plate absolute plate motion (PAC APM). The backarc spreading centre and trench are sketched with thick red and black lines, respectively. Spreading direction is perpendicular to the spreading ridges, as indicated by small double-headed red arrows. Several islands are named in bold text. Blue square on Guam is GSN station GUMO. Earthquakes (circles) are colour coded by increasing depth: <100 km are red, 100-200 km orange, 200-300 km yellow, 300-400 km green, 400-500 km blue and >500 km violet. spreading centre, in addition to significant along-strike changes in slab dip and convergence angles, suggest complicated patterns of mantle flow throughout the region. However, the Marianas are commonly cited as having mantle flow parallel to present-day Pacific Plate APM directions based on prior studies at Guam (Xie 1992;Fouch & Fischer 1998). Here we analyse shear wave splitting results from two well-distributed seismograph deployments to characterize mantle flow patterns throughout the Mariana arc system. We find that several mechanisms are required to explain the observed shear wave splitting patterns and that the conventional corner flow model in a subduction zone is not a comprehensive description of mantle flow in the Mariana system.

Regional setting
The Mariana arc system encompasses a wide variety of tectonic settings, with active serpentinite seamounts in the forearc, an active island arc and backarc spreading centre, and an extinct fossil arc on the overriding Philippine Sea Plate (Fig. 1). Geochemical variations have been observed across and along the arc (Kelley et al. 2003;Pearce et al. 2005) and several physical features change rapidly along strike of the arc. Near Pagan, the slab dip is nearly vertical and the slab appears to penetrate the 660 km discontinuity (van der Hilst et al. 1991). Convergence is highly oblique at a rate of ∼4 cm yr -1 and the half spreading rate at the Mariana trough is ∼1.6 cm yr -1 (Kato et al. 2003). The southern part of the arc near Guam and Rota is significantly different from the northern part of the arc. The slab dips at ∼55 • , seismicity extends only to depths of ∼250 km, and convergence is roughly perpendicular to the trench at a rate of ∼6.5 cm yr -1 (Stern et al. 2003). The half spreading rate of the southern Mariana trough is ∼4.5 cm yr -1 (Kato et al. 2003). The Mariana trench does not experience rollback (Stern et al. 2003), as both the trench and the overriding Philippine Sea Plate both move westward at rates of ∼2.5-5 and ∼3-8 cm yr -1 , respectively (Heuret & Lallemand 2005;Martinez et al. 2000). Upper plate and trench migration rates are maximum near Saipan and Tinian Islands and backarc spreading rates are maximum near Guam (Martinez et al. 2000). Plate motions suggest that backarc deformation is primarily controlled by upper plate motions (Heuret & Lallemand 2005).

Previous work
Despite the intrinsic value of shear wave splitting measurements for addressing questions regarding the flow pattern in the Mariana mantle wedge, prior observations are sparse with data obtained only at the Global Seismic Network (GSN) station on Guam (GUMO) and a few ocean-bottom seismographs (OBSs) (Xie 1992;Fouch & Fischer 1998;Volti et al. 2006). These studies found fast orientations roughly parallel to the Pacific plate APM orientation and attributed their observations to the conventional interpretation of mantle wedge corner flow. The one exception to this is near the backarc spreading centre, where Volti et al. (2006) found fast orientations subparallel to the spreading direction.
The present study provides a detailed and thorough analysis of shear wave splitting patterns observed along and across the entire Mariana arc system. We use a relatively dense deployment of land and ocean bottom seismographs ( Fig. 1) and analyse local S phases with multiple data processing schemes to describe mantle flow patterns throughout the region. We find that while the commonly cited APM-parallel fast orientations are present in certain areas, arc-parallel fast orientations dominate at shallow and intermediate depths. Results show that mantle flow patterns vary along-strike of the arc and in the across-arc direction, indicating that a simple corner flow model for the region is not ubiquitous.

Data acquisition
Two seismic deployments provide a high-resolution data set for this study. Most data is from the 2003-2004 MultiScale Seismic Imaging Experiment in the Mariana Subduction System consisting of an 11-month deployment of 20 land broadband seismographs and 58 semi-broadband OBSs ( Fig. 1 and Table 1). The 78 stations were deployed during 2003 May-June and recovered during 2004 April-May. The land stations used Streckheisen STS-2 and Guralp CMG-40T sensors and were deployed on each island between Guam and Agrihan, with Reftek 72A-08 dataloggers and GPS timing. The 58 OBSs surround the deepest earthquake locations near Pagan Island and traverse the trench, forearc, island arc, and backarc spreading centre, extending across the West Mariana Ridge. Fifty OBSs used three-component Mark Products L4 sensors with 1 Hz natural period and modified amplifiers to extend long-period performance  . Fifteen of these OBSs were an older 16-bit model and 35 of were a new 24-bit design, and they were operated by Lamont Doherty Earth Observatory (these instruments are denoted by the 'OBS' prefixes in Table 1). The remaining eight OBSs used Precision Measuring Devices (PMD -WB2023LP) sensors with a low frequency corner at 0.03 Hz and were built by H. Shiobara at the University of Tokyo (denoted by the 'PMD' prefixes in Table 1). The 35 new U.S. OBSs stopped recording data ∼50 days after deployment due to a firmware error, eight U.S. OBSs were not recovered, and the Anatahan Island station had several power failures due to ash from the eruption covering the solar panels (see Pozgay et al. 2005). Several of the U.S. OBSs also failed to properly deploy the sensor to the seafloor. All other stations operated smoothly for the duration of the 2003-2004 deployment. The smaller data set is from nine OBSs with PMD sensors deployed during 2001-2002 across the arc system (denoted by the 'MAR' prefixes in Table 1) (Shiobara et al. 2005).

Data selection and processing
We use local earthquakes with depths greater than 80 km from the USA National Earthquake Information Center (NEIC) global catalog. Most events are between 100 and 300 km depth, but several events are between 300 and 600 km. We also investigated earthquakes located by the local deployment, but not identified at the NEIC due to their small magnitudes; only a few of these earthquakes were added to the final event list. Earthquake locations and depths were checked against relocations made using the local network. Unlike Volti et al. (2006) who used earthquakes as shallow as 19 km, we use only intermediate and deep earthquakes to eliminate complex ray propagation effects associated with shallow regional events and to provide ray paths that are generally propagating vertically.
Most earthquakes have S arrivals within the shear wave window only at stations located fairly close to the events, which limits our sampling range. In addition, S arrivals west of the spreading centre experience higher attenuation than at stations east of the trough axis, such that only the larger events have high quality S arrivals in the far backarc. We also investigated SKS arrivals, but no reliable splitting measurements were obtained due to low signal-to-noise ratio and poor distribution of events in the proper distance ranges.  Table 1). Filtering is often necessary to eliminate noisy portions of the spectrum, particularly the microseism peak near 0.2 Hz and/or high frequencies that may result from near receiver scattering. Some frequency dependence has previously been reported in results observed at the GUMO station (Fouch & Fischer 1998). Therefore, we analyse each shear wave with three filters (a 1-Hz lowpass filter and bandpass filters at 0.3-0.7 and 0.5-1.5 Hz) and visually inspect them to determine which frequency band produced the best result. In many cases, only one filter (usually 0.3-0.7 Hz) is appropriate. For larger events with very high signal-to-noise ratio across the entire frequency band, the three filters produce nearly the same result.
We orient OBSs using polarization data from air-gun shots and several Rayleigh waves and average 7-15 high-quality measurements for each OBS to ensure accuracy of the final orientations. Standard errors for the orientations range from 3 • to 7 • , except for four OBSs with standard errors of 13 • -17 • . In total, we compute 27 OBS orientations and use two orientations determined from Pwave polarizations by Volti et al. (2006). Note that we only orient OBSs with quality-A splitting results (see below).

Shear wave splitting analysis
Shear wave splitting is the process by which an S-wave travelling through a seismically anisotropic medium is split into a fast component and an orthogonal slow component. Two parameters describe the effects of anisotropy on the waveform: the polarization orientation of the fast component (φ) measured in degrees clockwise from north and the time lag (δt) between the two components (see Savage 1999 for a review). (Although 'fast direction' is the predominant term throughout the literature, we in most cases refer to a 180 • -ambiguous direction, and so use the term 'orientation'.) Several methods, including a covariance matrix method (Silver & Chan 1988;Silver & Chan 1991) and a cross-correlation method (Fukao 1984;Bowman & Ando 1987) have been developed to determine the optimum splitting parameters. The covariance method solves for the most linear particle motion after correcting for arbitrary splitting by minimizing the second eigenvalue of the covariance matrix. In contrast, the cross-correlation method finds the maximum crosscorrelation between corrected components.
There are significant differences in how shear wave splitting studies are carried out. Many studies use only one method without intercomparison with other methods, and some studies numerically select high quality results with a small amount of visual inspection. In addition, frequency dependence, noise, cycle skipping, and improper window lengths can bias shear wave splitting measurements. To make our measurements as systematic as possible and to reduce the effect of possible biases, we employ both above-mentioned analysis techniques, solve for the optimal window length and position and grade each result visually with numerical verification.
For each waveform, we compute splitting parameters with the covariance matrix method utilizing an automatic windowing scheme to solve for the most stable analysis window (Teanby et al. 2004). Initial window position begins 0.5 s before the arrival and ends 1.5 s after. The start and end times of the window are incremented by 0.35 and 0.3 s, respectively, such that absolute window lengths range from 2 to 10 s in length. Splitting parameters are computed for the covariance matrix method at each of the 175 window positions. The optimal solution is found using a hierarchical cluster analysis technique (Everitt 1993;Teanby et al. 2004). Subsequent to automatic windowing computation of the covariance matrix solution, we compute φ and δt again with the cross-correlation method (Fukao 1984;Bowman & Ando 1987;Smith et al. 2001). For this method, we solve for the maximum cross-correlation coefficient over rotations (0 • -175 • ) and time-shifts (±3 s) between the two split waves. We use the previously determined optimum window position for the calculation. However, since the optimal window lengths determined are preferentially short due to possible interference with subsequent phases (Teanby 2005), we also repeat the cross-correlation analysis with 3 s added to the end of the window. Fig. 3 shows the cross-correlation solution of −60 • and 0.5 s for one waveform, which is essentially identical (at these periods) to the covariance matrix solution of −59 • and 0.487 s. Fig. 4 shows the automatic windowing scheme for the same waveform.
After processing each waveform with the above methods, we manually assign grades of A, B, C, D and NULL to the covariance matrix solution and both cross-correlation solutions based on the signal-to-noise ratio, particle motion analysis, waveform clarity, correlation coefficient or energy contours, and cluster identification match to error contours (see Fig. 4 and Teanby et al. 2004). In this case, although the optimum automatic window appears short (Fig. 3), the result is stable over more than 50 per cent of the analysis  Fig. 3. In all panels, the optimal solution is shown with the red cross. Top: clusters of all 175 φ and δt solutions. Middle: error surface for covariance matrix solution. Note that the error contours are representative of the cluster plot. Bottom: φ and δt solutions for each window in the AUTOWIN regime. Note that the optimal solution (red cross) is where the minimal errors occur over the longest range of windows. windows (Fig. 4 bottom) and is considered a robust solution (Teanby et al. 2004). The cross-correlation method usually required a longer analysis window, in agreement with Teanby (2005). The application of both splitting methods also ensured accurate identification of null splitting results (see Wüstefeld & Bokelmann in review), which occur when the shear waves appear to be unsplit. This would arise in the absence of anisotropy along the ray path or because the natural polarization is nearly aligned with the fast or slow orientation. The two analysis methods give a different pattern of splitting results when analysing a null measurement. For an unsplit shear wave, the cross-correlation δt will be less than 0.1 s and φ will be 45 • off of the initial polarization angle, while the covariance matrix δt will be either large (3 s) or small (<0.1 s) and φ will be 0 • or 90 • off of the initial polarization.
The application of multiple filters sometimes resulted in more than one quality-A solution per waveform. When more than one filter for a given waveform produced a quality-A result, we systematically chose the 0.3-0.7 Hz bandwidth filter consistent with the dominant frequency of most S arrivals. Overall, we computed splitting parameters for 1908 waveforms (1304 event-station pairs). After selecting from multiple filters when appropriate, we found 308 quality-A cross-correlation measurements and 268 quality-A covariance matrix measurements recorded at 48 different stations (see Results column in Table 1 for stations with some measurements but no quality-A results).  Processing each waveform with two different methods yields some inconsistent results. Individually, 194 measurements have consistent cross-correlation and covariance matrix results within <20 • Similarity in φ and <0.25 s difference in δt (Fig. 5). Although the overall splitting patterns of these results are similar to handpicked quality-A results, we thoroughly investigate any user-bias on manual grading. The median difference between the 308 quality-A cross-correlation results and corresponding covariance matrix measurements is 11 • and 0.04 s. Similar splitting patterns are obtained when comparing handpicked quality-A cross-correlation solutions to those with correlation coefficient >0.75 and to solutions with excellent agreement between the two methods (<10 • similarity in φ and <0.25 s difference in δt). However, nearly 30 per cent of quality-A cross-correlation results have corresponding covariance matrix results of quality-D or null (Fig. 5). Since all tests show that handpicked quality-A cross-correlation results are robust and provide the largest number of measurements, we henceforth report only quality-A cross-correlation results. This observation of fewer usable covariance matrix solutions than cross-correlation solutions for noisy data has been reported previously (Restivo & Helffrich 1999;Long & van der Hilst 2005a).

Spatial averaging
One difficulty with the interpretation of shear wave splitting measurements is that the observed splitting can theoretically occur at any location along the ray path. For complex anisotropy, inadequate interpretation of the results may arise when plotting the splitting measurement at the station or event location or at the epicentral midpoint. Although shear wave splitting tomography based on the Christoffel equation is a possibility (Abt et al. 2006), the number of crossing rays in this region is insufficient to obtain a well-constrained result. Therefore, we employ a spatial averaging technique after the method of Audoine et al. (2004) to help characterize the overall anisotropy for regions of high complexity.
In this method, we associate each observation with the spatial volume sampled by the ray path and average the splitting results from all the ray paths sampling a given region. In practice, we assign the best-fitting cross-correlation splitting parameters to a set of nodes on a circular grid of concentric rings surrounding (and oriented perpendicular to) each ray path. These nodes define the region sampled by the ray path and are henceforth referred to as the 'ray path grid'. In order to spatially average the splitting parameters from all the individual ray paths, we subsequently superimpose a less-dense global coordinate grid (henceforth referred to as the 'global grid'). This global grid is oriented parallel to the geographic boundaries of our study area and is equally discretized in each of the X , Y , and Z directions. For each node of the global grid, we calculate the weighted average of all φ's and δt's from all the ray path nodes that reside within a box centred on the global coordinate node. Averages are computed at the surface global coordinate node (Z = 0 km) for all ray path segments shallower than 250 km (see Audoine et al. 2004 for details). Inverse distance weighting for each ray path node ensures that splitting parameters from longer ray paths will have smaller influence on many individual global nodes and shorter ray paths have larger influence on a lesser number of global nodes.
Averages are computed only for the northern part of the island arc, where we have the densest ray path coverage. However, we include ray paths from northern earthquakes to southern stations and from southern earthquakes to northern stations if they traverse the averaging region. The averaging scheme breaks down in areas of poor ray path coverage. Therefore, we de-emphasize any results based on single measurements (grey vectors in Fig. 6).
We performed several tests to ensure that the averaging results are robust with respect to the particular choice of averaging parameters. The only free parameters in the averaging scheme are the 'ray path grid' node spacing, 'global grid' node spacing, and starting location of the 'global grid'. Variations of the 'global grid' node spacing between 14 and 40 km showed some difference in individual average fast orientations, but the overall pattern remained the same. Averaging results were very similar when varying the 'ray path grid' node spacing (left column, Fig. 6). Similarly, varying the starting location of the 'global grid' altered individual averaged fast orientations (right column, Fig. 6), but again the overall pattern remained the same. After several trials, we found a 'ray path grid' spacing of 10 km and an overlying 'global grid' spacing of 25 km to provide the visual best match of averaged results to the raw splitting results. However, throughout the free parameter variations, we emphasize the small change of averaged splitting orientations and magnitudes in regions with a high density of crossing ray paths.

Northern region
Rose diagrams of all quality-A fast orientations plotted at each station show several trends (Fig. 7). The dominant fast splitting orientation is arc-parallel for most stations in the arc and for stations between the arc and backarc spreading centre, whereas the pattern becomes more arc-perpendicular in the far backarc region. We subset individual quality-A splitting results by event depth and plot them at their midpoint to provide further clarity. Results from earthquakes shallower than 250 km (Fig. 8) show roughly arc-parallel splitting orientations for ray paths between the island arc and spreading centre. Most measurements at the OBSs immediately north of Pagan are slightly oblique to arc-parallel with a few arc-perpendicular  measurements. Mid-latitude island arc stations (Anatahan through Guguan) show variable φ ranging between subparallel to APM and subparallel to the arc (Figs 7 and 8). West of the spreading centre and on the West Mariana Ridge, fast orientations are roughly parallel to APM. Shallow event results recorded at stations along the spreading centre show principally arc-parallel fast orientations near the main OBS line. OBSs near the northern part of the spreading centre show fast orientations roughly parallel to the spreading orien-tation. OBSs in the forearc show variable fast orientations ranging from arc-parallel to very oblique.
Events deeper than 250 km (Fig. 9) show somewhat different patterns compared to shallower events. Fast orientations are subparallel to APM along most of the arc, except at Pagan, where φ is subparallel to the arc. Stations at the spreading centre exhibit variable fast orientations, ranging from arc-parallel, APM-parallel, and APMperpendicular. Stations in the forearc record φ subparallel to the arc  Spatial averages for the depth range 0-250 km show roughly arcparallel fast orientations in the island arc and east of the backarc spreading centre (Fig. 10). We observe APM-subparallel φ near the spreading centre and towards the West Mariana Ridge. Some indication of fast orientations subparallel to APM is detected at mid-latitude island arc regions and the northernmost forearc region. These patterns clearly elucidate trends observed in the raw data (Figs 7-9).

Southern region
Nearly all measurements in the southern region have focal depths <250 km and two patterns of fast orientations predominate (Figs 11 and 12). Between Rota and Saipan, we record splitting measurements subparallel and oblique to the arc. In this region, the strike of the arc is ∼20 • and the average φ of these measurements is 3.1 • . We observe variable splitting orientations near Guam, but can make several general observations from the splitting map (Fig. 12). For example, events in the forearc show variable split-ting orientations, events north of the island show fast orientations subparallel to APM, and events southwest of the island show dominantly arc-parallel orientations. However, results are complicated and data coverage is sparse.

Depth extent of anisotropy
We investigate the possibility of depth-dependent lag times, as would be expected if anisotropy extends throughout the depth range studied, and has been noted for several other subduction zones (Yang et al. 1995;Fouch & Fischer 1996). Delay times throughout the northern region range from 0.1 to 2.1 s for focal depths <250 km and from 0.1 to 1.5 s for events ≥250 km. For both depth ranges, average lag times are 0.55-0.56 s. Southern region delay times range from 0.1 to 1.2 s with a 0.36 s average. There is a lack of systematic variation of δt with hypocentral depth (Fig. 13, right). However, there is a slight suggestion of a small increase of δt with path length that is more apparent for the southern stations (Fig. 13, lower left), but also noticeable at northern stations (upper left).
Percent anisotropy (k) is calculated for the ith ray path by where V S is an assumed S velocity of 4 km s −1 and d is the hypocentral distance. Arithmetic mean k AVE and individual k MAX are detailed in Table 2. Mean percent anisotropy for northern and southern events is 1.0 and 1.4 per cent, respectively. Nearly all southern events are between 80 and 200 km, therefore, we cannot discuss the depth extent of anisotropy there, except to mention that the mean value of 1.4 per cent must be characteristic of the upper 80-100 km. For the northern region, mean anisotropy for events shallower than 250 km is 1.3 per cent, which reduces to 0.5 per cent for events 250-600 km depth. The reduction in average anisotropy with depth is consistent with a small amount azimuthal anisotropy below 250-300 km in the mantle wedge (Fischer & Wiens 1996). We note that, since there is little or no variation in splitting time with depth or path length, the maximum percent anisotropy can be calculated with the minimum path length, which gives a maximum percent anisotropy of 1.6 per cent, which is slightly higher than the mean calculations. Figure 12. Comparison of results near Guam. All splitting results are plotted at the epicentre, oriented by φ, and scaled to δt. Solid black lines are results from this study and dashed lines are local S results from Fouch & Fischer (1998). Diamonds show qualitative NW (pink) or NE (yellow) fast orientations from Xie (1992). Circles are events used in this study. See text for discussion and Fig. 1 for base figure explanation.

Prior Mariana studies
Shear wave splitting fast orientations parallel to APM recorded at the GUMO GSN station on Guam are commonly cited as evidence for a corner flow dominant mantle flow regime. Previously reported fast splitting orientations at Guam range from −20 • to −80 • . With splitting times of less than ∼0.4 s (Fouch & Fischer 1998;Fig. 12). Xie (1992) also observed δt < 0.4 s and NNW-SSE fast orientations for events near Guam. The earlier GUMO study additionally found NNE-SSW φ for events southwest of the island. Fast orientations in this study from events southwest of Guam are roughly arc-parallel (Fig. 12), in agreement with the comparable qualitative measurements of Xie (1992) (yellow diamonds in Fig. 12); the later GUMO study did not report measurements from earthquakes in this area. Events located closer to the island show predominantly APM-subparallel fast orientations for all three studies and events in the forearc have variable orientations. Fig. 12 shows agreement of all Guam splitting measurements, with fast orientations dependent on earthquake location, such that the dominant pattern is of APMparallel φ near the island and arc-parallel φ southwest of the island. Fouch & Fischer (1998) Volti et al. (2006), using three of the same 2001-2002 OBS records used in this study and three OBSs from an earlier 1999-2000 deployment, found fast orientations subparallel to APM in the northernmost part of the arc and along the West Mariana Ridge and φ quasi-perpendicular to the trough axis near the spreading centre. We image fast orientations subparallel to APM on the West Mariana Ridge and to the west of the Mariana spreading centre in agreement with the previous OBS study. At the same time, our larger aperture array with greater station density allows us to image the rotation of fast orientations into the arc-parallel orientation at the backarc spreading centre, which was not observed in the previous studies.
In addition to shear wave splitting and seismic anisotropy investigations, other studies suggest preliminary interpretations about mantle flow patterns. Preliminary velocity and attenuation tomography results from the 2003-2004 deployment show separate low velocity and high attenuation regions beneath the volcanic arc and backarc spreading centre (Barklage et al. 2006;Pozgay et al. 2006). A small higher velocity and lower attenuation region separates these two zones at shallow depth ( < ∼ 50 km). This region might represent the boundary between different flow mechanisms-arc-parallel flow surrounding the island arc and APM-parallel flow beyond the backarc spreading centre. In addition, detailed geochemical analyses of the Mariana arc and trough suggest arc-parallel mantle flow between several centres of passive upwelling directly beneath the spreading centre (Pearce et al. 2005).

Other subduction zones
Several other subduction zones exhibit varying fast orientations that can be compared to observations in the Marianas. In Tonga, percent anisotropy is estimated at ∼1.5 per cent for the upper ∼200 km (Fischer & Wiens 1996), similar to calculations in this study, and fast orientations at stations near the arc are arc-parallel, while APMparallel values are observed farther into the backarc near the Fiji plateau (Bowman & Ando 1987;Smith et al. 2001). Southward infiltration of the Samoan plume above the retreating Pacific slab is one likely interpretation for this pattern, in agreement with geochemical evidence for plume infiltration (Turner & Hawkesworth 1998). Our observations in the Mariana arc show a similar pattern of arcparallel fast orientations near the arc and APM-parallel φ in the far backarc. In both subduction zones, the rotation from arc-parallel to APM-parallel occurs near the backarc spreading centre. The Mariana region does not show a slab window or the slab rollback thought to contribute to the Tonga dynamics, however the similarity of shear wave splitting measurements in the two different subduction zones may be indicative of a general pattern of upper-mantle flow in arcs with active backarc spreading.
Local S splitting results in Kamchatka show a maximum of 2.6 per cent anisotropy and highly variable fast splitting orientations (Levin et al. 2004). Fast orientations behind the arc are roughly arcparallel throughout the study area. However, the northeastern region shows oblique and scattered φ near the trench and arc, compared to relatively uniform arc-perpendicular φ in the southwest. Levin et al. (2004) conclude that corner flow is highly localized near the southern region and suggest that complex slab morphology and temporal slab dip variations might explain the heterogeneous fast orientations in the north. In our study, we also observe highly variable splitting orientations and time lags in some regions and observe similar percent anisotropy. However, we observe nearly opposite fast orientations compared to Kamchatka, with arc-parallel fast orientations dominating near the Mariana arc, but APM-parallel φ in the far backarc.
Splitting orientations from teleseismic S phases recorded in southern Japan are dominantly trench-parallel at the arc and trenchperpendicular in the backarc (Long & van der Hilst 2005b), while local S observations in central Honshu and Hokkaido show arc-parallel fast orientations in the forearc and APM-parallel φ beneath the arc (Nakajima & Hasegawa 2004;Long & van der Hilst 2005b). Fast orientations in the northern part of the Izu arc are oriented subparallel to both the convergence direction and the trench-parallel orientation (Anglin & Fouch 2005), whereas φ is consistently trench-parallel in the Ryukyu arc (Long & van der Hilst 2006). Fast orientations in this study are similar to that observed throughout parts of Japan and the Izu and Ryukyu islands and our estimate of percent anisotropy for northern events agrees well with similar calculations in Japan (Fouch & Fischer 1996).
In addition to those mentioned above, other subduction zones are also dominated by arc-parallel fast splitting orientations with some nearby APM-parallel component. Fast splitting orientations in the mid-latitude regions of the South American subduction zone are dominantly trench-parallel in regions of 'normal' slab behaviour and are subparallel to the strike of slab contours in regions where the slab flattens (Anderson et al. in review). Fast orientations are dominantly trench-parallel throughout the Central American arc and backarc and are only trench-perpendicular in the forearc (Abt et al. 2006), whereas the north island of New Zealand shows trenchparallel φ in the arc and forearc with trench-perpendicular φ only in the backarc (Marson-Pidgeon et al. 1999;Audoine et al. 2004;Styles et al. 2006). Geodynamic modelling results typically suggest APM-parallel mantle flow, however most arcs show arc-parallel fast orientations implying arc-parallel mantle flow. This study shows that although APM-parallel fast orientations are present in the far backarc and in a small region surrounding the island of Guam, arc-parallel fast splitting orientations predominate throughout the Mariana subduction system.

Fossil anisotropy in the subducting slab
Anisotropy might be controlled by fossil sea floor spreading (FSS) (Hess 1964) from ray paths travelling through the slab. Based on 45 • magnetic isochrons prior to subduction (Nakanishi et al. 1992), slab dips of ∼88 • and ∼55 • , and arc strikes of ∼0 • and ∼20 • for the northern and southern regions, respectively, and assuming the dominant fast anisotropy orientation in the oceanic lithosphere is isochron-perpendicular (Forsyth 1992), we would expect fast orientations approximately N-S or NNW-SSE in the northern region and roughly −20 • in the south. Although anisotropy in the slab could be distinguished by time lag differences between events in the upper or lower plane of the double seismic zone , data for such investigation is insufficient. Similarly, any anisotropy due to in-slab strain resulting from downdip extension or compression is not resolvable with this study.
We can, however, infer the presence or absence of slab anisotropy by other means. Along the arc and east of the trough, arc-parallel φ dominates for events <250 km (Figs 7-10). The arc-parallel orientation is similar to the FSS direction in the northern region, however we do not attribute the majority of these measurements to anisotropy in the slab because most ray paths travel primarily through the mantle wedge and sample a minimal region of the slab itself. Several deep northern events recorded at Pagan and at southern stations show fast orientations parallel to predicted FSS orientations with small δt and may result from slab anisotropy (dashed vectors in Fig. 9). These ray paths sample a significant amount of the slab and the solutions are oriented in the expected FSS directions and are in agreement with prior slab anisotropy solutions with average lag times of ∼0.25 s (Volti et al. 2006). APM-parallel fast orientations are recorded at Agrihan and Guguan.

Anisotropy in the crust
Considering the apparent lack of depth dependence on our splitting measurements, crustal anisotropy might seem a plausible source. However, we first note that most crustal anisotropy measurement are on the order of 0.05-0.2 s (Savage 1999), which is much smaller than our delay times (average δt is ∼0.5 s). Second, we note that crustal thickness beneath the Mariana Trough is only ∼5-7, ∼20 km thick beneath island arc stations, and ∼10 km thick beneath forearc stations (Takahashi et al. 2007). If the crust was the main source of anisotropy at OBS stations, where we observe an average 0.5 s delay time, percent anisotropy would be on the order of >28 per cent (based on 10 km thick crust and shear velocity of 4 km s −1 ), which is significantly larger than has been observed in crustal rocks (Babuska & Cara 1991;Crampin 1994). Although we cannot completely rule out a small effect of crustal anisotropy on our measurements, it cannot be a major source of splitting in our measurements.

Effect of water on the anisotropic slip system
Significant amounts of water and high shear stress may align the olivine a-axis perpendicular to the direction of maximum strain, resulting in a shear wave splitting fast orientation perpendicular to mantle flow (Jung & Karato 2001). Recent modelling results show that such changes in olivine slip systems are likely for slow subduction rates (< ∼3-4 cm/yr) and in areas with extremely high stresses ( 50 MPa), low temperatures (700-1000 • C), and presence of water, for example in forearcs (Kneller et al. 2005). Geodynamic modelling (Currie & Hyndman 2006) and heat flow (Blackwell et al. 1982) results suggest that the upper mantle beneath arc and backarc regions is characterized by temperatures that are too high (>1200 • C) to produce this alignment. Thus, 'B-type fabric' is not a good explanation for the extensive region of along-strike fast orientations extending from the arc to the backarc spreading centre in our results.
B-type fabric remains a possibility for the interpretation of forearc measurements, where the mantle is cold and has a high water content, as indicated by widespread serpentinization (Fryer 1996). However, if B-type fabric exists in the Mariana forearc, we would expect a small magnitude of anisotropy and the arc-parallel φ measurements would be interpreted as APM-parallel flow. Since island arc fast orientations are predominantly arc-parallel and that area is too hot for B-type fabric, this scenario would then require different physical mechanisms to invoke APM-parallel flow beneath the forearc and arc-parallel flow beneath the arc. However, due to the predominance of arc-parallel φ throughout much of the arc system and the significant magnitude of splitting observed in this region, we suggest that the forearc fast orientations are probably due to similar mechanisms that produce the arc-parallel fast orientations beneath the arc and are probably not due to B-type fabric.

Oriented melt pockets
Experimental studies show that a small amount of melt can weaken the overall a-axis alignment, such that abundant melt bands aligned in the direction of stress would rotate the a-axis 90 • to the dominant flow direction (Holtzman et al. 2003). If the flow direction in the vicinity of the backarc spreading centre were parallel to the spreading direction as inferred from a-axis orientations obtained by ocean-bottom electromagnetic studies (Baba et al. 2004), we would expect fast orientations aligned 90 • to the extension direction. We observe fast orientations that are dominantly ridge-parallel near the arc and switch to ridge-perpendicular immediately to the west of the spreading centre axis (Figs 7-10). This might suggest dominantly ridge-perpendicular mantle flow along the entire spreading centre with the along-strike fast orientation between the arc and spreading centre produced by aligned melt bands.
Although this hypothesis requires further study, there are several difficulties. This model would warrant significant melt porosity in the upper mantle over a wide depth range and over a large spatial extent; it is not clear whether sufficient in situ melt exists in the upper mantle to cause this effect on a wide scale. In fact, U-series disequilibria studies suggest that the porosity of the upper mantle is very low, both beneath oceanic spreading centres (Spiegelman & Elliott 1993;Lundstrom et al. 1998) and island arcs (Turner & Foden 2001). In addition, if widespread melt exists beneath the arc and between the arc and backarc spreading centre where we see alongstrike fast orientations, it is not clear why this mechanism would not also operate somewhat west of the spreading centre where we see APM-parallel fast orientations.

Slab-driven flow
A long-standing model for mantle flow in a subduction zone consists of mantle wedge material coupled to the downgoing slab producing APM-parallel fast orientations (e.g. van Keken 2003). The Mariana arc is often cited as a model for APM-parallel corner flow. However, we observe APM-parallel fast orientations only for deep events, at stations near the West Mariana Ridge, and close to Guam. We record a slight rotation of fast orientations towards APM-parallel at mid-latitude island arc stations, but fast orientations are dominantly arc-parallel. We conclude that APM-parallel corner flow is likely at 250 km distance and depth away from the trench, but is not dominant throughout the arc system.

Arc-parallel mantle flow due to slab morphology and convergence angle variations
One mechanism that may be important for producing arc-parallel flow in the Mariana wedge is arc-parallel flow in the asthenosphere. Pressure gradients that drive flow may be generated in a variety of ways, such as flow through slab tears or windows due to spatial (Smith et al. 2001) or temporal variations (Levin et al. 2004), slab rollback (Russo & Silver 1994;Buttles & Olson 1998), variations in downgoing slab morphology (Hall et al. 2000), or variations in convergence angle along-strike of the arc (Honda & Yoshida 2005).
Analytical models show that a change in slab dip of one degree over ∼100 km of along-strike distance will produce arc-parallel pressure gradients in the wedge corner that may, in turn, elicit arcparallel flow above the slab (Hall et al. 2000). Numerical modelling results based on this suggest that the magnitude of arc-parallel flow must be >25 per cent of the plate convergence rate to noticeably rotate olivine a-axes away from an APM orientation (Blackman & Kendall 2002). Such arc-parallel flow would decrease in magnitude with distance away from the wedge corner. In the Marianas, slab dip decreases by ∼40 • between Agrihan and Guam, a distance of ∼555 km. In accordance with modelling results of Hall et al. (2000), we would expect arc-parallel flow near the wedge corner resulting from large changes in slab dip and waning effects of arc-parallel flow with increasing distance away from the corner. Such a model is consistent with our observations of arc-parallel φ east of the spreading centre for events shallower than ∼250 km and APM-parallel φ for deeper events and at stations west of the spreading centre. If arc-parallel flow is due to slab dip variations, effects may only be imaged near the arc for shallow events (<250 km) with waning influence to the spreading centre (∼250-300 km laterally). Additionally, laboratory experimental results modelling a steeply dipping slab suggest that a-axis alignment may only predominate near the slab (Buttles & Olson 1998), further reinforcing a confined spatial extent of arc-parallel flow.

Arc-parallel pressure gradients-geodynamic modelling
To illustrate the qualitative viability of arc-parallel pressure gradients controlling the LPO structure within the mantle wedge, we create a 2-D finite-element model of mantle wedge flow using a stressand temperature-dependent viscosity structure. First, we calculate the viscosity structure of the mantle wedge, with the known inverse exponential relationship to temperature (Karato & Wu 1993), using prior thermal modelling methods (Conder et al. 2002;Conder 2005). The viscosities are lowest just outside of the 'cold nose' forearc, in the uppermost corner where temperatures and corner flow strain rates are highest. We impose an arc-parallel pressure gradient of a few kPa km -1 to that resultant viscosity structure and calculate the induced arc-parallel velocity structure (Fig. 14a) and the subsequent strain rate from the induced flow (Fig. 14b). The pressure gradient is an a priori assumption in this simple modelling, the purpose of which is to investigate whether a pressure gradient can produce the observed spatial distribution of splitting orientations. Thus, the magnitude of both the pressures and the flow velocity are meant to be representative and should not be as a precise measure.
We approximate conditions for a cross-section along our main OBS line near Pagan in our model with a slab dip of ∼75 • and a convergence rate of 4 cm yr -1 . To keep the methodology simple, we ignore the spreading centre in these models. While pressure gradients for either flow through slab tears or pressure-gradients from changes in slab morphology will be strongest near the corner, our 2.5-D assumption requires a constant pressure gradient over the entire wedge, noting that a variable pressure gradient would be possible with full 3-D modelling, but is outside the scope of this paper for our illustrative purposes. As such, arc-parallel flow rates far from the corner are likely overestimated relative to the arc-parallel flow rates within the corner. However, even with this overestimation, clear differences are apparent between the backarc and the uppermost corner of the wedge.
Two localized regions of high strain rate are visible in Fig. 14(b). One region is along the top of the slab and extends to depths 150 km. The other region is just below the overriding Philippine Sea Plate with a maximum depth extent of ∼150 km and dissipates towards the backarc. In this light, the high strain rates in the wedge corner dissipate between the island arc and spreading centre, becoming close to the ambient background. These regions of high strain rate due to arc-parallel velocity gradients illustrate regions of expected contribution to arc-parallel fast orientations and are a good match to our observations. Similarly, 3-D models with curved slabs also show strong trench-parallel stretching in the wedge corner (Kneller & van Keken 2006). Ray paths that travel through these high strain regions result in dominantly arc-parallel fast orientations, whereas ray paths travelling mostly through regions of low arc-parallel strain rate have roughly APM-parallel φ. In contrast, deep events recorded at northern island arc stations are nearly all subparallel to APM. Ray paths from deeper hypocentres will spend more time at depths that are largely unaffected by arc-parallel pressure gradients (Fig. 14b), and are more affected by corner flow in a direction parallel to APM.

C O N C L U S I O N S
Local S splitting results show complex patterns along and across the Mariana subduction system. Dense seismic arrays, multiple data analysis techniques, and comprehensive manual and numerical grading schemes provide a high-resolution image of seismic anisotropy patterns throughout the region. Shallow and intermediate depth events exhibit predominantly arc-parallel fast orientations in the forearc, island arc and backarc regions. Splitting measurements recorded at stations west of the backarc spreading centre and measurements from deep events show fast orientations subparallel to APM. We observe larger average splitting times in the northern region (0.55 s) compared to the southern region (0.36 s), but there is no increase in lag time with depth. Arc-parallel fast orientations for shallow and intermediate depth events are likely controlled by arc-parallel flow due to trench-parallel pressure gradients induced by a combination of slab dip and convergence angle variations and by arc-parallel extension. These causes would dominate near the wedge corner and have waning influence at greater distances. Our results indicate a transition at the backarc spreading centre from arc-parallel fast orientation on the east side to APM-parallel on the western side, similar to the transition observed in other subduction zones (e.g. Tonga arc). In addition, modelling results suggest that any arc-parallel pressure gradient induced by tectonic variations would not affect the far backarc and we see no evidence of arc-parallel fast orientations in this region. We conclude that simple 2-D corner flow does not dominate throughout the Mariana arc system and the slab-mantle system is only strongly coupled at ≥ ∼250 km depth and at large distances from the wedge corner.

A C K N O W L E D G M E N T S
Thanks to Gideon Smith, Nick Teanby, and Martha Savage for helpful comments on shear wave splitting methods and implementation. Theodora Volti provided two OBS orientations. Martha Savage, Vadim Levin, Karen Fischer, Megan Anderson, Maureen Long, Garrett Euler, Erica Emry, Moira Pyle and one anonymous reviewer all provided helpful advice and suggestions that greatly improved this paper. We thank numerous people, particularly Patrick Shore, Spahr Webb, Allan Sauter, Patrick Jonke, Juan Camacho, Joe Kaipat and Ray Chong, as well as the captain and the crew of the R/V Kaiyo, the R/V Wecoma, and the Super Emerald for assistance with deploying and recovering the seismographs. Land seismic instrumentation was provided by the PASSCAL program of the Incorporated Research Institutions in Seismology (IRIS) and the Lamont Ocean Bottom Seismograph Facility provided ocean bottom seismographs. This research was supported by the MARGINS program under National Science Foundation grant OCE0001938.