Heterogeneity Within the Subducting Pacific Slab Beneath the Izu-Bonin-Mariana Arc: Evidence from Tomography Using 3D Ray Tracing Inversion Techniques
M.S. Miller, A. Gorbatov, B.L.N. Kennett
Abstract
Distinct zones of seismic heterogeneity along the Izu-Bonin-Mariana arc have been investigated using complimenting regional bulk sound, shear wave speed, and P-wave tomographic images. The distribution of seismic anomalies and the inferred geometry of the subducting Pacific plate have been modeled in unprecedented detail using joint tomography and new P-wave model using a 3D inversion algorithm. The use of 3D ray tracing techniques and smaller cell parameterizations have greatly enhanced the resolution of gradients, therefore, the models show much more detail about the structure and physical properties of the subduction zone. The well-defined features from the multiple wave speed images are used to elucidate the distinct morphology change between slab beneath the Izu-Bonin (horizontal slab) and Mariana (vertical) arcs and the distribution of physical properties in the mantle and the subducting oceanic lithosphere. Changes in physical properties within the slab tear at the southern end of the Izu-Bonin arc, identified as a ″gap or thinning″ in the tomographic images, could be the result of the distortion of the Pacific plate as its shape transforms between near horizontal to near vertical, a decrease in the rigidity and strength of the lithosphere, the subduction of the Marcus-Necker Ridge, change in subduction velocity, or a combination of all these factors. Strong, slow anomalies exist both, in the mantle wedge, above and below the slab beneath the Izu-Bonin arc but are not present beneath the Mariana arc. At the junction between the Izu-Bonin and Mariana arcs it appears that the there are two separate pieces of the Pacific plate: the torn slab north of 26° N and a buckled near vertically dipping slab south of the tear. The transition between the two distinct slab morphologies coincides with the location of the Ogasawara Plateau and the trench.
1. Introduction
Seismic tomography has developed into one of the most effective and significant sources of information in modern geophysical research. One of the objectives of seismic tomography has been to understand tectonic processes and the interior of the earth but it can also provide details on the physical properties of the mantle [1–3]. For example, the physical properties of subducting slabs have been revealed through the comparison of bulk sound speed and shear wave speed anomalies from joint P and S data inversions [2]. Such joint inversions provide sharp definition of subduction zone features through an inversion exploiting the selection of P and S wave arrival times that have similar ray paths and can be used to uncover seismic heterogeneity on regional and global scales.
While many features of subduction zones in the Western Pacific were first identified in the 1990s [1,3–7], new detailed images of the variation of bulk sound and shear wave speeds and a new P-wave inversion can reveal the heterogeneous structure and, moreover, physical and chemical properties. An area of specific interest in this study lies at the ″junction″ of the Izu-Bonin and Mariana arcs (between 20° and 34° latitude) in Fig. 1. A distinct change in morphology and seismic property beneath Izu-Bonin was identified by Miller et al. [8] and was suggested to be related to the distortion of the Pacific plate as its shape transforms from near horizontal to vertical, coupled with the subduction of the Ogasawara Plateau. Farther south, along the Mariana arc the morphology of the subducting Pacific plate has a near vertical orientation [1,3,8–12]. We have imaged the discontinuous nature of the slab structure between the Izu-Bonin arc and the Mariana arc with the comparison of bulk sound, shear wave speed, and new P-wave inversions, which provides insight into both the change in physical properties and morphology of the subducting Pacific slab and the mantle wedge along Western Pacific margin.
2. Regional Setting
The Izu-Bonin arc forms the northern half of the Izu-Bonin-Mariana convergent margin in the Western Pacific, which has an extent of almost 1500 km from Japan to the end of Volcano Islands (Fig. 1). The Pacific plate, which is late-Cretaceous to early-Jurassic in age at the trench, is subducting beneath the Philippine Sea plate at a rate of approximately 9 cm⁄yr along the Izu-Bonin arc and 5 cm⁄yr along the Mariana arc [13,14]. Besides the relatively rapid movement of the Pacific plate towards the northwest, the Izu-Bonin trench has migrated from south to northeast with clockwise component of rotation since 48 Ma. In the last 17 my the trench-trench-trench triple junction of the Philippine, Pacific, and Eurasian plate has changed direction to move westward by about 50 km [13]. It is inferred that relatively young lithosphere was subducted beneath the Izu-Bonin arc during fast trench migration in the Oligocene to Mid-Miocene [11,13,15]. van der Hilst and Seno [11] proposed that the combination of young oceanic lithosphere and rapid trench migration produced a shallow subduction angle, with the slab being laid down in the transition zone in the southern portion of the Izu-Bonin arc. It was suggested that temporal changes in slab morphology gradually propagate in space from south to north based on the change in dip of the Wadati-Benioff zone and P-wave tomographic images from the early 1990s.
The variation in geometry of the subducted slab along the arc is one of the most intriguing aspects of the region. A north-south traverse approximately parallel to the Izu-Bonin-Mariana arc reveals the angle of the subducted slab becomes increasingly steeper; from a dip of approximately 45°, to horizontal along the 670 km discontinuity, to nearly penetrating vertically into the lower mantle beneath the Mariana arc, as seen in the distribution of seismicity in Fig. 1. The variation in subduction angle has been clearly imaged using seismicity, residual sphere analyses, discontinuity topography and tomography [3,5—7,16–19]. However, the details of the slab structure in the region between these two different geometries are not well understood and a tectonic model that clearly describes the change in morphology has not been established.
3. Tomographic Controls
Both the regional body wave joint tomographic inversion of Gorbatov and Kennett [2] and the new P-wave tomographic inversion were produced from the same arrival-time data. The detailed joint inversions used an inversion algorithm introduced by Kennett et al. [1] and then adapted by Gorbatov and Kennett [2] to include 3D ray tracing. This includes the trajectory of seismic ray propagation between source and receiver through the three-dimensional structure of the earth, which improves the resolution of gradients and strong variations in wave speeds. This approach has been successfully utilized for global and regional tomography [2,10,20].
P and S arrival-time data with the same source and receiver from the global catalogue of Engdahl et al. [21] were used. Single rays with both P and S readings were picked for events and stations within the study area to optimize data coverage and to ensure that the final images could be directly compared. Mantle structure in the Western Pacific region was parameterized into a non-overlapping grid of 5° x 5° with 16 layers ranging from 35 to 200 km down to a total depth of 1600 km. The average velocity values in the layers in each of the datasets are given in Table 1. The study area consisted of a grid of 19 layers with cells 0.5° x 0.5° in the uppermost mantle, 1° x 1° in the transition zone and lower mantle, and 2° x 2° beneath continents to a depth of 1500 km. The final dataset, using the ak135 model [22] as a reference, contains 900,000 pairs of P and S ray paths. This nested iterative approach resulted in P and S models that were then used in a non-linear joint tomographic inversion for bulk sound and shear wave speed inversions [2].
A separate P-wave tomographic inversion was obtained with a similar nonlinear scheme using the same dataset as in the joint tomography and with the same inversion parameters, but using a standard tomographic formulation. Although the inversion schemes used for joint tomography and the P-wave tomography are different, the three inversions can be compared with careful consideration of the smoothing (damping) parameters.
The potential resolution of the images can be estimated by checking the data coverage using density of ray paths. This quantity is measured by the sum of the lengths of the ray segments within each 0.5° cell in the grid. Fig. 2(A–D) represents the density of ray paths in kilometers for P- and S-wave seismic rays in both cross section and plan view. As seen in the images, the Izu-Bonin-Mariana region and its surroundings have dense ray path coverage both beneath the arcs, from the transition zone to the lower mantle, and throughout most of the volume beneath the Philippine Sea plate.
4. Tomographic Image Anomalies
4.1. Previous Models
From P-wave imagery van der Hilst et al. [7] and Fukao and Obayashi [4] (among others) have shown that the Izu-Bonin region has a prominent fast anomaly corresponding to the subducted Pacific slab lying horizontally on the 660 km discontinuity. This seismic velocity anomaly has been identified in multiple tomography inversions, but Widiyantoro et al. [3] used P and S wave arrival time data with similar path coverage to investigate this anomaly. Significant disparities in the configuration of P and S images were suggested to be differences in the flexural rigidity between the horizontal and the steeper dipping portions of the slab.
Gorbatov and Kennett [2] undertook a regional study using joint tomography from P and S arrival times to produce bulk sound and shear wave speed images, which confirmed the results in Widiyantoro et al. [3]. Beneath the Izu-Bonin arc, the steeply dipping slab is imaged as a strong shear wave feature but a weak signature in bulk sound, whereas the ″stagnant″ portion in the transition zone has a clear bulk sound speed expression. Strong shear wave anomalies associated with subducting slabs are primarily evident in subduction zones with lithosphere older than 90 Ma entering the trench, as below the Mariana arc, but further investigation is needed to describe the physical characteristics that cause these differences.
4.2. Multiple Inversion Models
Bulk sound and shear wave speed inversions from Gorbatov and Kennett [2] are visualized in EarthVision to provide a more detailed analysis of the data. The joint inversion results are supplemented with a new P-wave speed inversion, using 3D ray tracing and the same set of P arrival times, to compare to past P-wave inversions and as a complimenting dataset. The three sets of tomographic images for the Izu-Bonin-Mariana region (10–50° latitude and 110–150° longitude) are sensitive to different aspects of the structure of the subduction system.
Cross sections through all three models sliced along 31° latitude illustrate the subducting Pacific slab, which can be traced along the fast velocities down dip to the 660 km discontinuity and then horizontally through the transition zone (Fig. 3). A prominent bulk sound speed feature for the structure laying in the transition zone, as detected by Gorbatov and Kennett [2] (Fig. 3A), and a weak signature in the shear wave speed image at the same position (Fig. 3B). The P-wave model has a strong signature in both the horizontally positioned slab and the slab in the upper mantle but becomes weak at approximately 410 km (Fig. 3C). Features present in the new P-wave model derived with 3D ray tracing is consistent with past inversions based on inversions from 1D reference models [3,5–7,17,23] but are sharper and stronger than the previous models.
4.3. Southern Izu-Bonin and Northern Mariana Arcs
Weak deviations from the ak135 [22] reference model in a P-wave model in the southern end of the arc were initially noted in Miller et al. [8] and are confirmed in the new P-wave inversion that utilizes 3D ray tracing to greatly improve gradient resolution by including the three-dimensional earth structure between source and receiver (Fig. 4C). This region plus the northern portion of the Mariana arc was studied in detail using the joint inversions to understand the change in physical properties as detected by seismic velocities. The models shown in Fig. 3 are sliced along a traverse perpendicular to the trench, A–A’ in Fig. 1, to produce the cross sections in Fig. 4 and show earthquake hypocenters taken from the NEIC⁄USGS catalogue for events with magnitude greater than 4.5 between 1967–1995. The bulk sound and Pwave images show a distinct decrease or gap in the fast velocity anomaly between 400–450 km, where the shear wave image displays a ″thinning″ of the slab at the same position (Fig. 4A–C). This thinning or gap as seen along A–A’ (Fig. 1) is interpreted as a slab tear as in Miller et al. [8]. In contrast, the mantle wedge above (and west of), as well as below (or east of) the subducting slab is imaged with a strong, slow anomaly in the P-wave and shear wave speed images (Figs. 3 and 5). The subducting plate, interpreted as fast velocity anomalies, appears to be bordered by slow anomalies beneath the Izu-Bonin arc.
South of the anomalous area of southern Izu-Bonin, along the northern extent of the Mariana arc, the subducting slab penetrates vertically into the lower mantle Fig. 4(D–F) and has much different characteristics. The transect along B–B’ (shown in Fig. 1) illustrates the steeply dipping subducting Pacific plate and Wadati-Benioff zone, which is in agreement with previous studies [3,6,7,10]. The narrow band of positive anomalies interpreted as the Pacific plate subducting beneath the Mariana arc is consistent in all three inversions. The difference in available resolution in the images between A–A’ and B–B’ (as muted colors) is due to the lack of stations in the area and therefore dependent on sources rather than receivers. The sharp definition of width of the slab is also constrained by the width of cells used in the inversion matrix (approximately 100 km thick). Despite the less pronounced anomalies beneath the Mariana arc, there are only slightly slow velocities surrounding the slab, unlike in the upper mantle farther to the north.
The layer anomaly maps of bulk sound, shear wave speed, and P-wave speed give a different perspective of the anomalies in the transition zone (Fig. 5). The differences between the patterns of bulk sound and shear wave speed heterogeneity at 425 km depth beneath the Izu-Bonin region are apparent in each of the models. The bulk sound image details a ″break″ in the fast velocity anomaly at approximately 27.5–31.5° N latitude and 138–141° E longitude with a similar feature in the P-wave image (Fig. 5A,C). A ″thinning″ of the fast velocity at this same position in shear wave speed compliments the two previous images (Fig. 5B). The dimensions of the gap can be estimated as 4° (N–S), 3° (E–W), and 175 km in depth from Figs. 4 and 5. The intrinsic resolution of the tomography in this region is in the order of 50 km, which is smaller than the estimated size of the anomaly. The striking difference in the two geometries, the fast horizontal anomaly on the transition zone, the near vertical dip of the slab beneath Mariana, and the ″slab tear″ beneath southern Izu-Bonin are apparent in all three inversions both in map and cross section.
4.4. Three-dimensional Visualization Model
The volume of interest, which includes the Izu-Bonin and Mariana arcs as shown in an inset in Fig. 1, is visualized as a solid model in Fig. 6. The 3D visual model was created by picking the maximum gradient of velocity in the new P-wave model to define the top and bottom of the subducting slab, which compensates for the lower resolution in the southern portion of the model. This 3D schematic representation of the subducting plate beneath the Izu-Bonin arc illustrates the anomalous area as a missing block and the varying slab morphology as it transforms from steeply dipping to horizontal within a small distance between 25° and 30° N latitude. The horizontal portion of the slab has a sharply defined terminus at approximately 26° N, which is located where the high velocity anomalies become weak before disappearing completely.
In contrast to the geometry of the Izu-Bonin slab, the subducted slab beneath the Mariana arc is nearly vertical along the length of the arc. In Fig. 6 the southern region of the model from 25–15° N illustrates the steeply dipping Pacific plate as arcuate in shape similar to the trace of the trench location (Figs. 1 and 5). Although the Pacific plate beneath the Marianas penetrates into the lower mantle at a steep dip, it appears to have a crumpled or buckled morphology in the lower mantle, unlike the stagnant slab of Izu-Bonin.
Earthquakes collected from the NEIC catalogue for events between 1967–1995 are plotted with the 3D model and color-coded by depth (Fig. 6B). The deep seismicity plots within the slab subducting beneath the Mariana and Izu-Bonin arcs, but a cluster of hypocenters is found within the slab tear beneath southern Izu-Bonin (Figs. 4 and 6). The cluster of earthquakes appearing to be at the point of flexure in the subducting slab may be the result of accommodation of considerable amount of strain due to the large variation in geometry along arc [8]. The focal mechanisms for these events are predominately extensional with a NNW–SSE strike, which is consistent with mechanical breakup of a non-elastic, bending slab [24]. The general pattern of seismicity in this region was initially characterized by Lundgren and Giardini [18] as extensional, producing a shearing of the deep slab to the south along strike (parallel to the trench), which is also in agreement with the bending of the slab in this region.
5. Discussion and Conclusions
Previous studies have outlined the differences between bulk sound and shear wave speed inversions [1] and the implications of the heterogeneity of physical properties in subducted slabs [2,12]. Inferences about tectonic history and cause of differences in the physical properties have been explained by variations in age, rigidity, and strength of the subducting plate in the mantle [3,14,25–28], but it is mostly likely to be a combination of factors.
Between the southern end of the Izu-Bonin arc and the northern portion of the Mariana arc there is a dramatic change in the geometry and physical properties within the slab and mantle transition zone imaged in the three tomographic inversions for different wave speeds. Results of the joint bulk sound and shear wave speed tomographic images confirm the existence of a change in properties at 325–500 km depth beneath the southern end of the Izu-Bonin arc. Miller et al. [8] illustrated this change using a P-wave inversion with cells 1° by 1° down to 1600 km using 1D ray tracing. The three tomographic inversions used in this study are greatly improved due to the use of 3D ray tracing techniques and smaller cell parameterizations. All of the inversions show a decrease in the fast velocity perturbations from reference model ak135 [22] at depths between 325–500 km beneath southern Izu-Bonin (Fig. 5A-C). The reduced velocities for all three inversions within this anomalous region are evidence for a distinct change in property within the subducting Pacific slab. As the anomaly is clear in bulk sound, shear wave speed, and P-wave models it can be inferred that heterogeneity within this region is the result of a true change in physical properties within the slab.
The combination of strong bulk sound signature, but a weak signature in shear wave speed along the horizontally lying deflected slab, and the strong shear anomaly above the hinge point suggests that there is a change in rigidity within this region. Reidel and Karato [28] suggested that weak zones can develop in old (and cold) slabs in association with phrase transformations and reduced grain size at the tip of a metastable olivine wedge. The slab could become less rigid as it bends and is deflected in the transition zone at approximately 350–400 km depth. This configuration is consistent with the area of weakening and region of anomalous rheology within subducted slabs argued by Riedel and Karato [28]. As the strength within the slab decreases and the conditions are suitable, the slab may begin to break and tear. This breaking is evident in the ″tear or gap″ in the P-wave, bulk sound, and shear wave speed images (Fig. 5A-C) and the position of an earthquake cluster at the point of flexure (Figs. 4 and 6).
Another feature that is evident in all of the inversions is a strong, negative anomaly to the west of the slab and a similar slow velocity to the east (Figs. 3B-C and 5B-C) in P-wave and shear wave speed. The bright, high amplitude feature in the mantle wedge could be related to a region of lower viscosity. Low viscosity, due to water and volatiles released from the subducted slab into the mantle wedge, has been proposed by geochemical studies to explain arc volcanism [29,30]. The region of slow velocity and therefore low viscosity could thus be due to the presence of water in the mantle wedge [31–35]. Low viscosity in the mantle wedge has also been used as a constraint in mantle convection and subduction zone initiation models for the Izu-Bonin-Mariana subduction system [25,26,36,37] but has also been considered to modify flow patterns and lead to larger velocities in the mantle wedge [38]. Models presented in Hall et al. [26] have a wide zone of extended lithosphere directly above the top of the slab allowing for hot mantle to be in direct contact with the slab, in which dehydration could produce large amounts of partial melting and dehydration into the mantle wedge. Since there is a similar low velocity feature under the slab beneath the Izu-Bonin arc, but is not present along the Mariana arc, viscosity changes may not be a suitable explanation.
Grain-size reduction or low viscosity in the mantle wedge may not be the sole possibility for this decrease in rigidity in the slab; in fact it is likely to be due to a combination of origins. The change in style of subduction and age of the subducting oceanic lithosphere could also be a factor. As proposed by van der Hilst and Seno [11] the subducting slab was laid down in the transition zone due to the large oceanward migration of the Izu-Bonin trench (approximately 1000 km) from 30 to 17 Ma. There has been little trench retreat since 17 Ma, which could remove the conditions for slab deflection onto the transition zone, leading to stress build up at the point of flexure and a reduction in the strength of the slab at 325–500 km depth, where the slab transitions from a 45° dip and horizontal on the 660 km boundary.
An important facet in the distortion and heterogeneity of the Pacific slab at depth is the collision of the Ogasawara plateau (and Marcus-Necker Ridge) with the Izu-Bonin trench and its probable continuation at depth. As the Pacific plate is subducting obliquely in a north-northwest orientation at about 9 cm⁄yr [11,14] the extension of the Marcus-Necker Ridge on the downgoing slab can be estimated. Burbach and Frohlich [39] suggested that material at 30° N would enter the trench at 24° N but failed to consider the impact of the aseismic ridge on the Pacific plate. Using their approach the position of the anomalous seismic velocities in the slab beneath the Izu-Bonin arc between 28° and 32° N latitude can be traced back to a place of entry at approximately 24° N, which is where the Ogasawara Plateau collides with the trench (Figs. 1, 5 and 6). Interestingly, the current location of the Ogasawara Plateau is at a transition point between the two distinct pieces of the Pacific plate (Fig. 6). Hsui and Youngquist [40] use an approximation of 15 Ma as the time of initial collision of the Ogasawara Plateau and the trench based on tectonic reconstructions [41] and paleomagnetic data [42]. These papers suggest a mid-Miocene collision of the aseismic ridge at the trench, which would also correspond to the slowed migration (or possible reversal) of the trench-trench-trench triple junction at 17 Ma that has been advocated in tectonic reconstructions of the Philippine Sea plate [13–15,43–45]. The tectonic history of events involving the mid-Miocene collision, the subsequent subduction of the Marcus-Necker Ridge, and change from trench retreat to trench advance suggests that it has had a significant role in the evolution of the morphology of the subducted Pacific plate.
The geology of the Ogasawara Plateau also contributes to the interpretation of the heterogeneity and morphology of the region. A seismic survey of the aseismic ridge categorized the feature as a topographic high (2000–3000 m) above the surrounding seafloor consisting of a series of seamounts (guyots) of Cretaceous age [46]. The flat summits of the guyots were interpreted to be filled with lagoon sediments of considerable thickness (1.2 s in two-way travel time), rimmed by reef. If similar sediment-covered seamounts extend down into the subduction zone, the carbonate-rich sediment could be the source of the volatiles that have been released into the mantle wedge [30] to create the low-velocity zone found beneath the southern Izu-Bonin arc [8].
The visualization of the Pacific plate subduction beneath the Izu-Bonin and Mariana arcs shows a geometry that differs somewhat from previous models that describe the shape of the subducting Pacific slab [2–4,6,7,19,20]. The anomalous area in the plate is well imaged, but another interesting feature is the crumpled slab beneath the northern part of the Mariana arc and the stagnant slab beneath Izu-Bonin (Fig. 6). The two different geometries appear to be two almost separate pieces of the Pacific plate. The Mariana slab appears to be anchored and has ″buckled″ in the lower mantle as it penetrated through the 660 km discontinuity as suggested in earlier studies [11]. In contrast to the near vertical geometry beneath the Marianas, the slab beneath Izu-Bonin is interpreted as having a relatively sharp edge, which is consistent with seismic structure interpretation from comparing regional distance waveforms with synthetic seismograms [47]. The extent of the fast velocities associated with the stagnant slab dissipates around 26° N, which has led to the interpreted geometry of the slab at depth. The slightly rounded shape of the edge of the stagnant slab could be due to erosion of the plate from mantle flow similar to the slab edge beneath Kamchatka [48].
The two pieces of the plate seem to have had different responses to the tectonic history of the Western Pacific margin, which could be due to the ongoing subduction of the extinct Marcus-Necker Ridge and the Ogasawara Plateau. The plateau appears to be in a position that could have ″sliced″ the two distinct morphologies. The collision of the large topography high of the Ogasawara Plateau in the mid-Miocene could be the final event that has affected both the Izu-Bonin trench motion and the morphology of the subducting Pacific plate at depth.
The combination of bulk sound, shear wave speed, and P-wave tomography is strong evidence for heterogeneity within the subducting slab beneath the southern Izu-Bonin arc. The change in seismic properties within the slab is related to the distortion of the Pacific plate as its shape transforms from steeply dipping to near horizontal, as seen in the 3D visualization model, and the decrease in rigidity and strength of the oceanic lithosphere. The combination of a weak, old slab, the change of subduction style from fast trench rollback to a slow, opposite rate of motion, and the collision of the Marcus-Necker Ridge has produced the heterogeneity seen in the tomographic images.
Acknowledgments
We would like to thank Scott King, Mike Gurnis, and an anonymous reviewer for their constructive comments and suggestions that greatly improved the paper. We would also like to thank Christian Noll and Jim Dunlap for their support and reviewing early versions of the manuscript.
References
For complete list of References see PDF download.
