The Mid‐Lithospheric Discontinuity Caused by Channel Flow in Proto‐Cratonic Mantle

Global geophysical observations show the presence of the enigmatic mid‐lithospheric discontinuity (MLD) at depths of ca. 80–150 km which may question the stability and internal structure of the continental lithosphere. While various mechanisms may explain the MLD, the dynamic processes leading to the seismic observations are unclear. Here we present a physical mechanism for the origin of MLD by channel flow in the cratonic mantle lithosphere, triggered by convective instabilities at cratonic margins in the Archean when the mantle was hot. Our numerical modeling shows that the top of the frozen‐in channel flow creates a shear zone at a depth comparable to the globally observed seismic MLD. Grain size reduction in the shear zone and accumulation of percolated melts or fluids along the channel top may reduce seismic wave speeds as observed below the MLD, while the channel flow itself may explain radial anisotropy of seismic wave speeds and change in direction of the seismic anisotropic fast axis. The proposed mechanism is valid for a broad range of physically realistic parameters and that MLD may have been preserved since its formation in the Archean. The intensity of the channel flow ceased with time due to secular cooling of the Earth's interior. The new mechanism may reshape our understanding of the evolution and stability of cratonic lithosphere.

the North American craton, joint inversion of the high-frequency P wave coda with long period surface wave dispersion data identified not only a velocity discontinuity with a 2.1% Vs decrease at depths of ca. 100-120 km beneath the cratonic stations, but also a positive intra-lithospheric discontinuity at around 150-170 km depth with a 3.6% Vs increase (Calo et al., 2016).
The intra-lithosphere seismic discontinuity with a negative velocity converter at depths ca. 70-120 km globally under stable continents has also been identified by high-frequency seismic receiver function (RF) studies Rychert & Shearer, 2009). Similar studies, largely based on S-RF (S wave RF), reported a negative converter at ca. 60-85 km depth beneath the Precambrian Western Australia (Sun et al., 2018) and typically at ca. 80-120 km depth beneath the cratonic North America (Chen et al., 2018). Beneath the Kalahari craton, some S-RF studies identified a shallow (80-90 km depth) negative converter found in other S-RF studies (Rychert & Shearer, 2009), while other S-RF studies found a prominent interface with a 4.5% reduction in seismic velocity at ca. 150-160 km depth (Savage & Silver, 2008). Notably, the deep converter corresponds to the top of melt-metasomatized mantle beneath the Kaapvaal craton (Griffin et al., 2003;Figure 1b).
Vertical stratification of the cratonic lithosphere with the transition between compositionally distinct layers at ca. 120-150 km depth is revealed by mantle xenolith studies in the Archean-Paleoproterozoic provinces of the Slave craton and the Baltic Shield (Kopylova & Caro, 2004). Vertical stratification of the cratonic lithospheric mantle is also imaged in seismic azimuthal anisotropy studies of the North American and the Australian cratons, which show a boundary between mantle layers with different fast directions at ca. 100-150 km depth (Sun et al., 2018;. Likewise, strong electrical anisotropy has been reported for the Superior Province with the top of the anisotropic layer at ca. 70-100 km depth (Mareschal et al., 1995), while in the Slave craton a pronounced decrease in electrical conductivity by 0.5 log unit at 135-150 km depth (Jones et al., 2003) closely correlates with the compositional boundary established from mantle xenoliths (Kopylova & Russell, 2000).
This brief overview demonstrates the presence of an intra-lithosphere heterogeneity at a depth of around 100 km in stable continental lithosphere globally, although with regional depth variations (Figure 1a). It is, however, still debated if the top of an intra-lithospheric seismic refraction LVZ and negative seismic convertors imaged by RFs (later referred to as the mid-lithospheric discontinuity, MLD) image the same intra-lithospheric transition zone, and if anisotropic layering and metasomatic imprint on composition of the cratonic lithospheric mantle correspond to the same feature.
Various models have been brought forward to explain the MLD. They include: (a) rheological layering at temperature close to the solidus in the presence of melt/fluids (Thybo, 2006), (b) compositional layering due Figure 1. Compilation of depth to intra-lithosphere discontinuities in cratons (a) and thickness of the anomalous intra-lithosphere layer below these discontinuities (b) from seismic, MT, and xenolith data (see references in Supporting Information S1).
These mechanisms are under discussion. First, it is debated if solidus temperature may be reached at mid-lithosphere depths below stable continents even for a fertile and volatile-rich composition (Karato et al., 2015). Second, compositional layering may not necessarily relate to the seismic MLD. Minerals rich in carbon or amphibole may produce a >5% seismic wave velocity drop (Selway et al., 2015), but mantle xenoliths suggest that such minerals may not occur ubiquitously in cratonic lithosphere (Rader et al., 2015;Selway et al., 2015). Third, an analysis of feasible mechanisms for a seismic velocity change in S-RF suggests that radial anisotropy could be a strong candidate to explain the seismic MLD (Selway et al., 2015). However, this model requires a global operation of a geological process which produces the same anisotropic layering pattern at the same depth range in every craton, and such process is yet unknown. Azimuthal anisotropy may either increase or decrease seismic velocity depending on the back azimuth of seismic waves, but this is not revealed by long-range controlled-source seismic profiles in the cratons of North America, Siberia, and Baltica (Thybo, 2006). Elastically accommodated grain-boundary sliding and mantle hydration, for example, by melt metasomatism, seem to predict a velocity drop similar to observed at the MLD (Karato et al., 2015). However, the model indicates no relationship between a reduction in seismic velocity at the MLD and anisotropy variations within the lithospheric mantle.
In summary, mechanisms for the presence of the MLD should explain the following features: (a) a seismic velocity interface observed both in seismic controlled-source and RF models, (b) a transition between layers with different directions of seismic, and possibly electrical, anisotropy, and (c) the ubiquitous appearance of the MLDs in cratons.
Here, we propose a new mechanism for the origin of the MLD by frozen-in channel flow formed at the early cratonic evolution. This mechanism may resolve the controversy between various MLD models. We use a 2D state-of-the-art geodynamic modeling software "Underworld" (Moresi et al., 2007) to demonstrate that channel flow may explain the effects observed in various geophysical data and interpreted as the MLD. The numerical model predicts the formation of a rheologically weak mid-lithosphere zone in the cratons globally, with the associated changes in seismic velocity and formation of anisotropic layering. We propose that a global mid-lithosphere weakness zone may have formed in the Precambrian at the early stages of cratonization by channel flow in the mantle. We argue that the intensity of the channel flow mechanism ceased with time due to secular cooling of the Earth's interior.

Thickness of Chemical and Thermal Boundary Layers
We perform numerical experiments (see Methods in Supporting Information S1) to examine the long-term evolution of Archean lithosphere for a range of possible thicknesses of the proto-lithospheric mantle thickness ( Figure 2a and details below). Edge-driven convection around the transition between thick and thin lithosphere triggers a downwelling mantle flow, which creates a horizontal pressure gradient in the lower part of cratonic lithosphere. This, in turn, causes a Poiseuille-type channel flow which produces features similar to the MLD.
The model time zero corresponds to the Archean formation of a compositionally low-density continental lithospheric root below the craton, which forms the chemical boundary layer (CBL) composed of a buoyant material (Figure 2b), left behind after an extensive melt extraction from the mantle (Carlson et al., 2005). The base of the CBL defines the chemical LAB (CLAB). Besides, the lithosphere of both the thick and thin proto-CBL domains makes the conductive thermal boundary layer (TBL) atop the convecting mantle, and the boundary between the conductive and convective mantle defines the thermal LAB (TLAB).
Initially, CLAB is deeper than TLAB (Figure 2c). This situation corresponds to the start of cratonic evolution in the Archean (time zero), when a high concentration of radioactive elements in the chemical lithosphere of cratons and a high Archean geothermal gradient (Herzberg et al., 2010) cause partial melting of lower lithospheric mantle above the CLAB, so that the base of the conductive layer TLAB in the hot Earth is shallower than CLAB (Figure 2c). This creates a gravitationally unstable situation with a thick depleted, low-density, and at the same time warm cratonic root above undepleted convective mantle ( Figure 2b). The situation may have analogy with modern oceanic plateau when intensive mantle melting in regions with hot mantle upwellings may produce a low-density/low-viscosity layer. Density inversion in Archean mantle may also be caused by subduction-induced viscous underplating of depleted oceanic sublithospheric mantle (Perchuk et al., 2020). The heterogeneity of mantle upwelling or proto subduction can cause thickness variation of the CBL in the mantle.
Our model setup at time zero differs from a xenolith-constrained model for cratonic lithosphere with shallower CLAB than TLAB based on the presence of fertile peridotites in the lower portion of the conductive lithosphere (Lee et al., 2005) as commonly used in geodynamic models of cratonic stability (Cooper et al., 2004;Lenardic et al., 2003). In our interpretation, the situation with shallower CLAB than TLAB develops during cratonic evolution through thinning of the CBL and thickening of the TBL. Continuous interaction between mantle convection and the lithospheric root leads to metasomatic modification of the lower lithospheric portion (Aulbach et al., 2017;Griffin et al., 2003), thereby reducing the thickness of the CBL. At the same time, secular cooling of the mantle and radiogenic depletion of the lithospheric root increase the thickness of the TBL.

Model Geometry and Mantle Thermal Structure
The model (4,000 km × 660 km) includes a 40 km thick crust, lithospheric mantle with different thicknesses in the left domain (x = 0-2,400 km) and the right domain (x = 2,600-4,000 km) of the model, and sublithospheric mantle down to 660 km depth. At model time zero (Archean formation of a craton), the thickness of the CBL  (Table S1 in Supporting Information S1), which includes a thick proto-lithosphere (x = 0-2,400 km) and a thin proto-lithosphere (x = 2,600-4,000 km) area with a linear transition between the two domains. (a) The model includes three layers: crust-top 40 km; mantle chemical boundary layer (CBL) with the base (chemical LAB, CLAB) at a depth of 235 km in the thick domain and 95 km in the thin domain, underlain by asthenospheric mantle down to 660 km depth. Circular markers in the upper half of the model (at x = 500-3,000 km) trace mantle deformation. (b) Initial in situ density profiles at thick (x = 2,300 km, blue) and thin proto-lithosphere domain (x = 2,800 km, red). In the sublithospheric mantle, a temperature-controlled density decrease with depth leads to convective instability (gray arrows). The lithostatic pressure difference between the two profiles is illustrated by the orange line. (c) Initial temperature profile at time zero with a linear temperature increase from 0°C at the model top to a potential temperature of 1,500°C at depth 90 km (thermal LAB, TLAB), with further increase in deeper mantle with an adiabatic gradient of 0.3°C/km. is 195 km in the thick proto-lithosphere domain and 55 km in the thin proto-lithosphere domain, with a linear change in the transition zone between the two domains. This model is an "end-member model" which demonstrates the formation of the channel flow ( Figure 2, see Table S1 in Supporting Information S1 for all model parameters) due to a big difference (140 km) of CBL thickness in the model. Models with smaller differences (25-85 km) are investigated in the next sections.
Our assumption, that initially the CLAB is deeper than the TLAB, implies the presence of a partially molten material in the lower part of the Archean CBL (i.e., the layer between the TLAB and the CLAB) which may form a secondary convection system within the chemical lithosphere root with complex geotherms (Sleep, 2003a(Sleep, , 2003b. We simplify the geotherm shape by (a) assuming a continuous gradient of 0.3°C km −1 in all convecting mantle layers, including the continental root (Sleep, 2003b) and (b) initially placing the TLAB at 90 km depth (Herzberg et al., 2010), that is maintaining CBL > TBL in both thick and thin proto-CBL domains ( Figure 2c). Therefore, initially the thick domain has a 145 km vertical difference between the CLAB (at 235 km depth for a 40 km thick crust) and the TLAB, while in the thin domain this vertical difference is 5 km only. Latent heat due to partial melting in the continent root gives rise to the complexity of the temperature distribution in the continent root. The effect of latent heat is neglected in the reference model but is discussed in Section 3.2.4.
We assume that the CBL and sublithospheric mantle have different density and heat production ( Figure 2b and Table S1 in Supporting Information S1): depleted cratonic lithospheric mantle with density of 3,300 kg m −3 is 30 kg m −3 (1%) less dense than sublithospheric mantle (Carlson et al., 2005), and has much higher radioactive heat production (0.2 versus 0.02 μW/m 3 ) than present values. The presence of a compositionally depleted, low-density cratonic lithospheric root is required to initiate the process, since the CLAB perturbation at the thick-to-thin proto-CBL transition triggers convective instabilities in the sublithospheric mantle (Section 3.1).
At model time zero the mantle potential temperature is 1,500°C (Figure 2c), that is ca. 200°C higher than at present (Nisbet et al., 1993). Due to secular cooling of the Earth and high rate of radioactive heat generation in Archean lithosphere, a 200°C lithosphere temperature anomaly converges to the same geotherm after 1 Gyr (Michaut & Jaupart, 2007). We incorporate this effect through a time-dependent heat generation model. At time , the rate of radioactive heat generation is = 0 − , where 0 is the concentration of heat producing elements at model time zero, and is the mean lifetime. The four major radioactive nuclides 40 K, 232 Th, 235 U, and 238 U with different half-life times are present in different proportions which leads to isotopic heat production rates with one order of magnitude difference (Van Schmus, 1995). Since the total radioactive heat production in the young Earth was dominated by 40 K and 235 U (see Figure S1 in Supporting Information S1), we adopt a half-life time of 1 Gyr for the radioactive nuclides. Therefore, at time zero the heat generation rate is about 4 times higher than present cratonic values (Artemieva & Mooney, 2001). The radioactive heat sources are modeled as Lagrangian particles advecting through time. Temperature at the lower boundary decays with time at a rate of 0.1°C/Myr from the initial temperature of 1,698°C at time zero (Figures 2a and 2c).

Formation of the Channel Flow
The density decreases with depth due to increasing temperature creates density inversion in the mantle and triggers the development of convective instabilities in the sublithospheric mantle ( Figure 2b). Their intensity is different below the thick and the thin proto-CBL domains, despite at time zero the model has a horizontally homogeneous temperature field. The only lateral heterogeneity in the mantle is associated with the transitional zone (x = 2,400-2,600 km) from the thick-to-thin proto-CBL domain, where the depth of the CLAB has a steep change of 140 km over a 200 km distance.
Due to its larger thickness and larger temperature difference between the top (CLAB) and the bottom (660 km depth), the sublithospheric mantle outside the thick proto-CBL domain (x > 2,600 km) has a higher Rayleigh number than below the thick proto-CBL domain (x < 2,400 km) ( Figure 2c). This triggers edge convection beneath the transitional zone (x = 2,400-2,600 km) in the sublithospheric mantle (Figures 2b and 3a), which reduces the pressure in the mantle around the thick proto-CBL domain edge as compared to the pressure within the thick proto-lithosphere far from the transition zone (x = 500-2,200 km) (Figure 3a, top). As a result, the horizontal pressure gradient in the lower chemical lithosphere creates channel (Poiseuille) flow at depths between ca. 100 and 230 km (Figures 3a and 3a′).  ; note that the scale varies with time due to mantle cooling and viscosity increase. Orange lines-temperature isotherm of 1,400°C, which closely follows the TLAB. Small circles in the upper mantle placed with a 500 km distance step illustrate mantle deformation. (a) Convective instability at the thick proto-CBL domain edge produces a low-pressure zone to the right of the thick domain. It drives the rheologically weak, high-temperature chemical lithosphere root toward the thin proto-CBL region, thus forming a continuous channel flow over a few thousand kilometers distance. The insert (white box on the left) shows a typical Poiseuille flow for stress exponent n = 1 (blue), 3 (orange), and 5 (red). A higher n produces a less deformed central part of the channel. (b) In the thick proto-CBL domain, the channel flow forms one strong shear zone (labeled s1) within the chemical lithosphere at a depth of ∼100 km and another (labeled s2) at the top of the sublithospheric mantle. (c) The driving pressure gradient in the channel decreases in magnitude with time, reducing the intensity of the channel flow. (d) After ∼2.5 Gyr of model evolution, TLAB has deepened to ∼220 km depth; shear zone s1 is preserved as frozen-in strain at mid-lithosphere depths, while the sublithospheric shear zone s2 is destroyed by the underlying mantle convection. The temporal evolution of the channel flow is further illustrated by depth profiles of material velocity, strain rate, effective viscosity, and temperature (Figures 4a-4d) at the thick proto-CBL part of the model (at distance = 1, 500 km). The top of the channel flow (s1) deepens at a rate of 1.5 km/Myr from its initial depth of 90 km depth (TLAB) at 0.4 Myr to a depth of 120 km at time = 19.8 Myr (Figures 4a-4c). While strain rate in CBL decreases by two orders of magnitude, CBL viscosity increases by one order of magnitude in ca. 20 Myr (Figure 4c). These changes are caused by a reduction of the pressure gradient in the channel (Figures 3a-3c, top panels), despite a ∼40°C temperature increase in the channel (Figure 4d) which would otherwise decrease the viscosity.

Testing Parameter Space for the Channel Flow Formation
We demonstrate a broad stability of the channel flow mechanism by testing a series of model setup parameters (see Table S2 in Supporting Information S1 for overview). It is outside the scope of the present study to investigate the full parameter space for channel flow formation, since our goal is to demonstrate that this new mechanism for the MLD origin is plausible.

Difference Between CBL Thickness
Continent-scale channel flow forms in the lower part of Archean CBL by a horizontal pressure gradient at the craton margin caused by a downwelling mantle flow near the thick-to-thin proto-CBL transition. Mantle downwelling is triggered by a density inversion in the mantle (Figure 2b) which leads to differences in intensity of mantle flow below the thick and thin proto-CBL domains. Therefore, the initial contrast in the CLAB depth in two domains may play a critical role in the model evolution. We test model results for initial CLAB depth in the thin proto-CBL domain of 210 , 170, and 150 km, and 95 km depth (95 km depth in the end-member model), while keeping the CLAB of the thick proto-CBL fixed at 235 km depth as before ( Figure 5).
The thermal expansion coefficient of mantle material affects both the buoyancy of the continental root and the vigor of mantle convection. Laboratory measurements on forsterite show that the thermal expansion coefficient reduces from 2.4 × 10 −5 to 1.3 × 10 −5 K −1 for a pressure reduction from 1 GPa (top of the upper mantle) to 18 GPa (bottom of the upper mantle) (Chopelas, 1990). The value used in the end-member model (1.0 × 10 −5 K −1 ) better represents the deep part of the mantle. With all other parameters fixed as in the end-member model, numerical experiments with 3.0 × 10 −5 K −1 (Figure 5a) produce channel flow like Figure 3b. A thermal expansion coefficient of 3.0 × 10 −5 K −1 is used in all other the studies of CBL thickness differences across the transition zone in Figure 5.
Decreasing the CBL thickness contrast across the transition zone from 140 km (Figure 5a) to 25 km (Figure 5d) reduces the pressure difference in the channel from 20 to 5 MPa, thus reducing the mantle material flow motion rate by one order of magnitude (Figures 5a-5d). Although the mantle material motion rate is deceased due to a smaller CBL thickness contrast than that in the end-member model, a continuous channel flow at depths of ∼100-200 km still develops in the CBL. After ∼20 Ma, material in the s1 shear zone (the channel top) of the model with a 65 km contrast in CBL thickness across the transition zone ( Figure 5c) generates >10 km dislocation at x = 1,500 km, which is about 1,000 km away from the convection center at x = 2,500 km (see circles deformation at x = 1,500 km in Figure 5c). However, in the model with a chemical lithosphere thickness of 210 km in the thin domain, there is almost no shear zone developed at x = 1,500 km ( Figure 5d). If the viscosity in the sublithospheric mantle is reduced by a factor of 5, with all other parameters unchanged, the model in Figure 5d with a thickness contrast of 25 km produces a continuous the shear zone s1 at a depth of ∼100 km as well due to lower viscosity and higher Rayleigh number in the sublithospheric mantle, where there is a more rigorous mantle convection (Figures 6a, 6b versus Figure 5d).

TLAB Depth and Density Contrast
We test how the initial TLAB depth (and therefore, initial cratonic geotherm, Figure 2c) affects the model evolution. Changing the initial TLAB depth from 90 km in the end-member model to 150 km does not prevent formation of the channel flow ( Figure S3c in Supporting Information S1). However, in this case the top shear zone (s1) forms at ∼150 km, that is, deeper than in the reference model. If the TLAB depth is too shallow, that means there are higher temperatures in the CBL than the end-member model, we do not observe the channel flow, but instead convection within the continent root. Whether or not the convection in the continent root interacts or merges with the sublithospheric mantle convection, depends on the density contrast between CBL and the mantle below and the Rayleigh number of the sublithospheric mantle Sleep, 2003b) (see more discussions in next sections). If the density contrast in the end-member model is reduced to 5 kg/m 3 , the continent root will be destroyed, and no channel flow is observed.

Width of the Thick-To-Thin Proto-CBL Transition
The width of the thick-to-thin proto-CBL transition does not affect the maximum pressure drop at the craton edge ( Figure 7) and therefore it does not control the conditions for the formation of a channel flow in the CBL. However, the width of the transition affects the style of mantle flow pattern. We illustrate this by a model with a wide CLAB vertical step of 145 km over 1,800 km distance (at x = 800-2,600 km) as compared to the end-member model with a narrow 145 km CLAB step over a 200 km distance at x = 2,400-2,600 km (Figure 3). A very broad transition restricts most of the deformation to the transitional zone (Figures 7a and 7b). After 6.0 Myr (Figure 7c), a large-scale convection cell near the right end of the transition (x = 2,300 km) starts splitting into two smaller cells with the opposite propagating directions, while the channel flow continues.

Mantle Melting
Our model does not consider melting effects, although melt extraction may consume latent heat and thereby lower temperature in the lithospheric root. We test this in two experiments ( Figure S2 in Supporting Information S1) with wet and dry solidus. The compensation of latent heat in a melting system is a complex process. If partial melts and solid grains are in equilibrium, that is, the temperature in the melting system is equal to dry (or wet) solidus of olivine, the continental lithospheric root can still accommodate large-scale channel flow at a temperature lower than in the end-member model ( Figure S2 in Supporting Information S1).

Delayed Formation of Channel Flow in the Archean
Model geometry and variations in mantle physical properties control a complex flow pattern, which may be different from those discussed above (Figures 3, 5-7 and Figures S2, S3 in Supporting Information S1). In particular, the formation of the channel flow depends on the initial conditions, which are poorly constrained for the Archean Earth. Among many tested models, we found a few cases, when a large-scale channel flow does not form at the early stages of model evolution, but develops after tens of million years. The model, named as M_den15_exp2, with the delayed channel flow formation has the same parameters as the end-member model (Tables S1 and S2  Figure 3. The model has a thermal expansion coefficient of 3 × 10 −5 K −1 , which is triple value of that in the end-member model. Increasing thermal expansion does not affect the channel flow pattern (a). Decreasing the thickness contrast across the transition zone reduces the channel flow motion rate due to lower pressure gradient in the channel.
in Supporting Information S1 and Figures 2 and 3), except that the density contrast across the CLAB is reduced from 30 to 15 kg/m 3 , and the thermal expansion coefficient of 2.0 × 10 −5 K −1 is double of that in the end-member model. Due to the higher value of thermal expansion and higher Rayleigh number (see Sections 3.2.6), the evolution of this model (Figure 8a) creates a very vigorous edge convection at x = 2,000-3,000 km.
The initial channel flow formed during the first 1.0 Myr is <500 km long and most of the deformation localized at the transitional part (Figures 8a and 8b). The leftward movement of this convection cell from the thick proto-CBL domain edge toward its interior deforms the shape of the CLAB. Vigorous convection below the thick proto-CBL domain (x = 700-1,200 km; Figure 8c) increases pressure in the directly overlying channel segment (Figure 8c, top). The rising pressure creates a large-scale channel flow in the chemical (depleted) lithospheric mantle to the right from the convection cell center. After 48.5 Myr of model evolution, the channel flow extends to the whole model length with a pressure contrast of 15 MPa between the maximum and minimum values across the model domains (Figure 8d). The pressure bulge above the convection cell center in the sublithospheric mantle does not form in the end-member model (Figures 3a-3c), which instead shows a monotonic pressure decay from the interior of the thick proto-CBL domain (model left) to its edge (model right). This difference is caused by the more vigorous convection in model M_den15_exp2 (Figure 8), which causes deformation of the CLAB geometry during the leftward propagation of the sub-CLAB convecting cells. If the compositional density contrast across the CLAB is further reduced, for example, by metasomatism of depleted cratonic mantle, the cratonic root will be destroyed by mantle convection, as reported in numerous studies of craton stability and longevity (Cottrell et al., 2004;Lenardic et al., 2003).

Mantle Rayleigh Number, Convection Vigor, and Channel Flow
Formation of a depleted mantle may be a progressive process, which may delay the channel flow formation. Early stages of the lithospheric mantle depletion may have led to an initially small contrast in density and lithospheric mantle thickness between cratonic and noncratonic domains. These small contrasts result in a very slow channel flow (Figures 6 and 8). Importantly, an experiment with a 25 km difference in CBL thickness demonstrates that even small perturbations at the CLAB can trigger channel flows in cratonic roots during or after the proto-craton formation. While mantle heterogeneities (e.g., mantle fluids) tend to perturb the CLAB shape, the resultant channel flow in cratonic roots tends to flatten it. Figure 6. Snapshots of the model evolution after 3.7 Myr (a) and 28.6 Myr (b). Compared to the model in Figure 5d (Figure 3), the viscosity in the sublithospheric mantle is reduced by a factor of 5. Notations as in Figure 3. The model develops a channel flow in the chemical lithosphere root, and the flow is more rigorous than that in Figure 5d, though they have the same geometry in the initial setup.
Model M_den15_exp2 (Figure 8) demonstrates that the early stage evolution, which is strongly controlled by the initial conditions, does not define whether or not channel flow may form at later stages. Therefore, we track changes in model key physical properties through time, especially at stages when convection cells begin to deform the CLAB geometry from below and/or above. We adopt a standard convection modeling approach to describe model evolution by dimensionless numbers, which incorporate most critical model parameters, given their large number other approaches are not justified.

The dimensionless Rayleigh number = ∆ 3
and lithosphere Buoyancy number = ∆ ∆ characterize mantle flow patterns. Here is thickness of the continent root, its chemical density, ∆ the chemical density contrast across the CLAB, ∆ the maximum temperature difference within the continent root, thermal diffusivity, thermal expansion coefficient, and average viscosity in the cratonic root.
We calculate Ra and B for a set of models with various ∆ , and , preserving the same geometry as in the end-member model (Section 2.2). As ∆ and change both in space and time, we limit observations to the area within the thick proto-CBL domain interior. For example, at 1.0 Myr of the M_den15_exp2 model evolution (Figure 8a), the mantle below its interior (x = 0-700 km) remains stable, while below its edge (x > 2,000 km) the mantle is involved in strong convection. Since averaging over the whole cratonic lithosphere root hides local heterogeneities, we only average the properties over ∼500 km distance within the thick proto-CBL interior.
The results fall into two major convection domains, the stable domain and the whole-layer mixed convection domain (Figure 9). In the stable domain, the thick continental root is neither affected by convection from below nor develops small-scale convection inside the continent root, which allow channel flow to develop. The other end-member case, with mantle convection in the layer between the CLAB and 660 km depth, includes delamination of the continent root into the hot, deep mantle (Figure 9b). The boundary between the two domains is defined by a critical Rayleigh number (Ra c ) (dashed line in Figure 9a). The Ra c linearly increases from ∼200 to ∼1,000 (with some outliers) for the Buoyancy number between 0.0 and 0.4-0.5 and is almost constant (Ra c ∼ 1,000) for B > 0.5. The latter value is close to Ra c ∼ 1,100 as reported for a typical Rayleigh-Benard convection with one rigid and one free slip boundary (Pellew & Southwell, 1940) and in comprehensive theoretical instability studies . Similar value of Ra c ≈1, 100 was reported in laboratory studies of convection instability for  (Figure 2), the model has the half density contrast across the chemical LAB (CLAB) (15 kg/m 3 ) and double value of thermal expansion coefficient of 2 × 10 −5 K −1 (see details in Tables S1 and S2 in Supporting Information S1). Notations as in Figure 3. At the early stage (a, b), the model does not develop a channel flow in the chemical lithosphere root; a channel flow forms at later stages (c, d) due to a pressure increase directly above the convection cell centered around x = ∼ 800-1,000 km. two-layered fluids; this critical value is almost independent on the viscosity change across the material interface (Cottrell et al., 2004), that corresponds to the CLAB in our study.
The exact position of the critical Ra c -B curve has some uncertainty (Figure 9a). Close to the critical boundary, the two mantle layers (depleted lithospheric and fertile asthenospheric) are stratified, but the CLAB interface oscillates (Figure 9c). The convective overturning instability essentially drives the dynamics of the model (Figure 2b), while viscous forces and thermal diffusion stabilizes the density contrast across the CLAB and work as a restoring force. The efficiency of heat diffusion in the depleted cratonic root affects the mantle flow pattern. Large-scale channel flow in the cratonic root is less efficient for heat diffusion than mantle convection, which leads to accumulation of extra heat in the cratonic root when radioactive heat generation rate is high (Figure 4d). It increases the Rayleigh number and causes oscillation of the CLAB interface (Figure 9c).  Tables S1 and S2 in Supporting Information S1 for details) is the same as the end-member model except for a smaller density contrast across the chemical LAB (CLAB) (5 kg/m 3 ) and triple value of thermal expansion coefficient (3 × 10 −5 K −1 ). (c) Snapshot of model M_den15_ exp2 evolution at time = 294.7 Myr, which shows a stratified oscillatory CLAB interface. Dark blue-chemical boundary layer (depleted lithosphere), light bluesublithospheric mantle. White arrows-velocity of material motion with the scale atop the snapshots.

Discussion
Our numerical model provides a plausible dynamic process for the formation of MLD by channel flow within the cratonic lithosphere during the Archean. Channel flow typically forms below 80-120 km depth (deeper at later times due to cooling of the Earth and decay of the driving pressure gradient) with the bottom at depths >200 km. The top (s1) and bottom (s2) of the channel flow experience maximum shear deformation, whereas its middle part (120-200 km) is less deformed (Figures 3-8). Therefore, the main signature of channel flow is frozen-in deformation at the channel top. It will manifest itself in a tens of km thick layer by changes in for example, grain size, seismic velocities, seismic and electrical anisotropy (Figures 1a and 1b). The bottom of the channel flow is likely to be destroyed by mantle convection. It may be preserved in rare cases of an exceptionally deep depleted cratonic roots and will manifest itself similar to the top shear zone.

Channel Flow and Seismic Low-Velocity Zones
Intensive shearing at the top of the channel flow (layer s1, depth 90-120 km, Figures 3 and 7) corresponds to a low-viscosity zone (Figure 4c) at depths consistent with the depth to seismic velocity discontinuities (Rader et al., 2015;Thybo & Perchuć, 1997) (Figure 1a). Although a low viscosity does not necessarily require a low seismic velocity (Karato, 1995), the seismic low-velocity zone is generally regarded as a rheologically weak upper mantle zone (Thybo, 2006). Mineral physics experiments on melt-free polycrystalline aggregates of Fo90 olivine indicate that shear modulus decreases with reducing grain size, for example, at 1,300°C, the shear modulus is decreased by 3-4% if the grain size is reduced from 10 to 1 mm at a period of 10-100 s (Jackson & Faul, 2010). Grain size reduction my be caused by intensive shearing in the mantle (Bystricky et al., 2000), as evidenced by mantle xenoliths from the Kaapvaal craton with a drop in olivine grain size from 12 to 5 mm at 130-160 km depth (Mercier, 1980). However, at 1,100°C, the shear modulus is decreased only by ∼2%, corresponding to ∼1% reduction in shear-wave velocity, for the grain size reduced from 10 to 1 mm (Jackson & Faul, 2010). The secular cooling of the earth may make the seismic velocity discontinuity more difficult to be detected at present than at the Archean, if the grains size in the shear zone has not changed significantly since the formation of cratons (Speciale et al., 2020). Therefore, the grain size reduction caused by channel flow alone is difficult to explain current observations.

Channel Flow, Phase Transitions, and Fluid Migration
Volatile-rich melts percolating through the lower lithosphere may concentrate at some depths, leading to accumulations of minerals with low seismic velocities (pyroxenes, phlogopite, amphibole, carbonates) (Aulbach et al., 2017). The presence of these minerals at mid-lithosphere depths was proposed as an explanation of MLD observed in seismic S-RFs (Selway et al., 2015). However, our compilation shows little correlation between the depth of seismic MLD and the petrologically constrained top of the metasomatically reworked lithosphere mantle in nearly all cratons (Figure 1a). The top of the melt metasomatism is everywhere deeper than the geophysical observations of MLD (Figure 1a). Instead, seismic observations of MLDs are likely to be related to spinel-garnet transition (Figure 1a).
Both shear zones (s1 at 80-120 km and s2 at 200 km) form below low-deformed, high-viscosity lithospheric mantle. Above the channel top (s1) high viscosity is due to low temperature; above s2 it is due to high strength from low strain rate for dislocation creep in the central part of the channel (Figure 3a′). Volatile-rich melts from asthenosphere (>200 km) may be blocked by the central part of the channel flow, while at shallower depth (60-140 km) volatile fluids may be released by chemical reactions (Scambelluri et al., 2006;Skirrow et al., 2018) and can be capped by high-viscosity lithosphere above the top shear zone (s1). The top shear zone (s1) at ∼100 km coincides with the depth of pargasite dehydration solidus, 1,100°C at ∼90-100 km for old cratons, and the fluids rising from underlying asthenosphere may accumulate at the base of the pargasite stable zone at a depth around 100 km, preventing the rest of the mantle at a shallower depth than 100 km from rehydration and chemical alteration (Kovács et al., 2021). Both mechanisms may work together to create a barrier to rising fluids and accumulate fluids at ∼100 km.
The thickness of the less deformed middle part of the channel is comparable to the thickness of intra-lithosphere anomalous layers in cratons from seismic, MT and xenolith data (Figure 1b), which varies from 20 to 70 km with a global average of ∼50 km. Our model indicates that when conditions are suitable to produce channel flow (i.e., sufficient pressure gradient), the depth of shear zone s1 depends on mantle temperature (Figure 3 and Figure S3c in Supporting Information S1). If the craton is later reheated, there could be multiple s1 shear zones.

Channel Flow, Seismic Anisotropy, and Contrasting Seismic Signals From MLD and LAB
Channel flow causes strong radial anisotropy at the top (s1) and base (s2) of the channel, but not in-between, where strain rate is very low as typical of Poiseuille flow (Figure 3a′). This is consistent with observations of stronger radial anisotropy in the shallow mid-lithosphere mantle than at around the LAB in the cratons in western Australia (Sun et al., 2018). Radial anisotropy may produce a seismic MLD as observed in S-RF (Selway et al., 2015). Dislocation creep in non-Newtonian material with the stress exponent n > 1 (Figure 3a′) would enhance a localized shear zone s1 and associated seismic anisotropy (Karato, 1995). Our channel flow model predicts two distinct motion patterns for mantle material in the young and present Earth (Figures 3a-3c and 9). In the young Earth, a continuous channel flow forms beneath the TLAB at a depth of ∼100 km, while after billions of years of mantle cooling, multiple small convection cells underlie the TLAB of present earth at a depth of ∼200 km (Figures 3d and 9). Therefore, in the Archean a nearly horizontal shear below the TLAB due to the channel, while the mantle flow below the present TLAB is characterized by alternating upwelling and downwelling (Figure 3d). Higher horizontal strain and less vertical shear developed in the channel top (s1) than the channel bottom (s2) should produce a stronger radial anisotropy (V SH > V SV ) in the top than the bottom.

Implications for Mantle Dynamics
Channel flow cannot form in modern cratonic settings: cold and rigid cratons with high-viscosity lithosphere roots tend to resist the development of a channel flow, even in cases of significant horizontal pressure gradient, for example, from edge convection (King & Ritsema, 2000). We note that during the evolution of the end-member model, at ca. 20 Myr, geotherms develop an intriguing curvature at ∼200 km depth (Figure 4d), when temperature in the CBL is higher than in sublithospheric mantle. Such kinked, hot lithospheric geotherm may exist during a period of very high rate of radioactive heat generation in the Precambrian lithosphere (Michaut & Jaupart, 2007). High temperature in the cratonic root reduces its viscosity and facilitates formation of channel flow once a sufficient pressure gradient develops. In our numerical model, we assume that this thermal condition is satisfied in the Archean with mantle temperatures ca. 200°C higher than at present. However, large-scale mantle thermal anomalies may introduce additional heat at the base of cratons which may increase temperature at the lithosphere base by 100-150°C, possibly over distances of 500 km away from the impact center (Sleep, 2003a). If the increased temperature at the continental lithospheric root reaches the wet solidus, channel flow may be triggered in the lower cratonic lithosphere too ( Figure S2 in Supporting Information S1). This may contribute to explaining some enigmatic features, such as high (2-7%) shear-wave velocity anomalies which appear to continue from the ancient continental cratonic roots of the southern African cratons to beneath the oceanic basins 500-1,300 km away from the cratonic margin (Grand, 2002). If the mantle thermal anomaly associated with the separation of the South American and African continents increased the temperature of the cratonic root by 100-150°C, it would have reduced the mantle viscosity by approximately two orders of magnitude. Mantle edge convection tends to create a low-pressure zone in the upper mantle at cratonic margins with respect to cratonic interior, which may cause a weakened cratonic root to flow toward the ocean like the pattern observed in the reference model (Figures 3a-3c). If a continental plate moves in the direction opposite to the channel flow, it may enhance shearing at depths of the top of the channel flow (s1), and this process may explain observations of high seismic velocity anomalies such as cratonic continent roots beneath the oceans around southern Africa (Wang et al., 2017).
As reviewed in Section 1, seismic, MT, and xenolith data report variable depths to the top of the intra-lithosphere discontinuities in cratons (Figure 1). Our numerical experiments generate the top shear zone (s1) in a relatively narrow depth range (100-130 km), which fits observations of seismic anisotropy and the seismic low-velocity zone from controlled-source seismic experiments better than other intra-lithospheric geophysical features. Our model provides a generic framework for understanding the origin of seismic radial anisotropy; it also provides a physical mechanism to explain the major MLD features observed from S-RF (Selway et al., 2015). It is possible that, during later cratonic evolution, various local tectonic events can significantly modify preexisting frozen-in shear zones. For example, a high-resolution teleseismic S wave study imaged two MLDs around the Midcontinent Rift and the southern Appalachians, that is, in the regions that have experienced significant tectonomagmatic modifications since cratonization (Liu & Shearer, 2021). Phanerozoic rifting processes in cratons can raise temperature in the continental root and create a CBL thickness contrast, thus supporting the formation of channel flow in the continent root as well (Huismans & Beaumont, 2011). The reheating event may cause deeper and smaller scale shear zones than the original channel flow and produce multi MLDs, but generally over much shorter distances than the originally formed MLD.

Conclusions
We propose that the MLD was formed by channel flow in the Archean Earth and test this hypothesis by two-dimensional thermomechanical numerical modeling. The critical aspects of our model are that at time zero, the thermal LAB (TLAB) is shallower than the chemical LAB (CLAB), and the CLAB transition between thick and thin proto-lithosphere triggers mantle flow in the lithospheric (above the CLAB) and sublithospheric mantle ( Figure 10).
1. High mantle temperature in the Archean Earth is critical for the formation of the channel flow in a thick buoyant cratonic lithosphere. Continent-scale channel flow forms in the lower chemical lithosphere due to a horizontal pressure gradient originating from mantle convective instabilities around the craton edge. 2. The channel flow has a Poiseuille-type deformation style and creates deformation zones at 80-120 km and 200-250 km depths. The top deformation zone freezes-in at mid-lithosphere depths, while the bottom zone may be destroyed by mantle convection. The less deformed central part of the channel may become a barrier to fluids rising from the sublithospheric mantle. 3. The mid-lithosphere channel flow predicts and explains the presence of seismic low-velocity zones and seismic MLD in the cratons. Intensive shearing near the top of the channel flow may reduce olivine grain size, promote accumulation of volatile-rich melts and create frozen-in anisotropy to develop a characteristic seismic discontinuity. 4. The interplay between the mantle Rayleigh number and the lithosphere chemical buoyancy controls the formation of either a mid-lithosphere channel flow or mantle convection cells with possible destruction of cratonic roots. 5. Secular cooling affects the flow structure beneath the TLAB, which switches from the early mode with a channel flow to the formation of small multiple convection cells at later stages. These multiple convection cells with alternating upwelling and downwelling may reduce the strength of radial anisotropy caused by large-scale channel flow in shallow proto-cratonic lithosphere.

Data Availability Statement
This study is based on published data, which are available as referenced in the article. The numerical open-source code Underworld2 used for the calculations is available at https://doi.org/10.5281/zenodo.1436039.