Numerical Modeling of Melting Processes During Slab Break-off: Insights Into Tectonic Setting for Massif-Type Anorthosites

The concept that lithosphere detachment or break-off has long been conceived as a viable mechanism to explain prominent geological phenomena in Earth’s crust and the surface. One of the strengths of slab delamination mechanism is that it can account for the extensive magmatism in active orogenic belts due to the upwelling of the asthenosphere after the slab break-off. However, in the last 20 years, geodynamic simulations show that the inflow of the asthenosphere upon slab break-off is insufficient to cause significant melting of the overriding lithosphere adjacent to the slab window. The primary reasons include the occurrence of slab break-off at a location that is too deep to effectively heat the overriding lithospheric mantle. Another factor is the presence of a thin film of crustal material that is retained during the slab break- off, inhibiting a significant thermal perturbation within the lithosphere. In this work, we couple petrological–thermomechanical simulations with magmatic melting processes to examine the lithospheric melting and surface lithological expression associated with slab break-off. Our work shows that in the early Earth when the mantle temperature is relatively higher, shallow slab break-off can give rise to significant lithospheric melting during the development of slab break-off. Moreover, because the slab becomes weaker in the earlier hotter mantle, it may break-off prior to the stage of continental collision, thus the magmatism it induced may not give a direct constraint on the time of continental collision. Our study has implications for the interpretation of geological and tomography studies in orogenic belts. It also provides insights into reconciling conflicts between geodynamic and geological studies regarding slab break off-induced melting and magmatism.One of the most peculiar lithologies in Earth’s middle age is Proterozoic massif-type anorthosites (PMAs), a plutonic batholith-forming rock type temporally restricted to the Proterozoic [1-4]. Formally, PMAs are composed of at least 90% plagioclase feldspar accompanied by minor mafic silicates and Fe-Ti oxides [5]. PMAs are areally and volumetrically extensive, with the largest PMA being the Kunene Complex in SW Angola, which covers an area of 18,000 km2 [6]. A consensus has been largely reached on the mechanism by which anorthosites were concentrated. They formed through the accumulation of magmatic plagioclase at the top of a magma chamber due to the low density of plagioclase compared to coexisting melt [7]. However, despite their simple mineralogy and have been studied for over a century, the geodynamic setting accounting for PMAs remains hotly debated [1, 2, 8, 9].A variety of tectonic settings have been proposed for PMAs, including Andean-type continental arc, post-orogenic, anorogenic, and continental rift settings [2, 4, 10, 11]. In recent years, a new tectonic regime—slab break-off—has been adopted accounting for the origin of several PMAs in Asia based on detailed geological, geochronological, and geochemical analyses [12-14]. Their work attributes the formation of PMAs to slab break-off–related thermal perturbation. This thermal perturbation leads to the melting of the metasomatized overriding mantle lithosphere, generating basaltic magmatism that gives rise to the mafic-type PMAs at the surface. The slab break-off model, widely adopted for post-collision magmatism [15], is supported by improved seismic tomography, geological observation, and tectonic evolution models [16, 17]. Although the geodynamic processes of slab break-off have been well examined [18-27], whether or not slab break-off can generate significant lithospheric melting is still controversial. Earlier studies by Davies and von Blanckenburg (1995) using analytical models show that the ascent of asthenospheric material during break-off events can induce a thermal disturbance. This disturbance can lead to the melting of the overlying metasomatised mantle lithosphere, but it requires break-off to occur at depths shallower than the overriding plate [28]. Recent numerical work utilized integrated petrological–thermomechanical models to examine the dynamics of break-off and post-collisional magmatism [18, 25, 29], however, their simulations have presented findings inconsistent with the aforementioned break-off–induced melting scenario [25, 30]. Their results reveal that most slab break-off events take place at a deep depth that is insufficient to induce substantial thermal perturbation to melt the overriding lithosphere. Moreover, it is found that the upwelling of asthenospheric mantle does not sustain long enough to elevate the temperature of the metasomatised overlying lithosphere to cause melting [25]. Thus, how to reconcile geodynamic studies with existing geological and geochemical observations that support slab break-off–induced magmatism [13, 14, 17, 28, 31] remains a significant subject warranting further investigation.Because PMAs formed in Proterozoic, a time the mantle potential temperature could be 150–250°C higher than today [32]. A hotter asthenosphere may weaken the subducted slab, which may lead to earlier and shallower slab break-off compared to the present asthenosphere [33]. Furthermore, an increased temperature in the asthenospheric mantle can not only increase the scale but also accelerate the rate of melting in the adjacent overriding lithosphere upon slab break-off. In this work, we coupled thermomechanical numerical models with melting processes to examine the slab break-off–associated melting process. We also utilized thermodynamic melting algorithm to model the magmatic chemical differentiation following the lithospheric melting. This investigation aimed to determine whether it can result in the crystallization of abundant plagioclase and reproduce the compositional range of plagioclases observed in PMAs. Our coupled petrological–thermomechanical simulations with magmatic melting process support slab break-off model as a potential tectonic setting accounting for the PMAs indicated by earlier geochemical studies.The three conservation equations of mass, momentum, and energy are solved under the extended Boussinesq approximation using finite element method with the Underworld package in a 2D cartesian geometry [34-36]. Distinct materials and properties are tracked as Lagrangian particles within a Eulerian finite element grid using particle-in-cell method.where u is velocity, σ′ is deviatoric stress, P is pressure, ρ is density, g is gravitational acceleration, z^ is the vertical unit vector, cp is heat capacity, k is thermal conductivity, and H is heat production.In our simulations, we incorporated additional heating terms in the above equations. For instance, radiogenic heating (⁠Qradiogenic⁠) and partial melting (⁠Qpartial melt⁠) is added into the energy equations as:where A is the radiogenic heat production rate, Lf is the latent heat of fusion which is the amount of energy consumed when a phase changes from solid to liquid, and Mf is the melt fraction. Therefore, in our simulations, the rate of temperature variation is a function of heat diffusion, heat advection, radiogenic heat production, and heat changes due to partial melting, expressed together as:where α is the thermal diffusivity.The density of a material in our experiments thus depends on composition, temperature, and the melt fraction:where ρr , Tr is reference density and reference temperature, respectively. MΔρr represents the density change during melting process. The intrinsic density of the continental crust is 600 kg/m3 smaller than the asthenospheric mantle (3300 kg/m3), while the continental lithospheric mantle has the same density compared to the asthenospheric mantle. Following previous work [37], we only consider thermal-related density variations for oceanic lithosphere, thus the oceanic lithosphere and crust have intrinsic densities the same as the asthenospheric mantle.The melt fraction (⁠Mf⁠) is dynamically calculated during our simulations, using the super-solidus dimensionless temperature given by [38]:where Tss is the super-solidus, and Ts and Tl are the solidus and liquidus, respectively. In this study, our primary focus is on lithospheric partial melting during slab break-off. Unless specified otherwise, we only consider melting within the lithosphere in our calculations to prevent potential overlap with melting formed during subduction process. The solidus and liquidus for the lithosphere mantle follow a polynomial relationship between temperature and pressure as:where P is pressure, and as , bs , cs , al , bl , and cl are defined in Table 1.We use a nonlinear, temperature-dependent, and strain rate-dependent viscoplastic rheology to model the subduction, continental collision, and followed slab break-off process. The viscous deformation is represented by the power law equation with dislocation and diffusion creep in the form as:where η is the viscosity, A is the preexponential factor, n is the stress exponent, d is the grain size, p is the grain-size exponent, f is water fugacity, r is the water fugacity exponent, E is the activation energy, P is the lithostatic pressure, V is the activation volume, R is the gas constant, T is the temperature, and ε˙II is the second invariant of strain rate tensor. Values of these variables are listed in Table 1.The Drucker–Prager pressure-dependent criteria is used for plastic deformation:where τy is the yield stress, C is the cohesion at the surface, ϕ is the internal friction angle, and sin(⁠ϕ⁠) is the friction coefficient. Similar to previous studies [39, 40], C and sin(⁠ϕ⁠) linearly decrease with accumulation of plastic strain, and afterward τy remains constant. To simulate plasticity in areas of large stresses, we impose an upper limit on the stress for lithosphere mantle to mimic the Peierls mechanisms (low‐temperature plasticity) [41]. The effective viscosity η for each material is then determined through the minimum viscosity of each deformation mechanism:In our simulations, the rheological effect related to melting process is determined by an effective viscosity reduction by 2 orders of magnitude when melt fraction is equal or larger than 0.2. It is important to note that there is no melt segregation from the host rocks in our experiments, but we perform thermodynamic calculations to explore its subsequent magmatic process. The viscosity range is limited between 1019 and 1024 Pa s in this study.The numerical experiments were performed in a two-dimensional domain with a width of 800 km and depth of 400 km, and some cases have a dimension of 1200×600 km. Most models use a uniform grid resolution of 513 × 257 nodes, which gives a cell width of ~1.5 km. Each cell contains 30 Lagrangian particles, and more than 3 million particles were used to track material properties in each experiment. An oceanic plate with a width of 250 km was initially placed between two 225 km-long continental plates. All plates are composed of an upper and lower crust and lithospheric mantle (Figure 1). To simulate the topography evolution, a free surface is created by using a 30-km thick “sticky air” layer with a minimum viscosity (1019 Pa s) and low density (1 kg/m3) [42]. Following previous numerical work on slab break-off [20], the initial setup also includes a weak zone near the transition from the oceanic plate and the right continental plate to initiate subduction. The weak zone cuts through the whole lithospheric mantle with a dipping angle of 50°, and is characterized by a minimum viscosity. We fix the dipping angle for the weak zone, even though it may influence the dynamics during the development of slab break-off. However, previous studies suggest that the slab dip angle plays a minor role compared to viscous necking during slab break-off [21]. The viscosity of oceanic crust, when subducted below a depth of 50 km, is reduced by a factor of 103 to mimic the viscosity reduction caused by sediment dehydration of subducting slab [43] and deformation along the mega-thrust faults [44].The boundary condition is a free-slip condition on the top and right boundaries. We apply a constant, convergent velocity (1–10 cm/year) to the left wall in each experiment until the entire 250-km oceanic plate is subducted into the mantle. Afterward, subduction and continental collision are solely driven by slab pull. The bottom boundary is designed as a permeable boundary allowing for inflow and outflow of materials as used in previous studies [45].A constant temperature of 0°C and 1300°C (1450°C) were used for the top and bottom boundary, respectively. The thermal structure of the oceanic lithosphere follows that from the half-space cooling model for a prescribed plate age [46]. For the initial internal temperature of continental plates, it follows a geothermal gradient of 25°C per kilometer for the first 10 km, and then decreases to 19°C (21.8°C) per kilometer until it reaches a temperature of 1300°C (1450°C) at the lithosphere-asthenosphere boundary.To model the magmatic evolution after lithospheric melting, we carry out thermodynamic calculations using alphaMELTS software, version 1.9 [47], a menu-driven interface of MELTS [48, 49], pMELTS/pHMELTs [50]. The entire thermodynamic modeling process consists of two components. The first part involves modeling the isentropic melting of lithospheric mantle during the decompression melting process, followed by the cooling of the magma chamber in an isobaric crystallization fractionation process. We examine the magmatic system under two conditions: one is anhydrous (using pMELTs) and the other with water using pHMELTs, which can calculate water solubility in nominally anhydrous minerals. Calibrated on a large range of terrestrial bulk compositions, pMELTS is ideally suited for modeling magmatic processes like partial melting and crystallization fractionation processes at pressures between 1 and 4 GPa [50]. Major elements of the primitive mantle [51] are used as the initial bulk mantle composition for decompression melting process from 3 to 1GPa, with 0.01 GPa pressure steps and assuming continuous melting with a residual porosity of 0.2%. The depth of the first decompression melting is constrained based on the mantle potential temperature in our thermochemical convection models (1300°C or 1450°C) and the nominally volatile-free peridotite solidus using pMELTS calculations. We stop the decompression melting at the Moho depth (1 GPa) and retrieve the major and trace elements of aggregated liquids at that depth to model the crystal fractional process under isobaric conditions and fayalite—magnetite—quartz oxygen buffer to explore whether it will reproduce the compositional range of plagioclases observed in PMAs.Here, we restrict ourselves to shallow slab break-off scenarios as a minimum requirement for slab break-off–induced magmatism, because if shallow slab break-off does not show any lithospheric melting, we do not need to consider it under intermediate or deep slab break-off. To that end, we start exploring conditions that facilitate a shallow slab break-off by changing plate velocity (1–10 cm/year), plate yield stress (60, 200 MPa), length, and age of the oceanic lithosphere. During the development of slab break-off, we particularly focus on the melting process of the overriding lithosphere. As will be shown below, under present mantle temperature even shallow slab break-off will not be able to generate melting of the lithospheric mantle adjacent to the subducting slab. Next, we conduct a shallow break-off model under warmer mantle temperatures in the early Earth, which clearly shows that upwelling of asthenosphere during slab break-off can cause strong thermal perturbations that melt the overriding mantle lithosphere. Last, because mantle convection code is not designed to model melt segregation and chemical differentiating processes, we perform pMELTs/pHMELTs thermodynamic modeling to model the subsequent magmatic process.Previous numerical investigations indicate faster plate velocity and young oceanic plate promote shallow slab break-off [20]. Here we describe a shallow break-off reference model using this particular setting in detail (Figure 2). Due to the prescribed high plate velocity (10 cm/year) and weak zone between the oceanic and continent plate, subduction of the 20 Ma-old oceanic lithosphere is initiated shortly (Figure 2(a) and (e)). After 2.5 Myr, the entire oceanic plate undergoes subduction, leading to a collision between the two continents (Figure 2(b) and (f)). Right after the continental collision and before the continental crust begins to subduct, the subducting slab rapidly narrows at a depth of 50 km and there is a clear asthenospheric inflow developed under the subduction plate (Figure 2(c) and (g)). The rupture then grows larger and causes the detachment of the subducting slab along the subduction plane (Figure 2(d) and (h)). In line with prior research [25], we find the persistence of a thin crustal film during the viscous creep process subsequent to slab break-off. Considering water-dependent rheology did not substantially change the results. In summary, despite the reference model exhibiting a shallow break-off mode with robust mantle inflow through the slab window, no lithospheric mantle melting associated with slab break-off is observed.In another model, we increase the yield stress of the lithospheric mantle from 80 to 150 MPa, closer to higher estimates [52], we observe a similar shallow break-off dynamics as in the reference model (Figure 3). Due to a higher maximum yield stress, it takes longer for the occurrence of slab break-off, but it still exhibits strong asthenospheric upwelling during slab break-off process after continental collision (Figure 3(b) and (e)). Nevertheless, once more, no lithospheric mantle melting is observed attributed to the slab break-off processes.As shown above under present-day mantle temperature, we do not observe any lithospheric mantle melting associated with slab break-off. However, in the early Earth when PMAs were formed the mantle potential temperature could be 150–250°C higher than present-day value [32], we thus model the shallow slab break-off under a hotter mantle potential temperature. Here, model 3 is characterized with the same setup with respect to the reference case except its mantle temperature is 150°C higher and we correspondingly increase its crustal geothermal gradient at the beginning. Similar to the reference case, subduction of oceanic slab occurs shortly after the model start (Figure 4(a)). However, due to a hotter mantle temperature, the slab becomes weaker and quickly starts to get narrowed (Figure 4(b)). As the subducting slab continues to narrow, a strong inflow of asthenospheric mantle begins to replace the narrowed slab region under the overriding continental plate. Meanwhile, because of the hotter asthenosphere and faster heating rate, the asthenospheric upwelling indeed causes partial melting of the overriding lithospheric mantle (Figure 4(b) and inserted). As the slab necking continues, it leads to more significant melting when the asthenospheric mantle upwells (Figure 4(c) and inserted), and the whole melting process can last for nearly 0.2 Myr (Figure 4(d) and inserted). To further illustrate this lithospheric melting is directly caused by the asthenospheric upwelling during slab break-off, rather than a basal melting, we slightly lower the solidus of the overriding lithosphere by 20°C (1100°C). As shown in Figure 5, significant melting of the overriding lithosphere forms exactly along the curved region due to asthenospheric upwelling, and no melting is observed at other locations of the overriding lithosphere.Above, we illustrate that warmer mantle temperature allows lithospheric melting caused by asthenospheric inflow during the slab break-off. Next, we parameterize the magmatic melting process using alphaMELTS algorithm [47] following the lithospheric melting (Figure 6(a)). We aim to explore whether the simple magmatic evolution pathways can result in the crystallization of abundant plagioclases and if it could have formed the compositional range of plagioclases observed in PMAs. We first model isentropic decompression melting process of the lithospheric mantle without water from 3 to 1 GPa. We find that garnet is completely exhausted in the restite as the fraction of the melted mantle increases, followed by clinopyroxene close to 1 GPa (online Supplementary Table S1). Meanwhile, the extracted melt composition becomes more evolved with higher SiO2 and lower MgO (online Supplementary Table S2). When the decompression melting ceases at 1 GPa, the equilibrated liquid melt is basaltic, which is then used to model isobaric fractional crystallization at the Moho pressure. Our fractional crystallization calculations begin at an initial superliquidus temperature (1500°C), and run until the model hits subsolidus conditions (800°C). Output data of phase proportions and compositions are calculated using a temperature interval of 10°C (online Supplementary TableS3).Figure 6(b) and (c) illustrates the fractional crystallization results, which show olivine first appears, followed by clinopyroxene and spinel, and when temperature drops to 1220°C plagioclase starts to crystallize which is consistent with previous studies [13, 53]. The fractionation trendlines also show that after saturation of plagioclase it remains as the dominant fractionated phase until about 1100°C. Moreover, the composition of the plagioclase is characterized by anorthite content within 50%–30% (Figure 6(c)) (online Supplementary Table S3), consistent with plagioclases in typical PMAs worldwide [1]. We repeat our thermodynamic calculations but for hydrous primitive mantle, and the results remain similar (Figure 7) (online SupplementaryTables S5 and S6). Thus, the results of our simple magmatic evolution following lithosphere melting suggest it can roughly reproduce the plagioclases and its composition known in PMAs. Due to the buoyancy instability of plagioclase among coexisting minerals/melt at the Moho depth, plagioclase-rich mushes have the propensity to undergo descent and form anorthosite diapirs, which will subsequently emplace along major fault zones as PMAs [54].Since the publication of two seminar articles on slab break-off model almost three decades ago [28, 55], it has increasingly gained popularity as an explanation for post-collisional magmatism in many locations where a sudden shift in compositional changes occurs in associated magmatic rocks [56]. But geodynamic models typically reveal slab break-off at a depth that is too deep to cause significant lithospheric melting due to the asthenospheric mantle inflow through the slab window [25]. Another factor is that the upwelling mantle is not sustained enough to sufficiently heat the lithospheric mantle to melt. Moreover, numerical studies consistently show that slab break-off results from viscous creep, leading to the preservation of a thin layer of crustal material (less than 10 km in thickness) beneath the overriding lithosphere [21, 25]. This thin layer acts as an additional shield, limiting the extent of thermal perturbation experienced by the lithospheric mantle. Our results are consistent with previous observations that slab break-off–induced melting is difficult in the modern Earth (Figures 2 and 3). However, due to the secular cooling, the mantle temperature in early Earth could be 300°C higher than today [57]. Here, our new thermomechanical simulations, incorporating melting processes, reveal that during the early stages of Earth’s evolution, elevated mantle temperatures are conducive to the occurrence of shallow slab break-off (Figure 4). Additionally, the increased mantle heat facilitates the melting of the adjacent lithospheric mantle during the necking process, preceding the actual slab detachment (Figure 5), although slab rollback origin of the mantle upwelling cannot be completely ruled out [30]. Nevertheless, the higher mantle temperature in the early Earth extends both the extent of melting and duration of the lithospheric melting associated with slab break-off. Despite that the melting extent is still relatively limited compared to the proposed extent inferred from many geological and geochemical studies [58], it is important to note that post-collisional magmatism is typically a smaller phenomenon compared to magmatism associated with subduction [59]. Furthermore, we do not take into account that hydration (metasomatism) of the overlying lithosphere resulting from oceanic subduction can significantly lower the solidus [60], thereby increasing the extent of melting in the overriding lithosphere due to the thermal perturbation induced by slab break-off.While the present work demonstrates increased lithospheric mantle melting due to slab break-off during the early Earth, it’s crucial to emphasize a distinct aspect of this break-off process. It is observed that these slab break-off events may take place before the actual collision of continents (Figure 4(g) and (h)). As a result, surface magmatism timing records may not necessarily indicate that continental collision occurred before the magmatism induced by lithospheric mantle melting during slab break-off. Because multiple slab break-off events and associated magmatism can occur in the same subduction zone. This may help to explain the cause of multiple exhumation pulses as suggested by previous studies [15]. Similarly, as mantle temperature gradually cools, slabs get stronger and multiple slab break offs should be less common compared to the early Earth. Our work also agrees with a preceding investigation by van Hunen and van den Berg [33] which shows that hotter mantle temperature will lower viscosity leading to more frequent slab break offs [33]. Thus, it may also explain why the post-collisional mafic magmas are ubiquitous since the late Archean [61]. However, due to the possibility of multiple slab break-off events in the same subduction zone, future studies should exercise caution when attempting to interpret the timing of continental collisions solely based on post-collisional magmatic activities.The partial melting and petrological expression during slab break-off are inherently three-dimensional and contain a lot of complexities, which can vary the results and thermochemical evolution [19, 62]. Additionally, it’s important to note that we have not incorporated models for melt segregation and emplacement in our analysis, primarily because coupling multiphase flow in large-scale mantle convection models remains computationally expensive. Due to these factors, it remains challenging to quantifiably constrain how the interconnected melts create their pathways and ascend to the surface, or what controls the surface distributions of PMAs. Nevertheless, the gravity-driven diapiric uprise of anorthosite plutons is favored [63], especially in areas where crustal-scale zones of weakness and terrane boundaries extend to Moho offsets [54]. These aspects hold particular significance in the context of comprehending the formation, migration, and emplacement of anorthositic magma within the shallow crustal settings. Nonetheless, our objective is to furnish a preliminary, comprehensive knowledge of how lithospheric melting may occur during the process of slab break-off and exert an influence on surface magmatic events.We perform coupled petrological–thermomechanical numerical models to explore lithospheric melting and magmatic evolution processes associated with slab break-off. Our simulations show that under present-day mantle temperature conditions, slab break-off is unlikely to generate lithospheric melting. However, within the context of an early Earth characterized by increased mantle temperatures, our findings elucidate that slab break-off possesses the capability to induce substantial lithospheric mantle melting, consequently giving rise to surface magmatism. Nevertheless, it is imperative to exercise caution and employ additional geological and geochemical examination when assessing whether these magmatic events are formed within post-collision scenarios or not as slab break-off can occur prior to continent collision.The data for this study are available in this manuscript and supplementary material. The mantle convection code Underworld2 version 2.10.1b is available at https://zenodo.org/record/3975252. All other data are included in the article and/or supplementary materials. Thermodynamic code AlphaMELTS is free to download at https://magmasource.caltech.edu/alphamelts/.The author declares that they have no conflicts of interest.The author thanks Mingming Li, Yida Li, Michael Gurnis, Paul Asimow, Wei Mao for discussions. Thanks are extended to the editor Dr. J.Y Yin for handling of the manuscript and to reviewers for their insightful comments. Q. Yuan acknowledges support from the O.K. Earl Postdoctoral Fellowship at Caltech and partial support by NSF EAR-2330810. The numerical models were mostly performed on the Agave cluster at Arizona State University.