Can metasomatic weakening result in the rifting of cratons?

Accepted Manuscript Can metasomatic weakening result in the rifting of cratons?

Stefanie Wenker, Christopher Beaumont PII: DOI: Reference:

S0040-1951(17)30259-7 doi: 10.1016/j.tecto.2017.06.013 TECTO 127525

To appear in:

Tectonophysics

Received date: Revised date: Accepted date:

30 January 2017 12 May 2017 8 June 2017

Please cite this article as: Stefanie Wenker, Christopher Beaumont , Can metasomatic weakening result in the rifting of cratons?, Tectonophysics (2017), doi: 10.1016/ j.tecto.2017.06.013

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT Can metasomatic weakening result in the rifting of cratons? Stefanie Wenker1 and Christopher Beaumont2 1. Earth Sciences Department, Dalhousie University, Halifax, NS, B3H 4R2, Canada

CR

IP

T

2. Oceanography Department, Dalhousie University, Halifax, NS, B3H 4R2, Canada

US

Highlights

Cratons are cold and strong, yet they occasionally rift, e.g. North China craton



Can mantle metasomatism weaken cratons sufficiently that they will rift?



Numerical models with progressive craton melt metasomatism are rift tested



Combined metasomatic heating, refertilization and densification can lead to rifting



Once sufficiently weakened cratons are not protected by surrounding standard lithosphere

PT

ED

M

AN



Abstract

CE

Keywords craton weakening, metasomatism, rifting, numerical models

AC

Cratons are strong and their preservation demonstrates that they resist deformation and fragmentation. Yet several cratons are rifting now, or have rifted in the past. We suggest that cratons need to be weakened before they can rift. Specifically, metasomatism of the depleted dehydrated craton mantle lithosphere is a potential weakening mechanism. We use 2D numerical models to test the efficiency of simulated melt metasomatism and coeval rehydration to weaken craton mantle lithosphere roots. These processes effectively increase root density through a parameterized meltperidotite reaction, and reduce root viscosity by increasing the temperature and rehydrating the 1

ACCEPTED MANUSCRIPT cratonic mantle lithosphere. The models are designed to investigate when a craton is sufficiently weakened to undergo rifting and is no longer protected by adjacent standard Phanerozoic lithosphere. We find that cratons only become vulnerable to rifting following large-volume melt metasomatism (~30% by volume) and thinning of the gravitationally unstable cratonic lithosphere from >250 km to

IP

T

~100 km; at which point its residual crustal strength is important. Furthermore, our results indicate

CR

that rifting of cratons depends on the timing of extension with respect to metasomatism. An important effect in the large-volume melt models is the melt-induced increase in temperature which must have

US

time to reach peak values in the uppermost mantle lithosphere before rifting. Release of water stored in the transition zone at the base of a big mantle wedge may provide a suitable natural setting for both

AN

rehydration and refertilization of an overlying craton and is consistent with evidence from the eastern

M

North China Craton. An additional effect is that cratons subside isostatically to balance the increasing

ED

density of craton mantle lithosphere where it is moderately metasomatized. We suggest that this forms intracratonic basins and that their subsidence and subsequent uplift, and cratonic rifting

AC

CE

PT

constitute evidence of progressive metasomatism of cratonic mantle lithosphere.

2

ACCEPTED MANUSCRIPT 1. Introduction Cratons are distinctly different from Phanerozoic lithosphere (e.g. Griffin et al., 2009), in that the latter actively participates in tectonics and is recycled continuously, whereas cratons resist deformation. The high yield strength of cold cratonic mantle lithosphere (CML), combined with the

IP

T

secondary effects of high viscosity and neutral buoyancy (Jordan, 1975), has been used to explain the

CR

longevity of cratons; supported by studies of CML rheology (Griffin et al., 2009; Hirth and Kohlstedt, 1996; Karato, 2010) and geodynamics (Beuchert and Podladchikov, 2010; Doin et al., 1997; François et

US

al., 2012; Lenardic et al., 2003; Lenardic and Moresi, 1999; O’Neill et al., 2008; Sleep, 2003; Wang et al.,

AN

2014; Yoshida, 2012). Additionally, it has been shown that strong lithosphere is protected against deformation by surrounding weaker lithosphere which rifts more readily (Yoshida, 2010, 2012;

M

Lenardic et al., 2000, 2003; Gorczyk et al., 2013; Corti et al., 2007; Wenker and Beaumont, 2016.).

ED

Since cratons are strong and protected by surrounding weaker lithosphere, why do some thin

PT

and eventually rift? This certainly happened to the North China Craton (Menzies et al., 2007; Tang et al., 2013b), the North Atlantic Craton (Sand et al., 2009; St-Onge et al., 2009; Tappe et al., 2011, 2007)

CE

and is starting in the Tanzanian Craton (Aulbach et al., 2008).

AC

Geochemical data from xenoliths and xenocrysts have provided evidence that even stable cratonic roots may be altered by fluid and/or melt metasomatism (Carlson et al., 2005; Griffin et al., 2009; Lee et al., 2011; Rudnick et al., 1998). This widespread process can cause refertilization, rehydration, and increase the temperature, thereby counteracting the characteristics that promote craton survival (Lee, 2003; O’Reilly et al., 2001; Schutt and Lesher, 2010). Refertilization increases the density of originally depleted lithospheric mantle. Rehydration decreases its bulk viscosity by a factor of up to 140 at constant stress (Hirth and Kohlstedt, 1996), but only by a factor of 10 at constant strain 3

ACCEPTED MANUSCRIPT rate (e.g. Wang et al., 2014). Increasing temperatures decrease the temperature-dependent viscosity, and may also induce fractional melting, further decreasing the viscosity (Hirth and Kohlstedt, 1995). Overall, metasomatism alters the mantle lithosphere and may destabilize the cratonic root (Schutt and Lesher, 2010). We propose that this type of weakening is required in order to rift an originally strong

IP

T

craton.

CR

We use 2D-numerical models to investigate the rifting of cratonic lithosphere following varying degrees of simulated metasomatism. The models comprise welded adjacent regions of ‘standard’ and

US

‘cratonic’ lithosphere, which is depleted and dehydrated. Both regions contain equal small inherited

AN

weaknesses representing vestiges of previous tectonic processes, namely; weak sutures, shear zones, preferred foliations (Manatschal et al., 2015; Tommasi and Vauchez, 2001) and may serve to localize

M

deformation. The models are used to investigate where extension localizes in the model and under

ED

what circumstances rifting shifts from the standard lithosphere to the craton as the CML is

PT

metasomatically weakened. Specifically, we consider the effects of basaltic-melt metasomatism and coeval rehydration the CML. We parametrically approximate the melt-peridotite reaction, thereby

CE

transforming hartzburgite to lherzolite, which increases the density. This contrasts with the effect of

AC

crystallization of the basaltic melt alone. We also include rehydration as part of the metasomatism. Because basaltic- melt metasomatism alone will not rehydrate the CML we assume that hydrous fluids from the same or a different source rehydrate the CML and reduce its viscosity up to a factor of 5. In order to simplify the calculations and the number of variables the parametric rehydration occurs on the same timescale as the refertilization in the models, that is, they are coeval. Although this is an approximation it allows us to investigate the combined effects of both processes. As we show in the results, large volume melt metasomatism of approximately 30% by volume is required to weaken 4

ACCEPTED MANUSCRIPT cratonic lithosphere sufficiently to rift. We therefore initially focus on large volume metasomatism models, to determine the maximum effect, and later analyze lesser amounts and the individual contributions of the refertilization, heating and reduced viscosity. Metasomatism starts under neutral boundary conditions and may cause convective removal of the lower CML. The models are then rift

IP

T

tested near the end of, or after, metasomatism.

CR

In what we term M-type models melt is added uniformly to the whole CML over 3, 10, 30 and 50 Ma, respectively, with lesser amounts close to the edge of the craton. These models demonstrate

US

convective removal of the lower CML and the effects of heat diffusion as a function of the duration of

AN

metasomatism. We next employ more geologically realistic progressive upward migrating metasomatism (Le Roux et al., 2009, 2007; Lenoir et al., 2001; Soustelle et al., 2009) using a reactive

M

front model, termed RF. We test a range of models to see whether rifting localizes in the cratonic or

ED

standard lithosphere. We then assess the results in terms of published modeling and observations and

PT

offer a mechanism whereby the necessary weakening could be achieved naturally. Finally, we note the possible link between refertilization of CML and the formation of intracratonic basins.

CE

2. What makes cratons vulnerable to rifting?

AC

Many modeling studies that investigated the long term stability of cratons focused on thermal erosion by convective mantle entrainment (Doin et al., 1997; Shapiro et al., 1999; Sleep, 2003; Wang et al., 2014) and delamination (Gorczyk et al., 2013, 2012). This process has been suggested for the Colorado Plateau (Bashir et al., 2011; Levander et al., 2011), North China Craton (Gao et al. 2009), and the North Atlantic Craton (Sand et al., 2009; Tappe et al., 2011, 2007). Here we take the view, supported by studies of craton longevity noted above, that these processes may thin the CML, but not weaken it sufficiently to cause rifting when surrounded by weaker lithosphere. Instead, we consider 5

ACCEPTED MANUSCRIPT the effects of melt and hydrous fluid interaction with depleted CML peridotite, processes that change the internal characteristics that make cratons strong. Widespread refertilization of CML is supported by major and trace element data, and Sr, Nd, and Os isotopic compositions from xenoliths (e.g. Griffin et al., 2009; O’Reilly and Griffin, 2013; Tang et

IP

T

al., 2013). Mantle metasomatism is the interaction of fluids (by definition) and melts (by extension)

CR

with the CML (Foley, 2008; O’Reilly and Griffin, 2013; Zheng et al., 2015). As noted by O’Reilly and Griffin (2013), “The lithospheric mantle is a palimpsest recording the multiple fluid events that have

US

affected each domain since it formed. These events, involving different fluids and compositions, have

AN

repeatedly overprinted variably depleted original mantle wall-rocks. This produces a complex, essentially ubiquitously metasomatised lithospheric mantle, heterogeneous on scales of microns to

M

terranes and perhaps leaving little or no “primary” mantle wall-rock”. We note, however, that in some

ED

cases the observed metasomatism may have developed after the cratonic lithosphere had already

PT

been thinned, so caution should be used in applying all of the evidence to thick cratons. Here we take a simplified approach and model the effects of refertilization that results from the

CE

reaction between infiltrating basaltic melt, which converts harzburgites to lherzolites and increases the

AC

density by increasing the FeO concentration (Schutt and Lesher, 2006). Metasomatism by hydrous fluids can also rehydrate the CML thereby decreasing its viscosity, making it easier to deform, and may produce minor melts (Foley, 2008; Lee et al., 2009; Zheng et al., 2015). Additionally, delivery of melts to the CML increases the temperature (O’Reilly and Griffin, 2013), which, owing to thermal expansion, partly offsets the metasomatic increase in the density, but also weakens the CML by decreasing the viscosity. We also include these effects in the modeling.

6

ACCEPTED MANUSCRIPT 3. Summary of modeling methods Calculations (details Appendix A) are made with the arbitrary Eulerian-Lagrangian finiteelement thermo-mechanical software SOPALE-nested (Beaumont et al., 2009; Butler et al., 2014; Fullsack, 1995). Deformation is described by frictional-plastic and thermally activated power-law

IP

T

viscous rheologies based on laboratory experiments. The latter are scaled by a factor f to change the

CR

relative strength among models and as the models evolve as explained below. The upper-mantle-scale 2D models (Figure 1a) comprise adjacent welded regions of thermally

US

stable standard (left) and cratonic lithosphere (right) containing small weak regions that may serve to

AN

initiate necking instabilities when the model extends under velocity boundary conditions. We include the weak regions on the basis that continental lithosphere is extensively reworked and not pristine,

M

and will therefore contain one or more weak regions. Extensional boundary conditions (‘rifting tests’)

ED

are applied in the late- to post-metasomatism phase.

PT

The initial thickness, thermal regime, and material properties change across the lithospheric boundary at the center of the model (Figure 1). Crust is uniformly 35 km thick. Standard and cratonic

CE

lithospheres are 140 and 250-km thick, respectively. Part of the base of the CML is inclined to reflect

AC

prior flow and thinning of the CML and to represent an initially stable system. The thermal parameters within each lithosphere differ (Table 1) resulting in Moho temperatures of 600°C (standard lithosphere) and 450°C (craton). The base of both lithospheres is initially at 1350°C. The CML is initially depleted (low density) and dehydrated (strong). We investigate what degree of weakening makes cratons vulnerable to rifting and we parameterize the weakening of the CML in terms of refertilization by basaltic melt metasomatism and rehydration. There are ten 20-km thick layers in the CML and one 15-km thick layer just below the Moho (Figure 1 blue/ pink layers). This 7

ACCEPTED MANUSCRIPT layering allows us to metasomatise the whole lithosphere concurrently (M-type models), or a layer at a time (RF-type models). The amount of metasomatism decreases from the central craton to its boundary over a distance of 200 km, with no metasomatism in the 50 km adjacent to the standard lithosphere, Figure 1. This design represents metasomatism of craton interiors, thereby avoiding

IP

T

complexities caused by metasomatism spanning the boundary with the standard lithosphere. An alternative would be to metasomatise the mantle lithosphere in columns. However, when

CR

viewed in 3D the metasomatised column would need to extend across most of the craton in the 3rd

US

dimension to achieve sufficient bulk weakening. This requires metasomatism that is effectively craton wide in at least one horizontal dimension. Our models are 2D so we only have the choice of the in-

AN

plane dimension of the model. Therefore we choose the layer approach to metasomatism. We believe

M

that equivalent necking of the lithosphere would result if the scale of metasomatism were restricted to

ED

the horizontal length that gives the same necking style (approximately 2-3 times the thickness of the lithosphere). However, both scaling arguments require testing, particularly with fully three dimensional

PT

models.

CE

Only two laboratory-determined flow laws are used in these models, ‘Wet’ Quartzite (WQz) (Gleason and Tullis, 1995) for the crust and ‘Wet’ Olivine (WOl) (Karato and Wu, 1993) for the mantle

AC

(Table 1). These flow laws are scaled by scaling factor f to increase/decrease the relative strength (Butler et al., 2014, Appendix). This simplification allows for a relative comparison between materials of different strengths without having to incorporate other flow laws (Beaumont et al., 2006; Jamieson et al., 2007; Karato, 2010), and allows the juxtaposition of the standard and craton lithospheres with different properties in the rifting tests. The scaling approach is also valid for a considerable range of f values, even for the reference flow laws, owing to the significant experimental uncertainties associated with the laboratory-determined flow-law parameters (Beaumont et al., 2006; Butler et al., 2014; 8

ACCEPTED MANUSCRIPT Jamieson et al., 2007; Karato, 2010; see Butler et al., 2014 for a more complete explanation and figure showing scaled effective viscosities versus temperature for a range of f values ). The models presented have a constant scaling factor of WQz x 3 for the standard crust (SL1) and WQz x 4 (SL2), and WQz x 5 for the cratonic crust in C1 and WQz x 4 in C2 (Table 2). Scaling factors

IP

T

of f = 3-5 are an appropriate approximation for a mixed quartz-feldspar controlled crust, or ‘dry’ quartz

CR

with increasing feldspar content. We scale the flow laws for the mantle lithospheres to represent dehydration of peridotite controlled by the olivine flow law (WOl). In SL1 the standard lithosphere has

US

WOl x 3, and in SL2 it is WOl x 4, both accounting for moderate dehydration expected in Proterozoic or

AN

Phanerozoic mantle lithosphere (Griffin et al., 2009). The CML starts with an effective viscosity of WOl x 5 corresponding to low water content (<0.001 wt% H2O), and close to the ‘dry’ olivine flow law

M

(Karato, 2010, Butler et al. 2014, Appendix). As the CML is metasomatized the scaling factor decreases

ED

to 1 (as explained below) so that it reaches WOl x 1 by the end of simulated metasomatism. This

flow law for wet olivine.

PT

reproduces a decrease in the viscosity expected with rehydration of the CML such that WOl x 1 is the

CE

4. Simulated craton metasomatism in the numerical models

AC

4.1 Refertilization and rehydration Parameterized melt metasomatism, in the form of uniformly distributed basaltic melt at 1300°C, is added from the 1350 °C sublithospheric mantle at steady fractional volumetric rate, RBM, to selected layers in the CML for a given time interval, tBM . We assume that melt infiltrates rapidly into the layer that is being metasomatized, with velocities of 20-50 m a-1 ( O’Reilly and Griffin, 2013; Bodinier et al., 1990; Le Roux et al., 2007) or faster by diking (Johnson and Jin, 2009), thereby making it nearinstantaneous for our model time scales (McKenzie, 1989; Ranalli et al., 2007; Vernières et al., 1997) 9

ACCEPTED MANUSCRIPT and without significant loss of heat during delivery. We do not specify where the melt is produced, which could be from a mantle plume (e.g. Sobolev et al., 2011), melting in the mantle wedge above a subducting slab located beyond our right boundary (e.g. Zhu et al., 2011), or most likely from water release and melting in the mantle transition zone (Karato, 2011; Liu et al., 2016). We discuss these

IP

T

possible sources in Section 8.4. In most models the melt volume added is 30% and it is added at a

CR

uniform rate determined by the duration of metasomatism. We choose 30% to approximate melt addition required to invert melt extraction required to form CML from Archean primitive mantle

US

(Bernstein et al., 2007; Griffin et al., 2003; Lee, 2003). For each time step the heat and mass balance of melt added to each layer is treated in a self-consistent manner and we use a simple simulated reaction

AN

model between the depleted peridotite and the basaltic melt in which both the heat and mass of the

M

metasomatized layer are conserved.

ED

As pieces of detached CML sink, the material properties are changed to those of the

PT

sublithospheric mantle (density 3378 kg m-3, WOl x 1) after they reach a depth of 500 km. These pieces of CML can only reach this depth after delamination and this is normally close to the end of

CE

metasomatism. The 500 km depth was chosen to minimize interference with current delamination by

AC

previously delaminated CML. 4.2 Conservation of heat

Basaltic melt is injected during each time step at a specified temperature, TBM, and is normally at higher temperature than the surrounding lithospheric layer and therefore heats this layer. We take into account cooling of the basaltic melt to the local ambient temperature of the lithosphere, which may include the release of latent heat when melt solidifies. For each Eulerian element in the finite

10

ACCEPTED MANUSCRIPT element model the heat released by a fractional volume

of the melt, referred to as material 2, is

equal to:

is the specific heat capacity of the basaltic melt,

the difference in temperature between

T

where

IP

the melt and the ambient temperature, ratio between the actual temperature change and the total the melt fraction and

the latent heat of

CR

temperature change of the region (see below),

US

crystallization. The other variables are the same as described above.

There are four specified cooling intervals based on Annen et al. (2006) (Figure A.1): 1) specific

AN

heat released during cooling to the liquidus temperature of material 2; 2) specific and latent heat

M

released between the liquidus and the intermediate temperature (Annen et al., 2006); 3) specific and

ED

latent heat released between the intermediate temperature and the solidus temperature of material 2, and; 4) specific heat released below the solidus temperature. The intermediate temperature reflects

PT

changes in melt fraction with temperature that has two linear segments (Figure A1).

CE

The total heat released by cooling of material 2 is absorbed locally by material 1, the peridotite, and is

AC

described by a similar equation:

is the heat absorbed by material 1, equation,

is the change in volume calculated from the mass balance

the initial volume of the element at the beginning of the timestep,

capacity of material 1,

the density of material 1,

the specific heat

the melt fraction for either a low or high

temperature region between the solidus and liquidus, and

the latent heat of fusion for material 1.

The temperature regions are defined in the same manner as for the heat released, but with 11

ACCEPTED MANUSCRIPT parameters specific to material 1. However, in these models the temperature never goes above the solidus temperature of peridotite and so the temperature of material 1 only changes because specific heat is absorbed (Hirschmann, 2000). The thermal parameters of material 1 are constant, although heat is technically absorbed by the evolving peridotite. This simplification is justified because the

IP

T

density of the peridotite increases by a very small amount.

CR

Strictly speaking, this calculation should be iterated because the final temperature at the end of the time step is higher than the ambient temperature assumed in the cooling of the melt. However,

US

the amount of melt added each time step is sufficiently small that this temperature difference is

AN

negligible. 4.3 Conservation of mass and rehydration

M

The original material is the depleted CML, a peridotite with WOl x 5 and a density of 3335 kg m-3

ED

(corresponding to Mg# 92). The basaltic melt is injected at a constant rate calculated by the total

PT

volume change required divided by the length of time of metasomatism. Simulated metasomatism occurs over specified time intervals and will, when complete, increase the initial volume by a specified

CE

amount, approximately 30% in most cases for the models used here.

AC

Although the temperature at which material 2 is injected dictates all of it to be melt, it is injected with the density of the solid. We justify this by assuming that the melt crystallizes in less time than a model time step and therefore has a solid density by the end of the timestep. The basaltic melt (material 2) is assumed to react with the peridotite (material 1) to make it less depleted. The effect of this reaction is simulated by progressively increasing the peridotite density from 3335 kg m-3 (Mg# 92) to a specified ‘fertile’ maximum of 3378 kg m-3 (Mg# 89), based on Lee (2003). By the time 30% by volume of melt has been added to the CML its material density is equal to 12

ACCEPTED MANUSCRIPT that of the sublithospheric mantle (Bernstein et al., 2007). This range of Mg# may be a conservative choice given Mg#’s as low as 86 are observed (Zheng et al., 2015). However, such low values observed in xenoliths may be local and not ubiquitous. In addition, the important effects are the change in density, not the absolute values, and the difference between peridotite density and that of the

IP

T

underlying sublithospheric mantle (3378 kg m-3).

CR

The volume change of material 1, required by the mass balance during the simulated reaction

is the constant volume of material 2 added per time step,

AN

Where

US

is described by the following equation:

M

the volume of material 1 at the beginning of the timestep,

the initial density and final

the fractional change in volume of material 1. As

ED

density for material 1 during this time step, and

and,

the density of material 2,

is explained in the Appendix the finite element calculation is incompressible (Equations A.1-A.3) but of each model layer is updated at the end of each time step by

PT

the fractional change in volume

CE

geometrically dilating the layer .

We assume that rehydration of the CML also occurs during metasomatism by the addition of

AC

hydrous fluids. These fluids may have a different source from the melt. We do not model this process per se, but instead linearly decrease the scaling factor f from WOl x 5 to WOl x 1 to simulate rehydration over the range of metasomatism. We assume the rehydration is coeval with the refertilization. The CML becomes more deformable as the viscosity of the mantle lithosphere decreases. In nature rehydration would have the additional effect of lowering the solidus temperature. This is not taken into account in the models presented because the decrease in melting temperature of

13

ACCEPTED MANUSCRIPT the peridotite is not sufficient to produce melt in the CML (Hirschmann, 2000). There are no model volume or mass changes related to hydration. In summary, metasomatism has four effects in the models. 1) As hot melt cools all heat released, including latent heat, is absorbed by the CML, increasing its temperature and decreasing its

IP

T

viscosity. 2) The melt reacts with the peridotite to make it less depleted thereby increasing its density.,

CR

3) The CML evolves from a dry to a hydrated rheology (Appendix A) by linearly changing the scaling factor in the viscous flow law from WOl x 5 to WOl x 1 as metasomatism proceeds assuming hydration

US

parallels the 30% melt addition. 4) The volume of the CML increases by the volume change of the

AN

peridotite after reaction with the melt. This volume increase is balanced by isostatic subsidence of the CML which pumps sublithospheric mantle through the side boundaries of the model to equilibrate the

M

model basal pressure with an external reference pressure. Thus the volume of the sublithospheric

PT

5. Rifting cratons in the models

ED

mantle, the melt source, decreases by the melt volume and source and sink volumes are balanced.

CE

5.1 Stability of the initial model configuration All models were computed for an initial 1 Ma phase to allow for isostatic adjustment of both

AC

lithospheric regions prior to metasomatism. This produces a stable system with weak edge-driven convection (Figure 1b), which when computed for 100 Ma remains close to stable, showing only minimal buoyant lower lithosphere counterflow (<100 km) at the vertical lithospheric boundary and minor smoothing of the inclined base of the CML (Figure 1c). We conclude that the starting model is stable, unlike initially unstable states in some models criticized by Sleep (2007, 2003), and that any additional destabilization of the CML is the result of the metasomatism.

14

ACCEPTED MANUSCRIPT 5.2 Rifting tests The models are computed through the metasomatism phase and beyond with zero velocity lithospheric boundary conditions; that is no rifting. The effect of rifting is then tested by ‘restarting’ the models from output files generated during the model computation. In restarts the applied horizontal

IP

T

uniform boundary velocities are increased from 0 cm a-1 to +/- 0.5 cm a-1. By performing a series of

CR

these rifting tests for two reference standard lithospheres, SL1 and SL2, juxtaposed against the reference craton, C1 (combinations SL1-C1, SL2-C1, Section 6.1, Table 2), we determine if the

US

localization of rifting relocates from the initially weaker standard lithosphere to the craton. This approach is to understand how natural systems may behave if lithospheric extension starts during

AN

metasomatism or after it has finished. To investigate the effect of crustal strength we also rift test the

M

combination SL2-C2 in which the rheological composition of both crustal regions is the same (Section

ED

6.1, Table 2). The times of rift tests, reported with respect to start of metasomatism, are chosen to correspond to maximum CML weakening, including factors such as the time the Moho reaches its

PT

maximum temperature, or when the CML first reaches its final reduced thickness (Table 3).

CE

Decompression melts are not added as cratonic lithosphere thins during rifting because the models are designed to test whether weakening by precursor metasomatism is sufficient. The vertically integrated

AC

horizontal tensile tectonic boundary forces on the lithosphere are monitored during the rift tests and were found to have maximum values of 900-2400 GNm-1, less than or equal to estimates of ridge push (2500 GNm-1) or slab pull (>3000 GNm-1).

6. Results: Simultaneous metasomatism of the craton mantle lithosphere M-type models (Table 2) investigate the effects of simultaneous metasomatism of the whole CML (all 11 layers) over 3 Ma (M3), 10 Ma (M10), 30 Ma (M30), and 50 Ma (M50) intervals. All models 15

ACCEPTED MANUSCRIPT receive the same 30% total melt volume in the center of the craton but the volume decreases toward the edge of the craton as shown in Figure 1. The volume of the CML is increased by the volume of melt that is added. M-type models are the simplest archetypes and are used to establish the basic behavior. We discuss M30 as an example, followed by a comparison with the others. The primary difference

IP

T

among these models is the rate at which the density and volume of the CML increase and viscosity

CR

decreases, versus the time available for the heat from the melt to diffuse.

In M30, by the time that 25% melt has been added (25 Ma) a gravitational instability develops

US

at the bottom left edge of the metasomatically thickened CML (Figure 2a). CML removal caused by this

AN

instability thins approximately half of the craton by the end of metasomatism (30 Ma, Figure 2b), and propagates rapidly (Figure 2c), leaving a residual thickness of 120 km.

M

Melt heating increases the Moho temperature from 450°C to 600°C by 30 Ma (Figure 3c).

ED

Sufficient time has elapsed for heat to diffuse upward into the crust increasing the temperature at 25

PT

km depth from 345°C to 440°C. However, the 30 Ma timescale is not long enough for significant heat loss through the surface. The evolving temperature exhibits a maximum increase at 160 km and a

CE

minor increase in the surface geothermal gradient beginning at 25 Ma and peaking at 37 Ma (Figure

AC

3c). The convective removal of lower CML affects the temperature of the lower 20 km of the residual lithosphere (~100-120 km depth) In models M3, M10 and M50 (Figure 4), thinning and removal of the lower CML evolves in a similar manner to that in M30 (Figure 2), except that initiation of the instability scales approximately with the duration of metasomatism. In M3 and M10 they initiate at the end of metasomatism or soon after. In contrast, in M50 the instability develops 5 Ma before the end of metasomatism and approximately half of the CML has delaminated by the end of metasomatism (Figure 4e). 16

ACCEPTED MANUSCRIPT The Moho temperature at the end of metasomatism is dictated by heat diffusion during metasomatism (Figure 3). Accordingly, the shortest interval of metasomatism (M3) achieves the highest Moho temperature of all M-type models as soon as metasomatism finishes (635°C; Figure 3a). The crustal geothermal gradient only starts increasing at this time and reaches its maximum by 10 Ma.

IP

T

In contrast, in M50 the crustal geothermal gradient has already reached its maximum by 50 Ma (Figure

CR

3e). The temperature profiles of M30 and M50 have a similar character at the end of and after metasomatism (Figure 3d, e) when scaled by the duration of metasomatism. Thinning of the CML is

US

indicated by the sharp increase to 1350°C in the M30 and M50 profiles at 140 km depth. The deeper temperature decrease (250 km) stems from delamination and sinking of the cold CML. For all M-type

AN

models, when the CML has thermally stabilized after delamination, the thermal depth profile

M

resembles the starting one, but with a vertically condensed distribution (Figure 3).

ED

6.1 M-series rifting tests

PT

M-type models were subjected to rifting tests using two different strengths for the standard lithosphere juxtaposed against the reference craton; combinations SL1-C1 and SL2-C1. SL1 has crustal

CE

and mantle rheologies described by WQz x 3 and WOl x 5, whereas SL2 has WQz x 4 and WOl x 4 crust

AC

and mantle rheologies. The reference craton, C1, has WQz x 5 crust and WOL x 5 initial mantle rheologies (Table 2). In Table 3 the columns on the right list the results of the rifting tests in the following way: time(S) or time (C) for rifting in the standard (S) or cratonic (C) lithosphere, respectively. We also rift tested the combination SL2-C2 where the C2 cratonic lithosphere has crustal WQz x 4 and mantle WOl x 5 (Table 2, Table 3), thereby making the S and C crusts the same composition, but not the same temperature.

17

ACCEPTED MANUSCRIPT Model M30 rifts in the standard lithosphere in all tests. SL2-C1 exhibits typical results, shown for SL2-C1 at 30 and 37 Ma (Figure 5a-d) with high SL strain rates resulting in necking and final breakup. For SL1-C1 and SL2-C1, all M-type models localized extension in the standard lithosphere except

IP

T

M3 (Table 3). This establishes that the standard lithosphere cannot be significantly weaker than the

CR

initial cratonic lithosphere if rifting is to localize in the craton for the parameter values used in these models. We choose to focus on results using SL2, unless otherwise indicated. The M-type results for

US

SL2-C1 and SL2-C2 are the same in regard to rift localization (Table 3).

AN

In M3 at 3 Ma, immediately after metasomatism ends, the Moho temperature is higher in the craton lithosphere than in the standard lithosphere (635°C vs. 600°C, Figure 3). Although the cratonic

M

lithosphere has only just started to delaminate (Figure 4a), the refertilization combined with increased

ED

temperature weakens the cratonic lithosphere enough to accommodate rifting in SL2-C1 (Figure 5e, f).

PT

Additional cratonic lithosphere thinning to ~100km by 10 Ma (Figure 4b) allows even higher strain rates

CE

in the cratonic lithosphere, also resulting in craton rifting (Figure 5g, h). In the other M-type models with more protracted metasomatism, refertilization effects are the

AC

same as M3 at the end of metasomatism, but the melt heat has had longer to diffuse, resulting in a cooler Moho (Figure 3). This reduction in temperature is sufficient for rifting to localize in the standard lithosphere in all of these cases. In summary, the M series results (Table 3) show that rifting localizes in the craton only when its Moho temperature is higher than in the standard lithosphere, or where the craton has thinned substantially and the Moho is also hot.

18

ACCEPTED MANUSCRIPT 7. Results: Propagating reaction front models It is improbable that the whole mantle lithosphere gets metasomatized at the same time and at the same rate as assumed in the M-type models. A closer approximation to progressive upward migration of metasomatism is the reaction front model, which we term RF, where layers are

IP

T

metasomatized one at a time, evolving from the bottom to top of the CML over 33 Ma. Metasomatism

CR

is active for 3 Ma in each layer and the layer receives the full 30% melt by volume during this time. In contrast to the comparable M30 model lower regions of the CML are maximally weakened long before

US

upward propagating metasomatism starts in the upper layers. This evolution results in the

AN

development of multiple smaller-scale convective instabilities in RF, which peel away successive layers of the CML (Figure 6a-d) more slowly than in M30 where there is only one major thinning event (Figure

M

2). The upward propagating metasomatism front also means that the critical heating and weakening of

ED

the upper CML just below the Moho does not occur until close to 33 Ma. The final thickness of the CML in RF (Figure 6d) is similar to that in M30 (Figure 2). The Moho

PT

temperature does not start increasing until after 27 Ma and has increased from 470°C to 540°C by 30

CE

Ma (Figure 6e). After metasomatism of the last layer has finished the Moho temperature reaches 665°C. This is significantly hotter than 600°C in M30, owing to additional heat diffusing from below

AC

during the 33 Myr of metasomatism. Even after the 4th layer of the CML has metasomatized completely (12 Ma) there is only a small temperature perturbation at 160-200 km depth (Figure 6e; orange line). The character of this deviation marks the local increase (over a 20 km depth slice) in temperature during the metasomatism interval and a similar deviation can be seen for the 27, 30, and 33 Ma profiles, which represent the end of metasomatism for the 9th to 11th layers, respectively.

19

ACCEPTED MANUSCRIPT 7.1 RF rifting tests RF model rifting tests were performed starting at 30, 33, and 36 Ma for both SL1-C1 and SL2-C1. Only in the SL2-C1 case does rifting localize in the cratonic lithosphere (Table 3). Comparison of the effective strain rates for the SL2-C1 cases (Figure 8) shows that in all cases deformation initially

IP

T

localizes at the SL2 weak seed and no localization develops at the cratonic weak seed, although there is

CR

distributed cratonic extension in the 33 and 36 Ma cases. In the former, cratonic rifting progressively prevails (Figure 8c, d), whereas in the latter the rift in SL2 wins (Figure 8e, f). This behavior can be

US

understood in terms of the temperature in the cratonic lithosphere. In the RF model metasomatism thins the CML to 120-100 km by 30 Ma, but the Moho temperature only increases to 540°C (Figure 6e),

AN

less than that of the standard lithosphere, which probably accounts for rifting in the standard

M

lithosphere when rift tested beginning at 30 Ma (Figure 8a, b). Later, when metasomatism finishes at

ED

33 Ma, and the craton Moho temperature is higher than that in the standard lithosphere, extension leads to localization in the cratonic lithosphere (

PT

Figure 8c, d). Interestingly, only 3 Ma later at 36 Ma, the roles reverse, rifting reverts to the SL2

CE

standard lithosphere because the mantle just below Moho has now cooled (Figure 6e). The background effective stress profile plots confirm that the integrated strength of the C1 craton is lower at 33 Ma

AC

than at 30 Ma (Figure 7a, b), but increases by 36 Ma as the Moho temperature decreases. In addition to the SL2-C1 rift tests the role of crustal strength of the craton was considered in SL2-C2 rifting models at the same times (Figure 7). The C2 crust has the same rheology as SL2, WQz x 4, making it slightly weaker than the C1 crust, WQz x 5 (Table 2). The stress profiles (Figure 7c) demonstrate that late in the metasomatism phase of RF much of the residual cratonic lithospheric strength is in the mid and upper crust, primarily owing to it being significantly colder than the standard crust. This difference allows C1 (WQz x 5) to retain enough strength, even when its mantle is weakened 20

ACCEPTED MANUSCRIPT to WOl x 1 and the Moho heated to > 600 °C, to remain stronger than SL2 lithosphere with C1 (WQz x 4) crust. However, when both crusts have the same WQz x 4 rheologies in SL2-C2 the strain rate is slightly increased in the cratonic lithosphere ~ 300 km from the lithospheric boundary (Figure 7d). This allows the necking instability to grow faster in the craton than at the weak seed in the standard

IP

T

lithosphere (Wenker and Beaumont, 2016).

CR

In summary, the strength contrast between the standard lithosphere and the craton at the time of extension (Figure 7) and the growth rate that of necking instabilities combine to control the location

US

of rifting. All RF models with SL1 standard lithosphere rift in this region demonstrating that

AN

metasomatic weakening is not sufficient for craton rifting. Only when the models have an SL2 (WQz x 4 and WOl x 4) rheology is craton rifting possible. These results indicate that to rift a metasomatized

M

craton the growth rate of its necking instability must be faster than that in the adjacent standard

ED

lithosphere. We also note that adding melt underplating/metasomatism into the lower crust of the

PT

models would not modify these results significantly because this region is already very weak (Figure 7). The region that controls the rifting is the strong uppermost mantle.

CE

In the RF models we have only investigated metasomatism lasting for 33 Ma because this is

AC

probably more realistic than the 3 Ma in the M3 model. An RF model with 3 Ma of metasomatism would likely result in craton rifting as in M3 (Table 1). However, it is also clear that the RF model with 33 Ma metasomatism can result in craton rifting. 7.2 Sensitivity of reaction front models to the components of metasomatism The RF and M-type models demonstrate the combined effects of metasomatism in regard to rendering the craton vulnerable to rifting, but we are also interested in each of the three contributing effects, namely; rehydration, refertilization, and increasing temperature, acting alone. The series of 21

ACCEPTED MANUSCRIPT models M3-( ) and RF-( ) (Table 2) is used to investigate which effect dominates. These models are end members in which melt addition affects only one of temperature, density by refertilization, or hydration through the f-scaling. In model RF-T only the effect of melt heating and the associated volume increase is retained.

IP

T

The results (Figure 9a) show an increase in Moho temperature to 685°C, slightly higher than in RF. By

CR

33 Ma buoyant lower-lithosphere counterflow has started to smooth the original inclined base of the CML. In model RF-F only the refertilization is retained, which changes the density but not the viscosity.

US

The instability therefore takes a long time (100 Ma) to develop (Figure 9b). In model RF-H only the

AN

rehydration is included. The results (Figure 9c) are similar to RF-T because the CML viscosity is reduced in both cases, resulting in lower-lithosphere counterflow. Finally, Model RF-H&F demonstrates the

M

combined effect of rehydration and refertilization. As the CML becomes denser and less viscous small

ED

Rayleigh-Taylor instabilities form (Figure 9d) similar to the RF model (Figure 6d), but the extent of the

PT

removal is much less. Only 50 km of the lower CML is removed because the colder upper CML remains stable, demonstrating the importance of heating by the melt in RF.

CE

A subset of the RF-( ) models was rift tested at t = 33 Ma using the SL1-C1, SL2-C1 and SL2-C2

AC

combinations (Table 3). These were chosen on the basis that they were the most likely to rift in the craton. Neither hydration nor refertilization acting alone or in combination renders the craton vulnerable to rifting and the craton remains protected by the standard lithosphere. The temperature increase in RF-T does weaken the cratonic lithosphere sufficiently that the initial strain rate is higher in the cratonic mantle lithosphere in SL2-C1 and SL2-C2 at 33 Ma, immediately after metasomatism ends. However, as heat diffuses the integrated strength of the cratonic lithosphere increases and rifting once again localizes in the standard lithosphere. 22

ACCEPTED MANUSCRIPT Equivalent tests were performed for select M-( )3 series models (Table 2). In all cases rifting localizes in the SL1 and SL2 standard lithosphere (Table 3). 7.3 Sensitivity to the volume of melt added Addition of 30% melt by volume in RF results in large-scale removal of the basal CML and a

IP

T

temperature increase throughout the CML. Two equivalent tests with 20% (RF-20%) and 25% (RF-25%)

CR

melt addition during 3 Ma per layer were also computed. For RF-20% the reduced melt volume results in a smaller increase in density (28.7 kg m-3 vs. 43 kg m-3), a smaller degree of rehydration, and a

US

smaller increase in temperature (Table 3). No Rayleigh-Taylor instability develops in RF-20% (Figure

AN

A2a).

When 25% melt is added, the density increase of 35.8 kg m-3 does result in a Rayleigh-Taylor

M

instability (Figure A2b). However, the convective removal is more limited than in RF with the full 30%

ED

melt addition (Figure 6c). In RF-25% the delamination front takes 70 Ma to propagate to the center of

PT

the cratonic lithosphere (Figure A2c). These results indicate that when other conditions remain the same the total amount of melt added determines the amount of convective thinning of the CML. These

CE

models will rift in the standard lithosphere because the metasomatic weakening is significantly less

8. Discussion

AC

than in RF with the 30% melt addition.

8.1 Weakening by metasomatism; Does it allow cratons to rift? The abundance of lherzolite in the lower mantle lithosphere of many cratons has recently been interpreted as evidence of widespread silicate melt metasomatism of depleted harzburgites (Carlson et al., 2005; Griffin et al., 2009, 2003; Lee et al., 2011; Tang et al., 2013a). Here we suggest that in nature basalt-melt metasomatism weakens cratonic mantle by increasing CML density and temperature, and 23

ACCEPTED MANUSCRIPT in addition rehydration from the same or a different source reduces the effective viscosity. Removal of CML and its replacement by hot fertile sublithospheric mantle further enhances the weakening. The M and RF models, with strength distributions SL2-C1, illustrate that the timing of extension relative to metasomatism is very important. Heat diffusion is the controlling factor in this respect. All M-type

IP

T

models (Section 6) have a similar evolution when scaled to the duration of metasomatism. The time

CR

available for heat to diffuse is the remaining control that makes the model behaviors different. In M3 the exceptionally short duration of metasomatism precludes significant heat diffusion during

US

metasomatism and results in a generally hot upper cratonic mantle lithosphere by the time significant convective removal ends. By comparison, in M10, M30, and M50 heat diffusion into the crust has

AN

decreased the Moho temperature sufficiently so that it is no longer much hotter than the standard

M

lithosphere Moho and rifting takes place in the standard lithosphere.

ED

The RF, SL2-C1 combination, supports this interpretation. The CML has thinned to 125-150 km

PT

at 27 Ma (Figure 6c) but rifting still localizes in the standard lithosphere at 30 Ma. At 33 Ma, when the CML is fully metasomatized and the Moho temperature is maximized (Figure 6e), and the stress in the

CE

upper CML decreases (Figure 7b), rifting switches to the craton (Figure 8c, d).

AC

In all models that were metasomatized by 30% melt and included all the effects of metasomatism, the CML was thinned by ~100-150 km. This suggests that metasomatism can be responsible for thinning of the cratonic mantle lithosphere in natural settings. Rifting of the cratonic lithosphere was achieved only in a limited number of models and only when a large amount of melt and heat were added to the CML in a short/relatively short period of time (e.g. M3, RF). Sensitivity tests where the strength of the standard lithosphere is altered to decrease the strength contrast with the CML, SL2-C1 versus SL1-C1 (Table 3), illustrate that properties of any 24

ACCEPTED MANUSCRIPT lithosphere surrounding natural cratons most likely plays a major role in the amount of weakening needed to rift the craton. Cratons surrounded by similar strength lithosphere may rift following lesser amounts of metasomatism than required in the models presented here. For example, a cooler Moho, 500-550°C, makes the standard lithosphere stronger, therefore, less metasomatism of the CML would

IP

T

be required for it to rift. The results of the SL2-C2 tests show that the strength of the cratonic crust is

CR

also important, in this case with craton rifting where the craton and standard crusts have the same rheological properties but not the same temperatures.

US

Overall, rifting mostly localized in the standard lithosphere. This is true, even for models with

AN

large amounts of convective removal of the CML and with Moho temperatures similar to that in the standard lithosphere. This result underlines how difficult it is to rift an initially strong craton, and is

M

consistent with the observations that imply few cratons have rifted in the last 2-3.5 Ga.

ED

It is also possible that under some circumstances melt is added in sufficient amounts to cause

PT

melt-induced weakening of the craton mantle lithosphere (e.g. Gerya et al., 2015; Sizova et al., 2010), which will certainly assist craton rifting. We have not considered this effect because the small amounts

CE

of melt added each time step in our models cool rapidly and solidify such that the amount of melt

AC

present during any time step is negligible. This is because we consider the added melt to be distributed in the model layers, which leads to the rapid cooling, and not intruded in thick sills which would cool more slowly.

8.2 Comparison with other modeling research Many modelling studies attribute craton stability to a higher viscosity and lower density than adjacent standard lithosphere. Here we investigate the opposite problem. Even though cratons are cold, highly viscous, strong, buoyant, and surrounded by weaker lithosphere, there is evidence that 25

ACCEPTED MANUSCRIPT they rift. In the case of the North China Craton and North Atlantic Craton, rifting did not initiate until late in the Mesozoic (Menzies et al., 2007; Tappe et al., 2007), and was preceded by extensive thinning of the CML (Griffin et al.,1998; Tappe et al., 2007; Wu et al., 2006). The North China Craton was thinned at least locally by 60 - 160 km between the Ordovician and the Mesozoic (Griffin pers. comm.)

IP

T

Mantle plume-continent interaction leading to continental breakup has been well studied with

CR

geodynamical models (review in Buiter and Torsvik, 2014), but only a few models consider plumes as a mechanism to thin the CML (Guillou-Frottier et al., 2012; Wang et al., 2015). Wang et al. (2015) results

US

suggest that plume thinning is limited when a chemically distinct, strong, and depleted cratonic root is

AN

present. In contrast, when rheological weakening and a decrease in the buoyancy of the cratonic root are included, representing the effects of melt metasomatism, more dramatic thinning is seen,

M

especially around the edges of the craton. However, their models assume the effects of

ED

metasomatism, and do not include its progression, nor the effect of heating by melt. He (2014)

PT

considers the effects of mantle wedge convection combined with peridotite-melt interaction to thin a cratonic lithosphere. Thinning of >50 km was seen in these models but they lack tests to determine

CE

whether the craton rifts. Koptev et al. (2015) model rifting of standard lithosphere bordering the

AC

Tanzanian craton above a mantle plume, but do not investigate processes that may weaken and rift the craton per se. A model that strongly contrasts with our relatively slow metasomatism of the mantle lithosphere assumes rapid addition of large amounts of melt to the Siberian craton from a plume in a few hundred thousand years (Sobolev et al., 2011). Their numerical modelling results indicate this amount of melt and associated thermomechanical diapirism can destroy the craton. This model represents an end member in which melt was delivered in a shorter timespan than our fastest model, M3. 26

ACCEPTED MANUSCRIPT A lateral boundary between strong, cratonic, and weaker lithosphere is included in some analogue (Chemenda et al., 2002; Corti et al., 2013) and numerical models (Huerta and Harry 2007, Corti et al. 2007, Gorczyk et al. 2013, Wenker and Beaumont ,2016.). These studies agree that rifting of laterally heterogeneous lithosphere localizes in the weaker lithosphere (Corti et al., 2013), with the

IP

T

possible exception when the two lithospheres have similar integrated strengths (Corti et al., 2013;

CR

Wenker and Beaumont, 2016). None of these studies includes any weakening effect in the CML. Our modeling approach is novel in that we consider: 1) the heat and mass balance of melt that

US

is progressively added, with associated volume change, heating, rehydration and refertilization; 2) the

AN

effect of applying extensional boundary conditions during and after metasomatism; 3) the competition between rifting in adjacent standard lithosphere versus the craton and the effect of embedded weak

M

regions.

ED

8.3 Comparison with natural examples

PT

The presence of thin (80 km thick) ‘oceanic’-type lithospheric mantle is recognized under the eastern North China Craton, but the way the CML was removed is still debated (Menzies et al., 2007).

CE

Alteration of the CML was likely linked to tectonics including the North and South China Blocks collision

AC

(Late Permian-Triassic), southeast-directed subduction followed by collision of the Yangtze block with the North China Craton, and recent Pacific plates subduction (Tian et al., 2009; Wei et al., 2012). More specifically, xenoliths from Ordovician kimberlites in the eastern North China Craton indicate that cratonic depleted mantle lithosphere existed to a depth of 180 km and that heat flow was ~40 mW m-2 (Griffin et al., 1998; Menzies et al., 2007). In contrast, Mesozoic and Cenozoic xenoliths and volcanics suggest a more fertile melt source (Menzies and Xu, 1998). The timing of CML removal with respect to extension is not well constrained. Magmatism spanned a 100 Ma time interval with two 27

ACCEPTED MANUSCRIPT main episodes (190-155 Ma and 135-115 Ma), with the main phase of craton thinning and rifting probably related to the second phase (He, 2015; Wu et al., 2006; Xu et al., 2009). Heat flow peaked in the early Cenozoic at 80 mW m-2. Our results predict a lag between CML heating by melt infiltration and increased surface heat flow, which, if correct, would be consistent with melt infiltration during the

IP

T

Cretaceous.

CR

That the North China Craton ‘rifted’ is evident from extensional basins overlying the craton (Menzies et al., 2007; Menzies and Xu, 1998), with Triassic-Late Jurassic compression changing to Early

US

to mid-Cretaceous tension. Thus, North China craton rifting and the evidence for metasomatism (e.g.

AN

Zheng et al., 2015) is consistent with our proposed need for metasomatism but does not prove it. Other examples of thinned and rifted cratons (Wang et al., 2016) are the North Atlantic Craton

M

(Tappe et al., 2007), which is crosscut by the Labrador Sea, and the Tanzanian Craton, which has been

ED

metasomatically altered but only recently started to deform (Koornneef et al., 2009; Le Gall et al.,

PT

2008).

CE

That large volumes of melt can be injected into cratonic lithosphere receives support from observations from the Main Ethiopian Rift and particularly the adjacent ‘cratonic’ Ethiopian Plateau

AC

which both have low seismic velocity lower crust and uppermost mantle (e.g. Keranen et al., 2009) attributed to the presence of melt, as are the corresponding inferred high temperatures. The volume of melt inferred depends on the geometry of the melt pockets but may be as high as several percent (e.g. Hammond and Kendall, 2016). Widespread Oligocene flood basalts on the Ethiopian Plateau also indicate that melt injection started 20 Ma before the onset of rifting.

28

ACCEPTED MANUSCRIPT 8.4 Can large-volume metasomatism develop naturally? In the models we have considered large-volume metasomatism with a short duration of up to 3050 Myr that achieves both rehydration and refertilization. This rate corresponds to a vertical net flux of melt of approximately 2 x 10-3 m a-1 per unit length of the base of the craton, a total volumetric rate

IP

T

2000 m3 a-1 per unit strike length for the whole of the craton. This rate is approximately a factor of

CR

three larger than the equivalent volumetric rate of oceanic crust produced at fast spreading midoceanic ridges and is therefore not as large as it might seem. It is comparable to that supplied to the

US

crust of large igneous provinces (LIP), for example in oceanic settings where the crust is three or more times thicker than standard oceanic crust, but is larger than arc crustal growth rates above mantle

AN

wedges which range from 20 to 200 m3 a-1 (values summarized in Gerya and Meilick, 2011). The rate is,

M

however, less than that estimated for the Siberian traps (Sobolev et al., 2011) where several million

ED

cubic kilometers of magma were delivered in less than a million years. The models also include metasomatism of the shallow mantle to a degree that is larger than

PT

observed in the chemical tomography of cratons (e.g. O’Reilly et al., 2001). Clearly if our type of models

CE

were to mimic observed chemical tomography, metasomatism will deliver insufficient weakening to the upper mantle lithosphere for craton rifting. We propose that our large-volume metasomatism may

AC

develop in two phases. The RF model (Figure 6) shows phase 1 instabilities beginning at 12 Ma that thin the 250 km thick lithosphere to ~190 km and finally thinning to 150 km in places by 24 Ma. Such thinning will render the residual craton vulnerable to phase 2, which is metasomatism by fluids including basaltic and other silicic melts, seen as metasomatising fluids in off-craton styles of metasomatism where the continental lithosphere has a normal thickness. Given the availability of sufficient melts, phase 2 melt metasomatism is expected to propagate to higher levels than currently observed in cratons. Even if this metasomatism has a small effect on density it will weaken the residual 29

ACCEPTED MANUSCRIPT craton root by heating and rehydration. We suggest preserved cratons are still in phase 1. The difference between these cratons and the eastern North China Craton, for example, appears to be the need for larger volumes of melt metasomatism and perhaps rehydration to complete phase 1 and continue through phase 2 as in the RF model.

IP

T

What natural setting may produce the additional melt and water? On the basis of the

CR

‘headspace’ argument (Zheng et al., 2015) it is unlikely that sufficient anhydrous decompression melting will develop below cratons even in the presence of a mantle plume. Zheng et al. (2015), for

US

example, estimate only small volume melt will develop where the craton is thicker than 230 km, certainly insufficient for the volume required by our models. However, off-craton plumes associated

AN

with major flood basalts apparently metasomatise craton margin lithosphere (e.g. Griffin et al., 2005;

M

Sobolev et al., 2011). Fluid release from subducted slabs at active continental margins produces mantle

ED

wedge melting. However, in the case of an overlying thick craton fluid release needs to be at 6 GPa which could derive from antigorite or the 10 Angstrom phase if the subduction velocity is high and the

PT

subducted lithosphere is cold (Currie and Beaumont, 2011 and references therein), but the associated

CE

melt volumes are probably too small. The most likely natural setting for rapid, large-volume metasomatism is for cratons located above the mantle transition zone (MTZ) of a ‘Big Mantle Wedge’

AC

(BMW) (He, 2014; Zhao et al., 2009), particularly where subducted oceanic slabs have stagnated, flattened and stacked near the upper to lower mantle boundary. Water release from the MTZ (e.g. Richard et al., 2006; Richard and Bercovici, 2009; Sheng et al., 2016) as modelled by Wang et al. (2016), results in hydrous material that either plumes or percolates upward into the upper mantle beneath the craton, thereby lowering the melting point and enhancing melting as water is transferred to adjacent rocks, and lowering the effective viscosity. The hydrous basaltic melts pond at the base of the lithosphere causing weakening and then, under our assumptions, would ascend into the lithosphere by 30

ACCEPTED MANUSCRIPT diking or percolation. This mechanism, termed hydroweakening, provides both water and melts (Wang et al., 2016). Hydroweakening has recently gained favor as a potential mechanism for thinning the cratonic lithosphere beneath the eastern North China craton (see for example review Wang et al., 2016; Windley et al., 2010; Xia et al., 2013 and references therein). Although there are many

IP

T

uncertainties in regard to the Wang et al. (2016) models, the results are roughly consistent with our RF

CR

model (Figure 6) in that 2-6% partial melting in the asthenosphere leads to 1-3% partial melt in the lithosphere at any given time and thinning of the lithosphere from 200 to ~100 km in approximately 30

US

Myr after the onset of hydroweakening. Although evaluating the BMW hydroweakening mechanism is in its infancy, we suggest it holds the best potential for explaining rare large-volume widespread rapid

AN

metasomatism of cratonic roots, sufficient to render them vulnerable to rifting. It is also likely that

M

some degree of refertilization and rehydration has accumulated in cratonic mantles over time as is

ED

seen in cratons that have not been destroyed, i.e. they are still in phase 1, as discussed above. Such partially metasomatized cratons would require a smaller ‘large-volume metasomatic event ’of the type

PT

we propose to render them vulnerable to rifting.

CE

8.5 Metasomatism and intracratonic basins

AC

We note that melt metasomatism also affects the isostatic balance of CML by increasing its density. This can lead to subsidence and formation of intracratonic basins (Wenker and Beaumont, in prep.). Increased density alone produces subsidence, but thermal expansion from melt-heating counteracts subsidence and can result in temporary surface uplift followed by subsidence on cooling. Dynamical uplift above a mantle plume similarly uplifts the surface. In addition, large amounts of metasomatism produce lithospheric instabilities as in the M and RF models, resulting in a phase of isostatic uplift during the replacement of lower lithosphere by hot, less dense sublithospheric mantle. 31

ACCEPTED MANUSCRIPT We mention these effects here because initiation, subsidence, and subsequent uplift of intracratonic basins prior to rifting of the craton may be key observational evidence that the CML is being metasomatized and by how much.

T

9. Conclusions

IP

We have used 2D numerical models to investigate the effects of metasomatic weakening of

CR

cratons and whether this increases their propensity to rift. Weakening can occur through refertilization (increase in density) and rehydration (decrease in viscosity), but temperature increases from local

US

metasomatic melt cooling and solidification can also be an important model component in regard to

AN

decreasing the strength of cratonic mantle lithosphere (CML).

M

The model results indicate that significant melt infiltration and reaction is needed for cratons to rift. Even after melt metasomatism has heated the CML so that the Moho temperature is higher than in

ED

the adjacent standard lithosphere, and has thinned the CML by 100-120 km, rifting still localized in the

PT

standard lithosphere in most of our SL2-C1 models. This result is consistent with the long-time survival of most cratons; they are difficult to rift. Attempted craton rifting is, however, most likely to succeed

CE

when the temperature of the CML in the vicinity of the Moho is a maximum. At this time the strength

AC

of cratonic crust versus that of surrounding lithosphere is important because much of the residual lithospheric strength remains in the mid and upper crust following extensive CML metasomatism. In summary, given that most, if not all, cratons have been subjected to rifting forces this means they are capable of resisting rifting. What we have shown is that their ability to resist rifting is primarily the result of protection afforded by surrounding standard, but weaker, lithosphere. This is a point not commonly considered in association with rifting. It follows that the strength of the craton must be reduced such that it is weaker than the standard lithosphere if it is to rift. If this weakening is the result 32

ACCEPTED MANUSCRIPT of metasomatism it requires large volume melt metasomatism. Such large volume metasomatism is probably extremely rare, and not just the result of a mantle plume or location above a normal mantle wedge. We suggest that a ‘Big Mantle Wedge’ may be capable of large volume metasomatism, but admit that more research is needed on the destabilization of water in the mantle transition zone in

IP

T

association with a big mantle wedge.

CR

Refertilization densification also causes isostatic subsidence of cratons thereby forming intracratonic basins. Conversely, uplift and erosion followed by protracted subsidence, as identified by

US

unconformities in some cratonic basins, may be explained by removal of lithospheric mantle and

AN

subsequent cooling of hot replacement sublithospheric mantle. In summary, we suggest subsidence and uplift of cratons may be one key indicator of CML metasomatism which, if sufficiently intense, can

AC

CE

PT

ED

M

lead to craton rifting provided rifting occurs when the CML is weak and vulnerable.

33

ACCEPTED MANUSCRIPT Appendix A Thermo-Mechanical Modeling Methods A.1 Overview of SOPALE-Nested software We use a 2D upper-mantle-scale model comprising adjacent regions of standard and cratonic lithosphere subject to extension under specified velocity boundary conditions following simulated

IP

T

metasomatism of the cratonic mantle lithosphere. We calculate the model evolution using an

CR

Arbitrary-Lagrangian-Eulerian (ALE) finite element method using the software SOPALE-nested. It solves the thermo-mechanically coupled, incompressible viscous-plastic creeping (Stokes) flow equations and

US

has frictional-plastic and thermally activated power-law viscous rheologies.

AN

For each time step, the force balance equations, for an incompressible creeping flow, are solved on a 2D Eulerian grid (Equations A.1 and A.2). These are coupled to the energy balance

M

equation (Equation A.3) through the temperature-dependent viscosity and density and constrained by

AC

CE

PT

ED

mechanical and thermal boundary conditions.

where σij is the deviatoric stress tensor, density, heat,

spatial coordinates,

the gravitational acceleration acting vertically, the absolute temperature, time,

production per unit volume, and

is the pressure (mean stress),

a component of velocity,

the thermal conductivity,

the

the specific

radioactive heat

the volumetric thermal expansivity. We use an extended

Boussinesq approximation that includes adiabatic heating (last term Eqn A.3) but we purposely omit 34

ACCEPTED MANUSCRIPT shear heating because we do not want to include its effects in strain localization during rifting. This also means that the STP densities we quote (Table 1) only vary with thermal expansion and not pressure. Models have a high resolution nest embedded in the large scale domain (Figure 1). The solution

IP

T

for the velocity and thermal fields are first found for the 2200 x 600 km domain. The solution is then

CR

used to set boundary condition for the nested domain in the center of the model from 100-2100 km and 360 km deep. A two-way connection is maintained by the Lagrangian tracking particles, which are

US

shared between the two solutions. Particles in the nested domain always obey the higher resolution

AN

solution.

M

A.2 Model Design

ED

The large-scale model has basal and side boundaries with free slip in the horizontal and vertical directions, respectively, and no basal vertical velocity. Extension during rift tests is defined by

). A small horizontal flux of material into the sub-lithospheric mantle allows for conservation of mass,

CE

1

PT

horizontal velocities of +/- 0.5 cm a-1 on the right and left vertical lithospheric boundaries (total 1 cm a-

and there is no material flux through the base. The model (Figure 1) juxtaposes ‘standard‘ 140 km

AC

thick lithosphere (SL1 or SL2) with a 250 km thick ‘cratonic’ lithosphere (C1 or C2) (Table 1) (Artemieva and Mooney, 2001; Griffin et al., 2009). Two 12 x 5 km Mantle Weak Zones (MWZ) are embedded in the upper mantle lithosphere equidistant from the lithospheric boundary. The MWZ are 500 km away from the lithospheric boundary and have a reduced internal angle of friction of 2˚ to represent the effect of inherited plastic strain softening. The size and position of the weak zones was chosen to represent simplified residual mantle weaknesses; elsewhere, annealing in the hotter viscous parts of the lithosphere (Yamasaki et al., 2006), where grain regrowth is efficient, is assumed to have removed 35

ACCEPTED MANUSCRIPT strength contrasts associated with grain size reduction in lithospheric-scale sutures, although some form of heterogeneity (e.g. compositional) may remain (Vauchez et al., 2012).

The base of the cratonic mantle lithosphere (CML) is tapered over a horizontal distance of 300

T

km to represent the effect of a history of destabilization of the edge of the craton by various processes

IP

and to reduce the tendency for vigorous edge driven convection (EDC), which is not the subject of this

CR

investigation. In addition, the models are computed for an initial spin-up phase of 1 Ma to achieve an

US

isostatic equilibrium between the lithospheres and the sublithospheric mantle. The model was also tested for an initial 100 Ma to demonstrate stability against EDC and convective CML removal (Figure

AN

1). There is only minimal deformation of the lower boundary of the lithospheres during this phase

M

which demonstrates initial stability against EDC erosion of the CML owing to disequilibrium in the initial conditions, a requirement emphasized by Sleep (2007). Destabilization of the CML during the

ED

main phase of the computations is therefore the result of the metasomatism, not the initial conditions.

PT

A.3 Material Properties

CE

The deformation mechanism is determined by the state of stress in the flow. Flow is viscous (Equation A.1) where stress is lower than the frictional-plastic yield stress. In contrast, where stress is

AC

on the frictional-plastic yield envelope, described by the pressure-dependent Drucker-Prager yield criterion (Equation 5), material deforms in a plastic (stiff) manner.

In power-law flow (Equation A.1), f is the viscosity scaling factor which modifies the relative strength of viscous materials with respect to the reference flow law, f =1, as explained below,

is the pre36

ACCEPTED MANUSCRIPT exponential scaling factor, is the power law exponent,

gas constant. In Equation (2),

is the activation energy, =

invariant of the deviatoric stress

is the activation volume, R is the universal

, the effective stress, is the square root of the second ,

is the cohesion,

is the effective internal angle of friction

T

,

is the second invariant of the deviatoric strain rate tensor

`~ 15˚

where

where

is the

). Parameter values are given in Table 1.

CR

pore-fluid pressure), which gives

IP

taking account of near-hydrostatic fluid pressure (using

US

Strain-softening is included in the brittle regime to represent the formation of weak faults. The effective internal angle of friction is decreased linearly from 15˚ to 2˚ as the effective strain increases

AN

from 50% to 150%. The models are not particularly sensitive to the range of strain over which materials

M

undergo strain softening (Huismans and Beaumont, 2007). The requirement for a low strain-softened angle of internal of friction is a consequence of the resolution of the computational finite element grid,

PT

A.4 Thermal properties

ED

such that with higher resolution grid a larger value would give equivalent results.

CE

The initial 2D steady-state temperature of the system is determined by the insulated side boundaries, the basal heat flux, the radioactive heat production in the crust (Ar), and the thermal

AC

conductivity (K(TK)). These parameters are different for the standard and cratonic lithospheres (Table 1) and are selected so that the temperature at their base is 1350˚C. The standard lithosphere has a higher geothermal gradient (surface heat flow is 62.5 mW m-2) than the cratonic lithosphere (surface heat flow of 38.6 mW m-2). The basal heat flux changes from 20 to 22 mW m-2 at the center of the model. The standard lithosphere has a radioactive heat production of 1.5 µW m -3 in the upper crust and 0.5 µW m-3 in the lower crust (Hasterok and Chapman, 2011) and the Moho temperature is 600˚C.

37

ACCEPTED MANUSCRIPT The cratonic crust has a lower radiogenic heat production of 0.6 µW m-3 and 0.2 µW m-3 for the upper and lower crust, respectively, and the Moho temperature is 450˚C. The choice of thermal expansivity of the mantle lithospheres and sublithospheric mantle is important because it determines the thermal buoyancy which directly competes with the effects of

IP

T

depletion on material density, particularly in the CML which is designed to be approximately isopycnic.

CR

Uniform values of thermal expansivity for the mantle in the range 3.0-3.2 x 10-5 deg-1 are commonly used, but there is increasing evidence that thermal expansivity depends on temperature (Bouhifd et

US

al., 1996). We approximate this mantle temperature dependence such that for T < 500 K α(TK) = 3.1 x

AN

10-5 K-1 and for T> 2000 K α(TK) = 3.9 x 10-5 K-1 and varies linearly between these temperatures. The high temperature expansivity has been reduced by comparison with the Bouhifd et al. (1996) value to

M

account for the secondary effect of pressure on thermal expansivity (Katsura et al., 2009). During time

ED

stepping, the heat conservation equation A.3 is solved with the same boundary conditions as those

PT

described above. The methods used for the effects of simulated metasomatism are described in the

AC

CE

main text.

38

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

AC

CE

PT

ED

M

Figure A1. Melt Fractions. (Annen et al., 2006)

39

PT

ED

M

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

AC

CE

Figure A2. Model sensitivity to the volume of melt added in Reaction Front, RF, models (Table 2) in which the CML is progressively metasomatized by addition of melt to successively shallower layers over a 33 Ma timespan . a) RF-20%, total melt added to each CML layer by 33 Ma is 20%. No CML thinning has developed by 80 Ma. b,c) RF-25%, total melt added to each CML layer by 33 Ma is 25% in. A Rayleigh-Taylor instability has started, but by 70 Ma the CML has not thinned to the same extent as in the RF model with 30% melt addition (Figure 6). Metasomatism decreases to the left of the vertical black line in the CML (see Figure 1).

Acknowledgements

This research was funded by an NSERC Discovery Grant and Canada Research Chair in Geodynamics to CB. We thank Cin-Ty A. Lee for helpful advice on the melt-peridotite reaction. We thank Douglas Guptill for his computing support in regard to the development of SOPALE-nested and the associated graphics programs and, in particular, coding of the metasomatism. Reviews by colleagues Rebecca Jamieson, Charlotte Keen, Nicholas Culshaw and Barrie Clarke are gratefully acknowledged. SOPALE was originally developed at Dalhousie University by Philippe Fullsack. We also thank the journal reviewers, W. Griffin, J. van Hunen, and

40

ACCEPTED MANUSCRIPT T. Gerya, and one anonymous reviewer, and the editor P. Agard for comments that resulted in significant

AC

CE

PT

ED

M

AN

US

CR

IP

T

improvements to the paper.

41

CE

PT

ED

M

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

AC

Figure 1. a) Model design. Standard lithosphere (left) is 140 km thick with Moho temperature 600°C. Cratonic lithosphere (right) is 250 km thick with Moho temperature 450°C. In all figures temperature contours (blue) are in °C. The boundary between the standard and cratonic lithosphere is at 1100 km. Two mantle weak zones (12 x 5 km) with identical properties are embedded at equal distances from the boundary, one in each lithosphere. The standard lithosphere (Table 2) has a crust governed by the WQz x 3 flow law (see Appendix A for explanation of f-scaling) and a mantle lithosphere rheology of WOl x 5 in the SL1 models. The SL2 models (Table 2) have WQz x 4 and WOl x 4 flow laws. Cratonic lithosphere, C1, has WQz x 5 crust and mantle lithosphere (CML) that starts at WOl x 5. As metasomatism proceeds (pink), the CML becomes more hydrated; this weakening corresponds to a linear change in model mantle lithosphere rheology to WOl x 1. Inclined base of the CML is designed to reflect any prior flow and thinning of the CML and to represent an initially stable system. Higher resolution domain is outlined by the dashed box. ML = mantle lithosphere. Temperature distributions presented in figures 3 and 6 are at x = 500 km and 1700 km. Vertical black lines in the CML at 1150 and 1350 km show region of reduced metasomatism of the CML in M and RF-type models. From 1100 to 1150 km there is no metasomatism and the amount of melt added increases linearly from zero to the maximum between 1150 and 1350 km. Note that variably metasomatised CML will be advected across these regions during deformation. b) Reference model (no metasomatism) after the initial phase of 1 Ma, this is t = 0 Ma for all the models in the results section. c) Reference model after 100 Ma. Only a small amount of buoyant counterflow has developed at the vertical lithospheric boundary. Velocities in the convecting mantle are -1 less than 1 cm a and no CML is removed by basal shear. The boundary between the standard lithosphere (SL, green) and the cratonic lithosphere (CL, blue) is shown.

42

PT

ED

M

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

AC

CE

Figure 2. Convective removal of the CML in M30 (Table 2). a) At 25 Ma the CML is starting to become unstable. b) At the end of metasomatism (30 Ma) the delamination front is propagating and thinning the craton. c) Lithospheric thinning is complete at 36 Ma and in some places leaves the craton lithosphere thinner than the standard lithosphere. Metasomatism decreases to the left of the vertical black line in the CML (see Figure 1).

43

ACCEPTED MANUSCRIPT b) M10

1500

0

50

50

100

100

150

0 Ma 3 Ma

200

25 Ma

250

17 Ma SL

d) M50

1500

0

Temperature (°C) 500 1000

1500

M

0

ED

50

150

CE

0 Ma

PT

100

30 Ma

250

37 Ma

AC

200

Depth (km)

0

Temperature (°C) 500 1000

AN

300

0

50 100 150 200 250

0 Ma 50 Ma 59 Ma SL

SL

300

10 Ma

US

300

Depth (km)

0 Ma

250

SL

c) M30

150 200

10 Ma

1500

T

0

0

Temperature (°C) 500 1000

IP

Temperature (°C) 500 1000

Depth (km)

Depth (km)

0

CR

a) M3

300

Figure 3. Temperature evolution for M-type models (Table 2). a-d) Temperature evolution at x = 1700 km for M3, M10, M30 and M50, with metasomatism spanning 3, 10, 30 and 50 Ma, respectively. Standard lithosphere (SL) temperature distribution at x = 500 km. Moho depth is indicated by the dashed black line. Note that metasomatism alone does not heat the CML to the temperature of SL. It is only when accompanied by convective removal of the lower CML that CML temperatures approach SML temperatures.

44

ED

M

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

AC

CE

PT

Figure 4. Model evolution for the cratonic parts of M3, M10, and M50 (Table 2). a) A convective instability has started to form in M3 by the end of metasomatism (3 Ma). b) Maximum thinning of the CML in M3 has finished by 10 Ma. c) At the end of metasomatism in M10, part of the CML has already delaminated. d) Thinning of the CML in M10 has reached the model boundary by 17 Ma. e) By the end of metasomatism (50 Ma) approximately half of the CML in M50 has delaminated. f) Maximum thinning in M50 is complete by 59 Ma. Metasomatism decreases to the left of the vertical black line in the CML (see Figure 1).

45

AC

CE

PT

ED

M

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

Figure 5. Rifting tests for M30 and M3 for the SL2-C1 (Table 2) combination showing effective strain rates (effective strain 1/2 rate = (second invariant deviatoric strain rate) ). t values give the start time of the rifting tests. Temperature contours in °C. a) Strain rate M30 at t1 = 30 Ma. Although the cratonic mantle lithosphere is weakened by metasomatism rifting initiates in the standard lithosphere, shown by higher strain rates. b) 50 km of extension, rifting localizes in the standard lithosphere. c) Strain rate M30 at t2 = 37 Ma. Thinning of the cratonic mantle lithosphere does not change the strain rate distribution in favor of rifting in the cratonic lithosphere. d) Rifting again localizes in the standard lithosphere. e) Strain rate M3 at t1 = 3.0 Ma (end of metasomatism), CML has not thinned by the end of metasomatism. f) The strain rate is higher in cratonic lithosphere after 50 km of extension. g,h) Strain rate M3 at t2 = 10 Ma. Heat has diffused into the crust and the cratonic lithosphere also rifts in this case. The red areas in the lower parts of the strain rate plots reflect rapid deformation at and below the delamination front. White boxes show positions of the weak seeds and M marks the initial depth to the Moho.

46

AC

CE

PT

ED

M

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

Figure 6. Reaction front, RF, metasomatism (Table 2). Layers are sequentially metasomatized for 3 Ma at a time (pink layers) over a 33 Ma interval. a) The instability starts after the 4th layer has completed metasomatism. b,c) Delamination is characterized by smaller instabilities than in the M30 model. d) Final thinning results in a 100-120 km thick CML. e) th Temperature evolution of RF at x= 1700 km. A small temperature deviation from the 0 Ma profile can be seen after the 4 layer of the CML has finished metasomatism at 12 Ma. The temperature at the Moho (black dashed line) has not discernably increased by 27 Ma and only reaches its maximum after the top layer metasomatises (33 Ma). Metasomatism decreases to the left of the vertical black line in the CML (see Figure 1).

47

ACCEPTED MANUSCRIPT Stress (MPa)

0

10

10

20

20

30

30

40

40

60

30 Ma

90

SL 33 Ma

100

AN

200

300

M

0 10

ED

20

PT

40

60

90 100

CE

30 Ma 33 Ma

AC

Depth (km)

30

50

33 Ma 36 Ma SL 33 Ma

d) Strain Rate SL2-C1 vs. SL2-C2

Stress (MPa) 100

36 Ma

5 X 10 0

-17

Strain Rate (s-1) 1.05 X 10

-15

10

CL

20

1700 km

2.05 X 10

-15

30

Depth (km)

0

80

60

80

36 Ma

c) SL2-C2

70

50

70

33 Ma

100

300

30 Ma

70

90

200

IP

50

100

T

0

300

0

80

Stress (Mpa)

b) SL2-C1

200

CR

100

Depth (km)

Depth (km)

0

US

a) SL1-C1

40 50 60 70

SL

CL 1400 km

80 90

SL 33 Ma

100

Figure 7. Background vertical effective deviatoric stress profiles for Reaction Front, RF, models (Table 2) at 500 and 1700 km ½ during rift testing (effective stress = (second invariant deviatoric stress) ) . ‘Background’ means that these profiles represent regional stress and not stress influenced by the weak seeds. Model properties are: SL1 (crust, WQz x 3, mantle WOl x 5); SL2 (crust, WQz x 4, mantle WOl x 4); C1 (crust, WQz x 5, mantle WOl x 5); C2 (crust, WQz x 4, mantle WOl x 5). Other properties are given in Table 1. a) SL1-C1 rift tested at 30, 33, 36 Ma. b) SL2-C1 rift tested at 30, 33, 36 Ma. c) SL2-C2 rift tested at 30, 33, 36 Ma. The bolder lines in a-c indicate models that rifted in the cratonic lithosphere. Note that the background integrated strength of SL is greater than CL, but SL strength is locally much less at the CML weak seed, 30 vs 220 MPa. d) Effective strain rate plots vs depth for RF extended at 36 Ma for SL2-C1 (dashed) and SL2-C2 (solid). The profiles are labeled with their locations relative to the left boundary of the model. The profiles at 1400 km are closest to the location of breakup, when breakup occurs in the CL. The strain rate differs most between SL2-C1 and SL2-C2 in this region. Although

48

ACCEPTED MANUSCRIPT

CE

PT

ED

M

AN

US

CR

IP

T

lithospheric stresses are similar between SL2-C1 and SL2-C2, the rate at which the necking instability grows in the SL2-C2 model is indicated by the higher strain rate.

AC

Figure 8. Rifting tests for Reaction Front, RF, model for the combination SL2-C1 (Table 2) showing effective strain rates 1/2 (effective strain rate = (second invariant deviatoric strain rate) ). t values give the start time of the rifting tests. Temperature contours in °C. a,b) Strain rate at 30.0 Ma. Although the CML has thinned to less than150 km and the Moho temperature has reached 540°C at 30 Ma, the strain rate is higher in the standard lithosphere. Rifting localizes in the standard lithosphere. c,d) After metasomatism of the layer immediately beneath the Moho, rifting now develops in the cratonic lithosphere. e,f) Only 3 Ma after the end of metasomatism, heat diffusion has strengthened the craton enough for rifting to shift back to the standard lithosphere. White boxes show positions of the weak seeds and M marks the initial depth to the Moho.

49

CE

PT

ED

M

AN

US

CR

IP

T

ACCEPTED MANUSCRIPT

AC

Figure 9. Sensitivity tests for the RF-( ) series (Table 2). a) RF-T. A temperature increase in the CML results in buoyant counterflow at the vertical lithospheric boundary. The Moho temperature has increased to 685°C. b) RF-F. CML is refertilized, but the viscosity does not change, delamination is delayed. c) RF-H. Hydration of the CML results in extensive counterflow. Lower layers, which are hotter than the upper layers, are incorporated into the counterflow more readily than the upper layers. d) RF-H&F. Result resembles RF model, as small parts of CML delaminate. However the final thinning is not as extensive as in RF. Metasomatism decreases to the left of the vertical black line in the CML (see Figure 1).

50

ACCEPTED MANUSCRIPT Property Frictional-plastic parameters Internal angle of friction (˚) Over strain range of (%) Crust cohesion (Pa) Mantle cohesion (Pa)

Symbol

Value

CCrust CMantle

15 → 2 50-150 7 1·10 6 2·10

Crust: Wet Quartzite Power law exponent -1 Activation energy (J mol ) -n -1 Initial constant (tensor invariant) (Pa s ) Scaling factor -3 Crustal density (kg m )

WQz nCrust Q A fCrust

Mantle: Wet Olivine Power law exponent -1 Activation energy (J mol ) 3 -1 Activation Volume (m mol ) -n -1 Initial constant (tensor invariant) (Pa s )

WOl nML = nSML Q V* A

T IP

US

CR

S

AN

Scaling factor -3 Standard mantle lithosphere density (kg m ) -3 Cratonic mantle lithosphere density (kg m ) -3 Sub-lithospheric mantle density (kg m )

CE

AC

3 3 430·10 -6 12·10 -14 1.7578·10

fMantle

3 (SL1), 4 (SL2), 5-1 (C1, C2) 3360 3335 3378

TMoho TCMoho TL q qs qc

600 450 1350 21-22 62.5 38.6

PT

ED

M

Standard Moho temperature (˚C) Cratonic Moho temperature (˚C) Base of lithosphere temperature (˚C) -2 Basal heat flux (mW m ) -2 Standard lithosphere surface heat flux (mW m ) -2 Cratonic lithosphere surface heat flux (mW m ) Radiogenic Heat Production -3 Standard upper crust (µW m ) -3 Standard lower crust (µW m ) -3 Cratonic upper crust (µW m ) -3 Cratonic lower crust (µW m ) -1 Coefficient of thermal expansion (˚C ) Over temperature range (˚C) -1 -1 Crust thermal conductivity (W m ˚C ) -1 -1 Standard mantle lithosphere thermal conductivity (W m ˚C ) -1 -1 Cratonic mantle lithosphere thermal conductivity (W m ˚C ) -1 -1 Sub-lithospheric mantle thermal conductivity (W m ˚C ) Over temperature range (˚C) Melt Metasomatism -1 -1 Heat capacity peridotite (J Kg ˚C ) -1 -1 Latent heat of crystallization peridotite (J Kg ˚C ) Peridotite liquidus, intermediate temperature and solidus (˚C) Peridotite melt coefficient high, low for Temperature of the melt (˚C) -1 -1 Heat Capacity melt (J Kg ˚C ) -1 -1 Latent heat of crystallization melt (J Kg ˚C ) Melt liquidus, intermediate temperature and solidus (˚C) Melt coefficient high, low Table 1. Model Parameters

4 3 223·10 -28 8.574·10 3 (SL1), 4 (SL2), 5 (C1), 4 (C2) 2800

KCrust KMLS(T) KMLC KSLM(T)

Cp1 L1 T1L, T1i, T1S, X1h, X1l Tm Cp2 L2 T2L, T2i, T2S, X2h, X2l

1.5 0.5 0.6 0.2 -5 -5 3.2·10 → 3.9·10 500-2000 2.25 2.7 5.3 3-40 1323-1423 1100 5 4.0·10 1900, 1750, 1600 0.6, 04 1300 1480 5 4.0·10 1250, 1075, 725 0.6, 04

51

ACCEPTED MANUSCRIPT

CR

Sensitivity tests for components of metasomatism

T

RF-Type

Properties Melt addition of 30% added to all CML layers over a specified length of time (with reduced amounts at the edge of the craton). Models are named by the length of time specified to add the melt, M3 (3 Ma), M10 (10 Ma), M30 (30 Ma) and M50 (50 Ma). Progressive upward addition of melt to each CML layer (with reduced amounts at the edge of the craton) (Figure 1) such that 30% melt addition takes 3 Ma for each layer and the total time for complete melt addition is 33 Ma (RF-33). Models are named by the volumetric amount of melt added: RF (30%), RF-25% (25%), RF20% (20%). M and RF-type models that test the component effects of metasomatism are named according to the active component: M-T (temperature effect), M-F(refertilization effect), M-H (rehydration effect) and M-H&F (rehydration and refertilization effects). Same notation for RF models but with prefix RF.

IP

Model Type M-Type

Model material properties during rift tests. Rift tests were conducted in which the standard (S) and cratonic (C) lithospheres were assigned specific power-law flow properties determined by the f-scaling factor (Section 3) for the wet quartz (WQz) and wet olivine (WOl) flow laws. All other properties do not vary among models and are given in Table 1. Section 4 has an explanation of f-scaling and the reasons for choosing the f-values. SL1 Standard Lithosphere. Crust WQz, f = 3, mantle lithosphere, WOl, f = 3 SL2 Standard Lithosphere. Crust WQz, f = 4, mantle lithosphere, WOl, f = 4 C1 Cratonic Crust. Crust WQz, f = 5, mantle lithosphere, WOl, f = 5 to 1 (change results from metasomatism) C1 Cratonic Crust. Crust WQz, f = 4, mantle lithosphere, WOl, f = 5 to 1 (change results from metasomatism) Table 2. Summary of model properties that vary among models

AC

CE

PT

ED

M

AN

US

Power Law Properties

52

ACCEPTED MANUSCRIPT Model

M-Models M3 M10 M30 M50 Sensitivity tests M3-T M3-H M3-F M3-H&F RF-Models RF Sensitivity tests RF-T RF-H RF-F RF-H&F RF-25% RF-20%

Max Moho Temp. (°C)

Time of Max thinning (Ma)

Moho Temp. at Max Thinning (°C)

SL1-C1

635 625 600 580

10 17 37 59

605 595 585 565

3(S) 10(S) 30(S) 50(S)

630 450 450 450

-

-

4(S) 20(S)

665

28

450

30(S)

685 450 450 450 640 615

100* 100* 60 -

450 450 550 -

E C

T P

D E

6(C) 16(S) 37(S) 55(S)

Rifting tests and position of localization Ma (S vs C) SL2-C1 10(C)

3(C) 10(S) 30(S) 50(S)

6(C) 16(S) 37(S) 55(S)

T P

10(C)

I R

C S

4(S) -

33(S)

M

A

33(S) 33(S) 33(S) 33(S)

U N

36(S)

30(S)

33(C) 33(S) 33(S) 33(S)

36(S)

3(C) 10(S) 30(S) 50(S)

SL2-C2 6(C) 16(S) 37(S) 55(S)

10(C)

33(C)

36(C)

4(S) 30(S)

33(S) -

-

Table 3. Summary of model results. M-models have simultaneous metasomatism throughout the cratonic mantle lithosphere. RF (Reaction Front) models have metasomatism that propagates upward from the base of the cratonic lithosphere. Max Moho temperature, time of max thinning, and Moho temperature at max thinning refer to conditions during metasomatism without rift testing. * Models were computed to 100 Ma, but more thinning is expected after this time (Figure 9). In the rifting test columns, the number is the time of testing (Ma) and S and C stand for localization in the standard and cratonic lithospheres, respectively. SL1 models have a standard lithosphere with WQz x 3 and WOl x 5 for the crust and mantle lithosphere respectively. SL2 models have standard lithosphere with WQz x 4 and WOl x 4 for the crust and mantle lithosphere. C1 and C2 models have cratonic lithosphere with crust WQz x 5 and WQz x 4, respectively and mantle WOl x 5. The models presented in the results section are of the SL2 type unless otherwise indicated. Where rift test results are not shown models are predicted to rift in the standard lithosphere on the basis of other rift test results.

C A

53

ACCEPTED MANUSCRIPT

Need to add references from Response to Reviewers Reference for mantle wedge magma production Gerya, T.V., and Meilick, F.I., 2011. Geodynamic regimes of subduction under active margins:

T

effects of rheological weakening by fluids and melts. J. Metamorphic Geol. 29, 7-31.

CR

IP

doi:10.1111/j.1525-1314.2010.00904.x

There is also a typo in He 2014 ref …should be big (not bit) mantle wedge

US

References

AN

Annen, C., Blundy, J.D., Sparks, R.S.J., 2006. The genesis of intermediate and silicic magmas in

M

deep crustal hot zones. J. Pet. 47, 505–539. doi:10.1093/petrology/egi084

ED

Artemieva, I.M., Mooney, W., 2001. Thermal thickness and evolution of Precambrian lithosphere: A global study. J. Geophys. Res. 106, 16387–16414.

PT

Aulbach, S., Rudnick, R.L., McDonough, W.F., 2008. Li-Sr-Nd isotope signatures of the plume

CE

and cratonic lithospheric mantle beneath the margin of the rifted Tanzanian craton

AC

(Labait). Contrib. to Mineral. Pet. 155, 79–92. doi:10.1007/s00410-007-0226-4 Bashir, L., Gao, S.S., Liu, K.H., Mickus, K., 2011. Crustal structure and evolution beneath the Colorado Plateau and the southern Basin and Range Province: Results from receiver function and gravity studies. Geochemistry Geophysics Geosystems 12, Q06008. doi:10.1029/2011GC003563 Beaumont, C., Jamieson, R.A., Butler, J.P., Warren, C.J., 2009. Crustal structure: A key constraint

54

ACCEPTED MANUSCRIPT on the mechanism of ultra-high-pressure rock exhumation. Earth Planet. Sci. Lett. 287, 116–129. doi:10.1016/j.epsl.2009.08.001 Beaumont, C., Nguyen, M.H., Jamieson, R.A., Ellis, S., 2006. Crustal flow modes in large hot orogens, in: Law, R.D., Searle, M.P., Godin, L. (Eds.), Channel Flow, Ductile Extrusion and

IP

T

Exhumation in Continental Collision Zones. Geological Society, London, Special

CR

Publications, pp. 91–145. doi:10.1144/GSL.SP.2006.268.01.05

Bernstein, S., Kelemen, P.B., Hanghøj, K., 2007. Consistent olivine Mg# in cratonic mantle

US

reflects Archean mantle melting to the exhauston of orthopyroxene. Geology 35, 459–462.

AN

doi:10.1130/G23336A.1

M

Beuchert, M.J., Podladchikov, Y.Y., 2010. Viscoelastic mantle convection and lithospheric

ED

stresses. Geophys. J. Int. 183, 35–63. doi:10.1111/j.1365-246X.2010.04708.x Bodinier, J.L., Vasseur, G., Vernieres, J., Dupuy, C., Fabries, J., 1990. Mechanisms of mantle

PT

metasomatism: Geochemical evidence from the Lherz peridotite. Jour. Pet. 31, 597–628.

CE

Bouhifd, M.A., Andrault, D., Fiquet, G., Richet, P., 1996. Thermal expansion of forsterite up to

AC

the melting point. Geophys. Res. Lett. 23, 1143–1146. doi:10.1029/96GL01118 Buiter, S.J.H., Torsvik, T.H., 2014. A review of Wilson Cycle plate margins: A role for mantle plumes in continental break-up along sutures? Gondwana Res. 26, 627–653. doi:10.1016/j.gr.2014.02.007 Butler, J.P., Beaumont, C., Jamieson, R.A., 2014. The Alps 2: Controls on crustal subduction and (ultra)high-pressure rock exhumation in Alpine-type orogens. J. Geophys. Res. 119, 5987–

55

ACCEPTED MANUSCRIPT 6022. doi:10.1002/2013JB010799 Carlson, R.W., Pearson, D.G., James, D.. E., 2005. Physical, chemical, and chronological characteristics of continental mantle. Rev. Geophys. 43, RG1001.

T

doi:10.1029/2004RG000156

IP

Chemenda, A., Déverchère, J., Calais, E., 2002. Three-dimensional laboratory modelling of

CR

rifting: application to the Baikal Rift, Russia. Tectonophysics 356, 253–273.

US

doi:10.1016/S0040-1951(02)00389-X

Corti, G., Iandelli, I., Cerca, M., 2013. Experimental modeling of rifting at craton margins.

AN

Geosphere 9, 138–154. doi:10.1130/GES00863.1

M

Corti, G., van Wijk, J., Cloetingh, S., Morley, C.K., 2007. Tectonic inheritance and continental rift

ED

architecture: Numerical and analogue models of the East African Rift system. Tectonics 26.

PT

doi:10.1029/2006TC002086

Currie, C.A., Beaumont, C., 2011. Are diamond-bearing Cretaceous kimberlites related to low-

CE

angle subduction beneath western North America? Earth Planet. Sci. Lett. 303, 59–70.

AC

doi:10.1016/j.epsl.2010.12.036 Doin, M.-P., Fleitout, L., Christensen, U., 1997. Mantle convection and stability of depleted and undepleted continental lithosphere. J. Geophys. Res. 102, 2771. doi:10.1029/96JB03271 Foley, S.F., 2008. Rejuvenation and erosion of the cratonic lithosphere. Nat. Geosci. 1, 503–510. doi:10.1038/ngeo261 François, T., Burov, E., Meyer, B., Agard, P., 2012. Surface topography as key constraint on 56

ACCEPTED MANUSCRIPT thermo-rheological structure of stable cratons. Tectonophysics 602, 106–123. doi:10.1016/j.tecto.2012.10.009 Fullsack, P., 1995. An arbitrary Lagrangian-Eulerian formulation for creeping flows and its application in tectonic models. Geophys. J. Int. 120, 1–23. doi:10.1111/j.1365-

IP

T

246X.1995.tb05908.x

CR

Gerya, T. V., Meilick, F.I., 2011. Geodynamic regimes of subduction under an active margin: Effects of rheological weakening by fluids and melts. J. Metamorph. Geol. 29, 7–31.

US

doi:10.1111/j.1525-1314.2010.00904.x

AN

Gerya, T. V., Stern, R.J., Baes, M., Sobolev, S. V., Whattam, S.A., 2015. Plate tectonics on the

M

Earth triggered by plume-induced subduction initiation. Nature 527, 221–225.

ED

doi:10.1038/nature15752

Gleason, G.C., Tullis, J., 1995. A flow law for dislocation creep of quartz aggregates determined

PT

with the molten salt cell. Tectonophysics 247, 1–23. doi:10.1016/0040-1951(95)00011-B

CE

Gorczyk, W., Hobbs, B., Gerya, T., 2012. Initiation of Rayleigh–Taylor instabilities in intra-

AC

cratonic settings. Tectonophysics 514–517, 146–155. doi:10.1016/j.tecto.2011.10.016 Gorczyk, W., Hobbs, B., Gessner, K., Gerya, T., 2013. Intracratonic geodynamics. Gondwana Res. 24, 838–848. doi:10.1016/j.gr.2013.01.006 Griffin, W.L., Andi, Z., O’Reilly, S.Y., Ryan, C.G., 1998. Phanerozoic evolution of the lithosphere beneath the Sino-Korean craton, in: Flower, M.F.J., Chung, S., Lo, C., Lee, T. (Eds.), Mantle Dynamics and Plate Interactions in East Asia. pp. 107–126.

57

ACCEPTED MANUSCRIPT Griffin, W.L., Natapov, L.M., O’Reilly, S.Y., van Achterbergh, E., Cherenkova, A.F., Cherenkov, V.G., 2005. The Kharamai kimberlite field, Siberia: modification of the lithospheric mantle by the Siberian Trap event. Lithos 81, 167–187. doi:http://dx.doi.org/10.1016/j.lithos.2004.10.001

IP

T

Griffin, W.L., O’Reilly, S.Y., Abe, N., Aulbach, S., Davies, R.M., Pearson, N.J., Doyle, B.J., Kivi, K.,

CR

2003. The origin and evolution of Archean lithospheric mantle. Precambrian Res. 127, 19– 41. doi:10.1016/S0301-9268(03)00180-3

US

Griffin, W.L., O’Reilly, S.Y., Afonso, J.C., Begg, G.C., 2009. The composition and evolution of

AN

lithospheric mantle: a re-evaluation and its tectonic implications. J. Pet. 50, 1185–1204.

M

doi:10.1093/petrology/egn033

ED

Guillou-Frottier, L., Burov, E., Cloetingh, S., Le Goff, E., Deschamps, Y., Huet, B., Bouchot, V., 2012. Plume-induced dynamic instabilities near cratonic blocks: Implications for P-T-t

PT

paths and metallogeny. Glob. Planet. Change 90–91, 37–50.

CE

doi:10.1016/j.gloplacha.2011.10.007 Hammond, J.O.S., Kendall, J.-M., 2016. Constraints on melt distribution from seismology: a case

AC

study in Ethiopia. Geol. Soc. London, Spec. Publ. 420. doi:10.1144/SP420.14 Hasterok, D., Chapman, D.S., 2011. Heat production and geotherms for the continental lithosphere. Earth Planet. Sci. Lett. 307, 59–70. doi:10.1016/j.epsl.2011.04.034 He, L., 2015. Thermal regime of the North China Craton: Implications for craton destruction. Earth-Science Rev. 140, 14–26. doi:10.1016/j.earscirev.2014.10.011

58

ACCEPTED MANUSCRIPT He, L., 2014. Numerical modeling of convective erosion and peridotite-melt interaction in big mantle wedge: Implications for the destruction of the North China Craton. J. Geophys. Res. 119, 3662–3677. doi:10.1002/2013JB010657 Hirschmann, M.M., 2000. Mantle solidus: Experimental constraints and the effects of peridotite

IP

T

composition. Geochemistry Geophysics Geosystems 1, 1–26. doi:10.1029/2000GC000070

CR

Hirth, G., Kohlstedt, D.L., 1996. Water in the oceanic upper mantle: implications for rheology, melt extraction and the evolution of the lithosphere. Earth Planet. Sci. Lett. 144, 93–108.

US

doi:10.1016/0012-821X(96)00154-9

AN

Hirth, G., Kohlstedt, D.L., 1995. Experimental constraints on the dynamics of the partially

M

molten upper mantle: Deformation in the diffusion creep regime. J. Geophys. Res. 100,

ED

1981–2001. doi:10.1029/94JB02128

Huerta, A.D., Harry, D.L., 2007. The transition from diffuse to focused extension: Modeled

PT

evolution of the West Antarctic Rift system. Earth Planet. Sci. Lett. 255, 133–147.

CE

doi:10.1016/j.epsl.2006.12.011

AC

Huismans, R.S., Beaumont, C., 2007. Roles of lithospheric strain softening and heterogeneity in determining the geometry of rifts and continental margins, in: Karner, G.D., Manatschal, G.D., Pinheiro, L.M. (Eds.), Imgaging, Mapping and Modelling Continental Lithosphere Extension and Breakup. Geological Society, London, Special Publications, pp. 111–138. doi:10.1144/SP282.6 Jamieson, R.A., Beaumont, C., Nguyen, M.H., Culshaw, N.G., 2007. Synconvergent ductile flow

59

ACCEPTED MANUSCRIPT in variable-strength continental crust: Numerical models with application to the western Grenville orogen. Tectonics 26, TC5005. doi:10.1029/2006TC002036 Johnson, S.E., Jin, Z.H., 2009. Magma extraction from the mantle wedge at convergent margins through dikes: A parametric sensitivity analysis. Geochemistry Geophysics Geosystems 10.

IP

T

doi:10.1029/2009GC002419

CR

Jordan, T.H., 1975. The continental Tectosphere. Rev. Geophys. Sp. Phys. 13, 1–12.

US

Karato, S., 2010. Rheology of the deep upper mantle and its implications for the preservation of the continental roots: A review. Tectonophysics 481, 82–98.

AN

doi:10.1016/j.tecto.2009.04.011

M

Karato, S., Wu, P., 1993. Rheology of the upper mantle: A synthesis. Science 260, 771–778.

ED

Karato, S.I., 2011. Water distribution across the mantle transition zone and its implications for

PT

global material circulation. Earth Planet. Sci. Lett. 301, 413–423.

CE

doi:10.1016/j.epsl.2010.11.038 Katsura, T., Shatskiy, A., Manthilake, M.A.G.M., Zhai, S., Fukui, H., Yamazaki, D., Matsuzaki, T.,

AC

Yoneda, A., Ito, E., Kuwata, A., Ueda, A., Nozawa, A., Funakoshi, K. ichi, 2009. Thermal expansion of forsterite at high pressures determined by in situ X-ray diffraction: The adiabatic geotherm in the upper mantle. Phys. Earth Planet. Inter. 174, 86–92. doi:10.1016/j.pepi.2008.08.002 Keranen, K.M., Klemperer, S.L., Julia, J., Lawrence, J.F., Nyblade, A.A., 2009. Low lower crustal velocity across Ethiopia: Is the Main Ethiopian Rift a narrow rift in a hot craton?

60

ACCEPTED MANUSCRIPT Geochemistry Geophysics Geosystems 10. doi:10.1029/2008GC002293 Koornneef, J.M., Davies, G.R., Döpp, S.P., Vukmanovic, Z., Nikogosian, I.K., Mason, P.R.D., 2009. Nature and timing of multiple metasomatic events in the sub-cratonic lithosphere beneath

T

Labait, Tanzania. Lithos 112, 896–912. doi:10.1016/j.lithos.2009.04.039

IP

Koptev, a., Calais, E., Burov, E., Leroy, S., Gerya, T., 2015. Dual continental rift systems

CR

generated by plume–lithosphere interaction. Nat. Geosci. 8, 1–5. doi:10.1038/ngeo2401

US

Le Gall, B., Nonnotte, P., Rolet, J., Benoit, M., Guillou, H., Mousseau-Nonnotte, M., Albaric, J., Deverchère, J., 2008. Rift propagation at craton margin. Distribution of faulting and

AN

volcanism in the North Tanzanian Divergence (East Africa) during Neogene times.

M

Tectonophysics 448, 1–19. doi:10.1016/j.tecto.2007.11.005

ED

Le Roux, V., Bodinier, J.L., Alard, O., O’Reilly, S.Y., Griffin, W.L., 2009. Isotopic decoupling during porous melt flow: A case-study in the Lherz peridotite. Earth Planet. Sci. Lett. 279, 76–85.

PT

doi:10.1016/j.epsl.2008.12.033

CE

Le Roux, V., Bodinier, J.L., Tommasi, A., Alard, O., Dautria, J.M., Vauchez, A., Riches, A.J. V, 2007.

AC

The Lherz spinel lherzolite: Refertilized rather than pristine mantle. Earth Planet. Sci. Lett. 259, 599–612. doi:10.1016/j.epsl.2007.05.026 Lee, C.-T.A., 2003. Compositional variation of density and seismic velocities in natural peridotites at STP conditions: Implications for seismic imaging of compositional heterogeneities in the upper mantle. J. Geophys. Res. 108, 2441. doi:10.1029/2003JB002413

61

ACCEPTED MANUSCRIPT Lee, C.-T.A., Luffi, P., Chin, E.J., 2011. Building and Destroying Continental Mantle. Annu. Rev. Earth Planet. Sci. 39, 59–90. doi:10.1146/annurev-earth-040610-133505 Lee, C.-T.A., Luffi, P., Plank, T., Dalton, H., Leeman, W.P., 2009. Constraints on the depths and temperatures of basaltic magma generation on Earth and other terrestrial planets using

IP

T

new thermobarometers for mafic magmas. Earth Planet. Sci. Lett. 279, 20–33.

CR

doi:10.1016/j.epsl.2008.12.020

Lenardic, A., Moresi, L., Muhlhaus, H., 2000. The role of mobile belts for the longevity of deep

US

cratonic lithosphere: the crumple zone model. Geophys. Res. Lett. 27, 1235–1238.

AN

Lenardic, A., Moresi, L.N., 1999. Some thoughts on the stability of cratonic lithosphere: Effects

M

of buoyancy and viscosity. J. Geophys. Res. 104, 12747–12758. doi:10.1029/1999JB900035

ED

Lenardic, A., Moresi, L.N., Muhlhaus, H., 2003. Longevity and stability of cratonic lithosphere: Insights from numerical simulations of coupled mantle convection and continental

PT

tectonics. J. Geophys. Res. 108, 2303. doi:10.1029/2002JB001859

CE

Lenoir, X., Garrido, C.J., Bodinier, J.L., Dautria, J.M., Gervilla, F., 2001. The recrystallization front

AC

of the Ronda peridotite: Evidence for melting and thermal erosion of subcontinental lithospheric mantle beneath the Alboran basin. J. Pet. 42, 141–158. doi:10.1093/petrology/42.1.141 Levander, A., Schmandt, B., Miller, M.S., Liu, K., Karlstrom, K.E., Crow, R.S., Lee, C.-T.A., Humphreys, E.D., 2011. Continuing Colorado plateau uplift by delamination-style convective lithospheric downwelling. Nature 472, 461–465. doi:10.1038/nature10001

62

ACCEPTED MANUSCRIPT Liu, Z., Park, J., Karato, S., 2016. Seismological detection of low-velocity anomalies surrounding the mantle transition zone in Japan subduction zone. Geophys. Res. Lett. 43, 2480–2487. doi:10.1002/2015GL067097.Received Manatschal, G., Lavier, L., Chenin, P., 2015. The role of inheritance in structuring hyperextended

CR

Gondwana Res. 27, 140–164. doi:10.1016/j.gr.2014.08.006

IP

T

rift systems: Some considerations based on observations and numerical modeling.

McKenzie, D., 1989. Some remarks on the movement of small melt fractions in the mantle.

US

Earth Planet. Sci. Lett. 95, 53–72.

AN

Menzies, M.A., Xu, Y., 1998. Geodynamics of the North China craton, in: Flower, M.F.J., Chung,

M

S., Lo, C., Lee, T. (Eds.), Mantle Dynamics and Plate Interactions in East Asia. Wiley Online

ED

Library, pp. 155–165.

Menzies, M., Xu, Y., Zhang, H., Fan, W., 2007. Integration of geology, geophysics and

PT

geochemistry: A key to understanding the North China Craton. Lithos 96, 1–21.

CE

doi:10.1016/j.lithos.2006.09.008

AC

O’Neill, C.J., Lenardic, a., Griffin, W.L., O’Reilly, S.Y., 2008. Dynamics of cratons in an evolving mantle. Lithos 102, 12–24. doi:10.1016/j.lithos.2007.04.006 O’Reilly, S.Y., Griffin, W.L., 2013. Mantle Metasomatism, in: Harow, D.E., Austrheim, H. (Eds.), Metasomatism and the Chemical Transformation of Rock. Springer, pp. 471–533. doi:10.1007/978-3-642-28394-9 O’Reilly, S.Y., Griffin, W.L., Djomani, Y.P., Morgan, P., 2001. Are lithospheres forever? Tracking

63

ACCEPTED MANUSCRIPT changes in subcontinental lithospheric mantle through time. GSA Today 11, 4–10. Ranalli, G., Piccardo, G.B., Corona-Chávez, P., 2007. Softening of the subcontinental lithospheric mantle by asthenosphere melts and the continental extension/oceanic spreading

T

transition. J. Geodyn. 43, 450–464. doi:10.1016/j.jog.2006.10.005

IP

Richard, G.C., Bercovici, D., 2009. Water-induced convection in the Earth’s mantle transition

CR

zone. J. Geophys. Res. 114, 1–11. doi:10.1029/2008JB005734

US

Richard, G.C., Bercovici, D., Karato, S.-I., 2006. Slab dehydration in the Earth’s mantle transition

AN

zone. Earth Planet. Sci. Lett. 251, 156–167. doi:10.1016/j.epsl.2006.09.006 Rudnick, R.L., McDonough, W.F., O’Connell, R.J., 1998. Thermal structure, thickness and

M

composition of continental lithosphere. Chem. Geol. 145, 395–411. doi:10.1016/S0009-

ED

2541(97)00151-4

PT

Sand, K.K., Waight, T.E., Pearson, D.G., Nielsen, T.F.D., Makovicky, E., Hutchison, M.T., 2009. The lithospheric mantle below southern West Greenland: A geothermobarometric

CE

approach to diamond potential and mantle stratigraphy. Lithos 112, 1155–1166.

AC

doi:10.1016/j.lithos.2009.05.012 Schutt, D.L., Lesher, C.E., 2010. Compositional trends among Kaapvaal Craton garnet peridotite xenoliths and their effects on seismic velocity and density. Earth Planet. Sci. Lett. 300, 367– 373. doi:10.1016/j.epsl.2010.10.018 Schutt, D.L., Lesher, C.E., 2006. Effects of melt depletion on the density and seismic velocity of garnet and spinel lherzolite. J. Geophys. Res. 111, 1–24. doi:10.1029/2003JB002950

64

ACCEPTED MANUSCRIPT Shapiro, S.S., Hager, B.H., Jordan, T.H., 1999. Stability and dynamics of the continental tectosphere. Lithos 48, 115–133. Sheng, J., Liao, J., Gerya, T., 2016. Numerical modeling of deep oceanic slab dehydration: Implications for the possible origin of far field intra-continental volcanoes in northeastern

IP

T

China. J. Asian Earth Sci. 117, 328–336. doi:10.1016/j.jseaes.2015.12.022

CR

Sizova, E., Gerya, T., Brown, M., Perchuk, L.L., 2010. Subduction styles in the Precambrian:

US

Insight from numerical experiments. Lithos 116, 209–229. doi:10.1016/j.lithos.2009.05.028 Sleep, N.H., 2007. Edge-modulated stagnant-lid convection and volcanic passive margins.

AN

Geochemistry Geophysics Geosystems 8. doi:10.1029/2007GC001672

ED

doi:10.1029/2001JB000169

M

Sleep, N.H., 2003. Survival of Archean cratonal lithosphere. J. Geophys. Res. 108, 1–29.

PT

Sobolev, S. V, Sobolev, A. V, Kuzmin, D. V, Krivolutskaya, N. a, Petrunin, A.G., Arndt, N.T., Radko, V. a, Vasiliev, Y.R., 2011. Linking mantle plumes, large igneous provinces and

CE

environmental catastrophes. Nature 477, 312–316. doi:10.1038/nature10385

AC

Soustelle, V., Tommasi, A., Bodinier, J.L., Garrido, C.J., Vauchez, A., 2009. Deformation and reactive melt transport in the mantle lithosphere above a large-scale partial melting domain: The Ronda peridotite Massif, Southern Spain. J. Pet. 50, 1235–1266. doi:10.1093/petrology/egp032 St-Onge, M.R., Van Gool, J.A.M., Garde, A.A., Scott, D.J., 2009. Correlation of Archaean and Palaeoproterozoic units between northeastern Canada and western Greenland:

65

ACCEPTED MANUSCRIPT constraining the pre-collisional upper plate accretionary history of the Trans-Hudson orogen, in: Cawood, P.A., Kroner, A. (Eds.), Earth Accretionary Systems in Space and Time. Geological Society, London, Special Publications, pp. 193–235. doi:10.1144/SP318.7 Tang, Y.-J., Zhang, H.-F., Ying, J.-F., Su, B.-X., 2013a. Widespread refertilization of cratonic and

IP

T

circum-cratonic lithospheric mantle. Earth-Science Rev. 118, 45–68.

CR

doi:10.1016/j.earscirev.2013.01.004

Tang, Y.-J., Zhang, H.F., Santosh, M., Ying, J.F., 2013b. Differential destruction of the North

US

China Craton: A tectonic perspective. J. Asian Earth Sci. 78, 71–82.

AN

doi:10.1016/j.jseaes.2012.11.047

M

Tappe, S., Foley, S.F., Stracke, A., Romer, R.L., Kjarsgaard, B.A., Heaman, L.M., Joyce, N., 2007.

ED

Craton reactivation on the Labrador Sea margins: 40Ar/39Ar age and Sr–Nd–Hf–Pb isotope constraints from alkaline and carbonatite intrusives. Earth Planet. Sci. Lett. 256, 433–454.

PT

doi:10.1016/j.epsl.2007.01.036

CE

Tappe, S., Pearson, D.G., Nowell, G., Nielsen, T., Milstead, P., Muehlenbachs, K., 2011. A fresh isotopic look at Greenland kimberlites: Cratonic mantle lithosphere imprint on deep source

AC

signal. Earth Planet. Sci. Lett. 305, 235–248. doi:10.1016/j.epsl.2011.03.005 Tian, Y., Zhao, D., Sun, R., Teng, J., 2009. Seismic imaging of the crust and upper mantle beneath the North China Craton. Phys. Earth Planet. Inter. 172, 169–182. doi:10.1016/j.pepi.2008.09.002 Tommasi, A., Vauchez, A., 2001. Continental rifting parallel to ancient collisional belts: an effect

66

ACCEPTED MANUSCRIPT of the mechanical anisotropy of the lithospheric mantle. Earth Planet. Sci. Lett. 185, 199– 210. Vauchez, A., Tommasi, A., Mainprice, D., 2012. Faults (shear zones) in the Earth’s mantle.

T

Tectonophysics 558–559, 1–27. doi:10.1016/j.tecto.2012.06.006

IP

Vernières, J., Godard, M., Bodinier, J.-L., 1997. A plate model for the simulation of trace

CR

element fractionation during partial melting and magma transport in the Earth’s upper

US

mantle. J. Geophys. Res. 102, 24771–24784. doi:10.1029/97JB01946 Wang, H., van Hunen, J., Pearson, D.G., 2015. The thinning of subcontinental lithosphere: The

AN

roles of plume impact and metasomatic weakening. Geochemistry Geophysics Geosystems

M

16, 1156–1171. doi:10.1002/2015GC005784

ED

Wang, H., van Hunen, J., Pearson, D.G., Allen, M.B., 2014. Craton stability and longevity: The roles of composition-dependent rheology and buoyancy. Earth Planet. Sci. Lett. 391, 224–

PT

233. doi:10.1016/j.epsl.2014.01.038

CE

Wang, Z., Kusky, T.M., Capitanio, F.A., 2016. Lithosphere thinning induced by slab penetration

AC

into a hydrous mantle transition zone. Geophys. Res. Lett. 567–577. doi:10.1002/2016GL071186 Wei, W., Xu, J., Zhao, D., Shi, Y., 2012. East Asia mantle tomography: New insight into plate subduction and intraplate volcanism. J. Asian Earth Sci. 60, 88–103. doi:10.1016/j.jseaes.2012.08.001 Wenker, S., Beaumont, C., 2016. Effects of lateral strength contrasts and inherited

67

ACCEPTED MANUSCRIPT heterogeneities on necking and rifting of continents. Tectonophysics. doi:http://dx.doi.org/10.1016/j.tecto.2016.10.011 Windley, B.F., Maruyama, S., Xiao, W.J., 2010. Delamination/thinning of sub-continental lithospheric mantle under eastern China: The role of water and multiple subduction. Am. J.

IP

T

Sci. 310, 1250–1293. doi:10.2475/10.2010.03

CR

Wu, F.Y., Walker, R.J., Yang, Y.H., Yuan, H.L., Yang, J.H., 2006. The chemical-temporal evolution

5013–5034. doi:10.1016/j.gca.2006.07.014

US

of lithospheric mantle underlying the North China Craton. Geochim. Cosmochim. Acta 70,

AN

Xia, Q.K., Liu, J., Liu, S.C., Kovács, I., Feng, M., Dang, L., 2013. High water content in Mesozoic

M

primitive basalts of the North China Craton and implications on the destruction of cratonic

ED

mantle lithosphere. Earth Planet. Sci. Lett. 361, 85–97. doi:10.1016/j.epsl.2012.11.024 Xu, Y., Li, H., Pang, C., He, B., 2009. On the timing and duration of the destruction of the North

PT

China Craton. Chinese Sci. Bull. 54, 3379–3396. doi:10.1007/s11434-009-0346-5

CE

Yamasaki, T., O’Reilly, B., Readman, P., 2006. A rheological weak zone intensified by post-rift

AC

thermal relaxation as a possible origin of simple shear deformation associated with reactivation of rifting. Earth Planet. Sci. Lett. 248, 134–146. doi:10.1016/j.epsl.2006.05.019 Yoshida, M., 2012. Dynamic role of the rheological contrast between cratonic and oceanic lithospheres in the longevity of cratonic lithosphere: A three-dimensional numerical study. Tectonophysics 532–535, 156–166. doi:10.1016/j.tecto.2012.01.029 Yoshida, M., 2010. Preliminary three-dimensional model of mantle convection with deformable,

68

ACCEPTED MANUSCRIPT mobile continental lithosphere. Earth Planet. Sci. Lett. 295, 205–218. doi:10.1016/j.epsl.2010.04.001 Zhao, D., Tian, Y., Lei, J., Liu, L., Zheng, S., 2009. Seismic image and origin of the Changbai intraplate volcano in East Asia: Role of big mantle wedge above the stagnant Pacific slab.

IP

T

Phys. Earth Planet. Inter. 173, 197–206. doi:10.1016/j.pepi.2008.11.009

CR

Zheng, J.P., Lee, C.-T.A., Lu, J.G., Zhao, J.H., Wu, Y.B., Xia, B., Li, X.Y., Zhang, J.F., Liu, Y.S., 2015. Refertilization-driven destabilization of subcontinental mantle and the importance of initial

US

lithospheric thickness for the fate of continents. Earth Planet. Sci. Lett. 409, 225–231.

AN

doi:10.1016/j.epsl.2014.10.042

M

Zhu, R.X., Chen, L., Wu, F.Y., Liu, J.L., 2011. Timing, scale and mechanism of the destruction of

AC

CE

PT

ED

the North China Craton. Sci. China Earth Sci. 54, 789–797. doi:10.1007/s11430-011-4203-4

69

ACCEPTED MANUSCRIPT Highlights Cratons are cold and strong, yet they occasionally rift, e.g. North China craton



Does mantle metasomatism weaken cratons sufficiently that they will rift?



Numerical models with progressive craton melt metasomatism are rift tested



Combined metasomatic heating, refertilization and densification lead to rifting



Once weakened sufficiently cratons are not protected by surrounding standard lithosphere

AC

CE

PT

ED

M

AN

US

CR

IP

T



70