Dynamics of subduction initiation with different evolutionary pathways

Changes of plate motion may have induced subduction initiation (SI), but the tectonic history of SI is different from one subduction zone to another. Izu-Bonin-Mariana (IBM) SI, accompanied by strong backarc spreading and voluminous eruption of Boninites, contrasts with the Aleutians which shows neither. Using finite element models, we explore visco-elasto-plastic parameters and driving boundary conditions for SI evolution. With an imposed velocity, we find three different evolutionary modes of SI: continuous without backarc spreading, continuous with backarc spreading and a segmented mode. With an increase in the coefficient of friction and a decrease in the rate of plastic weakening, the amount of convergence needed for SI increases from ∼20 to ∼220 km, while the mode changes from segmented to continuous with backarc spreading and eventually to continuous without backarc spreading. If the imposed velocity boundary condition is replaced with an imposed stress, the amount of convergence needed for SI is reduced and backarc spreading does not occur. These geodynamic models provide a basis for understanding the divergent geological pathways of SI. First, IBM evolution is consistent with subduction of an old strong plate with an imposed velocity which founders causing intense backarc spreading and Boninitic volcanism. Second, the New Hebrides SI is in the segmented mode due to its weak plate strength. Third, the Puysegur SI is in the continuous without backarc spreading mode with no associated volcanic activities. Fourth, the Aleutians SI has neither trench rollback nor backarc spreading because the slab is regulated by constant ridge-push forces

Simultaneous inversion of mantle properties and initial conditions using an adjoint of mantle convection

Through the assimilation of present-day mantle seismic structure, adjoint methods can be used to constrain the structure of the mantle at earlier times, i.e., mantle initial conditions. However, the application to geophysical problems is restricted through both the high computational expense from repeated iteration between forward and adjoint models and the need to know mantle properties (such as viscosity and the absolute magnitude of temperature or density) a priori. We propose that an optimal first guess to the initial condition can be obtained through a simple backward integration (SBI) of the governing equations, thus lessening the computational expense. Given a model with known mantle properties, we show that a solution based on an SBI-generated first guess has smaller residuals than arbitrary guesses. Mantle viscosity and the effective Rayleigh number are crucial for mantle convection models, neither of which is exactly known. We place additional constraints on these basic mantle properties when the convection-induced dynamic topography on Earth's surface is considered within an adjoint inverse method. Besides assimilating present-day seismic structure as a constraint, we use dynamic topography and its rate of change in an inverse method that allows simultaneous inversion of the absolute upper and lower mantle viscosities, scaling between seismic velocity and thermal anomalies, and initial condition. The theory is derived from the governing equations of mantle convection and validated by synthetic experiments for both one-layer viscosity and two-layer viscosity regionally bounded spherical shells. For the one-layer model, at any instant of time, the magnitude of dynamic topography is controlled by the temperature scaling while the rate of change of topography is controlled by the absolute value of viscosity. For the two-layer case, the rate of change of topography constrains upper mantle viscosity while the magnitude of dynamic topography determines the temperature scaling (lower mantle viscosity) when upper-mantle (lower-mantle) density anomaly dominates the flow field; this two-stage scheme minimizes the tradeoff between temperature and lower mantle viscosity. For both cases, we show that the theory can constrain mantle properties with errors arising through the adjoint recovery of the initial condition; for the two-layer model, this error is manifest as a tradeoff between the temperature scaling and lower mantle viscosity

Evolving force balance during incipient subduction

Nearly half of all active subduction zones initiated during the Cenozoic. All subduction zones associated with active back arc extension have initiated since the Eocene, hinting that back arc extension may be intimately associated with an interval (several tens of Myr) following subduction initiation. That such a large proportion of subduction zones are young indicates that subduction initiation is a continuous process in which the net resisting force associated with forming a new subduction zone can be overcome during the normal evolution of plates. Subduction initiation is known to have occurred in a variety of tectonic settings: old fracture zones, transform faults, and extinct spreading centers and through polarity reversal behind active subduction zones. Although occurring within different tectonic settings, four known subduction initiation events (Izu-Bonin-Mariana (IBM) along a fracture zone, Tonga-Kermadec along an extinct subduction boundary, New Hebrides within a back arc, and Puysegur-Fiordland along a spreading center) were typified by rapid uplift within the forearc followed by sudden subsidence. Other constraints corroborate the compressive nature of IBM and Tonga-Kermadec during initiation. Using an explicit finite element method within a two-dimensional domain, we explore the evolving force balance during initiation in which elastic flexure, viscous flow, plastic failure, and heat transport are all considered. In order to tie theory with observation, known tectonic settings of subduction initiation are used as initial and boundary conditions. We systematically explore incipient compression of a homogeneous plate, a former spreading center, and a fracture zone. The force balance is typified by a rapid growth in resisting force as the plate begins bending, reaching a maximum value dependent on plate thickness, but typically ranging from 2 to 3 × 1012 N/m for cases that become self-sustaining. This is followed by a drop in stress once a shear zone extends through the plate. The formation of a throughgoing fault is associated with rapid uplift on the hanging wall and subsidence on the footwall. Cumulative convergence, not the rate of convergence, is the dominant control on the force balance. Viscous tractions influence the force balance only if the viscosity of the asthenosphere is >1020 Pa s, and then only after plate failure. Following plate failure, buoyancy of the oceanic crust leads to a linear increase with crustal thickness in the work required to initiate subduction. The total work done is also influenced by the rate of lithospheric failure. A self-sustaining subduction zone does not form from a homogeneous plate. A ridge placed under compression localizes subduction initiation, but the resisting ridge push force is not nearly as large as the force required to bend the subducting plate. The large initial bending resistance can be entirely eliminated in ridge models, explaining the propensity for new subduction zones to form through polarity reversals. A fracture zone (FZ) placed in compression leads to subduction initiation with rapid extension of the overriding plate. A FZ must be underthrust by the older plate for ~100–150 km before a transition from forced to self-sustaining states is reached. In FZ models the change in force during transition is reflected by a shift from forearc uplift to subsidence. Subduction initiation is followed by trench retreat and back arc extension. Moderate resisting forces associated with modeled subduction initiation are consistent with the observed youth of Pacific subduction zones. The models provide an explanation for the compressive state of western Pacific margins before and during subduction initiation, including IBM and Tonga-Kermadec in the Eocene, and the association of active back arcs with young subduction zones. On the basis of our dynamic models and the relative poles of rotation between Pacific and Australia during the Eocene, we predict that the northern segment of the Tonga-Kermadec convergent margin would have initiated earlier with a progressive southern migration of the transition between forced and self-sustaining states

Slabs in the lower mantle and their modulation of plume formation

Numerical mantle convection models indicate that subducting slabs can reach the core-mantle boundary (CMB) for a wide range of assumed material properties and plate tectonic histories. An increase in lower mantle viscosity, a phase transition at 660 km depth, depth-dependent thermal expansivity, and depth-dependent thermal diffusivity do not preclude model slabs from reaching the CMB. We find that ancient slabs could be associated with lateral temperature anomalies ~500°C cooler than ambient mantle. Plausible increases of thermal conductivity with depth will not cause slabs to diffuse away. Regional spherical models with actual plate evolutionary models show that slabs are unlikely to be continuous from the upper mantle to the CMB, even for radially simple mantle structures. The observation from tomography showing only a few continuous slab-like features from the surface to the CMB may be a result of complex plate kinematics, not mantle layering. There are important consequences of deeply penetrating slabs. Our models show that plumes preferentially develop on the edge of slabs. In areas on the CMB free of slabs, plume formation and eruption are expected to be frequent while the basal thermal boundary layer would be thin. However, in areas beneath slabs, the basal thermal boundary layer would be thicker and plume formation infrequent. Beneath slabs, a substantial amount of hot mantle can be trapped over long periods of time, leading to “mega-plume” formation. We predict that patches of low seismic velocity may be found beneath large-scale high seismic velocity structures at the core-mantle boundary. We find that the location, buoyancy, and geochemistry of mega-plumes will differ from those plumes forming at the edge of slabs. Various geophysical and geochemical implications of this finding are discussed

Adjoint models of mantle convection with seismic, plate motion, and stratigraphic constraints: North America since the Late Cretaceous

We apply adjoint models of mantle convection to North America since the Late Cretaceous. The present-day mantle structure is constrained by seismic tomography and the time-dependent evolution by plate motions and stratigraphic data (paleoshorelines, borehole tectonic subsidence, and sediment isopachs). We infer values of average upper and lower mantle viscosities, provide a synthesis of North American vertical motions (relative sea level) from the Late Cretaceous to the present, and reconstruct the geometry of the Farallon slab back to the Late Cretaceous. In order to fit Late Cretaceous marine inundation and borehole subsidence, the adjoint model requires a viscosity ratio across 660 km discontinuity of 15:1 (reference viscosity of 10^(21) Pa s), which is consistent with values previously inferred by postglacial rebound studies. The dynamic topography associated with subduction of the Farallon slab is localized in western North America over Late Cretaceous, representing the primary factor controlling the widespread flooding. The east coast of the United States is not stable; rather, it has been experiencing continuous dynamic subsidence over the Cenozoic, coincident with an overall eustatic fall, explaining a discrepancy between sea level derived from the New Jersey coastal plain and global curves. The east coast subsidence further constrains the mantle viscosity structure and requires an uppermost mantle viscosity of 10^(20) Pa s. Imposed constraints require that the Farallon slab was flat lying during Late Cretaceous, with an extensive zone of shallow dipping Farallon subduction extending beyond the flat-lying slab farther east and north by up to 1000 km than previously suggested

The effects of chemical density differences on convective mixing in the Earth's mantle

Subducted lithosphere and crust are chemically distinct from ambient mantle, and some have hypothesized that the resulting anomalous buoyancy may affect mixing and, in particular, may cause chemical heterogeneity to persist. To explore this effect, two-dimensional convective mixing calculations have been carried out. Boundary conditions mimicking features of plate kinematics have been imposed and negatively buoyant tracers placed beneath converging plate margins to simulate subduction. The principal parameter of interest is the density anomaly ratio: Δρ/(ρ_o α ΔT), the ratio of the chemical to the thermal buoyancy. For a density anomaly ratio appropriate to complete separation of the oceanic crust from the rest of the lithosphere, it is marginally possible for the tracers to segregate into the warm region near the bottom and neutralize the thermal buoyancy there, for internally heated flows, and to segregate to the base of the uprising limbs of convection, for bottom heated flows. However, for the more likely scenario of oceanic crust and lithosphere maintaining their layering after subduction, there is neither segregation nor neutralization of thermal buoyancy. The average residence time is the average time tracers remain in the box from subduction to sampling beneath ridges and could be increased by 4% from the passive case at the very most if the crust and lithosphere remain layered. Even in the more extreme case of separation of crust from lithosphere, the residence times are only increased by 15% from the passive case

On the location of plumes and lateral movement of thermochemical structures with high bulk modulus in the 3-D compressible mantle

The two large low shear velocity provinces (LLSVPs) at the base of the lower mantle are prominent features in all shear wave tomography models. Various lines of evidence suggest that the LLSVPs are thermochemical and are stable on the order of hundreds of million years. Hot spots and large igneous province eruption sites tend to cluster around the edges of LLSVPs. With 3-D global spherical dynamic models, we investigate the location of plumes and lateral movement of chemical structures, which are composed of dense, high bulk modulus material. With reasonable values of bulk modulus and density anomalies, we find that the anomalous material forms dome-like structures with steep edges, which can survive for billions of years before being entrained. We find that more plumes occur near the edges, rather than on top, of the chemical domes. Moreover, plumes near the edges of domes have higher temperatures than those atop the domes. We find that the location of the downwelling region (subduction) controls the direction and speed of the lateral movement of domes. Domes tend to move away from subduction zones. The domes could remain relatively stationary when distant from subduction but would migrate rapidly when a new subduction zone initiates above. Generally, we find that a segment of a dome edge can be stationary for 200 million years, while other segments have rapid lateral movement. In the presence of time-dependent subduction, the computations suggest that maintaining the lateral fixity of the LLSVPs at the core-mantle boundary for longer than hundreds of million years is a challenge

Subduction initiation at relic arcs

Although plate tectonics is well established, how a new subduction zone initiates remains controversial. Based on plate reconstruction and recent ocean drilling within the Izu-Bonin-Mariana, we advance a new geodynamic model of subduction initiation (SI). We argue that the close juxtaposition of the nascent plate boundary with relic oceanic arcs is a key factor localizing initiation of this new subduction zone. The combination of thermal and compositional density contrasts between the overriding relic arc, and the adjacent old Pacific oceanic plate promoted spontaneous SI. We suggest that thermal rejuvenation of the overriding plate just before 50 Ma caused a reduction in overriding plate strength and an increase in the age contrast (hence buoyancy) between the two plates, leading to SI. The computational models map out a framework in which rejuvenated relic arcs are a favorable tectonic environment for promoting subduction initiation, while transform faults and passive margins are not