Abstract

Subduction initiation may unfold via different pathways in response to plate strength, plate age, and driving mechanism. Such pathways influence volcanism on the overriding plate and may be preserved in the sequence of erupted volcanic products. Here, we parameterize melting in a mechanical model to determine the volcanic products that form in response to different subduction initiation modes. We find that with a mode of continuous initiation with infant-arc spreading, the foundering of the subducting slab and water release from the slab govern a succession from basalts with compositions similar to mid-ocean-ridge basalts (MORB) to boninites. The modeled transition from MORB-like to boninite composition typically occurs within a few million years. When plate strength is reduced, the subducting slab tends to segment, with extensive melting occurring when the slab breaks; most melting occurs close to the trench. When plate strength increases, subduction initiation becomes continuous without infant-arc spreading; such a mode leads to a limited, very low degree of melting occurring during a long interval of plate convergence before subduction initiation starts, although extensive melting near the trench is still possible when subduction initiation starts after a protracted period of plate convergence (∼10 m.y.). If the subduction initiation is driven by constant stresses, such as through ridge push, the slab subducts rapidly in response to continuous acceleration of the plate under action of the far-field push; significant melting, including boninite eruption, can be generated within a few million years with no trench migration. Based on the tectonic and volcanic evolution, these different modes may be applicable to the initiation of the Izu-Bonin-Mariana arc (infant-arc spreading and a sequence from MORB-like to boninites), the New Hebrides arc (slab segments in the upper mantle), the Puysegur Trench in New Zealand (scarce distribution of volcanism and no infant-arc spreading), and the Aleutian Trench (strong volcanism and no infant-arc spreading).

INTRODUCTION

Subduction initiation occurs at the interfaces of adjacent plates, such as transform faults between two oceanic plates, possibly induced by a change in plate motions (McKenzie, 1977; Gurnis et al., 2004; Stern, 2004). The detailed initiation processes can be deciphered from studies of sedimentary and volcanic rock sequences of specific nucleating margins (House et al., 2002; Sutherland et al., 2009; Reagan et al., 2010), and then placed within the context of the regional kinematics of the major plates (e.g., Sharp and Clague, 2006). Computational models have been used to investigate the mechanics of initiation, including the tectonic evolution of the overriding and subducting plates, in response to induced subduction initiation by external forces (Toth and Gurnis, 1998; Hall et al., 2003; Gurnis et al., 2004; Thielmann and Kaus, 2011; Leng and Gurnis, 2011). However, the volcanic expression on the overriding plate in response to unfolding subduction initiation remains unexplored.

The volcanic expression on the overriding plate shows distinctive patterns for the several known subduction initiation events. Extensive magmatism characterized by depleted (tholeiitic) and ultradepleted (boninitic) magmas occurred during the initiation of the Izu-Bonin-Mariana subduction zone (Stern and Bloomer, 1992). Detailed, in situ sampling of the Mariana forearc near Guam has revealed a volcanic stratigraphy in which basalts similar to mid-ocean-ridge basalts (MORB) were the first to erupt in the nascent arc and were quickly followed by boninites and then by normal arc lavas within several million years (Reagan et al., 2010). The same volcanic sequence is found ∼1500 km to the north near the Bonin Islands (Ishizuka et al., 2011). In contrast, the formation of the Aleutian island arc has been accompanied by a high rate of magma production (Jicha et al., 2006), but no infant-arc spreading is observed. For the Puysegur subduction zone, i.e., the northern segment of the Macquarie Ridge complex between the Australian and Pacific plates, plate convergence has been under way since 15 Ma, and the downdip tip of slab has reached 150 km depth (Sutherland et al., 2009), but the evidence for volcanism is scarce, with only small isolated seamounts found on the Pacific plate, notably Solander Island (Reay and Parkinson, 1997).

These different volcanic patterns may have been caused by divergent evolutionary pathways in response to different mechanical conditions. Driven by an imposed constant velocity, subduction initiation may occur catastrophically after a finite amount of plate convergence (Hall et al., 2003; Gurnis et al., 2004). In this case, the older and denser plate founders with the development of infant-arc spreading and rapid trench rollback, all occurring within a few million years. Once the foundering starts, the theoretical model of induced subduction is similar to Stern and Bloomer’s (1992) conceptual model of self-nucleating subduction. The difference between the two is the prediction of compression and rapid uplift of the overriding plate in the induced model just before slab foundering, while the self-nucleating model has no uplift and compression. Gurnis et al. (2004) argued that relic structures on the overriding plate showed uplift just before the Izu-Bonin-Mariana arc formed, but this possibility needs further testing. The mechanical conditions under which this mode will occur have been clarified by Leng and Gurnis (2011). When we systematically increase the coefficient of internal friction or decrease the rate of plastic weakening of the plates, the mode of subduction initiation induced by imposed velocity is found to change from (1) segmented initiation to (2) continuous initiation with infant-arc spreading and eventually to (3) continuous initiation without infant-arc spreading (Leng and Gurnis, 2011). Here, the terms “segmented” versus “continuous” describe whether or not the slab breaks at depth during initiation. These three different modes have respectively been applied to the initiation of the New Hebrides, Izu-Bonin-Mariana, and Puysegur subduction zones (Leng and Gurnis, 2011). Subduction initiation can also be induced by a constant stress, such as by a ridge push force. The constant stress case may apply to the Aleutian subduction zone, initiated by subduction of the small Kula plate (Seton et al., 2012). In this case, the mode of subduction initiation is similar to continuous initiation without infant-arc spreading, but it requires a much shorter plate convergence time (Leng and Gurnis, 2011).

The relationship between volcanism and the evolutionary pathway of a nascent subduction zone can be addressed through coupling between a quantitative melting model and a mechanical model; such models can provide insights into the connections between the chemical composition of magmatic products and the evolution of temperature, pressure, mantle water contents, and surface plate kinematics. Using thermodynamic principles, petrological models such as pHMELTS (Asimow et al., 2004) have been developed to compute the extent of melting and major- and trace-element compositions of mantle-derived melts. This approach may not be efficient for application to large-scale geodynamics models (Hebert et al., 2009a), as pHMELTS needs to compute melt extent and composition for tens of thousands of points for every time step in a geodynamic model, which typically runs for tens of thousands of time steps at significant computational cost.

Alternatively, computations of extent of melting and magma compositions can be parameterized (Katz et al., 2003; Kinzler and Grove, 1992a, 1992b; Parman and Grove, 2004), with several controlling parameters constructed to describe complicated melting processes. These parameters are calibrated to match laboratory experiments on melting. From the parameterized approach, extent of melting (Katz et al., 2003), major-element composition of melts for both lherzolite (Kinzler and Grove, 1992a, 1992b) and harzburgite melting (Parman and Grove, 2004), and trace-element compositions can be determined. Although the parameterizations are largely empirical and may not reflect the exact physical and chemical process of melting, they are simple and computationally efficient. More importantly, they are able to largely reproduce the melting results of laboratory experiments and provide us with a first-order understanding of melting processes during subduction initiation. Here, we collect the relevant parameterizations from the literature, build a geodynamic melting model, and obtain a complete picture of the way in which subduction initiation evolves with time and how this process can influence volcanic sequences.

A COUPLED GEODYNAMIC MELTING MODEL FOR SUBDUCTION INITIATION

Geodynamic Model

We used the two-dimensional (2-D) finite-element model Ellipsis (Moresi et al., 2003) to study the visco-elasto-plastic behavior of the lithosphere and mantle during subduction initiation. Ellipsis uses Lagrangian particles to accurately track material properties and deformational history. Mechanically, the geodynamic model is the same as the subduction initiation model in Leng and Gurnis (2011), but here we add melting. Governing equations and model details can be found in Moresi et al. (2003) and Leng and Gurnis (2011). We only provide an overview of mantle rheology and the model domain.

We used an incompressible Maxwell material for visco-elasticity. The deformation rate is the sum of an elastic part and a viscous part: 
graphic
where forumla, forumla, sij, G, and η are the strain rate, time rate of change of deviatoric stress, deviatoric stress, shear modulus, and viscosity, respectively, and i and j are spatial indices. The viscosity is non-Newtonian and temperature dependent: 
graphic
where η0, forumla, forumla, n, E, R, T, and T0 are the reference viscosity, the second invariant of the deviatoric strain rate tensor, reference strain rate, strain exponent, activation energy, gas constant, and absolute and reference temperatures (both in Kelvin) (Table 1). In addition, the maximum and minimum viscosities are limited with ηmax and ηmin, respectively (Table 1). We use a Drucker-Prager yielding model for the plasticity: 
graphic
where τy, μ, P, and C are the yield stress, coefficient of friction, pressure, and cohesion. Cohesion is the strength of the material at zero normal stress. The values of μ and C are linearly reduced with accumulated plastic strain to reproduce the weakening of mantle material with increasing plastic strain (Buck and Poliakov, 1998; Poliakov and Buck, 1998): 
graphic
 
graphic
where μ0, ɛp, ɛf, C0, and Cf are the initial coefficient of friction, accumulated plastic strain, reference plastic strain, initial cohesion, and minimum cohesion. Here, ɛp represents the total accumulated plastic strain after the yield stress is reached, and ɛf represents a threshold value of plastic strain. When the accumulated plastic strain is larger than ɛf, no further weakening is accounted for. Table 1 shows the invariant model parameters used. It has been suggested that μ0 and ɛp determine the mode of subduction initiation (Leng and Gurnis, 2011). We vary these two parameters herein to produce different types of subduction initiation.

The computational domain is 1050 km in width by 300 km in depth (Fig. 1). A 50-km-thick compressible “sticky air layer” is inserted above the normal solid domain to simulate a free surface, which may be important for subduction (e.g., Crameri et al., 2012) (Fig. 1). The air layer has zero density but has a viscosity ηmin. In nominal computations, 128 × 64 elements are used with 16 tracking particles initially inserted uniformly in each element; additional higher resolutions were used by Leng and Gurnis (2011), and the resulting dynamics were resolved with this resolution. The subducting and overriding plates have different initial thermal ages and are divided at the box center with a weak zone (17.5 by 17.5 km) (Fig. 1). The overriding plate age is 2 m.y., and the subducting plate age tp is varied according to different tectonic settings. There is no heat flux at the left and right boundaries. The initial temperature field is computed from the half-space cooling model with a fixed reference temperature T0 at the bottom boundary. The top, bottom, and left boundaries are free slip. We use two different boundary conditions to induce subduction initiation. The first is an imposed constant velocity at the right boundary. The vertical velocity is zero, and an imposed inward velocity Vx is applied from 0 to 100 km depth. To conserve mass, an outward Vx is applied from 200 to 300 km depth. We use a tanh function from 100 to 200 km depth for velocity change from inward to outward. The second condition uses imposed stresses in the subducting plate. We put a constant stress at a fixed position, 700 km from the left boundary, from the surface down to 50 km depth (Fig. 1). With this stress boundary condition, the right boundary is set to be free slip, and another weak zone (dashed rectangle in Fig. 1, 50 by 100 km) is placed at the right end of the subducting plate to detach it from the boundary.

Parameterized Model of Mantle Melting Process

We constructed a parameterized melting model that can be conveniently implemented into Ellipsis. Each of the 16 tracking particles initially inserted in each element represents a mantle parcel and tracks the extent of melting, phase proportions, and major- and trace-element compositions of the parcel. The initial major- and trace-element compositions (in weight percent) for these particles are obtained from the depleted MORB mantle (DMM) composition (Workman and Hart, 2005). This DMM composition is slightly depleted with respect to bulk silicate earth, as required by the complementary trace elements and isotopic signatures of continental crust and the residual upper-mantle MORB source, but it is actually a fertile lherzolite composition. We only consider eight major elements and replace others with SiO2 (Table 2). The initial mantle phase partitions are (in weight percent): olivine 57%, orthopyroxene 28%, clinopyroxene 13%, and spinel 2% (Workman and Hart, 2005). Each time a particle moves to a new position, we compute the extent of melting from the temperature, pressure, water content, and weight percent of modal clinopyroxene based on the parameterized method of Katz et al. (2003). If the extent of melting is larger than that tracked by the particle, new melt is generated and extracted. The extent of melting and residual composition of the particles are then updated.

We consider that newly generated melts are instantaneously extracted from the mantle and vertically transferred to the surface to form volcanic rocks. At the plate surface, between 175 and 525 km in horizontal direction (Fig. 1), we uniformly place 40 sampling points, which store the thickness and compositions of erupted volcanic rocks and subsequently translate them with the surface (plate) velocity. The effective widths of these points are the same as the element widths beneath the sampling point. These sampling points move in response to plate motion and stretching. If the distance between two adjacent sampling points becomes too large, >87.5 km, 10 more sampling points are uniformly added between these two. For the newly generated melts, the trace-element compositions are computed from constant partition coefficients for each trace element (Table 3) (Workman and Hart, 2005) with the batch melting equation (Schilling, 1966): 
graphic
where CL, CB, D, and F are the element concentration in liquid, bulk composition of the source mantle (updated in time, so this formulation in fact describes near-fractional melting), partition coefficient, and melt extent. A nonmodal trace-element partitioning model would be superior, but for the present exercise, this simple model is sufficient; note that the bulk D values used are appropriate to describe the average melting process observed in mid-ocean-ridge settings. The major-element compositions are computed with either the lherzolite multiple saturation surface (when clinopyroxene is present) based on the parameterized method of Kinzler and Grove (1992a, 1992b) or with the harzburgite multiple saturation surface (when clinopyroxene has been completely consumed) based on the parameterized method of Parman and Grove (2004). Therefore, the thickness and compositions of the surface volcanic rocks associated with mantle melting processes can be computed from our model. We output these results every 0.25 m.y. After the melt extraction, the phase proportions of the source mantle are also computed and updated according to lherzolite melting (Kinzler and Grove, 1992a, 1992b) or harzburgite melting (Parman and Grove, 2004).

We computed a test of the major-element compositions for different extents of melting. With pressure at 1.0 GPa and no water content, we increased the temperature from 1500 to 1680 K with an interval of 15 K. New melts were extracted each time temperature increased. The composition of melts varied when the extent of melting increased from 0.8% to 26.6% (Table 2). MgO continuously increased with extent of melting. When the extent of melting was larger than 13%, all clinopyroxene was consumed, and the melting process changed from lherzolite to harzburgite melting (Table 2). For lherzolite melting, SiO2 and Al2O3 decreased while CaO increased with extent of melting. However, for harzburgite melting, CaO and Al2O3 decreased while SiO2 increased with extent of melting. We used the oxide weight percents of SiO2 and MgO to divide the resulting volcanic rocks into four compositional groups: Group A represents the very low-degree-of-melting rocks with SiO2 > 52% and MgO < 10%; group B is the low-degree-of-melting (normal) MORB-like rocks with SiO2 < 52% and MgO < 13%; group C is the high-degree-of-melting MORB-like rocks with SiO2 < 52% and MgO > 13%; and finally, group D represents boninites with SiO2 > 52% and MgO > 13%.

For our melting model, the reference temperature T0 plays a key role for computing the extent of melting and compositions. We determined T0 using the following approach. Since the age of the overriding plate in our model is just 2 m.y., strong melting and a layer of surface volcanic rock are generated at the first time step when we start the model. We consider that this layer resembles the oceanic crust formed from MORB. When T0 is 1355 °C, the thickness of this volcanic layer is 7 km, the same as the typical thickness of oceanic crust (Orcutt et al., 1984). We therefore use this T0 in our melting model. Note that in the following results and discussion, this initial non-subduction-related basement volcanic layer is removed from the melting results.

The subduction of oceanic crust brings water into the mantle wedge, which significantly alters the melting process. The modeling of chemical signatures of subduction fluid and water transport in the mantle wedge can be complicated (e.g., Iwamori, 1998; Hebert et al., 2009b). Here, we use a simple method to model the effects of water and slab-derived fluids on volcanism associated with subduction initiation. We consider that there is a thin wet layer at the top of the subducting plate. The mantle zone extending vertically 60 km above the subducting plate is considered “wet,” with a water content of 0.3%. To simulate refertilization of the mantle wedge by slab-derived fluids, the concentrations of trace elements in this wet mantle zone are updated at each computation time step to a mixture of 99.7% residual source composition and 0.3% of slab-derived fluid composition. The slab-derived fluid composition (Table 3) was compiled from the results of Stolper and Newman (1994). With this method, we obtain a wet mantle zone that moves with the subducting plate and has the chemical signature of subduction fluid.

VOLCANIC SEQUENCES ASSOCIATED WITH DIFFERENT SUBDUCTION INITIATION MODES

Mode of Continuous Subduction Initiation with Infant-Arc Spreading

We start with case A01, for which we set the subduction plate age tp and the imposed velocity Vx to be 82 m.y. and 4 cm/yr, respectively. With plate strength parameters μ0 = 0.2 and ɛp = 0.5, we obtain a subduction initiation mode of continuous initiation with infant-arc spreading (Leng and Gurnis, 2011). It takes several million years to underthrust the old plate beneath the young one (Fig. 2A). After the negative buoyancy of the subducting plate overcomes the resistance from plate bending and friction, it begins to founder rapidly, leading to significant asthenospheric upwellings (up to several tens of centimeters per year), strong infant-arc spreading, and trench retreat (Fig. 2B).

To illustrate the effects of water on the volcanism, we first computed melting without the inclusion of water. During the first several million years of plate convergence, there is no melt generated, and the mantle melting region is confined to the initial melting region induced below the thin overriding plate (Figs. 2A and 2D). When the slab founders and subduction begins, infant-arc spreading causes adiabatic melting beneath the spreading center and forms volcanic rocks (Figs. 2B and 2E). Initially, these rocks are mostly within group B, the low-degree-of-melting (normal) MORB-like rocks, which move from the spreading center toward the trench with plate motion (Fig. 2E). As the slab continues to founder, the spreading center moves with trench retreat, and more extensive melting generates rocks of group C, high-degree-of-melting MORB-like rocks (Figs. 2C and 2F). The accumulated thickness of the volcanic rocks reaches ∼8 km.

When the effects of water are included, the melting extent is enhanced in the mantle wedge due to water injection (Figs. 3A, 3B, and 3C). The wet mantle zone with 0.3% of water content occurs above the slab, extending 60 km upward (Fig. 3D). The enhanced melting leads to much larger thicknesses of volcanic rocks above the wet zone, up to 15 km (Fig. 3E). The horizontal distribution of erupted compositional groups at different time indicates that all the eruptions start at ∼4.2 m.y. when slab begins to founder (Fig. 3F). Notably, the mantle source entering the mantle wet zone is residual to the melt extraction beneath the infant-arc spreading center and is remelted due to water injection into this zone. This process characteristically generates boninites, which are erupted at the surface near the trench (Figs. 3E and 3F). The sites with boninite eruption move with trench retreat. On closer inspection of the rock sequence at the sampling point that ends up at position 500 km in Figure 3E (between the spreading center and the trench), we find a prominent stratigraphic sequence for the volcanic rocks, changing from group B to C and then to D within a few million years (Fig. 4A). The major-element compositions of these rocks show the transition from normal MORB-like rocks to boninites (Fig. 4B). The Mg# increases with this transition from 70.2 to 76.5 (Fig. 4B). The trace-element compositions for group B, which are normalized by primitive mantle (PM) compositions (Sun and McDonough, 1989), show similar patterns to the normal MORB rocks (Fig. 4C; compared with Hofmann, 1997). For groups C and D, the trace elements are more depleted for compatible elements but are enriched for elements that are rich in slab-derived fluid (Fig. 4C).

Starting with case A01 and holding plate strength parameters μ0 and ɛp constant, we vary the subduction plate age tp between 82 and 42 m.y. and the imposed velocity Vx between 4 and 2 cm/yr and obtain eight additional cases. The results show that the subduction initiation modes for all nine cases are continuous initiation with infant-arc spreading. The surface volcanic sequences typically include boninite and show the transition between different groups (Fig. 5). However, two cases with Vx = 2 cm/yr and tp = 62 and 82 m.y. generate no boninites. To investigate the missing boninite formation, we place a tracking particle initially at 35 km in depth and 350 km in horizontal distance for both case A01 and case A01a, which has Vx = 2 cm/yr and tp = 62 m.y. The paths of the tracking particle during the whole process for these two cases (black curves in Figs. 6A and 6B) show that mantle material beneath the spreading center moves toward the slab when subduction initiates with slab foundering and trench retreat. For case A01, the depleted material beneath the spreading center reaches the surface of the slab. The injection of water from the slab remelts the depleted material, and boninite is generated. For case A01a, however, the combined effects of decreased plate velocity and increased plate age delay the starting of infant-arc spreading. The subducting plate reaches the bottom of the model domain shortly after infant-arc spreading starts. The depleted material beneath the spreading center does not reach the wet zone above the slab before the slab is stopped by the bottom boundary. Therefore, no boninites are observed (Fig. 5).

Segmented Subduction Initiation Mode

With plate strength parameters μ0 = 0.1 and ɛp = 0.2, we have a weak plate, which leads to a segmented mode of subduction initiation (Leng and Gurnis, 2011). We set the subduction plate age tp and the imposed velocity Vx to be 42 m.y. and 2 cm/yr, respectively (case A02). After several million years of plate convergence, the subducting plate bends and infant-arc spreading occurs (Fig. 7A). However, unlike in the mode of continuous subduction initiation with infant-arc spreading, the weak plate cannot sustain the slab pull forces, and the subducting plate breaks to form a slab segment (Fig. 7A). As plate convergence continues, another slab segment breaks from the plate (Fig. 7B). During the whole process, significant melting occurs, and large volumes of volcanic rocks form. However, the weak plate leads to rapid vertical descent of the broken slab segments, and the wet zone due to water release from the slab segments only exists in a small region (Figs. 7A and 7B). Therefore, extensive melting and thick volcanic sequences including boninites only occur close to the trench; this zone moves with the trench as it rolls back (Figs. 7C and 7D). Starting with case A02 and holding plate strength parameters μ0 and ɛp constant, we vary the subduction plate age tp between 82 and 42 m.y. and the imposed velocity Vx between 4 and 2 cm/yr, and the resulting subduction initiation modes are all segmented, with similar volcanic expression as shown for case A02.

Mode of Continuous Subduction Initiation Without Infant-Arc Spreading

High plate strength leads to continuous subduction initiation with neither infant-arc spreading nor trench retreat (Leng and Gurnis, 2011). With tp = 42 m.y. and Vx = 2 cm/yr, we increase plate strength to μ0 = 0.6 and ɛp = 0.4 (case A03). Until 6.8 m.y. after plate convergence begins (Fig. 8A), very few volcanic rocks form (less than 0.5 km in thickness from 0 to 6.8 m.y.) (Fig. 8C). Most of these rocks are of group A, very low-degree-of-melting rock (Fig. 8C). The typical major-element concentrations of group A are characterized by high SiO2 and Al2O3, but low MgO (Fig. 4B). After several more million years of convergence, the subducting plate bends, with relatively weak asthenospheric upwelling beneath the overriding plate (Fig. 8B). The upwelling causes melting, but the melting region is confined to a small mantle region where water is released from the slab surface (Fig. 8B). Since there is no trench migration, the site of extensive volcanic eruption does not move, which causes thick accumulation of volcanic rocks at this site after eruption starts (Fig. 8D). We also hold plate strength parameters μ0 and ɛp constant for case A03 and vary the subduction plate age tp between 82 and 42 m.y. and the imposed velocity Vx between 4 and 2 cm/yr. For the case with tp = 42 m.y. and Vx = 4 cm/yr, the resulting subduction initiation mode is the same as case A03, which leads to similar volcanic expression, except that the corresponding time scales are approximately halved due to increased plate velocity. For the cases with tp = 82 m.y. and Vx = 2 and 4 cm/yr, due to the increased subducting plate age, the subduction initiation modes change to continuous with infant-arc spreading (Leng and Gurnis, 2011) and have similar volcanic expression as shown for case A01.

The subduction initiation in cases A01 through A03 is induced by an imposed constant velocity on the subducting plate. Subduction initiation may also be induced by a constant stress applied on the subducting plate (Toth and Gurnis, 1998; Leng and Gurnis, 2011). With tp = 42 m.y., and the same plate strength parameters as we chose for case A01, μ0 = 0.2 and ɛp = 0.5, we change from the imposed velocity boundary condition to the imposed stress boundary condition (Fig. 1) for case A04. A tectonic stress of 80 MPa is applied within the old plate from the surface down to 50 km depth (Fig. 1), which is equivalent to a tectonic force of 4 × 1012 N/m. This stress boundary condition leads to a subduction initiation mode similar to case A03. With such a force, the subducting plate accelerates during plate convergence. It only takes 1.5 m.y. for the plate to reach the bottom of the model domain (Figs. 9A and 9B). Note that here the mantle melting extent due to the initially thin overriding plate is larger than in previous cases, because the mantle has less time to cool. There is no infant-arc spreading associated with this subduction initiation process. The surface volcanic rock patterns are similar to case A03 (Figs. 9C and 9D). Moderate melting occurs beneath the overriding plate, while extensive melting and boninites form close to the trench (Fig. 9D). There is little horizontal transport of depleted mantle beneath the overriding plate for this case. However, since the mantle has less time to cool, the water released from the slab causes large melting extents close to the trench, which can remelt the depleted mantle due to the initially thin overriding plate (cf. Figs. 9A and 9B for the melting contour). Therefore, a large volume of boninite is generated.

The typical ridge push force is estimated to be ∼4 × 1012 N/m (McKenzie, 1977). Here, we apply such a force uniformly between 0 and 50 km depths, which leads to a constant stress of 80 MPa. Applying this force uniformly between 0 and 67 km depths leads to a stress of 60 MPa. This stress causes subduction initiation with volcanism similar to case A04, except that the time for the slab to reach the bottom of the model domain becomes slightly longer, ∼3.0 m.y. If this force is applied uniformly between 0 and 100 km depths, the resulting stress of 40 MPa is smaller than the initial cohesion C0 and cannot induce subduction initiation.

GEOLOGICAL APPLICATIONS OF DIFFERENT MODES OF SUBDUCTION INITIATION

The Izu-Bonin-Mariana subduction zone initiated at ca. 50 Ma with strong infant-arc spreading and voluminous boninite eruptions (Stern and Bloomer, 1992; Reagan et al., 2010) and may be associated with the change in direction of the Pacific plate from northward to westward (Sharp and Clague, 2006). Case A01, having a mode of continuous subduction initiation with infant-arc spreading, provides an evolutionary pathway similar to the sequence of tectonism and volcanism of the juvenile Izu-Bonin-Mariana subduction zone. After several million years of plate convergence, the subducting plate founders and causes strong infant-arc spreading and volcanic eruptions. In particular, this mode of subduction initiation produces volcanic rocks that are similar to the observed sequences, including the order, thicknesses, and major- and trace-element compositions (Fig. 4) (Reagan et al., 2010; for observed units, see fig. 4 inIshizuka et al., 2011). The concentrations of highly incompatible trace elements in boninite obtained from our model are larger than those observed at Izu-Bonin-Mariana forearc region (Fig. 4C), indicating that mantle source in Izu-Bonin-Mariana region may be more depleted than the averaged mantle composition considered in our model. In the volcanic sequences obtained (Figs. 4 and 5), the rock composition typically changes from MORB-like to boninite (if any) within a few million years. The thickness of these sequences varies from several to more than 10 km. We consider such volcanic sequences to be characteristic of the volcanic expression related to the mode of continuous subduction initiation with infant-arc spreading. They indicate the occurrence of slab foundering, infant-arc spreading, and remelting of depleted mantle due to water release from the slab. Such sequences are observed at both the Mariana forearc (Reagan et al., 2010) and the Bonin Islands (Ishizuka et al., 2011), suggesting that a significant portion (more than 1500 km in width) of Pacific plate foundered during Izu-Bonin-Mariana subduction initiation and caused infant-arc spreading and strong volcanism along the whole Izu-Bonin-Mariana forearc region. One limitation is that the age of the modeled overriding plate is only a few million years, which would be consistent with oceanic lithosphere near a spreading center. However, the Mariana and Bonin sections are about ∼1500 km from each other, and there is no evidence that both sections are adjacent to a spreading center perpendicular to a nascent trench.

The New Hebrides subduction zone may have initiated with a segmented subduction initiation mode similar to case A02. The New Hebrides subduction zone initiated through a polarity reversal at ca. 10–12 Ma (Greene et al., 1994), but a slab segment is found beneath the North Fiji Basin (Hamburger and Isacks, 1987; Richards et al., 2011). Richards et al. (2011) suggested that this slab segment occurred by the tearing away of the deeper portion from the shallower portion of the slab at ca. 5 Ma, or ∼5–7 m.y. after the subduction initiation event. The subducted plate at the New Hebrides Trench is continuous with the crust of the South Fiji Basin (Seton et al., 2012). Therefore, it is possible that this plate formed within a backarc basin and remains mechanically weak due to high volatile concentrations. Our model suggests that if the New Hebrides subduction initiation is in the segmented mode, we expect to observe thick volcanic products close to the trench that possibly include boninites.

Plate convergence has been occurring continuously since ca. 15 Ma at the Puysegur subduction zone south of New Zealand, which nevertheless shows neither infant-arc spreading nor significant volcanism. Therefore, this subduction zone may have been initiated with the mode of continuous subduction initiation without infant-arc spreading, as in case A03. This mode can explain the observed scarce distribution of volcanism in the Puysegur subduction zone so far. Moreover, the volcanic rocks found in Solander Island in this zone contain mostly adakite and andesite (Reay and Parkinson, 1997). No basalts are reported. There is little to no garnet involved in the signature. Therefore, these rocks are probably derived from a low degree of partial melt from the wedge (Foley et al., 2011). This is consistent with results from case A03, which produces rocks mostly in group A, very low-degree-of-melting rocks, for the first several million years. Nevertheless, the volcanism in case A03 will eventually be extensive after ∼10 m.y. of plate convergence. We thus suspect that this indicates possibly strong volcanic eruptions in the Puysegur subduction zone in the future, although the region of eruption would be confined to a narrow region close to the trench.

The Aleutian subduction zone initiated at approximately the same time as Izu-Bonin-Mariana, ca. 50 Ma. The subducting plate in the Aleutians is the Kula plate, a small plate confined by the Kula-Pacific ridge a few thousand kilometers to the south (Seton et al., 2012). This plate geometry can lead to a constant ridge push force from the Kula-Pacific ridge during the formation of Aleutian subduction zone. Therefore, Aleutian subduction initiation may be driven by a constant stress from mid-ocean-ridge push. Relatively strong volcanic eruptions have been reported during the formation of the Aleutian subduction zone (Jicha et al., 2006), while no infant-arc spreading has been inferred for the Bering Sea. There are no reports of boninite along the Aleutian arc. However, if case A04 is applicable, this suggests that a large volume of boninite may exist in the forearc, on the overriding plate close to the Aleutian trench, and might be found by drilling in this area. With neither slab foundering nor infant-arc spreading, a volcanic sequence that changes from MORB-like to boninite is not found for this mode (Fig. 9D).

DISCUSSION AND CONCLUSION

Our model indicates that the initial coefficient of friction plays an important role in determining the mode of subduction initiation. The initial coefficient of friction may not vary much for different rocks. However, in Equation 3 for yielding, we did not consider the effect of pore fluid pressure, which can significantly reduce the yield stress. Therefore, the coefficient of friction in our model represents an effective value that could vary significantly between 0.0 and 0.6 for different locations due to the possible existence of different pore fluid pressures.

Water transport, melt generation, and melt migration in subduction zones can each be complicated (e.g., Iwamori, 1998; Hebert et al., 2009b). Here, we consider a wet zone with uniformly distributed water content above the slab and simplify the treatment of slab-derived trace elements in the melting region. We ignore the consequences of melting for the energy budget, the effect of latent heat on the thermal field, and the feedback to melt generation. Also, we consider that melts are extracted vertically from the mantle to the surface immediately after they are generated. In short, the detailed migration and fractionation processes of melts and water transport in the mantle wedge are ignored. So, our results only reflect a first-order approximation of magma compositions, spatial distribution, and eruption timing during subduction initiation.

In summary, we find that with a mode of continuous initiation with infant-arc spreading, the foundering of the subducting slab, and water release from the slab govern a succession from MORB-like composition to boninites. The modeled composition transition typically occurs within a few million years. When plate strength is reduced, the subducting slab tends to segment, with extensive melting occurring when the slab breaks; most melting occurs close to the trench. When plate strength increases, subduction initiation becomes continuous without infant-arc spreading; such a mode leads to limited, very low degrees of melting occurring during a long interval of plate convergence before subduction initiation starts, although extensive melting near the trench is still possible when subduction initiation starts after a protracted period of plate convergence (∼10 m.y.). If the subduction initiation is driven by constant stresses such as through ridge push, the slab subducts rapidly in response to continuous acceleration of the plate under action of the far-field push; significant melting, including boninite eruption, can be generated within a few million years with no trench migration. Based on the tectonic and volcanic evolution, these different modes may be applicable to the initiation of the Izu-Bonin-Mariana arc, the New Hebrides arc, the Puysegur Trench in New Zealand, and the Aleutian Trench, respectively.

We thank Bob Stern and two anonymous reviewers, whose comments help to improve the manuscript. Leng was supported by the O.K. Earl Fellowship at Caltech. Additional support was provided by the National Science Foundation (EAR-0810303) and the Caltech Tectonics Observatory (by the Gordon and Betty Moore Foundation) with contribution number TO 209.