Convective Self-Compression of Cratons and the Stabilization of Old Lithosphere

Despite being exposed to convective stresses for much of the Earth's history, cratonic roots appear capable of resisting mantle shearing. This tectonic stability can be attributed to the neutral density and higher strength of cratons. However, the excess thickness of cratons and their higher viscosity amplify coupling to underlying mantle flow, which could be destabilizing. To investigate the stresses that a convecting mantle exerts on cratons that are both strong and thick, we developed instantaneous global spherical numerical models that incorporate present-day geoemetry of cratons within active mantle flow. Our results show that mantle flow is diverted downward beneath thick and viscous cratonic roots, giving rise to a ring of elevated and inwardly-convergent tractions along a craton's periphery. These tractions induce regional compressive stress regimes within cratonic interiors. Such compression could serve to stabilize older continental lithosphere against mantle shearing, thus adding an additional factor that promotes cratonic longevity.

2 of 10 found a similar amplification of tractions, but also showed that the strain-rates at the cratonic base diminish as lithospheric roots get thicker. This inverse relation between tractions and the strain-rates may slow the deformation of a cratonic root, and therefore might be an important factor for the long-term survival of cratons. Cooper and Conrad (2009) attributed elevated tractions at the base of cratons to greater coupling to mantle flow, which has been noted in models with thick cratonic roots (Becker, 2006;Zhong, 2001). However, more recent models, especially those employing free-slip surface boundary conditions that more closely resemble Earth's own conditions, show that tractions are primarily amplified along the periphery of cratons ( Figure  3 from Paul et al., 2019). Although Paul et al. (2019) speculated that cratonic edges might more effectively absorb mantle stresses compared to cratonic interiors, a proper quantitative analysis of such a phenomenon is lacking.
Here, we explore the origin of higher tractions along craton boundaries and consider their implications for the stability of cratons. We build instantaneous global models of mantle convection and examine how mantle flow is modified due to the presence of thick and viscous cratons. We hypothesize that the diversion of mantle flow by the thick and highly viscous root of a craton can generate strong and inwardly-convergent tractions at the craton's periphery. We test our hypothesis using various models with different viscosity combinations for cratons and asthenosphere. We consider how large convergent tractions, which are generated by the cratons themselves, may support cratonic stability against mantle shearing, and therefore could be essential for cratonic longevity.

Mantle Convection Models
We use the finite element code CitcomS to develop instantaneous spherical models of mantle convection (Zhong et al., 2000). The code assumes the mantle to be a viscous and incompressible fluid. It solves the conservation of mass, momentum, and energy equations with the Boussinesq approximation and infinite Prandtl number. The smallest resolution of our models in the horizontal direction is ∼0.7° × 0.7°. The vertical resolution in the top 300 km is 24 km, and from 300 km to the core-mantle boundary (CMB) it is ∼50 km. Mantle flow is driven by the density anomalies obtained from SMEAN2 seismic tomography (Jackson et al., 2017), which is a combination of S40RTS (Ritsema et al., 2011), GyPSuM-S (Simmons et al., 2010) and SAVANI (Auer et al., 2014). Following earlier, similar efforts (Becker, 2006;Paul & Ghosh, 2020), a velocity-density scaling value of 0.25 is used to convert velocity anomalies into density anomalies. Higher velocity regions under the continents were removed down to 300 km to impose neutrally buoyant cratons. We keep a free-slip boundary condition at the surface and at the CMB. Reference viscosity, Rayleigh number, thermal expansivity and thermal diffusivity values are kept at η ref = 10 21 Pa.s, Ra = 4 × 10 8 (considering Earth radius as the length scale), α = 3 × 10 −5 K −1 , and κ = 10 −6 m 2 /s, respectively.
In our models, the mantle is divided into four layers based on their relative viscosity with respect to the upper mantle (300-600 km) reference viscosity (η ref = 10 21 Pa.s). The top 100 km is assigned as the lithosphere with a radial viscosity of 30 × η ref (30 × 10 21 Pa.s). The radial viscosity of the asthenosphere (100-300 km) is varied between 0.01 (10 19 Pa-s), 0.1 (10 20 Pa-s) and 1 (10 21 Pa-s) times the reference upper mantle viscosity. The radial viscosity of the lower mantle (660-2,900 km) is made 50× larger than the reference viscosity (50 × 10 21 Pa.s). On top of this radially-varying viscosity structure, we impose lateral viscosity variations. In the top 300 km, we approximate temperature-dependent viscosity using a linearized Arrhenius law η = η R × exp(E(T 0 − T)), where η R is the radial viscosity of any layer, T 0 is the non-dimensionalized reference temperature, and T is the non-dimensionalized actual temperature, where the maximum temperature corresponds to 1,300°C. E is a dimensionless quantity that controls the strength of the temperature dependence. We have tested several models to find suitable values for E (cf. Paul et al., 2019) and use a value of 5, which produces 10× weak plate margins compared to the continental interiors. Weak plate margins originate due to slow velocity anomalies inherent within the SMEAN2 tomography model. Stronger continental interiors with weaker plate margins enhance plateness and produce plate velocities comparable to observations ( Figure S1 of Paul et al., 2019). We also incorporate high viscosity cratons in our models, where the locations of cratons are taken from the 3SMAC model (Nataf & Ricard, 1996). Cratons are made 10×, 100×, and 1,000× more viscous than the surrounding lithosphere, making their actual viscosities 30 × 10 22 Pa.s, 30 × 10 23 Pa.s and 30 × 10 24 Pa.s, respectively. Cratons have uniformly viscous keels up to a depth of 300 km. Our reference models omit cratons and only incorporate temperature-dependent viscosity to create lateral viscosity variations.
PAUL ET AL.

Tractions Within Cratons
We analyze the rϕ and rθ components of stress tensor (σ ij ; i,j = r: radial component, ϕ: longitudinal component, θ: co-latitudinal component) from model outputs and calculate traction vectors ( ⃗0 ) from the reference model ( Figure 1). In the reference model, the magnitudes of traction vectors are less than ∼5 MPa, and their orientations are guided by density anomalies within the model (Figure 1a). Incorporating 100× viscous cratons in the same model significantly affects traction vectors ( ⃗ ) along the edges of cratons (Figure 1b). A few enlarged maps near the cratonic regions show this effect more prominently (Figures 1c-1i). Most cratons show rings of high traction magnitude along their periphery, where traction directions become inwardly convergent. Elevated inwardly convergent tractions appear prominently along the western margin of the North and South American cratons (Figures 1c and 1d), the eastern, western and southern margins of the Siberian and Australian cratons (Figures 1e and 1f), the northern and southern margins of the Scandinavian craton (Figure 1g), and the eastern margins of African cratons ( Figure 1h). The Indian craton, being very small in size, experiences convergent tractions all around its periphery ( Figure 1i). The southernmost part of the African craton shows an outwardly directed traction, which is the only exception ( Figure 1h). We have tested a model with high lithospheric viscosity (150×) and similarly found large traction ratios along cratons' periphery ( Figure S1 in Supporting Information S1).
To quantify the increase in traction magnitudes caused by the presence of cratons, we normalize the traction magnitudes from models with cratons using those from the reference model ( |⃗|∕| ⃗0| ) . In the presence of cratons that are 100× more viscous than the rest of the lithosphere, the maximum traction ratio increases by up to 80-100 times at ∼120 km depth along the edges of cratons (Figures 1c-1i). The magnitude of the traction ratio along the craton edges can be influenced by the viscosity structure imposed in our models ( Figure 2). To investigate the dependence of the traction ratio on viscosity structure and depth, we calculate the average traction ratio at various depths along the edges of cratons ( Figure 2a). The edges of cratons are identified by regions with traction ratio ( |⃗|∕| ⃗0| ) more than 5 at 120 km depth (Figures 1c-1i). The general trend shows that the average traction ratio varies between 10 and 15 within the top 100 km of craton edges (Figure 2a), which are proximal to viscous non-cratonic lithosphere. The average traction ratios increase with depth, reaching peak values in the mid-cratonic depth range of ∼160 km. The highest traction ratio occurs near the depth of peak horizontal velocity, which occurs in the mid-asthenosphere. With increasing depth, the traction ratio gradually falls before reaching another smaller peak near the base of cratons at ∼270 km depth ( Figure 2a). The magnitude of the traction ratio depends on the combination of the craton and asthenosphere viscosity. Higher viscosity contrast between a craton and its surroundings can enhance traction ratio. Indeed, highly viscous (1,000×) cratons exhibit the largest traction ratios, which exhibit peak values of 35-45 for mid-asthenospheric depths. Models with smaller viscosity contrasts (e.g., stronger asthenosphere with relative viscosity 1×) exhibit relatively smaller traction ratios near the craton edges.
We use centroid moment tensor (CMT) type symbols ( Figure 1) to quantify the state of stress within cratons due to inwardly convergent tractions. CMT symbols are colored by the ratio of mean horizontal stress ( h = 1 2 ( + ) ) and the second invariant of the deviatoric stress . A negative ratio (σ h /σ II < 0) indicates a compressive stress regime and vice-versa. The deformation states shown imply that the model without cratons (Figure 1a) has compression only along the convergent plate boundaries, that is, along the margins of the Pacific and the Indo-Eurasia collision zones. The same model with 100× viscous cratons (Figures 1b-1i) acquires a compressive stress regime within all cratons, except in South Africa (near the Kalahari and Kaapvaal cratons) (Figure 1h). This compressive nature is consistent throughout the cratonic root at greater depths ( Figure S2 in Supporting Information S1).

Origin of Compression Along Craton Edges
Our models demonstrate an amplification of tractions ( ⃗ ) along craton edges that induce a highly compressive state within viscous cratons. To understand the origin of this regional compressive stress regime, we calculate the traction vector ( ⃗ ) from the σ rϕ and σ rθ components of the deviatoric stress tensor that relate to horizontal shear, 10.1029/2022GL101842 5 of 10 where v ϕ , v θ , and v r are the horizontal and vertical components of the velocity vector and η is the viscosity. If these shear components dominate the stress tensor, then the magnitude of horizontal traction is given by The presence of a thick and highly viscous craton obstructs horizontal asthenospheric flow, and deflects it downward near the craton edges. Such velocity diversion can make the v r and v r components small near the craton edges, implying that the first terms in Equations 1 and 2 can be neglected. With stronger downward diversion, the vertical velocity (v r ) increases approaching a craton edge. Thus, the horizontal gradient of the vertical velocity component, that is, the second term in Equations 1 and 2, becomes the controlling factor for the origin of high tractions along craton boundaries. A small change in velocity gradients near cratons can thus induce higher tractions around them as tractions originate from velocity gradients multiplied by the high viscosity of cratons (Equation 3).
Large horizontal gradients of vertical velocities, induced by viscosity heterogeneity associated with cratons, thus amplify tractions. We calculate such gradients as:  6 of 10 To highlight the impact of thick cratons on the gradient, we compute the ratio of the horizontal velocity gradient (RHG) as: no_craton are the horizontal gradient of vertical velocities from models with and without cratons, respectively. RHG can quantify the concentration of downward flow due to the presence of viscous cratons, where RHG ≫ 1 indicates strong vertical velocity deflection.
Similar to the rings of high traction zones, we find rings of elevated RHG along the craton periphery (Figure 3). Elevated RHG values can be interpreted as horizontal velocities converting into vertical velocities near the craton boundary due to cratons' excess thickness and viscosity. Horizontal gradients of vertical velocity should amplify tractions, and indeed contours of traction ratio >5 typically lie next to regions of high RHG values (Figure 3). The reduction of horizontal velocity at craton edges is clearly visible underneath North and South America (Figures 3a and 3b). The strong velocity decrease arises because slabs underneath these two cratons force a rapid asthenospheric flow that is impeded by stiff cratons. The velocity gradient variations expressed by RHG are also controlled by the angle between the craton edge and the direction of horizontal flow. In our density-driven flow models, the mantle flows from west to east under the North American plate, remaining almost perpendicular to the western face of the craton (Figure 3a). Hence, the maximum velocity diversion, or highest RHG value, occurs along the western margin of the North American craton. On the contrary, the southeastern margin of the craton, being almost parallel to flow, shows no significant change in RHG values. This pattern resembles the change of traction vectors along the western and eastern margins of the North American craton (Figure 1d

of 10
To investigate how downwelling on the craton edges varies with depth and viscosity structure, we calculate variations of average RHG within the region where RGH value > 5 (Figure 2b). In the top 100 km, the average RHG varies within 15-17. In the mid-cratonic depth range (100-250 km), the average RHG value decreases to slightly less than 12.5. Deeper than 250 km depth, RHG increases again, reaching a peak near the base of cratons. These two peaks near the top and bottom of craton may appear due to the most significant change of velocity gradients occurring above and below the asthenosphere, giving rise to a "z" type velocity profile, considering left to right horizontal flow (e.g., Figure 4d).
We compare the velocity cross-sections from our models with and without cratons (Figure 4) to investigate the actual nature of flow diversion along craton edges. Downward mantle flow near craton edges has previously been attributed to lateral temperature variations (i.e., edge-driven convection, e.g., King and Ritsema, 2000), but our results suggest that such flow diversion is a natural consequence of global mantle convection operating in the presence of lithospheric viscosity heterogeneity. Cross-sections underneath the South American and the North American cratons show the most notable changes in velocity along their western margins (Figures 4a-4d). In both cases, the mantle flows from west to east in the absence of a craton due to density heterogeneity present in our models (Figures 4a and 4c). Convergent flow west of the South American craton occurs due to the subducting Nazca slab. In the presence of a thick and viscous craton, the convergent flow velocity is diverted along the craton margin and gets concentrated below it (Figure 4b). Similar velocity diversion is also visible around the western margin of the North American craton, where the flow gets diverted downwards and is concentrated below the craton (Figure 4d). Flow diversion by the Scandinavian craton occurs along a north-south orientation (Figures 4e  and 4f). However, the diversion is relatively weaker, most likely due to the absence of nearby mantle slabs to drive the flow in the model. Weaker velocity diversion is also reflected in less elevated traction magnitudes compared to the American cratons. The size of the Western African craton is significantly smaller than the rest, but the change of velocity field is considerably pronounced (Figures 4g and 4h). A downward flow along the eastern margin of the craton denotes the change in RHG (Figure 4h). The Australian craton also shows velocity diversion (Figures 4i and 4j) along an east west profile, leading to amplified tractions. The South African craton  (Figure 4j). In this scenario, the horizontal velocities get diminished due to the craton, and vertical upward velocities on the eastern side become stronger along the craton boundary. Therefore, the traction magnitudes increase along the South African craton's southeastern margin near the Kalahari and Kaapvaal cratons, as they do for the other cratons, but the tractions are outwardly directed and the stress regime becomes extensional (Figure 1h). Such extension could be a potential reason for recent thinning of the Kaapval craton (cf. Mather et al., 2011).

Role of Self-Compression in Craton Stabilization
Understanding cratonic survival has remained a long-standing problem in the geoscience community. Bedle et al. (2021) noted three significant geodynamical properties of a stable craton: (a) thick and buoyant cratonic roots, (b) highly viscous roots, and (c) integrated high yield strength that minimizes deformation. However, depending on their evolution, cratons can become unstable or partially destroyed (Bedle et al., 2021;Lee et al., 2011). For example, rapidly thickened lithosphere (e.g., Beall et al., 2018) can be subjected to basal erosion, subsequently leading to destabilization (Lenardic & Moresi, 1999). Thus, a self-driven and sustained process of gradual thickening may be essential to craton stabilization. Wang et al. (2018) has previously attributed such self-thickening to tectonic shortening stabilized by gradual gravitational thickening. However, they did not explore the nature of the stresses and tractions acting within the cratons, which may underpin slow and gradual thickening. We infer that such slow thickening could be controlled by self-compression within cratons and may be crucial for craton stabilization. Recently, a study suggested that the Slave craton may have regrown with time after being destroyed by the McKenzie plume . Self-compression could support such recratonization.
The shape of cratonic roots can influence the diversion of flow along the craton boundary, which subsequently deforms the craton interior Currie & van Wijk, 2016). Cooper et al. (2021) showed that a vertical craton margin can resist such deformation compared to margins that slope downward toward the craton interior. Our models consider roots with a sharp vertical viscosity contrast between the craton and the surrounding asthenosphere. In the future, it will be interesting to investigate the nature of flow diversion for cratons of different root geometries and more gradual viscosity contrasts with their surroundings. However, the horizontal length scale of the mantle flow diversion is on the order of 1,000 s of km (e.g., Figure 4). Hence, the sharpness of the viscosity contrast may have a smaller effect on cratonic self-compression compared to the actual magnitude of the lateral viscosity variations. Slow and continuous thickening induced by self-compression may also help to maintain steeper edges for cratonic roots, enhancing their stability.
Geologically, cratons are not individual single units; instead, they are composed of multiple protocratons that together form a larger continental mass (Bleeker, 2003). For example, the North American craton is composed of the Superior, Slave, Wyoming, Hearne, Rae, and several other small blocks (Canil, 2008); the Indian craton is assembled with five smaller units, Dharwar, Bastar, Singhbhum, Bundelkhand, and Aravalli (Pandey, 2020). Since their formation and amalgamation, larger continental units have remained together for more than a couple of billion years. There are some instances of delamination or partial destruction of cratons Menzies et al., 1993), but none of them were completely split apart. Self-compression could help to keep smaller continental blocks together within larger cratonic units. It also may be a key reason that older continental units did not split away during supercontinental break-up events. In the future, time-dependent numerical models should be developed to study the effect of self-compression in the craton stabilization process.

Conclusions
The diversion of mantle flow by thick and viscous cratonic lithosphere induces self-compression within the cratons themselves. Traction magnitudes increase along the craton periphery, and their directions become convergent toward craton interiors. Traction magnitudes depend on the viscosity structure of the craton, asthenosphere, and lithosphere. In the presence of a 100× viscous craton, traction magnitude increases to 15-20 MPa (Figures 1b-1i), more than an order of magnitude larger than cases without cratons (Figure 1a). The inward-directed orientation of tractions along the craton boundary appears to be a universal phenomenon (Figure 1b, Figure S1 in Supporting Information S1), except for the southernmost part of the African craton (Figure 1h). We infer that such convergent tractions originate from the diversion of (typically downward) mantle flow due to thick and viscous cratonic roots. We test our hypothesis by calculating the ratio of the horizontal gradient of vertical velocity (RHG, 9 of 10 Equation 5). Our calculations demonstrate that large velocity gradients along craton margins amplify tractions along the craton periphery. For most cratons the downward diversion of mantle flow produces inward-directed trac tions that induce a compressive stress regime within all cratons. The South African craton presents the only exception, where upwelling flow generates extension. We conclude that self-compression could be a key mechanism that drives the slow and gradual thickening of cratons, enhancing their stability. Such compression may also hold multiple smaller cratons together, merging them into larger blocks. Cratonic self-compression thus may be an essential stabilizing component that allows cratons to resist the destructive forces of mantle convection over billion years.

Data Availability Statement
The latest version of CitcomS code is freely available for download on GitHub (https://github.com/geodynamics/ citcoms). Example input files and the model output can be downloaded from JP's personal GitHub repository: https://jyotirmoyp.github.io/research/craton/ or https://doi.org/10.5281/zenodo.7264900. Detailed mathematical calculations and formulations are given in text which can be used to reproduce the results.