Continental Mid‐Lithosphere Discontinuity: A Water Collector During Craton Evolution

Continental mid‐lithosphere discontinuity (MLD) seems to be ubiquitous beneath cratons around the world with the dominant depth of 70–100 km, and is characterized by a shear‐wave velocity drop of 2%–7%, according to geophysical observations. The MLD is often considered to be related to a rheologically weak layer; however, the mechanism of MLD formation is widely debated. In this study, we have conducted systematic numerical modeling with a new and simplified deep hydration method to study the dynamics of craton evolution and MLD formation. The results indicate that the MLD may be induced by slow hydration processes within the mantle lithosphere during craton evolution. The top boundary of this hydrated layer is characterized by high water content and low shear‐wave velocities, and is consistent with the depth and properties of natural MLD observations. Thus, we propose that the MLD may act as a water collector during craton evolution.

• The dynamics of craton evolution and mid-lithosphere discontinuity (MLD) formation is investigated by numerical models with a new and simplified deep hydration method • The MLD may be induced by slow hydration within the mantle lithosphere and acts as a water collector during long-term craton evolution • The MLD represents the top boundary of a hydrated layer with high water content and low shear-wave velocity within cratonic lithosphere

Supporting Information:
Supporting Information may be found in the online version of this article.
The MLD may be corresponding to a layer with high water content within the mantle lithosphere produced by fluid metasomatism. The hydration process could be induced by a series of tectonic processes including plate subduction and upwelling of mantle plumes during the evolution (Zhu et al., 2021). If enough water is released during dehydration, excess water can migrate upward and further interact with the overlying lithosphere  Table S1 in Supporting Information S1. (a) Global distribution and depth of the MLD. (b) Frequency histogram of the MLD depths, with the peak values of 70-100 km (gray dashed box). (c) Frequency histogram of the shear-wave velocity (Vs) drops of the MLD, with the dominant values of 2%-7% (gray dashed box). (Li, 2020(Li, , 2022. The released water may hydrate the lithospheric mantle and then produce a high-water-content layer, due to the high water capacity of mid-lithospheric mantle with low temperature and pressure. In the 200-km-thick continental lithosphere, the water capacity of mantle rock increases significantly above the depth of ∼100 km due to the appearance of hydrous minerals, whereas it is nearly zero below this depth with limited water only in the nominally anhydrous minerals (NAMs) ( Figure S1 in Supporting Information S1). The high water content will lead to a significant drop of shear-wave velocity within the mantle lithosphere ( Figure S2 in Supporting Information S1), which may explain the MLD detected by seismic methods.
In this study, we aim to explain the global occurrence of the MLD based on a systematic thermomechanical model with thermodynamic fluid activity. A new and simplified deep hydration method is applied in the present model. The factors influencing the formation of the MLD are investigated, and the depth and shear-wave velocity drop of the MLD in the model are further compared with natural observations.
Numerical models are configured in a 2-D Cartesian box of 600 km × 420 km. The initial model contains a 10-km thick "sticky air" layer, a 200-km-thick continental lithosphere and a sublithospheric mantle layer ( Figure S3 in Supporting Information S1). The detailed model setup and boundary conditions are shown (Text S2 and Figure  S3 in Supporting Information S1).

Hydrated Block Generation Routine
When dehydration occurs in the deep mantle, excess water is released and moves upward (Figures 2a and 2b). The propagation of released water is modeled by the migration of fluid markers (Figures 2c and 2d), each of which stores a certain amount of water, and moves with a calculated velocity . When a fluid marker enters the grid of dry mantle rocks or hydrated ones but still below water capacity ( Figure S1 in Supporting Information S1), hydration will occur in this grid ( Figure 2e). The water content of hydrated mantle depends on its water capacity and the water content of the fluid marker. In the meantime, with the occurrence of hydration, the water content of the fluid marker decreases gradually during the migration, and finally becomes zero. At this moment, the fluid marker is eliminated, which means that there is no excess water.
A new and simplified hydrated block generation routine ( Figure 2) is applied instead of simulating the whole processes of (de)hydration during subduction or mantle plume, because the complicated and regional dehydration over a long geological history cannot be well constrained. During the present model evolution, hydrated blocks are generated at random horizontal locations with a constant depth in the model domain and with random time steps, representing the spatial and temporal uncertainty of the occurrence of dehydration in the deeper mantle. The water content of hydrated blocks represents the amount of released water, and the migration of blocks is modeled by generating fluid markers with a certain water content in the block area, where the number and location of fluid markers are the same as those of rock markers ( Figure 2f).

Model Results
In the reference model, hydrated blocks are generated at random horizontal locations with a constant depth of 404 km in the model domain. The size (S h ) of hydrated block is 2 km × 2 km, with the water content (W h ) of 6 wt%, and the average generation timestep (t h ) of hydrated block is 0.05 Myr. Before the model time of 323 Myr, the continental lithosphere thickens gradually ( Figure S4 in Supporting Information S1). A layer with high water content forms within the mantle lithosphere, due to the successive generation of hydrated blocks (Figures 3a and 3b). Meanwhile, the water content of mantle rock along the migration pathway of fluid markers increases (Figures S5, S6 in Supporting Information S1), leading to a reduction of the viscosity by about one order of magnitude (Figures 3a', 3b', Figure S7 in Supporting Information S1). At the model time of 323 Myr, the layer with high water content is almost horizontally continuous, and the top boundary of this layer moves upward gradually. The continental lithosphere is no longer thickening; instead, convective dripping develops at the base of lithosphere, due to the effect of Rayleigh-Taylor instability ( Figure 3c). Correspondingly, the mantle effective viscosity decreases rapidly at the location where the convective dripping occurs (Figure 3c'). At the final model stage of 551 Myr, a continuous high-water-content layer is formed within the lithosphere. The top boundary of this layer, representing the MLD, has a horizontally averaged depth of 83.8 km (Figure 3d). The detailed model evolution is discussed (Text S3 and Figures S4-S7 in Supporting Information S1, Movies S1-S4).
During model evolution, the horizontally averaged depth of the MLD decreases slightly, that is, the depth of the MLD becomes a bit shallower with the increase of water content in the lithosphere. The depth of the MLD in the numerical model is in good agreement with the dominant depth of natural MLDs (70-100 km) given by seismic observations. In addition, the lithospheric mantle is difficult to maintain a cooling-induced steady thickening, but exhibits a "thickening-dripping-thickening" pattern and keeps a dynamically balanced lithospheric thickness ( Figure 3e, Figures S8, S9 and Text S3 in Supporting Information S1).
The MLD is characterized by obviously low shear-wave velocity (Vs) within the mantle lithosphere. In order to directly compare with the geophysical observations, the water content and its corresponding Vs at the model time of 551 Myr are plotted (Figure 4). The Vs interpolated from the Vs database ( Figure S2 in Supporting Information S1) shows that high water content causes the reduction of Vs at the depth of the MLD (Figure 4b), and the horizontally averaged Vs in the model is significantly lower than the normal Vs without MLD (Figure 4d). In this study, the database of Vs at different temperatures, pressures and water contents is calculated using Perple_X    Figure 1). (Connolly, 2005(Connolly, , 2009Holland & Powell, 1998), with the rock composition of mean Archean subcontinental lithospheric mantle given by Griffin et al. (2003) (Figure S2 and Table S4 in Supporting Information S1). Compared with the normal value, the Vs at the depth of the MLD is reduced by about 1%-6% ( Figure 4c). The maximum horizontally averaged Vs drop of the whole model domain is ∼3%, which is consistent with the dominant range of Vs drop (2%-7%) obtained from geophysical observations (Figure 4e).

Discussion
The time-dependent depth evolution of the hydrated layer in mid-lithosphere is shown in Figure 5a. The top boundary of this layer may represent the occurrence of the MLD, with the accumulation of hydrous minerals by metasomatism. However, it is still controversial about the specific minerals that can accumulate at the MLD depth. The presence of amphibole, phlogopite or carbonates has been reported in mantle xenoliths at the depth of the MLD, and can cause shear-wave velocity drop (Eeken et al., 2018;Rader et al., 2015). Amphibole is most likely to form MLD, due to its P-T condition of stability (<1,070°C, <3 GPa) corresponds to the dominant depth of natural MLDs (Selway et al., 2015). Nevertheless, the presence of phlogopite can limit the amphibole-stability, and explain the MLD with a wide depth range (Aulbach et al., 2017;Saha et al., 2018). However, the xenoliths from the Siberian craton show that no rock rich in amphibole or phlogopite (Liu et al., 2022). Amphibole and phlogopite may be not ubiquitous within cratons based on mantle xenoliths data (Selway et al., 2015). Due to such kind of complexity, the minerals, associated with the formation of the MLD, require further constraints from geological observations. In the present study, we focus on the dynamic mechanism of MLD formation, and do not aim to constrain the specific minerals for the MLD.
The numerical models indicate that the MLD is formed gradually over a long geological history, and acts as a water collector during craton evolution. In the present model, there are two variable parameters: (a) the rate of hydrated block generation, which represents the amount of water that migrates into the asthenosphere from below per unit time, and corresponds to the dehydration processes in the deeper mantle ( ℎ = ℎ × ℎ ℎ , where S h is the size, W h the water content, and t h the average generation timestep of hydrated blocks); (b) the full time of model evolution (t f ). In the reference model, R h = 4.8, which is hard to be well constrained in nature due to the complicated and spatially variable tectonic events of dehydration beneath cratons. The entire model time, t f = 551 Myr, indicating that only the Phanerozoic period is considered in the reference model. On the other hand, a series of comparative models are conducted (Figure 5), in which the rate (R h ) multiplied by time (t f ) is constant, indicating that the total water content is identical in all the models. As a result, the MLD in these models are consistently resulting at ∼80 km depth (Figure 5c), with similar shear-wave velocity and velocity drop ( Figure S10 in Supporting Information S1). It indicates that these two parameters (i.e., R h and t f ) can compensate each other, for example, if the hydration rate (R h ) is lower, the MLD can still be achieved by increasing the model evolution time (t f ). The water distribution of early Earth is unclear (Brown et al., 2020), and the hydration rate (R h ) may be greatly promoted by subduction process in the plate tectonics regime. However, the exact times of global cratonization of continent (Cawood et al., 2018;Zhai & Peng, 2020) and the onset of plate tectonics (Korenaga, 2013) are both debating. Thus, the time for the MLD formation cannot be constrained from the current study; but the numerical models do indicate that the MLD could be gradually formed during the slow and nonsteady state hydration process.
Recent studies also suggest that multiple MLDs exist in the cratonic lithosphere Sodoudi et al., 2013;Wirth & Long, 2014), and the detected MLD may not be a single discontinuity but consist of multiple small-scale thin layers (Kind et al., 2015;Sun et al., 2020). Calò et al. (2016) distinguished two discontinuities in the middle of cratonic lithosphere from seismic data, which are interpreted as the top and bottom boundaries of the low-velocity layer, respectively. The low-velocity layer within the lithosphere has a high-water-content layer in the numerical model, and the top boundary, that is, MLD, represents the front of fluid infiltration. In addition, the bottom boundary of the hydrated layer represents the maximum depth with high water capacity in the lithosphere. This depth varies with the thickness and thermal structure of cratons ( Figure S11 in Supporting Information S1). Below this depth, the rather limited water only exists in NAMs. During model evolution, the thickness of high-water-content layer is 22.2 km at the model time of 551 Myr, which increases with the amount of water entering the lithosphere. Therefore, the MLD may have a depth range of 50-160 km, and the thickness of the low-velocity layer reflects the variable thickness of cratonic lithosphere and effects of regional tectonic events during craton evolution.
The MLD forms gradually during craton evolution. In this study, slow hydration processes within the mantle lithosphere lead to the formation of the MLD, during which the cratonic lithosphere fails to thicken steadily, and convective dripping occurs at the base of lithosphere. On the other hand, the occurrence of the MLD, which is often considered as a rheologically weak layer, may further affect the stability of cratons. In the previous study, the MLD could decouple the overlying and underlying mantle, facilitate the delamination of cratonic lithosphere, and even lead to the craton destruction Shi et al., 2020;Wang et al., 2017Wang et al., , 2018. However, similar to the results of present numerical models (e.g., Figure 3, Figure S7 in Supporting Information S1), the strength of the MLD is generally not low enough to induce lithospheric delamination by itself (Liao et al., 2013;Wang & Kusky, 2019). This may account for the long-term stability of most cratons, although the MLD exits within them. Thus, the roles of the MLD in the craton destruction require further evaluations.

Conclusions
In this study, systematic numerical models have been conducted with a new and simplified deep hydration method to study the dynamics of craton evolution and MLD formation. We propose that the MLD may be formed by slow lithospheric hydration processes during the upwelling of aqueous fluid, released from plate subduction or mantle plume in the deep mantle. The MLD represents the top boundary of the hydrated layer within the mantle lithosphere, and thus acts as a water collector during long-term craton evolution. The depth of hydrated layer represents the favorable P-T condition for the occurrence of hydrous minerals, that is, hydrous mineral absent below this depth in the craton lithosphere due to high P-T condition. The variable depth of the MLD and the thickness of low-velocity layer are affected by the thickness and tectonic history of cratons. However, the MLD cannot induce lithospheric delamination by itself, and may require additional perturbations to facilitate craton destruction.

Data Availability Statement
The water capacity database of Figure S1 in Supporting Information S1 and shear-wave velocity (Vs) database of Figure S2 in Supporting Information S1 are computed by Perple_X (Connolly, 2005(Connolly, , 2009. The figures of numerical models are produced by Matlab (https://ww2.mathworks.cn/products/matlab.html) and further compiled by Adobe Illustrator (https://www.adobe.com/cn/products/illustrator.html). The related data are provided in the public repository of Zenodo (https://doi.org/10.5281/zenodo.7309267).