Towards scaling laws for subduction initiation on terrestrial planets: constraints from two-dimensional steady-state convection simulations

Abstract

The strongly temperature-dependent viscosity of rocks leads to the formation of nearly rigid lithospheric plates. Previous studies showed that a very low yield stress might be necessary to weaken and mobilize the plates, for example, due to water. However, the magnitude of the yield stress remains poorly understood. While the convective stresses below the lithosphere are relatively small, sublithospheric convection can induce large stresses in the lithosphere indirectly, through thermal thinning of the lithosphere. The magnitude of the thermal thinning, the stresses associated with it, and the critical yield stress to initiate subduction depend on several factors including the viscosity law, the Rayleigh number, and the aspect ratio of the convective cells. We conduct a systematic numerical analysis of lithospheric stresses and other convective parameters for single steady-state convection cells. Such cells can be considered as part of a multi-cell, time-dependent convective system. This allows us a better control of convective solutions and a relatively simple scaling analysis. We find that subduction initiation depends much stronger on the aspect ratio than in previous studies and speculate that plate tectonics initiation may not necessarily require significant weakening and can, at least in principle, start if a sufficiently long cell develops during planetary evolution.

Background

Plate tectonics is central to many aspects of the geology and evolution of terrestrial planets. While Earth is the only planet where plate tectonics is observed, its driving mechanism and timing of initiation are still poorly understood. Subduction is thought to be the fundamental process for plate tectonics initiation because the slab pull of subducting slab contributes most to the forces that drive plate movements. On the Earth, initiation of subduction is greatly facilitated by tectonic forces associated with plate motions already occurring elsewhere (Mueller and Phillips 1991; Hall et al. 2003). Various models for subduction initiation has been proposed (e.g., McKenzie 1977; Turcotte 1977; Ogawa 1990; Mueller and Phillips 1991; Kemp and Stevenson 1996; Toth and Gurnis 1998; Stern 2004; Solomatov 2004a; Ueda et al. 2008; Nikolaeva et al. 2010), many of which involve existing plate boundaries or weak zones. Incipient subduction zones are often found near transform faults or fracture zones because of their physical weakness (e.g., Mueller and Phillips 1991; Gurnis et al. 2004).

On one-plate planets such as Venus and Mars, the absence of plate tectonics is likely to be due to the difficulty of subduction initiation in the absence of forces due to plate motions. In other words, the problem of plate tectonics initiation can be viewed as the problem of the very first occurrence of subduction. Due to the high sensitivity of viscosity to temperature, the lithosphere acts as the cold rigid thermal boundary layer that has a very high strength. On these planets, mantle convection is likely to be in the stagnant lid regime (e.g., Morris and Canright 1984; Fowler 1985; Solomatov 1995). One possible mechanism for the very first episode of subduction is due to the lithospheric stresses generated by mantle convection (Ogawa 1990; Fowler and O’Brien 2003; Solomatov 2004a). The magnitude of these stresses is relatively small compared to the lithospheric strength suggested by laboratory and field observations (e.g., Kohlstedt et al. 1995; Gurnis et al. 2004), and thus, it is usually believed that to initiate subduction, some weakening mechanisms must be present in the lithosphere.

Much effort has been devoted to understand the weakening mechanisms of the lithosphere. Several studies showed that the frictional shear stress resisting subduction at transform faults and fracture zones have to be less than 10 MPa for subduction to occur (Toth and Gurnis 1998; Hall et al. 2003; Gurnis et al. 2004). Stress drop estimates from earthquakes also indicate that fault strength may be approximately 10 MPa (Kanamori 1994; Kanamori and Brodsky 2004). Models are able to describe global reduction in the lithospheric strength, as well as localized weak zones such that plate-like features can be generated from mantle convection in a self-consistent manner (e.g., Trompert and Hansen 1998; Moresi and Solomatov 1998; Tackley 2000a; Bercovici 2001; Branlund et al. 2001; Regenauer-Lieb et al. 2001; Regenauer-Lieb and Kohl 2003; Regenauer-Lieb et al. 2006; Korenaga 2007; Landuyt et al. 2008). Various approaches have been used to deal with the creation of weak zones (Bercovici 2001; Bercovici and Ricard 2005; Landuyt et al. 2008; Branlund et al. 2001; Regenauer-Lieb and Kohl 2003). The two-phase damage theory with a grain-size dependent rheology was developed to explain the formation of weak plate boundaries and track the evolution of deformation (e.g., Bercovici and Ricard 2005; Landuyt et al. 2008; Bercovici 2012). Some studies suggested that water might play an important role in the localization of deformation (Regenauer-Lieb et al. 2001; Regenauer-Lieb and Kohl 2003; Regenauer-Lieb et al. 2006). Water also weakens the lithosphere by lowering the activation energy (Regenauer-Lieb et al. 2001; Regenauer-Lieb and Yuen 2004; Regenauer-Lieb et al. 2006) and increasing the pore fluid pressure (Kohlstedt et al. 1995).

One approach to quantify the weakening of lithosphere is to set a yield value to the rheology of the lithosphere to simulate brittle behavior (Fowler 1993; Trompert and Hansen 1998; Moresi and Solomatov 1998; Richards et al. 2001; Tackley 2000a,b; Fowler and O’Brien 2003; Solomatov 2004a; Stein et al. 2004; O’Neill et al. 2007; Stein and Hansen 2008). The yield stress can be regarded as a simplification of mechanisms that describe the strength of the lithosphere. Convection with yield stress is usually categorized into three regimes: mobile lid regime, transitional regime with some episodic failure, and stagnant lid regime (Moresi and Solomatov 1998; Tackley 2000b; Stein et al. 2004). Stein and Hansen (2008) further subdivided the transitional regime into episodically mobile and stable plate mobilization regimes. To access the conditions of a planet to have plate tectonics, some researchers presented regime diagrams in terms of Rayleigh number, viscosity contrast, and yield stress (e.g., Stein et al. 2004; O’Neill and Lenardic 2007).

A number of studies attempted to derive scaling relations for convective stresses and yield stress to extrapolate to planetary conditions (e.g., Moresi and Solomatov 1998; Fowler and O’Brien 2003; Solomatov 2004a,b; O’Neill et al. 2007; Valencia and O’Connell 2009; Korenaga 2010b; Van Heck and Tackley 2011; Stamenkovic and Breuer 2014). Yet the accurate description of these convection-induced stresses inside the lithosphere and thus the yield stress is lacking.

This study seeks to understand the stress distribution of the steady-state convecting cell with respect to various convective parameters using the pseudoplastic rheology as a first step. The goal of this study is to find a scaling law for the lithospheric stress (hereafter referred as lid stress) and the critical yield stress, which is the highest yield value at which the stagnant lid could be mobilized. Note that an alternative and, perhaps, more intuitive approach would be to assume that the yield stress is known and to try to figure out under what dynamic conditions it can be reached. However, (a) the ‘normal’ yield stress is so high that it is nearly impossible to reach and (b) given the uncertainties in the weakening mechanisms and thus the actual magnitude of the yield stress, it should be treated as an unknown.

In this study, we first examine the stress structure in steady-state stagnant lid convection and explore scaling relationships between convective parameters especially in relation to aspect ratio to develop a scaling theory for lid stress and critical yield stress. We then compare the theoretical scaling laws with numerical results. In addition, we investigate the accuracy of the Frank-Kamenetskii approximation for modeling the initiation of plate tectonics.

Rheology

Viscous creep governs the flow in the mantle as it has high temperatures and low stresses. It can be described by a constitutive relation (Hirth and Kohlstedt 2003), which is an Arrhenius function of temperature T, activation energy E, pressure P, and activation volume V with power law dependences on stress τ (second invariant of stress tensor), grain size d, water fugacity $f_{\text {H}_{2}\text {O}}\phantom {\dot {i}\!}$, and an exponential function of melt fraction ϕ:

$$ \eta = A\tau^{1-n}d^{m}f_{H_{2}O}^{r}\exp(-\alpha\phi)\exp\left(\frac{E+PV}{RT}\right), $$

((1))

where A and α are constants, R is the gas constant, and m,n, and r are exponents for grain size, stress, and water fugacity, respectively. Depending on the temperature, grain size, stress, pressure, and composition, the dominating deformation mechanism in the mantle would be different (Karato and Wu 1993; Karato et al. 1995; Hirth and Kohlstedt 2003). In the lithosphere, the major factor controlling the viscosity is temperature. Thus, the viscosity function to investigate subduction initiation is often written as:

$$ \eta = A'\exp\left(\frac{E}{RT}\right). $$

((2))

Frank-Kamenetskii approximation

Many numerical studies used use a relatively low viscosity contrast to observe plate behavior, which has limited applications to realistic planetary situations. Moresi and Solomatov (1998) investigated the convective regimes with viscosity contrast ranging from 3×10⁴ to 3×10⁷, and in Tackley (2000a), the viscosity contrast was limited to 10⁴, whereas Richards et al. (2001) and Stein and Hansen (2008) used viscosity contrast on the order of 10⁵. The viscosity contrast across the terrestrial lithosphere is many orders of magnitude higher.

The low-viscosity contrast is used because high-viscosity contrasts are difficult to treat in numerical calculations (Moresi and Solomatov 1995). Thus, the Arrhenius function is often approximated by the Frank-Kamenetskii function, which reduces the viscosity contrast by many orders of magnitude compared to Arrhenius viscosity function. This makes the problem of convection with strongly temperature-dependent viscosity more computationally tractable.

Frank-Kamenetskii approximation originated from the combustion theory. Frank-Kamenetskii pointed out that since the activation energy E is large, we can consider the rate of reaction only in a narrow range of temperature around the combustion temperature (Frank-Kamenetskii 1969). The equation for the rate of reaction is similar to the strain rate in the constitutive relations, which also has an Arrhenius form exp(−E/R T). Since convection mostly takes place in the interior of the cell where the temperature is close to the interior temperature T _i , we use the same approximation by expanding the exponent E/R T in the Arrhenius form so that the viscosity can be expressed as an exponential function of temperature only:

$$ \eta = B\exp(-\gamma T), $$

((3))

where B and γ are constants. In the interior:

$$\begin{array}{*{20}l} \eta_{i,\text{Arr}} &= \eta_{i,\text{exp}}, \end{array} $$

((4))

$$\begin{array}{*{20}l} \left(\frac{d\eta_{i,\text{Arr}}}{dT}\right)_{T=T_{i}} &= \left(\frac{d\eta_{i,\text{exp}}}{dT}\right)_{T=T_{i}}, \end{array} $$

((5))

where η _Arr and η _exp are the interior viscosities of the Arrhenius function and that of Equation 3 (hereafter referred as exponential viscosity), respectively. Equations 4 and 5 give γ in terms of activation energy and interior temperature:

$$ \gamma = \frac{E}{R{T_{i}^{2}}}. $$

((6))

This method of expanding the terms in the exponent preserves the interior viscosity and the change of viscosity with temperature close to T _i , where convection actively takes place. Some studies expanded the terms inside the exponents differently (e.g., King 2009). However, it is important to use Equations 4 and 5 to ensure the asymptotic accuracy of Frank-Kamenetskii approximation (Morris 1982; Morris and Canright 1984; Fowler 1985; Frank-Kamenetskii 1969).

Frank-Kamenetskii approximation was shown to be sufficiently accurate for the interior of the convective layer with large viscosity contrast (Solomatov and Moresi 1996; Ratcliff et al. 1997; Reese et al. 1999). Recent studies have examined convection with Arrhenius rheology and suggested slightly different scaling laws compared to convection with Frank-Kamenetskii viscosity (Korenaga 2009; Stein and Hansen 2013). Here, we assess the accuracy of the Frank-Kamenetskii approximation in predicting the values of critical yield stress.

Pseudoplastic rheology and plastic yielding

The brittle behavior of the lithosphere can be simplified with a viscoplastic rheology that causes yielding when the convective stresses exceed a plastic yield stress τ _y (Moresi and Solomatov 1998; Trompert and Hansen 1998; Tackley 2000b; Fowler and O’Brien 2003). The yield stress can be defined by Byerlee’s law (Byerlee 1978):

$$ \tau_{y} = \tau_{0} + \mu\rho gz, $$

((7))

where τ ₀ is the yield stress at the zero hydrostatic pressure, μ is the frictional coefficient, and ρ g z is the hydrostatic pressure. Viscous deformation occurs according to Equation 3 when stresses are less than the yield stress. Above the yield stress, the deformation follows a plastic flow law defined by a non-linear effective viscosity:

$$ \eta_{\text{eff}} = \frac{\tau_{y}}{\dot{e}}, $$

((8))

where $\dot {e}$ is the second invariant of the strain rate tensor. The yield stress defines a change on deformation mechanism based on the second invariant of the deviatoric stress tensor, which corresponds to the Von Mises yield criterion. In this study, we consider two types of yield stress: a constant yield stress τ _y or a depth-dependent yield stress with a constant gradient τ y′.

Methods

Formulation of the problem

Equations of thermal convection

The equations of thermal convection of an incompressible fluid in Boussinesq approximation and infinite Prandtl number are as follows:

$$\begin{array}{*{20}l} \frac{\partial u_{i}}{\partial x_{i}} &= 0, \end{array} $$

((9))

$$\begin{array}{*{20}l} \alpha\rho g_{i} T'-\frac{\partial p'}{\partial x_{i}}+\frac{\partial\tau_{ij}}{\partial x_{j}} &= 0, \end{array} $$

((10))

$$\begin{array}{*{20}l} \frac{\partial T'}{\partial t}+u_{i}\frac{\partial T'}{\partial x_{i}} &= \kappa\frac{\partial^{2}T'}{\partial {x_{i}^{2}}}, \end{array} $$

((11))

where ρ is density, p ^′ and T ^′ are pressure and temperature perturbations, g _i is the gravity vector, α is the thermal expansivity, $\kappa = \frac {k}{\rho c_{p}}$ is the thermal diffusivity, k is the thermal conductivity, and c _p is the isobaric specific heat. τ represents the elements of the stress tensor according to the following equation:

$$ \tau_{ij} = \eta\left(\frac{\partial u_{i}}{\partial x_{j}}+\frac{\partial u_{j}}{\partial x_{i}}\right), $$

((12))

where i and j, are indices of the coordinate axes.

The boundary conditions are as follows. For a cell with only base heating, the top and bottom surfaces are isothermal. The temperature of the top surface as T ₀ and that of the bottom surface as T ₁. The temperature difference is Δ T=T ₁−T ₀. The vertical boundaries are thermally insulated. All surfaces are free-slip. The velocity normal to a cell boundary is zero.

Non-dimensionalization

The above equations are non-dimensionalized as follows:

$${} {\fontsize{9.4pt}{9.6pt}\selectfont{\begin{aligned} \bar{x_{i}} = \frac{x_{i}}{d},\; \bar{u_{i}} = u_{i}\frac{d}{\kappa},\; \bar{\eta} = \frac{\eta}{\eta_{1}},\; \bar{\tau} = \tau\frac{d^{2}}{\kappa\eta_{1}},\; \bar{t} = t\frac{\kappa}{d^{2}},\; \bar{T} = \frac{T}{\Delta T}, \end{aligned}}} $$

((13))

where d is the layer depth, t is time, η ₁ is the reference viscosity (at the bottom of the convective layer), and Δ T is the temperature drop across the layer. The Rayleigh number can then be used to characterize the system:

$$ Ra = \frac{\alpha\rho g\Delta T d^{3}}{\kappa\eta_{1}}. $$

((14))

The Arrhenius viscosity is non-dimensionalized as:

$$ \bar{E} = \frac{E}{R\Delta T}, $$

((15))

$$ \bar{\eta} = \eta_{r, \text{Arr}} \exp\frac{\bar{E}}{\bar{T_{0}}+\bar{T'}}, $$

((16))

while the dimensionless exponential viscosity is as follows:

$$ \bar{\eta} = \eta_{r, \text{exp}} \exp(-\theta\bar{T'}), $$

((17))

where the pre-factors η _r,Arr and η _r,exp are chosen to ensure that the viscosity is equal to unity at $\bar {T'}=1$ and the Frank-Kamenetskii parameter θ is the non-dimensionalized form of the constant γ (Equation 6):

$$ \theta = \gamma\Delta T. $$

((18))

In this case, the viscosity contrast is characterized by only one parameter (θ):

$$ \Delta\eta = e^{\theta}, $$

((19))

The non-dimensional yield stress is as follows:

$$ \bar{\tau_{y}} = \bar{\tau_{0}} + \bar{\tau'}_{y}\bar{z}, $$

((20))

where

$$ \bar{\tau}_{0} = \frac{d^{2}}{\kappa \eta_{1}} \tau_{0} $$

((21))

is the non-dimensional yield stress at the surface and

$$ \bar{\tau'}_{y} = \frac{\rho g d^{3}}{\kappa \eta_{1}} \mu $$

((22))

is the non-dimensional yield stress gradient. The dimensionless pseudoplastic viscosity is as follows:

$$ \bar{\eta}_{\text{eff}} = \frac{\bar{\tau}_{y}}{\bar{\dot{e}}},\; \bar{\dot{e}} = \frac{d^{2}}{\kappa}\dot{e}. $$

((23))

In the following discussion, all parameters are assumed to be non-dimensionalized and the bar sign will be dropped. The non-dimensional forms of Equations 9 to 11 are as follows:

$$\begin{array}{*{20}l} \frac{\partial u_{i}}{\partial x_{i}} & = 0, \end{array} $$

((24))

$$\begin{array}{*{20}l} Ra T'\textbf{e}_{i} - \frac{\partial p'}{\partial x_{i}} + \frac{\partial\tau_{ij}}{\partial x_{j}} & = 0, \end{array} $$

((25))

$$\begin{array}{*{20}l} \frac{\partial T'}{\partial t} + u_{i}\frac{\partial T'}{\partial x_{i}} & = \frac{\partial^{2} T'}{\partial {x_{i}^{2}}}. \end{array} $$

((26))

where e _i is a unit vector in the direction of gravity.

Matching Arrhenius viscosity and exponential viscosity in non-dimensional form

To compare the two viscosity laws, the Arrhenius viscosity and the exponential viscosity are matched according to Equations 4 and 5:

$$\begin{array}{*{20}l} {}\eta_{\mathrm{r, Arr}}\exp\frac{E}{T_{i}+T_{0}} & = \eta_{\mathrm{r, exp}}\exp{(-\theta T_{i})}, \end{array} $$

((27))

$$\begin{array}{*{20}l} {}\frac{E}{(T_{i}+T_{0})^{2}}\eta_{\mathrm{r, Arr}}\exp\frac{E}{T_{i}+T_{0}} & = \theta\eta_{\mathrm{r, exp}}\exp{(-\theta T_{i})}. \end{array} $$

((28))

Equations 27 and 28 would yield:

$$ \theta = \frac{E}{(T_{i}+T_{0})^{2}}. $$

((29))

T _i differs from the bottom temperature T ₁ by a rheological temperature difference Δ T _rh, which is on the order of θ ⁻¹.

Equation 29 shows that there are various combinations of E and T ₀ that would give the same θ, and they would result in different Arrhenius viscosity contrasts:

$$ \Delta\eta_{\text{Arr}} = \exp\frac{ET_{i}}{T_{0}(T_{i}+T_{0})}. $$

((30))

Thus the ratio of Arrhenius viscosity contrast to exponential viscosity contrast expθ is as follows:

$$ \frac{\Delta\eta_{\text{Arr}}}{\exp(\theta)} = \exp\frac{ET_{i}}{T_{0}(T_{i}+T_{0})^{2}}. $$

((31))

Results

Steady-state convection

We use the finite element code CITCOM (Moresi and Solomatov 1995) to simulate convection in a 64a×64 box, where a is the aspect ratio. Several high viscosity cases were ran with 128a×128 resolution for more accurate results. All cases were run until they reached a steady state. We consider the range of parameters in which convection is in the stagnant lid regime (Solomatov 1995).

The structure of a steady-state convection cell is shown in Figure 1. Due to the temperature-dependent viscosity, the top part of a convective cell forms a stagnant lid and convection only penetrates into the lid by length of δ _rh - the rheological layer thickness (e.g., Solomatov 1995). A cold rigid lid, which is often defined by an isotherm, is naturally developed in the top part of the cell sloping downward to the downwelling end. The stress field (Figure 1, right) shows a stress boundary layer near the surface. This is consistent with the analytical solutions of Fowler (1985).

Figure 1

Temperature (left) and stress fields (right) of a steady-state convective cell. Color scale goes from high (red) to low (blue). The lid is defined by an isotherm T _L , and the interior temperature T _i is found by averaging the temperature of the convecting interior excluding boundary effects. For scaling purposes, the lid slope λ and rheological sublayer thickness δ _rh is taken at mid-width, whereas lid thickness δ _lid is extrapolated to the edge from the lid slope in the middle.

To compare the stresses in exponential and Arrhenius viscosity, we choose a range of T ₀ and calculate their corresponding E that gives the same θ according to Equation 29 with T _i ≈1.

Aspect ratio

The horizontally averaged profiles (Figure 2) show that the lid thickness slightly depend on the aspect ratio of the convective cell but the bulk of the temperature and viscosity profile does not vary much with the aspect ratio. However, the stress profile in small aspect ratio cells (a=0.25) is distinctly different from that in larger cells (a = 0.5 to 1).

Figure 2

Comparison of (a) viscosity, (b) temperature, (c) stress profiles of exponential viscosities for θ=16, R a=3×10⁷, and varying a. The somewhat different stress profile of a=0.25 suggests that the surface stress boundary layer is not the region with highest stress.

This difference is more apparent in the 2-D plots (Figure 3). In wider cells, the layer with highest stress (red) is approximately symmetric along the half-width of the cell, increasing in depth towards both edges and greater towards the downwelling edge. Below the surface stress boundary layer, the stresses in the middle of the lid are highest as they are not affected by the free-slip boundary conditions at the vertical edges. Although surface stress boundary layer is obvious in horizontally averaged stress profiles, the 2-D stress fields reveal that the surface stress are not always greater than that at depth. Figure 3 shows that at mid-width, it is possible that the surface stresses are lower than the interior. The high-stress region (orange to lime) roughly correspond to the cold lid shown in the temperature distribution, both having slopes towards the downwelling edge. In the narrowest cell however, this high-stress slope deviates from the thermal lid slope. There is a high-stress ‘core’ within the lid where the magnitude of stresses is close to that of the surface stress boundary layer. This implies that steady-state convection in small aspect ratios may be in another regime of plastic failure where the interior stresses reach the yield stress first, such that the plastic zone could initiate at depth while the surface may or may not be plastic, depending on the yield stress.

Figure 3

Stress (top row) and temperature fields (bottom row) of convecting cells with R a=3×10⁷, θ=16, and various a.

For cases with high-viscosity contrasts and larger aspect ratios, convection is localized to a small part of the cell, or there could be multiple convective cells (Figure 4). These cases might belong to the subcritical convective regime (Solomatov 2012), which has a different behavior from supercritical convection and therefore not in our scope of study. To stabilize one-cell flow, smaller aspect ratio cases are considered to investigate the dependence of critical yield stress on aspect ratio.

Figure 4

Localized convection and multiple subcells in steady-state convection.(a) θ=19,T ₀=0.6,a=1,R a=1×10⁷. (b) θ=16, T ₀=0.6, R a=3×10⁷.

Rayleigh number

Increasing Ra reduces the lid thickness as well as the lid slope (Figure 5). The stresses are larger with higher Ra, as is expected with more vigorous convection.

Figure 5

Stress structure of steady state cases with θ=16,a=0.75,T ₀=0.8.(a) R a=10⁷. (b) R a=3×10⁷. (c) R a=10⁸. The area of convecting interior becomes larger as Ra increases and the thickness of the stagnant lid decreases. Scale shows the values of logτ.

Viscosity contrast

The effects of viscosity contrast on the interior profile are illustrated in the plots in Figure 6. The conductive lid becomes thicker and the interior temperature is closer with the bottom temperature with increasing θ. In the stagnant lid regime (θ>10), the stress boundary layer near the surface is more pronounced than in the transitional regime.

Figure 6

Temperature, viscosity, and stress profiles of R a=3×10⁷, a=0.75, and varying Frank-Kamenetskii parameter θ.

The interior viscosity, temperature, and stress for Arrhenius viscosity and exponential viscosity are close. We investigate a range of different Arrhenius viscosity values by varying T ₀, noting that at T ₀=2.0 is a rather high surface temperature. As T ₀ decreases, the difference in viscosity contrasts calculated by the two viscosity laws becomes larger. However, the temperature and the stress profiles are similar, as shown in the horizontally averaged profiles in Figure 7. There is only a slight decrease in thickness of the stress boundary layer as the Arrhenius viscosity contrast increases. The 2-D stress distributions in Figures 8 and 9 reflect the small differences in the stress distribution as viscosity contrast increases, both in exponential viscosity and Arrhenius viscosity. This is contrary to the findings of Stein and Hansen (2013), which observed a distinctly thinner lid in Arrhenius viscosity and slightly different temperature and viscosity profiles. The small difference that we observe may affect the scaling laws for subduction. The vertical variation in stresses at the downwelling edge may be particularly important. Although the lid thickness is about the same, the contrast in stresses at x=a seems to be greater at low T ₀ or high-viscosity contrast (Figure 9). Since subduction occurs at the downwelling edge, this may influence the scaling of yield stress.

Figure 7

Comparison of temperature, viscosity, and stress profiles of Arrhenius and exponential viscosities for θ=16, R a=3×10⁷, a=0.75, and various T ₀.

Figure 8

Stress structure of steady state cases with R a=3×10⁷, T ₀=0.8, a=0.75.(a-d) Arrhenius viscosity with (a) θ=22, (b) θ=19, (c) θ=16, (d) θ=13, and (e) exponential viscosity with θ=13. Color scale shows the values of log(τ). White lines represent streamlines.

Figure 9

Stress structure of steady state cases with R a=3×10⁷, θ=16, a=0.75. Cases (a-d) use Arrhenius viscosity with (a) T ₀=0.6, (b) T ₀=0.8, (c) T ₀=1.2, and (d) T ₀=2.0. Case (e) uses exponential viscosity. Color scale shows the values of log(τ). White lines represent streamlines.

Lid stress scaling theory

Fowler (1985) obtained a polynomial expression for the stress in the lid below the surface stress boundary layer in large lid slope approximation, which allows a comparatively simple scaling relation for stress. In order to solve the equations of convection, Fowler took Ra and θ to be asymptotically large and assumed the magnitude of lid slope to be either on the order of lid thickness or rheological sublayer thickness. For the Earth and some smaller terrestrial planets as well as most numerical simulations, Ra may not be as high as the asymptotic theory require, and thus, some other theory may be needed for lid stress scaling. Moreover, as we find in this study, the lid slope does not follow either of these two end-member cases and thus needs to be scaled based on numerical simulations.

Fowler also found that the interior flow can be uncoupled from the rheological sublayer, which makes the problem setup akin to a viscous lid gravitationally sliding along a slope. We can therefore estimate the shear stress in the rigid lid τ _lid by considering the force balance on the lid (Figure 1):

$$ \frac{d\tau}{dy} = \Delta\rho g_{x}. $$

((32))

For density changes due to temperature variations, the lid stress can be integrated from Equation 32:

$$ \tau_{\text{lid}} = -\alpha\rho_{0}\frac{dT}{dy}g\sin \lambda\frac{y^{2}}{2} + \tau_{i}, $$

((33))

where τ _i is the stress at the interior temperature T _i .

There are two assumptions that allow us to simplify Equation 33. The first is that τ _i is negligible since τ _i ≪τ _lid. The second is the small lid slope approximation. For small λ, sin(λ) is approximately equal to λ. The non-dimensional form of Equation 33 becomes:

$$ \tau_{\text{lid}} = -\text{Ra}\frac{dT}{dy}\frac{y^{2}}{2}\lambda. $$

((34))

Thus, the lid stress is determined by Ra, d T/d y, λ, and y, which will be defined in the following discussion. This scaling is similar to that obtained from the analytical solutions of Fowler (1985).

Lid base temperature

The lid base is often defined by an isotherm:

$$ T_{L} = T_{i} - C\theta^{-1} $$

((35))

where C is a constant. We determine the lid base from velocity profile (Solomatov and Moresi 2000) to find the constant C. We first find the greatest velocity gradient at a specific distance x from the upwelling edge of the cell. This velocity gradient is then extended to the depth at which velocity is zero, as shown in Figure 10. This depth defines the lid thickness. This process is repeated for all x (from 0 to a) to obtain the shape of the lid base across the convecting cell. This velocity gradient-defined lid base is then used to find the thermal lid base defined by Equation 35. The interior temperature T _i is found by averaging the temperature in the middle part of the interior to exclude the boundary effects. The constant C is determined by matching the lid thickness at mid-width (x=0.5a) given by T _L and that defined by the velocity gradient (Figure 11). We choose the value at mid-width because in most cases, the two lid bases are closest around the middle of the cell for a large lateral extent. In some cases, especially for narrower cells, the temperature lid base and the velocity lid base may not match, so the mid-width serves as a reference point for consistency in defining the lid base.

Figure 10

Velocity profile taken at mid-width of the cell (x=0.5a). Dotted line is the linear extrapolation of the maximum velocity gradient, and the lid base is marked at the depth at which this line intersect with the vertical axis.

Figure 11

Top part of a convective cell showing various definitions of the lid base. The lid base defined by the velocity gradient is shown in dashed line, and it matches the temperature lid base defined at T _L =T _i −3.2θ ⁻¹ at around the mid-point. The lid slope is estimated at the mid-point.

Equation 35 allows us to determine the rheological temperature difference Δ T _rh=T _i −T _L =C θ ⁻¹. It varies with Ra, θ, and aspect ratio, and the scaling exponents are summarized in Table 1.

Table 1 Numerical results of power law coefficients in scalings of different parameters with Ra, aspect ration a , and Frank-Kamenetskii parameter θ

Lid slope and lid thickness

Gravitational sliding requires a downward dipping slope λ as indicated in Equation 34. In larger aspect ratio cells, although the lid thickness varies horizontally, the lid slope is approximately constant in the middle portion of the cell. This is different for smaller cells, where the lid base could be some function of x instead of a straight line (Figure 12). For example, Fowler suggests that the lid base varies with x ^0.4 (Fowler 1985). For consistency in our scaling analysis, the lid slope is taken to be the slope of the thermal lid base at mid-width.

Figure 12

The slope of the lid base at various aspect ratios. The lid slope deviates from the linear approximation as the cell aspect ratio and viscosity contrast increase; therefore, a non-constant value of the lid slope may affect the scalings.

To check whether the lid slope scales with lid thickness or rheological sublayer thickness (Fowler 1985), we look at the vertical drop δ _x (Figure 1) and define the lid slope as δ _x/a. If the lid slope scales with the lid thickness, then δ _x∼δ _lid∼Nu⁻¹, and thus

$$ \delta_{\mathrm{x}}/\delta_{\text{lid}} \sim \text{Nu}\delta_{\mathrm{x}} \sim \text{constant}. $$

((36))

If the lid slope scales with δ _rh, then δ _x∼δ _rh. As δ _rh/δ _lid∼θ ⁻¹, thus δ _x/δ _lid∼θ ⁻¹, or

$$ \theta\delta_{\mathrm{x}}/\delta_{\text{lid}} \sim \text{Nu}\delta_{\mathrm{x}}\theta \sim \text{constant}. $$

((37))

We plot Nu δ _x and Nu δ _x θ in Figure 13. Neither combination remains constant with θ, with Nuδ _x increases with θ and Nuδ _x θ decreases with θ. This suggests that the lid slope is somewhere in between these two extreme cases. Therefore in deriving the scaling laws for stresses, we need to determine the dependence of lid slope on various convective parameters (Table 1).

Figure 13

Plots of Nu δ _x (or δ _x/δ _lid) and Nu δ _x θ (or θ δ _x/δ _lid) as a function of θ.

Thermal gradient

The temperature is approximately a linear function of depth and the thermal gradient in the lid is about constant with depth. At mid-width (x=0.5a), it is approximately equal to the Nusselt number, which is the non-dimensional horizontally averaged surface temperature gradient. Since we are looking at temperature changes from the interior to the bottom of the lid which includes the rheological sublayer, we also check the thermal gradient in the rheological sublayer Δ T _rh/δ _rh to note any difference in the scaling relations. As before, we choose the values Δ T _rh/δ _rh at the mid-width to exclude boundary effects for scaling purposes.

In previous theories, Δ T _rh/Δ T∼θ ⁻¹ and δ _rh/δ _lid∼θ ⁻¹. The determination of C follows the description in the previous section on lid base temperature, and it is found to be dependent on aspect ratio and θ. Therefore, Δ T _rh/Δ T and δ _rh/δ _lid will also have a dependence on a and θ, and their scaling relations are summarized in Table 1.

Stress scaling

Equation 34 can now be expressed in a non-dimensional form with the definitions of various parameters in the previous discussion. At the bottom of the lid at δ _rh from the interior, where the temperature difference is Δ T _rh, the stresses are as follows:

$$ \tau_{\text{lid}} \sim \text{Ra}\frac{\Delta T_{\text{rh}}}{\delta_{\text{rh}}}\lambda y^{2}. $$

((38))

Further into the lid, the non-dimensional temperature gradient is the Nusselt number. Therefore Equation 34 can be alternatively scaled as:

$$ \tau_{\text{lid}} \sim \text{Ra}\text{Nu}\lambda y^{2}. $$

((39))

As shown in Figure 10, there is a slight difference in the thermal gradient in the lid and in the rheological sublayer; therefore, Equations 38 and 39 may result in slightly different scaling exponents. We check both scaling relations to see whether the stresses at the lid base and those in the lid can be scaled similarly.

We plot the stress profile according to Equations 38 and 39 and compare with that from numerical solutions (Figure 14). The prefactor of the stress as a function of y calculated from the thermal gradient in the rheological sublayer (T _rh/δ _rh) is 5.9 ×10⁶, whereas that from the Nusselt number is 9.6 ×10⁶. This demonstrates that the theoretical stress profiles match fairly closely with the numerical one, and the best fit can be obtained with some small adjustments in the coefficient.

Figure 14

Comparison of stress profiles at mid-width obtained from numerical calculations and theory. Blue line represents the best fit to numerical solutions of stresses. The stress profile is taken at mid-width, where the stresses in the surface boundary layer is lower than the interior as shown in the 2-D stress fields.

All the above parameters depend on Ra, aspect ratio a, and θ; thus, they can be expressed as Ra ^β a ^ζ θ ^α where β, ζ, and α are scaling exponents. The results are summarized in Table 2.

Table 2 Power law coefficients in scaling laws for stresses: numerical results vs theory

Convection with yield stress

We use the steady-state solutions as the starting point before imposing a yield stress to simulate plastic yielding. For the yield stress gradient, a small cohesion term (surface yield stress) was introduced to stabilize the solution.

Regimes of convection with constant yield stress or constant yield stress gradient

When a yield stress is present, the regions with stresses that reach the yield value would have an approximately constant stress close to τ _y . These plastic zones develop first at the corners of the cell where the stresses are highest. As τ _y decreases, the plastic zones extend both in depth and horizontally, narrowing the width of the high-viscosity part of the lid (Figure 15). If yield stress is too high, the depth of the plastic zone is small or the plastic zone is entirely absent; thus, the stagnant lid does not fail. If the yield stress is sufficiently low, the plastic zone extends suffcieintly deep so that the stagnant lid is mobilized.

Figure 15

Stress (top) and viscosity (bottom) fields of case R a=3×10⁷, θ=16, a=0.75, T ₀=0.8.(a) τ _y =6.7×10⁵, (b) τ _y =7×10⁵, (c) τ _y =∞. Failure occurs at τ _y =6.6×10⁵, while the stagnant lid remains at higher τ _y . Color bars show log τ values.

We examine the stress and viscosity profiles at various locations x to see how they change in the presence of a yield stress (Figures 16 and 17). At high yield stress or high yield stress gradient, the plastic zone only occurs at shallow depths and the bulk of the stress and viscosity profiles are unaltered from the stagnant lid state. The plastic zone extends deeper as the yield stress or yield stress gradient decreases. As the yield stress approaches the critical value, a small change in yield stress induces a change in plastic depth that is comparable to the change caused by an order of magnitude change in yield stress when it is far from critical. Although subduction does not occur, this implies a change in convection regime from stagnant lid to some sort of transitional regime.

Figure 16

Stress profiles at different widths at various yield stress or yield stress gradient.

Figure 17

Viscosity profiles at different widths at various yield stress or yield stress gradient.

In this transitional regime, the yield stress slightly changes the interior dynamics as can be seen in various convective parameters of the interior region of the cell. The yield stress increases the lid slope whereas the lid thickness remains approximately the same. When the yield stress is slightly above the critical value, these changes caused by yield stress are negligible and convection remains in the stagnant lid regime. Therefore in deriving the critical yield stress scalings, we refer to the steady-state structure that has a yield stress just above the critical value, so that the scaling relations for various convective parameters (Table 1) from steady-state stagnant lid convection can still be used.

Time evolution of lid weakening and failure

When the lid fails, the surface velocity continuously increases and overturn occurs (Figure 18). Figure 19 shows the time sequence of stress structures before and during failure. When the surface velocity is still low compared to the bottom velocity (Figure 18, left), the variation in stress structure is not obvious. It is not until the velocity begins to increase drastically that the plastic yield zones from the two corners start to connect in the middle of the cell to form a plastic lid. The weak lid then becomes unstable and starts to subduct.

Figure 18

Surface velocities over time for R a=3×10⁷, θ=16, a=0.75, (a) T ₀=0.8, τ _y =6.6×10⁵, (b) T ₀=0.4, τ _y =6.1×10⁵, bottom velocity ≈650. The surface velocity in (b) increases slowly. Although it has evolved for approximately 25 times as much as the period for (a) to fail, it is far from reaching the bottom velocity.

Figure 19

Snapshots of stress fields before and at the point of failure for case R a=3×10⁷, θ=16, a=0.75, T ₀=0.8, τ _y =6.6×10⁵. Time sequence goes from left to right and top to bottom. White arrows show velocities. (a) and (b) are close to beginning of simulation, (c) is at the mid-point between the start and failure, and (d-i) are right before overturning (depicted in (j-l)) occurs.

A possible contribution to the uncertainty in determining the critical yield stress is that at the vicinity of the critical value, the behavior may be difficult to interpret. In some cases that as the yield stress gets close to a critical value, the surface velocity increases slowly and it may take more than 10⁵ timesteps to reach a point of overturn, whereas typically, it takes less than 10⁴ timesteps to a drastic increase in surface velocity (Figure 18 right). This may be due to the behavior of dynamic system near a critical point.

Critical depth of plastic failure zone

For subduction to occur, the lithosphere has to be sufficiently weak so that it can be mobilized by stresses arisen from mantle convection. As seen in Figures 16 and 17, the lid remains stagnant if the plastic yield zone is small. Therefore, the question is how deep does this plastic zone have to penetrate for the lid to be mobilized? Previous theories proposed that the plastic zone has to penetrate through some critical depth δ _pl defined by a critical temperature T _c :

$$ \frac{\delta_{\text{pl}}}{\delta_{0}} = \frac{T_{c} - T_{0}}{T_{i} - T_{0}}, $$

((40))

where δ ₀=δ _lid+δ _rh. The model of Fowler and O’Brien (2003) predicts that δ _pl is defined by the temperature that gives the interior viscosity. For Newtonian rheology, this means that the plastic zone has to extend through the base of the lid. Solomatov (1995,2004a) suggested that δ _pl only has to penetrate to the isotherm at which the viscosity contrast with the interior viscosity is e ⁴⁽ⁿ⁺¹⁾, where n is the stress exponent for non-Newtonian viscosity.

To examine these hypotheses, we look at the stress profile at the downwelling edge (x=a) to determine the depth of the plastic zone, as the stresses at this edge are the highest and this is where subduction starts (Figure 20). The depth of the plastic zone is defined by extent of the stress modified by the plastic flow law in Equation 7. The exact value of δ _pl and the stress at this depth are found by the intercept of the stress calculated from linear extrapolation from the top (where stress is determined by the yield stress) and exponential extrapolation from the creep regime just below the plastic depth. The values of δ _pl for constant yield stress cases and that for constant yield stress gradient cases are close. It is noted that the stress at depth δ _pl in constant yield stress gradient cases is about twice as much as that in constant yield stress cases. The force on the lid is δ _pl τ _{y,c r} for constant yield stress and 0.5δ _pl τ y,c r′ for constant yield stress gradient. This implies that the force on the lid is about the same for both constant yield stress and yield stress gradient.

Figure 20

Depth of plastic zone determined by the yield stress (right) and the drop in viscosity due to the yield stress (right). Profiles taken at the downwelling edge (x=a) of the convecting cell.

As we see from the effects of aspect ratio on steady-state stress distribution, the stresses in the interior become comparable or even exceed the surface stresses as the aspect ratio and Ra decrease and viscosity contrast increases. This means that the plastic zone may not propagate from the top but also develop at depths in the lid, and the plastic zone does not span the whole top part of the cell. (Figure 21). This may represent another regime of lid failure. These cases are thus excluded from our scaling analysis.

Figure 21

Stress field (left) and stress profile at the downwelling edge. R a=3×10⁶, a=0.s x5, θ=19, τ _y =3.2×10⁴, exponential viscosity.

The depth of the plastic zone is determined from drop in stress by the yield stress. We note that the zone of reduced viscosity due to the yield stress may correspond to the zone of reduced stress as shown in Figure 20. However, this does not always hold, especially for cases with higher θ. Figure 22 shows that the transition of the plastic zone to creep flow may not correspond to the sharp change in the viscosity. This means that the viscosity at the plastic depth and the maximum viscosity are different, since the reduction in viscosity is not only determined by the yield stress but also by the strain rate.

Figure 22

Strain rate, stress, and viscosity profile at the edge with R a=3×10⁶, a=0.75, θ=19, exponential viscosity. The point of maximum viscosity may not correspond to the brittle-plastic transition.

To find out whether there is a critical viscosity contrast in lithospheric failure, we examine both the viscosity contrast at plastic depth Δ η _pl and the maximum viscosity contrast Δ η _max. The maximum viscosity needs to be determined by extrapolation as the resolution near the point where the viscosity is maximum is not high enough to resolve sharp changes in stress and viscosity. The point of maximum viscosity is found by extrapolating the values from both above and below the maximum point (e.g., Figure 22 right). The viscosity is extrapolated linearly from the two points above the maximum point, and below the maximum point, the viscosity is calculated from temperature. The intersection of these two curves determines the maximum viscosity and its depth. We found that these two viscosity contrasts are mostly within the same order of magnitude. Figure 23 shows that the maximum viscosity seems to depend on the original non-yielding viscosity contrast Δ η. For exponential viscosity and Arrhenius viscosity cases at low non-yielding viscosity contrast, the maximum viscosity contrast Δ η _max increases with non-yielding viscosity contrast. This may be because they are close to transitional regime. At higher Δ η, the increase in Δ η _max seems to decline with increasing Δ η, but the spread of data prevents us from concluding that Δ η _max converges towards higher viscosities. As with the maximum viscosity contrast, Δ η _pl also does not display linearity or convergence clearly with either the non-yielding viscosity contrast or θ (Figure 24).

Figure 23

Maximum apparent viscosity in the convective cell at the point of lithospheric failure as a function of the viscosity contrast in the absence of yield stress. Black, exponential viscosity with constant τ _{y,c r} ; red, Arrhenius viscosity with constant τ _{y,c r} ; green, exponential viscosity with constant τ y,c r′; blue, Arrhenius viscosity with constant τ y,c r′.

Figure 24

Viscosity contrast at plastic depth, fraction of lid that experience plastic failure, and the distance between the interior and plastic depth as a function of θ. Solid symbols represent constant τ _{y,c r} cases, while open symbols are constant τ y,c r′.

Since the critical viscosity contrast is neither a constant or a function of θ, we look at the depth of the plastic zone to derive scaling relations for the yield stress. We investigate the plastic depth δ _pl as a fraction of the lid thickness. For scaling purposes, the lid thickness was previously defined at the middle of the convecting cell. Here, since δ _pl is defined at the downwelling edge, we have to determine a lid thickness at the edge δ _lid,max. This is done by extrapolating the mid-width lid slope to the downwelling edge (Figure 10). As shown is Figure 24, the plastic depth is approximately 0.3 to 0.5 of the lid thickness. We take approximate values for our scaling relations rather than scaling these properties with convective parameters because the trends observed in Figure 24 maybe due to insufficient viscosity contrasts which place convection on the boundary of transitional regime, especially for θ=13 in which the viscosity is reduced to approximately 10⁴ by the critical yield stress.

To find the lid stress at δ _pl using Equation 34, we also need to determine the distance of the base of the plastic zone from the convective interior y _pl. The lid base is at δ _rh from the interior, so we express y _pl in terms of δ _rh to give a sense of distance in relation to rheological sublayer thickness. While there is a general trend of increasing of y _pl/δ _rh with θ, it is difficult to discern a correlation as y _pl/δ _rh fluctuates; thus, we take y _pl≈3δ _rh for our scaling relations.

The dependence of δ _pl, y _pl , and Δ η _pl on Ra and aspect ratio are very weak and therefore assumed negligible.

Scaling for critical yield stress and critical yield stress gradient

Critical yield stress scaling theory

In deriving a theoretical scaling for the critical yield stress, we assume an approximate balance between the force generated by the shear stresses acting at the base of the lid and the normal stress acting on the side of the lid and we assume that the latter are largely dominated by the stresses in the plastic zone δ _y (Solomatov 2004a and Figure 25). Thus, we can express the yield stress as:

$$\begin{array}{@{}rcl@{}} \tau_{y} & \sim & \tau_{\text{lid}}\frac{a}{\delta_{y}}\\ & = & \text{Ra}\frac{dT}{dy}\lambda\frac{y_{\text{pl}}^{2}}{2}\frac{a}{\delta_{y}}. \end{array} $$

((41))

Figure 25

Schematic diagram of surface stresses on the plastic zone in the lid. The shear stress τ _lid acting on the base of the lid of horizontal length a is balanced by the normal stress τ _y acting on the side with depth δ _y , developed under free-slip boundary conditions.

In the case of a constant yield stress, gradient (τ _y =τ y′z), τ y′ can be scaled as

$$\begin{array}{@{}rcl@{}} \tau'_{y} & \sim & \tau_{\text{lid}}\frac{a}{{\delta_{y}^{2}}}\\ & = & \text{Ra}\frac{dT}{dy}\lambda\frac{y_{\text{pl}}^{2}}{2}\frac{a}{{\delta_{y}^{2}}}. \end{array} $$

((42))

In Equations 41 and 42, the yield stress is treated as the normal stress whereas the lithosphere stresses are shear stress. However, the magnitude of stresses is expressed in second invariant, and the yield stress in von Mises criterion also put a limit the second stress invariant. To see which stress component contributes to the second invariant, we plot the normal stress and shear stress profiles in Figure 26. Inside the lid where plastic failure occurs, normal stress dominates, whereas shear stress exceeds normal stress below the plastic depth.

Figure 26

Horizontally averaged stress components and second invariant of stress as a function of depth.

The critical yield stress is the stress at δ _y =δ _pl . The critical plastic depth δ _pl is taken to be 0.3 to 0.5 δ _lid, and y _pl is about 2 to 4 δ _rh .

$$\begin{array}{*{20}l} \tau_{\text{lid}} & \sim C_{1}\text{Ra}\Delta T_{rh}\lambda\delta_{rh}, \end{array} $$

((43))

$$\begin{array}{*{20}l} \tau_{y,cr} & \sim C_{2}\text{Ra}\Delta T_{rh}\lambda a\frac{\delta_{rh}}{\delta_{\text{lid}}}, \end{array} $$

((44))

$$\begin{array}{*{20}l} \tau'_{y,cr} & \sim C_{3}\text{Ra}\Delta T_{rh}\lambda a\frac{\delta_{rh}}{\delta_{\text{lid}}^{2}}. \end{array} $$

((45))

For τ _lid expressed in terms of Nu (Equation 39), noting that Nu $\sim dT/dy \sim \delta _{\text {lid}}^{-1}$,

$$\begin{array}{*{20}l} \tau_{\text{lid}} & \sim C_{1}\text{Ra}\text{Nu}\lambda\delta_{rh}^{2}, \end{array} $$

((46))

$$\begin{array}{*{20}l} \tau_{y,cr} & \sim C_{2}\text{Ra}\text{Nu}^{2}\lambda a\delta_{rh}^{2}, \end{array} $$

((47))

$$\begin{array}{*{20}l} \tau'_{y,cr} & \sim C_{3}\text{Ra}\text{Nu}^{3}\lambda a\delta_{rh}^{2}. \end{array} $$

((48))

where C ₁, C ₂, C ₃ ranges from 4 to 16, 8 to 53, and 16 to 178, respectively.

From the previous sections, since Δ T _rh/Δ T, δ _rh/d, δ _lid/d, and λ are all scaled in terms of Ra, a, and θ with scaling exponents summarized in Table 1. τ _lid, τ _{y,c r} , and τ y,c r′ can be scaled in terms of Ra, a, and θ. The results are listed in Table 2.

Numerical results for critical yield stress and yield stress gradient: Arrhenius vs. exponential viscosity

Figures 27 and 28 show that both τ _{y,c r} and τ y,c r′ decrease with increasing θ and total viscosity contrast and might converge to asymptotic values at high viscosity contrasts, although it is difficult to tell from our limited data. To estimate the accuracy of Frank-Kamenetskii approximation in the prediction of critical yield stress, we express the ratio of yield stress for Arrhenius viscosity to that for exponential viscosity R _τ =τ _{y,c r,Arr}/τ _{y,c r,exp} and similarly for critical yield stress gradient with $R_{\tau ^{\prime }}=\tau ^{\prime }_{y,cr,\text {Arr}}/\tau ^{\prime }_{y,cr,\text {exp}}\phantom {\dot {i}\!}$ to look at the dependence of these ratios on θ (Figures 29 and 30). Both R _τ and $R_{\tau ^{\prime }}$ increase with θ and the values of yield stress for Arrhenius viscosity and that for Frank-Kamenetskii approximation get closer as the viscosity contrast increases, assuming that R _τ <1 and R τ′<1 at all θ.

Figure 27

Critical yield stress τ _{y,c r} as a function of Frank-Kamenetskii parameter θ. The Arrhenius viscosity is calculated with the non-dimensional surface temperature T ₀ and the activation energy E that gives the corresponding θ. Lower T ₀ gives a higher viscosity contrast Δ η _Arr.

Figure 28

Critical yield stress gradient τ y,c r′ as a function of θ.

Figure 29

The ratio of Arrhenius yield stress to exponential yield stress R _τ as a function of Frank-Kamenetskii parameter θ.

Figure 30

The ratio of Arrhenius yield stress gradient to exponential yield stress gradient $R_{\tau ^{\prime }}$ as a function of θ.

For the cases tested at resolution that is doubled, we find that the values for R _τ and $R_{\tau ^{\prime }}$ are within 5% difference. Therefore, 64×64a resolution is sufficient for our single-cell steady-state convection analysis.

To find the dependence of R _τ and $R_{\tau ^{\prime }}$ on the Arrhenius viscosity contrast Δ η _Arr, we consider the ratio of Arrhenius viscosity contrast to exponential viscosity contrast (Δ η _Arr/ expθ). This ratio reflects the difference between Arrhenius viscosity contrast and exponential viscosity contrast: the larger the ratio, the greater the difference. Figure 31 shows that values of the two τ _{y,c r} approach each other as the Arrhenius viscosity gets closer to the exponential viscosity. It also suggests that R _τ and $R_{\tau ^{\prime }}$ are mainly determined by the Frank-Kamenetskii parameter θ but less sensitive to Ra and a. As the difference between the Arrhenius viscosity and exponential viscosity increases, both R _τ and $R_{\tau ^{\prime }}$ appear to approach some asymptotic value.

Figure 31

Top figures: ratio of yield stress for Arrhenius viscosity to that for exponential viscosity R _τ and $R_{\tau ^{\prime }}$ as a function of Arrhenius viscosity contrast normalized to exponential viscosity contrast Δ η _Arr/ exp(θ). Δ η _Arr/ exp(θ)≥1, and it can go up to many orders of magnitude. Bottom figures: ratio of yield stress multiplied by θ. For cases with the same θ but various Ra and a, these ratios do not differ much, meaning that θ is the controlling factor for the difference in yield stress predicted by Arrhenius viscosity and exponential viscosity. Asymptotically towards high Δ η _Arr/ exp(θ), R _τ and R τ′ are approximately proportional to θ.

We also investigate the dependence of critical yield stress and yield stress gradient on aspect ratio (Figures 32 and 33). The positive power law coefficient (see Table 1) implies that with smaller aspect ratio, the critical yield value is lower and thus more difficult to reach the yielding criterion. Therefore, smaller cells are more stable. This may explain the phenomenon where overturning of the cold lid is observed once or a few times but then stabilized after reconfiguring into smaller cells.

Figure 32

Critical yield stress τ _{y,c r} as a function of aspect ratio. Circles represent Ra = 3 ×10⁷, squares Ra = 10⁷, diamonds Ra = 3 ×10⁶.

Figure 33

Critical yield stress gradient τ y,c r′ as a function of aspect ratio. Legends are the same as in Figure 32.

Figures 34 and 35 and Table 1 show that τ _{y,c r} and τ y,c r′ are approximately proportional with Ra.

Figure 34

Critical yield stress τ _{y,c r} as a function of Rayleigh number.

Figure 35

Critical yield stress gradient τ y,c r′ as a function of Rayleigh number.

The plots in Figures 27, 28, 32, 33, 34, and 35 show that the scaling exponents could have a range of values, and the cases using Arrhenius viscosity may have a bit different values from those with exponential viscosities. Figure 36 shows that while the scaling exponent for Ra and a for exponential viscosity cases lie between the range obtained from Arrhenius viscosity cases, the exponent for θ varies with surface temperature T ₀, Ra, and a. In general, the scaling exponent of θ increases with the viscosity contrast (i.e., lower T ₀), a, and moderately with Ra. This may be related to the difference in the stress distribution in the lid in Arrhenius cases, especially towards the downwelling edge, as discussed previously. In the future, it will be worthwhile to look in more detail at the stress variations for the Arrhenius viscosity and to see if there are any relationships between R _τ , $R_{\tau ^{\prime }}\phantom {\dot {i}\!}$, and T ₀ at even lower T ₀ (i.e., at higher viscosity contrasts.)

Figure 36

Scaling exponents in Ra^β a ^ζ θ ^α for τ _y (left three figures) and τ y,c r′ (right three figures) with varying T ₀ in Arrhenius viscosity. The grey stripes represent the value obtained from exponential viscosity calculations, with the width determined by error bars.

Discussion and conclusions

Comparison with other studies of stress scaling laws

The studies dealing with convective stress scaling often aim to provide an expression for stress in terms of radius and mass of a planet to predict the likelihood of plate tectonics. They reached various conclusions (e.g., O’Neill et al. 2007; O’Neill and Lenardic 2007; Valencia and O’Connell 2007;2009; Korenaga 2010a; Van Heck and Tackley 2011; Stamenkovic and Breuer 2014). One of the main difficulties in deriving convincing scaling laws for plate tectonics initiation was a poor understanding of lid stresses and how they are related to lid failure. In the present study, we have addressed these issues using two-dimensional steady-state convective cell simulations. This is the simplest system to analyze, and yet even for this system, the derivation of scaling laws proved to be complicated and not well described by the existing asymptotic theories. Below, we discuss some differences between our study and previous studies and summarize our scaling laws in a dimensional form.

In some studies (e.g., Moresi and Solomatov 1998; Trompert and Hansen 1998; Tackley 2000b; Fowler and O’Brien 2003), the authors assumed that subduction occurs when the stresses in the convective interior exceed the yield stress. This means that subduction begins when not only the lid but also the interior of the convective cell fails. However, subduction initiation may not necessarily require the failure of interiors but instead may only require failure of just a small portion of the lid. The stresses in the lid are several orders of magnitude higher than the stresses in the interior, and also, they scale differently. Thus, the assumption regarding what part of the convective cell must fail in order for subduction to begin is critically important. In this study, we have investigated this assumption quantitatively, based on a detailed analysis of stresses and other parameters in the convective cell, and then formulated the critical conditions for subduction initiation.

In agreement with Fowler (1985), we have shown that the lid slope is a key factor in determining the stresses in the lid. However, our model has several important differences from Fowler (1985). The theoretical solution in Fowler (1985) is a similarity solution and does not take into the finite horizontal extent of the lid. Our model has vertical boundaries, and thus, the structure of the lid in our model is more complex. Also, the solution in Fowler (1985) is an asymptotic solution requiring very high values of parameters, such as Ra and θ, and a satisfaction of certain asymptotic conditions, which are not reached in our simulations and may not necessarily be reached on planets. Thus, our scaling laws are not asymptotic in this sense. Also, solutions in Fowler (1985) are obtained for two end-member cases, the large lid slope case, and the small lid slope case. We find that the lid slope behaves in a more complex way and is between these two end-member cases. We have determined a scaling law for the lid slope numerically and used it to derive the scaling law for the stresses in the lid.

Our analysis suggests that the stresses in the lid increase approximately as a square of the distance from the bottom of the lid (Equation 34 and Figure 14). This agrees with the asymptotic analysis of Fowler (1985) but is different from the stress distribution in Solomatov (2004a). In Solomatov (2004a), the stress distribution was more complex because the convective cell was heated from within rather than from the bottom and the internal heating affected the temperature-induced density distribution in the lid. At Rayleigh numbers higher than those reached in Solomatov (2004a), the lid is expected to become sufficiently thin so that the heat production inside the lid would be negligible compared to the heat flux at the base of the lid. Thus, we expect that for convection with internal heating, the stress distribution in the lid should approach the quadratic distribution that we observe for convection with bottom heating.

We find that subduction initiation requires that only a part of the lid undergoes plastic failure, roughly 0.3 to 0.5 of the total lid thickness. This generally agrees with the analysis in Solomatov (2004a) and confirms that the plastic failure does not have to extend all the way to the bottom of the lid as was assumed in Fowler and O’Brien (2003). However, unlike Solomatov (2004a), we determine the distance to the boundary of the plastic failure zone by measuring it from the base of the lid and scaling it in terms of the rheological boundary layer thickness. We find that such an approach is more appropriate because the mobility of the lid is largely controlled by the viscosity contrast between the zone of failure and the convective interior of the cell, which in turn is scaled with the rheological boundary layer thickness.

Estimates for the Earth

To compare our results with those obtained in Solomatov (2004a,b), we convert the critical yield stress τ _{y,c r} and critical yield stress gradient τ y,c r′ into their dimensional forms (Equations 20 and 22) and estimate the critical yield strength and the critical coefficient of friction μ for subduction initiation on the Earth.

The interior viscosity cannot be reliably estimated from the viscosity law alone and is usually determined from better constrained properties such as lithospheric thickness. Therefore, following Solomatov (2004a), we use the scaling law for δ ₀ (Table 1) and present the results in terms of the thickness of the thermal boundary layer δ ₀∼100 km instead of the mantle viscosity η ₁.

The scaling law for the critical yield stress depends strongly on aspect ratio a. Previous studies have scaled the aspect ratio from half-space cooling of lithosphere (Korenaga 2010b; Stamenkovic and Breuer 2014) or estimated from numerical simulations (Solomatov 2004a,b), whereas it was assumed to be on the order of 1 in Valencia and O’Connell (2009). We use the horizontal width of the convective cells as l _hor=a d∼100 km as a very rough value to compare our estimates with those in Solomatov (2004a,b). This value was inferred from observational constraints on the present-day horizontal scale of sublithospheric convective structures (Solomatov 2004a).

Using the results in Table 1, we obtain that the dimensional critical yield stress for subduction initiation is as follows:

$$ { \fontsize{9}{6}\begin{aligned} \tau_{y,cr}\sim 1.95\alpha\rho g\left(\frac{E}{R{T_{i}^{2}}}\right)^{-1.03}\Delta T^{-0.03}\delta_{0}^{-0.41}l_{\text{hor}}^{1.78}d^{-0.37}. \end{aligned}} $$

((49))

For comparison, Solomatov (2004a,b) give:

$$ \tau_{y,cr}\sim 13\alpha\rho g\left(\frac{E}{R{T_{i}^{2}}}\right)^{-2}\Delta T^{-1}l_{\text{hor}}. $$

((50))

The critical coefficient of friction μ for subduction initiation is as follows:

$$ \mu\sim89\alpha\left(\frac{E}{R{T_{i}^{2}}}\right)^{-1.74}\Delta T^{-0.74}\delta_{0}^{-1.55}l_{\text{hor}}^{1.87}d^{-0.32}, $$

((51))

and in Solomatov (2004a,b):

$$ \mu\sim50\alpha\left(\frac{E}{R{T_{i}^{2}}}\right)^{-2}\Delta T^{-1}\delta_{0}^{-1}l_{\text{hor}}. $$

((52))

Using the typical values of various physical parameters (Table 3), we estimate that the yield strength for the Earth is 5 MPa which is of the same order of magnitude as 3 MPa obtained by Solomatov (2004a,b). To see how variations in various parameters may affect these estimates, it is useful to present the estimates in a different form. Our estimate (Equation 49) can be written as:

$${} {\fontsize{9.4pt}{9.6pt}\selectfont{\begin{aligned} \tau_{y,cr}\sim 5 \left(\frac{100\ \text{km}}{\delta_{0}}\right)^{0.41} \left(\frac{l_{\text{hor}}}{100\ \text{km}}\right)^{1.78} \left(\frac{500\ \text{km}}{d}\right)^{0.37}\ \text{MPa}, \end{aligned}}} $$

((53))

Table 3 Parameters used to estimate τ _y,cr and τy,cr′ for Earth as in Solomatov ( 2004a ,b)

and the estimate from Solomatov (2004a,b) (Equations 50) is as follows:

$$ \tau_{y,cr}\sim 3 \left(\frac{l_{\text{hor}}}{100\ \text{km}}\right)\ \text{MPa}, $$

((54))

from Equation 50.

Our estimates of the critical friction coefficient μ is 8 ×10⁻³, which is a factor of 3 larger than 3 ×10⁻³ obtained in Solomatov (2004a,b). Our estimate can be written as

$${} \mu\sim 0.008 \left(\frac{100\ \text{km}}{\delta_{0}}\right)^{1.55} \left(\frac{l_{\text{hor}}}{100\ \text{km}}\right)^{1.87} \left(\frac{500\ \text{km}}{d}\right)^{0.32}, $$

((55))

and the estimate from Solomatov (2004a,b) (Equation 50) is as follows:

$$ \mu\sim 0.003 \left(\frac{l_{\text{hor}}}{100\ \text{km}}\right) \left(\frac{100\ \text{km}}{\delta_{0}}\right). $$

((56))

If we take into account the fact that Frank-Kamenetskii approximation that we used to derive the scaling laws overestimate the critical yield stress and the critical friction coefficient (Figure 31), then both ours and the estimates in Solomatov (2004a,b) should be further reduced by a factor of 2 (Figure 31), depending on the values of the viscosity parameters and the Rayleigh number.

One major difference between our scaling laws and the scaling laws obtained in Solomatov (2004a) is a much stronger dependence of the critical yield stress and critical friction coefficient on the width of the convecting layer l _hor - they scale roughly as approximately $l_{\text {hor}}^{2}$ as opposed to the previous scaling approximately l _hor. This means that the critical values of the yield stress and friction coefficient would increase by 2 to 4 orders of magnitude if the width of the convective cells increased by 1 to 2 orders of magnitude (for example, in the past history of the Earth),and thus, at least in principle, could reach the experimentally observed values that are on the order of 1,000 MPa for τ _{y,c r} and μ∼0.6 to 0.85 (e.g., Byerlee 1978; Goetze and Evans 1979; Kohlstedt et al. 1995; Mei et al. 2010) and values constrained by loading models with in situ stress measurements of Hawaiian Islands (Zhong and Watts 2013) which are 0.25 to 0.7 for μ and 100 to 200 MPa for lithospheric stress. This implies that the chances of plate tectonics might be higher than we thought before. Time-dependent calculations and a more realistic formulation of the problem are required to better understand the implications of these results for plate tectonics initiation.

Uncertainties in stress scaling

The scaling laws derived here are applicable to Newtonian rheology; therefore, the activation energy for diffusion creep is used in our calculations. However, it should be noted that dislocation creep is probably the dominant mechanism in the lithosphere (Karato and Wu 1993). For the Earth, wet dislocation creep may be preferable (Solomatov and Moresi 2000) while for other terrestrial planets such as Venus might have dry lithosphere. To apply on a wider range of planets including icy bodies, scaling laws based on non-Newtonian rheology will be required.

In previous scaling theories, the lid slope is often considered to be small because the lid thickness is assumed to be relatively small. Even in the large lid slope end member case in Fowler’s theory, δ _lid is assumed to be small relative to the thickness of the convecting layer.

However, our simulation indicates that the slope may be significant, so the derivations may need to be modified to take this into account.

Free-slip boundary conditions are often used in solving equations for thermal convection, but this restricts the vertical motion of the surface. Recent studies have used the free-surface boundary conditions, which is closer to natural surface condition as both normal and shear stress on the surface is reduced to zero (Zhong et al. 1996; Schmeling et al. 2008; Kaus et al. 2010; Crameri et al. 2012; Kramer et al. 2012). It maybe computationally expensive to implement this for the time being, but it could be worthwhile to explore its effect on scaling relations for stresses in the future.

Our numerical results show that τ _{y,c r} of Arrhenius viscosity approaches that of exponential viscosity as the Frank-Kamenetskii parameter θ increases. This enables us to use exponential viscosity law to extrapolate to high Arrhenius viscosity contrast conditions. Besides the Frank-Kamenestskii approximation, the viscosity contrasts can be reduced in other ways, one of which is to set a cutoff viscosity. The stress structure resulting from the cutoff viscosity will have to be examined. We can then compare the accuracy of these approximations and apply them to extrapolate the results to planetary parameters.

Our results generally support previous conclusions that in order for the convective regime on the terrestrial planets in the inner Solar System to change from stagnant lid convection to plate tectonics, the yield stress of the lithosphere should be much smaller (several MPa) than that predicted by laboratory experiments on rock deformation (hundreds of MPa as predicted by Byerlee’s law). However, our results suggest a much stronger dependence of the critical yield stress on the horizontal width of the convective cells. This opens a possibility of subduction initiation even for the large, experimentally measured, lithospheric strength provided that a sufficiently long convective cell forms in a time-dependent mantle convection. In the future, it would be important to investigate the role of initiation conditions and statistical fluctuations of convective cells for the initiation of subduction in time-dependent convection.