Geochemical observations suggest that cratonic xenoliths originate from depths no greater than 250 km, which implies a maximum craton thickness. What determines this thickness? In general, the stability and longevity of cratons depend on their ability to resist deforming forces exerted from below by the flowing, evolving mantle. We employ an analytical approach to relate the lithospheric viscosity structure of a craton to shear tractions exerted by asthenospheric shear flow. As the net thickness of the chemical and thermal components of a craton increases, these tractions increase dramatically for non-Newtonian rheology and begin to deform the cratonic base. Thus, overly thick lithosphere is subjected to large basal stresses that weaken the thermal boundary layer and diminish its ability to buffer the craton from destructive mantle flow. This feedback prevents cratonic lithosphere from growing thicker than a maximum value. However, we show that this maximum thickness increases slightly with increasing vigor of mantle convection but decreases rapidly as convective vigor decreases. Thus, we predict relative stability of cratonic thickness during most of Earth's history but instability in the future.

Cratons are regions of long-term stability on Earth's surface. This long-term stability has been attributed to the intrinsic chemical buoyancy, elevated viscosity, and finite strength of cratonic lithosphere (Jordan, 1978; Lenardic and Moresi, 1999; Sleep, 2003; King, 2005). Previous studies have used buoyancy arguments to define the maximum and minimum lithospheric thickness required for craton stability (Lenardic and Moresi, 1999; Cottrell et al., 2004; King, 2005; Cooper et al., 2006). If a craton is too thin, then it may not possess substantial positive chemical buoyancy to resist foundering within convective downwellings (Cooper et al., 2006). Similarly, an overly thick craton may also become unstable due to the increased role of downwelling convective instabilities (Cottrell et al., 2004; Conrad, 2000; Lenardic and Moresi, 1999). However, buoyancy alone cannot control the craton's ability to resist deformation over large geologic time scales because the mechanical strength of the craton serves to resist the deformation that may lead to thinning of the craton (Lenardic and Moresi, 1999). Furthermore, the stresses exerted on the base of a craton by the convecting mantle will interact with cratonic rheology to affect stability. Several authors have shown that deep continental roots increase the coupling of the mantle to the cratonic base (e.g., Zhong, 2001; Becker, 2006). Thus, because basal tractions increase with lithospheric thickness (Conrad and Lithgow-Bertelloni, 2006), they may place an upper bound on the maximum thickness that cratons can achieve.

Thickness estimates for cratonic lithosphere range from 150 to 300 km based on geochemical and seismic observations (Gung et al., 2003; Carlson et al., 2005; Lee et al., 2005). The large range of thicknesses is due in part to the discrepancy between the two methods. Seismic imaging tends to give larger estimates for craton thickness, while geochemical estimates of lithospheric thickness based on xenolith data rarely exceed 200 km (Gung et al., 2003; Carlson et al., 2005; Lee et al., 2005). The xenolith data could be marking the boundary of a chemically distinct lithosphere (or chemical boundary layer), while seismic imaging could be distinguishing a thermally distinct lithosphere (or thermal boundary layer). Regardless, if cratonic lithosphere is rigid, and thus unable to participate in convective heat transfer, the thermal boundary layer will reside between the conducting chemical boundary layer and the convecting asthenosphere to gradationally guide the transition between the two heat-transfer methods. Therefore, seismological estimates of cratonic lithosphere thickness, which are based on the thermal signature of the lithosphere, will include this transitional boundary layer in addition to the chemically distinct rigid lithosphere (Fig. 1).

This transitional layer also could be acting as a mechanical buffer between the convecting mantle and the rigid cratonic lithosphere. In this case, the transition layer may absorb some of the shear deformation that the convecting mantle would otherwise impart on to the rigid cratonic lithosphere, and, in this way, it contributes to the overall stability of the craton. However, Conrad and Lithgow-Bertelloni (2006), using a global mantle-flow model, observed that as lithospheric thickness increases, the shear tractions exerted by the mantle on the lithospheric base also increase. Thus, the thickest cratonic lithosphere experiences the greatest mantle traction stress yet remains resistant to deformation. The thickness estimates used in calculations by Conrad and Lithgow-Bertelloni (2006) are based on a seismically determined thermal lithospheric thickness (Gung et al., 2003; Ritsema et al., 2004) that includes both the rigid cratonic lithosphere and the transitional boundary layer. Therefore, the deformation regime they infer for the base of the thermal lithosphere may not reflect that which occurs at the base of the rigid cratonic lithosphere.

Numerical models that include both chemical and thermal lithospheric layers show that the thickness of the transitional layer depends on the thickness of the rigid cratonic lithosphere (Cooper et al., 2004). In particular, Cooper et al. (2004) showed that a thicker rigid chemical lithosphere leads to a thinner transitional boundary layer. This occurs because as the rigid cratonic lithosphere thickness increases, it dominates more of the combined boundary layer and thus limits the portion that can be subjected to convective stresses, and thus deformation. Therefore, this potentially deformable region (i.e., the transitional boundary layer) controls the deformation regime that is transmitted to the base of the rigid cratonic lithosphere. If this transitional boundary layer is thick, then there is a greater amount of deformable material between the convecting mantle and the rigid, cratonic lithosphere. Such interplay between the rheological sublayer (i.e., transitional sublayer) and deformation in the lithosphere has also been applied to subduction initiation (Solomatov, 2004); again, the rheological sublayer acts as a modulator for the stress environment experienced by the lithosphere.

Because the deformable sublayer buffers the cratonic lithosphere from potential destruction by convective stresses, craton stability should be highly sensitive to the thickness and rheology of the sublayer. Thus, the maximum thickness of the cratonic lithosphere should depend on a combination of the thickness of the “rigid” craton (which controls the depth of the deforming sublayer within the mantle), the rheology and thickness of the sublayer beneath the craton (which controls the deformability of that layer), and the convective vigor of the mantle (which controls the amplitude of basal tractions on the craton). Here, we examine the relationship between stress and deformation within the cratonic sublayer in an effort to constrain the role of basal tractions from mantle flow in determining the thickness of cratonic lithosphere. In doing so, we answer the question: is there a maximum thickness beyond which cratons are rendered unstable? If so, does this maximum lithosphere thickness change or remain constant with time as the mantle cools and convective stresses change?

Resistance to deformation depends on the balance between the mechanical strength of the material and the stresses available for deformation. Therefore, to define a maximum cratonic lithosphere thickness based on its resistance to deformation, we relate the asthenospheric convective stresses and associated strain rates to lithospheric thickness. We begin by estimating the strain rate within the asthenosphere as
where v is the mantle velocity relative to the surface plate, ha is the thickness of the asthenosphere, which we assume to be the difference between the depth to the base of the asthenosphere (Da) and the thickness of the upper thermal boundary layer (i.e., the thermally defined lithospheric thickness; hTBL). Assuming a Newtonian rheology, this strain rate can be converted to a convective stress that applies to the base of the craton:
where ηa is the asthenospheric viscosity. We can simplify this relationship by normalizing by a reference case with a lithospheric thickness of hr. This normalization removes the dependence on asthenospheric viscosity and velocity and reduces the relationship to one solely between convective stress and lithospheric thickness:

This nondimensionalization allows for efficient comparison between the two key parameters investigated (convective stress and lithospheric thickness), and it predicts that convective stresses increase with increasing thermal lithospheric thickness, in agreement with previous studies (e.g., Conrad and Lithgow-Bertelloni, 2006). The inherent dependence of stress on velocity and viscosity (i.e., two additional “unknowns”) is maintained within the reference case, which will be useful, because Equation 3 does not provide any additional information about the lithosphere's deformability. Note that we have assumed a simplified (Newtonian) rheology for the asthenosphere, as have others (e.g., Richards et al., 2001; Conrad and Lithgow-Bertelloni, 2006), in order to provide a framework for the stresses available for deformation within the overlying lithosphere. In other words, at this point within our analysis, we are more concerned about the effect of basal stresses on the lithosphere than we are with the nature and style of deformation within the asthenosphere. Next, we employ differing rheologies for lithosphere and asthenosphere, assuming that the two layers are acting as distinct materials. Beneath cratonic regions, this may be valid because the asthenosphere is actively shearing, whereas the cratonic lithosphere has remained relatively stable for much of Earth's history. In addition, the textures of xenoliths sampled within cratonic regions also point to a possible shift in material properties between the lithosphere and asthenosphere (Harte, 1983; Menzies, 1990).

To address potential response of the cratonic lithosphere to convective stresses, we assume a power-law rheology that relates these stresses to the local strain rate of the transitional sublayer (TSL) beneath the craton's chemical lithosphere (Fig. 1),
where A is a pre-exponential term related to the magnitude of viscosity, and n is a power-law exponent that we later assume to be 3 (Kohlstedt et al., 1995). Again, we normalize by a reference case (ε̇r = Aσrn) using the corresponding convective stress (σr) of the reference case. This provides a nondimensional relationship between the local strain rate of the craton's transitional boundary layer (which defines its deformability), and the mantle convective stresses:
Then substituting Equation 3 for σ/σr, the equation becomes:

This relationship then relates ε̇TSL, an indicator for lithospheric deformation rate, to the thermal lithosphere thickness of the craton. As we found for mantle tractions (Eq. 3), the local strain rate of the cratonic lithosphere increases with increasing total lithosphere thickness. However, it does so much more dramatically because of the power-law relationship between stress and strain rate. Note that if the asthenosphere were also represented by a power-law relationship rather than a Newtonian rheology, the mantle tractions and accompanying lithospheric deformation would still increase with lithospheric thickness, although not as dramatically.

Using our derived relationships, we can now further explore the potential for craton deformation as determined by its thickness. To begin, we compare our scaling relationships against numerical simulations based on the global mantle-flow models of Conrad and Lithgow-Bertelloni (2006). These finite-element (CitcomS model of Moresi et al. [1996] and Zhong et al. [2000]) simulations use mantle density heterogeneity inferred from seismic tomography (the S20RTSb model of Ritsema et al. [2004]) to drive mantle flow beneath a stationary lithosphere with laterally varying thickness. To assign a viscosity structure for the lithospheric and asthenospheric layers in continental regions, Conrad and Lithgow-Bertelloni (2006) estimated a characteristic length scale for the lithosphere thickness using lithospheric ages for oceanic regions and the thickness of fast seismic velocity anomalies for continental regions. A characteristic thickness of 100 km was imposed as a maximum value for the oceans and a minimum value for the continents. Conrad and Lithgow-Bertelloni (2006) then used these thicknesses to assign error function temperature profiles above 300 km depth, and they used a temperature-dependent viscosity to impose elevated viscosity in the lithosphere. We reproduce the flow models of Conrad and Lithgow-Bertelloni (2006) here, using an activation energy of 300 kJ/mol to produce a rapid viscosity increase between the asthenosphere, which has a viscosity (ηa) of 0.1ηum (one-tenth the upper-mantle viscosity), and the upper lithosphere, which has a viscosity that cannot exceed 1000 ηum.

We measured the magnitude of basal shear tractions at the depth at which the lithospheric viscosity is 10 times larger than the asthenospheric viscosity (i.e., at η = 10ηa). For a reference model, we performed a calculation in which the entire lithosphere was assigned a characteristic thickness of 100 km; this produced a lithospheric viscosity structure where η = 10ηa occurred at 87.1 km depth. Therefore, to test our prediction of the dependence of basal shear stress on lithospheric thickness (Eq. 3), we measured σ and σr from the two flow calculations and plotted their ratio as a function of (Fig. 2). When grouped into 10-km-depth bins, the average value hTBL of these estimates (Fig. 2, gray line) matches the analytical prediction from Equation 3 (Fig. 2, black line, which uses hrTBL = 87.1 km, and Da = 300 km) within the error bars. (The size of the error bars represents the standard deviation of all ratio values, and it gets larger with lithospheric thickness because the number of measurements decreases with thickness.) Basically, we were able to validate Equation 3 by using it to make a new analytical prediction of the results of Conrad and Lithgow-Bertelloni (2006).

To express the deformation of the transitional sublayer as a function of total lithospheric thickness, we applied Equation 6 to the traction ratios presented in Figure 2. The result, shown in Figure 3 (solid black line and gray curve and dots), shows that lithospheric deformation, as expressed by the strain rate, increases with lithospheric thickness. However, it is important to note the dramatic manner in which the strain rate increase occurs: there is more than a threefold increase in strain rate for lithosphere thicknesses between 150 and 200 km, and there are much steeper increases in strain rate for even thicker lithosphere. This also is the range where the thickness of the chemical lithosphere begins to dominate the thermal lithosphere (Fig. 3, dashed line), which reduces the thickness of the transitional boundary layer. To determine chemical lithosphere thickness relative to the thermal lithosphere, we used the ratio between these quantities,
which was constrained by Cooper et al. (2004) using numerical models.

The dramatic increase in strain rate within our calculations (Fig. 3, solid and dashed black lines) occurs over the range of observed maximum craton thicknesses obtained from xenoliths. This suggests an explanation for the lack of xenolith samples coming from depths greater than 200 km (Carlson et al., 2005; Lee et al., 2005). If the xenoliths are sampling the chemical lithosphere of a craton, then a craton's chemical lithosphere that becomes thicker than ~200 km could experience a potentially destructive environment due a lack of a sufficient buffer zone (Fig. 3). Thus, convecting mantle could modulate the maximum thickness of a rigid craton because the stress, and especially the resulting strain rate, that the mantle imparts on the cratonic base increases dramatically as the chemical lithosphere thickness increases beyond 150 km (Fig. 3). This is the region where the chemical lithosphere dominates the thermal lithosphere of a craton, which causes the transitional sublayer between the two to thin significantly, decreasing its ability to protect a cratonic lithosphere from the deforming mantle. The decline in the buffer zone thickness may explain the paucity of xenolith data from depths >200 km because cratonic material at thicknesses greater than this could be deformable.

Finally, our analysis is static, i.e., it depends strongly on an assumed constant asthenospheric flow velocity and viscosity, which is the source of the destructive deformation and thus sets the stage for the background mantle strain rate. Since the values of these parameters will be influenced by changes in the mantle's interior temperature (which certainly would occur over the history of a craton as Earth cools), one might predict that the maximum craton thickness might also change with time in response.

We can explore the stability of cratons through the past and future by scaling our method above to the Rayleigh number, Ra, which is a measure of convective vigor, and which thus affects the flow rates in the mantle and the tractions exerted by the mantle on the lithospheric base. Because Ra is inversely proportional to the mantle viscosity, it should decrease with time as Earth cools and becomes more viscous. In our original approach, we assumed constant Ra equivalent to the Rayleigh number of our reference case, Ra0. In other words, our reference case describes a single instant during a craton's long history. Relaxing the assumption of a constant Rayleigh number allows us to expand our analysis to the past, when a hotter mantle and increased convective vigor could have operated. We can also make inferences about the future, when the mantle will presumably become cooler, and convection will be more sluggish.

A changing Rayleigh number will affect several parameters within the original analysis. For instance, the viscosity scales inversely with the Rayleigh number (η~ Ra−1), while the velocity of the mantle flow scales as ~Ra, where β is a power-law exponent that is defined by the relation Nu ~ Raβ, which describes how convective heat transport (as expressed by the Nusselt number, Nu) scales with the Rayleigh number (e.g., Conrad and Hager, 2001). This relationship can be used to describe Earth's thermal evolution because it serves as an indicator of the efficiency of convection to cool Earth's interior. The β parameter expresses the sensitivity of that efficiency and is typically equal to 1/3 for standard boundary layer theory (Turcotte and Oxburgh, 1967). However, β may be closer to zero if the thermal boundary layer thickness (and convective heat flow) does not change with convective vigor, as may be required to explain the thermal history of Earth (e.g., Christensen, 1985). Within boundary layer theory, the thermal boundary layer thickness also scales with the Rayleigh numbers according to h ~ Ra−β. While this scaling does describe the oceanic lithosphere appropriately, it does not hold for continental (and cratonic) lithosphere, which it is not consumed by subduction. Therefore, within our following analysis, we only consider the variations in mantle flow (i.e., viscosity, velocity, strain rate, etc.) with the Rayleigh number, since these are global features that will influence cratonic lithosphere. In doing so, we are implicitly assuming that the boundary layer analysis that describes mantle convection associated with the oceanic lithosphere also describes convective flow beneath the cratons, even though the cratons themselves do not participate in that convection (e.g., Lenardic et al., 2003).

To determine how the maximum craton thickness might vary with mantle dynamics, we need to step through our approach accounting for the dependence of the Rayleigh number on relevant parameters. Since we are concerned with the influence of deviations from the reference case, we normalize the Rayleigh number by Ra0. Beginning with Equation 1, the asthenospheric strain rate scales as
because the velocity does (e.g., Conrad and Hager, 2001); note that hTBL does not scale with the Rayleigh number because this is the local thermal lithosphere of the craton, which does not adhere to boundary layer theory. We carry this strain rate scaling to Equation 2 for convective stresses. Since convective stresses depend on both viscosity (which scales as η ~ Ra−1) and strain rate (~Ra), Equation 2 scales as
(as in Lenardic et al., 2008, when β = 1/3). Finally, when we incorporate the power-law rheology in Equation 4, our scaling relationship for strain rate versus maximum craton thickness becomes:

Note that the exponent for the Rayleigh number term is equal to −1 when β = 1/3 and n = 3. Thus, values of Ra larger than Ra0 will cause the strain rate for a given value of hr to decrease (Fig. 3). As a result, the sharp increase of strain rate with increasing depth occurs at slightly larger lithospheric thicknesses for higher Ra. This increase in thickness occurs because the tractions exerted by the mantle on the base of the craton depend on both the mantle velocity and the mantle viscosity, as shown in Equation 2. The former increases as Ra increases, but the latter decreases more rapidly, causing a net decrease in the amplitude of basal tractions for increasing Ra. In effect, the greater convective vigor associated with larger Ra causes mantle-flow velocities to increase, but this faster flow is more decoupled from the craton because of a smaller asthenospheric viscosity. The net decrease in deformation rates at the cratonic base will be amplified as β decreases (which prevents increasing mantle-flow rates for higher Ra) or as n increases (which increases the sensitivity of craton deformation rates to increasing tractions).

We estimate the maximum craton thickness as a function of Rayleigh number by measuring the point on individual lithosphere thickness versus strain-rate curves (e.g., Fig. 3), where the slope becomes greater than one dimensionless strain-rate unit per kilometer (~225 km for the Ra = [solid black] curve in Fig. 3). The result (Fig. 4) shows that for all Ra0 values of β, the maximum thickness that a craton can sustain decreases with decreasing Rayleigh number, because the deforming basal tractions on the craton increase (until the h ~ Ra−β scaling factor that describes the thickness of an oceanic-style thermal boundary layer begins to dominate, as shown by the kinks in the curves of Fig. 4). Higher values of β maintain a more constant craton thickness over variable Rayleigh numbers because craton deformation rates are less sensitive to Ra for higher β (Eq. 7). However, regardless of the value of β, during times of higher Rayleigh number (during the past with higher mantle temperatures), craton thickness does not vary greatly, and the variations are even smaller if β is large.

We have shown that craton thickness should have remained relatively stable during much of Earth's history, despite the greater convective vigor (higher Ra) of the past mantle (Fig. 4). In the future, as the mantle cools, and the Rayleigh number decreases further to values smaller than our reference case, this relative stability should begin to deteriorate because of increasingly large tractions that are exerted by a more highly viscous mantle on the base of thick cratons. These tractions tend to deform and thin the transitional sublayer and dramatically decrease the maximum thickness that a craton can sustain as mantle convection becomes less vigorous, especially for small βvalues. Indeed, there is evidence for recent removal of cratonic material in North China (Gao et al., 2002).

Our analysis suggests that the nearly constant craton thickness that is thought to have persisted over billions of years based on xenolith data, and that is suggested by geodynamical calculations (Boyd et al., 1985; Pearson et al., 1995; Shirey et al., 2004; Sleep and Jellinek, 2008), must have occurred in a convection environment with a relatively high value of β or a high Rayleigh number, or both. Thus, a large value of β, which is associated with a slow mantle cooling rate and thus small changes in mantle temperatures (and thus Rayleigh number) since the Archean (Christensen, 1985), tends to promote craton stability. In contrast, a small value of β, which has been proposed by several authors (e.g., Christensen, 1985; Conrad and Hager, 2001), and especially a negative value (e.g., Korenaga, 2003), is associated with larger changes in mantle temperatures and Rayleigh number and thus should eventually cause larger changes in craton thickness (Fig. 4). Thus, to explain the relatively constant thickness of cratons, β should be positive and at least greater than ~0.1. Finally, our study also explains why cratons formed early in Earth's history, when the mantle was hotter and convecting more vigorously. In our view, vigorous convection associated with a hot mantle origin actually promotes the generation of thick cratons because the convective stresses that are destructive to cratons are smaller.

This work was supported by National Science Foundation (NSF) grant EAR-0914712 (to Conrad). We thank James Conder and an anonymous referee for helpful reviews that improved the manuscript.