<figure><img src="image_1.png" /></figure>

HAL
open science

# Numerical models of the magmatic processes induced by slab breakoff

Rebecca Freeburn, Pierre Bouilhol, Ben Maunder, Valentina Magni, Jeroen van Hunen

► To cite this version:

Rebecca Freeburn, Pierre Bouilhol, Ben Maunder, Valentina Magni, Jeroen van Hunen. Numerical models of the magmatic processes induced by slab breakoff. Earth and Planetary Science Letters, 2017, 478, pp.203-213. 10.1016/j.epsl.2017.09.008. hal-01785194

HAL Id: hal-01785194

https://uca.hal.science/hal-01785194v1

Submitted on 23 Jul 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Numerical models of the magmatic processes induced by slab breakoff

Rebecca Freeburnª, Pierre Bouilholª,<fnref n="1" />, Ben Maunderª,<fnref n="2" />, Valentina Magniª,<fnref n="3" />, Jeroen van Hunenª*

**Abstract**

After the onset of continental collision, magmatism often persists for tens of millions of years, albeit with a different composition, in reduced volumes, and with a more episodic nature and more widespread spatial distribution, compared to normal arc magmatism. Kinematic modelling studies have suggested that slab breakoff can account for this post-collisional magmatism through the formation of a slab window and subsequent heating of the overriding plate and decompression melting of upwelling asthenosphere, particularly if breakoff occurs at depths shallower than the overriding plate.

To constrain the nature of any melting and the geodynamic conditions required, we numerically model the collision of two continental plates following a period of oceanic subduction. A thermodynamic database is used to determine the (de)hydration reactions and occurrence of melt throughout this process. We investigate melting conditions within a parameter space designed to generate a wide range of breakoff depths, timings and collisional styles.

Under most circumstances, slab breakoff occurs deeper than the depth extent of the overriding plate; too deep to generate any decompressional melting of dry upwelling asthenosphere or thermal perturbation within the overriding plate. Even if slab breakoff is very shallow, the hot mantle inflow into the slab window is not sustained long enough to sufficiently heat the hydrated overriding plate to cause significant magmatism. Instead, for relatively fast, shallow breakoff we observe melting of asthenosphere above the detached slab through the release of water from the tip of the heating detached slab. Melting of the subducted continental crust during necking and breakoff is a more common feature and may be a more reliable indicator of the occurrence of breakoff. We suggest that magmatism from slab breakoff alone is unable to explain several of the characteristics of post-collisional magmatism, and that additional geodynamical processes need to be considered when interpreting magmatic observations.

**Keywords**

slab breakoff; magmatism; continental collision; numerical modelling; subduction

**1 Introduction**

The nature and distribution of post-collisional magmas differ markedly from those corresponding to the pre-collisional history of convergence zones. Oceanic subduction creates linear belts of arcs dominated by calc-alkaline magma series, whereas post-collisional magmatism is more diverse chemically, and more widespread spatially. It can extend 100s of kilometres from the suture zone and can be sustained for millions of years after initial collision of the continental plates, although short-lived pulses of more intense melting and longer quiescent periods have been observed, as well as systematic variations in space and time (e.g. Chung et al., 2005; Pearce et al., 1990). Magmatism across

these areas appears to originate from a variety of different sources, displaying signatures characteristic of crustal, subcontinental lithospheric and asthenospheric material (e.g. Lee et al. 2012).

Many different mechanisms have been proposed to explain magmatism in continental collisional areas. These include processes involving solely the continental crust: shear heating along crustal faults (Harrison et al., 1998); melting of thickened crust (England and Thompson, 1986); and the exhumation of crust back to the surface (Harris and Massey, 1994). Additionally, processes may involve an influx of heat or melt from the underlying mantle: large scale lithospheric delamination (Bird, 1978); slab breakoff (Davies and von Blanckenburg, 1995); edge-driven convection (Kaislaniemi and van Hunen,

<footnote>ª Durham University, Stockton Road, Durham, DH1 3LE, UK</footnote>

<footnote n="1">¹ Present address: Centre de Recherches Pétrographiques et Géochimiques, 15 rue Notre Dame des Pauvres, 54500 Vandœuvre les Nancy, France.</footnote>

<footnote n="2">² Present address: Department of Earth Science and Engineering, Imperial College London, London, UK.</footnote>

<footnote n="3">³ Present address: The Centre for Earth Evolution and Dynamics (CEED), University of Oslo, Sem Sælands vei 24, PO Box 1048, Blindern, NO-0316 Oslo, Norway.</footnote>

<footnote>* Corresponding author: jeroen.van-hunen@durham.ac.uk</footnote>

<footnote># This paper is published as:</footnote>

<footnote>Rebecca Freeburn, Pierre Bouilhol, Ben Maunder, Valentina Magni, Jeroen van Hunen (2017). Numerical models of the magmatic processes induced by slab breakoff, *Earth and Planetary Science Letters*, Volume 478, Pages 203-213, doi: 10.1016/j.epsl.2017.09.008</footnote>
2014); and small-scale convection at the base of the lithosphere (Kaislaniemi et al., 2014). Slab breakoff after collision is often invoked to account for the changes in the nature and composition of magmas in collisional areas (e.g. Coulon et al., 2002; Mahéo et al., 2002; von Blanckenburg and Davies, 1995). Upon continental collision, the subducted slab may weaken due to thermal heating and/or increased extensional stresses between the more buoyant continental and denser oceanic lithospheres, resulting in slab breakoff and the formation of a slab window, and possibly leading to magmatism (van Hunen and Allen, 2011).

The dynamics of slab breakoff has been investigated extensively, and it has been shown that the strength of the subducting oceanic and continental lithospheres, in part influenced by the oceanic slab age, convergence velocity, continental crustal and lithospheric thicknesses, and the mechanism of detachment, all have a control on the depth of breakoff (Andrews and Billen, 2009; Duretz et al., 2011; Gerya et al., 2004). Models have shown a wide range in this depth, from 40 to over 500 km (Baumann et al., 2010; Duretz et al., 2011), and have shown the potential importance of a 3D geometry which allows for lateral progression of the break-off process through slab tearing (van Hunen and Allen, 2011; Magni et al., 2014b; Menant et al., 2016; Capitanio and Replumaz, 2013).

Slab breakoff may be fairly well constrained dynamically, but its consequences in terms of magmatism remain unclear. Previous studies have suggested that hot asthenosphere, flowing up through the slab window, may have the potential to melt through decompression and/or heat both the subducting crustal material and the overriding continental lithosphere (Davies and von Blanckenburg, 1995; van de Zedde and Wortel, 2001). Melting of the overriding lithosphere is viable where breakoff occurs at depths shallower than the base of this, but, asthenospheric melting may be less likely, as it requires very shallow breakoff at depths less than 50 km, unless the mantle is volatile rich.

Many previous slab breakoff studies used either fully dynamic models to study the complex mechanical process of breakoff and the associated mantle flow (e.g. Gerya et al., 2004), or used kinematic models of subduction to focus on the hydration state of the system and magmatic consequences (e.g. van de Zedde and Wortel, 2001). Recent numerical studies have started to use coupled petrological-thermomechanical numerical models to study the dynamics of collision and breakoff (Li et al., 2013; Menant et al., 2016); but a systematic investigation of the conditions required to generate breakoff-induced melting has not yet been conducted, which is undertaken in this study.

Continental collision involves a succession of complex, interacting dynamical and petrological processes, of which the contribution of slab breakoff to post-collisional magmatism is the least well understood and constrained. Therefore, we use coupled petrological-thermomechanical numerical models to examine conditions for post-collisional magmatism, focusing on the magmatic consequences of slab breakoff in particular, so that its effect can be better understood and dissociated from other processes. We will do so by varying a number of key parameters which are likely to affect the breakoff dynamics, and by tracking the occurrence of melting of a number of different materials over time. The dynamics in our models are in line with previous numerical studies and replicate many of the features that have been previously observed during collision and breakoff (Duretz and Gerya, 2013), and we thus use these models to explore the magmatic processes. We show that slab breakoff by itself is unlikely to explain all of the characteristics of post-collisional magmatism, and that additional processes may play an important role.

## 2 Method

### 2.1 Governing Equations

We use the Cartesian version of the finite element code Citcom to model oceanic and continental subduction (Moresi and Gurnis, 1996; Zhong et al., 2000). This solves for the conservation of mass, momentum, energy and composition in a fluid, assuming incompressible flow and adopting the Boussinesq approximations. A standard non-dimensionalisation is applied to the governing equations (Table 1):

$$x = x'h \quad t = t'h^2/\kappa \quad \eta = \eta'\eta_0 \quad T = T'\Delta T$$

Giving:

$$\nabla \cdot u = 0 \\ -\nabla P + \nabla \cdot (\eta(\nabla u + \nabla u^T)) + (RaT - \sum_{i=1}^{N} Rb_i C_i)e_z = 0 \\ \frac{\partial T}{\partial t} + u \cdot \nabla T = \nabla^2 T + Q_0 \\ \frac{\partial C}{\partial t} + u \cdot \nabla C = 0$$

The thermal Rayleigh number, Ra, is defined as:

$$Ra = \frac{\alpha \rho_0 g \Delta T h^3}{\kappa \eta_0}$$

and similarly the compositional Rayleigh number of a given composition i, $Rb_i$:
$$Rb_i = \frac{\delta \rho_i g h^3}{\kappa \eta_0}$$

where $\delta \rho_i$ is the intrinsic density of composition i relative to mantle material.

The conservation of mass, momentum and energy are solved using a finite element technique, whilst the conservation of composition is solved using a CVI particle-tracking technique (Wang et al., 2015). Model parameters are described in Table 1.

## 2.2 Rheology

Stress, temperature and pressure dependent rheology is applied throughout the model. The strength of material is reduced through three deformation mechanisms; dislocation creep, diffusion creep and brittle yielding. Furthermore, a maximum viscosity limit of $\eta_{max}$ is applied. The viscosity resulting from diffusion and dislocation creep are both described by:

$$\eta = A^{\frac{1}{n}} \dot{\epsilon}^{\frac{1-n}{n}} \exp \left( \frac{E^* + pV}{nRT} \right)$$

The composite viscosity caused by these viscous deformation mechanisms is described as:

$$\frac{1}{\eta_{visc}} = \frac{1}{\eta_{diff}} + \frac{1}{\eta_{disl}}$$

The effect of brittle yielding is also included:

$$\eta_y = \frac{\sigma_y}{\dot{\epsilon}}$$

where the yield stress, $\sigma_y$, is described as:

$$\sigma_y = \min(\sigma_0 + \mu p, \sigma_{max})$$

in which the effect of fluid pressure on the yield strength is ignored.

The effective viscosity is defined as the minimum of the viscosity components from viscous deformation, brittle yielding and $\eta_{max}$:

$$\eta_{eff} = \min(\eta_{visc}, \eta_y, \eta_{max})$$

A number of commonly used flow laws were used to determine a representative flow law for partially hydrated and hydrated mantle (Maunder et al., 2016). Given the expected profound influence of rheology on the dynamics of the continental collision process (e.g. van Hunen and Allen, 2011; Duretz and Gerya, 2013), mantle and crustal rheology are part of the

parameter investigation in this study, by investigating different mantle rheological flow laws (partially hydrated and hydrated, see Maunder et al., 2016), and crustal rheological flow laws (plagioclase & dry quartzite, see Ranalli (1995)). Flow laws for oceanic crust and sediment are adopted from Ranalli, (1995). Used values of rheological parameters for each composition are shown in Table 2.

## 2.3 Model setup

The development of continental collision and breakoff of the oceanic lithosphere is modelled through the introduction of a continental lithospheric block into a subduction zone (Fig. 1). As an initial condition, a subducting and overriding lithosphere are juxtaposed. The subducting lithosphere has an oceanic part, which, at one end, enters the subduction zone as a slab that is initially reaching down to a depth of 250 km underneath the continental overriding lithosphere. At its other end, a passive margin connects the oceanic lithosphere to a continental block, which extends to the left hand side boundary. To decouple the subducting and overriding plates, a 1.3 km thick layer of sediments with a fixed viscosity of 2 x 10²⁰ Pa s is positioned at the top of the subducting oceanic and continental lithospheres, which gradually increases to the background viscosity over a length of 13 km (see Maunder et al., 2016). Both oceanic and continental lithospheres are composed of crust and lithospheric mantle. The continental crust (assigned a typically observed thickness of either 35 or 40 km (Fowler, 2005) in the models), is split into an upper and lower crust, which have different rock compositions (affecting its dehydration and melting, see below), but otherwise have the same physical properties. Continental crust and sediments have a fixed relative density of -600 kg/m³ and oceanic crust -300 kg/m³ with respect to a reference mantle density of 3330 kg/m³. At depths greater than 80 km, oceanic crust is prescribed a new density of +200 kg/m³ to represent a simplified transition from basalt to eclogite. Oceanic lithosphere initially obeys the half space cooling thermal profile for a given plate age of 50 Myr. Continental lithosphere is assigned a commonly observed thicknesses of either 100 or 125 km (Fowler, 2005) in the presented models (see section 2.5) with a thermal profile from 0°C at the surface, to 1350°C at the base of the lithosphere. A constant initial mantle temperature of 1350°C is defined elsewhere. The rheological effect of adiabatic heating, which is absent in our models, is partly compensated for by using reduced values of the activation volume, V, in order to give a background viscosity profile in line with glacial isostatic adjustment data (Mitrovica and Forte,
2004). To investigate the effect of continental crust radiogenic heating, a uniform crustal value of 0.7 µWm⁻³ (Fowler, 2005) is included in some model calculations, along with an initially higher crustal temperature profile.

As boundary conditions, top and bottom boundaries have fixed temperatures of 0°C and 1350°C respectively, whilst side boundaries have a fixed continental geotherm (Fig. 1). To ensure trench retreat (i.e. an oceanward velocity of the overriding plate relative to the mantle underneath, as observed for the majority of subducting slabs worldwide (Funiciello et al., 2008)) during the initial slab development, a horizontal mantle flow of between 3 and 8 cm yr⁻¹ out of the model domain is applied initially through the right boundary (Fig. 1) until the slab reaches a depth of 600 km. Mantle material is allowed to flow into the model domain through the left boundary, where a zero horizontal deviatoric stress (σ'xx=0 and vz=0) is imposed. After this initial period, the right hand side boundary adopts this same zero horizontal deviatoric stress boundary condition, which removes any external forcing. This means that, after initialisation, subduction progresses unforced and flow is able to develop naturally with material free to enter/leave the model domain on both sides (Chertova et al., 2012). Due to a combination of the endothermic phase transition and viscosity jump between the upper and lower mantle and slab roll-back, many slabs in nature temporarily stagnate in the transition zone before penetrating into the lower mantle (Goes et al., 2017). Therefore, slabs are modelled to stagnate at the base of the upper mantle, and hence a closed bottom boundary is applied, with a horizontal velocity that matches the imposed right hand side mantle flow during the initial period, and no-slip afterwards. The top boundary is free slip to allow the natural horizontal (but no vertical) movement of the subduction trench.

The model domain represents an area of 1980 x 660 km and consists of 512 x 128 finite elements. These have a varying grid size to achieve the highest resolution at the top of the model, with an element size of 3.3 x 3.3 km² in this area, which is sufficient for the slab dynamics and melting features studied here. Compositional information is stored and advected using particles.

## 2.4 Hydration and melting implementation

We use a Gibbs free energy minimisation strategy to determine the metamorphic reactions, including (de)hydration reactions and melting reactions, occurring in the subducting slab and mantle. Equilibrium mineral assemblages are calculated beforehand using Perple_X (Connolly, 2005) for oceanic, upper and lower continental crustal, sediment, and primitive and depleted mantle lithologies, for a range of pressures, temperatures and water contents, and are stored in lookup tables (see Magni et al. (2014a) and Supplementary Material for full details). Since Boussinesq approximations are used, a correction for the mantle adiabat of 0.5°C/km is applied to the temperature profiles when determining the stable mineral assemblages. We use passive tracers that hold information about the composition to compute the stable phases at each time step. The presence and amount of melt under any given P-T-H₂O conditions can hence be determined and the variations in water content tracked. An initial hydration profile, with water stored in hydrous minerals, is defined for the oceanic and continental lithosphere within the subducting plate, which is advected with the flow field (Fig. 1). Oceanic crust and sediments are uniformly hydrated with 2.0 wt% H₂O, continental crust with 0.2 wt%, and the upper 8 km of oceanic lithospheric mantle with 4.0 wt% (Faccenda, 2014; Bouilhol et al., 2015). During subduction, dehydration might occur due to a change in the mineral assemblage. If free water becomes present, this water is allowed to percolate upwards, potentially rehydrating overlying material (see Bouilhol et al. (2015) and Magni et al. (2014a) for full details). During each timestep any free water is either re-absorbed into new hydrous minerals or is assumed to escape through the top of the model; hence we assume free water is capable of migrating significantly faster than the mantle flow.

To investigate the nature of any magmatism, we record the melting of six different compositions throughout the model calculations; the mantle lithosphere, asthenosphere, sediments, oceanic crust, and upper and lower continental crust. If melt is present and the amount exceeds a threshold, the melt is extracted instantaneously to the surface. Extraction threshold values of 2 and 8 wt% are chosen for mantle and crustal material respectively, based upon the amount required to form an interconnected melt pathway and ascend to the surface (Brown, 2007). Upon melt extraction, the residue will change in composition making it less likely to melt again. To account for this, we assume that material that has undergone melt extraction is unable to re-melt by marking it as “depleted” and removing it from further thermodynamic calculations for the
remainder of the model run. We record post-collisional melting, which is defined as any melting occurring after subduction-related melting has ceased. We also record breakoff-related melting, which is defined as melting which can be clearly linked to the necking and breakoff process, and is not a result of other post-collisional processes.

## 2.5. Model parameter investigation

We aim to analyse if and how breakoff-induced melting is generated, and constrain the geodynamic conditions required. To do so, we initially performed a pilot study on the dynamic effects of a large number of model parameters, including viscous strength (of the upper crust, lower crust, and mantle), brittle strength (surface yield strength, maximum yield strength, and friction coefficient), thickness of the overriding and subducting continental lithosphere, oceanic plate age, radiogenic heating, and the strength of the interplate weak zone. This list is unlikely to be exhaustive, and other model parameters, such as the width of the passive margin or slab penetration into the lower mantle, might affect slab breakoff dynamics and/or the magmatic consequences.

From this pilot study, we chose the five parameters that have the most profound effect on the collision dynamics. Here, we systematically investigate those parameters in terms of their control on the breakoff depth and/or timing, and record any changes that these dynamical differences produce in the style of melting. Models are conducted within this five-dimensional parameter space, with two values chosen for each of the following parameters; the ratio of continental crust to lithospheric thickness (35:125 km & 40:100 km), the radiogenic heating of the continental crust (0 & 0.7 μWm⁻³), the mantle rheological flow law (partially hydrated and hydrated), the crustal rheological flow law (plagioclase & dry quartzite) and the overriding lithospheric thickness (100 km & 60 km) (see Table 3). As will be further discussed below, these values are not intended to cover the full observed parameter range, but instead focus on conditions that are potentially relevant for slab breakoff related melting, which is the key objective of this study.

## 3 Results

We first describe the development of the dynamical process of collision and breakoff. To quantify these results, the timing of collision is defined as the moment when the continental lithosphere first reaches the subduction trench and breakoff as the moment when, due to a necking process, the 1000°C isotherm ceases to be continuous in the slab. We discuss any dehydration and melting for a reference

model (Model 1, Section 3.1), which will be shown to produce only small amounts of melt, and then for a model with a shallower breakoff depth that produces a significantly larger amount of melt (Model 13, Section 3.2). Next, we summarise the variation in breakoff dynamics and melting observed across the full parameter range of models within this study (Sections 3.3 and 3.4).

### 3.1 Reference model: breakoff dynamics and melting

Before examining and describing the effect of various model parameters that affect the dynamics and melting processes during continental collision, we first describe the reference model (Model 1) in more detail. Model 1 (Table 3, Fig. 2) adopts a plagioclase crustal flow law (Table 2), a mantle flow law that is appropriate for a partially hydrated mantle (Table 2), and has a roughly globally-averaged continental crust and plate thicknesses. An initial period of oceanic subduction (Fig. 2) and roll-back creates a fully developed subduction system, in which the mantle wedge and overriding lithosphere become hydrated. Next, the continental lithosphere enters the subduction zone (Fig. 2). The convergence velocity decreases due to the positive buoyancy of the continental crust, until the former passive margin stagnates at a depth of 150 km (Fig. 2A.ii, 2B.ii). The temperature of the stagnated slab increases, which weakens both the buoyant subducted continental crust and oceanic lithosphere, and results in necking within the oceanic lithosphere in this model (Fig. 2A.iii, 2B.iii). The subducted continental crust begins to flow upwards, whilst necking proceeds through viscous deformation and eventually results in breakoff at 44 Myr after initial collision at a depth of 240 km. In total, the process from initial necking to slab breakoff lasts for 14 Myr. Buoyant crustal flow continues and leads to a slight upwelling of asthenospheric mantle underneath the continental crust (Fig. 2A.iv, 2B.iv).

Throughout oceanic subduction, breakdown of hydrous minerals within the sediments, oceanic crust and lithospheric mantle result in mantle wedge melting and hydration of the overriding lithosphere (Fig. 2B.i), as expected from oceanic subduction. As the continental lithosphere enters the subduction zone, melting of the mantle wedge corner declines and then stops completely 10 Myr after initial collision (Fig. 2C). The oceanic crust and lithosphere continue to dehydrate throughout the early stages of collision, due to the continuously increasing slab temperature. However, the reduction in convergence velocity causes a steepening of the subducted slab, which results in the flux of water moving upwards into the subducted mantle lithosphere and continental crust (Fig. 2B.ii). Another consequence of this
steepening is the appearance of saturated melting in the oceanic crust at depths > 200 km due to fluid advection (Bouilhol et al. 2015), as the subducted slab begins to neck and crustal temperatures exceed 1000°C, but does not reach sufficient volumes to be extracted (i.e. > 8 wt%, Figs. 2B.ii and 2C). As the continental crust subsequently begins to be exhumed, small amounts of the deepest part of the subducted continental crust and sediments melt due to decompression melting (Figs. 2B.iv and 2C). The continental crust can also become significantly hydrated throughout the collisional process (~ 2 wt%), resulting in small amounts of upper crustal melting at shallow depths of 50-100 km (Fig. 2B.iv). This progresses for 15 Myr before the crust ceases exhuming and crustal temperatures decrease.

## 3.2 Shallow breakoff model: breakoff dynamics and melting

Model 13 (Table 3, Fig. 3) is characterised by a weaker crustal rheology, thicker continental crust and thinner continental lithosphere with respect to the reference model, and represents the other end of the spectrum of the explored parameter space, as it shows significant differences in breakoff depth, timing and location, and provides a clear example of the breakoff-related melting that we observe across the parametric study. The continental crust is subducted to similar depths after collision as compared to model 1. At the passive margin, necking and weakening of the crust and mantle lithosphere occurs, which progresses through both simple shearing and viscous necking (Figs. 3A.i and 3A.ii). Slab breakoff in this model occurs within the continental lithosphere at a significantly shallower depth of 170 km and with a shorter breakoff time of 15 Myr after initial collision (Figs. 3A.ii and 3A.iii).

A continuum of melting is observed, affecting different lithologies at different depth, whilst slab breakoff proceeds after collision. Melting occurs along the edge of the stagnated continental crust during necking of the slab (Fig. 3A.ii), and includes small amounts of sediment melting (Fig. 3B). This melting results from heating of continental crustal material within the tips of the detached slab as it descends into the mantle (Fig. 3A.iii). This crustal and sediment melting is most profound during the breakoff process and a few Myrs after (Fig. 3B), but small amounts (<10⁶ m³/m/Myr) of lower continental crust melting continues for another 30 Myrs. This heating also results in the breakdown of stable hydrous phases within crustal material in this area leading to deep (~ 150 km) asthenospheric melting, directly above the tip of the detached slab (Fig. 3A.iii). Asthenospheric melting begins 4 Myr after the cessation of subduction related magmatism, at

the time of breakoff (Fig. 3B). The locus of this melting progresses into the foreland region over time, due to the movement of the detached slab towards the subducting plate as it descends deeper into the mantle (Fig. 3A.iii). Asthenospheric melting can last only a few Myrs, since it is caused by the flux of water released from the detached part of the slab, which rapidly sinks into the mantle. A more prolonged period of continental crustal melting begins during crustal exhumation, at around the same time as breakoff, and is sustained for over 30 Myr, increasing as the exhumation advances and then decreasing as the plates become stable and cool (Figs. 3A.iv and 3D).

## 3.3 Breakoff dynamics across the parametric study

So the earlier, faster, and shallower break-off of the slab of Model 13 compared to Model 1 increases the chances of breakoff related melting. With this in mind, we will explore a range of model parameters to investigate their effect on post-collisional magmatism. We decided to focus on the conditions that lead to relatively shallow and rapid slab breakoff, which has been shown above to favour slab breakoff associated melting. As a consequence, most models in our chosen parameter range exhibit shallower and more rapid breakoff than the reference Model 1 (Fig. 4, Table 3). Breakoff localises where the tensile stresses, between the dense oceanic lithosphere and buoyant continental lithosphere, exceed the local strength of the slab. This condition depends upon a number of factors. Firstly, whilst within the oceanic lithosphere, these tensile stresses increase towards the surface, as the slab at that point supports more of the dense, subducted slab underneath. Secondly, the strength of the slab is in part controlled by the thickness of the slab and the strength of the mantle lithosphere and continental crust (e.g. compare continental slabs in Figs. 2B.ii and 3A.ii). Thirdly, the thickness of the cold slab core also determines the location of breakoff, which increases by thermal diffusion whilst the slab is stagnated. This is reduced by a thinner overriding plate and a shorter time between collision and breakoff, allowing the slab to start warming up at a shallower depth. Fourthly, the thickness of the continental crust plays a complex role. A thicker continental crust reduces the integrated strength of the lithosphere; in part because crustal material is intrinsically weaker than mantle material, but also because the strength of the continental crust can be further reduced by radiogenic heating. However, counterintuitively a weaker crustal rheology (and in some cases an increase in radiogenic heating) generally leads to an increase in breakoff depth. This results from a greater crust-mantle decoupling and earlier exhumation of the crust, locally reducing the
buoyancy of the continental part of the slab, and therefore also reducing the tensile stresses in the slab underneath (c.f. Baes et al., 2011). Full delamination of the mantle lithosphere from the continental crust can ensue in these cases (Fig. 5A), leading to slab roll-back before a final breakoff event that occurs when the slab has weakened sufficiently through heat diffusion.

Despite the wide range of breakoff depths that are observed (Fig. 4), breakoff does not occur at depths shallower than the base of the overriding lithospheric mantle for any of the conditions applied within the parametric study. To explore the consequences of even shallower breakoff and its possible magmatic effects, we add another model that lies outside the systematic parameter space, which induces breakoff at a depth of 83 km, by using a 50 km thick continental crust within a very thin, 70 km, continental lithosphere, with a crustal radiogenic heating of 0.7 µWm⁻³ (Fig. 5D, Model 33). However, other conditions, such as the subduction of a very young subducting oceanic plate, or using models where a low-temperature plasticity deformation mechanism, such as Peierl's mechanism is at work, may also lead to such shallow breakoff (see Duretz et al., 2011).

### 3.4 Melting across the parametric study

Modelled post-collisional melting occurs due to various dynamical processes, including slab breakoff (exposing lithosphere of the detached slab to hot mantle temperatures, leading to dehydration and saturated melting), decompression melting resulting from exhumation, and melting linked to delamination or extension. Below we consider melting patterns for each material that is directly related to slab breakoff and that reaches sufficient amounts to be extracted.

Continental crust melting from slab breakoff occurs across most of the parametric study. It shows no clear correlation with breakoff depth or location within oceanic or continental lithosphere (Fig. 4). Instead, it depends upon the maximum depth of the subducted continental crust relative to the thickness of the overriding plate. For breakoff within continental lithosphere, melting can occur around the necking and breakoff, which causes a local slab temperature increase (Fig. 3A.ii). Continental crust can still melt in some instances where breakoff occurs within the oceanic lithosphere from crustal heating within the upper detached slab after breakoff (Fig. 5B). This takes place when the overriding lithosphere is relatively thin, and hence the crust is in close proximity to hotter asthenosphere. For the same reason, melting of sediments also depends upon the thickness of the overriding plate. This melting is

widespread throughout our models, as some sediments are always carried down atop the subducting slab, independent of breakoff dynamics.

Asthenospheric melting is observed in a few cases (e.g. Models 12 - 16, Fig. 6). This melting is linked to the breakoff depth and timing (Figs. 4 and 6). Hydration of the inflowing mantle is essential for asthenospheric melting to occur. A shorter, faster breakoff increases the likelihood that sufficient water will be retained within the detached oceanic lithosphere after breakoff. The tip of detached continental crust is also able to store water which could potentially dehydrate. However, when melting occurs in continental crust, all the water is consumed, preventing free water flow into the inflowing asthenospheric mantle and induce asthenospheric melting.

Oceanic crustal melting is not a common feature across our parametric study. This only occurs in two cases (Models 9 and 10), for which breakoff occurs within the oceanic lithosphere (Fig. 4). After breakoff, the upper, stagnated portion of the slab continues to heat up, weakening the oceanic crust that is still present. As this eclogitised crust is denser than the surrounding mantle, it drips off and melts as it descends (Fig. 5B).

We do not observe any post-collisional melting in the overriding mantle lithosphere, except for Model 33, which shows small amounts of melt extraction ($<$5x10⁵ m³/m). This occurs at the very base of the overriding plate, to the right of the suture (Fig. 5D.ii) from hydration by water fluxing from the underlying slab. Melting at this location occurs before and after continental collision, and also during the process of slab necking and breakoff.

Melting is also generated through mechanisms that are not directly resulting from slab breakoff. Small amounts of decompression melting of the continental crust can occur during exhumation, which is not necessarily initiated by slab breakoff (Figs. 2B.iv and 3A.iv). Larger amounts of melting can be generated at the base of the continental crust from heating by upwelling asthenosphere during post-collisional extension (e.g. Model 11, Figs. 5C and 6). Finally, significant melting of continental crust and inflowing asthenosphere can occur during mantle lithospheric delamination in models with a weaker crustal rheology (e.g. Model 24, Figs. 5A and 6 (also for Models 19-23 and 27-32)).

## 4 Discussion

In convergence zones, changes in magmatic compositions are often interpreted in terms of changes in the dynamics of the system, such as a
change in dip angle (Gutscher et al., 2000), or onset of continental collision (e.g. Bouilhol et al. 2013), and after collision, abrupt shifts in composition are commonly hypothesised to represent slab breakoff (e.g. Coulon et al., 2002). Observed post-collisional melts are usually small in volume and distributed fairly sparsely compared to subduction-related melting. We have performed a series of 2D numerical models in which dehydration and melting processes are studied in the context of collision dynamics. The models are all two-dimensional, so the effect of processes like lateral slab tearing cannot be examined, but this does not significantly affect the timing and thermal conditions of the breakoff process (van Hunen and Allen, 2011). Our study enables us to better constrain the source and conditions required for such melting. Our results show that post-collisional magmatism due to slab breakoff is feasible, but perhaps not as easy to achieve as previously suggested by kinematic models. Below we discuss the conditions required to melt the different components of the system.

## 4.1 Lithospheric mantle melting

The models of Davies & von Blanckenburg (1995) show that the upwelling of asthenosphere upon breakoff can cause a thermal perturbation and melting of the overriding metasomatised mantle lithosphere (Fig. 7A). Subsequently, many post-collisional magmatic observations, with geochemical signatures typical of metasomatised mantle lithosphere, have been interpreted to be induced by the process of slab breakoff-induced melting (e.g. von Blanckenburg & Davies 1995; Ahmadzadeh et al. 2010). However, we do not observe any significant lithospheric mantle melting adjacent to the slab window upon slab breakoff in our study. Below we consider a number of possible reasons why breakoff, under the conditions imposed here, might be unable to induce a sufficient thermal perturbation to generate melting.

For melting to occur at the subduction interface (or resulting suture), breakoff needs to be at depths shallower than the base of the overriding lithosphere (Fig. 7A). Our results, and those of previous numerical studies, suggest breakoff at these depths is not common, and perhaps only occurs under unusual conditions (Fig. 4) (e.g. Duretz et al., 2011; van Hunen and Allen, 2011). Hence, hot asthenospheric material may not normally be brought into direct contact with the colder overriding lithosphere (Fig. 7B). So lithospheric melting near the suture, directly as a result of breakoff, may not be a common phenomenon, and the breakoff depth should be considered before using lithospheric melts as an indicator of direct exposure of overriding lithosphere to hot asthenosphere due to a slab window. However,

even for breakoff at depths shallower than the base of the overriding lithosphere, no significant melting of the lithosphere adjacent to the suture is observed (Fig. 5D). The model of lithospheric mantle melting relies on hydration (metasomatism) of the overlying lithosphere from previous oceanic subduction, which will lower the solidus, and enhances melting of the overlying lithosphere during a thermal perturbation. Although our models predict a wet solidus that is very close to those determined experimentally (e.g. Condamine and Médard, 2014) for metasomatic mantle, the thermal perturbation does not extend sufficiently into the overriding plate in order to significantly melt this metasomatised area (Figs. 5D.ii and 5D.iii). Davies and von Blanckenburg (1995) found that, after breakoff at 80 km depth, a 25-Myr period is required to heat a point only 15 km away from suture. Our model calculations agree with these findings, although the modelled rapid slab descent and subsequent plate cooling prevents the hydrated portion of the lithosphere from exceeding 800°C. Unless the hot asthenospheric flow is sustained for a longer period of time, it seems unlikely that the metasomatised mantle lithosphere would be able heat sufficiently to melt. Furthermore, the necking process during breakoff results in a thin film (< 10 km) of crustal material and sediments remaining below the overriding lithosphere, which further shields the lithospheric mantle from a more significant thermal perturbation.

Hence slab breakoff may not commonly be shallow enough to expose the overriding lithosphere to hotter mantle asthenosphere, and, other mechanisms may be at work to melt the mantle lithosphere in post-collisional areas. We do not rule out the general mechanism of slab breakoff-induced melting of the overriding lithosphere, but suggest that slab breakoff and subsequent slab window heating may not always be the most likely cause for melting of the mantle lithosphere. This is consistent with some observations from post-collisional areas, for example the Iranian plateau, where magmatism is prolonged, widespread and occurs hundreds of kilometres from the suture (e.g. Kaislaniemi et al., 2014).

## 4.2 Asthenospheric mantle melting

No dry decompressional melting of asthenosphere occurs in our study. Davies and von Blanckenburg (1995) also suggest breakoff depths <50 km are required, thus making dry asthenospheric melting unlikely. However, melting of asthenospheric mantle can occur where the mantle is percolated by metamorphic water derived from the slab tip. Indeed, slab breakoff provides a mechanism to heat marginally hydrated slab areas, releasing water that otherwise is likely to remain stable (Fig. 7B). In some cases the tip of the detached slabs melts, and
may contribute significantly to this metasomatic agent. This slab fluid/melt triggers asthenospheric melting above the slab tip. This melting is short-lived (a few Myrs), and synchronous with slab breakoff. Although the resulting magmas might be chemically distinct, with a strong slab affinity, (e.g. enrichment in potassium), it might be difficult to identify those magmas in the geological record due to their short-lived appearance. Nevertheless, potassic magmatism observed in the Mediterranean would fit both the spatial, temporal and compositional occurrence observed in our models (e.g. De Astis et al., 2006).

## 4.3 Crustal and sediment melting

Slab breakoff may result more easily in melting of the continental crust and sediments than mantle melting. The increase in temperature of this subducted material during slab necking and breakoff is sufficient to generate significant amounts of melt. This melting is not constrained to the more extreme cases of fast and shallow breakoff, such as for asthenospheric melting, but is widespread across our study. Instead it is primarily dependent upon the maximum depth of continental crustal subduction and overriding plate thickness. We would expect crustal melts to be extruded in localised, linear belts of magmatism near the suture, which may be a more reliable indicator of breakoff than observations of lithospheric mantle melts.

Some studies suggest that the observations of mafic melts may be due to a heating of oceanic lithosphere upon slab breakoff (e.g. Omrani et al., 2008). This is a rare feature in our models, and instead, melting of lower continental crust during breakoff is more likely, and may produce a similar geochemical signature.

## 4.4 Other post-collisional melting processes observed across the parametric study

Our models show several melting mechanisms not directly related to slab breakoff. Separation of crust and mantle lithosphere in subducting slabs occurs for most of our weaker (dry quartzite) continental crust models (Fig. 5A), due to reduced coupling between the crust and mantle lithosphere, in agreement with other numerical studies (e.g. Duretz and Gerya, 2013). This generates significant melting at the base of the delaminated crust (e.g. Göğüş and Pysklywec, 2008), and sometimes some asthenospheric melting, through crustal dehydration within the delaminating slab (e.g. Models 19-24 and 27-32, Fig. 6). Exhumation of the continental crust and in some cases subsequent extension can also lead to crustal melting through both decompression melting and heating around areas of hot inflowing asthenosphere

(Figs. 2B.iv, 3A.iv and 5C). These processes are not strictly breakoff generated, but may occur during or be accentuated by breakoff. In some collisional areas today (Arabia-Eurasia, India-Eurasia), convergence continues long after the initial collision (Molnar and Tapponnier, 1975), potentially suggesting additional compressional forces at work (e.g. Magni et al., 2017), that might suppress the extension observed in the presented models. However, delamination and exhumation have been suggested to occur and account for melting in some areas (e.g. Mohammadi et al., 2016). Additionally, melting may also have been generated through post-collisional extensional processes within older continental collisional settings, such as the Western Gneiss Region, Norway (e.g. Gordon et al., 2013).

## 5 Conclusions

The occurrence of, or changes in, post-collisional magmatism are often attributed to slab breakoff. In this study, we used numerical modelling to quantify any melting processes during post-collisional slab breakoff. Results suggest that the classical style of slab breakoff related magmatism, i.e. lithospheric melts from the exposure in a slab window created by slab breakoff (Fig. 7A), may not commonly occur, since breakoff depths are typically deeper than the overriding plate thickness, and may not induce a significant thermal perturbation in the overriding lithosphere. However, melting of the continental crust and sediments provides an alternative, viable mechanism for generating magmatism through breakoff (Fig. 7B). Crustal melts in collisional areas may commonly result from breakoff, and may be a more reliable indicator of slab necking and breakoff than mantle derived melts, provided these are not long-lived and located within narrow belts parallel to the suture. Inflowing asthenosphere can also melt through hydration from the detached, descending slab if breakoff occurs relatively shallowly, and <20 Myr after initial collision.

Slab breakoff cannot induce the observed prolonged periods of post-collisional magmatism, or those observed 100s of kms away from the suture zone. Therefore, it is unlikely that the breakoff process by itself can account for all observations of magmatism in post-collisional areas, and additional, alternative mechanisms are required to account for the diverse nature of magmatic observations and to explain some of the post-collisional magmatism that has previously been attributed to breakoff.

## Acknowledgments

The authors would like to thank Fabio Capitanio and Rob Govers for thorough and constructive reviews. This work has been supported by the European Research Council (ERC Starting Grant 279828) and
the FP7 Marie Curie Initial Training Network
Topomod, contract 264517. PB acknowledges his
Auvergne Fellowship (French Government
Laboratory of Excellence initiative no. ANR-
10LABX-0006, ClercVoc contribution no. 254). VM
also acknowledges support from the Research
Council of Norway through its Centres of Excellence
funding scheme, Project Number 223272.

**Appendix A. Supplementary material**

Supplementary material related to this article can be
found on-line at
https://doi.org/10.1016/j.epsl.2017.09.008.

**References**

Ahmadzadeh, G., Jahangiri, A., Lentz, D., Mojtahedi, M., 2010. Petrogenesis of Plio-Quaternary post-collisional ultrapotassic volcanism in NW of Marand, NW Iran. J. Asian Earth Sci. 39, 37–50.

Andrews, E.R., Billen, M.I., 2009. Rheologic controls on the dynamics of slab detachment. Tectonophysics 464, 60–69.

Baes, M., R. Govers, and R. Wortel (2011), Switching between alternative responses of the lithosphere to continental collision, Geophys. J. Int., 187(3), 1151–1174.

Baumann, C., Gerya, T. V., Connolly, J. a. D., 2010. Numerical modelling of spontaneous slab breakoff dynamics during continental collision. Geol. Soc. London, Spec. Publ. 332, 99–114.
doi:10.1144/SP332.7

Bird, P., 1978. Initiation of intracontinental subduction in the Himalaya. J. Geophys. Res. 83, 4975–4987.

Bouilhol, P., Jagoutz, O., Hanchar, J.M., Dudas, F.O., 2013. Dating the India–Eurasia collision through arc magmatic records. Earth and Planetary Science Letters 366, 163–175.

Bouilhol, P., Magni, V., van Hunen, J., Kaislaniemi, L., 2015. A numerical approach to melting in warm subduction zones. Earth Planet. Sci. Lett. 411, 37–44.

Brown, M., 2007. Crustal melting and melt extraction, ascent and emplacement in orogens: mechanisms and consequences. J. Geol. Soc. London. 164, 709–730. doi:10.1144/0016-76492006-171

Capitanio, F. A., and A. Replumaz (2013), Subduction and slab breakoff controls on Asian indentation tectonics and Himalayan western syntaxis formation, Geochem. Geophys. Geosyst., 14, 3515–3531, doi:10.1002/ggge.20171.

Chertova, M. V., Geenen, T., Van Den Berg, A., Spakman, W., 2012. Using open sidewalls for modelling self-consistent lithosphere subduction dynamics. Solid Earth 3, 313–326. doi:10.5194/se-3-313-2012

Chung, S.L., Chu, M.F., Zhang, Y., Xie, Y., Lo, C.H., Lee, T.Y., Lan, C.Y., Li, X., Zhang, Q., Wang, Y., 2005. Tibetan tectonic evolution inferred from spatial and temporal variations in post-collisional magmatism. Earth-Science Rev. 68, 173–196.

Condamine, P., Médard, E., 2014. Experimental melting

of phlogopite-bearing mantle at 1 GPa: Implications for potassic magmatism. Earth Planet. Sci. Lett. 397, 80–92. doi:10.1016/j.epsl.2014.04.027

Connolly, J.A.D., 2005. Computation of phase equilibria by linear programming: A tool for geodynamic modeling and its application to subduction zone decarbonation. Earth Planet. Sci. Lett. 236, 524–541. doi:10.1016/j.epsl.2005.04.033

Coulon, C., Megartsi, M., Fourcade, S., Maury, R.C., Bellon, H.B., Louni-Hacini, A., Cotten, J., Coutelle, A., Hermitte, D., 2002. Post-collisional transition from calc-alkaline to alkaline volcanism during the Neogene in Oranie (Algeria): magmatic expression of a slab breakoff. Lithos 62, 87–110.

Davies, J.H., von Blanckenburg, F., 1995. Slab breakoff: A model of lithosphere detachment and its test in the magmatism and deformation of collisional orogens. Earth Planet. Sci. Lett. 129, 85–102.

De Astis, G., Kempton, P.D., Peccerillo, A., Wu, T.W., 2006. Trace element and isotopic variations from Mt. Vulture to Campanian volcanoes: constraints for slab detachment and mantle inflow beneath southern Italy. Contrib. to Mineral. Petrol. 151, 331–351.

Duretz, T., Gerya, T. V., 2013. Slab detachment during continental collision: Influence of crustal rheology and interaction with lithospheric delamination. Tectonophysics 602, 124–140.
doi:10.1016/j.tecto.2012.12.024

Duretz, T., Gerya, T. V., May, D.A., 2011. Numerical modelling of spontaneous slab breakoff and subsequent topographic response. Tectonophysics 502, 244–256.

England, P.C., Thompson, A., 1986. Some thermal and tectonic models for crustal melting in continental collision zones. Geol. Soc. London, Spec. Publ. 19, 83–94.

Faccenda, M. 2014. *Water in the slab: A trilogy*, Tectonophysics 614, 1-30

Fowler, C.M.R., 2005. The solid Earth, 2nd ed. Cambridge University Press, Cambridge.

Funiciello, F., Faccenna, C., Heuret, A., Lallemand, S., Di Giuseppe, E., Becker, T.W., 2008. Trench migration, net rotation and slab-mantle coupling. Earth Planet. Sci. Lett. 271, 233–240.
doi:10.1016/j.epsl.2008.04.006

Gerya, T. V., Yuen, D.A., Maresch, W. V., 2004. Thermomechanical modelling of slab detachment. Earth Planet. Sci. Lett. 226, 101–116.

Goes, S., Agrusta, R., van Hunen, J. & Garel, F. 2017. *Subduction-transition zone interaction: A review*. Geosphere 13(3): 644-664.

Göğüş, O.H., Pysklywec, R.N., 2008. Near-surface diagnostics of dripping or delaminating lithosphere. J. Geophys. Res. Solid Earth 113, 1–11.
doi:10.1029/2007JB005123

Gordon, S.M., Whitney, D.L., Teyssier, C., Fossen, H., 2013. U-Pb dates and trace-element geochemistry of zircon from migmatite, Western Gneiss Region, Norway: Significance for history of partial melting in continental subduction. Lithos 170–171, 35–53.
doi:10.1016/j.lithos.2013.02.003

Gutscher, M.A., Maury, F., Eissen, J.P., Bourdon, E., 2000. Can slab melting be caused by flat subduction? Geology 28, 535–538.
doi:10.1130/0091-7613(2000)28<535:CSMBCB>2.0.CO;2

Harris, N., Massey, J., 1994. Decompression and anatexis of Himalayan metapelites. Tectonics 13, 1537–1546.

Harrison, T.M., Grove, M., Lovera, O.M., Catlos, E.J., 1998. Greater Himalayan Crystallines. J. Geophys. Res. 103, 27017–27032.

Kaislaniemi, L., van Hunen, J., 2014. Dynamics of lithospheric thinning and mantle melting by edge-driven convection: Application to Moroccan Atlas mountains. Geochemistry, Geophys. Geosystems 15, 3175–3189. doi: 10.1002/2014GC005414

Kaislaniemi, L., van Hunen, J., Allen, M.B., Neill, I., 2014. Sublithospheric small-scale convection-A mechanism for collision zone magmatism. Geology 42, 291–294

Lee, H.Y., Chung, S.L., Ji, J., Qian, Q., Gallet, S., Lo, C.H., Lee, T.Y., Zhang, Q., 2012. Geochemical and Sr-Nd isotopic constraints on the genesis of the Cenozoic Linzizong volcanic successions, southern Tibet. J. Asian Earth Sci. 53, 96–114.

Li, Z.-H., Xu, Z., Gerya, T., Burg, J.P., 2013. Collision of continental corner from 3-D numerical modeling. Earth Planet. Sci. Lett. 380, 98–111.

Magni, V., Bouilhol, P., van Hunen, J., 2014a. Deep water recycling through time. Geochemistry, Geophys. Geosystems 15, 4203–4216.

Magni, V., Faccenna, C., van Hunen, J., Funiciello, F., 2014b. How collision triggers backarc extension: Insight into mediterranean style of extension from 3-d numerical models. Geology 42, 511–514. doi:10.1130/G35446.1

Magni, V., Allen, M.B., van Hunen, J., Bouilhol, P., 2017. Continental underplating after slab break-off. Earth Planet. Sci. Lett. 474, 59–67. doi: 10.1016/j.epsl.2017.06.017.

Mahéo, G., Guillot, S., Blichtert-Toft, J., Rolland, Y., Pêcher, A., 2002. A slab breakoff model for the neogene thermal evolution of South Karakorum and South Tibet. Earth Planet. Sci. Lett. 195, 45–58. doi:10.1016/S0012-821X(01)00578-7

Maunder, B., van Hunen, J., Magni, V., Bouilhol, P., 2016. Relamination of mafic subducting crust throughout Earth’s history. Earth Planet. Sci. Lett. 449, 206–216. doi:10.1016/j.epsl.2016.05.042

Menant, A., Sternai, P., Jolivet, L., Guillou-Frottier, L., Gerya, T., 2016. 3D numerical modeling of mantle flow, crustal dynamics and magma genesis associated with slab roll-back and tearing: The eastern Mediterranean case. Earth Planet. Sci. Lett. 442, 93–107.

Mitrovica, J.X., Forte, A.M., 2004. A new inference of mantle viscosity based upon joint inversion of convection and glacial isostatic adjustment data. Earth Planet. Sci. Lett. 225, 177–189. doi:10.1016/j.epsl.2004.06.005

Mohammadi, A., Burg, J.P., Bouilhol, P., Ruh, J., 2016. U-Pb geochronology and geochemistry of Zahedan and Shah Kuh plutons, southeast Iran: Implication for closure of the South Sistan suture zone. Lithos 248–251, 293–308. doi:10.1016/j.lithos.2016.02.003

Molnar, P., Tapponnier, P., 1975. Cenozoic tectonics of Asia: Effects of a continental collision. Science (80-

.), 189, 419–426.

Moresi, L., Gurnis, M., 1996. Constraints on the lateral strength of slabs from three-dimensional dynamic flow models. Earth Planet. Sci. Lett. 38, 15–28.

Omrani, J., Agard, P., Whitechurch, H., Benoit, M., Prouteau, G., Jolivet, L., 2008. Arc-magmatism and subduction history beneath the Zagros Mountains, Iran: A new report of adakites and geodynamic consequences. Lithos 106, 380–398.

Pearce, J.A., Bender, J.F., De Long, S.E., Kidd, W.S.F., Low, P.J., Guner, Y., Saroglu, Y., Yilmaz, Y., Moorbath, S., Mitchell, J.G., 1990. Genesis of collision volcanism in Eastern Anatolia, Turkey. J. Volcanol. Geotherm. Res. 44, 182–229.

Ranalli, G., 1995. Rheology of the Earth. Springer.

van de Zedde, D.M.A., Wortel, M.J.R., 2001. Shallow slab detachment as a transient source of heat at midlithospheric depths. Tectonics 20, 868–882.

van Hunen, J., Allen, M.B., 2011. Continental collision and slab break-off: A comparison of 3-D numerical models with observations. Earth Planet. Sci. Lett. 302, 27–37. doi:10.1016/j.epsl.2010.11.035

von Blanckenburg, F., Davies, J.H., 1995. Slab breakoff: A model for syncollisional magmatism and tectonics in the Alps. Tectonics 14, 120–131.

Wang, H., Agrusta, R., Van Hunen, J., 2015. Advantages of a conservative velocity interpolation (CVI) scheme for particle-in-cellmethods with application in geodynamic modeling. Geochemistry Geophys. Geosystems 16, 2015–2023. doi: 10.1002/2015GC005824

Zhong, S., Zuber, M.T., Moresi, L., Gurnis, M., 2000. Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection. J. Geophys. Res. 105, 11063–11082.
<figure><img src="image_3.png" /><figcaption>Figures</figcaption></figure>

**Figure 1** - Initial setup of the reference Model 1, showing initial positioning of materials, boundary conditions and the initial hydration state of both the continental and oceanic lithospheres. Two right hand side and bottom boundary conditions are used: one during the initial forced slab roll-back period, and one for after this period. Legends apply to all subsequent figures as well. Abbreviations used: sed = sediments, UCC = upper continental crust, LCC = lower continental crust, OC = oceanic crust, LM = mantle lithosphere, AM = asthenospheric mantle.
<figure><img src="image_4.png" /></figure>

**Figure 2** - Dynamics and melting of Model 1 (reference model). Panels A.i-A.iv show the dynamics of Model 1, showing oceanic subduction, continental subduction, slab stagnation and final breakoff. Colours and black contours show temperature, and white contours outline continental upper and lower crust. The area inside orange boxes are detailed in panels B.i-B.iv. For the materials, hydration and melt colour legend, see Fig. 1. Times are given respective to the time of initial collision. The horizontal length for each inset plot is 200 km. Grey lines are 200-K temperature contours. The extracted melt through time is shown below for panel C. The vertical dashed lines show the times given in panels A.ii - A.iv, whilst the vertical solid lines defines the point of breakoff. Pale colours show where any melt is extracted within 1 Myr, whilst bright colours represent melt that exceeds a limit of 1x10⁶ m³/m/Myr. The volumes of melt are based on our thermodynamic database and defined in m³/m/Myr, but are only provided to compare the amount of melting observed through time between different models. For reference, if extrusive melt would spread evenly over a distance of 100 km, a volume of 1x10⁸ m³/m would correspond to a melt thickness of 1 km.
Model 13

<figure><img src="image_5.png" /></figure>

Figure 3 – Dynamics and melting of Model 13. Panels A.i to A.iv: as Figure 2B.i-B.iv. Panel B: as Figure 2C.
<figure><img src="image_6.png" /><figcaption>**Figure 4** – Breakoff depth versus breakoff time across the parametric study (see Table 3 for details of each model). Models which show melt extraction over $1\times10^6 \text{ m}^3/\text{m}$ as a direct result of breakoff are highlighted for the asthenospheric mantle, upper and lower continental crust, sediments, and oceanic crust. Squares refer to models with breakoff in the oceanic part of the slab and circles indicate breakoff in the continental part of the slab. Models 1-32 are those within the parametric study, whilst Model 33 is an anomalous case, designed to induce very shallow breakoff.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 5** - Compilation of other key model dynamics (see Table 3 for details of each model). A - Model 24 (which uses a weaker (dry quartzite) continental crustal flow law): crust-mantle delamination, resulting in both asthenospheric mantle melting and melting of the delaminating lower continental crust. B - Model 10: melting of the oceanic crust if breakoff occurs within the oceanic lithosphere, due to the later dripping away of denser, eclogitised oceanic crust. C - Model 11: Post-collisional extension, which leads to upwelling of hot asthenospheric mantle and melting of the continental crust. D - Model 33: Shallow breakoff where break-off occurs at ~80 km, showing no significant melting of the lithospheric mantle adjacent to the suture, nor any significant thermal perturbation in the lithosphere at the location of breakoff. Large areas of melting (although with quantities less than the extraction threshold) are observed towards the base of the crust within the subducting continental lithosphere. This is a result of the higher crustal geotherm created by the imposed model conditions, and is not linked to the collision or breakoff process.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6** – A comparison of the amount of melt extraction for all models. Bright colours represent breakoff related melting, whilst colours with white hatching represent the remaining total post-collisional magmatism, from processes such as crustal exhumation, extension and delamination (see Fig. 5 for examples of these). Significant crust-mantle delamination occurs in Models 19-24 and 27-32, which can induce large volumes of continental crustal and asthenospheric melting.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7** - Schematic diagram of slab breakoff. A –possible melting induced by very shallow breakoff (similar to that suggested by Davies and von Blanckenburg (1995)), and B – possible melting induced by deeper breakoff, as observed in this study.</figcaption></figure>
<table><caption>Tables</caption><thead><tr><th>Symbol</th><th>Parameter and units</th><th>Default value</th></tr></thead><tbody><tr><td>u</td><td>velocity [m s<sup>-1</sup>]</td><td>-</td></tr><tr><td>P</td><td>pressure [Pa]</td><td>-</td></tr><tr><td>η</td><td>viscosity [Pa s]</td><td>-</td></tr><tr><td>T</td><td>temperature [°C]</td><td>-</td></tr><tr><td>Q<sub>0</sub></td><td>radiogenic heat production [µW m<sup>-3</sup>]</td><td>see Table 3</td></tr><tr><td>Ra</td><td>thermal Rayleigh number [-]</td><td>2.8 x 10<sup>6</sup></td></tr><tr><td>Rb</td><td>compositional Rayleigh number [-]</td><td>-</td></tr><tr><td>C</td><td>compositional function [-]</td><td>-</td></tr><tr><td>e<sub>z</sub></td><td>vertical unit vector [-]</td><td>-</td></tr><tr><td>t</td><td>time [s]</td><td>-</td></tr><tr><td>α</td><td>thermal expansion coefficient [K<sup>-1</sup>]</td><td>3.5 x 10<sup>-5</sup></td></tr><tr><td>ρ<sub>0</sub></td><td>reference density [kgm<sup>-3</sup>]</td><td>3330</td></tr><tr><td>g</td><td>gravitational acceleration [m s<sup>-2</sup>]</td><td>9.8</td></tr><tr><td>h</td><td>model height [km]</td><td>660</td></tr><tr><td>κ</td><td>thermal diffusivity [m<sup>2</sup> s<sup>-1</sup>]</td><td>8 x 10<sup>-7</sup></td></tr><tr><td>η<sub>0</sub></td><td>reference viscosity [Pa s]</td><td>2 x 10<sup>20</sup></td></tr><tr><td>δρ</td><td>density difference</td><td>see Table 2</td></tr><tr><td>A</td><td>pre-exponential exponent [Pa<sup>n</sup> s<sup>1</sup>]</td><td>see Table 2</td></tr><tr><td>x</td><td>distance [m]</td><td>-</td></tr><tr><td>ΔT</td><td>maximum temperature difference [°C]</td><td>1350</td></tr><tr><td>n</td><td>power-law exponent [-]</td><td>see Table 2</td></tr><tr><td>ε̇</td><td>second invariant of strain rate [s<sup>-1</sup>]</td><td>-</td></tr><tr><td>E*</td><td>activation energy [Jmol<sup>-1</sup>]</td><td>see Table 2</td></tr><tr><td>V<sub>diff</sub></td><td>activation volume [cm<sup>3</sup> mol<sup>-1</sup>]</td><td>0.7 (diff.), 8.0 (disl.)</td></tr><tr><td>R</td><td>gas constant [J K<sup>-1</sup> mol<sup>-1</sup>]</td><td>8.3</td></tr><tr><td>η<sub>max</sub></td><td>maximum viscosity [Pa s]</td><td>1 x 10<sup>24</sup></td></tr><tr><td>σ<sub>y</sub></td><td>yield stress</td><td>-</td></tr><tr><td>σ<sub>0</sub></td><td>surface yield stress [MPa]</td><td>40</td></tr><tr><td>σ<sub>max</sub></td><td>maximum yield stress [MPa]</td><td>400</td></tr><tr><td>μ</td><td>friction coefficient</td><td>0.2</td></tr><tr><td>N</td><td>number of materials with differing densities</td><td>4</td></tr></tbody></table>

<table><caption>Table 1 - Parameters, symbols, units and default values used within this study</caption><thead><tr><th rowspan="2">Material</th><th colspan="6">Rheological Values</th><th rowspan="2">Flow law</th><th rowspan="2">Density relative to mantle (kg m<sup>-3</sup>)</th></tr><tr><th>A<sub>diff</sub> (Pa s)</th><th>E*<sub>diff</sub> (Jmol<sup>-1</sup>)</th><th>n<sub>diff</sub></th><th>A<sub>disl</sub> (Pa<sup>n</sup> s)</th><th>E*<sub>disl</sub> (Jmol<sup>-1</sup>)</th><th>n<sub>disl</sub></th></tr></thead><tbody><tr><td>Mantle</td><td>3.79x10<sup>10</sup><br>1.17x10<sup>11</sup></td><td>3.19 x 10<sup>5</sup><br>2.96 x 10<sup>5</sup></td><td>1<br>1</td><td>1.47 x10<sup>16</sup><br>2.30 x10<sup>16</sup></td><td>4.96 x 10<sup>5</sup><br>4.73 x 10<sup>5</sup></td><td>3.48<br>3.47</td><td>PH<sup>1</sup><br>Hydrated<sup>1</sup></td><td>-</td></tr><tr><td>Oceanic crust</td><td>-</td><td>-</td><td>-</td><td>1.26 x10<sup>24</sup></td><td>2.60 x 10<sup>5</sup></td><td>3.40</td><td>Diabase<sup>2</sup></td><td>-300 /+200<br>(basalt/eclogite)</td></tr><tr><td>Continental crust</td><td>-</td><td>-</td><td>-</td><td>4.80 x10<sup>22</sup><br>3.75 x10<sup>19</sup></td><td>2.38 x 10<sup>5</sup><br>1.56 x 10<sup>5</sup></td><td>3.20<br>2.40</td><td>Plagioclase<sup>2</sup><br>Dry Quartzite<sup>2</sup></td><td>-600</td></tr><tr><td>Sediments</td><td>-</td><td>-</td><td>-</td><td>1.97 x10<sup>17</sup></td><td>1.54 x 10<sup>5</sup></td><td>2.30</td><td>Wet Quartzite<sup>2</sup></td><td>-600</td></tr></tbody></table>

<caption>Table 2 - Rheological values, flow laws and densities used for materials in this study. PH = Partially Hydrated. 1 - (Maunder et al., 2016); 2 - (Ranalli, 1995).</caption>
<table><caption>**Table 3** - Parameters and breakoff information for all models within the parametric study. Non-reference values are highlighted in grey. For the breakoff location, O represents within the subducted oceanic lithosphere and C, within the subducted continental lithosphere. Abbreviations: CC=continental crust, DQ=Dry Quarzite, PH=Partially Hydrated</caption><thead><tr><th>Model</th><th>CC flow law</th><th>Mantle flow law</th><th>Continental crustal/ lithospheric thickness (km)</th><th>Radio- genic heating (µWm<sup>-3</sup>)</th><th>Overriding lithospheric thickness (km)</th><th>Breakoff depth (km)</th><th>Breakoff time post- collision (Myr)</th><th>Breakoff location (O/C)</th></tr></thead><tbody><tr><td>1</td><td>Plagioclase</td><td>PH</td><td>35/125</td><td>0</td><td>100</td><td>241</td><td>44.3</td><td>O</td></tr><tr><td>2</td><td>Plagioclase</td><td>PH</td><td>35/125</td><td>0</td><td>60</td><td>224</td><td>26.7</td><td>O</td></tr><tr><td>3</td><td>Plagioclase</td><td>PH</td><td>35/125</td><td>0.7</td><td>100</td><td>191</td><td>34.0</td><td>C</td></tr><tr><td>4</td><td>Plagioclase</td><td>PH</td><td>35/125</td><td>0.7</td><td>60</td><td>172</td><td>25.3</td><td>C</td></tr><tr><td>5</td><td>Plagioclase</td><td>PH</td><td>40/100</td><td>0</td><td>100</td><td>175</td><td>24.2</td><td>C</td></tr><tr><td>6</td><td>Plagioclase</td><td>PH</td><td>40/100</td><td>0</td><td>60</td><td>159</td><td>16.6</td><td>C</td></tr><tr><td>7</td><td>Plagioclase</td><td>PH</td><td>40/100</td><td>0.7</td><td>100</td><td>178</td><td>17.6</td><td>C</td></tr><tr><td>8</td><td>Plagioclase</td><td>PH</td><td>40/100</td><td>0.7</td><td>60</td><td>149</td><td>12.3</td><td>C</td></tr><tr><td>9</td><td>Plagioclase</td><td>Hydrated</td><td>35/125</td><td>0</td><td>100</td><td>241</td><td>24.6</td><td>O</td></tr><tr><td>10</td><td>Plagioclase</td><td>Hydrated</td><td>35/125</td><td>0</td><td>60</td><td>195</td><td>14.7</td><td>O</td></tr><tr><td>11</td><td>Plagioclase</td><td>Hydrated</td><td>35/125</td><td>0.7</td><td>100</td><td>165</td><td>17.0</td><td>C</td></tr><tr><td>12</td><td>Plagioclase</td><td>Hydrated</td><td>35/125</td><td>0.7</td><td>60</td><td>135</td><td>10.5</td><td>C</td></tr><tr><td>13</td><td>Plagioclase</td><td>Hydrated</td><td>40/100</td><td>0</td><td>100</td><td>169</td><td>14.8</td><td>C</td></tr><tr><td>14</td><td>Plagioclase</td><td>Hydrated</td><td>40/100</td><td>0</td><td>60</td><td>145</td><td>9.2</td><td>C</td></tr><tr><td>15</td><td>Plagioclase</td><td>Hydrated</td><td>40/100</td><td>0.7</td><td>100</td><td>168</td><td>10.0</td><td>C</td></tr><tr><td>16</td><td>Plagioclase</td><td>Hydrated</td><td>40/100</td><td>0.7</td><td>60</td><td>129</td><td>5.6</td><td>C</td></tr><tr><td>17</td><td>DQ</td><td>PH</td><td>35/125</td><td>0</td><td>100</td><td>221</td><td>44.8</td><td>C</td></tr><tr><td>18</td><td>DQ</td><td>PH</td><td>35/125</td><td>0</td><td>60</td><td>241</td><td>33.1</td><td>O</td></tr><tr><td>19</td><td>DQ</td><td>PH</td><td>35/125</td><td>0.7</td><td>100</td><td>208</td><td>51.7</td><td>C</td></tr><tr><td>20</td><td>DQ</td><td>PH</td><td>35/125</td><td>0.7</td><td>60</td><td>231</td><td>28.0</td><td>C</td></tr><tr><td>21</td><td>DQ</td><td>PH</td><td>40/100</td><td>0</td><td>100</td><td>191</td><td>23.2</td><td>C</td></tr><tr><td>22</td><td>DQ</td><td>PH</td><td>40/100</td><td>0</td><td>60</td><td>188</td><td>17.9</td><td>C</td></tr><tr><td>23</td><td>DQ</td><td>PH</td><td>40/100</td><td>0.7</td><td>100</td><td>162</td><td>16.7</td><td>C</td></tr><tr><td>24</td><td>DQ</td><td>PH</td><td>40/100</td><td>0.7</td><td>60</td><td>155</td><td>11.5</td><td>C</td></tr><tr><td>25</td><td>DQ</td><td>Hydrated</td><td>35/125</td><td>0</td><td>100</td><td>215</td><td>27.6</td><td>C</td></tr><tr><td>26</td><td>DQ</td><td>Hydrated</td><td>35/125</td><td>0</td><td>60</td><td>155</td><td>15.0</td><td>C</td></tr><tr><td>27</td><td>DQ</td><td>Hydrated</td><td>35/125</td><td>0.7</td><td>100</td><td>231</td><td>17.6</td><td>C</td></tr><tr><td>28</td><td>DQ</td><td>Hydrated</td><td>35/125</td><td>0.7</td><td>60</td><td>215</td><td>11.4</td><td>C</td></tr><tr><td>29</td><td>DQ</td><td>Hydrated</td><td>40/100</td><td>0</td><td>100</td><td>178</td><td>15.2</td><td>C</td></tr><tr><td>30</td><td>DQ</td><td>Hydrated</td><td>40/100</td><td>0</td><td>60</td><td>162</td><td>7.1</td><td>C</td></tr><tr><td>31</td><td>DQ</td><td>Hydrated</td><td>40/100</td><td>0.7</td><td>100</td><td>139</td><td>8.9</td><td>C</td></tr><tr><td>32</td><td>DQ</td><td>Hydrated</td><td>40/100</td><td>0.7</td><td>60</td><td>125</td><td>5.4</td><td>C</td></tr><tr><td>33</td><td>Plagioclase</td><td>Hydrated</td><td>50/75</td><td>0.7</td><td>100</td><td>83</td><td>1.8</td><td>C</td></tr></tbody></table>