Roles of Continental Mid‐Lithosphere Discontinuity in the Craton Instability Under Variable Tectonic Regimes

The continental mid‐lithosphere discontinuity (MLD) is widely detected within cratons, with the dominant depth range of 70–100 km and a significant reduction of shear‐wave velocity of 2%–12%. However, the formation mechanism and corresponding strength of the MLD are widely debated, which may strongly affect the roles of MLD in craton evolution. The comparisons among variable mechanisms indicate that the strength of the MLD varies from the relatively high viscosity of wet olivine to the rather low viscosity of antigorite. Thus, systematic numerical modeling has been conducted with the MLD of contrasting strengths, that is, the wet olivine‐induced MLD or antigorite‐induced MLD, to investigate the roles of MLD in the craton instability under variable tectonic regimes (stable, extension, compression, mantle flow traction, or mantle plume). The models show that the cratonic lithosphere with wet olivine‐induced MLD maintains its stability under all the tectonic regimes. In contrast, the antigorite‐induced MLD with lowest viscosity could significantly promote the decoupling of lithosphere, and facilitate the lithospheric deformation. However, lithospheric delamination only occurs with the rather weak MLD interacting with the sub‐plate asthenosphere upwelling during craton extension or mantle plume activity. The sufficient amount of melts is essential for this process, which requires a large amount of extension or a mantle plume with rather high temperature anomaly and large size. Therefore, craton destruction is still difficult, and requires additional strict conditions. This may explain the general stability of most cratons with widespread MLDs.

Plain Language Summary Geophysical observations have revealed a seismic discontinuity within cratons, which is defined as continental mid-lithosphere discontinuity (MLD).The MLD is often considered to be rather weak, which could facilitate the lithospheric delamination, and even the craton destruction.However, most cratons preserve the stability during long-term evolution.Thus, the viscosity of the MLD and its roles in craton evolution require further investigations.In this study, we have conducted systematic numerical modeling to investigate the effects of the MLD on craton instability under variable tectonic regimes.The wet olivine-induced MLD and antigorite-induced MLD are two end-members of MLD with contrasting strengths, which are compared in the present models.The model results indicate that only the connection and interaction between the weak antigorite-induced MLD and the sub-plate asthenosphere could cause craton destruction during tectonic extension or upwelling of mantle plume.The sufficient amount of melts is essential for this process to destroy the lower part of lithosphere.Thus, craton destruction requires additional strict conditions, which may explain the stability of most cratons on the Earth, although with the presence of MLD.water capacity (Figure S1 in Supporting Information S1) could absorb sufficient amount of water from dehydration process in the deep mantle, such as subduction process (Li, 2020(Li, , 2022)).Further on, a high-water-content layer, accompanied by a significant Vs reduction (Figure S2 in Supporting Information S1), may form gradually within cratonic lithosphere during the long-term evolution, which would be recognized as the MLD revealed by geophysical observations.Based on the "water collector" hypothesis, the viscosity of MLD would decrease due to its high water content.The upper mantle is mainly composed of olivine (>60%) and pyroxene, and the strength of olivine is generally lower than that of pyroxene (Bystricky & Mackwell, 2001).Thus, the rheology of olivine largely controls the deformation of the upper mantle (Bürgmann & Dresen, 2008;Wang, 2016).The viscosity of olivine decreases with the increase of water content (Hirth & Kohlstedt, 2003).Thus, for the wet olivine-induced MLD, the viscosity is relatively lower than that of the surrounding mantle.However, based on the comparison between the stability field of serpentinite (Schmidt & Poli, 1998) and the P-T conditions of MLD, serpentinization may occur at the shallow part of the MLD, which has been proposed as a high-water-content layer (Figure S3 in Supporting Information S1).At the dominate depth of the MLD, the viscosity of serpentinite may be strongly controlled by antigorite, that is, the only stable serpentine mineral under the P-T conditions of up to 1.5 GPa and 650°C (Hilairet et al., 2006;Ulmer & Trommsdorff, 1995).In addition, the viscosity of antigorite is much lower than that of wet olivine (Hilairet et al., 2007).Thus, the antigorite-induced MLD may represent the case with fully serpentinized MLD and lowest viscosity.We proposed that the viscosity of MLD varies dramatically depending on the degree of serpentinization, with the wet olivine-induced MLD and antigorite-induced MLD as the two end-member cases of MLD.
In the previous models, the MLD is usually considered to be rather weak and thus may strongly influence the craton instability.The interaction between a rather weak MLD and subduction or mantle plume could promote the delamination and even destruction of cratonic lithosphere (Liu et al., 2018;Shi et al., 2020Shi et al., , 2021;;Wang et al., 2018).In addition, a weak MLD could decouple the lithosphere, and foster the lateral displacement between the overlying and underlying lithosphere during the plate motion (Wang et al., 2017).However, the viscosity range of MLD is not well constrained, and systematic comparisons of the roles of MLD in craton evolution under variable tectonic regimes are still lacked.On the other hand, the MLD exists within most cratons, whereas only very few cratons have undergone destruction (Figure 1a).Thus, the effects of MLD with contrasting strengths on the craton instability require further evaluations.
In this study, we aim to investigate the effects of MLD on the craton instability under variable tectonic regimes (stable craton, tectonic extension, tectonic compression, mantle flow traction, as well as the upwelling of mantle plume), using a thermomechanical model with fluid migration.Two end-member cases of the MLD with contrasting viscosities, that is, wet olivine-induced MLD and antigorite-induced MLD, are applied.In addition, the rheological strength of MLD induced by variable formation mechanisms is systematically compared to constrain the viscosity range of the MLD.

Numerical Model Setup
Numerical models in a 2-D Cartesian box of 1,200 km × 800 km are conducted with the program I2VIS (Gerya, 2010).The general algorithms and implementations are shown in Section S1 in Supporting Information S1 and the detailed model setup and boundary conditions shown in Section S2 of the Supporting Information S1.

Mantle Viscous Flow Law With Fluid/Melt Effects
The previous numerical models generally used dry or wet olivine rheology separately, for the normal and hydrated mantle rocks, respectively.In addition, the effective viscosity of hydrated rock is not further dependent on the water content.In this study, the rheology of mantle rock depends on water content and melt fraction according to Hirth and Kohlstedt (2003): (2) 10.1029/2023JB028022

of 27
Where f H2O is water content, d grain size (d = 10 mm), ϕ melt fraction.A H (pre-exponential factor), n (creep exponent), p (grain size exponent), r (water content exponent), α (pre-melt-fraction factor), E (activation energy) and V (activation volume) are rheological parameters (Hirth & Kohlstedt, 2003).In the present model, when the water content of mantle rock is lower than 3 ppm, the rheological parameters of dry olivine are used (C* and E* in Table 1); otherwise, the parameters of wet olivine are used (D* and F* in Table 1).This algorithm is applied in order to keep the effective viscosity of hydrated mantle rocks lower than that of dry ones, because the calculated viscosity of wet mantle rock with rather low water content of <3 ppm (Equation 1) is even higher than the dry mantle (Figure S4 in Supporting Information S1) which is thus not reasonable.

Effective Viscosity of MLD
Similar to the normal mantle rocks, the effective viscosity of the MLD is integrated with the ductile viscosity (η ductile ), the plastic equivalent (η plastic ), and the Peierls viscosity (η peierls ).However, variable viscous flow laws and rheological parameters are applied for the MLD resulting in different rheological strengths.For the wet olivine-induced MLD, the ductile viscosity is calculated by the olivine flow law of Hirth and Kohlstedt (2003) (Equations 1 and 2) with the rheological parameters of wet olivine (D* and F*) in Table 1.This ductile viscosity of the MLD is not low enough, and even exceeds the cut-off value of 10 25 Pa•s at the shallow part of the lithosphere (blue lines in Figures 2e and 2f).In contrast, for the antigorite-induced MLD, the ductile viscosity is calculated by the flow law of Ranalli (1995) (Equation S10 in Supporting Information S1) with the rheological parameters of antigorite (G1*, G2* or G3*) in Table 1 (Agard et al., 2016;Hilairet et al., 2007).The viscosity of the antigorite-induced MLD is much lower than that of the wet olivine one at the dominant depth of the MLD.Finally, the flow law parameters (G3*) are applied for the antigorite-induced MLD in the numerical models, since its viscosity is the lowest (red lines in Figures 2e and 2f).The effects of contrasting rheological strengths of the MLD are systematically compared in the present study.

Aqueous Fluid Generation and Migration
The water could be carried into the deep Earth by hydrous minerals, which could be released during the dehydration process (Figures 3a and 3b) (Ohtani, 2020).Further on, the excess water would migrate upward and hydrate the surrounding mantle.In the numerical model, the aqueous fluid markers are generated, which store a certain amount of water, when the dehydration process occurs (Figures 3c and 3d).These fluid markers migrate with a calculated velocity, and hydrate the water-unsaturated mantle, whereas its water is gradually consumed until the water content becomes zero (Figure 3f).
Besides the subduction, upwelling of hydrous mantle plume could also carry water from the deep mantle into the continental lithosphere (Kuritani et al., 2019).The complicated and regional (de)hydration processes during craton evolution cannot be well constrained.Thus, a new generation routine of hydrated blocks as shown in Figure 3 is applied to simplify the simulation of the whole process of (de)hydration in the deep mantle (Fu ) is ∼15 km 2 /Myr.
The detailed generation routine of hydrated blocks is shown in Fu et al. (2022) and the Supporting Information S1 (Section S1).
In addition, the water capacity of nominally anhydrous minerals (NAMs) and the actual water content of mantle transition zone (MTZ) are still controversial (Li, 2020;Peslier et al., 2017).Thus, hydrated blocks are generated above the MTZ in the present model, and a constant water capacity of NAMs of 50 ppm is applied for the whole upper mantle.The rate of hydrated block generation (R h ) is a critical parameter that cannot be well constrained due to the complex tectonic events associated with dehydration beneath the craton.Meanwhile, the time for the formation of MLD is still controversial, which may be related to the time of global cratonization and the onset of plate tectonics.Variable generation rates and formation time have been discussed in Fu et al. (2022).Therefore, in this study, the constant generation rate of R h = 15 km 2 /Myr and model time of 550 Myr are applied in the numerical model to simulate the formation of the MLD.

Reference Model With Wet Olivine-Induced MLD
In this model, the MLD is formed gradually with wet olivine rheology (Hirth & Kohlstedt, 2003) (Equations 1 and 2) with the rheological parameters shown in Table 1 (D* and F*).
In the numerical model, hydrated blocks are generated at the depth of 410 km.These blocks hydrate the surrounding mantle until the water content of mantle reaches its capacity, and then move upward, which lead to the reduction of viscosity along the migration pathway (Figures 4a and 4aʹ).During the model evolution, the hydrated blocks gradually move into the lithosphere, and lead to the formation of a high-water-content layer within cratonic lithosphere (Figures 4b-4d).Finally, at the model time of 550 Myr, a continuous layer with the maximum water content of ∼3% is formed at the depth range of 70-100 km, which further leads to a significant Vs drop with the maximum of ∼9% (Figure 4e, Figure S5 in Supporting Information S1).The detailed model evolution and calculation of Vs are discussed in Section S2 of the Supporting Information S1.
We propose that the observed MLD may represent this layer with high water content and significant Vs drop, which is formed gradually during craton evolution.The water content of the MLD increases gradually, whereas the Vs at the depth of the MLD decreases significantly, due to the (de)hydration process in the deep mantle.However, for the wet olivine-induced MLD, the viscosity of the MLD is not low enough, which is only lower than that of the surrounding mantle by about half an order of magnitude at ∼100 km (Figure 4eʹ, Figure S6 in Supporting Information S1).
With the formation of the MLD, the thermal lithosphere cannot maintain a cooling-induced steady thickening.
After the model time of 250 Myr, convective dripping occurs at the bottom of the lithosphere, accompanied by the significant reduction of viscosity at these locations (Figures 4cʹ and 4dʹ).The thickness of thermal lithosphere decreases rapidly at 250 Myr, and then remains stable at around 200 km (i.e., identical to the initial model), whereas the thickness of compositional lithosphere reduces slightly and continuously after 250 Myr (Figure S6 in Supporting Information S1).
In the numerical model, the upper boundary of MLD may still move upward, and even exceed the dominant depth of natural MLD of 70-100 km, during the continuous generation of hydrated blocks (Figure S7 in Supporting Information S1).This is partly due to the limited water capacity (3 wt%) of mantle rock applied in this study, (e-f) The viscosity profiles of olivine and antigorite, calculated using the rheological parameters in Table 1, with the strain rate (   ) of 10 −19 s −1 , 10 −14 s −1 , respectively.The gray rectangles represent the dominant MLD depth range of 70-100 km as shown in Figure 1.
whereas the maximum water capacity of serpentinite can be as high as ∼13% (Hilairet et al., 2006).On the other hand, the occurrence of grain size reduction due to brittle/ductile transition may form a barrier zone, which could further contribute to the accumulation of fluid and melt within the MLD of continental lithosphere (Liu et al., 2022).Collectively, the MLD may form and absorb a large amount of water in the middle of lithosphere, instead of migrating upward continuously.In the previous numerical model (e.g., Li et al., 2019), when the dehydration process occurs (c), the aqueous fluid markers are generated in the grid (d), then migrate upward and hydrate the surrounding rocks (f).In the present model, a new generation routine of hydrated blocks (e) is applied to simplify the simulation of whole process of (de)hydration in the deep mantle.The gray dots represent mantle rocks, while the blue dots represent the addition of released water.The yellow and green grids represent the dry and hydrated mantle domains, respectively.

Reference Model With Antigorite-Induced MLD
In this model, the MLD is formed gradually with a low viscosity, calculated by the flow law of Ranalli (1995) (Equation S10 in Supporting Information S1) with the rheological parameters of weakest antigorite (G3*).
For the composition field, the evolution of the model is similar to the model with the wet olivine-induced MLD.The MLD forms gradually with the random and continuous generation of hydrated blocks.Meanwhile, for the thermal lithosphere, cooling-induced thickening does not last for 550 Myr, that is, the final stage of model evolution.At the model time of 250 Myr, the convective dripping develops at the base of the lithosphere, due to the Rayleigh-Taylor instability (Figures 5a-5e).The strength of the antigorite-induced MLD is much weaker than that of wet olivine-induced MLD (cf.Figures 5 and 4).The viscosity of the MLD decreases significantly during model evolution, and even reduces to ∼10 19 Pa•s at 550 Myr (Figures 5aʹ-5eʹ, Figure S8 in Supporting Information S1).In addition, a weak MLD could lead to strain localization, and further contribute to the reduction of the strain rate dependent viscosity.However, the lithospheric delamination does not occur along the MLD, even with such a low viscosity.

Effect of Contrasting MLD Strengths on Craton Instability
In the above reference models, all the other parameters are identical, except the flow law parameters of the MLD.At the model time of 550 Myr, the MLDs are located at the same depth range of 70-100 km, with the same water content of ∼3% (Figures 6a and 6e).The MLDs are accompanied by the same Vs and Vs drop, which is interpolated from the Vs database based on the temperature, pressure and water content (Figures 6c and 6d).The viscosity of the MLDs shows a large contrast.The viscosity of antigorite-induced MLD is much lower than that of wet olivine-induced MLD, with a viscosity contrast of about five orders at 550 Myr (Figure 6b).However, the evolution of the lithosphere exhibits the same pattern in both models.The thermal lithosphere keeps a cooling-induced steady thickening before 250 Myr, and experiences lithospheric thinning rapidly due to the occurrence of convective dripping.Afterward, the craton maintains a dynamically balanced lithospheric thickness of about 200 km (Figure 6f).Therefore, for a craton without boundary forces, the presence of the MLD, even the weakest MLD, has little effect on the craton instability, and could not lead to lithospheric delamination or craton destruction (Figure 5).
The long-term evolution of cratons involves complicated tectonic processes.Thus, the roles of the MLD in craton evolution are further investigated under variable tectonic regimes, including tectonic extension, tectonic compression, mantle flow traction and the upwelling of mantle plume.

MLD Within Craton Under Tectonic Extension
In the tectonic extension regime, constant velocities with the same value but opposite directions are prescribed at the entire left and right boundaries (Liao et al., 2013;Liao & Gerya, 2014).The velocity at the lower boundary is calculated based on the mass conservation.The velocities are applied from the model time of 500 Myr, when the high-water-content layer comparable to the natural MLD is formed within craton (Figure S9 in Supporting Information S1, two reference model results at 500 Myr without boundary forces).The boundary velocities will be canceled after the model time of 520 Myrs.Both representative models in Figure 7 have the same extension rate of 0.5 cm/yr (half rate).There is no weak seed to localize the strain in the model domain.Thus, the position of continental break-up is relatively random.
In the model with wet olivine-induced MLD (Figure 7, left column), diffused deformation dominates at the early extension stage, until strain localization occurs (Figures 7a and 7aʹ).Subsequently, the lithosphere above the MLD has been broken during extension process.Meanwhile, the MLD migrates upward to shallower depth.At the model time of 510 Myr, decompression melting occurs at the base of the lithosphere.The melts with low viscosity (Equation 1) and low density form gradually during continuous extension (Figures 7b-7c, 7bʹ-7cʹ).However, lithospheric delamination does not occur during continental break-up (Figures 7d and 7dʹ).The average depth of compositional LAB decreases steadily during extension process until the final stage of 520 Myr.The evolution of continental lithosphere with variable extension rates exhibits the same pattern, in which the thickness of lithosphere does not drop rapidly (Figure 7e).The wet olivine-induced MLD is not weak enough and thus could not lead to lithospheric delamination during craton extension, even with the longer time of extension (e.g., Figure S10 in Supporting Information S1 with the extension time of 50 Myr).In the model with antigorite-induced MLD (Figure 7, right column), the model evolution in early extension stage is similar to the previous model with wet olivine-induced MLD (Figures 7f and 7fʹ).A MLD with rather low viscosity and relatively low density, prefers to move upward during craton extension, which further hinders the deformation of the underlying lithosphere (Figures 7g and 7gʹ).At the model time of 517 Myr, the underlying lithosphere has been destructed due to the formation and accumulation of melts at the location of continental break-up.The melts with low viscosity (Equation 1) migrate along the weak MLD (Figure 7h), and further lead to the lithospheric delamination, which corresponds to a significant decrease of average thickness of lithosphere (red line in Figure 7j).Although the melt fraction during craton extension is similar to that with wet olivine-induced MLD (Figure S11 in Supporting Information S1), the model results are quite different.Such a weakest MLD plays an important role in the decoupling of overlying and underlying lithosphere (Figure S12 in Supporting Information S1), which facilitates the delamination and even destruction of craton under tectonic extension.
Model results with variable extension rates indicate that the two parameters (i.e., extension rate and extension time) could compensate each other, for example, craton destruction could also occur by increasing the time of extension, if the extension rate is lower (Figure 7j).Thus, extension rate multiplied by extension time, that is, the amount of extension, has a significant effect on the craton stability.Meanwhile, the amount of extension strongly controls the melt fraction of decompression melting.The sufficient amount of melts is essential for craton destruction during craton extension, which requires large amount of extension.
The cratonic lithospheric mantle may be depleted with lower reference density than the asthenospheric mantle at the same pressure and temperature; the density reduction could be ∼2.5%-1%(Lee et al., 2011;Zhu et al., 2021).Thus, an additional set of models with depleted lithospheric mantle (ρ 0 = 3,250 kg/m 3 ) is conducted (Figure S13 in Supporting Information S1).In these models, a weak seed is set in the lithospheric mantle in order to localize the position for initial break-up, and all the other parameters are identical as the representative models in Figure 7.In these additional models, the depleted lithospheric mantle with relatively low density significantly prevents the lithospheric delamination along MLD, and further improves the stability of craton (Figure S14 in Supporting Information S1).

MLD Within Craton Under Tectonic Compression
In this regime, the configuration of velocity boundary is similar to the model of craton extension in Section 3.2, but with compression.In the representative models, convergence rates are applied in the opposite directions with the same value of 0.5 cm/yr at the entire left and right boundaries from the model time of 500 Myr, which will be canceled after the model time of 520 Myrs (Figure 8).
In the model with wet olivine-induced MLD (Figure 8, left column), the lithosphere thickens extensively with diffused deformation at the early stage of craton compression.At the model time of 510 Myr, strain localization occurs, and the left part of the lithosphere underthrusts beneath the right part.The MLD with relatively high viscosity could not decouple the lithosphere.The lithosphere as a whole localizes the strain along the pathway of underthrusting (Figure S15 in Supporting Information S1), which further leads to the reduction of viscosity along this pathway (Figures 8bʹ-8cʹ).However, the lithosphere does not undergo delamination or destruction along the MLD due to the relatively high viscosity.
In the model with antigorite-induced MLD (Figure 8, right column), the underthrusting of the lithosphere also occurs during model evolution.It is worth noting that there is no weak seed prescribed in the initial model; thus, the direction of underthrusting is opposite to that of the above model, due to the uncertainty of strain localization (Figures 8f and 8g).In the present model, the weak MLD localizes large strain, and decouples the overlying and underlying lithosphere.The strain localization mainly occurs at the MLD and the overlying lithosphere during compression process, which further hinders the deformation of the underlying lithosphere (Figure S15 in Supporting Information S1).Thus, the strength of underlying lithosphere is only slightly affected by the underthrusting.The evolution of viscosity field is obviously different from the model with wet olivine-induced MLD.In addition, the underlying lithosphere with relatively high viscosity further inhibits the destruction of craton, although a rather weak MLD is formed within the cratonic lithosphere.
Model results with variable convergence rates indicate the averaged depth of compositional LAB increases monotonically during compression process, in the model with either wet olivine-induced MLD (Figure 8d) or antigorite-induced MLD (Figure 8h).There is no significant reduction in the average thickness of lithosphere.The occurrence of MLD with variable strength could not lead to lithospheric delamination during craton compression, although the MLD would facilitate the strain localization and promote the deformation of lithosphere.

MLD Within Craton Under Mantle Flow Traction
In this section, models with the sub-lithosphere mantle flow traction are further conducted.The mantle flow velocities are fixed with the same direction below the depth of 200 km on the left and right boundaries.Permeable boundary condition is applied for the left boundary during mantle flow, in which a ghost region is located to the left with additional rock markers.Once markers migrate out of the model domain from the right side boundary, additional markers could migrate across the left side permeable boundary into the model domain, in order to guarantee the mass conservation.In the representative models, the mantle flow velocity of 5 cm/yr is applied from the model time of 500 Myr, and will be canceled after the model time of 520 Myrs (Figure 9).
In the model with wet olivine-induced MLD (Figure 9, left column), the hydrated sub-lithospheric mantle moves horizontally from left to right, and leaves the model domain during mantle flow.Meanwhile, the dry sub-lithospheric mantle enters the model domain through the left boundary, and is hydrated with the generation and migration of hydrated blocks (Figures 9a-9c).The water content of these blocks is gradually consumed in sub-lithospheric mantle.Thus, the hydrated blocks are more difficult to migrate into the continental lithosphere during mantle flow.The water content of the MLD would increase more slowly than the reference model (Figure 4).On the other hand, the model evolution with antigorite-induced MLD is similar to that with wet olivine-induced MLD (Figures 9aʹ-9cʹ), although the viscosity differs significantly (Figures 9d and 9dʹ).Craton could remain long-term stability during mantle flow, and the MLDs with contrasting strengths have rare effect on craton evolution.
Models with variable velocities of mantle flow show that the lithosphere has similar viscosity structure and almost constant thickness, with the similar strain at the final model time of 520 Myr (Figure S16 in Supporting Information S1).The occurrence of mantle flow has rather limited effect on the cratonic lithosphere, because the LAB with low viscosity decouples the lithosphere and asthenosphere.In addition, the strain rate dependent viscosity at LAB decreases with the increasing mantle flow velocity, which further weakens the effect of mantle flow traction.

Interaction Between MLD and Mantle Plume
Mantle plume may play an important role in the lithospheric thinning and even craton destruction (Shi et al., 2021;Wu et al., 2014).In this section, the interactions between the MLD and mantle plume are studied.The mantle plume in a circular shape is configured at the model time of 500 Myr.The center of the plume is located at the depth of 410 km.In the representative models, the plume has a diameter of 200 km and the temperature of 2,123 K which is 400 K higher than that of surrounding mantle at 410 km (∼1723 K).The corresponding temperature field evolution is shown in Figure S17 in Supporting Information S1.
In the model with wet olivine-induced MLD (Figure 10, left column), the melts are generated during upwelling of the mantle plume due to the high temperature anomaly and rapid decompression.At the model time of 500.02Myr, the mantle plume impinges the base of the lithosphere, and extends laterally along the LAB.The lithospheric thinning occurs due to the thermal erosion of mantle plume.However, the instability and dripping only localize at the base of the lithosphere, but the cratonic keel preserves its general stability (Figures 10b-10d).In contrast, in the model with antigorite-induced MLD (Figure 10, right column), the weak MLD and underlying lithospheric mantle experience larger strain during the upwelling of plume.The lower part of lithosphere with low viscosity could be destructed by the mantle plume easily.Further on, the plume connects with the weak MLD, and leads to the delamination of lithosphere along the MLD (Figures 10f-10h).Thus, the interaction between the antigorite-induced MLD and the mantle plume could cause craton destruction during model evolution.
An additional set of models with depleted lithospheric mantle (ρ 0 = 3,250 kg/m 3 ) is further conducted (Figure 11); all the other parameters are identical to the model with antigorite-induced MLD in Figure 10 (right column).During the model evolution, the base of lithospheric mantle is disturbed by the upwelling mantle plume (Figures 11a and 11aʹ).However, the plume cannot destroy the lower part of lithosphere and contact with the MLD (Figures 11b and 11bʹ).The cratonic lithosphere with relatively low density remains stable, even with the weakest antigorite-induced MLD.It indicates that the depleted lithospheric mantle with relatively lower density could further improve the stability of craton.
Systematic numerical models are conducted with variable MLD rheology (wet olivine vs. antigorite), plume properties (diameter and temperature), as well as different degrees of lithospheric mantle depletion (fertile, depleted, or mixed).The model results, as summarized in Figure 12, indicate that the rheological strength of MLD plays a first-order role in craton destruction.In all the models with wet olivine-induced MLD, the lithospheric delamination is not predicted, even in the cases with an entire fertile lithospheric mantle (Figure 12a) and a large and hot plume.In the alternative regime with the weakest antigorite-induced MLD, the craton destruction occurs more easily with fertile lithospheric mantle (Figure 12b).In this situation, the large mantle plume with a diameter of 400 km and temperature anomaly of ≥150 K leads to craton destruction.In the cases with slightly decreased plume diameter (e.g., 300 km) or temperature anomaly (100 K), the transitional mode with local craton destruction is predicted (Figure S18 in Supporting Information S1).It is worth noting that the nature of such a plume is at the high value end-member or even beyond the typical natural plumes in the upper mantle of the present Earth (Koppers et al., 2021;Leng & Liu, 2023), which indicates that the plume-induced craton destruction is also not easy.In the cases with depleted lithospheric mantle (Figure 12f), the craton destruction is even more difficult, which only occurs with a plume diameter of 400 km and temperature anomaly of ≥250 K.In the intermediate cases with layered lithospheric mantle, that is, depleted above and fertile below MLD, the craton destruction requires a similar parameter range as the fertile cases (cf.Figures 12b and 12d).However, the layered lithospheric mantle does not produce transition mode with local craton destruction, which indicates its more stable nature than the fertile cases.
As a summary, the interaction between the weak antigorite-induced MLD and the mantle plume has significant effects on craton destruction, whereas the depletion of lithospheric mantle could reduce these effects.In addition, critical plumes with high temperature anomaly and large size are still required for craton destruction, which represent the strict conditions.
In the above models, variable boundary conditions are applied at the model time of 500 Myr, that is, the late stage of model evolution with the occurrence of MLD.However, the MLD within cratonic lithosphere might experience various tectonic processes at its different stages of formation and evolution.Thus, an additional model is conducted with upwelling of mantle plume at the model time of 100 Myr, that is, the early stage of MLD formation.All the other parameters are identical to the model with antigorite-induced MLD in Figure 10 (right column).In this model, the mantle plume could foster the convective dripping at the base of lithosphere, and further affect the later evolution of craton and MLD.However, the lithospheric delamination does not occur, because the MLD is not entirely formed (Figure S19 in Supporting Information S1).At the early stage of MLD formation, this layer would not play an important role in craton instability under variable tectonic regimes; however, the occurrence of tectonic process could affect the nature of MLD during the formation, such as its depth.

Formation Mechanism and the Corresponding Strength of MLD
The roles of MLD in craton instability are systematically investigated by numerical models in this study.The MLDs with contrasting viscosities, that is, wet olivine-or antigorite-induced, are applied based on the "water collector" hypothesis (Fu et al., 2022).However, the formation mechanism of MLD is still widely debated.Regime diagrams of model results with variable MLD rheology, plume diameter and temperature, as well as variable degree of lithospheric mantle depletion.The reference density of ρ 0 = 3,300 kg/m 3 is used for the fertile lithospheric mantle, whereas ρ 0 = 3,250 kg/m 3 is applied for the depleted lithospheric mantle.In the layered case, the depletion (ρ 0 = 3,250 kg/m 3 ) is only applied for the top lithospheric mantle above the MLD and the lithospheric mantle below MLD is still fertile with ρ 0 = 3,300 kg/m 3 .The white stars represent the typical models as shown in Figures 10 and 11.
the MLD would have a relatively low viscosity at the dominant depth of 70-100 km, based on the "partial melting" hypothesis.However, partial melting is difficult in the MLD depth of 70-100 km, but may only occur below the depth of ∼150 km with sufficient water content of >0.1% or 1,000 ppm (Figure 13a), based on the comparison between the solidus of peridotite with variable water content (Litasov, 2011) and the typical geothermal structure of craton with surface heat flux of 40 mW m −2 (Lee et al., 2011).Such a high water content significantly exceeds the water capacity of craton lithospheric mantle at the relevant depth (Figure S1 in Supporting Information S1).Thus, the MLD may not be caused by partial melting, although the presence of melts could lead to significant reduction of viscosity (Karato et al., 2015).2. Elastically accommodated grain-boundary sliding (EAGBS).At modest temperatures, the EAGBS, that is, one mechanism of anelasticity, may occur.It could lead to the reduction of elastic modulus, and further result in the decrease of seismic velocity within the lithosphere.Hence, EAGBS could explain the occurrence of the MLD with significant Vs drop (Karato et al., 2015;Karato & Park, 2018).During this deformation process, the strain is rather small, which is comparable to the order of elastic strain.Therefore, the occurrence of EAGBS may not influence the long-term viscosity of cratonic lithosphere, and could not play a significant role in the craton instability.3. Seismic anisotropy, especially the radial anisotropy.A radially anisotropic layer, independent of the back azimuth, could cause the consistent negative velocity change; thus, it has been proposed as a possible formation mechanism of the MLD (Sun et al., 2018).However, the radial anisotropy is related to the tectonic history of cratons, and thus may not explain the global distribution of MLD (Selway et al., 2015).On the other hand, the occurrence of seismic anisotropy, which is attributed to the presence of lattice preferred orientations (LPO), could lead to grain-boundary sliding, and further result in a significantly anisotropic viscosity in the dislocation creep regime (Castelnau et al., 2008;Hansen et al., 2012).Therefore, the viscosity of MLD may decrease due to the change of the direction of strain during craton evolution, but it is hard for the seismic anisotropy to explain the global MLD. 4. Channel flow.The channel flow could create a shear zone with radial anisotropy at the depth comparable to the natural MLD, which could lead to a negative change of seismic velocity (Yang et al., 2023).In addition, the channel flow could reduce the grain size, and further facilitate the accumulation of fluid and melts.The seismic low velocity minerals may be produced during this process, which would further reduce the seismic velocity.Thus, the "channel flow" hypothesis combines the effects of anisotropy and low velocity minerals.However, it is intuitively difficult for the channel flow occurring within the global cratonic lithospheres.5. Accumulation of seismic low velocity minerals.The MLD may be induced by low velocity minerals within cratonic lithosphere, which is produced by metasomatism, such as amphibole, phlogopite and carbonate minerals (Aulbach et al., 2017;Eeken et al., 2018;Rader et al., 2015;Selway et al., 2015;Smart et al., 2021).Viscosity of variable low velocity minerals are calculated under the cratonic thermal structure with a typical surface heat flux of 40 mW m −2 .The relevant flow law parameters determined by experiments are summarized in Table S5 in Supporting Information S1 (Agard et al., 2016;Getsinger & Hirth, 2014;Getsinger et al., 2013;Hilairet et al., 2007;Hirth & Kohlstedt, 2003;Holyoke et al., 2013Holyoke et al., , 2014;;Kirby & Kronenberg, 1987; Rybacki  et al., 2006).The viscosity profiles indicate that the presence of low velocity minerals could decrease the viscosity significantly, either at a low strain rate of 10 −19 s −1 or at a high strain rate of 10 −14 s −1 .At the dominant depth of MLD, the viscosities of variable minerals follow the rule that: wet olivine > amphibole > phlogopite > carbonate minerals > antigorite (Figures 13c and 13d).The viscosities of wet olivine and antigorite characterize the highest and lowest end-members, respectively.
The abundances of low velocity minerals might be limited at the depth of MLD.The presence of 5%-10% phlogopite or 10%-15% carbonate minerals could cause a Vs drop of 2%-7%, to the MLD observations of 2%-12% (Rader et al., 2015).Based on the mantle xenoliths, ≤∼10% amphibole, <∼2% phlogopite and ≤∼0.2% magnesite could be contained within the cratonic lithosphere, and further lead to a Vs drop of 2%-3% at the depth of the MLD (Saha et al., 2021).For the "water collector" hypothesis, the maximum water capacity of mantle rocks is about 3% in the present model, while the serpentine could contain a large amount of water up to 13% (Hilairet et al., 2006).Thus, the MLD might be a partially serpentinized layer, and the abundance of  Pollack and Chapman (1977).The black dashed lines represent the solidus of peridotite with variable water content (Litasov, 2011).(b) The viscosity profiles of olivine with variable melt fraction.The viscosity is calculated at the second invariant of strain rate (   ) of 10 −19 s −1 (dashed lines) and 10 −14 s −1 (solid lines) under the geothermal structure with surface heat flux of 40 mW•m −2 .(c-d) Comparison of viscosities of variable low velocity minerals.The viscosity is calculated using the rheological parameters in Table S5 in Supporting Information S1 with the value of   = 10 −19 s −1 and 10 −14 s −1 , respectively.The solid lines represent the depth range for the stability of these minerals, while the dashed lines represent the depth range out of the stability field of minerals.The viscosity of wet olivine is generally the highest, whereas that of antigorite is the lowest at the dominant depths of the MLD (gray rectangle).
antigorite at the depth of MLD is relatively low, although this value is not well constrained.The natural MLD could only contain a certain amount of low velocity minerals, and thus its viscosity may be a combined effects of the minerals shown in Figure 13.
As a summary, variable formation mechanism of the MLD does not always cause a rheologically weak MLD, such as EAGBS, whereas several mechanisms, which could lead to a weak MLD, may not explain the global presence of the MLD within cratons, such as seismic anisotropy or channel flow (Table 2).The accumulation of low velocity minerals is most likely to form the global MLD with relatively low viscosity.Based on the comparison among variable low velocity minerals, the viscosity of MLD may vary from the relatively high viscosity induced by wet olivine to the rather low viscosity induced by antigorite.Thus, we propose that the numerical model results, applying the MLD with contrasting strengths in the present study, represent the end-member cases of the roles of MLD in craton instability.

Effect of MLD on the Craton Instability
The effects of the MLD on craton destruction have been investigated in the previous studies (Liu et al., 2018;Shi et al., 2021;Wang et al., 2018;Wang & Kusky, 2019), with a rheologically weak MLD decoupling the overlying and underlying lithosphere.Once a weak MLD links to the hot asthenospheric mantle due to the lithospheric thinning by small-scale convective erosion, such a MLD could foster lithospheric delamination, and further cause craton destruction (Liu et al., 2018;Wang & Kusky, 2019;Wang et al., 2018).In these models, a rather weak MLD at constant depth is essential, and its connection with the sub-plate asthenosphere is strongly dependent on the occurrence of small-scale, edge-driven convection.However, the relatively weak MLD induced by the accumulation of seismic low velocity minerals, is mostly limited by the P-T conditions of stability field of minerals and the thermal structure of craton (Figures 13c and 13d).When the convective dripping occurs at the base of lithosphere, the MLD would migrate upward due to the increase of temperature condition during lithospheric thinning (e.g., Figure 6).The connection of weak MLD with sub-plate asthenosphere is still difficult, unless a certain section of the lower lithosphere has been destroyed.Therefore, the large amount of melt is essential in the present models, which contributes to destroy a lower part of lithosphere and trigger craton destruction (e.g., Figures 7 and 10), similar to the previous study (Shi et al., 2020(Shi et al., , 2021)).
In the present study, the stability field of hydrous minerals at the depth of MLD has been considered.Model results indicate that craton with the presence of the MLD, even the weakest antigorite-induced MLD, could preserve its stability under most tectonic regimes, such as tectonic compression and the mantle flow traction.Only the connection of weakest MLD with sub-plate asthenosphere could lead to lithospheric delamination and even craton destruction under tectonic extension regime or the upwelling of mantle plume.The rapid formation and accumulation of sufficient amount of melts are essential for this process to destroy the underlying lithosphere and connect the MLD with asthenosphere, which require the large amount of extension or the mantle plume with high temperature anomaly and large size.
It is worth noting that craton destruction prefers to occur with the fertile lithospheric mantle, whereas the depleted cratonic lithospheric mantle with low density could increase the intrinsic buoyancy of lithosphere, which hinders the delamination process, and enhances the stability of craton.Collectively, the strict conditions for craton destruction are not widely applicable in the natural cases.Thus, the role of MLD in the craton destruction may be not as significant as previous studies, which explains the general stability of natural cratons even with widespread MLDs.

Model Limitations and Future Perspectives
In the present model, the effects of MLD on craton instability during subduction process has not been investigated directly, due to the limitation of model domain.However, the interaction between the MLD and subduction process has been divided into two parts.The first one is the water activity and (de)hydration process during subduction, which contributes to the formation of the MLD.This process is simplified by the hydrated blocks generation routine in the present numerical model.The second part is the evolution of continental lithosphere with subduction and its induced fluid/melt activity, which is represented by the craton evolution under variable tectonic regimes in the numerical model.For example, the craton may be subjected to compressive stress during flat subduction (Yan et al., 2020), which corresponds to the craton compression as discussed in Section 3.3.Meanwhile, back-arc extension could occur during steep subduction, and is comparable to the tectonic extension regime as shown in Section 3.2.Thus, based on the model results, the MLD within overriding cratonic lithosphere may affect the craton instability during subduction.
On the other hand, we focus on the roles of MLD in the instability of normal cratons with constant thickness, that is, craton interiors.Previous numerical studies indicated that the lithospheric delamination along the MLD prefers to occur at the craton margins, where the MLD links with hot asthenospheric mantle more easily (Liu et al., 2018;Wang & Kusky, 2019).In addition, the MLD has been detected within both Archean and Proterozoic cratons (Figure 1a), and the thickness of craton varies with age (Artemieva, 2006;Artemieva & Mooney, 2001).Thus, the roles of the MLD in the craton instability with heterogeneous thickness require further studies.

Conclusion
Systematic numerical models are applied to investigate the roles of MLD in craton instability under variable tectonic regimes, including the stable craton, tectonic extension, tectonic compression, mantle flow traction, as well as the upwelling of mantle plume.The hydrated blocks generation routine is applied to simulate the formation of MLD.The MLDs with contrasting viscosities, that is, the wet olivine-induced MLD and antigorite-induced MLD, are compared in this study.The main conclusions are shown below: 1.The antigorite-induced MLD with lowest viscosity and wet olivine-induced MLD with relatively high viscosity represent two end-member cases of MLD with contrasting rheological strengths.2. The wet olivine-induced MLD could not lead to lithospheric delamination and craton destruction under all the tested tectonic regimes.3. The weak antigorite-induced MLD could localize large strain and decouple the overlying and underlying lithosphere significantly; despite this, the lithospheric delamination still requires additional conditions.4. Craton destruction only occurs with the connection of the weak antigorite-induced MLD and the sub-plate asthenosphere during craton extension or mantle plume activity.The partial melting process during large amount of extension or upwelling of mantle plume with high temperature anomaly and large size is a key condition.In addition, the depleted cratonic lithospheric mantle with low density would increase the intrinsic buoyancy of lithosphere, and inhibit the lithospheric delamination and craton destruction.
The systematic numerical studies indicate that the role of MLD in the craton destruction is not as significant as previously considered in the models, which requires additional strict conditions that are not widely satisfied on the Earth.Thus, it explains the generally stable cratons even with clearly observed MLDs.This work was supported by the NSFC project (42225403), National Key R&D Program of China (2022YFF0801001) and the Fundamental Research Funds for the Central Universities.Numerical simulations were run on the clusters of National Supercomputer Center in Guangzhou (Tianhe-II) and the GeoModeling cluster in UCAS.Yang Wang is greatly acknowledged for the helpful discussion.The constructive reviews by two anonymous reviewers improved the paper.

Figure 1 .
Figure 1.(a) Distribution and depth of the continental mid-lithosphere discontinuity (MLD) based on multiple geophysical observations, with the detailed data shown in Table S1.The distribution of cratons is shown in gray (modified from Bally et al. (2020)).(b) Frequency histogram of the depth of natural MLD, with the data grouped by cratons.

Figure 2 .
Figure2.Initial model design.(a) Composition field with the white dashed lines for the depth of 70 and 100 km, respectively.The yellow dashed lines represent the 410 and 660 km discontinuity, respectively.The colors in (a) are used for rock types as specified in the color grids: 1-sticky air (not shown in the figure); 2-sediment; 3-continental upper crust; 4-continental lower crust; 5-dry lithospheric mantle; 6-dry sub-lithospheric mantle; 7-hydrated lithospheric mantle; 8-hydrated sub-lithospheric mantle; 9-serpentinized mantle; 10-partially molten mantle.The temperature profile at the location of X = 600 km is shown in (c).(b) The effective viscosity field, with a vertical profile at the location of X = 600 km as shown in(d).(e-f) The viscosity profiles of olivine and antigorite, calculated using the rheological parameters in Table1, with the strain rate (   ) of 10 −19 s −1 , 10 −14 s −1 , respectively.The gray rectangles represent the dominant MLD depth range of 70-100 km as shown in Figure1.

Figure 3 .
Figure 3. Tectonic background and conceptual sketch of the hydrated blocks generation routine in the numerical model.The water could be carried into the deep mantle during subduction, and be released from subducting slab due to the continuous heating from surrounding mantle (a, b).In the previous numerical model (e.g.,Li et al., 2019), when the dehydration process occurs (c), the aqueous fluid markers are generated in the grid (d), then migrate upward and hydrate the surrounding rocks (f).In the present model, a new generation routine of hydrated blocks (e) is applied to simplify the simulation of whole process of (de)hydration in the deep mantle.The gray dots represent mantle rocks, while the blue dots represent the addition of released water.The yellow and green grids represent the dry and hydrated mantle domains, respectively.

Figure 4 .
Figure 4. Reference model evolution with the wet olivine-induced MLD.(a-e) Composition field evolution at the model time of 10, 100, 250, 400, 550 Myr, respectively.The white solid line represents the isotherm of 1,350°C.The black dots represent the hydrated blocks generated at the depth of 410 km during the model evolution.The colors in the composition field represent the rock types as shown in Figure 2a.(aʹ-eʹ) Effective viscosity field evolution at 10, 100, 250, 400, 550 Myr, respectively.The black solid line represents the viscosity contour of 10 22 Pa•s.

Figure 6 .
Figure 6.The horizontally averaged water content (a), effective viscosity (b), Vs (c) and Vs drop (d) as a function of depth at 550 Myr.(e) Evolution of the depth range of the MLD.The solid and dashed lines represent the top and bottom boundaries of MLD, respectively.(f) Evolution of the depth of thermal LAB.The blue and yellow lines represent the model evolution in Figures 4 and 5, respectively.The gray rectangle represents the dominant depth range of natural MLDs (Figure 1).

Figure 7 .
Figure 7. Two representative model results during craton extension.Both models have the same extension rate of 0.5 cm/ yr (half rate).(a-d and aʹ-dʹ) Composition and effective viscosity field with wet olivine-induced MLD.(f-i and fʹ-iʹ) Composition and effective viscosity field with antigorite-induced MLD at the model time of 500, 510, 517, 520 Myr, respectively.Before the model time of 500 Myr, no boundary extension is applied; the model evolves similarly as the reference models in Figures 4 and 5. (e, j) Time-dependent evolution of the depth of compositional LAB with variable extension rate.Craton destruction only occurs with the presence of a weaker MLD and larger amount of extension (i.e., extension rate multiplied by extension time).

Figure 8 .
Figure 8. Two representative model results during craton compression.Both models have the same convergence rate of 0.5 cm/yr (half rate).(a-c and aʹ-cʹ) Composition and effective viscosity field with wet olivine-induced MLD.(e-g and eʹ-gʹ) Composition and effective viscosity field with antigorite-induced MLD at the model time of 500, 510, 520 Myr, respectively.(d, h) Time-dependent evolution of the depth of compositional LAB with variable convergence rate.

Figure 9 .
Figure 9. Two representative models under mantle flow traction with a sub-lithosphere velocity of 5 cm/yr.The wet olivine-induced MLD (left column) and the antigorite-induced MLD (right column) are applied, respectively.(a-c, aʹ-cʹ) Composition field evolution at 500, 510, 520 Myr, respectively.Effective viscosity field are shown for the model time of 520 Myr (d, dʹ).

Figure 10 .
Figure 10.Two representative models with interaction between MLD and plume.The mantle plume is generated at the depth of 410 km at the model time of 500 Myr.The diameter of plume is 200 km and the temperature of plume is 2,123 K. (a-d and aʹ-dʹ) Composition and effective viscosity field with wet olivine-induced MLD.(e-h and e'-h') Composition and effective viscosity field with antigorite-induced MLD at 500.02, 500.2,500.5, 500.8Myr, respectively.

Figure 11 .
Figure 11.Model evolution with the depleted lithospheric mantle (ρ 0 = 3,250 kg/m 3 ) and the antigorite-induced MLD.The properties of mantle plume are identical as the model in Figure 10, with the diameter of 200 km and temperature of 2123 K. (a-c) Composition field evolution.(aʹ-cʹ) Effective viscosity field evolution at 500.05, 501, 510 Myr, respectively.

Figure
Figure12.Regime diagrams of model results with variable MLD rheology, plume diameter and temperature, as well as variable degree of lithospheric mantle depletion.The reference density of ρ 0 = 3,300 kg/m 3 is used for the fertile lithospheric mantle, whereas ρ 0 = 3,250 kg/m 3 is applied for the depleted lithospheric mantle.In the layered case, the depletion (ρ 0 = 3,250 kg/m 3 ) is only applied for the top lithospheric mantle above the MLD and the lithospheric mantle below MLD is still fertile with ρ 0 = 3,300 kg/m 3 .The white stars represent the typical models as shown in Figures10 and 11.

Figure 13 .
Figure 13.(a) Geothermal structures of cratons with variable surface heat flux as proposed byPollack and Chapman (1977).The black dashed lines represent the solidus of peridotite with variable water content(Litasov, 2011).(b) The viscosity profiles of olivine with variable melt fraction.The viscosity is calculated at the second invariant of strain rate (   ) of 10 −19 s −1 (dashed lines) and 10 −14 s −1 (solid lines) under the geothermal structure with surface heat flux of 40 mW•m −2 .(c-d) Comparison of viscosities of variable low velocity minerals.The viscosity is calculated using the rheological parameters in TableS5in Supporting Information S1 with the value of   = 10 −19 s −1 and 10 −14 s −1 , respectively.The solid lines represent the depth range for the stability of these minerals, while the dashed lines represent the depth range out of the stability field of minerals.The viscosity of wet olivine is generally the highest, whereas that of antigorite is the lowest at the dominant depths of the MLD (gray rectangle).

Table 1 Viscous Rheological Flow Laws Used in the Numerical Experiments Figure 2. et
al., 2022).In the present numerical model, hydrated blocks are generated at the constant depth of 410 km.The size (S h ) of hydrated blocks is 5 km × 3 km, with the water content (W h ) of ∼5%.The averaged generation timestep (t h ) of hydrated blocks is 0.05 Myr.The generation rate of hydrated block (

Table 2
Formation Mechanism of MLD and Its Rheological Strength Variation