Rapid subduction initiation and magmatism in the Western Pacific driven by internal vertical forces

The study addresses how subduction initiates, a fundamental yet elusive aspect of plate tectonics. Few modern examples exist and none at whole-plate scale, so understanding relies on numerical models and the fragmented geologic record. The Izu-Bonin-Mariana (IBM) system preserves the most complete magmatic record of subduction initiation, transitioning from fore-arc basalts (FAB) to boninites and then to typical arc magmas. High-precision dating from IODP Expedition 352 shows FAB magmatism preceded boninites by only ~0.6–1.2 Myr and occurred nearly synchronously along the entire IBM (~2400 km), implying a rapid, whole-plate-scale initiation event. The prevailing conceptual model invokes a buoyancy-driven, “vertically driven” initiation where an older dense plate founders adjacent to a younger buoyant plate across a large-offset transform, creating asthenospheric upwelling and decompression melting (FAB), followed by slab-tip dehydration/melting leading to boninites. The paper tests whether internal vertical forces versus external horizontal forcing better explain the observed timing and spatial distribution of magmatism and the progression to self-sustained subduction.

Prior work includes reviews of subduction initiation in nature and models (e.g., Stern & Gerya, 2018) and the IBM magmatic record (Reagan et al., 2010; Ishizuka et al., 2011; 2006). The FAB-to-boninite sequence and near-synchronous ages suggest rapid, system-wide initiation. Concepts of “spontaneous” or internally buoyancy-driven initiation (Stern & Bloomer, 1992; Stern, 2004) contrast with “induced” or externally forced initiation (e.g., Hall et al., 2003; Lallemand, 2016). Three-dimensional models have shown initiation may start at a point (e.g., ridge-transform junction) and propagate or “unzip” along strike (Zhou et al., 2018). Additional hypotheses invoke initiation at relic arcs (Leng & Gurnis, 2015) or plume rejuvenation (Ishizuka et al., 2013). Geochemical constraints, including Hf–Nd isotopes in earliest boninites, indicate slab melt contributions (Li et al., 2019). Global tectonic reorganizations near the Eocene (e.g., Hawaiian–Emperor bend timing) provide potential links or triggers (Whittaker et al., 2007; O’Connor et al., 2013).

The authors construct a 2D trench-perpendicular thermomechanical model of the IBM system using the open-source finite element code Fluidity on an adaptive mesh (<300 m local resolution). The domain is 1000 km wide by 300 km deep. Initial setup: a 50 Myr Pacific plate abuts a 5 Myr Philippine Sea plate across a 10 km wide, plastically weak, hydrated transform/fracture zone (to 40 km depth). The Philippine plate transitions to 50 Myr age in the far backstop (last 100 km). The top 1 km is sediment, underlain by 7 km crust; the remainder is mantle. A free surface enables development of an internally generated horizontal push analogous to ridge-push from the buoyant young plate; this internal force is ~2 TN/m. Boundary conditions: top free surface; bottom closed, no-slip; Philippine-side boundary closed/no-slip to 80 km and open below; Pacific-side boundary open with optional applied pressure above 50 km to impose external horizontal push. Forces: At baseline, an opposing horizontal push equal to 2 TN/m is applied at the boundary to balance internal forces. A small additional vertical pull is applied within a 20–40 km box beneath/near the transform on the Pacific plate by locally increasing density to emulate out-of-plane forces from along-strike unzipping. A suite of models varies vertical pull and external horizontal compression. Rheology: Ductile deformation via diffusion, dislocation, and Peierls creep combined in parallel; flow law parameters chosen to fall within commonly used dry peridotite laws; Peierls behavior approximated by Arrhenius fit to Katayama & Karato (2008) via parameter grid search. Plasticity uses a depth-dependent Byerlee yield law with cohesion and friction coefficient reduced with accumulated plastic strain (“damage”); sediments and fracture zone are initially maximally damaged; no healing is modeled. Effective viscosity combines ductile and plastic components in parallel, with bounds 1e19–1e24 Pa·s. Melting and magmatism: Mantle melting follows Katz et al. (2003) parameterization. Hydration: mantle above dehydrating slab crust is assigned 0.05 wt% bulk H2O (assumes vertical fluid migration only); otherwise dry. Slab crust solidus from Bouilhol et al. (2015), close to Schmidt & Poli (2014). Sediment melting tracked using Nichols et al. (1994) solidus. Decompression melting rate is computed from upwelling and integrated along streamlines tracking depletion. All melt is instantaneously extracted and transported vertically to the surface; erupted products are then advected horizontally using the surface velocity (upwind scheme) to predict spatial accumulation (fore-arc vs nearer future arc/back-arc). Model termination occurs when decompression melting and fore-arc surface velocities become negligible, prior to the slab reaching the bottom boundary. Parameter exploration includes 20 models spanning vertical pull and horizontal compression; a regime diagram summarizes outcomes (no initiation, vertically driven with FAB, horizontally forced without FAB).

- A small additional vertical pull of 9 TN/m (far below a typical slab pull ~30 TN/m) suffices to trigger rapid, vertically driven subduction initiation in the Pacific–Philippine Sea configuration with initial force balance of horizontal pushes. - Vertically driven initiation produces asthenospheric upwelling beneath the former transform, a lithospheric gap, and decompression melting that forms fore-arc basalts (FAB). The average FAB crust formation (eruption) rate is ~4 cm/yr. - The interval between onset of FAB magmatism and onset of boninitic magmatism is ~0.6 Myr, matching drill core constraints of 0.6–1.2 Myr from IODP Exp. 352. - Spatial distribution matches observations: FABs erupt and are advected to accumulate nearest the future trench (fore-arc); subsequent boninites erupt above the slab tip and are advected to accumulate landward, closer to the future arc. Drill sites nearest the trench contain only FAB; sites ~30 km toward the arc contain only boninite. - Early sediments are scraped off and delay entering the subduction zone; isotopic signatures of sediments appear only in younger boninites, consistent with model predictions. Earliest boninites show a slab melt component, consistent with shallow slab crust crossing its solidus as initiation proceeds. - As down-dip subduction establishes, decompression melting migrates away from the trench toward the future arc/back-arc; predicted timing aligns with IODP Exp. 351 back-arc decompression melts west of the Kyushu-Palau Ridge. - The brief pre-boninite FAB production period in a single 2D slice (~0.6 Myr) is consistent with near-synchronous FAB ages along the ~2400 km IBM, implying extremely rapid along-strike propagation. - In contrast, models with predominant external horizontal compression (“horizontally forced”) do not create a lithospheric gap or shallow upwelling; no FAB is produced. In a representative case with +4 TN/m extra push, initiation takes ~5 Myr before the Pacific plate begins to subduct beneath the Philippine plate; once initiated, subsequent kinematics are comparable to the vertically driven case. Such behavior resembles smaller, externally forced systems (e.g., Puysegur, Matthew-Hunter). - Regime analysis shows three outcomes depending on vertical pull and horizontal compression: no initiation; initiation without FAB (horizontally forced); and initiation with FAB (vertically driven). While precise boundaries shift with other parameters (fracture zone width, plate ages, rheology, fluid weakening), the trend that FAB uniquely signifies vertically driven initiation is robust.

The results demonstrate that the IBM magmatic sequence and its rapid timescale are best explained by initiation dominated by internal, vertical buoyancy forces rather than far-field horizontal compression. The model reproduces both the temporal evolution (FAB to boninite within ~0.6–1.2 Myr) and spatial segregation of magmatic products observed in IODP drill cores. The necessity of a lithospheric gap and decompression melting to form FAB explains why FAB is absent in primarily horizontally forced initiation. Therefore, the presence of FAB can serve as an indicator (“smoking gun”) of vertically/buoyancy-driven, whole-plate-scale subduction initiation events. Such events likely require strong plate age/buoyancy contrasts (or old buoyant arcs) and are rare, but when they occur, they may coincide with broader reorganizations of plate motions. The timing of IBM initiation relative to the Hawaiian–Emperor bend suggests a potential link, though causality remains debated. The similarity of fore-arc sequences and FAB-like rocks in Tethyan ophiolites implies that segments of the Tethyan margin may have experienced similar internally driven initiation.

This study integrates geologically and geochemically constrained numerical models with the unique in situ IBM record to show that subduction initiation driven primarily by internal vertical buoyancy forces can rapidly generate the observed FAB-to-boninite progression and evolve into self-sustained subduction. The models match both timing (~0.6–1.2 Myr from FAB onset to boninite) and spatial patterns (fore-arc FAB, landward boninite) and predict the absence of FAB when initiation is governed by external horizontal compression. Consequently, FAB occurrence is proposed as a diagnostic signature of vertically driven initiation, likely associated with whole-plate-scale processes and potentially global plate reorganizations. Future research directions include fully 3D modeling of along-strike nucleation and propagation of initiation; systematic exploration of parameter sensitivities (fracture zone properties, plate ages, rheologies, fluid weakening and healing); improved treatment of fluid and melt transport (non-vertical pathways, permeability evolution, two-phase flow); integration with additional geochronologic and geochemical datasets from IBM and analogous ophiolites; and testing whether similar signatures identify other past initiation events (e.g., along the Tethyan margin).

- Dimensionality: The primary models are 2D and cannot explicitly capture along-strike nucleation and propagation; out-of-plane forces are emulated by an imposed local vertical pull. - Force implementation: External horizontal pushes and localized vertical pulls are parameterized; their magnitudes and spatial extents are simplified proxies for complex 3D dynamics. - Parameter sensitivity: Regime boundaries depend on fracture zone width, plate ages, rheological parameters, and fluid weakening; only a subset was explored. - Rheology and damage: Plastic weakening uses a simplified damage formulation without healing; flow laws are chosen within experimental bounds but remain approximations (e.g., Arrhenius fit to Peierls creep). - Melt and fluid transport: Assumes instantaneous melt extraction and vertical ascent; mantle hydration is prescribed at 0.05 wt% H2O above dehydrating slab crust with vertical fluid migration only; these simplifications neglect lateral transport and two-phase coupling. - Thermal and compositional structure: Initial conditions (layer thicknesses, ages, hydration of fracture zone) and solidus choices influence results; alternative plausible structures may shift thresholds. - Boundary effects: Although simulations stop before interaction with the bottom boundary, finite domain size and boundary conditions may still influence late-stage evolution. - Resolution and model termination: While mesh is highly refined, sub-kilometer processes and long-term arc evolution (post-initiation) are not resolved or modeled.