Contributions of Grain Damage, Thermal Weakening, and Necking to Slab Detachment

We investigate the impact of three coupled weakening mechanisms on the viscous detachment of a stalled lithospheric slab: structural weakening due to necking, material weakening due to grain size reduction, using a two-phase grain damage model, and thermal weakening due to shear heating (thermal damage). We consider a combined flow law of dislocation and diffusion creep. To understand and quantify the coupling of these three nonlinear weakening processes, we derive a mathematical model, which consists of three coupled nonlinear ordinary differential equations describing the evolution of slab thickness, grain size, and temperature. With dimensional analysis, we determine the dimensionless parameters which control the relative importance of the three weakening processes and the two creep mechanisms. We derive several analytical solutions for end-member scenarios that predict the detachment time, that is the duration of slab detachment until slab thickness becomes zero. These analytical solutions are then tested against numerical solutions for intermediate cases. The analytical solutions are accurate for end-member scenarios where one of the weakening mechanisms and one of the creep mechanisms is dominant. Furthermore, we use numerical solutions of the system of equations to systematically explore the parameter space with a Monte Carlo approach. The numerical approach shows that the analytical solutions typically never deviate by more than 50% from the numerical ones, even for scenarios where all three weakening and both creep mechanisms are important. When both grain and thermal damage are important, the two damage processes generate a positive feedback loop resulting in the fastest detachment times. For Earth conditions, we find that the onset of slab detachment is controlled by grain damage and that during later stages of slab detachment thermal weakening becomes increasingly important and can become the dominating weakening process. We argue that both grain and thermal damage are important for slab detachment and that both damage processes could also be important for lithosphere necking during continental rifting leading to break-up and ocean formation.


INTRODUCTION
Since the establishment of plate tectonic theory and early seismicity observations at subduction systems (Isacks and Molnar, 1969), slab detachment, also often called slab break-off, is considered as a prominent process associated with subduction systems. During slab detachment, the lower part of a subducting lithospheric plate detaches from its upper portion and sinks into the mantle. There is overall agreement that buoyancy forces within subducting slabs induce extensional stresses which can result in slab detachment. This process has been suggested to cause seismicity patterns in subduction zones (Isacks and Molnar, 1969), magmatism during orogeny (Davies and von Blanckenburg, 1995), rapid surface uplift (e.g., Fox et al., 2015;Sternai et al., 2019), the transition from Flysch to Molasse sedimentation in the European Alps (Schlunegger and Kissling, 2015), and abrupt changes in plate motions (McKenzie, 1969;Bercovici et al., 2015). The abruptness of these events suggests that slab detachment has to occur in equally short timescales. The occurrence of slab detachment has been suggested for various regions, among them the Mediterranean (Wortel and Spakman, 2000), the Pamir-Hindu Kush region (Kufner et al., 2016), and the European Alps (Lippitsch et al., 2003). However, geodynamic interpretations that consider slab detachment as controlling process are frequently contentious (see Garzanti et al., 2018, for a review). Also, indirect observations of slab detachment based on seismic tomography are sometimes contradictory, such as in the Central and Western Alps (Kästle et al., 2020), where seismic tomography is interpreted as showing a detached slab (Lippitsch et al., 2003) but also a continuous slab (Zhao et al., 2016). Some of the controversies concerning slab detachment are also due to the fact that the thermo-mechanical process of slab detachment is still incompletely understood.
Using simplified models, Wong A Ton and Wortel (1997) found that the rheology of the subducting lithosphere is the most important parameter affecting slab detachment, while buoyancy and resistive forces are less significant. To investigate the impact of additional processes, such as brittle failure, mineral phase transitions, partial melting, grain size reduction, or shear heating, and to study the inherent three-dimensionality of slab detachment, complex numerical models have been conducted in both 2D (e.g., Gerya et al., 2004;Andrews and Billen, 2009;Duretz et al., 2011) and 3D (e.g., Burkett and Billen, 2010;van Hunen and Allen, 2011;Capitanio and Replumaz, 2013;Duretz et al., 2014;von Tscharner et al., 2014;Chen and Gerya, 2016;Magni et al., 2017;Pusok et al., 2018). In general, these models confirm the finding of Wong A Ton and Wortel (1997) that slab rheology is the most important parameter controlling slab detachment. However, when standard olivine rheologies are employed, slab detachment is usually slow and occurs on time scales of tens of millions of years, thus resulting in large detachment depths (e.g., Schott and Schmeling, 1998;van Hunen and Allen, 2011), contrary to what is suggested from tomographic images and seismicity patterns. However, additional mechanisms have been suggested to significantly accelerate slab detachment, such as weakening due to grain size reduction (Bercovici et al., 2015Bellas et al., 2018), elaborated viscoelastoplastic rheology models (Duretz et al., 2011(Duretz et al., , 2012, preexisting weak zones (Burkett and Billen, 2010) or thermal weakening due to shear heating (Gerya et al., 2004;Duretz et al., 2012). Shear heating can efficiently localize deformation and form large-scale shear zones under lithospheric deformation conditions (Thielmann and Kaus, 2012;Duretz et al., 2015Duretz et al., , 2019Kiss et al., 2019). In particular, the combined action of shear heating and grain size reduction may result in significant weakening of the slab and associated strain localization (Thielmann et al., 2015;Foley, 2018), which may even trigger earthquakes in slab detachment settings (Thielmann, 2018). Therefore, it seems important to include different ductile weakening effects in models of slab detachment. Shear heating and grain size reduction have been shown to result in comparable weakening (Duretz et al., 2012;Bercovici et al., 2015), but it is not yet clear which of these mechanisms ultimately control slab detachment under what conditions and whether these processes are competitive or collaborative processes (Kameyama et al., 1997;Thielmann et al., 2015;Foley, 2018;Thielmann, 2018).
While 2D and 3D numerical models capture the complex geometrical interactions during slab detachment, their computational expense makes it difficult to systematically investigate the impact of different weakening mechanisms and their mutual feedbacks. Furthermore, due to the many equations and parameters involved, such 2D and 3D models are often little transparent. In contrast, simple 0D and 1D models are usually more transparent and have proven to provide fundamental physical insight into the slab detachment process and to efficiently test the impact of different model parameters (Schmalholz, 2011;Duretz et al., 2012;Bercovici et al., 2015Bercovici et al., , 2019Ribe and Xu, 2020).
While slab detachment is likely dominated by viscous necking, it has also been invoked to explain deep seismicity in various regions on Earth, such as e.g., the Hindu-Kush (Kufner et al., 2017), the Vrancea region (Mitrofan et al., 2016), or the Alboran slab in the Mediterranean (Sun and Bezada, 2020). The mechanisms resulting in the observed seismicity are still contentious, as brittle failure at this depth seems to be unlikely due to the high pressures and temperatures. However, since rocks are intrinsically viscoelastic, a significant local increase of strain rates during slab detachment could shift the rock deformation mechanism from effectively viscous to effectively elastic, and hence potentially cause brittle failure. Currently, it is still unclear which mechanism(s) result in seismicity generation, but two mechanisms are frequently favored, namely dehydration embrittlement and thermal runaway (see Zhan, 2019, and citations therein). Some intermediate-depth earthquakes pose, despite their depth, a significant seismic hazard (Frohlich, 2006) and it is hence important to identify their rupture mechanisms and the conditions favoring their occurrence.
Here, we investigate viscous slab detachment for a combination of dislocation creep and diffusion creep flow laws. We quantify the coupling of three nonlinear weakening mechanisms during slab detachment: (1) structural weakening due to a necking instability (see Schmalholz and Mancktelow, 2016, and citations therein) causing localized thinning of the slab, (2) weakening due to grain size reduction and grain size sensitive diffusion creep, referred to here as grain damage, and (3) weakening due to shear heating and temperature sensitive effective rock viscosities, referred to here as thermal damage. We consider a simple slab detachment setting and use a simple 0D mathematical model that allows us to comprehensively examine the parameter space. We derive a system of three coupled nonlinear differential equations that describe the evolution of grain size, temperature and slab thickness. We apply dimensional analysis to determine the controlling nondimensional parameters. We further derive several analytical solutions for the duration of slab detachment, that is when the slab thickness is zero, for end-member scenarios where one of the weakening mechanisms or one of the two flow laws is dominant. Using a Monte Carlo simulation with 2 × 10 5 simulations, we are able to determine the most important parameters controlling slab detachment. Based on the results from these simulations, we derive empirical relationships that can be used to predict detachment times for a given set of parameters. Finally, we determine which mechanism is most likely to govern slab detachment.

METHODS
In this section, we introduce our detachment model (section 2.1), together with the equations governing the evolution of slab thickness. In sections 2.2-2.4, we introduce the governing equations for rheology, for the evolution of grain interface curvature (as a proxy of grain size) and for the temperature evolution. In section 2.5, we non-dimensionalize the governing equations to derive a set of non-dimensional equations and parameters that will be used to investigate the impact of different parameters on the detachment process. The meaning of the derived non-dimensional parameters is discussed in section 2.6.

Detachment Model
We consider a vertical slab with initially zero velocity attached to a rigid surface; a so-called stalled slab. The slab has a length H, width D 0 and density ρ 1 , which is larger than the density ρ 2 of the surrounding mantle (see also Figure 1). Due to this density difference ρ = ρ 1 − ρ 2 , a negative buoyancy force per unit length, F B = ρgD 0 H, is induced, which causes the slab to sink and to stretch. For this configuration, localized thinning and necking will develop in the upper region of the slab where the buoyancy-driven stresses are highest. The extensional strain ratė ε is related to the temporal evolution of the slab thickness D at the necking region by equation (Schmalholz, 2011): where t is time. As the buoyancy force remains constant and the slab is considered much stronger than the surrounding, the deviatoric stress in the necking region is half the total stress, F B /D, and is given by (Turcotte and Schubert, 2014)  The initial state of the system is thus given by the initial stress τ 0 = 1 2 ρgHD 0 D 0 and the mean temperature T 0 of the slab. We will link the above two equations by rheological flow laws, which are mathematical relations betweenε and τ and are described in the next section.

Rheology
The buoyancy force F B induces deformation in the stalled slab. We assume that deformation is governed by a viscous composite rheology consisting of grain size insensitive dislocation and grain size sensitive diffusion creep (e.g., Hirth and Kohlstedt, 2003;Billen and Hirth, 2007;Kohlstedt and Hansen, 2015): where the strain rateε is a function of stress τ , temperature T, and grain size R. A dis and A dif are rheological prefactors, n is the stress exponent of dislocation creep and m the grain size exponent (see also Table 1 for the definition of the different rheological parameters). Q dis and Q dif are the respective activation enthalpies divided by the gas constant R given by: where P is pressure and E i , V i are activation energy and volume, respectively. The index i denotes either diffusion or dislocation creep. It is useful to define a transition grain size R t at which the contributions from dislocation and diffusion creep equally contribute to deformation: Frontiers in Earth Science | www.frontiersin.org Creep parameters are taken from Hirth and Kohlstedt (2003), while interface coarsening parameters are taken from Bercovici and Ricard (2012).
Using this transition grain size, Equation (3) can be rewritten as: where the ratio R t R indicates whether dislocation creep ( R t R < 1) or diffusion creep ( R t R > 1) is the dominant deformation mechanism.

Grain Size Evolution
Grain size can significantly affect rock rheology if it is sufficiently small and if it changes with progressive deformation. The evolution of grain size is therefore crucial to determine the transient response of the rock to an applied stress. The reduction of the effective viscosity, quantifying viscous strength, due to a reduction in grain size is frequently termed grain damage. As mantle rocks are polymineralic, consisting of mostly olivine and pyroxene, we here treat grain size evolution using the twophase grain size evolution law of Bercovici and Ricard (2012). This evolution law takes into account that a rock does not consist of a single mineral phase, but of at least two phases, which in the case of the mantle lithosphere corresponds to an olivine-pyroxene mixture. Due to the presence of pyroxene, the grain boundary mobility of olivine is strongly reduced due to Zener pinning. Due to this effect, the curvature of the grain boundaries is also significantly distorted, as illustrated in Figure 2. If we consider a single grain with radius R 1 , its interface curvature without the effect of secondary particles would equal R 1 . However, due to the presence of secondary particles with radius R 2 , this interface is significantly distorted and its mean interface curvature, or roughness, is quantified by r. In the "pinned state" limit, the system is assumed FIGURE 2 | Sketch of the impact of Zener pinnning on grain interface roughness (after Bercovici and Ricard, 2012;Mulyukova and Bercovici, 2019). Shown is a sketch of a grain surrounded by secondary particles, with a zoom on the interface around one of the secondary particles.
to be in a state where grain size evolution is controlled by the evolution of interface roughness r. Grain damage is therefore solely controlled by grain interface damage in the pinned state limit. Although we focus here on grain interface damage, we will use the more general term grain damage throughout the manuscript. If the secondary phase is dispersed (which we assume here), this effect results in interface roughness depending on the grain size R i of each phase as (see Bercovici and Ricard, 2012, Appendix G): Frontiers in Earth Science | www.frontiersin.org where c = 3λ 4 160λ 2 = 3 160 e 6σ 2 is a factor which depends on the half width σ of the lognormal grain size distributions, which are assumed to have the same shape (see Bercovici and Ricard, 2012, Appendix F.4). Consequently, the grain sizes of both phases depend on r as (e.g., Bercovici and Ricard, 2012;Bercovici et al., 2015;Mulyukova and Bercovici, 2018): where φ i are the respective volume fractions of the two phases. With this definition, the mean grain size R can be computed as: Assuming a pyroxene content of 40 % and a half width of the lognormal grain size distributions of 0.8 (e.g., Bercovici and Ricard, 2013), the mean grain size is given by R ≈ π 2 r. Additionally assuming that both mineral phases have approximately the same rheological properties, we can-in analogy to the transition grain size R t (Equation 5)-define a transition interface roughness r t , which only differs by a constant prefactor π 2 m from the transition grain size R t : Therefore, the composite rheology can be re-written in the same manner as Equation (6), only replacing R R t with r r t : In the "pinned state", the grain size of either mineral phase is assumed to be solely determined by the interface curvature r. The interface curvature r evolves with time according to (Bercovici and Ricard, 2012): where q is an effective grain growth exponent, γ the interfacial energy, η = 3φ 1 φ 2 is a parameter that depends on the volume fraction φ i of both phases, G I is the temperature-dependent interface coarsening rate given by (e.g., Mulyukova and Bercovici, 2017): The deformational work rate , given by the product ofε and τ , is: Both the work done in dislocation and diffusion creep contribute to grain damage and hence interface roughness. Last, the partitioning factor λ I determines the relative amount of deformational work that is used for grain damage. This factor is relatively unconstrained and is most likely dependent on temperature (Rozel et al., 2011;Mulyukova and Bercovici, 2017), phase dispersion (Bercovici and Ricard, 2016;Bercovici and Skemer, 2017) as well as dislocation density (Chrysochoos and Belmahjoub, 1992;Holtzman et al., 2018). Given the missing constraints on most of these dependencies, we neglect these effects for simplicity. In steady state, when grain damage and recovery balance each other such that dr/dt = 0, we can define a steady state roughness r s that has to fulfill the following equation: This equation is not analytically solvable for r s for the general case. However, we can consider two particular cases, one where deformation is dominated by dislocation creep and another where diffusion creep is dominant. In the dislocation-dominated case (r s ≫ r t ), the steady state interface roughness is given by: whereas in the diffusion-dominated case (r s ≪ r t ), it is: Using Equations (16) and (17), we find that r s,dis and r s,dif are related by:

Thermal Evolution
Assuming adiabatic conditions, the temperature evolution in the necking region is given by: Only the deformational work that is not used for grain damage results in heating, which is controlled by the coefficient 1 − λ I (this coefficient is also frequently called the Taylor-Quinney coefficient). We neglect the surface energy that is released during interface coarsening as its contribution is likely negligible at typical interface roughnesses (Thielmann et al., 2015;Thielmann, 2018).

Reference State and Non-dimensionalization
The equations derived above contain a large number of material parameters. To identify potential trade-offs between different material parameters and to reduce their number it is useful to (i) define them in terms of a reference state and (ii) nondimensionalize them using appropriate scales.
We use the initial state as the reference state with initial stress τ 0 , initial temperature T 0 and initial grain interface transition roughness r t,0 . First, we express the temperature T as the sum of the initial temperature T 0 and a temperature perturbationT: Using this, an Arrhenius term of the form e can be rewritten as: Second, with the above form of the Arrhenius term the flow law given in Equation (11) can be rewritten in terms of the reference state: where we used Equation (10). Rewriting the rheology also allows to rewrite the grain size evolution Equation (12) in terms of the reference state: where, using Equation (16), we can identify: It is useful to define a normalized grain interface roughnessr = r r t,0 and introduce dr dt = r t,0 dr dt , to obtain an equation for the evolution of the normalized interface roughness (which does not have any units): Additionally, we can rewrite the relation between the steady state dislocation creep grain size and the steady state diffusion creep grain size given by Equation (18) Using Equation (20) and the definition of the normalized grain interface roughness, the temperature Equation (19) can be written as: Equating (1) and (22) yields: This leaves us with three coupled ordinary differential equations for the three unknowns r, T, and D. While this notation may seem more complicated, it is significantly simplified when nondimensionalizing Equations (26)-(29) with a set of characteristic scales for length, temperature, stress, and time, respectively, which are given by the initial slab thickness, initial slab temperature, initial slab stress, and the initial strain rate at the transition between diffusion and dislocation creep: Frontiers in Earth Science | www.frontiersin.org With the chosen characteristic temperature scale, the Arrhenius term can be rewritten as: where T ′ 0 = Q dis /T 0 . Replacing τ τ 0 = D 0 D and nondimensionalizing the respective quantities in Equations (24), (28), and (29) with the respective scales then yields for better readability. Asr is already a non-dimensional quantity, it does not have to be additionally non-dimensionalized. The superscript ′ indicates a non-dimensionalization with the characteristic scales in Equation (30), whereasˆindicates a normalization with the initial transition roughness r t,0 . We use two characteristic length scales because our model involves both the scale of the lithosphere thickness and the scale of a mineral grain, which can be more than ten orders of magnitude different. The two length scales are useful because they provide a more intuitive meaning of both ther and D ′ . Additionally, the numerical treatment of this system of equations is made much easier, as values ofr and D ′ do not attain extreme values prone to roundoff errors. The following non-dimensional numbers can be identified: as well as the three non-dimensional exponents n, m, and q.

Nondimensional Parameters
It is worthwhile to discuss the meaning of the different nondimensional parameters introduced in Equation (35) as well as the values they may take for Earth-like conditions. Due to their definition, the initial values for temperatureT ′ 0 and slab thickness D ′ 0 are 0 and 1, respectively. The meaning of the nondimensional initial temperature and interface curvature T ′ 0 andr 0 are straightforward, with values ofr 0 > 1 indicating that the initial deformation is mostly accommodated by dislocation creep and forr 0 < 1 mostly accommodated by diffusion creep. The values of δQ g and δQ dif indicate the respective temperature sensitivity of grain growth and diffusion creep with respect to dislocation creep. Values larger than 1 denote that dislocation creep is more sensitive to a temperature change than either grain growth or diffusion creep. This may, for example, be important during the necking process, as an increase in temperature due to shear heating may result in a switch from diffusion to dislocation creep. The value forr s,dis,0 is the normalized steady state interface curvature in dislocation creep at the initial conditions, which is closely linked to the normalized steady state interface curvature in dislocation creep (see Equation 27). Depending on the initial grain sizer 0 , the interface curvature will either increase (r 0 < r s,dis,0 ) or decrease (r 0 >r s,dis,0 ) to attain this steady state value. As for the initial interface curvaturer 0 , values ofr s,dis,0 > 1 indicate that the steady state interface curvature will favor dislocation creep. However, this is only true at the initial conditions. During the necking process, stress (being inversely proportional to the slab thickness D), temperature and interface curvature will change and thus also the transition curvature and steady state curvature. Whether deformation in these cases will be governed by dislocation or diffusion creep then depends on the values of n, m, and δQ dif .
Finally, the damage numbers I and S determine the importance and speed of grain damage and temperature induced damage, respectively. Large values of the respective parameter denote fast and efficient weakening. Both the grain damage number I and the thermal damage number S can be considered as expressions of the available energies in the system. The grain damage number I can be considered as the ratio between the specific deformation energy available for grain damage, given by the term λ I τ 0 and the energy stored in the grain interfaces, which is given by the term γ η r t,0 (Bercovici and Ricard, 2012). Hence, when the available specific deformation energy is also much larger than the energy stored in the interface (resulting in a large value of I), grain damage is more efficient. The same applies to the thermal damage number S, which can be considered as the ratio between the specific deformation energy available for thermal damage, (1 − λ I )τ 0 , and the specific thermal energy, ρcT 0 . Due to our choice of non-dimensionalization, this ratio is then additionally scaled by the factor Q dis /T 0 , also known as the Arrhenius number (e.g., Peters et al., 2016), which is in itself a measure of the ratio between activation energy and thermal energy. A large value of S thus requires that (1) the system is highly temperature dependent (represented by large values of the Arrhenius number) and (2) the available deformation energy is much larger than the thermal energy.
To determine Earth-like ranges of S, I,r s,dis,0 , and T ′ 0 , we assume slab lengths H of 100-500 km, mean slab temperatures T 0 of 900-1,200 K and an interface partitioning factor λ I within the range of 10 −6 − 10 −3 (e.g., Bercovici and Ricard, 2016;Mulyukova and Bercovici, 2017). There are at present no experimental constraints on the value of λ I , thus its value is usually assumed based on indirect evidence. With the material parameters given in Table 1, we then compute the respective nondimensional numbers for 10 5 random combinations of the input parameters to estimate their range. The resulting histograms for I and S,r s,dis,0 , and T ′ 0 are shown in Figure 3. For the chosen parameters, I ranges between 1 and 1,000, S shows much less variation and varies between 1 and 10.r s,dis,0 and T ′ 0 vary in the range of 10 −3 -10 −2 and 50-70, respectively (see also

Comparison to Previous Studies
Before applying these models to Earth, we will first have a general look at the system of equations given by (32)-(34). Unless otherwise noted specifically, the primes indicating nondimensional values are dropped for the remainder of the paper for better readability. We can identify several particular cases that have been studied in the literature, the first being the analytical necking model by Schmalholz (2011). In this model, only dislocation (powerlaw) creep was considered, without taking into account either shear heating or grain size evolution. In this case, only (34) has to be considered, where the grain size dependent diffusion creep term is dropped. This is equivalent to setting I = S = 0 andr 0 = ∞. AsT = 0, the remaining Arrhenius term equals one (34), which results in: The above equation corresponds to the equation derived in Schmalholz (2011). For n = 1, the detachment time t d , that is the duration until the slab thickness becomes zero, is given by the characteristic time scale t c . For n > 1, the detachment time then depends on the powerlaw exponent as t d = 1 n . Additionally considering shear heating (by setting S > 0) then yields the following system of equations: which corresponds to the system investigated by Duretz et al. (2012) (the thermal damage number S is equal to A in their Equation 6). Finally, Bercovici et al. (2015) investigated the effect of grain damage and a composite dislocation/diffusion creep rheology on slab detachment, but not shear heating (S = 0, I > 0, 0 <r 0 < ∞). In that case, (32)-(34) reduce to (due to the constant temperature, the respective Arrhenius terms equal 1): This system of equations is equal to the one used in Bercovici et al. (2015), except that they additionally included the effect of the viscous resistance of the mantle to the sinking slab. The damage number D used in their paper is equal to I used here, whereas thê r s,dis,0 used here equals their C qI . Both the studies of Duretz et al. (2012) and Bercovici et al. (2015) concluded that the respective weakening mechanism resulted in a significant reduction of the detachment time.

RESULTS
To determine the impact of different parameters on the necking process, we assessed their effect on the detachment time of the slab, t d . We follow a two-fold approach which consists of (1) analytical considerations of several simplified cases and (2) a numerical Monte-Carlo approach to test the prediction potential of the analytical approximations. We solved the set of ordinary differential Equations (32)-(34) numerically using a stiff ordinary differential equation solver of the software Matlab. The numerical code can be found in the online repository on figshare (Thielmann and Schmalholz, 2020a).
For the two non-damage cases, detachment times are longest. However, once a single or both damage mechanisms are activated, detachment times are significantly reduced. The reduction is largest when both interface and thermal damage are active. For all cases where damage is activated, we observe a very rapid, almost instantaneous detachment phase. As can be seen in Figure 4C, the interface curvature is significantly reduced if grain damage is active and rapidly approaches very small values. On the contrary, if thermal damage is active, temperatures significantly rise during the detachment  phase (Figure 4D), which in some cases results in a switch from diffusion to dislocation creep at late detachment stages ( Figure 4B).
While these examples qualitatively show the impact of both damage mechanisms on detachment times, we aim at quantifying these effects. We will first develop analytical approximations for t d for several extreme cases before turning to numerical modeling. Even though Earth-like values for different nondimensional parameters are restricted (see Figure 3), we will also investigate regimes that go beyond these parameter ranges.

Analytical Considerations
3.2.1. Low Damage Numbers (I ≪ 1, S ≪ 1) (32) and (33) can be neglected and the evolution of slab thickness is solely governed by Equation (34). Due to the low damage numbers, the interface curvaturer as well as the temperature perturbationT remain at their initial valuesr 0 and 0. Depending onr 0 , deformation will either be dominated by diffusion or dislocation creep and the evolution of the slab thickness D is governed by:

At low damage numbers, both Equations
which can be integrated to yield: from which we can compute the detachment time t d at which D = 0:

Shear Heating Dominated Weakening (I ≪ 1)
In the case of shear heating dominated weakening, only Equation (32) can be neglected. The system of equations is then given by Equations (33) and (34). As in section 3.2.1, we will also distinguish between a diffusion and a dislocation dominated case. To understand the necking evolution in this case, we now examine it at the initial stages where temperature changes can be assumed to be small. The evolution of slab thickness in each case is then given by Equation (42). The evolution of temperature (Equation 33) reads as: As we are only considering initial stages, we can additionally approximate the term e T using a Taylor series expansion around T = 0: where we neglect higher order terms. Rewriting Equation (44) and inserting Equation (42) then results in: which can be integrated using separation of variables to yield: As can be seen in Figure 4,T(t) goes to infinity as D approaches zero when thermal damage is important. For the expressions of T(t) given above, the critical time t crit at which T → ∞ is therefore an indicator of the detachment time t d . This occurs when:r which yields the following expressions for t crit : For small values of S, these expressions are equal to the lowdamage case derived in section 3.2.1, thus reproducing the prediction of Schmalholz (2011)

Grain Damage Dominated Weakening (S ≪ 1)
If grain damage dominates weakening during the detachment process, the necking process can be described by Equations (32) and (34). When deformation is governed by dislocation creep, grain damage does not affect the necking process. Necking within this regime can thus be described using Equation (42). In this case (and assuming that grain growth is negligible), grain size evolution is described by which can be integrated to yield: Using this equation, we can now determine the time t t at which the rheological transition between diffusion and dislocation creep is reached, which is the case when the dislocation and diffusion creep yield equal strain rates. If grain size evolution is sufficiently fast, this will be the case whenr ≈ 1. Inserting this value forr(t) in the above equation and solving for t results in: Frontiers in Earth Science | www.frontiersin.org TABLE 2 | Nondimensional model parameters and parameter ranges computed with the material parameters given in Table 1 and for slab lengths H of 100-500 km, average slab temperatures T of 900-1,200 K, pressures of 3-9 GPa, and partitioning factors λ I ranging from 10 −6 to 10 −3 .

Dislocation creep Diffusion creep
a In the dislocation dominated case deformation is initially controlled by dislocation creep, but eventually switches to diffusion creep.
Once this time is reached, the importance of diffusion creep will increase. To assess the necking behavior in this necking phase, we again consider its initial stage, where we assume that slab thickness D does not change significantly. In this case, only Equation (32) has to be considered. Atr 0 = 1, dislocation and diffusion creep yield the same strain rate so that Equation (32) can be rewritten as (additionally neglecting grain growth): which, upon integration, results in the following expression for the evolution ofr:r For this equation, we can now define a critical time t crit at whicĥ r approaches zero: The total detachment time for the case of grain damage dominated weakening with the initial deformation dominated by dislocation creep is therefore the sum of t t and t crit : If, on the other hand, the initial interface curvaturer 0 < 1, deformation at the initial stages will directly be governed by diffusion creep. In this case, Equation (32) can be written as (again neglecting grain growth): which, can again be integrated and used to compute a critical time t t , which only slightly differs from the one derived above: This critical time then represents the detachment time of the slab in the case of grain damage dominated weakening with the initial deformation dominated by diffusion creep. In sections 3.2.1-3.2.3, we have derived a set of analytical predictions for the detachment time t d for different deformation regimes using simple approximations. These predictions are summarized in Table 3.

Numerical Monte Carlo Study
As the analytical approximations derived above are based on several simplification, it is important to determine to which extent they are capable of predicting the detachment time t d . To further quantify the impact of different material parameters on the detachment time t d , we now turn to numerical modeling. Here we define t d as the time it takes for the neck to reach a thickness of 1% of its original value. We determined the detachment time of 2 × 10 5 different simulations where we randomly varied S, I,r s,dis,0 and T 0 within the ranges given in Table 2. Random numbers were drawn from uniform distributions using the Mersenne Twister pseudorandom number generator. As the number of nondimensional parameters is still large, we used some simplifying assumptions: (1) Based on existing flow laws and two-phase growth laws for olivine (Hirth and Kohlstedt, 2003;Hiraga et al., 2010;Bercovici and Ricard, 2012), we assume fixed values for n, m, and q of 3.5, 3, and 4, respectively, (2) we assume that the activation energies for diffusion creep and interface coarsening are equal (δQ dif = δQ g = 0.566, e.g., Bercovici et al., 2019), and (3) at the initial conditions, grain size is large enough so that necking is dominated by dislocation creep (r 0 = 2).
The results of this Monte Carlo simulation are shown in Figure 5, where we plot the detachment time of each simulation with respect to the grain and thermal damage numbers I and S as well as the initial dislocation dominated steady state interface roughnessr s,dis,0 . Although the results are quite scattered, some potential underlying trends can be observed: (1) for small values of I < 1 and S < 1, necking times seem to be insensitive to both values, (2) for values of I < 1 and S > 1, necking times seem to primarily depend on S (Figure 5A), (3) I > 1, andr s,dis,0 < 1 significantly reduce detachment times (Figure 5B), albeit not in all simulations and (4) ifr s,dis,0 > 1, necking times seem to be primarily affected by S. To further explore these dependencies and to verify the previously derived analytical solutions, we will now consider three separate regimes: (1) a low grain damage regime, (2) a high grain damage, low thermal damage regime, and (3) a regime where both damage mechanisms play a role. As the initial grain sizer 0 is fixed at a value of 2, thus favoring dislocation creep, our analysis will be concerned with the analytical solutions derived for cases where initial deformation is governed by dislocation creep.

Low Grain Damage (I < 1)
As derived in section 3.2.2, detachment times in this regime are independent of I. In Figure 6A, we show the computed detachment times for all simulations with I < 1 as a function of S. In addition, the theoretical prediction of t d using Equation (50) withr 0 > 1 is shown in gray. We see that, despite the simplifying assumptions made in deriving Equation (50), the analytical solution predicts the numerical result well and provides an accurate upper bound for t d . Only at the transition from small to large values of S, numerical results deviate from the analytical prediction by up to 35% (Figure 6B). Deviations are largest for simulations where a large value of I is employed, indicating that in those cases grain damage already has an impact on the necking process. The larger deviations at values of S around 1 reflect the approximation errors introduced while deriving the respective approximations. In general, at values of S < 10 −1 , detachment times are not influenced by thermal damage and thus are only dependent on the value of the stress exponent n, while at large values of S > 10 1.5 , detachment times scale with 1/S.

Large Grain Damage and Low Thermal Damage
In this regime, detachment times should be largely independent of S, but rather depend on I. Due to the chosen value of r 0 = 2, deformation in this regime is always dominated by dislocation creep at the initial stages, but may switch to diffusion creep depending on the value ofr s,dis,0 . In Figure 7A, the detachment times for all simulations are shown, with their color denoting the employed value of r s,dis,0 . The two solid black lines show the predictions for t d given by Equations (43) and (57). The first equation predicts the detachment time in the case of dislocation creep dominated deformation (which is grain size insensitive and thus independent of I), while the second equation predicts the detachment time in the case of deformation governed by dislocation creep at the initial stages, but switching to diffusion creep at later stages, thus essentially being governed by diffusion creep.
The numerical results are quite scattered, but exhibit two main characteristics: at values of ∼r s,dis,0 > 10 0.5 detachment times are independent of I and are well predicted by Equation (50) for dislocation creep dominated deformation. This is illustrated in Figure 7B, where we plotted the ratio t d,num /t d,ana between the numerically computed and the predicted detachment time as a function of S. Here, colors denote the value ofr s,dis,0 for each simulation. As expected, the analytical prediction for t d in this dislocation creep dominated regime fits best for small values of S, while detachment times deviate from the prediction when S approaches 1. For values of ∼ r s,dis,0 < 10 −0.5 (diffusion dominated deformation), detachment times are well predicted by Equation (57) at values of I > 10 1 and deviate from the analytical prediction as I approaches 1 (Figure 7C). The analytical solutions thus predict detachment times well, except for the transition from grain damage dominated weakening to structural weakening (at I = 1) or the transition between dislocation creep dominated deformation to diffusion creep dominated deformation (at r s,dis,0 = 1).

Large Interface and Thermal Damage (I > 1, S > 1)
This regime is the most interesting one, as both interface and thermal damage will influence the necking process. Additionally, as Figure 3 indicates, detachment processes on Earth are most likely to lie within this regime. In Figure 8A, the detachment times of all simulations within this subset are shown, with colors denoting the value ofr s,dis,0 . Although detachment times are very scattered, a clear pattern can be seen, with diffusion creep dominated simulations (r s,dis,0 < 10 −0.5 ) representing a lower bound to detachment times. From this subset, we therefore extracted two additional subsets, one of which contained dislocation creep dominated simulations (r s,dis,0 > 10 0.5 ), the other one containing diffusion creep dominated simulations (r s,dis,0 < 10 −0.5 ). In the case of dislocation creep dominated simulations, Equation (50) should predict the detachment time, since grain damage does not affect the rheology. In Figure 8B, the ratio t d,num /t d,ana between the numerically computed detachment times and the analytical prediction is plotted vs. thermal damage parameter S. It can be seen that the analytical prediction predicts the detachment time very well, with deviations becoming larger as S approaches 1. Similarly, for diffusion creep dominated simulations, Equation (57) provides a good prediction for t d . The comparison between this prediction and the numerically computed results is shown in Figure 8C, where we plot the ratio t d,num /t d,ana between the numerically computed detachment times and the analytical prediction vs. the thermal damage parameter I. Colors denote the ratio between the grain damage parameter I and the thermal damage parameter S. We can see that if I is much larger than S, detachment times are best predicted, with numerically computed detachment times becoming significantly decreased for (1) similar values of I and S  and (2) for I approaching 1. As could be already seen in Figure 4, this highlights the cumulative effect of grain and thermal damage.
In summary, the comparison between numerical model results and analytically derived detachment time predictions shows that the analytical predictions are capable to reasonably predict numerical model behavior for either dislocation creep dominated cases (r s,dis,0 > 10 0.5 ) or diffusion creep dominated cases (r s,dis,0 < 10 −0.5 ), when thermal and grain damage parameters attain extreme values. At the transition between dislocation creep and diffusion creep dominated regimes as well as when thermal and grain damage parameters take values around 1, numerically computed detachment times deviate strongest from the analytical predictions. However, these largest deviations typically never exceed 50%.

Impact of Thermal Diffusion
In the analysis presented above we have so far neglected thermal diffusion for simplicity, thus omitting the diffusion term in Equation (19) and its non-dimensional counterpart in Equation (34). During slab detachment, the cold slab will be heated due to thermal diffusion by the surrounding hotter mantle. This heating will weaken the slab and, hence, thermal diffusion may have a significant effect on the necking process (e.g., Gerya et al., 2004). However, thermal diffusion is only important, if the time scale of diffusive heating is not significantly longer than the time scale of slab detachment. To assess the time scale and importance of thermal diffusion during slab detachment and how it may affect detachment times, we here employ a simple analysis. First, we approximate the diffusion term κ ∂ 2 T ∂x 2 using a simple finite difference approximation: where T m is the ambient mantle temperature and κ the thermal diffusivity. This term has units of K/s, which implies that it  has to be multiplied with t c /T c for non-dimensionalization. In non-dimensional form, this approximated diffusion term can be written as: where t diff = D 2 0 /κ is the diffusion time scale and t c is the characteristic time scale chosen for non-dimensionalization, here the dislocation creep strain rate at the initial conditions (see Equation 30). δT m = T m /T 0 is the ratio between the ambient mantle temperature and the slab temperature. Considering temperature diffusion thus introduces the two additional nondimensional numbers t c /t diff and δT m .
For Earth-like parameters, the range of δT m is relatively small and lies between ∼ 1.3-1.8, while the range of t c /t diff is significantly larger and may range between ∼ 10 −4 and10 1 . Values for t diff are limited by the range of Earth-like values for thermal diffusivity and slab thickness. Both these properties cannot vary orders of magnitude, thus the variations in t d are a result of variations in t c . This characteristic time scale is determined by the nonlinear rheology (see Equation 30) and yields large values for cold and short slabs and small values for long and warm slabs.
A large value of t c /t diff implies that the slab may be significantly heated due to thermal diffusion. Due to this heating, the slab will be weakened and detachment times will be reduced. On Earth, thermal diffusion therefore may become important for cold and short slabs with long detachment times. This is also reflected in a value of 200 Ma for t diff (assuming a slab with thickness 80 km and a thermal diffusivity of 10 −6 m 2 s −1 ).

Application to Earth
To determine the impact of grain size reduction and shear heating on slab detachment on Earth, we now compute detachment times for Earth-like parameter ranges. As shown in Figure 3, parameter ranges for I and S include values that span from ∼ 10 −1 to 10 2 and ∼ 10 −1 to 10 1 , respectively, thus putting Earth-like detachment processes right within the transitional regime, where both grain and thermal damage are important. Employing the Earth-like parameters ( Table 1) and assuming a range of slab lengths H and temperatures T slab , we computed the resulting detachment times numerically for these transitional regimes. Additionally, as the value of the partitioning factor λ I is highly unconstrained, we varied its value from 10 −6 to 10 −4 . The results show that for the chosen parameters, the detachment times range between 10,000 and more than 10 billion years (Figures 9A-C), with the extremely long detachment times being a result of cold slab temperatures in short slabs. When thermal diffusion is considered, these extremely long detachment times are reduced to values less than ∼100 million years for cold and short slabs (Figures 9A-C). As the value of λ I is increased, detachment times decrease, but the effect is limited to very cold slab temperatures and short stalled slab lengths.
To better illustrate the effect of thermal diffusion and variations in the partitioning factor λ I , we show a comparison of six different cases in Figure 10, where we assumed a slab temperature of 1,100 K, a slab thickness of 80 km and an initial interface curvature of 1 cm. Increasing λ I from 10 −6 to 10 −4 results in more than a 9-fold reduction of the detachment time ( Figure 10A, thus highlighting the importance of this parameter. The effect of thermal diffusion can only be seen in simulations where λ I is small and hence detachment times large. When looking at the evolution of the deformation mechanism ratio (Figure 10B), which indicates dominant dislocation creep for > 1, it becomes apparent that λ I is crucial in determining the dominant deformation mechanism due to its strong impact on interface curvature ( Figure 10C). The significant reduction in detachment time for the cases where λ I = 10 −4 indicates that the initial stages of slab detachment may be governed by grain damage. Indeed, interface curvature experiences a rapid decrease from 10 −2 m to about 5·10 −4 m at initial stages where slab thickness is reduced from 80 to 70 km (yellow line in Figure 10C). On the other hand, temperatures do not increase significantly at the initial necking stages. However, when a slab thickness of ∼30 km is reached, temperatures start to increase significantly. This implies that while shear heating may not govern initial stages of slab detachment, it becomes important at later stages.
Seismic imaging has revealed several regions on Earth where high-velocity bodies have been identified in the upper mantle, among them the Northwestern US (Schmandt and Humphreys, 2011;Wang et al., 2013;Jiang et al., 2018), the Hindu-Kush (Lister et al., 2008;Kufner et al., 2017), the Mediterranean (Wortel and Spakman, 2000), the Carpathians (Sperner et al., 2001;Ismail-Zadeh et al., 2012), and the European Alps (Lippitsch et al., 2003;Kästle et al., 2020). These high velocity bodies are commonly interpreted as cold lithospheric material. However, their origin is often debated, as they could be interpreted as either subducted slabs or lithospheric drips (e.g., Sperner et al., 2001;Zandt et al., 2004;Lister et al., 2008;Lorinczi and Houseman, 2009;Kufner et al., 2017;Molnar and Bendick, 2019). Observations linked to slab detachment processes in these regions suggest slab detachment events occurring within a few million years after continental collision. Only in cases for which short slabs remain after a previous detachment event, longer "stalling times" in excess of tens of millions of years have been suggested (e.g., Jiang et al., 2018). This is in accordance with our results (see Figure 9), which also only predict detachment times in excess of 10 Ma for sufficiently short or cold slabs.
Despite the initial detachment process being mainly controlled by grain size, shear heating may have a strong impact on the "terminal stage of subduction" (Kufner et al., 2017), when elevated stresses result in increased viscous dissipation and thus heating. These elevated stresses together may not only trigger brittle failure and thus earthquakes, but may also trigger ductile instabilities at larger depths that have been proposed to be the cause for intermediate-depth earthquakes (Kelemen and Hirth, 2007;John et al., 2009;Deseta et al., 2014;Thielmann et al., 2015;Kufner et al., 2017;Ohuchi et al., 2017;Thielmann, 2018;Zhan, 2019). To further investigate this issue, it will be necessary to additionally take into account plastic failure, as the large stresses in the slab at final detachment stages may also result in brittle deformation.

Model Limitations
Owing to the simplicity of our model, computed detachment times should not be taken at face value, but rather as a first-order constraint. Due to the simplifying assumptions, our model most likely underestimates the time scales of the detachment process, as we do not account for additional forces such as (1) the viscous resistance of the mantle to the sinking slab (2) resistance of the asthenosphere to flow into the necking region of the slab and (3) the viscous resistance of the lower mantle when a slab tip reaches the bottom of the transition zone during detachment (e.g., Burkett and Gurnis, 2013). These effects could potentially slow down slab detachment (Bercovici et al., 2015). However, it is not clear to which extent these viscous forces would contribute to a delay in slab detachment. This is partly due to the nonlinear rheology of the mantle itself, which may be affected by dislocation FIGURE 9 | Application to Earth. Detachment times t d for a range of slab lengths H and mean slab temperatures T slab . Detachment times were computed using the parameters listed in Table 1 creep and/or grain size reduction. In both cases, the viscosity of the mantle would be significantly reduced, thus reducing its resistance to the sinking slab (see also discussion in Bercovici et al., 2019).
The force balance assumed in our model may also be modified by other factors, such as the presence of positively buoyant crustal material at the top of the slab or negatively buoyant material due to the phase transition to eclogite at greater depths. Also, if the slab were not stalled, but rather subducting at a certain speed, the downward buoyancy force would be governed by their interplay with a number of other forces such as the viscous coupling of the plate at the surface to the asthenosphere and the coupling between subducting and overriding plate (e.g., Forsyth, 1975;Bercovici et al., 2019). Nevertheless, the fundamental interplay between the necking instability and the two damage processes will remain unchanged, whereas absolute time scales may change due to different stress levels in the slab. Combining the present model with the recent model by Ribe and Xu (2020) could potentially account for several of these effects. One uncertainty for all mathematical models arises from the uncertainty of the different material parameters. Rheological parameters for dislocation and diffusion creep are to some extent constrained by laboratory experiments, but the extrapolation to tectonic time scales and the application to kilometer-scale heterogeneous rock units generates uncertainties. One of the largest parameter uncertainty of our model arises from the lack of constraints on the partitioning factor λ I , which was assumed to be constant in this study. However, there is evidence that this factor depends on several other variables such as grain size, temperature, dislocation density, and phase mixing (see section 2.3). Some of these effects could potentially mitigate the effect of grain size reduction, as for example λ I is thought to decrease at small grain sizes (Bercovici and Skemer, 2017) and larger temperatures . Due to the strong effect of λ I on model results, it is indispensable to obtain better constraints on this parameter, which could prove to be a challenge for experimental geoscientists. Additional constraints on the interplay between necking, grain damage and shear heating can also be obtained from structural geology (e.g., Linckens et al., 2015).
Despite the uncertainties of our model, it provides the means to analyze the fundamental physical relationship between necking, grain damage, and thermal weakening, all of which have a first-order impact on slab detachment. It also shows the importance of different material parameters (or combinations thereof), thus indicating which parameters are necessary to constrain accurately and which parameters are less important.
This not only has implications for slab detachment, but also for other necking processes (e.g., Brune et al., 2016), where continental plate breakup, following lithosphere necking (e.g., Chenin et al., 2018), may be significantly accelerated due to the action of ductile weakening.

CONCLUSIONS
We have developed a fundamental mathematical model of slab detachment, considering structural weakening due to necking, thermal weakening due to shear heating (thermal damage) and microstructural weakening due to grain size reduction (grain damage; here represented by changes in interface curvature between two phases). Our model extends previously developed models (Schmalholz, 2011;Duretz et al., 2012;Bercovici et al., 2015) and is, to the best of our knowledge, the first model to link these three coupled weakening processes in a system of three nonlinear differential equations. The derived system of equations allows investigating the fundamental interplay of the three weakening mechanisms and to determine the conditions for which one, two or all three mechanisms are important. With dimensional analysis we determined the controlling dimensionless parameters of this system and we derived analytical solutions for the duration of slab detachment for several end-member regimes. We verified the controlling parameters and analytical solutions with a systematic numerical Monte-Carlo study. The analytical solutions for the detachment time are useful for understanding the fundamental control of the three weakening mechanisms and the involved parameters on slab detachment. Numerical solutions show that the analytical solutions are inaccurate for transitional regimes where more than one weakening mechanism is important. However, in these transitional regimes the analytical solutions typically do not deviate by more than 50% from the numerical solution, thus lying well within the range originating from uncertainties in the applied material parameters.
Both grain damage, due to grain size reduction, and thermal damage, due to shear heating, result in rapid detachment, showing that both damage mechanisms have a first order impact on the slab detachment process. Models where both damage mechanisms are equally active exhibit the shortest detachment times as a result of the positive feedback between necking, grain size reduction and shear heating. This positive feedback shows that both damage mechanisms should be considered in models of lithosphere necking.
For lithospheric-scale conditions on Earth, all three weakening mechanisms are presumably important during slab detachment. We find that for these conditions, grain damage and thermal weakening act in a collaborative manner, thus minimizing the detachment time. Grain damage is likely dominant during the initial stages of slab detachment, whereas thermal weakening becomes more important and presumably dominant during the later stages. Our results suggest that the typical duration of slab detachment is between 1 and 10 million years. Detachment times in excess of 10 Million years are only predicted for short (<≈200 km) and cold slabs, which is in broad agreement with results from other theoretical studies and indirect estimates based on geological and geophysical observations.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author. The code to create the datasets and figures used in this study are available on figshare (Thielmann and Schmalholz, 2020a,b).

AUTHOR CONTRIBUTIONS
MT and SS designed the study, analyzed the data, and wrote the manuscript. MT derived the governing equations and performed simulations. All authors contributed to the article and approved the submitted version.

FUNDING
MT was funded by the Bayerisches Geoinstitut Visitor's program and SS by the University of Lausanne.