Relamination of mafic subducting crust throughout Earth's history

Earth has likely cooled by several hundred degrees over its history, which has probably affected subduction dynamics and associated magmatism. Today, the process of compositional buoyancy driven upwelling, and subsequent underplating, of subducted materials (commonly referred to as “relamination”) is thought to play a role in the formation of continental crust. Given that Archean continental crust formation is best explained by the involvement of mafic material, we investigate the feasibility of mafic crust relamination under a wide range of conditions applicable to modern and early Earth subduction zones, to assess if such a process might have been viable in an early Earth setting. 
 
Our numerical parametric study illustrates that the hotter, thicker-crust conditions of the early Earth favour the upward relamination of mafic subducting crust. The amount of relaminating subducting crust is observed to vary significantly, with subduction convergence rate having the strongest control on the volume of relaminated material. Indeed, removal of the entire mafic crust from the subducting slab is possible for slow subduction (∼2 cm/yr) under Archean conditions. We also observe great variability in the depth at which this separation occurs (80–120 km), with events corresponding to shallower detachment being more voluminous, and that relaminating material has to remain metastably buoyant until this separation depth, which is supported by geological, geophysical and geodynamical observations. Furthermore, this relamination behaviour is commonly episodic with a typical repeat time of approximately 10 Myrs, similar to timescales of episodicity observed in the Archean rock record. We demonstrate that this relamination process can result in the heating of considerable quantities of mafic material (to temperatures in excess of 900 °C), which is then emplaced below the over-riding lithosphere. As such, our results have implications for Archean subduction zone magmatism, for continental crust formation in the early Earth, and provide a novel explanation for the secular evolution of continental crust.


Introduction
The mechanisms that have led to the formation of the felsic continental crust (CC) from the mantle throughout Earth's history are still poorly understood, but the presence of subduction systems seems to be an essential component. Indeed, the composition of today's subduction related magmas shares many similarities with the CC (Rudnick, 1995;Jagoutz and Kelemen, 2015). This also applies to the most ancient CC (>2.5 Ga), however, there are systematic differences (Taylor and McLennan, 1995). In particular, certain geochemical features common of this early crust are widely recognized to be derived from mafic material at high pressure (Defant and Drummond, 1990;Foley et al., 2002;Alonso-Perez et al., 2009), but the particular processes that oper-ever, when looking at Archean granulites this figure falls to 6%  in accordance with the fact that the composition of Archean CC is best explained by the involvement of mafic, not just sedimentary, crustal material. Secular changes in crustal composition suggest that CC formation differed from the modern day mechanism, but exactly how is debated (Condie, 2005;Van Hunen and Moyen, 2012). Whether subduction even played a role in the Earth's earliest history has been thrown into question (Davies, 1992;Sizova et al., 2010), although there is structural evidence for Archean subduction (Calvert et al., 1995;Stevens and Moyen, 2007). If subduction did occur in the early Earth, conditions within the subduction system are expected to have differed systematically from conditions today. In particular, mantle temperatures were probably higher (Herzberg et al., 2010), and oceanic crust was likely thicker as a result (Hoffman and Ranalli, 1988;Abbott et al., 1994). This oceanic crust may have differed compositionally from today's MORB (Herzberg et al., 2010) and, as such, may have had different physical properties. Moreover, it is debated whether subduction in the Archean was slower (Korenaga, 2006), faster (van Hunen and van den Berg, 2008), or comparable to today (Sleep, 1992) and perhaps occurred at a shallower angle (Abbott et al., 1994) although this point is also debated (van Hunen et al., 2004). All these differences may have affected the dynamics of the subduction interface and, in particular, allowed for the relamination of more of the subducting crust.
This leads to the hypothesis: if relamination of subducted felsic materials occurs in modern subduction zones then perhaps, under different, early Earth conditions, such a process may have been viable not only for subducted sediments, but for more of the mafic part of the subducting crust as well, thus producing the characteristics of Archean CC. In this paper we aim to constrain the conditions required for the relamination of the mafic subducting crust by using a thermo-mechanical numerical model of subduction, varying key parameters across a range of values deemed reasonable for the modern and Archean Earth. We show that under Archean conditions (i.e., with hotter mantle and thicker mafic crust), relamination of the mafic subducting crust is indeed viable under a wider range of subduction conditions. The observed dynamic behaviour can also account for other characteristics of Archean felsic terranes and related rocks, such as their episodic emplacement history , variable mantle signature (Smithies, 2000), and variable depth of primary melt formation, as well as the secular change of these properties throughout the Archean (Martin and Moyen, 2003).

Model physics
A 2-D Cartesian version of the finite element code Citcom is used (Moresi and Gurnis, 1996;Zhong et al., 2000) to numerically solve the governing equations that ensure the conservation of mass, momentum, thermal energy, and composition: where the primed variables represent dimensionless variables, non-dimensionalised as follows: and the thermal and compositional Rayleigh numbers, Ra T and Ra C , are defined as follows: All variables used in the above equations are defined in Table 1. We use the Boussinesq approximation and assume incompressibility, with no shear, adiabatic or radiogenic heating, which reduces the complexity of the system in order to allow for the clearer investigation of its primary features.
We use a purely viscous, non-Newtonian rheology for all materials. Brittle behaviour is not included as it is deemed not to play a significant role in our region of interest: namely the subducting crust and mantle wedge, below the over-riding plate. Instead, we impose a maximum viscosity of 10 24 Pa s, which limits the maximum stresses achievable within the rheologically stronger regions of the model. During preliminary work it was apparent that dislocation creep dominates throughout the model region. This is in line with the findings of Karato and Wu (1993) who put the dislocation-diffusion creep transition between 200 km and 300 km depth, and with the fact that dislocation creep is more likely to dominate in high-stress regimes such as around subduction zones. As such, the viscosity (η) is calculated as follows: All variables in the above equation are defined in Table 1. The effect of pressure on viscosity is assumed to be negligible as the model domain is shallow, and is partly offset by the omission of adiabatic heating, which has an opposite effect on the rheology. It is expected that the inclusion of pressure effects and adiabatic heating, would increase the viscosity at the bottom of the model area by no more than a factor of 3.2 (based on an activation volume of 14 × 10 −6 m 3 mol −1 for dislocation creep (Karato and Jung, 2003) and an adiabatic gradient of 7.8 • C/GPa). This is small compared to the viscosity uncertainties due to variations in flow law parameters arising from uncertainty in the subducting crustal composition as well as discrepancies between different rheological studies (see Fig. 2).
Material is defined to be either mantle or crust, both with different flow laws ( A, E and n) and densities. The density of crust relative to mantle material ( ρ) is assigned a fixed negative value above a certain "eclogitisation depth" and a fixed positive value below (see Table 1). This is included as a first order effort to simulate the effect of densifying phase changes (eclogitisation) expected to take place in the subducting crust (see Section 2.5 for details).

Model setup
The model domain is 450 × 150 km divided up into 512 × 256 finite elements. These elements are 600 × 600 m in a region of interest around the subducting slab and mantle wedge and 1200 × 600 m outside. The thermal and mechanical boundary conditions are illustrated in Fig. 1. The mechanical boundary conditions are as follows: V x = V z = 0 (no slip) on the surface boundary; V x = V trench (see below) and dV z /dz = 0 (stress free) on the bottom boundary at 150 km; dV x /dx = 0 (stress free) and V z = 0 (no slip) on the left hand side; and V x is linearly interpolated between V x = 0 at z = d plate and V x = V trench at z = h with V z = 0 (no slip) on the right hand side. The thermal boundary conditions are: T = 0 (0 • C) at the surface; T = 1 ( T or mantle temperature, see Table 1) at 150 km; the left hand side boundary is given the thermal profile of a half-space which has been conductively Rheological pre-factor for mantle dislocation creep Pa −n s −1 Acceleration due to gravity ms −2 9.8 -α Thermal expansivity K −1 3.5 × 10 −5 --

Fig. 1.
Schematic representation of the model setup with thermal and mechanical boundary conditions. Green area marks the position of crustal material and is given a thickness d crust , the location of the weak zone (see main text for further details) is depicted in orange, and the overriding plate is defined to extend down to d plate . The artificial slab pull force F is applied over an area within the slab mantle which is free to move laterally (indicated in blue). V trench is applied to the bottom boundary to mimic the effect trench migration has on the subduction geometry. * The outflow flow profile is V z = 0 above d plate and V x = V trench at the bottom boundary with V x linearly interpolated in between. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) cooled from above for a certain period of time: the "slab age" (see Table 1); and the right hand side boundary is given the linear thermal profile of an older plate where we assume that the base of the plate has been kept at a constant temperature by convective processes and the heat-flow through the plate is purely conductive and in steady-state.
In order to avoid modelling the full slab extending throughout the upper mantle, we approximate its presence by applying a slab pull force, F , achieved by applying a body force within the slab mantle at the bottom boundary. The bottom V x boundary condition, V trench , approximately mimics the effect of slab roll-back in the reference frame of the overriding plate, and allows for the sub- duction angle to be a controllable variable. Throughout the runtime of the model, the magnitudes of F and V trench are automatically adjusted such that a fixed, controllable, incoming plate velocity and subduction dip angle are achieved (see Table 1). By using this technique we neglect the effects of feedbacks between behaviour in our model and larger scale slab dynamics.
Randomly distributed tracers at an average density of ∼60 per element are advected with the flow to carry compositional information. Tracer velocities are conservatively interpolated between nodes to prevent gaps in tracer density forming in areas of high shear (Wang et al., 2015). A weak zone, 1.2 km wide, is imposed to decouple the two plates. Its position is defined using these tracers and it has a fixed low viscosity (see Table 1), which is logarithmically interpolated to the background viscosity over 10 finite elements (6 km).
During the model run time, the subducting crust is initially subjected to artificial stresses, resulting from the system finding its steady state geometry. We first run each model calculation with a uniformly strong (10 23 Pa s) crust, during an "initialisation stage" until it reaches a quasi-steady state temperature profile in order to prevent these stresses from deforming the crust. After this initialisation stage, all parts of the model are given their proper rheological properties, as described above.

Parametric study
Four default models are constructed for our parametric study. The first two "modern-day defaults" have reasonable values for a typical modern day subduction zone with convergence rates of 4.5 and 5 cm/yr, respectively. The other two "Archean defaults" are identical, but with an increased mantle temperature (+200 • C) and thicker subducting crust (15 km rather than 7 km). A parameter sensitivity study of these four default model setups is then undertaken by systematically isolating and varying key parameters, about the default model but within the parameter's bounds (see Table 1 for the values of these bounds). The varied parameters are incoming slab age, subduction convergence rate, subduction angle, mantle temperature, over-riding lithosphere thickness, crustal thickness, crustal buoyancy, "eclogitisation" depth (see Section 2.5), plate coupling (via viscosity of the imposed weak zone, see Section 2.2), mantle rheology, and crustal rheology. The default and bounds for each parameter are chosen to be representative for realistic subduction zones (and assuming that Archean crustal composition is similar to today's (see Section 4.4)) as well as to capture the full range of behaviours.

Rheological parameter values
It is recognised that the behaviour of the system is particularly sensitive to the choice of flow laws and as such the default and bounds for the mantle and crustal rheologies are informed by a thorough literature study. See Supplementary Material for a list of the experimentally derived flow laws for mantle and crustal dislocation creep used in this study. Values for the viscosity were calculated for a range of temperatures (300-1700 K) and strain rates (10 −18 -10 −10 s −1 ) for each flow law, and the logarithmic mean viscosity and standard deviation (SD) found. Flow laws were then fit to the mean and its ±1 SD using least squares regression.
The ±1 SD flow laws for mantle material are used as the bounds in this parametric study and the +1 SD as the default as this best represents dry peridotite. All mantle flow laws give modern-day upper mantle viscosities of order 10 20 Pa s (Fig. 2), in line with the conclusion of Hirth and Kohlstedt (2003) and estimates from isostatic rebound observations (Walcott, 1973). The experimentally derived crustal flow laws show a large variability with the collection of 15 studies showing a spread of viscosities over 4 orders of magnitude at temperature conditions of 800 K (Fig. 2). This is in part due to the inclusion of a variety of crustal materials, an approach deemed appropriate as the uppermost subducting crust is likely to consist of a mixture or "melange" of mafic crust and subducted sediments (Marschall and Schumacher, 2012) with the extent of mixing and the hydration state of these materials uncertain. All hydrated mafic flow laws lie near to or below the mean, with dehydrated mafic crust lying above and felsic material around the −1 SD flow law. As such, the −1 SD and mean crustal flow laws are used as upper and lower bounds for crustal material, respectively, and a flow law exactly halfway between used as the default. The default model therefore has a rheology that well represents the behaviour of hydrated mafic crust (see Supplementary Material for details). The upper bound of the crustal rheology is a realistic upper bound for hydrated mafic crust while the lower bound is realistic for a more felsic crust.
The effects of both dehydration (Hirth and Kohlstedt, 2003;Karato and Jung, 2003) (which can occur at highly variable depths Magni et al., 2014;Bouilhol et al., 2015) and partial melting (Hirth and Kohlstedt, 2003) on the rheology are uncertain. In addition to this, inclusion of dehydration involves modelling processes beyond the scope of this study. For these reasons, the effect of dehydration is not considered.

Crustal density and eclogitisation depth
Subducting mid-ocean ridge basalts densify via metamorphism, a process that we refer to here as "eclogitisation". The thermokinematic models of Duesterhoeft et al. (2014) demonstrate that this eclogitisation leads to gradual densification of the subducting mafic crust, which becomes negatively buoyant in the mantle below 100 km. However, there is seismic (Garth and Rietbrock, 2014), geodynamic (Gutscher et al., 1999;van Hunen et al., 2004) and field evidence (Austrheim, 1987;Hacker, 1996) showing that the mafic crust may remain metastable and therefore buoyant at greater depths, particularly if the crust is dry (Hacker, 1996). In addition, the effects of partial melts (not considered in this study) can reduce crustal densities considerably. The reader is referred to Section 4.4 for a more thorough discussion on the eclogitisation depth. As such, ρ crust = −300 kg/m 3 is used above an eclogitisation depth of 140 km (as default) and ρ crust = +200 kg/m 3 below. This eclogitisation depth and the crustal density above this depth are both varied as part of the parametric study. The effect of eclogitisation depth is investigated separately in more detail by running a further set of models that systematically explore the convergence rate -eclogitisation depth parameter space, the results of which are presented in Section 3.4.

Results
We performed four sets of parametric studies, in each of which, physical parameters were varied about a specific default model. Two default models apply to a modern day subduction system and the other two apply to an Archean setting, in which the mantle is made 200 K hotter and the crust is 15 km thick (Abbott et al., 1994) rather than 7 km thick. In the modern day default models, "normal" steady state subduction (i.e. no break-off or relamination) is observed. Normal subduction is also observed in the Archean default model with convergence of 5 cm/yr. However, upward removal of a portion of the mafic slab crust and subsequent underplating, henceforth referred to as relamination, occurs in the Archean default model with convergence of 4.5 cm/yr. These particular default models, with convergence rates of 4.5 cm/yr and 5 cm/yr, were chosen such that one Archean default model exhibits plume behaviour while the other does not, in order to allow for the investigation of the switching on/off behaviour of all parameters. Models that exhibit relamination appear to lie on a spectrum between a smaller scale, episodic behaviour that we refer to as "crustal plume" behaviour (Fig. 3B, described in section 3.1) Table 2 Summary of cases where relamination behaviour is observed, with the depth of detachment of material from the slab and the degree of relamination (a quantitative indication of the volume of crustal material expected to be removed from the slab. The method used to calculate the degree of relamination is described fully in the main text), and the peak vertical velocity of the crustal plume for each case. The depth of detachment is calculated as the shallowest depth of crustal material moving net upwards at the point of the formation of the first instability (or point "A" in Fig. 5). LB stands for "lower bound". Some models don't show relamination as they exhibit complete crustal decoupling instead (see Section 3.3 and Fig. 3C).  . 4. A graphical summary of the amount of relamination (see text for definition). Red symbols represent the Archean model runs and blue, modern day model runs. Darker, smaller symbols represent models where the default has a convergence rate of 4.5 cm/yr and lighter, larger symbols, 5 cm/yr. The circles represent Archean models which display relamination behaviour that is morphologically distinct from the normal "crustal plume" mode of relamination (see Sections 3.2 and 3.3). Crosses represent model runs which displayed slab break off instead of subduction (see Section 4.3). "LB" stands for lower bound, "def" for default and "UB" for upper bound. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) and larger scale, continuous behaviour that we term "whole crust" behaviour ( Fig. 3D, described in Section 3.2). Table 2 and Fig. 4 summarise the results of the two parametric studies, listing the conditions required to achieve crustal relamination. A key result to be noted from these results is that relamination does occur under a wider range of conditions in Archean models in comparison to modern day models.

Crustal plume behaviour
In models which display crustal plume behaviour, the top-most 2-3 km of crust becomes weak enough to form a Rayleigh-Taylor instability that subsequently feeds a plume consisting of crustal material which rises through the mantle wedge, driven by its compositional buoyancy (Fig. 5).
Convergence velocity strongly influences the development of these instabilities, with lower values promoting crustal plume development. When using the modern-day default with convergence 5 cm/yr, crustal plumes can be switched on by using the lower-bound crustal (i.e. felsic) rheology or decreasing the subduction convergence rate to ∼2 cm/yr. For the modern-day default with convergence 4.5 cm/yr, raising the mantle temperature or decreasing the crustal density also produces crustal plumes. In the Archean default with convergence 5 cm/yr, crustal plumes form under all the above conditions and also when the angle of subduction is significantly steeper than 45 • . The default Archean model with convergence 4.5 cm/yr is just slow enough to form crustal plumes without adjustment. But decreasing the eclogitisation depth to shallower than 140 km, increasing the slab age, decreasing mantle viscosity, decreasing the mantle temperature, increasing the over-riding plate thickness, decreasing the crustal density, or increasing the crustal strength all have the effect of switching off the crustal plumes (Fig. 4).  5 summarises the evolution of a model which shows typical crustal plume type behaviour. Relamination starts with the formation of an instability on the subduction interface (stage A), which corresponds to the start of an increase in the area of hot crustal material within the model, as this material is separated from the slab and heated in the mantle wedge. As the plume develops (stage B), crustal material is fluxed into the mantle wedge, which causes the mantle wedge to cool. At stage C, the mantle wedge has cooled by typically 100-300 • C and along with it, the surface of the newly incoming slab. At this point the surface of incoming crust is too cold and strong to relaminate, the source of the crustal plume is shut off and the crustal plume self-terminates. This is accompanied by a rapid removal of some of the colder relaminating material and surrounding cooled mantle. There is then a period where the ponded crustal material is gradually heated from below by the mantle wedge flow, with some of this material being caught in this flow and cycled around the mantle wedge. This causes some undulations in the thermal evolution of the crust (see curves between points C and D in Fig. 5). Eventually, the subduction interface becomes hot enough for another instability to form, which feeds a new crustal plume (stage D). Usually the second plume forms deeper in the model than the first as the mantle wedge is cooler overall, which is particularly apparent in the modern day models. The period of the episodicity of this behaviour is defined as the time between the formation of the two crustal plumes (points A and D) and is typically ∼10 Myrs although it ranges from 3 Myrs to 16 Myrs. In Archean conditions, crustal plumes rise fast through the mantle wedge, with velocity peaks of tens of cm/yr (40-80 cm/yr), typically taking only 100-200 kyrs to reach the base of the overriding plate.
All models that show crustal plume type behaviour have a similar evolution with points A, B, C and D easily identifiable. However, there is great variability in the depth at which these instabilities form and the amount of crustal material that is removed from the slab and eventually ponds below the over-riding plate. We define the depth of detachment as the shallowest depth of crustal material moving upwards at the point in time of the formation of the first instability (stage "A" in Fig. 5), and the amount of relamination as the area of crust that has passed the 900 • C isotherm, at the peak of the first crustal plume event. We use this measure as it offers a proxy for the amount of crust removed from the slab: the more crust that is removed, the more crust heats up to these high temperatures in the mantle wedge. In addition, it gives us a first order handle on the relative amounts of melt one could expect. These values are listed in Table 2, for all models that show relamination. Archean models lie on a spectrum between normal subduction with no relamination (Fig. 3A), through models which display smaller, more short-lived crustal plumes that form deeper in the model (up to 120 km, Fig. 3B), to the "whole crust" end member, described below, for which the instability forms shallower, ∼10-40 km below the over-riding plate, and crustal input into the mantle wedge is continuous (Fig. 3D). We observe a negative correlation between the depth of detachment and the amount of relamination in both the Archean and modern-day models, implying that crustal plumes that detach from the slab shallower are usually larger.

Whole crust relamination
It is possible for almost the entire crust to detach from the slab upon entering the mantle with only a small portion of material leaving through the bottom of the box (see Fig. 3D). This only occurs under Archean conditions and if the convergence is slow (∼2 cm/yr). In this model, the initial Rayleigh-Taylor instability forms at a depth between 75 and 90 km, or ∼10-40 km below the over-riding lithosphere. Unlike the crustal plume behaviour, this behaviour is not transient, with continued input of crustal material into the mantle wedge upon continued subduction and a new mantle wedge flow developing below the relaminating crust. As such, the definition of the "amount of relamination" (see Section 3.1) has to be altered for this model as there is no peak associated with a plume event. Instead, the amount of crustal material that has crossed the 900 • C isotherm after 10 Myrs from point A is used. All mantle material above the relaminating crust is cooled to <900 • C and therefore has lithospheric strength. The new mantle wedge below the crust behaves as it does under normal steady state subduction only with most crust flowing laterally above it, with only a small amount of crust, subducted with the slab, underlaying it.

Complete crustal decoupling
In Archean models where the lower crust is comparable in strength to the imposed inter-plate weak-zone prior to subduction, either because the crust is thick (>22 km) or rheologically weak, the lower crust assumes the role of the decoupling interface. The result of this is that the crust no longer subducts (Fig. 3C).

The effect of eclogitisation
A further set of model calculations was performed to determine the required eclogitisation depth for a given convergence rate (Fig. 6). It was found that the formation of the smaller, deeper plumes at moderate convergence rates (4 cm/yr), requires an eclogitisation depth of 140 km. However, for the larger, shallower relamination events which occur at slower convergence rates (2 cm/yr) an eclogitisation depth of 100 km is sufficient. The plumes that form at the faster convergence rates are smaller and form deeper in comparison to the crustal plumes that form at slower convergence rates. The result that relamination occurs at 4.5 cm/yr but not at 4 or 5 cm/yr when the eclogitisation depth is 140 km in Archean models, is due to the fact that this region is at the boundary between behavioural regimes. Any small changes under these conditions can result in markedly differing behaviours due to the fact that this system is very non-linear. This chaotic behaviour is also demonstrated in Fig. 4, panel 1.

Conditions for crustal relamination in modern subduction zones
Crustal plume behaviour can occur, under modern-day subduction zone conditions, when a weak rheology (e.g. for sediments or felsic crust; see Section 2.4) is used. Indeed, modern-day felsic and sedimentary material can relaminate from the subducting slab in the form of mantle wedge plumes (Gerya and Yuen, 2003;Behn et al., 2011;Marschall and Schumacher, 2012;Castro et al., 2013;Vogt et al., 2013), also known as "cold plumes" (Gerya and Yuen, 2003). The only condition under which material, which is as strong as hydrated mafic crust, can relaminate, is when the plate convergence rate is slow enough (<5 cm/yr or as low as ∼2 cm/yr when in a mantle at modern day temperatures and with a crust of relative density < +500 kg/m 3 ). This is primarily because, in slower subduction systems, the crust has more time to form Rayleigh-Taylor instabilities above the eclogitisation depth. In addition, while the topmost crust may heat up more slowly in slow subduction zones (Magni et al., 2014), the bulk of the crust heats up and therefore weakens more at shallower depth, which helps instabilities to form. This is particularly important for whole-crust relamination in Archean models (see Section 4.2).

Crustal relamination in the Archean
As in the modern-day cases, the parameter with the strongest control on relamination behaviour appears to be the subduction convergence rate. However, when subduction is slow (2 cm/yr) under Archean conditions, "whole-crust" relamination is observed (see Section 3.2). A much larger amount of crust is involved in the relamination process when compared directly to the modern-day equivalent (amount of relamination ∼4 times higher). In addition to this, crustal plume relamination now occurs under a far wider range of conditions. These results give strength to our original hypothesis that the relamination of the subducting crust was more viable in the early Earth.
The fact that one Archean default model displays relamination and the other does not allows the effects of all parameters to be gauged. The result that increased crustal buoyancy promotes relamination is perhaps unsurprising as with higher buoyancy forces, instabilities form more readily. Subduction with a hotter mantle means the crust heats faster, which leads to faster weakening and promotes the formation of instabilities. Steepening subduction means that slabs are more aligned with the direction of gravity which leads to increased stress weakening.
Whether subduction operated more slowly in the Archean is unclear. More dehydration stiffening near Archaean MORs may have led to consistently slower plate tectonics in the early Earth (Korenaga, 2006) (though this has been argued against Davies, 2009) or perhaps subduction was episodic, interrupted by intermittent slab break-off events (van Hunen and van den Berg, 2008;Moyen and Van Hunen, 2012). Episodic subduction may have had a particularly dramatic effect on relamination as during periods of slab stagnation the convergence rate is effectively zero. Stagnating mafic slab crust would heat up in such a situation, ideal conditions for the formation of crustal plumes. Due to the way we impose an artificial slab pull force, our models were not designed to investigate slab breakoff accurately beyond whether it occurs or not (see Section 4.3). As such, an investigation of this behaviour is outside the remit of this study but could be a subject for further study.
The complete crustal decoupling observed when the crust is made weak (felsic-like in strength) or thick (22 km) (see Section 3.3 and Fig. 3C) has previously been suggested as a possible mode of plate tectonics in the Archean ("flake tectonics") (Hoffman and Ranalli, 1988). Based on proposed Archean conditions, it was argued that the thicker crust, which is hot and weak at its base, may have been more decoupled from the underlying mantle lithosphere, allowing it to fully separate from the subducting slab upon subduction leaving only the mantle lithosphere to sink deeper into the mantle. This mechanism has been shown to be dynamically viable in other modelling studies, provided that a weak enough layer underlays the crust (Davies, 1992;Vogt and Gerya, 2014).

Slab break-off
Previous studies have suggested that, with a thicker crust, subduction may not have been viable in the early Earth (Davies, 1992;Sizova et al., 2010). This is because, with a thicker crust, buoyancy forces resisting subduction are higher and the overall lithospheric strength weaker. Indeed, in Archean models we find that slab break-off is more likely to occur when the crust is thickened (to ∼29 km in our models) or the slab made too young (∼20 Myrs), shutting down subduction. However, subduction with episodic, spontaneous break-off has been shown to be a viable mode of subduction in the Archean and may also explain the absence of ultra-high-pressure metamorphic rock from this time (van Hunen and van den Berg, 2008;Sizova et al., 2014).

The composition and density of the subducting crust
We have undertaken an extensive parametric study to constrain the conditions required for crustal relamination and investigate the effect of various subduction parameters on this behaviour. In our default models, the depth to which the crust remains buoyant (i.e. the eclogitisation depth) is 140 km. Additional simula-tions demonstrated that, for larger crustal plumes to develop at slower velocities, an eclogitisation depth of 100 km is sufficient. However an eclogitisation depth of 140 km is required for the development of smaller, deeper crustal plumes (see Section 3.4). The reason for choosing 140 km as a default eclogitisation depth is motivated by the evidence that the subducting crust remains metastable during its prograde path. Indeed, the geological record of exhumed ultra-high-pressure terranes (e.g. Western Norwegian Caledonides, Sulu belt, Eastern China) show that mafic material has resisted the transformation to eclogite well within eclogite-facies conditions (some as deep as >90 km), save for along hydrated veins (Austrheim, 1987;Zhang and Liou, 1997). This petrological evidence for a delayed eclogitisation, has been confirmed by the geophysical observation that slab crust remains seismically slow compared to mantle to great depths (>150 km) (Abers, 2000;Garth and Rietbrock, 2014). Such significant delays in eclogitisation have also been invoked to explain flat slab subduction of oceanic buoyant plateaus, arguing that the bulk of this plateau crust must remain buoyant at the depth at which the slab flattens (typically ∼150 km) for timescales of several Myrs (Gutscher et al., 1999;van Hunen et al., 2004). We have therefore used these observations to postulate that if the subducting crust is subject to eclogitisation delays, then the bulk crust may remain metastably buoyant.
As highlighted by the results in Section 3.4, the deepest, smallest scale crustal plumes require a deep eclogitisation depth and are therefore most likely to be prevented by densification reactions. However, these small plumes only form from the topmost 1-2 km of subducting crust, which is likely to consist of a melange of mafic subducting crust and sediments (Gerya and Yuen, 2003;Marschall and Schumacher, 2012). In this case, the physical addition of sediments to the oceanic crust would significantly affect the bulk density of this heterogeneous layer, leading to an even more buoyant slab top. In addition, along the warm subduction gradients of an "Archean subduction zone" the oceanic slab crust may undergo partial melting (Bouilhol et al., 2015) leading to significantly lower density (Gerya and Yuen, 2003).
In this study, our modelled Archean crust is assumed to be thicker than its modern equivalent but similar in composition. The nature of the Archean oceanic lithosphere is much debated. At higher mantle potential temperatures, more mantle material is processed at the ridge, leading to a thicker, MgO-rich oceanic crust. This thicker crust might have been stratified, in a similar way to modern oceanic plateaus, being more basaltic (MORB-like) at the top and more ultramafic/komatiitic at depth (Kerr et al., 1997). As such, MORB properties might be appropriate for the topmost Archean crust but the lower crust may have been denser and stronger. This would affect neither the eclogitisation depth nor the crustal plume behaviour, but could hinder the removal of the deeper parts of the crust from the slab during whole-crust relamination.

Magmatic output and Archean continental crust formation
Although studying the magmatic output of this process is beyond the scope of this work and we do not compute it in our models, it is worth noting that during normal subduction, under Archean conditions, very little crust heats up to above 800 • C and no crustal material heats up to above 900 • C within our model region. However, as Fig. 4 illustrates, the area of mafic crustal material that heats up to above 900 • C when relamination occurs (the "amount of relamination") can be substantial (typically a few 100 km 2 , but ∼1000 km 2 after 10 Myrs in the whole-crust relamination case). At this temperature, partial melting of mafic material is likely to occur, particularly if the system is being percolated by metamorphic fluids (Magni et al., 2014;Bouilhol et al., 2015). In all cases crustal plume material reaches these melting temperatures within the mantle wedge (at >50 km depth in our models) and, as such, has the potential to explain the high pressure mafic signature, pervasive in Archean terranes (Martin, 1999;Foley et al., 2002).
Episodic crustal plumes may result in episodic magmatism (Marschall and Schumacher, 2012) and previous studies on Archean terranes have suggested episodicity in their formation with a typical timescale of a few Myrs , the same timescale under which crustal plumes operate in our models. The variability of the depth of formation of the crustal instability, the variability in the size of the subsequent relamination event and the fact that shallow instability formation correlates with more extensive relamination imply that the degree of crust-mantle interaction may also vary significantly. When crustal plumes form deep and are small, the degree of mantle interaction is likely to be high whereas in the whole-crust relamination case, the crustal relaminant does not have to pass through the hot mantle wedge at all and potential melts from the ponded crustal material will only pass through mantle lithosphere which is <900 • C.
Our results demonstrate that more extensive and shallow relamination is more likely under the hotter mantle conditions of the early Earth. As the mantle probably cooled through the Archean (Herzberg et al., 2010), crustal relamination may have gradually transitioned from shallow and extensive to small-scale and deep. The result of this may have been a gradual transition from crust derived melting at relatively shallow depths with minimal mantle interaction, to crust derived melting at greater depths with greater mantle interaction. Interestingly, this trend in depth and degree of mantle interaction is observed in the composition of Archean felsic terranes of different ages (Martin and Moyen, 2003). This has previously been cited as evidence for, on average, shallower subduction in the Archean which gradually steepens through Earth's history (Abbott et al., 1994), or a secular change in the depth of slab melting through time (Martin and Moyen, 2003). Here we offer secular changes in subducting crust dynamics as an alternative explanation.

Conclusions
We performed two sets of parametric studies to investigate how different key parameters can affect subduction dynamics and, in particular, the possibility for the detachment of the subducting mafic crust from the slab and subsequent upwelling and underplating (or "relamination"), in a modern day setting and early Earth setting (with a 200 K hotter mantle and 15 km oceanic crust). Under modern day conditions, crustal relamination is only observed when the crust has a "felsic" flow law (in agreement with previous studies) or when subduction is slow (<5 cm/yr or ∼2 cm/yr for modern mantle temperatures and a realistic crustal density). Slowing down subduction, raising mantle temperature and decreasing crustal density or strength also enhances relamination in the early Earth models. However, in these Archean models, so does making the over-riding plate thinner, making the subduction angle significantly steeper or shallower than 45 • , making the slab younger, or making the mantle stronger. Relamination of the mafic crust requires that the crust remains buoyant in the mantle to depths greater than the equilibrium eclogitisation depth in all but the most extreme (slow convergence) cases. This may be achieved by mixing with more buoyant sediment, partial melting or the slow kinetics of eclogitisation in the absence of water.
Crustal strength and subduction convergence rate have the strongest control on behaviour. Indeed, when subduction is slow enough (∼2 cm/yr), the entire mafic crust is able to detach from the slab and relaminate. We demonstrate that crustal relamination behaviour exhibits a spectrum between this whole-crust end member (which occurs at relatively shallow depths) through an episodic "crustal plume" mode (which occurs deeper and is typically episodic on a ∼10 Myr timescale), to no relamination. We suggest that this trend may explain secular changes in the degree of mantle interaction and depth of primary magma production of Archean felsic terranes over time. It is also observed that, under early Earth conditions, continuous subduction is less viable if the oceanic crust is too thick (> ∼20 km), the slab is too young (∼20 Myrs), or the inter-plate coupling high, due to spontaneous slab-breakoff. Finally, when drawing a direct comparison between the modern day and early Earth parameter studies, we find that relamination of subducting mafic crust is a viable mechanism under a wider range of subduction conditions in an early Earth subduction zone and has the potential to be more substantial. This has implications for the generation of magmatism in Archean subduction zones, Archean subduction dynamics, and the long term recycling of mafic crust.