Scaling laws for convection with temperature-dependent viscosity and grain-damage

Abstract

1 INTRODUCTION

Answering the major questions of geodynamics, from the thermal evolution of the Earth to the origin of plate tectonics, requires a detailed understanding of the physics of mantle convection. Mantle convection theory is best summarized in the form of scaling laws, which relate the important aspects of mantle convection, such as the heat flow and plate speed, to the key quantities governing convective circulation such as the Rayleigh number (Ra). There has been extensive work in developing scaling laws for simple cases like constant viscosity convection (e.g. Turcotte & Oxburgh 1967; McKenzie et al.1974; Sotin & Labrosse 1999), convection with temperature-dependent viscosity (e.g. Christensen 1984a; Morris & Canright 1984; Ogawa et al.1991; Davaille & Jaupart 1993; Moresi & Solomatov 1995; Solomatov 1995; Grasset & Parmentier 1998), and weakly non-Newtonian convection (e.g. Parmentier et al.1976; Parmentier & Morgan 1982; Christensen 1984b; Reese et al.1998; Solomatov & Moresi 2000). However, the rheology of the mantle is necessarily more complex than these simple cases; a strongly non-linear rheology is required in order to generate plate tectonics (e.g. Tackley 2000a; Bercovici 2003). While there have been many studies on plate generation with exotic rheologies (e.g. Weinstein & Olson 1992; Trompert & Hansen 1998; Tackley 2000b; van Heck & Tackley 2008; Foley & Becker 2009), only a few attempt to develop detailed scaling laws (e.g. Moresi & Solomatov 1998; Korenaga 2010). In particular, no study has developed scaling laws for convection with damage physics, a promising mechanism for generating plate tectonics (e.g. Bercovici et al.2001; Bercovici & Ricard 2003, 2005; Landuyt et al.2008; Landuyt & Bercovici 2009; Bercovici & Ricard 2012, 2013, 2014). Therefore, the purpose of this paper is to develop scaling laws for the heat flow, interior convective velocity, and plate velocity for convection with damage.

The damage physics we employ is a grain size feedback referred to as grain-damage (e.g. Bercovici & Ricard 2005; Landuyt et al.2008; Bercovici & Ricard 2012). Grain-damage is motivated by its effectiveness at causing weakening and shear localization in the lithosphere, its allowance of dormant weak zones, and by the geological observation that peridotitic mylonites are ubiquitous in lithospheric shear zones (White et al.1980). In particular, grain-damage's allowance for dormant weak zones means this mechanism has rheological memory, a property that other mechanisms, such as plasticity, lack (e.g. Moresi & Solomatov 1998). Such memory of past deformation is thought to be crucial for initiating subduction (Gurnis et al.2000).

Grain-damage relies on a feedback between grain size reduction and a grain size dependent viscosity; deformation reduces the grain size, making the material weaker and hence causing more deformation. Previous studies have demonstrated that grain-damage is an effective mechanism for shear localization (Landuyt et al.2008), and for producing significant lithospheric weakening (Foley et al.2012). However, neither study develops scaling laws for the heat flow or plate velocity. Scaling laws for convection with grain-damage could differ significantly from those for Newtonian or weakly non-Newtonian convection, due to the effects of damage throughout the lithosphere and mantle. We therefore perform a large suite of numerical models and use the results to derive scaling laws for convection with grain-damage from boundary layer theory.

Convection with grain-damage has not been studied extensively, so we start with the very simple case of a viscosity that depends solely on grain size, progressively adding the complexities of temperature-dependent viscosity and finally temperature-dependent grain growth in later sections. The paper is therefore organized in the following manner: grain-damage and the numerical methods employed are reviewed in Section 2; scaling laws for temperature-independent viscosity are derived and compared to numerical experiments in Section 3; scaling laws for the stagnant lid regime are derived and compared to numerical experiments with temperature and grain size sensitive viscosity in Section 4; temperature-dependent healing is added, and scaling laws for the the stagnant lid, transitional, and fully mobile lid regimes are derived and compared to numerical results in Section 5; the boundaries between the stagnant lid, transitional, and fully mobile regimes, are derived in Section 6; a curve for the onset of convection is derived in Section 7; an application of our scaling laws to Earth and Venus and other implications of this study for the thermal and tectonic evolution of planets are given in Section 8 and concluding remarks in Section 9.

2 BACKGROUND

Our grain-damage mechanism relies on the feedback between deformation induced grain size reduction and grain size dependent viscosity. This combination may seem problematic, because these processes are thought to occur via distinct microphysical mechanisms: grain size reduction nominally occurs through the propogation of dislocations in the dislocation creep regime, whose flow law is independent of grain size; meanwhile grain size sensitive flow occurs in a diffusional regime, which does not implicitly involve dislocations and grain reduction. Thus, the dislocation creep and diffusion creep deformation regimes occur in separate domains of deformation space (depending on differential stress, temperature and grain size) and therefore do not necessarily interact in a way that would cause a grain size feedback (Karato 2008). However, in a two-phase material like peridotite, deformation and damage to the interface between phases (e.g. olivine and pyroxene) combined with pinning effects allows damage, grain-reduction and diffusion creep to co-exist (Bercovici & Ricard 2012); this leads to a state of small grain permanent diffusion creep, which is observed in natural peridotitic mylonites (Warren & Hirth 2006).

The full theory for grain-damage with Zener pinning involves a two-phase material with a composite rheology, taking into account both dislocation and diffusion creep. To make the problem more numerically tractable, we use a simplified version by assuming that the permanent diffusion creep ‘pinned’ state prevails throughout the mantle. In the pinned state, the grain size of the primary phase is controlled by the curvature of the interface with the secondary phase, and thus damage to the interface leads directly to damage of the primary phase. We can therefore assume that the bulk grain size of the material is governed by the same equation as the evolution of the interface curvature [eq. 4d of Bercovici & Ricard (2012)], reducing the problem to a single phase (i.e. we no longer need to track the evolution of both the interface curvature and grain size, we only need to track the grain size). The assumption that the pinned state prevails throughout the mantle also allows us to assume grain size sensitive diffusion creep throughout the domain, and neglect non-Newtonian dislocation creep. Thus the composite rheology is reduced to a simple grain size sensitive rheology. However, in reality the rheology will be controlled by whichever mechanism allows for the easiest deformation (e.g. Rozel et al.2011); that is when grains are large dislocation creep should predominate. We discuss how the transition to dislocation creep would affect our scaling laws, and under what conditions this transition should occur, in Appendix A1.

2.1 Damage formulation

2.2 Governing equations

In addition, we define parameters to describe the variation of viscosity and healing across the mantle due to temperature dependence; |$\mu _l^{\prime } = \mu _l/\mu _m$|⁠, the viscosity ratio in the absence of damage, (with A = A_ref), and |$h_l^{\prime } = h_l/h_m$|⁠, the healing ratio, where the subscript l denotes the value in the lithosphere (i.e. at T^′ = 0). We also define the effective Frank-Kamenetskii parameter, |$\theta = E_v^{\prime }/(1+T_s^*)^2$|⁠, to describe the temperature dependence of viscosity (Korenaga 2009). The Frank-Kamenetskii parameter comes from the linear exponential viscosity law, an approximation to the full Arrhenius viscosity law, and appears in scaling laws for stagnant lid convection. Thus it is necessary to define θ for numerical experiments with the Arrhenius viscosity law in order to develop and fit scaling laws for stagnant lid convection (e.g. Korenaga 2009).

2.3 Numerical methods

We solve the coupled convection and damage equations using a 2-D Cartesian finite volume code, similar to that used in Foley et al. (2012). The code uses the SIMPLER algorithm to solve the momentum equations (Patankar 1980), employing a multigrid method for the diffusion terms in both the momentum and temperature equations. The temperature equation uses a Crank–Nicholson time discretization and the non-oscillatory version of MPDATA for the advection term (Smolarkiewicz 1984; Smolarkiewicz & Grabowski 1990). For the fineness evolution equation, we linearize the source terms using the technique laid out in Patankar (1980), and utilize a Crank–Nicholson time discretization and the non-oscillatory version of MPDATA for advection as in the temperature equation. A non-oscillatory advection scheme for fineness is essential to the stability of the numerical solution. Dispersive ripples in the fineness solution will cause ripples in the viscosity field, which can grow due to the feedback with the momentum equations, eventually causing the numerical solution to diverge. Most numerical experiments in this study were performed with a 4 × 1 aspect ratio domain, though some cases were run in larger aspect ratio domains, up to 16 × 1 (see Tables 2–5 for a compilation of all numerical results). The typical resolution used was 512 × 128, though a higher resolution of 1024 × 256 was used for models with large D/H and/or large Ra₀ (see Tables 2–5). Rerunning select cases at double the resolution typically only changes the results for the Nusselt number (Nu) by 1–3 per cent, with a maximum change of 5.5 per cent, and typically only changes the results for the plate speed by 3–5 per cent, with a maximum change of 6.5 per cent (see Appendix A2); thus the resolution used for the numerical models does not significantly impact the results.

3 TEMPERATURE-INDEPENDENT VISCOSITY

We first explore the case where mantle viscosity is sensitive to grain size only, and healing is constant (i.e. both viscosity and healing are insensitive to temperature). This simple set-up is not a good approximation of mantle convection, but allows us to understand the effects of grain-damage in isolation; temperature-dependent viscosity and healing are added in later sections (Sections 4 and 5). The numerical models generally show a convective planform where upwellings and downwellings take the form of blob-like drips from the top and bottom thermal boundary layers, while the core of convection cells are dominated by a simple horizontal shear flow between the mobile boundary layers. As seen in the scaling theory below (Section 3.1), the fineness is controlled by the ratio of damage to healing, D/H, thus we discuss model results in terms of this quantity. At low D/H, downwellings and upwellings are more sheet-like, becoming more drip-like as D/H increases; a similar trend occurs for varying Ra₀ (Fig. 1). As expected, both increasing D/H or increasing Ra₀ enhances damage in the mantle interior (i.e. the isothermal core of convection cells) and boosts the vigour of convection.

Figure 1.

Open in new tabDownload slide

Convection snapshots for a model with low D/H (D/H = 10^{− 8}, Ra₀ = 10⁶) (a), high D/H (D/H = 10^{− 4}, Ra₀ = 10⁶) (b), low Ra₀ (D/H = 10^{− 8}, Ra₀ = 10⁵) (c) and high Ra₀ (D/H = 10^{− 8}, Ra₀ = 5 × 10⁷) (d). All models use m = 2 and p = 4. Each snapshot shows the surface velocity, |$u^{\prime }_{{\rm surf}}$|⁠, the fineness field, A^′, the viscosity field, μ^′, and the temperature field, T^′.

3.1 Scaling theory

We derive scaling laws based on the idealized model in Fig. 2. We assume that convection with grain-damage behaves like constant viscosity convection, with a viscosity set by the average interior fineness of the mantle, A_i; A_i is determined by the typical stress scale in the mantle. As seen by the numerical models, the shear stress, τ_xz, which results from the horizontal shear flow between the top and bottom boundaries, is the dominant stress component driving damage in the mantle. We note here that τ_xz is not the dominant stress component throughout the convecting layer; the normal stress, τ_xx, (τ_xx = −τ_zz due to mass conservation) is generally the largest stress component in the lithosphere and in upwelling or downwelling regions. This difference is important when we develop scaling laws for plate-tectonic style convection, where we consider damage in the lithosphere (see Section 5.1). A scaling law for the convective shear stress with grain-damage is derived assuming τ_xz = (2μ_effv_l)/d, where v_l is the lithosphere (or plate) velocity, and the factor of two arises from the fact that the horizontal velocity profile goes from v_l at the surface to zero at z = d/2 (see Fig. 2). The effective interior mantle viscosity, μ_eff, is defined as μ_eff = μ_i(A_i/A₀)^−m, where μ_i is the undamaged viscosity at the average interior mantle temperature, T_i; with temperature-independent viscosity, μ_i = μ_m.

Figure 2.

Open in new tabDownload slide

Sketch of the idealized model for convection with grain-damage and temperature-independent viscosity and healing. The assumed horizontal velocity profile and resulting mantle shear stress, τ_xz, are shown (left-hand panel), as well as the vertical temperature profile (right-hand panel). The sketch also illustrates the definitions of δ_l, δ_m, v_l, v_m, T_i and A_i.

Eq. (12) demonstrates that increasing Ra₀ increases the shear stress while increasing D/H decreases shear stress. Comparing this scaling law to the shear stress scaling for uniform viscosity elucidates how grain-damage influences convection. For isoviscous convection τ_xz ∼ μ(κ/d²)Ra^2/3; this shows that τ_xz ∼ μ^1/3, and thus decreasing μ results in a net drop in shear stress. Therefore increasing D/H decreases the shear stress because higher damage reduces the effective interior mantle viscosity. Increasing Ra₀ increases shear stress, but by a smaller amount than for isoviscous convection (note that the exponent for Ra₀ in (12) is less than 2/3). The increase in strain-rate that larger Ra₀ causes by boosting velocity is somewhat offset by increased damage in the mantle.

3.2 Comparison to numerical experiments

There is generally a good agreement between the numerical results where D/H is varied (with constant Ra₀ = 10⁶, m = 2 and p = 4), and the theoretical scaling laws with C₁ = 0.15 and C₂ = 2.5 (Fig. 3). In particular the scaling laws for |$v_l^{\prime }$| (14) and Nu (19) match the numerical results well (Figs 3 a and b). In addition, |$v_{{\rm rms}}^{\prime }$| scales the same as |$v_l^{\prime }$|⁠, albeit smaller by a factor of ≈2 due to averaging over areas where velocities are low (i.e. in the core of convection cells).

Figure 3.

Open in new tabDownload slide

Plots of velocity (a), Nusselt number (b), stress (c) and fineness (d) versus D/H for temperature independent viscosity and healing (Section 3). Numerical results are plotted as symbols with theoretical scaling laws plotted as solid lines. Panel (a) shows numerical data for |$v_l^{\prime }$| compared to the scaling law for |$v_l^{\prime }$| (14); data for |$v_{{\rm rms}}^{\prime }$| is also shown. Panel (b) compares the numerical results for Nu to the scaling law for Nu (19). Panel (c) compares the data for |${\tau _{xz}^{\prime }}_{{\rm rms}}$| to the scaling law for |$\tau _{xz}^{\prime }$| (12); data for |${\tau _{xx}^{\prime }}_{{\rm rms}}$| is also shown. Panel (d) compares the data |$\bar{A}^{\prime }$| to the scaling law for |$A_i^{\prime }$| (20) and also shows data for |$A^{\prime }_{{\rm max}}$|⁠.

The scaling laws for |$\tau _{xz}^{\prime }$| (12) and |$A_i^{\prime }$| (20) fit the numerical results with a small offset (Figs 3c and d). Strictly speaking, the scaling laws for |$\tau _{xz}^{\prime }$| and |$A_i^{\prime }$| describe interior mantle values, or averages in the isothermal core of convection cells, and not whole mantle averages. However, the difference between the whole mantle average and the interior value is likely to be small, as the well-mixed interior is by far the largest region of the convecting mantle. This is confirmed by the numerical data which generally show a good agreement with the scaling laws for |$\tau _{xz}^{\prime }$| and |$A_i^{\prime }$| (Figs 3c and d); there is only a small offset between the numerical data and the scaling laws for both quantities. In addition, the numerical results show that |${\tau _{xx}^{\prime }}_{{\rm rms}}$| scales in the same manner as the shear stress, and is typically somewhat larger in magnitude. The |${\tau _{xx}^{\prime }}_{{\rm rms}}$| average is likely dominated by the boundary layers and upwelling and downwelling regions in the mantle where normal stresses are highest. Therefore the metric |${\tau _{xx}^{\prime }}_{{\rm rms}}$| does not reflect the typical value of τ_xx in the core of convection cells; τ_xz is larger than τ_xx in the convecting interior (Fig. 4b). The maximum fineness, which is confined to the boundary layers and is likely driven by the large normal stresses where downwellings and upwellings first go unstable, appears to scale differently than |$\bar{A}^{\prime }$| (Fig. 3d). As the driving force for |$A_{{\rm max}}^{\prime }$| is different than for |$\bar{A}^{\prime }$|⁠, it is not unexpected that they scale differently.

Figure 4.

Open in new tabDownload slide

Vertical profiles of horizontally averaged temperature (a), stress (b) and velocity (c) for a representative numerical result with temperature-independent viscosity and healing in statistical steady-state. The plot for stress (b) shows depth profiles of the shear stress, |${\tau ^{\prime }_{xz}}$|⁠, and normal stress, |$\tau ^{\prime }_{xx}$|⁠. The plot for velocity (C) shows depth profiles of the full velocity vector, |$v^{\prime }_{{\rm rms}}$|⁠, the horizontal velocity, u^′ and vertical velocity, w^′. All horizontal averages for stress and velocity are rms averages. Parameters for model result shown here: Ra₀ = 10⁶ and D/H = 10⁻⁸.

The numerical results for varying Ra₀ (with D/H = 10⁻⁸, m = 2 and p = 4) compare well to the theoretical scaling laws (Fig. 5). The results are similar to the experiments with varying D/H, in particular the plate velocities and Nusselt number, which are well fit by the theory. As in Fig. 3, the convective stresses and mantle fineness show small offsets between the numerical data and the scaling theory due to the way the numerical data is calculated.

Figure 5.

Open in new tabDownload slide

Plots of time averaged velocity (a), Nusselt number (b), stress (c) and fineness (d) versus Rayleigh number, Ra₀, as in Fig. 3.

4 TEMPERATURE-DEPENDENT VISCOSITY: STAGNANT LID REGIME

Adding a strongly temperature-dependent viscosity (with a constant healing rate) allows us to explore how grain-damage influences convection in the stagnant lid regime, the regime where the cold, high viscosity lithosphere no longer participates in convection (Ogawa et al.1991; Davaille & Jaupart 1993; Moresi & Solomatov 1995; Solomatov 1995). Developing scaling laws for convection with grain-damage in the stagnant lid regime is important because it allows us to establish an important baseline scenario, the fully stagnant lid end member, with which we can compare mobile and plate-like models to in later sections (see Section 5). In addition, our results for stagnant lid convection may be applicable to planetary bodies that do not exhibit plate tectonics, such as Mars or rocky and icy satellites in the solar system.

Stagnant lid convection occurs when the viscosity ratio across the top thermal boundary layer, μ_l/μ_i, reaches a critical value of about 3000, which corresponds to a viscosity ratio across the mantle, μ_l/μ_m of ≈10⁴ for bottom heated convection (Solomatov 1995). We thus use viscosity ratios at or above 10⁴ in this section (where |$\mu _l^{\prime } = 10^5$|⁠) and in Section 5. Numerical models show changes in the convective planform with D/H or Ra₀ that are similar to the temperature-independent viscosity case (Fig. 6). At low D/H or Ra₀, convection beneath the lid is sluggish with sheet-like upwellings and downwellings; we even observe the cessation of convection entirely at very low D/H or Ra₀ (Section 7). With increasing D/H or Ra₀, upwellings and downwellings become drip-like, and the lid becomes flat and relatively thin.

Figure 6.

Open in new tabDownload slide

Snapshots of stagnant lid convection (Section 4) for a model with low D/H (D/H = 10^{− 6}, Ra₀ = 10⁶) (a), high D/H (D/H = 10^{− 3}, Ra₀ = 10⁶) (b), low Ra₀ (D/H = 10^{− 5}, Ra₀ = 3 × 10⁵) (c) and high Ra₀ (D/H = 10^{− 5}, Ra₀ = 10⁷) (d). All models use m = 2, p = 4, |$E_v^{\prime } = 23.03$| and |$T_s^* = 1$|⁠. Each snapshot shows the fineness field, A^′, the viscosity field, μ^′ and the temperature field, T^′.

4.1 Scaling theory

Figure 7.

Open in new tabDownload slide

Sketch of the idealized model for convection with grain-damage in the stagnant lid regime. The locations of the immobile lid and rheological boundary layer (δ_rh) are shown in the sketch on the left-hand side, as well as the asymmetry between δ_l and δ_m. On the right-hand side, the temperature profile shows the rheological temperature scale, ΔT_rh, and the definition of T_i in the stagnant lid regime.

However, there is some disagreement over these scaling laws. Boundary layer theories and steady-state numerical results find that |$\delta _l \sim \theta {\rm Ra}_{\rm eff}^{1/5}$| (Morris & Canright 1984; Reese et al.1998); however, for time-dependent convection the steady-state boundary layer theory breaks down and (27) is the correct scaling (Solomatov & Moresi 2000; Korenaga 2009). There is also disagreement over the scaling for velocity. Solomatov & Moresi (2000) find that v ∼ (Ra_eff/θ)^1/2 for their internally heated results. However, they also state that (26) fits the bottom heated experiments of Dumoulin et al. (1999), and that the 1/2 power-law scaling may be a transitional regime found at low Rayleigh numbers. We find that (26) is the correct scaling for our results.

4.2 Comparison to numerical experiments

We compare our theory to numerical results in the stagnant lid regime where D/H and Ra₀ are varied separately (Figs 8 and 9). The theoretical curves use C₃ = 4.24, C₄ = 0.125 and a_rh = 1.3, while the numerical models varying D/H use Ra₀ = 10⁶, m = 2, p = 4, |$E_v^{\prime } = 23.03$| and |$T_s^* = 1$|⁠; those varying Ra₀ use D/H = 10⁻⁵ and the same parameters otherwise. The temperature-dependent viscosity relation used here results in an effective Frank-Kamenetskii parameter of θ = 5.76 and a viscosity ratio of |$\mu _l^{\prime } = 10^5$|⁠. The constant a_rh is determined by calculating |$\Delta T^{\prime }_{rh}$| from our numerical models, and using |$a_{rh} = \Delta T^{\prime }_{rh} \theta$| (see Section 4.1). Assuming symmetry between |$\delta ^{\prime }_m$| and |$\delta ^{\prime }_{rh}$|⁠, |$\Delta T^{\prime }_{rh} = 2(1-T_i^{\prime })$| (see Fig. 7), and the internal temperature, |$T_i^{\prime }$|⁠, is calculated from the numerical results by taking the sublid volume average of the temperature field. The numerical results consistently yield a value of a_rh ≈ 1.3. The sublid average, also used for fineness and stress, is defined as a volume average only within the convecting region (i.e. a volume average beneath the base of the lid). Computing averages only within the convecting region is important because the rigid lid has a significant effect on whole mantle averages, especially for stress which is large in the lid. We define the base of the stagnant lid as the depth where the horizontally averaged vertical velocity profile reaches a non-dimensional value of 10. As illustrated by depth profiles of horizontally averaged temperature, velocity, and stress, the definition of 〈w′〉_rms = 10 as the base of the stagnant lid is, if anything, a conservative choice (Fig. 10); that is some of the lid is included in the average. Using a different value of 〈w′〉_rms to define the start of the convecting region does not change how the quantities considered here scale with D/H and Ra₀. Finally, we calculate |$v_m^{\prime }$| from the numerical results using (21) with u^′ evaluated at z^′ = 0 instead of z^′ = 1.

Figure 8.

Open in new tabDownload slide

Plots of velocity (a), Nusselt number (b), stress (c) and fineness (d) versus D/H for convection with temperature-dependent viscosity and constant healing (Section 4). Numerical results are plotted as symbols with theoretical scaling laws plotted as solid lines. Panel (a) shows numerical data for |$v_m^{\prime }$| compared to the scaling law for |$v_m^{\prime }$| (35); data for |$v_{{\rm rms}}^{\prime }$| and |$v_l^{\prime }$| are also shown. Panel (b) compares the numerical results for Nu to the scaling law for Nu (37). Panel (c) compares the data for the sublid |${\tau _{xz}^{\prime }}_{{\rm rms}}$| to the scaling law for |$\tau _{xz}^{\prime }$| (32); data for the sublid |${\tau _{xx}^{\prime }}_{{\rm rms}}$| is also shown. Panel (d) compares the data for sublid |$\bar{A}^{\prime }$| to the scaling law for |$A_i^{\prime }$| (33) and also shows data for |$A^{\prime }_{{\rm max}}$|⁠.

Figure 9.

Open in new tabDownload slide

Plots of time averaged velocity (a), Nusselt number (b), stress (c) and fineness (d) versus Rayleigh number, Ra₀, as in Fig. 8.

Figure 10.

Open in new tabDownload slide

Vertical profiles of horizontally averaged temperature (a), fineness (b), viscosity (c), stress (d) and velocity (e) for a representative stagnant lid case in statistical steady-state. The plot for stress (d) shows depth profiles of the second invariant of the stress tensor, |$\langle \tau ^{\prime }_{{\rm II}} \rangle$|⁠, shear stress, |$\langle {\tau ^{\prime }_{xz}}\rangle _{{\rm rms}}$|⁠, and normal stress, |$\langle \tau ^{\prime }_{xx} \rangle _{{\rm rms}}$|⁠. The plot for velocity (e) shows depth profiles of the full velocity vector, 〈v′〉_rms, the horizontal velocity, 〈u′〉_rms, and vertical velocity, 〈w′〉_rms. All plots show the base of the stagnant lid, defined by 〈w′〉_rms = 10, as a horizontal black line. Parameters for model result shown here: Ra₀ = 10⁶, D/H = 10⁻⁴, |$E_v^{\prime } = 23.03$| and |$T_s^* = 1$|⁠.

The numerical results for both the experiments varying D/H and those varying Ra₀ fit the data well, similar to the results for temperature-independent viscosity. The basal mantle velocity, the sublid average shear stress, and the sublid average fineness all match the theory (Figs 8a, c, d and 9a, c, d). The Nusselt number deviates from the theory more than any of the other observables (Figs 8b and 9b), but the data appears to asymptote to the theoretical curves at high D/H and Ra₀.

5 TEMPERATURE-DEPENDENT VISCOSITY AND HEALING

With both temperature-dependent viscosity and healing, numerical convection results display a wide range of behaviour from stagnant lid convection to convection that resembles the temperature-independent viscosity case. We also observe a parameter range where convection ceases entirely; the cessation of convection is more readily observed with temperature-dependent viscosity because it occurs at higher D/H or Ra₀ (i.e. the region of parameter space where convection will not occur is larger). As shown by Landuyt & Bercovici (2009) and Foley et al. (2012), stagnant lid convection occurs when either the damage to healing ratio in the lithosphere is low, or when the Rayleigh number is low, while the cessation of convection occurs when the damage to healing ratio in the mantle is very low, or Ra₀ is very low, as is explained further in Section 7. The numerical models show that with increasing D/H, the lid becomes steadily more mobilized over a broad region of parameter space (Fig. 11). At low D/H, convection shows weak mobility with a ‘mushy lid’ as damage is just able to soften the high viscosity lithosphere. At higher D/H, more coherent downwellings can form leading to faster plate velocities; however, the plate velocity is still slower than the basal mantle velocity due to the resistance provided by the viscosity of damaged lithospheric shear zones. In addition, the wavelength of instability for the top boundary layer is large compared to the bottom boundary layer; this is a result of the long horizontal length scale required for the top boundary layer to go unstable. Finally, at large D/H, models resemble temperature-independent viscosity convection because damage is so effective at weakening the lithosphere. Convection becomes symmetric: the plate velocity converges to the basal mantle velocity and the wavelength of instability for the top boundary layer becomes equal to that of the bottom boundary layer. Similar trends occur for varying Ra₀.

Figure 11.

Open in new tabDownload slide

Snapshots of convection with temperature-dependent viscosity and healing (Section 5) in the transitional regime. Models use three different D/H values: D/H = 10⁻⁷ (a), D/H = 10⁻⁵ (b) and D/H = 10⁻³ (c). All models use Ra₀ = 10⁶, m = 2, p = 4, |$E_v^{\prime } = E_h^{\prime }=23.03$| and |$T_s^* = 1$|⁠. Each snapshot shows the surface velocity, |$u^{\prime }_{{\rm surf}}$|⁠, the fineness field, A^′, the viscosity field, μ^′ and the temperature field, T^′.

5.1 Scaling theory

The numerical results show three main regimes of convection: the fully stagnant lid regime, the transitional regime, and the fully mobile regime, with most models shown here falling in the transitional regime. The fully stagnant lid regime is defined by the assumption that convection takes place beneath a rigid lid, and that the effective viscosity of the mantle interior controls instability of drips off the base of this lid; thus the important characteristics of convection, such as |$\delta _l^{\prime }$| or |$v_m^{\prime }$|⁠, are set by the interior mantle viscosity (Section 4.1). The viscosity that effectively controls convection changes in the transitional regime, where the viscosity of damaged lithospheric shear zones determines instability of the top thermal boundary layer and thus the plate motion. In this regime, plates show a wide degree of mobility relative to the interior, from near stagnant lid behaviour when convection first enters the transitional regime, to near isoviscous convection, or full mobilization, when convection enters the fully mobile regime. Finally, the fully mobile regime is defined as the point when damage in the lithosphere is so effective that the viscosity of lithospheric shear zones no longer provides any significant resistance to flow. Thus the interior mantle viscosity again controls boundary layer thickness and convective velocity. As a result, the scaling laws for the fully stagnant lid and fully mobile regimes derived below provide end member limits to the numerical data. For the example case of varying D/H, convection follows the fully stagnant lid scaling law at low D/H, deviates from this limit toward the transitional regime scaling law at moderate D/H, and finally converges to the fully mobile limit at high D/H.

We derive new scaling laws for convection in the transitional regime, as there are no previous studies directly constraining this style of convection. Solomatov (1995) explored convection in a transitional regime at low viscosity ratios (e.g. |$\mu _l^{\prime } \approx 10^2$|–10³) that lies between constant viscosity convection and stagnant lid convection, however, it is unclear if this scaling theory is applicable to the transitional regime we observe, which sits between stagnant lid convection and fully mobile convection at high viscosity ratios.

We first derive an equation for the thickness of the top thermal boundary layer, δ_l, as scaling laws for Nu and the plate velocity are easily obtained from those for δ_l (see the conceptual sketch, Fig. 12). The scaling law for δ_l is determined semi-empirically; we use boundary layer stability analysis to constrain which non-dimensional parameters should appear in the scaling law, and empirically determine the scaling law exponents from our numerical experiments. We choose to use empirical fits because the simple boundary layer theory used in previous sections (Sections 3.1 and 4.1) fails for transitional regime convection, where plate motion is dominated by the viscosity within the boundary layer rather than in the mantle interior.

Figure 12.

Open in new tabDownload slide

Sketch of the idealized model for convection with grain-damage in the transitional regime. Sketch of the convective planform (bottom left-hand panel) shows the lithospheric thickness δ_l is controlled by the viscosity of the damaged shear zone, with a fineness A_l. The vertical temperature profile (right-hand panel) denotes T_i, and shows the asymmetry between the top and bottom boundary layers. The surface velocity (top) shows how the velocity goes from v_l to −v_l over the lithospheric shear zone of thickness ξ. This sets the lithospheric normal stress, τ_xx, which drives damage in the lithosphere.

The lithospheric normal stress, |$\tau _{xx}^{\prime }$|⁠, can be written in terms of the strain-rate in a lithospheric shear zone as |$\tau _{xx}^{\prime } \sim (\mu _l^{\prime } A^{\prime -m} v_l^{\prime })/\xi$|⁠, where ξ is the characteristic width of the shear zone (Fig. 12); there is no known scaling law for ξ. Since we assume that the lithospheric shear zone viscosity controls the dynamics of the lithosphere, the shear zone thickness must be a function of lithospheric quantities and properties, namely the lithospheric damage to healing ratio, |$D/(Hh_l^{\prime })$|⁠, the plate velocity, |$v_l^{\prime }$|⁠, and the lithospheric viscosity, |$\mu _l^{\prime }$|⁠. Therefore |$A_l^{\prime }$| is some unknown function of |$\mu _l^{\prime }$|⁠, |$v_l^{\prime }$|⁠, and |$D/(Hh_l^{\prime })$|⁠.

The plate length, L^′, arises from the convecting system, and can be calculated from boundary layer theory (e.g. Turcotte & Schubert 2002); however, as an alternative, we choose to exploit our numerical results and calculate L^′ directly from the models. We therefore treat the plate length, L^′, as an unknown in (54), similar to |$D/(Hh_l^{\prime })$| or Ra₀, and determine the influence of varying L^′ empirically from the numerical results. We compute L^′ from the numerical models using the relationship between plate speed and boundary layer thickness; this gives |$L^{\prime } = v_l^{\prime }\delta _l^{\prime 2}$|⁠, where |$\delta _l^{\prime } = T_i^{\prime }/{\rm Nu}$| from the definition of the Nusselt number, and |$v_l^{\prime }$|⁠, |$T_i^{\prime }$| and Nu are output from the numerical models (see Sections 3.2 and 4.2). The numerical results show that for a fixed aspect ratio numerical domain, L^′ is nearly constant in the transitional regime (L^′ ≈ 1.2 − 1.8, with an average value of L^′ ≈ 1.5 for the models with a 4 × 1 aspect ratio). We therefore assume that L^′ can be approximated as independent of damage to healing ratio, Rayleigh number, and viscosity ratio, and that the influence of L^′ on the top boundary layer thickness can be constrained using a set of models where the aspect ratio of the numerical domain is changed (increasing the aspect ratio of the numerical model from 4 × 1 to 16 × 1 causes L^′ to increase from ≈1.5 to ≈4). Determining the influence of plate length in this manner is analogous to how the roles of |$D/(Hh_l^{\prime })$| or Ra₀ are assessed using models where these parameters are varied.

We determine the scaling law exponents empirically. We first perform a least squares fit to numerical experiments where Ra₀ is varied, and find β_Ra ≈ −2/3 for m = 2 and p = 4. Using β_Ra = −2/3, we perform least squares fits to determine β_D ≈ −1/3, β_μ ≈ 1/4 and β_L ≈ 1/10 (see Table 1 for the scaling law constants with different combinations of m and p); all three fits give C₅ ≈ 20 for m = 2 and p = 4 (Fig. 13). Our constraint on β_L is only based on a limited number of numerical experiments. However, given that the influence of L^′ is significantly less than the influence of damage to healing ratio, Rayleigh number, and viscosity ratio, and that the plausible range of variation in L^′ for planetary mantle convection is also significantly less than these other parameters, which can vary many orders of magnitude, our simple model to incorporate the effect of varying plate length is sufficient for constraining the first order scaling laws.

Figure 13.

Open in new tabDownload slide

Plots of numerical data with least-squares fits for varying Ra₀ (a), |$D/(Hh_l^{\prime })$| (b), |$\mu _l^{\prime }$| (c) and L^′ (d). The internal temperature, |$T_i^{\prime }$|⁠, varies as Ra₀, |$D/(Hh_l^{\prime })$|⁠, |$\mu _l^{\prime }$|⁠, or L^′ are varied. We therefore take the variation of |$T_i^{\prime }$| into account by fitting |${\rm Ra}_0T_i^{\prime }$| to the data for |$\delta _l^{\prime }$| (a), and fitting |$D/(Hh_l^{\prime })$|⁠, |$\mu _l^{\prime }$|⁠, and L^′ to |$\delta _l^{\prime }T_i^{\prime 2/3}$|⁠. The results shown here use m = 2, p = 4, and: (a) D/H = 10⁻⁵, |$\mu _l^{\prime }=10^6$|⁠, |$h_l^{\prime } = 10^{-5}$|⁠, L^′ = 1.5 (for a 4 × 1 aspect ratio numerical domain); (b) Ra₀ = 10⁶, |$\mu _l^{\prime }=10^5$|⁠, L^′ = 1.5; (c) Ra₀ = 10⁶, D/H = 10⁻⁵, |$h_l^{\prime } = 10^{-5}$|⁠, L^′ = 1.5 and (d) Ra₀ = 10⁶, D/H = 10⁻⁵, |$h_l^{\prime } = 10^{-5}$|⁠, |$\mu _l^{\prime } = 10^5$|⁠. The results in (d) use numerical domains with aspect ratios ranging from 4 × 1 to 16 × 1.

The scaling laws for |$v_m^{\prime }$| and |$A_i^{\prime }$| have the same form as the fully stagnant lid and fully mobile laws, because they are controlled by the viscosity of the interior mantle. They only differ from the fully stagnant lid and fully mobile scaling laws by the internal temperature, which determines both the buoyancy of the bottom boundary layer and the viscosity and healing rates in the mantle. The scaling laws for Nu and |$v_l^{\prime }$| have entirely different forms than the stagnant lid and fully mobile laws, as the driving forces for damage in lithospheric shear zones are dominated by in-plane normal stresses, which are different from those in the interior of convection cells where shear stresses dominate. The scaling laws for these lithospheric quantities are stronger functions of D/H and Ra₀ than either end member limit, and will thus make the transition from the fully stagnant lid limit to the fully mobile limit as either D/H or Ra₀ increase.

5.2 Comparison to numerical experiments

The theoretical scaling laws for convection with temperature-dependent viscosity and healing compare well to the numerical data for |$v_l^{\prime }$|⁠, Nu, |$v_m^{\prime }$| and |$A^{\prime }_i$| (see Figs 14 –17). In addition, the fully stagnant lid and fully mobile limits (shown in Fig. 14) illustrate how convection results evolve through the transitional regime and converge to the fully mobile limit as D/H increases. To plot the fully stagnant lid scaling laws, we calculate a_rh from stagnant lid numerical results, in the same manner as outlined in Section 4.2, and find a_rh ≈ 1.82. All results shown here use m = 2, p = 4, |$E_h^{\prime } = 23.03$| which results in |$h_l^{\prime } = 10^{-5}$| , and, save for those shown in Fig. 17, a 4 × 1 aspect ratio domain resulting in L^′ ≈ 1.5.