Plate tectonics is fundamental to Earth's processes, influencing everything from continent formation to volcanism and the planet's habitability. Subduction, where one tectonic plate slides beneath another, is the primary driver of plate tectonics, yet its initiation remains an enigma. While dynamic numerical models and fragmented geological records have been used to investigate this process, a complete understanding remains elusive due to the scarcity of ongoing, whole-plate scale subduction initiation events. This paper leverages the first in situ record of subduction initiation from IODP Expedition 352 at the Izu-Bonin-Mariana (IBM) subduction zone, a region exhibiting a well-preserved magmatic progression. This unique dataset allows for the construction of geologically and geochemically constrained dynamic models to investigate the driving forces behind subduction initiation. The IBM system's magmatic record, characterized by a sequence of fore-arc basalts (FAB) followed by boninites and then normal arc lavas, provides crucial temporal constraints on this process. High-precision dating reveals that FAB magmatism, preceding boninite magmatism, lasted only 0.6–1.2 Myr. The near-simultaneous occurrence of FAB magmatism along the entire IBM arc points to a single, rapid initiation event with significant implications for global plate reorganization at the beginning of the Cenozoic. This study explores the 'vertically driven' subduction initiation model, contrasting it with the 'horizontally forced' model, to determine the dominant forces involved in IBM subduction initiation. The vertically driven model emphasizes the role of internal buoyancy forces, contrasting with the horizontally forced model which relies on external, far-field forces.

Previous research on subduction initiation has relied heavily on dynamic numerical models and fragmented geological evidence. Stern and Gerya (2018) provide a comprehensive review of these efforts, highlighting the challenges in definitively establishing the dominant driving mechanisms. The Izu-Bonin-Mariana (IBM) subduction zone has emerged as a key area of study due to its relatively complete magmatic record (Reagan et al., 2017). The transition from fore-arc basalts (FAB) to boninites is a well-documented feature of this system (Reagan et al., 2010; Hickey & Frey, 1982), and high-precision dating has provided crucial temporal constraints on this transition (Ishizuka et al., 2006, 2011; Reagan et al., 2019). Conceptual models, such as the one proposed by Stern and Bloomer (1992), suggested a vertically driven process involving buoyancy differences and the formation of a lithospheric gap, resulting in decompression melting and FAB formation. However, dynamic testing of this model has been incomplete, and the role of external forces, such as far-field compression, remains debated (Leng & Gurnis, 2011, 2015; Leng et al., 2012; Dymkova & Gerya, 2013; Zhou et al., 2018). The isotopic composition of the initial boninites suggests slab melting contributions (Li et al., 2019), further complicating the picture. Existing studies have explored both spontaneous (Stern, 2004) and induced (Lallemand, 2016) subduction initiation, with varying degrees of success in capturing the observed magmatic sequences and timescales. The debate around the relative importance of internal versus external forces in driving subduction initiation is central to this study.

The authors employed a 2D numerical model of a trench-perpendicular cross-section through the IBM subduction system using the thermomechanical code Fluidity (Davies et al., 2011; Kramer et al., 2012). The model incorporated a 50 Myr old Pacific plate and a 5 Myr old Philippine Sea plate separated by a 10 km wide pre-damaged transform zone. The model included a free surface, generating an internal force akin to 'ridge-push,' which was balanced by an equal and opposite force applied at the boundary. The initial horizontal forces were thus balanced at the start of the model run. To investigate the effect of far-field compression, a larger external force was applied. Two-dimensional models can not capture the effects of along-strike propagation; thus, to mimic out-of-plane forces, a small additional vertical pull force was applied to the Pacific plate. The model simulated various scenarios by adjusting the external compressive force, allowing the observation of different initiation behaviors: no initiation, vertically driven initiation, and horizontally forced initiation. The model tracked temperature, melt generation, and the movement of the plates over time. The rheology of the model materials was described using ductile mechanisms (diffusion, dislocation, and Peierl's creep) and plastic deformation. The effective viscosity was calculated by combining the viscosities from the ductile and plastic mechanisms, incorporating a damage model to account for changes in the rock's properties over time. Mantle melting was calculated using the parameterization of Katz et al. (2003), incorporating hydration of the mantle via dehydration of the subducting slab. All melt was assumed to be fully extracted and travel vertically. The model output tracked the spatial and temporal distribution of magmatic products, providing predictions of FAB and boninite distribution. The model was run until negligible decompression melting and surface velocities in the fore-arc were observed, allowing for the analysis of magmatic product distribution at the time of subduction establishment. The initial conditions, rheological parameters, and melt calculations were based on existing literature and experimental data (Hirth & Kohlstedt, 2003; Jaoul et al., 1980; Houlier et al., 1990; Kirby & Kronenberg, 1987; Karato & Wu, 1993; Katayama & Karato, 2008; Bouilhol et al., 2015; Schmidt & Poli, 2014; Nichols et al., 1994). Specific parameters for the Arrhenius formulation of Peierl's creep were determined via a parameter grid search to minimize the difference with experimentally determined flow laws (Katayama & Karato, 2008).

The model successfully reproduced the observed magmatic sequence of the IBM, both temporally and spatially. When net initial horizontal push was set to zero, only a small additional vertical pull force was needed to initiate subduction and progress to fully down-dip subduction. This process generated a lithospheric gap, resulting in asthenospheric upwelling and decompression melting, forming FAB at a rate of ~4 cm/yr. The timing of the transition from FAB magmatism to boninite magmatism (~0.6 Myr) closely matched the drill core data. As subduction progressed, the sinking plate's tip crossed its solidus, initiating slab melting and boninite formation. The spatial distribution of FAB and boninites in the model also aligned with the observations, with FAB accumulating closer to the trench and boninites further away. The model correctly predicted the initial scraping of sediments, delaying their entry into the subduction zone, consistent with isotopic data. Importantly, the model demonstrated that FAB generation did not occur when far-field horizontal compressive forces dominated. In a suite of models investigating the influence of external compressive force, three types of behavior were observed: no initiation, vertically driven initiation (with FAB formation), and horizontally forced initiation (without FAB). A regime diagram was generated, summarizing the conditions required for each behavior. The key finding was that FAB is uniquely produced in the vertically driven type of initiation; therefore, the presence of FAB might serve as an indicator for this style of subduction initiation where internal forces dominate over external forces. The simulations suggest that vertically driven subduction initiation events are likely large-scale and thus rare, but their impact might extend to influencing wider plate motions. The presence of FAB-like rocks in Tethyan ophiolites further supports the significance of this internally-driven initiation mechanism on a global scale.

The findings strongly support the hypothesis that subduction initiation at the IBM was vertically driven by internal buoyancy forces. The ability of the model to accurately reproduce both the temporal and spatial distribution of magmatic products, particularly the presence and timing of FAB formation, provides compelling evidence for this mechanism. The finding that FAB formation is unique to the vertically driven initiation suggests that the existence of FAB in a geological record could be considered a diagnostic indicator of this process. The model also sheds light on the scale of these events, suggesting that they are necessarily whole-plate scale phenomena, implying rarity but potentially significant impact on global plate motions. The discovery of FAB-like rocks in Tethyan ophiolites suggests that similar self-driven, whole-plate scale events may have played a significant role in the formation of other subduction zones. The model's predictions about sediment melting timing and location further support its validity. The observed close match between model predictions and geological observations, specifically concerning the timing and location of magmatic products, strengthens the support for the vertically driven mechanism. The fact that FAB is not formed in horizontally forced scenarios highlights the importance of internal buoyancy forces in driving this specific type of subduction initiation.

This study provides strong evidence that the Izu-Bonin-Mariana subduction initiation was primarily driven by internal, vertical forces, as opposed to external, horizontal forces. The numerical models successfully replicated the observed magmatic progression (FAB to boninite) both temporally and spatially, highlighting the critical role of internal buoyancy. The unique association of FAB formation with this type of initiation suggests its potential use as a diagnostic indicator. The inferred whole-plate scale of these events explains their relative rarity yet significant implications for broader tectonic reorganizations. Future research should focus on further investigating the interplay between internal and external forces in various geological contexts, refining the model to encompass three-dimensional complexities, and exploring the implications of these findings for understanding other subduction zones and global plate tectonics.

The study utilizes a 2D model, which inherently simplifies the complexities of a three-dimensional system. Along-strike variations and the propagation of subduction initiation along the transform fault are not fully captured. The model also makes certain simplifying assumptions regarding melt extraction and migration, potentially impacting the quantitative aspects of the results. While the model accounts for damage accumulation, the lack of a healing mechanism might affect the model's long-term behavior. Despite these limitations, the model successfully reproduces the key temporal and spatial features of the magmatic sequence, providing robust qualitative insights into the underlying driving forces of subduction initiation.