# Stagnant-lid convection: comparison of viscosity laws and uniform scaling approach for temperature and heat flux prediction

Cécile Grigne

► To cite this version:

Cécile Grigne. Stagnant-lid convection: comparison of viscosity laws and uniform scaling approach for temperature and heat flux prediction. Geophysical Journal International, 2023, 235, pp.2410 - 2429. 10.1093/gji/ggad375. hal-05269518

HAL Id: hal-05269518

https://hal.univ-brest.fr/hal-05269518v1

Submitted on 19 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Stagnant-lid convection: comparison of viscosity laws and uniform scaling approach for temperature and heat flux prediction

Cécile Grigné*

Geo-Ocean, Univ. Brest, CNRS, Ifremer, UMR6538, F-29280 Plouzane, France. E-mail: cecile.grigne@univ-brest.fr

Accepted 2023 September 18. Received 2023 July 25; in original form 2023 March 8

## SUMMARY

Numerical simulations of stagnant-lid convection in a 2-D Cartesian fluid heated from below are carried out in order to study how the temperature dependence of the viscosity affects the vertical profile of temperature in the model. I test different viscosity laws, including the Arrhenius law with realistic parameters for the Earth’s mantle or for icy bodies. No approximation is made, which leads to extremely high viscosity contrasts. Results are compared to different approximations, in particular the Frank–Kamenetskii (FK) one. I propose a new approach for the scaling of the temperature drop across the convective part of the layer beneath the stagnant lid. The vertical profile of temperature as a function of the viscosity law is predicted, with a uniform scaling approach for all formulations of the temperature-dependent viscosity. The predicted profiles are in very good agreement with results of 2-D numerical simulations in Cartesian geometry. The complete scaling given here provides a rapid way to compare viscosity laws and to check how approximations affect the results, in terms of interior temperature, stagnant lid thickness and heat flux, compared to the real Arrhenius law for rocky mantles and for the icy outer shells of satellites. In particular, in the context of 2-D Cartesian convection heated from below, in the stagnant-lid regime, I propose a new approach to properly scale the FK formulation when it is used as an approximation of the Arrhenius law.

**Key words:** Heat flow; Rheology; mantle; Stagnant-lid convection; Scaling laws; Viscosity.

## 1 INTRODUCTION

The mantle of terrestrial planets and rocky bodies has a viscosity that strongly depends on temperature. The outer ice I layer of large icy satellites also has a rheology that depends on temperature. For a Newtonian rheology, that is to say for a linear relation between strain rate and stress, and with a viscosity that depends only on temperature, thermal convection can be in a so-called stagnant-lid regime: the upper cold part of the layer is so viscous that it is effectively rigid (e.g. Nataf & Richter 1982; Richter et al. 1983; Morris & Canright 1984; Fowler 1985; Ogawa et al. 1991; Davaille & Jaupart 1993; Moresi & Solomatov 1995). Only the warmer and less viscous medium underneath is mobile. This mode of convection is exhibited by telluric planets and rocky satellites other than the Earth, and probably also in the outer shells of icy moons (Deschamps & Sotin 2001; Yao et al. 2014; Hussmann et al. 2015). Scaling laws for the heat flux as a function of the Rayleigh number in this stagnant-lid regime have been proposed (e.g. Christensen 1984; Moresi & Solomatov 1995; Dumoulin et al. 1999; Deschamps & Sotin 2000). These scaling laws obtained with a Cartesian geometry

and without internal heating have been used as the base for more complex cases, for example with internal heating, in a spherical geometry and/or with a variable thermal conductivity (Yao et al. 2014; Guerrero et al. 2018; Thiriet et al. 2019; Harel et al. 2020; Deschamps 2021; Deschamps & Vilella 2021).

The temperature dependence of the viscosity of rocky mantles and the outer shells of icy moons can be described by an Arrhenius law (Kirby 1983; Karato & Wu 1993; Weertman 1983):

$$ \eta_{\text{Arr}}(T) = \eta_{\text{ref}} \exp \left( \frac{E_a}{R} \left( \frac{1}{T} - \frac{1}{T_{\text{ref}}} \right) \right), \quad (1) $$

where $E_a$ is the activation energy, $R = 8.314 \text{ J(K.mol)}^{-1}$ is the gas constant derived from the Boltzmann distribution, and $\eta_{\text{ref}}$ is a reference viscosity at the reference temperature $T_{\text{ref}}$. For estimated terrestrial values of $E_a$ (a few hundreds of kJmol$^{-1}$, Kirby 1983; Karato & Wu 1993) and for the estimated temperature drop across the Earth’s mantle for instance, the ratio between the surface and interior viscosity is larger than $10^{50}$, which is very challenging for numerical simulations (Stein & Hansen 2013). Therefore, a linearized form for the temperature term in the exponential part of eq. (1) is often used (e.g. Morris & Canright 1984; Fowler 1985; Moresi & Solomatov 1995; Solomatov & Moresi 1996; Grasset

<footnote>*Now at: Laboratoire Geo-Ocean, Technopôle Brest-Iroise, rue Dumont d’Urville, F-29280 Plouzané, France.</footnote>
& Parmentier 1998; Dumoulin et al. 1999). As is done by previ-
ous authors, it will be referred to as the ‘Frank–Kamenetskii (FK)
approximation’ hereafter:

$$
\eta_{\text{FK}}(T) = \eta_{\text{ref}} \exp\left(\theta (T_{\text{ref}} - T)\right), \quad (2)
$$

where θ is a parameter (in K⁻¹) that controls the temperature de-
pendence of the viscosity. For a total temperature drop ΔT = T_b -
T_s, where T_b and T_s are the bottom and surface temperatures, the
logarithm of the ratio between the maximum and minimum viscos-
ity ln(η(T_s)/η(T_b)), denoted hereafter by γ and referred to as the
'total viscosity contrast', is

$$
\gamma_{\text{Arr}} = \frac{E_a}{R} \frac{\Delta T}{T_b T_s} \qquad (3)
$$

for the Arrhenius viscosity law (eq. 1), and

$$
\gamma_{\text{FK}} = \theta \Delta T \qquad (4)
$$

for the FK approximation (eq. 2). The same total viscosity contrast ($\gamma_{\text{Arr}} = \gamma_{\text{FK}}$) is therefore obtained for $\theta = E_a/(R T_b T_s)$.

Because the FK approximation is mostly used to prevent high
viscosity contrasts that are problematic for numerical simulations,
the parameter θ is often chosen so that the total viscosity contrast
remains moderate. In the stagnant-lid regime, most of the viscosity
variations are accommodated by the stagnant lid in 2-D Carte-
sian geometry (Solomatov & Moresi 1996). Below the lid, in the
well-mixed interior, the variations of temperature, and therefore of
viscosity, are less important. The value of θ may be chosen so that
the variations of viscosity within the convective part of the model
are close to the ones that would be obtained with the real Arrhenius
law. For instance, to obtain the same viscosity gradient with the
two laws at the base of the model (for T = T<sub>b</sub>), θ must be cho-
sen as θ = E<sub>a</sub>/(R T<sub>b</sub><sup>2</sup>) (Korenaga 2009; Ferrick & Korenaga 2023).
Denoting T<sub>i</sub> the temperature of the well-mixed convective interior,
θ = E<sub>a</sub>/(R T<sub>i</sub><sup>2</sup>) can also be used (e.g. Solomatov & Moresi 1996;
Reese et al. 1999; King 2009; Harel et al. 2020). However for this
definition of θ, the value of T<sub>i</sub> may not be known a priori.

Comparisons between the results obtained with the real Arrhenius
law and the FK approximation were done in the stagnant-lid regime
and/or with a non-Newtonian rheology by Reese et al. (1999), King
(2009), Korenaga (2009), Stein & Hansen (2013) and Harel et al.
(2020). Most studies focus on the heat flux, written as the dimen-
sionless Nusselt number Nu, as a function of the Rayleigh number
Ra of the fluid. With similar viscosity gradients through the con-
 vective part of the model and for a similar interior temperature T<sub>i</sub>,
the Nu–Ra relationship is shown to differ significantly between the
Arrhenius law and the FK approximation, with a heat flux about
20 per cent higher with the FK approximation (Reese et al. 1999;
Harel et al. 2020).

Few studies focus on the scaling for the interior temperature
Tᵢ as a function of the parameters that control the temperature
dependence of the viscosity. Davaille & Jaupart (1993) proposed
that the temperature drop across the active thermal boundary layer
(hereafter denoted as TBL), below the stagnant lid, is proportional
to a viscous temperature scale ΔTᵥ defined as

$$
\Delta T_v = - \frac{\eta}{d\eta/dT}. \qquad (5)
$$

This reasoning was henceforth extensively used in studies about
the relation Nu–Ra in the stagnant-lid regime (e.g. Dumoulin et al.
1999; Reese et al. 1999; Solomatov & Moresi 2000; Deschamps
& Sotin 2000; Korenaga 2009; Deschamps & Lin 2014; Thiriet

et al. 2019; Harel et al. 2020; Deschamps 2021). With the FK approximation (eq. 2), the temperature scale is ΔT<sub>v</sub> = 1/θ and does not depend on T. One can then simply replace the temperature drop across the active TBL with θ<sup>-1</sup>. Given the temperature drop T<sub>b</sub> – T<sub>i</sub> across the bottom TBL, T<sub>i</sub> can be predicted from T<sub>b</sub> – T<sub>i</sub> ~ θ<sup>-1</sup>, which fits satisfyingly numerical simulations in the stagnant-lid regime (Deschamps & Sotin 2000; Deschamps & Lin 2014).

With the real Arrhenius law (eq. 1), the temperature scale ΔT<sub>v</sub>
depends on T. To my knowledge, no proper tests were carried out
so far to verify that the aforementioned approach can be extended
to all viscosity laws. In this study, I show that the scaling relying
on the viscous temperature scale ΔT<sub>v</sub> is actually not appropriate for
the Arrhenius law with a high viscosity contrast in 2-D Cartesian
geometry. The objective of this paper is to compare the effect of the
chosen formulation for the temperature dependence of the viscosity
on the thermal state of the layer, in the context of stagnant-lid
convection only. I do not address the conditions controlling the
transition from a mobile surface to a stagnant one. It was shown
by Yao et al. (2014) and Guerrero et al. (2018) that the viscosity
contrasts required to obtain stagnant-lid convection depended on
the geometry of the models, with an important influence of the
curvature of the system in 2-D spherical annulus and 3-D spherical
geometry. It must be kept in mind that the results presented here are
valid only for 2-D Cartesian geometry, for an incompressible fluid
heated from below.

First, I present the numerical model used in this study, with the
different viscosity laws that have been tested. Results of numerical
simulations are then discussed in comparison to previous proposed
scaling laws for the interior temperature beneath the stagnant lid,
namely predictions arising from eq. (5). I then propose a different
approach to predict the vertical temperature profile for stagnant-lid
convection. Finally, I present possible applications for this method.
In particular, I indicate how to properly scale parameters for the FK
approximation so that the same temperature and viscosity profiles
are obtained as with the Arrhenius law.

2 MODEL

2.1 Physical model

The governing equations for basally heated thermal convection at
infinite Prandtl number for a Newtonian fluid in the Boussinesq approximation are (i) mass conservation:

$$
\vec{\nabla} \cdot \vec{v} = 0, \tag{6}
$$

(ii) momentum conservation:

$$
\vec{0} = -\vec{\nabla} P + \vec{\nabla} \cdot (\eta (\vec{\nabla} \vec{v} + \vec{\nabla} \vec{v}^T)) + \rho_b \alpha (T - T_b) \vec{g} \quad (7)
$$

and (iii) heat conservation:

$$
\frac{\partial T}{\partial t} + \vec{v} \cdot \vec{\nabla} T = \kappa \nabla^2 T, \qquad (8)
$$

where $\vec{v}$ is the velocity vector, $P$ is the non-hydrostatic pressure, $\eta$ is the dynamic viscosity, $T$ is the temperature, $\rho_b$ is the reference density for the reference temperature $T_b$ (here defined at the bottom of the model), $\alpha$ is the thermal expansion, $\vec{g}$ the gravity acceleration and $\kappa$ the thermal diffusivity. Physical parameters with their symbols are in Table 1, along with their values when they are fixed in this study. I carry out a first set of numerical experiments with parameters appropriate for the Earth’s mantle. Parameters for
<table><caption>**Table 1.** Parameters used in this study, with the range of tested values for varying parameters. The layer is considered as incompressible and the total temperature difference $\Delta T = T_b - T_s$ is the superadiabatic jump between the bottom and the surface of the model.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Layer thickness</td><td><i>D</i></td><td>2890 km</td></tr><tr><td>Gravity acceleration</td><td><i>g</i></td><td>10 m s<sup>-2</sup></td></tr><tr><td>Thermal expansion</td><td><i>α</i></td><td>3 × 10<sup>-5</sup> K<sup>-1</sup></td></tr><tr><td>Thermal diffusivity</td><td><i>κ</i></td><td>10<sup>-6</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Reference density</td><td><i>ρ</i><sub>b</sub></td><td>3300 kg m<sup>-3</sup></td></tr><tr><td>Reference (bottom) temperature</td><td><i>T</i><sub>ref</sub> = <i>T</i><sub>b</sub></td><td>2000 K</td></tr><tr><td>Surface temperature</td><td><i>T</i><sub>s</sub></td><td>300 K</td></tr><tr><td>Reference (bottom) viscosity</td><td><i>η</i><sub>ref</sub> = <i>η</i><sub>b</sub></td><td>10<sup>20</sup> to 10<sup>22</sup> Pa.s</td></tr><tr><td>Reference (bottom) Rayleigh number</td><td><i>Ra</i><sub>ref</sub> = <i>Ra</i><sub>b</sub></td><td>3.985 × 10<sup>6</sup> to 3.985 × 10<sup>8</sup></td></tr><tr><td>Activation energy</td><td><i>E</i><sub>a</sub></td><td>23–1413 kJ mol<sup>-1</sup></td></tr><tr><td>Total viscosity contrast</td><td><i>γ</i></td><td>8–100</td></tr></tbody></table>

convection in the outer shell of icy satellites will be tested in Section 4.3. But the dynamics of the system, be it a rocky mantle or an ice I layer, is controlled by dimensionless parameters. A large range of parameters is covered in this study and results that are obtained here for Earth's mantle are actually also relevant for icy moons (see Section 4.3). A 2-D Cartesian geometry is used, with free-slip and isothermal conditions at the top and bottom, and with periodic boundary conditions along the vertical walls.

## 2.2 Numerical method: SFS2D code

Mass (eq. 6) and momentum conservation (eq. 7) are solved using the new code SFS2D (Stream Function Solver in 2-D). Equations are written in a dimensionless form (e.g. Schubert et al. 2001), introducing the reference Rayleigh number at the bottom of the layer:

$$Ra_b = \frac{\alpha \rho_b g (T_b - T_s) D^3}{\kappa \eta_b} \quad (9)$$

(see Table 1 for symbols). The horizontal velocity $u$ and the vertical velocity $w$ are written as functions of the stream function $\varphi$: $u = -\partial_y \varphi$ and $w = \partial_x \varphi$. This ensures that mass conservation (eq. 6) is satisfied. For a uniform viscosity, the conservation of momentum (eq. 7) can be written as a simple biharmonic equation for the stream function $\varphi$ (e.g. McKenzie et al. 1973; Jarvis 1984; Olson 1987):

$$\nabla^4 \varphi + Ra_b \frac{\partial T}{\partial x} = 0, \quad (10)$$

which can be easily solved through two successive Poisson equations (Hewitt et al. 1980). For a variable viscosity, the biharmonic equation is

$$\left(\frac{\partial}{\partial x^2} - \frac{\partial}{\partial y^2}\right) \left(\eta \left(\frac{\partial^2 \varphi}{\partial x^2} - \frac{\partial^2 \varphi}{\partial y^2}\right)\right) + 4 \frac{\partial^2}{\partial x \partial y} \left(\eta \frac{\partial^2 \varphi}{\partial x \partial y}\right) + Ra_b \frac{\partial T}{\partial x} = 0 \quad (11)$$

(Schubert et al. 2001; Nakakuki et al. 2008; Gerya 2019). This equation is solved with a finite-difference approach on a staggered grid, with the stream function $\varphi$ defined at cells' vertices. Other parameters are positioned on the staggered grid as is done in the code StagYY (Tackley 2008) (temperature and viscosity at cells' centres, velocities in the middle of cells' faces). Eq. (11) requires a 13-points stencil in 2-D. Multigrid cycles are used, with either Jacobi under-relaxations or Gauss-Seidel over-relaxations, the latter yielding a slightly faster convergence. Solving eq. (11) for the stream function avoids convergence problems that are often encountered for

viscosity contrasts higher than 20 (i.e. an increase of $10^8$ between bottom and surface viscosity) when solving conjointly pressure and velocity (eqs 6 and 7) through corrective iterations (e.g. Tackley 2008).

For heat conservation (eq. 8), advection-diffusion of temperature is solved with an alternating-direction implicit method (Peaceman & Rachford 1955), using a total variation diminishing (TVD) scheme for the advective terms. Various flux limiters (Roe 1986) are implemented in SFS2D and were tested, with no significant differences observed in the results (see the Appendix). The minmod scheme was used in the rest of the study.

The SFS2D code has been benchmarked against the data sets from Tosi et al. (2015) (see the Appendix), showing a good agreement with other codes.

## 2.3 Viscosity as a function of temperature

A Newtonian rheology is used, with a viscosity $\eta$ varying only as a function of the temperature $T$, according either to the Arrhenius law (eq. 1) or to the FK approximation (eq. 2). The total temperature drop $\Delta T = T_b - T_s$ is used to write temperatures in a dimensionless form: $T^* = (T - T_s)/\Delta T$, with the asterisk indicating dimensionless values. The viscosity $\eta^*$ is written as $\eta^* = \eta/\eta(T_b)$, which corresponds to $\eta^* = 1$ for $T = T_b$, at the bottom of the model where $T^* = 1$. Eq. (1) becomes:

$$\eta_{\text{Arr}}^*(T^*) = \exp\left(E_a^* \left(\frac{1}{T^* + T_{\text{off}}^*} - \frac{1}{1 + T_{\text{off}}^*}\right)\right), \quad (12)$$

where $E_a^* = E_a/(R\Delta T)$ is the dimensionless activation energy, and $T_{\text{off}}^* = T_s/\Delta T$ is the surface temperature normalized by $\Delta T$ if eq. (1) is directly used. To avoid too strong viscosity contrasts, $T_{\text{off}}^*$ has sometimes been chosen equal to 1 (e.g. Tackley 1996; Foley & Becker 2009; Rolf et al. 2017; Langemeyer et al. 2021), while for a surface temperature $T_s = 300$ K and a total temperature difference $\Delta T$ around 1500 K, it should be approximately 0.2. To study how strongly this affects the results in the stagnant-lid regime, two Arrhenius models are presented hereafter: ‘real Arrhenius’ with $T_{\text{off}}^* = T_s/\Delta T$, and ‘Arrhenius with offset’, for $T_{\text{off}}^* = 1$. Other values of the offset are studied in Section 4.2.

For the FK approximation, eq. (2) is non-dimensionalized as:

$$\eta_{\text{FK}}^*(T^*) = \exp(\theta^* (1 - T^*)) \quad (13)$$

where $\theta^* = \theta \Delta T$.

SFS2D can deal with extremely large viscosity contrasts in the stagnant-lid regime. Realistic values of $E_a$ for the mantle [a few
hundreds of kJ mol⁻¹, Karato & Wu (1993) are tested, leading to a total viscosity contrast γ up to 100 [i.e. η(T_s) = 2.7 × 10⁴³ η(T_b)]. The purpose here is to compare very different viscosity laws at high viscosity contrasts, in order to propose a uniformly defined scaling for the interior temperature T_i. Therefore, as a first approach, I do not try to find a way of writing θ (FK approximation) as a function of E_a (Arrhenius law) that would be appropriate in the convective part of the model, but rather consider models that give that same total viscosity contrast γ = ln(η(T_s)/η(T_b)) = ln(η*(0)).

## 2.4 Range of viscosity parameters

A total of 263 simulations were carried out: 60 in the FK approximation, 128 with the real Arrhenius law and 75 with the Arrhenius law with offset. Complete tables of numerical results are given in the Supplementary Material. The bottom (reference) viscosity is varied from 10²⁰ to 10²² Pa.s, which corresponds to a Rayleigh number varying from 3.985 × 10⁶ to 3.985 × 10⁸ with the fixed parameters given in Table 1. Models are 2-D in Cartesian geometry with an aspect ratio 2×1, with periodic boundary conditions along the vertical walls. A 256 × 128 regular grid is used for all the models, except for cases with a bottom viscosity equal to 10²⁰ Pa.s (bottom Rayleigh number equal to 3.985 × 10⁸), for which a 512 × 256 grid is used. The total viscosity contrast γ = ln(η*(0)) with the FK approximation ranges from γ = 8 to 35, but not all Rayleigh numbers were tested for γ ≥ 25, as cases where the stagnant lid occupies more than 60 per cent of the total thickness of the layer were not studied. In order to obtain the same range of interior temperature with the real Arrhenius law compared to the FK approximation, the activation energy E_a is varied from 23.5 to 293 kJ mol⁻¹, which corresponds to viscosity contrast going from 8 to 100. With the Arrhenius law with offset, in order to have total viscosity contrasts going from 8 to 50, E_a is varied from 226 to 1413 kJ mol⁻¹.

# 3 RESULTS AND SCALING LAWS

## 3.1 Interior temperature

The interior temperature Tᵢ* is estimated using either one temperature field when the model is in perfect steady state, or stacking together the temperature fields of at least 2000 time steps once statistical steady state is reached. Then, on this stacked temperature file, the horizontally averaged profile of temperature (see Fig. 1) is considered and Tᵢ* corresponds to the purely vertical part of the profile. Quantitatively, to obtain the observed value of Tᵢ*, I extract the median value over the part of the profile where the dimensionless thermal gradient is less than 0.05. To illustrate this, Fig. 1 presents six cases, with the FK approximation (eq. 13), the real Arrhenius law (eq. 12 with Tᵢ* = 0.1765), and the Arrhenius law with offset (eq. 12 with Tᵢ* = 1) and the corresponding vertical temperature profiles and values of Tᵢ*. The two columns in Fig. 1 are (left) for cases with the highest tested reference viscosity η_b = 10²² Pa.s, and therefore the lowest Rayleigh number Ra_b, and a moderate viscosity contrast (γ = 15) and (right) for the lowest reference viscosity η_b = 10²⁰ Pa.s and a high viscosity contrast (γ = 30). Cases in the left-hand column (a1–c1) are in steady state, while the ones in the right-hand column (a2–c2) are time-dependent. The grey lines in the vertical profiles of temperature are the predicted profiles, obtained with the method presented hereafter in Section 3.3.

As expected, the total viscosity contrast γ is not the parameter controlling the interior temperature Tᵢ*: in Fig. 1, two values of γ

are presented (γ = 15 in the left-hand column and γ = 30 in the right one), and the three viscosity laws for the same value of γ clearly give different thermal profiles. When using the FK approach as an approximation to the real Arrhenius law, it was shown (e.g. Solomatov & Moresi 1996; Reese et al. 1999; King 2009) that the parameter θ should be written as θ = E_a/(R Tᵢ²), i.e. in a dimensionless form θ* = E_a*/(Tᵢ* + Tᵢ*²)², and not as the value that gives the same total viscosity contrast: θ* = E_a*/(Tᵢ* (1 + Tᵢ*)). The appropriate way to scale θ* will be discussed in Section 4.4.

Deschamps & Sotin (2000) proposed that the temperature jump across the bottom thermal boundary layer (TBL) ΔTᵢ*bot = 1 - Tᵢ* is proportional to the viscous temperature scale (eq. 5). This is also implicitly used by Harel et al. (2020) when they fit the Nusselt number as a function of Raᵇ θ*ᶜ, where b and c are two exponents. In the FK approximation (eq. 13), ΔTᵢ*bot ~ ΔTᵢ*ν can directly be written

$$1 - T_i^* = \frac{a}{\theta^*} \qquad (14)$$

where a is a constant. For the Arrhenius law (eq. 12), when writing ΔTᵢ*ν for the temperature Tᵢ*, this becomes:

$$1 - T_i^* = a \frac{(T_i^* + T_{\text{off}}^*)^2}{E_a^*} \qquad (15)$$

In Fig. 2, ΔTᵢ*bot = 1 - Tᵢ* is plotted against ΔTᵢ*. The purpose here is to check whether eqs (14) and (15) can predict observed values of Tᵢ*. In Fig. 2 and in the rest of the study, I separate cases in the stagnant-lid regime from cases where the surface is partially deformed, with the following criterion: the RMS of the strain rate at the surface ε̇*surf*ms must be less than 5 per cent of the RMS of the velocity V*ms within the model. Cases for which ε̇*surf*ms > 0.05 V*ms are denoted as ‘non-stagnant’ hereafter. When considering only cases in the stagnant-lid regime, relation (14) is well observed for the FK approximation (Fig. 2a), with a slope a close to 1.1. But the linear relation in eq. (15) for the Arrhenius law (Figs 2b and c) is clearly not appropriate. For the Arrhenius law, a new approach to predict Tᵢ* must be proposed.

## 3.2 Scaling law for the temperature

I consider a convective loop model (Fig. 3(a)), as is done by Turcotte & Oxburgh (1967), Turcotte & Schubert (2002) and Grigné et al. (2005). I denote by uᵣₐₜ and uᵣₐ₉ the horizontal speeds at the bottom and top of the convective layer (i.e. uᵣₐ₉ is the horizontal speed below the stagnant lid), and by wₕ and wₐ₇ the vertical speeds of the hot upwelling and cold downwelling, respectively. The buoyancy forces of the vertical plumes are driven by temperature differences δTₕ/ₐ₉ (the subscripts h and c stand for ‘hot’ or ‘cold’), that are considered proportional to the temperature drops across the active thermal boundary layers, shown by their thicknesses δtₐ₉ and δtᵣₐ₉ in Fig. 3: δTₕ ~ Tₐ₇ - Tᵢ for the hot upwelling, and δTₐ₇ ~ Tᵢ - Tₐ₉ for the cold downwelling. Hereafter, the top and bottom thermal boundary layers (TBL) of the convective part of the model (i.e. not including the stagnant lid) will be named ‘top TBL’ and ‘bottom TBL’.

In this loop model, I do not consider the convective part below the stagnant lid as isoviscous, in contrast to previous authors who based their scaling on this assumption (e.g. Morris & Canright 1984; Fowler 1985; Solomatov 1995; Solomatov & Moresi 1996; Deschamps & Lin 2014; Yao et al. 2014; Thiriet et al. 2019). Indeed, when examining the vertical velocities of the hot and cold plumes,
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Snapshots of the temperature field and corresponding vertical profile of temperature, with the obtained temperature $T_i^*$ indicated in boxes above. Left-hand panels (a1, b1 and c1): reference viscosity is $\eta_b = 10^{22}$ Pa.s, corresponding to $Ra_b = 4 \times 10^6$, and viscosity contrast is $\gamma = 15$. These three cases are steady-state. Right-hand panels (a2, b2, and c2): reference viscosity is $\eta_b = 10^{20}$ Pa.s ($Ra_b = 4 \times 10^8$) and viscosity contrast is $\gamma = 30$. First line (a1 and a2) is for the FK approximation (eq. 13), second line (b1 and b2) for the real Arrhenius law (eq. 12 with $T_{off}^* = 0.1765$), and third line (c1 and c2) for the Arrhenius law with offset (eq. 12 with $T_{off}^* = 1$). Light grey lines in the vertical profiles of temperature are the predicted ones, with the method presented in Section 3.3.</figcaption></figure>

a non-negligible asymmetry appears, that may be attributed to different viscosities controlling the behaviours of the top and bottom TBLs, from where the plumes originate. I did not precisely analyse the vertical velocities within the plumes, which would require a tracking method to follow mobile plumes when the dynamics is not steady-state, but I studied the maximum upwards and downwards vertical speeds within the model. These can be seen as approximate measurements of the dimensionless values of $w_h$ (upwards hot plumes) and $w_c$ (downwards cold plumes). Fig. 4 shows that $w_h^*$ is almost always larger than $w_c^*$. Values of the ratio $w_h^*/w_c^*$ for all the experiments done here are given in the tables of numerical results in the Supplementary Material. For simplicity, the notations $w_h^*$ and $w_c^*$ are used in Fig. 4, but the values measured here (maximum upwards and downwards speeds) are only rough estimates of the real plume velocities when experiments are not steady-state. For comparison, the ratio between the maximum upwards and downwards speeds was also checked for isoviscous convection, with $Ra_b$ going from $10^3$ to $10^8$, and it was always equal to $1 \pm 0.005$. Fig. 4 shows

the clear asymmetry between cold and hot plumes. The scaling I propose hereafter relies on the fact that the heat flux through the top colder and more viscous TBL must be equal to the heat flux through the bottom hotter and less viscous TBL. Only certain combinations of $T_c$ and $T_i$ as a function of the viscosity parameters are then possible.

The thickness $\delta_h$ of the hot upwelling is approximately the same as $\delta_{bot}$ (hot plumes are initiated in the bottom boundary layer), and in the same way $\delta_c \sim \delta_{top}$. The amplitudes of the buoyancy forces of the plumes are therefore written as:

$$ \begin{cases} f_c \sim \alpha \rho g (T_i - T_c) \delta_{top} \\ f_h \sim \alpha \rho g (T_b - T_i) \delta_{bot} \end{cases} \qquad (16) $$

For the shear stresses experienced by the horizontal top active boundary layer and by the cold descending plume, which are both within the temperature range [$T_c, T_i$], I consider the highest possible viscous resistance, that is to say the highest viscosity $\eta(T_c)$. Similarly for the shear stresses on the bottom boundary layer and on the hot
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Bottom TBL temperature jump $\Delta T_{\text{bot}}^* = 1 - T_i^*$ as a function of the viscous temperature scale $\Delta T_v^*$, for the three viscosity laws tested here. The legend indicates the symbols used for real stagnant cases as opposed to non-stagnant ones, with different shapes and colours for the varying bottom Rayleigh number $Ra_b$. Other parameters are given in Table 1.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Loop model for stagnant-lid convection (a) and associated vertical temperature profile (b). D is the total thickness, and l is the thickness of the stagnant lid. $T_i$ is the interior temperature and $T_c$ the temperature at the base of the stagnant lid. The thicknesses of the TBL are denoted by $\delta_{\text{top}}$ and $\delta_{\text{bot}}$. The convective part of the model, below the stagnant lid, is not considered as isoviscous: the active top TBL is colder and more viscous than the bottom TBL, which implies $\delta_{\text{top}} > \delta_{\text{bot}}$ and $u_{\text{top}} < u_{\text{bot}}$. The cold plume originating in the top TBL is also thicker and slower than the hot plume originating in the bottom TBL: $w_c < w_h$. The grey line in the vertical temperature profile (b) is the conductive profile, shown here to indicate the definition of $\delta_{\text{bot}}^{\text{max}}$ (eq. 33).</figcaption></figure>

ascending plume, within the temperature range [$T_i$, $T_b$], I use the viscosity $\eta(T_i)$. Denoting by $\tau_{w_c/l}$ the shear stress on the cold/hot plumes and by $\tau_{u_{\text{top}}/\text{bot}}$ the shear stresses on the horizontal boundary layers, the following estimates can be written:

$$ \left\{ \begin{aligned} \tau_{w_c} &= \eta(T_c) \frac{w_c}{D/2} \\ \tau_{w_h} &= \eta(T_i) \frac{w_h}{D/2} \end{aligned} \right. \qquad (17) $$

and

$$ \left\{ \begin{aligned} \tau_{u_{\text{top}}} &= \eta(T_c) \frac{u_{\text{top}}}{(D-l)/2} \\ \tau_{u_{\text{bot}}} &= \eta(T_i) \frac{u_{\text{bot}}}{(D-l)/2} \end{aligned} \right. \qquad (18) $$

Mass conservation imposes

$$ \left\{ \begin{aligned} u_{\text{top}}(D-l) &= w_c D \\ u_{\text{bot}}(D-l) &= w_h D. \end{aligned} \right. \qquad (19) $$

Note that the use of $D/2$ as the length scale for the velocity change for both the upwelling and the downwelling, while allowing $w_h$ and $w_c$ to be slightly different, implies that mass conservation is not strictly respected here with the simplified velocity profiles of Fig. 3(a) [same for the difference between $u_{\text{top}}$ and $u_{\text{bot}}$, with the same length scale $(D-l)/2$].

The top and bottom heat flux across the boundary layers $\delta_{\text{top}}$ and $\delta_{\text{bot}}$ are written using either the half-space cooling model (Turcotte
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Observed dimensionless maximum upwards speed (considered as an estimate of the hot plume speed $w_h^*$) as a function of the maximum downwards speed (considered as an estimate of the cold plume speed $w_c^*$). The different symbols are for the different viscosity laws.</figcaption></figure>

& Oxburgh 1969) or a linear profile of temperature:

$$ \begin{cases} q_{\text{top}} = 2k (T_i - T_c) \left( \frac{u_{\text{top}}}{\pi \kappa D} \right)^{1/2} = k \frac{T_i - T_c}{\delta_{\text{top}}} \\ q_{\text{bot}} = 2k (T_b - T_i) \left( \frac{u_{\text{bot}}}{\pi \kappa D} \right)^{1/2} = k \frac{T_b - T_i}{\delta_{\text{bot}}} \end{cases} \quad (20) $$

First, to estimate the temperature $T_i$ and $T_c$, the local balance for momentum conservation for the two plumes is written as $f_c \sim \tau_{w_c}$ and $f_h \sim \tau_{w_h}$ (eqs 16 and 17). Using mass conservation (eq. 19) and noting that eq. (20) imply $u_{\text{top}} \sim \kappa D/\delta_{\text{top}}^2$ and $u_{\text{bot}} \sim \kappa D/\delta_{\text{bot}}^2$ (e.g. Turcotte & Oxburgh 1967; Solomatov 1995; Turcotte & Schubert 2002), it can be shown that

$$ \left\{ \begin{aligned} \alpha \rho g (T_i - T_c) \delta_{\text{top}}^3 &\sim 2\kappa \eta(T_c) (D-l)/D \\ \alpha \rho g (T_b - T_i) \delta_{\text{bot}}^3 &\sim 2\kappa \eta(T_i) (D-l)/D. \end{aligned} \right. \quad (21) $$

The equality of the top and bottom heat flux (eqs 20) implies $\delta_{\text{top}}/\delta_{\text{bot}} = (T_i - T_c)/(T_b - T_i)$. The ratio of the two lines in eqs (21) then gives

$$ \frac{T_i - T_c}{T_b - T_i} = \left( \frac{\eta(T_c)}{\eta(T_i)} \right)^{1/4} \quad (22) $$

This equation, which is key to the scalings presented hereafter, arises from the fact that the top and bottom TBL do not have the same velocities and thicknesses, with a colder, more viscous and therefore slower and thicker top TBL compared to the bottom one. The heat fluxes through the two TBLs, however, must be the same once thermal equilibrium is reached. From this reasoning, as appears in eq. (22), the temperature drop through the top TBL has to be larger than the one through the bottom TBL.

For the Arrhenius law (eq. 12) and in a dimensionless form, this becomes:

$$ \frac{T_i^* - T_c^*}{1 - T_i^*} = \exp\left(\frac{E_a^*}{4}\left(\frac{1}{T_c^* + T_{\text{off}}^*} - \frac{1}{T_i^* + T_{\text{off}}^*}\right)\right) \quad (23) $$

For the FK approximation (eq. 13), eq. (22) is

$$ \frac{T_i^* - T_c^*}{1 - T_i^*} = \exp\left(\frac{\theta^*}{4}\left(T_i^* - T_c^*\right)\right) \quad (24) $$

I hereafter use the notations $\Delta T_{\text{top}}^* = T_i^* - T_c^*$ and $\Delta T_{\text{bot}}^* = 1 - T_i^*$. For the FK approximation, eq. (24) becomes $\Delta T_{\text{bot}}^* = \Delta T_{\text{top}}^* \exp(-\theta^* \Delta T_{\text{top}}^*/4)$, showing that $\Delta T_{\text{bot}}^*$ can easily be computed as a function of $\Delta T_{\text{top}}^*$. For the Arrhenius law (eq. 23), in contrast, $T_i^*$ cannot be isolated as a function of $T_c^*$ or of $\Delta T_{\text{top}}^* = T_i^* - T_c^*$, but a numerical solution can be found.

Computed values of $\Delta T_{\text{bot}}^*$ as a function of $\Delta T_{\text{top}}^*$ using either eq. (23) or eq. (24) are plotted in Fig. 5, for the three viscosity laws tested here. For the FK approximation (Fig. 5a), the value of $\theta^*$ is varied from 10 to 40. For the real Arrhenius law (Fig. 5b), $E_a$ is varied from 100 to 400 kJ mol$^{-1}$, and for the Arrhenius law with offset (Fig. 5c), unrealistic values of $E_a$ (up to 1000 kJ mol$^{-1}$) are tested in order to reach similar interior temperatures as for the other laws. Fig. 5 shows that for a given value of $\theta^*$ (FK approximation) or $E_a$ (Arrhenius law), there is a maximum value of $\Delta T_{\text{bot}}^* = 1 - T_i^*$ (i.e. a minimum value for the interior temperature $T_i^*$) to satisfy eq. (22). This value of $\Delta T_{\text{bot}}^*$, denoted as $\Delta T_{\text{bot}}^{\text{max}}$ hereafter, decreases with the viscosity contrast, that is with increasing values of $\theta^*$ or $E_a$.

For the bottom TBL, combining eq. (20) and eq. (21) to eliminate $\delta_{\text{bot}}$, and using $Ra_b$ defined by eq. (9), it can be shown that in a
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Values of ($T_i^*$, $T_c^*$) respecting eq. (22), represented as $\Delta T_{\text{bot}}^* = 1 - T_i^*$ as a function of $\Delta T_{\text{top}}^* = T_i^* - T_c^*$, for (a) the FK approximation with different values of $\theta^*$, for (b) the real Arrhenius law and for (c) Arrhenius with the offset $T_{\text{off}}^* = 1$, with varying $E_a$. The same scale is used for the three frames, and values of $\theta^*$ and $E_a$ are chosen to exhibit similar temperatures for the three laws (which implies very high values of $E_a$ for the 'Arrhenius with offset' case). The perfectly isoviscous case is represented by the black dashed line $\Delta T_{\text{bot}}^* = \Delta T_{\text{top}}^*$. The maximum value of $\Delta T_{\text{bot}}^*$ is indicated in each case by a coloured circle.</figcaption></figure>

dimensionless form:

$$q_{\text{bot}}^* \sim \text{Ra}_b^{1/3} \frac{(1 - T_i^*)^{4/3}}{\eta^*(T_i^*)^{1/3}}. \qquad (25)$$

For each viscosity laws $\eta^*(T_i^*)$, one can derive the value of $\Delta T_{\text{bot}}^* = 1 - T_i^*$ for which the maximum heat flux $q_{\text{bot}}^*$ is obtained. It is always larger than $\Delta T_{\text{bot}}^{\text{max}}$. Therefore, among the possible values of $\Delta T_{\text{bot}}^* = 1 - T_i^*$ that can satisfy eq. (22) (Fig. 5), $\Delta T_{\text{bot}}^{\text{max}}$ is the one yielding the most efficient heat transfer. One could propose that the interior temperature $T_i^*$ is then such that $\Delta T_{\text{bot}}^*$ corresponds to this maximum heat flux, but observed values of $\Delta T_{\text{bot}}^*$ are systematically smaller than $\Delta T_{\text{bot}}^{\text{max}}$. In Fig. 5, the dashed line $\Delta T_{\text{top}}^* = \Delta T_{\text{bot}}^*$ corresponds to an isoviscous convective interior below the stagnant lid, that is to say a perfect symmetry between the top and bottom TBL, or equivalently, between the cold downwards and hot upwards plumes. Fig. 4 shows that there is actually a non negligible asymmetry in the convective cell below the stagnant lid, but one can suggest that this imbalance between hot and cold plumes should not be too strong for mass to be conserved within the convective loop (see Fig. 3). I propose that the observed temperature $T_i^*$ is such that $\Delta T_{\text{bot}}^*$ lies between the isoviscous case (exact balance between the top and bottom TBL: $\Delta T_{\text{top}}^* = \Delta T_{\text{bot}}^*$) and the value that yields the highest possible heat flux, that is to say, the maximum value of $\Delta T_{\text{bot}}^*$. This is written as:

$$1 - T_i^* = a \Delta T_{\text{bot}}^{\text{max}} \qquad (26)$$

where a is constant that will be determined numerically.

Fig. 6 shows $\Delta T_{\text{bot}}^* = 1 - T_i^*$ as a function of $\Delta T_{\text{bot}}^{\text{max}}$. For the three viscosity laws, the relation is indeed approximately linear when only real stagnant cases are considered (symbols with black edges in Fig. 6). Furthermore, the slope a is almost the same for the three viscosity laws. The best fit was computed separately for the three laws, as well as for all the cases together, each time for real stagnant cases only. The best-fitting slope for all cases together is 0.714, and it is 0.714, 0.713 and 0.692 for the FK approximation, the real Arrhenius law and the Arrhenius law with offset, respectively (see Fig. 6).

The following method can thus be applied to predict the interior temperature $T_i^*$:

(i) Find the maximum bottom TBL temperature difference $\Delta T_{\text{bot}}^{\text{max}} = 1 - T_i^{\text{min}}$ that allows eq. (22) to be satisfied by solving eq. (23) or eq. (24).

(ii) Compute $T_i^*$ from $1 - T_i^* = a \Delta T_{\text{bot}}^{\text{max}}$, with $a = 0.714$ (I here choose to use a unified scaling for all viscosity laws).

For the FK approximation, eq. (24) can be rewritten as

$$\Delta T_{\text{bot}}^* = \Delta T_{\text{top}}^* \exp\left(-\frac{\theta^*}{4} \Delta T_{\text{top}}^*\right), \qquad (27)$$

which can be used to search for the value of $\Delta T_{\text{top}}^*$ that gives $\Delta T_{\text{bot}}^* = \Delta T_{\text{bot}}^{\text{max}}$. It is such that

$$\left(\frac{d\Delta T_{\text{bot}}^*}{d\Delta T_{\text{top}}^*}\right) = \left(1 - \frac{\theta^*}{4} \Delta T_{\text{top}}^*\right) \exp\left(-\frac{\theta^*}{4} \Delta T_{\text{top}}^*\right) = 0, \qquad (28)$$

that is to say $\Delta T_{\text{top}}^* = 4/\theta^*$, which yields $\Delta T_{\text{bot}}^{\text{max}} = 4e^{-1/\theta^*}$, and therefore

$$\Delta T_{\text{bot}}^* = 1 - T_i^* = a \Delta T_{\text{bot}}^{\text{max}} = \frac{4 e^{-1} a}{\theta^*} \approx \frac{1.051}{\theta^*} \qquad (29)$$

when using $a = 0.714$. Many previous authors (e.g. Morris & Canright 1984; Davaille & Jaupart 1993; Solomatov 1995; Solomatov & Moresi 1996; Dumoulin et al. 1999; Reese et al. 1999; Solomatov & Moresi 2000; Deschamps & Sotin 2000; Korenaga 2009; Deschamps & Lin 2014; Thiriet et al. 2019; Harel et al. 2020; Deschamps 2021) have considered that the temperature jumps across the active TBLs for stagnant-lid convection were proportional to a viscous temperature scale (eq. 5), which is in the case of the FK approximation is proportional to $1/\theta^*$ (eq. 14). The approach proposed here, though derived from a different physical reasoning, corresponds to the same scaling for the FK approximation.

However, for the Arrhenius law, it was shown (Fig. 2) that the temperature jump $\Delta T_{\text{bot}}^*$ across the bottom TBL was not proportional to the viscous temperature scale (Fig. 2). With the new approach proposed here, $\Delta T_{\text{bot}}^*$ cannot be solved for analytically, but a numerical solution can easily be obtained as a function of the dimensionless activation energy $E_a^*$ and of the temperature $T_{\text{off}}^* = T_s/\Delta T$, using for instance Newton-Raphson method.

### 3.3 Thermal profile and heat flux

Now that $T_i^*$ can be predicted from the parameters controlling the temperature-dependent viscosity, an equation of the form of eq. (25) can be used to propose a scaling for the dimensionless heat flux $\text{Nu}$. Plotting the observed heat flux as a function of the right-hand side
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Bottom TBL temperature jump $\Delta T_{\text{bot}}^{\text{max}} = 1 - T_{\text{bot}}^{*}$ as a function of the maximum temperature jump $\Delta T_{\text{bot}}^{\text{max}}$, for the three viscosity laws. The best fit line is plotted as a coloured thick line for each viscosity law, with the coefficient a given in the box.</figcaption></figure>

(RHS) of eq. (25) does not show a perfect fit (Fig. 7). I infer that the aspect ratio $D/(D-l)$ of the convective cell (see Fig. 3), which depends on the thickness $l$ of the stagnant lid, cannot be neglected for the scaling of the heat flux (Turcotte & Schubert 2002; Grigné et al. 2005).

As is done by Grigné et al. (2005), the balance between the work of buoyancy forces and the work of shear stresses is written as

$$ (f_c w_c + f_h w_h)(D-l) = (\tau_{w_c} w_c + \tau_{w_h} w_h)(D-l) + (\tau_{u_{\text{top}}} u_{\text{top}} + \tau_{u_{\text{bot}}} u_{\text{bot}}) D. \quad (30) $$

Using the expressions of $f_{c/h}$, $\tau_{w_c/h}$ and $\tau_{u_{top}/bot}$ given by eqs (16), (17) and (18), and eqs (19) and (20) to express of the velocities $w_{c/h}$ and $u_{top/bot}$ as function of temperature jumps and of TBL thicknesses, one can show that eq. (30) can be written

$$ 2 \frac{\alpha \rho g}{\pi \kappa} \left( \frac{T_i - T_c}{\delta_{\text{top}}} + \frac{T_b - T_i}{\delta_{\text{bot}}} \right) = \left( \frac{D-l}{D} + \left( \frac{D}{D-l} \right)^3 \right) \left( \frac{\eta(T_c)}{\delta_{\text{top}}^4} + \frac{\eta(T_i)}{\delta_{\text{bot}}^4} \right). \quad (31) $$

The equality of the top and bottom thermal gradients $(T_i - T_c)/\delta_{\text{top}} = (T_b - T_i)/\delta_{\text{bot}}$, combined with eq. (22) to show that $\eta(T_c)/\delta_{\text{top}}^4 = \eta(T_i)/\delta_{\text{bot}}^4$, gives the following expression for the thickness $\delta_{\text{bot}}$ of the bottom boundary layer:

$$ \frac{1}{\delta_{\text{bot}}} = \left( \frac{2 \alpha \rho g (T_b - T_i)}{\pi \kappa \eta(T_i)} \right)^{1/3} \left( \frac{D-l}{D} + \left( \frac{D}{D-l} \right)^3 \right)^{-1/3} . \quad (32) $$

Fig. 7 shows that the best fits (coloured lines) obtained over the complete range of Nusselt numbers Nu, for the simple relation given by eq. (25), systematically underpredicts low values of Nu. This is especially visible in the inset in Fig. 7, for Nu values smaller than 4. One possible explanation is that the TBL thickness $\delta_{\text{bot}}$ predicted by an equation of the form of eq. (32) is too large, yielding a too low heat flux $q_{\text{bot}} = k(T_b - T_i)/\delta_{\text{bot}}$. With high viscosity contrasts, very low values of Nu $= q_{\text{bot}}^*$ can be obtained even for well developed convection, at high values of the Rayleigh number $Ra_b$, because the stagnant lid is very thick. For instance, case (a2) in Fig. 1, with the FK approximation with $\theta^* = 30$ and $Ra_b = 3.99 \times 10^8$, has a dimensionless heat flux Nu = 2.81. For comparison, the isoviscous case for this value of $Ra_b$, with the same boundary conditions and geometry, yields Nu = 125. In stagnant-lid convection, values of Nu can be close to the minimum value Nu = 1 obtained for conduction,

even for vigorous convection below the lid. For isoviscous convection, a simple relation of the form $Nu \sim Ra_b^{1/3}$ is not appropriate for low values of Nu, and I infer that in the stagnant-lid regime, some corrective terms might also be needed to account for low values of Nu.

When considering a loop model (Fig. 3), some important simplifications are made, among which considering that thermal boundary layers are well described by their mean thicknesses $\delta_{\text{top/bot}}$, proportional to $u_{\text{top/bot}}^{-1/2}$. As mentioned above and visible in Fig. 7, this possibly means an overprediction of $\delta_{\text{bot}}$ (i.e. underprediction of Nu). I suggest that there exists a maximum value of $\delta_{\text{bot}}$ at the onset of convection, when the Rayleigh number is just above the critical one, and that $\delta_{\text{bot}}$ should always be lower than this value $\delta_{\text{bot}}^{\text{max}}$. Indeed, using $\delta_{\text{bot}} \sim u_{\text{bot}}^{-1/2}$, as is done here to obtain eq. (32), implies that $\delta_{\text{bot}} \to \infty$ when $u_{\text{bot}}$ is close to zero, while for isoviscous convection, the maximum possible thicknesses of TBLs is half the thickness of the whole layer. For stagnant-lid convection, at the onset of convection when the thermal vertical profile slightly deviates from the conductive one (grey line in Fig. 3b), motion is possible only in a warm enough part of the whole layer, where the viscosity is low enough. Considering this ‘warm enough’ temperature to be close to $T_i^*$, this implies that

$$ \frac{T_b - T_i}{\delta_{\text{bot}}^{\text{max}}} = \frac{T_b - T_s}{D} \quad (33) $$

(see Fig. 3b for the vertical gradient of temperature along the conductive profile). To ensure that $\delta_{\text{bot}}$ is not overestimated for cases with a low heat flux, I write an effective TBL thickness $\delta_{\text{bot}}^{\text{eff}}$ as the harmonic mean between this maximum value and the value estimated through the loop model (eq. 32):

$$ \frac{1}{\delta_{\text{bot}}^{\text{eff}}} = \frac{1}{\delta_{\text{bot}}^{\text{max}}} + \frac{1}{\delta_{\text{bot}}} = \frac{1}{D} \left( \frac{T_b - T_s}{T_b - T_i} \right) + \frac{1}{\delta_{\text{bot}}} \quad (34) $$

I then use $\delta_{\text{bot}}^{\text{eff}}$ to write the prediction of the bottom heat flux $q_{\text{bot}}$ as

$$ \begin{aligned} q_{\text{bot}} &= k \frac{T_b - T_i}{\delta_{\text{bot}}^{\text{eff}}} \\ &= k \frac{T_b - T_s}{D} + k (T_b - T_i)^{4/3} \left( \frac{2 \alpha \rho g}{\pi \kappa \eta(T_i)} \right)^{1/3} \left( \frac{D-l}{D} + \left( \frac{D}{D-l} \right)^3 \right)^{-1/3} . \end{aligned} \quad (35) $$
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Observed dimensionless heat flux (Nusselt number Nu) as a function of the right-hand side (RHS) of eq. (25), with a numerically obtained prefactor c: Nu<sub>obs</sub> = c Ra<sub>b</sub><sup>1/3</sup>(1 - T<sub>i</sub><sup>*</sup>)<sup>4/3</sup> η<sup>*</sup>(T<sub>i</sub><sup>*</sup>)<sup>-1/3</sup>. Linear regressions show the obtained prefactor c for the best fit, for each viscosity law and for all cases together.</figcaption></figure>

Using the notations $\Delta T = T_b - T_s$, $l^* = l/D$ for the dimensionless
stagnant lid thickness, $\eta^* = \eta/\eta_b$ for the dimensionless viscosity,
and $T_i^* = (T_i - T_s)/\Delta T$ for the interior temperature, this can be
written as

$$
q_{\text{bot}} - k \frac{\Delta T}{D} = k \frac{\Delta T}{D} \left(\frac{2}{\pi}\right)^{1/3} \frac{(1 - T_i^*)^{4/3}}{\eta^*(T_i^*)^{1/3}} \left(\frac{\alpha \rho g \Delta T D^3}{\kappa \eta_b}\right)^{1/3} \left(1 - l^* + \frac{1}{(1 - l^*)^3}\right)^{-1/3} . \quad (36)
$$

With the definition of Ra<sub>b</sub> given by eq. (9), the dimensionless heat
flux Nu = q<sub>bot</sub><sup>*</sup> = q<sub>bot</sub>/(k ΔT/D) is then such that:

$$
\text{Nu} - 1 = c \text{Ra}_b^{1/3} \frac{(1 - T_i^*)^{4/3}}{\eta^*(T_i^*)^{1/3}} \left(1 - l^* + \frac{1}{(1 - l^*)^3}\right)^{-1/3}, \quad (37)
$$

where the term $(2/\pi)^{1/3}$ has been replaced by the prefactor c to take
into account the fact that the convective cell model used here (see
Fig. 3) is a simplification of the real geometry of the flow. This
prefactor c will be determined numerically. To do so, values of the
stagnant lid thickness $l^*$ are needed in the RHS of eq. (37). Some
measurements of $l^*$ from numerical simulations could be drawn, but
some criterion must then be chosen. To check for a fully predictive
scaling, where the only inputs are Ra<sub>b</sub> and the viscosity parameters
controlling $\eta^*(T_i^*)$, and the only output is Nu, I prefer to use the
following approach:

(i) Compute the predicted values of the interior temperature $T_i^{*\text{pred}}$ and of the temperature below the stagnant lid $T_c^{*\text{pred}}$ from the viscosity parameters, as is described in Section 3.2.

(ii) Compute the corresponding stagnant lid thickness $l^{*\text{pred}}$, using
the fact that the surface heat flux can be written as $q_{\text{top}} = k(T_c$

$-T_s)/l$. Since $q_{\text{top}} = q_{\text{bot}}$, this writes in a dimensionless form as
$\text{Nu} = T_c^{*\text{pred}}/l^{*\text{pred}}$. Here, for this first step conducted to evaluate
the prefactor $c$, Nu will be its observed value $\text{Nu}^{\text{obs}}$, so that the
stagnant lid thickness is evaluated as

$$
l^{*''\text{pred''}} = \frac{T_c^{*\text{pred}}}{\text{Nu}^{\text{obs}}} \quad (38)
$$

(iii) Evaluate the prefactor c with comparing

$$
(Nu - 1)_{\nu^{\text{pred}}} = c \text{Ra}_b^{1/3} \frac{(1 - T_i^{*\text{pred}})^{4/3}}{\eta^*(T_i^{*\text{pred}})^{1/3}} \left( 1 - \frac{T_c^{*\text{pred}}}{\text{Nu}^{\text{obs}}} + \frac{1}{\left(1 - \frac{T_c^{*\text{pred}}}{\text{Nu}^{\text{obs}}}\right)^3} \right)^{-1/3} \quad (39)
$$

to (Nu − 1)<sub>obs</sub> (which uses only the observed value of Nu). Note
that (Nu − 1)<sub>pred</sub> here is not yet a fully predicted value, since it
uses Nu<sup>obs</sup>. But once the prefactor c is known, eq. (39) can be used
to numerically estimate Nu directly from Ra<sub>b</sub> and from the viscosity
parameters only.

Comparison of (Nu − 1)<sub>obs</sub> and (Nu − 1)<sub>pred</sub> is done in Fig. 8,
showing that the prefactor c is close 0.5 for the three viscosity laws
tested here. The unified value c = 0.489 will be used hereafter. Now
that c is determined, eq. (39) can be used to numerically solve Nu
(making it effectively Nu<sup>pred</sup>). For simpler computation, eq. (39) is
written in terms of the stagnant lid thickness l<sup>*pred</sup> and not Nu<sup>pred</sup>
as:

$$
\left(\frac{T_c^*}{l^*} - 1\right)^3 = c^3 \text{Ra}_b \frac{(1 - T_i^*)^4}{\eta^*(T_i^*)} \left(1 - l^* + \frac{1}{(1 - l^*)^3}\right)^{-1}, \quad (40)
$$
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Observed values of the Nusselt number, written as $(Nu - 1)_{obs}$, as a function of the semi-predicted values $(Nu - 1)_{pred}$ (computed with eq. 39). Linear regressions show the estimates for the prefactor $c$ in eq. (39) (only stagnant cases are considered for these linear regressions).</figcaption></figure>

where the superscript ‘pred’ was dropped for simplicity. Isolation of $l^*$ cannot be achieved from eq. (40), but a numerical solution can be found once $T_i^*$ and $T_c^*$ are predicted from the approach given in Section 3.2.

This approach is checked in Fig. 9, for the three viscosity laws tested in this study (three columns), and for different combinations of the Rayleigh number $Ra_b$ and of the total viscosity contrast $\gamma$. The predicted thermal profile is given with black segments. The upper and lower thermal gradients are computed using the predicted values of $T_i^*$ and $Nu = T_c^*/l^*$. The horizontally averaged vertical profiles of temperature obtained with the simulations are given as thick coloured curves, showing a very good fit between the predicted and observed profiles, except of course for the ‘corners’ at the transitions between TBLs and the well-mixed convective part of the model. A wide range of parameters is presented in Fig. 9, with the interior temperature $T_i^*$ going from 0.84 to 0.96, and for thin and thick stagnant lids ($l^*$ going from 0.029 to 0.55). Additional cases can be seen in Fig. 1, where the predicted vertical temperature profile is given in grey. Note that the definition for the base of the ‘stagnant lid’ here is simply the depth at which the temperature is $T_i^*$ (see Fig. 3 and Section 3.2). It is only a thermal definition, useful for the scaling approach, but it does not rely on a direct measurement based on the velocity field as was done by previous authors (e.g. Davaille & Jaupart 1993; Deschamps & Lin 2014).

## 3.4 Summary of the complete scaling approach

The interior temperature, the thickness of the stagnant lid and the heat flux can be predicted with a uniform approach, valid for the different expressions of the temperature dependence of the viscosity.

The method proposed here cannot be expressed as simple equations, but numerical resolution is possible and offers a quick way to predict the vertical profile of temperature in the stagnant-lid regime, for incompressible convection heated from below and in 2-D Cartesian geometry.

The method is summarized here:

(i) With denoting by $T_c^*$ the temperature below the stagnant lid and by $T_i^*$ the interior temperature in the well mixed part of the layer, eq. (22) written in a dimensionless form is:

$$ \frac{T_i^* - T_c^*}{1 - T_i^*} = \left( \frac{\eta^*(T_c^*)}{\eta^*(T_i^*)} \right)^{1/4} \quad (41) $$

It is used to constrain the possible combinations $(T_c^*, T_i^*)$ that give an equal top and bottom heat flux.

(ii) These combinations are such that there is a maximum value of $\Delta T_{bot}^* = 1 - T_i^*$, and the observed value of $T_i^*$ is such that:

$$ 1 - T_i^* = a \Delta T_{bot}^{\text{max}} \quad (42) $$

with $a = 0.714$ for the 2-D Cartesian models analysed here. With eq. (41), this gives a prediction of both $T_i^*$ and $T_c^*$.

Note that for the FK approximation, eq. (42) can be simplified to $1 - T_i^* = 1.051/\theta^*$ (eq. 29), and corresponds to the usual scaling of temperature using a viscous temperature scale $\Delta T_\nu$ (eq. 5).

(iii) Using the reference viscosity (in my case, at the bottom of the model for temperature $T_b$) to compute the reference Rayleigh number $Ra_b$, the thickness of the stagnant lid $l^*$ can be deduced from eq. (40), rewritten here as

$$ \left( \frac{T_c^*}{l^*} - 1 \right)^3 \left( 1 - l^* + \frac{1}{(1 - l^*)^3} \right) = c^3 Ra_b \frac{(1 - T_i^*)^4}{\eta^*(T_i^*)} \quad (43) $$
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Predicted (black segments) and obtained (coloured curves) for six cases. Each row is for each viscosity law tested here: 'FK approximation' (eq. 13), 'real Arrhenius' (eq. 12 with $T_{off}^* = T_s/\Delta T$) and 'Arrhenius with offset' (eq. 12 with $T_{off}^* = 1$). Each column is for a different combination of bottom Rayleigh number $Ra_b$ and total viscosity contrast $\gamma$. The values of the predicted and observed interior temperatures $T_i^*$ are given in a coloured box in each frame, as well as the predicted and observed Nusselt numbers $Nu$, with the relative difference given in percent under each box. The thickness $l^*$ of the stagnant lid is also given.</figcaption></figure>
using $c = 0.489$.

(iv) The heat flux, written as the dimensionless Nusselt number Nu, is then computed as

$$Nu = T_c^* / l^* \qquad (44)$$

and the thicknesses of the top and bottom TBL can be deduced:

$$\left\{ \begin{array}{l} \delta_{\text{top}}^* = (T_i^* - T_c^*) / \text{Nu} \\ \delta_{\text{bot}}^* = (1 - T_i^*) / \text{Nu} \end{array} \right. \qquad (45)$$

A Python script that solves eqs (41)–(45), considering an Arrhenius viscosity law or the FK approximation, written either with dimension or in a dimensionless form, is given as Supporting Information. The limitation of this approach is that it is valid only for real stagnant lid cases, where the surface strain rate can be neglected.

## 4 POSSIBLE APPLICATIONS IN NUMERICAL MODELLING AND DISCUSSION

This study provides a method to predict the vertical profile of temperature corresponding to a given set of parameters for stagnant-lid convection in 2-D Cartesian geometry, for an incompressible layer heated from below. Possible applications include the easy comparison of formulations for the temperature dependence of the viscosity. If numerical approximations are necessary in a context where the actual Arrhenius law should be used, the proposed scaling offers a rapid way to check how much the thermal profile will deviate from the expected one in the stagnant-lid regime. Hereunder, I first discuss two common approximation methods used to avoid too high viscosity contrasts: a limitation through a viscosity cut-off and the choice of a non-physical value for $T_{\text{off}}^*$. I then apply the scaling approach of Section 3.4 to the ice I outer shell of a generic icy moon, to check that the scaling approach is valid in a context different from Earth’s mantle. Finally, I propose a new approach to properly scale the FK formulation when it is used as an approximation to the Arrhenius law.

### 4.1 Effect of a viscosity cut-off

One approach to avoid too high viscosity contrasts is to impose a maximum value of the viscosity. The effect of this cut-off has been studied by King (2009) and Stein & Hansen (2013), both finding that results are not affected as long as the cut-off viscosity is at least $10^4$ larger than the bottom viscosity. I did not extensively study the effect of a cut-off here, but still made a few tests for three cases with the real Arrhenius law: (1) $E_a = 300 \text{ kJ mol}^{-1}$ ($\gamma = 102$) with a bottom viscosity $\eta_b = 3.16 \times 10^{20} \text{ Pa.s}$ or (2) with a bottom viscosity of $\eta_b = 10^{21} \text{ Pa.s}$ and (3) $E_a = 100 \text{ kJ mol}^{-1}$ ($\gamma = 34$) with a bottom viscosity $\eta_b = 10^{21} \text{ Pa.s}$. Other parameters are given in Table 1. I found that the Nusselt number and interior temperature were deviating significantly from results obtained with the real Arrhenius law when the cut-off was equal or less than $7 \times 10^4$ times the bottom viscosity for case 1, $2 \times 10^5$ for case 2 and $10^4$ for case 3. The thermal profiles were unchanged in the three cases for a cut-off of $10^6$. This needs to be further explored, but the value $10^4$ for the cut-off seems a bit low for high viscosity contrasts (here for $E_a = 300 \text{ kJ mol}^{-1}$, i.e. $\gamma = 102$ for the Earth’s mantle).

<figure><img src="image_12.png" /><figcaption>Figure 10. Effect of $T_{\text{off}}^*$ for the Arrhenius law (eq. 12), for $E_a = 300 \text{ kJ mol}^{-1}$ and $\eta_b = 3.162 \times 10^{20} \text{ Pa.s}$, with $T_s = 300 \text{ K}$ and $T_b = 2000 \text{ K}$. (a) Temperature and Nusselt number, predicted (lines) and observed (symbols). (b) Total viscosity contrast $\gamma$ as a function of $T_{\text{off}}^*$ and value of $E_a$ that would give the same $\gamma$ for $T_{\text{off}}^* = T_s/\Delta T = 0.1765$ (‘Equivalent $E_a$’).</figcaption></figure>

### 4.2 Effect of the offset temperature $T_{\text{off}}^*$ in the Arrhenius law

To check that the scaling approach proposed here was valid for all formulations of the temperature-dependence of the viscosity, two cases for the Arrhenius law were studied: one with the normalized temperature surface $T_{\text{off}}^* = T_s/\Delta T$ corresponding to the actual choice made here for $T_s$ and $T_b$ (see Table 1), denoted as ‘real Arrhenius’, and one where $T_{\text{off}}^* = 1$, as is sometimes chosen for simulations of mantle convection with plate tectonics (e.g. Tackley 1996; Foley & Becker 2009; Rolf et al. 2017), denoted as ‘Arrhenius with offset’. The scaling for temperature and heat flux proposed in this study is valid for both cases, and to further check this, the value of $T_{\text{off}}^*$ was systematically varied. Results of simulations were compared to the predictions of the scaling proposed in Section 3.4 (see Fig. 10).

The obtained interior temperature $T_i^*$ and Nusselt number Nu are in very good agreement with the scaling laws. Results with ‘Arrhenius with offset’ in Section 3 (e.g. Fig. 9) already pointed to the important effect of $T_{\text{off}}^*$, but they were compared to ‘real Arrhenius’ for the same total viscosity contrast $\gamma$, which implied that very high values of $E_a$ had to be used when $T_{\text{off}}^* = 1$. Here, the activation energy $E_a$ is fixed to $300 \text{ kJ mol}^{-1}$, and the viscosity contrast $\gamma$ decreases with increasing $T_{\text{off}}^*$. This results in a thinner stagnant lid, lower interior temperature $T_i^*$ and larger heat flux Nu (Fig. 10b). With $T_{\text{off}}^* = T_s/\Delta T = 300/1700 = 0.1765$, $E_a = 300 \text{ kJ mol}^{-1}$ corresponds to a total viscosity contrast $\gamma = 102.24$ (upper horizontal dotted line in Fig. 10a). With $T_{\text{off}}^* = 1$, the viscosity contrast is only $\gamma = 10.61$, which would be obtained for $E_a = 31.14 \text{ kJ mol}^{-1}$ if $T_{\text{off}}^* = 0.1765$ were used (lower horizontal dotted line in Fig. 10a).
It is therefore essential to keep in mind that the use of a temperature $T_{off}^*$ larger than $T_s/\Delta T$ should only be used with the understanding that it is an important modification of the viscosity law and not merely a way to avoid very high viscosity values at the surface. The actual quantitative effect in the stagnant-lid regime can be predicted with the approach proposed here.

### 4.3 Further verification of the scaling approach: icy moons

So far in this study, the choices for the activation energy $E_a$ were made to study Earth's mantle, with parameters given in Table 1. Direct applications to Earth's dynamics, for instance to refine models of Earth's thermal history, are not yet possible with the scaling approach proposed here, as the effect of radioactive internal heating, curvature and compressibility are not considered. More importantly, Earth's mantle is not in a stagnant-lid regime.

A more straightforward possible domain of application concerns stagnant-lid convection for the ice I outer layer of icy bodies, when this layer is thin compared to the size of the whole body, so that curvature can be neglected. To check that the scaling approach proposed here can be valid in this context, I carried out a few experiments with parameters for a generic icy satellite of the size of Europa and with a similar surface temperature. Parameters are given in Table 2 and are similar to the ones chosen by Deschamps & Sotin (2001). The reference viscosity for ice I is at least 7 orders of magnitude smaller than the one of a silicate mantle, but with much colder temperatures, a smaller temperature jump $\Delta T$, and considering a 20-fold thinner layer, the Rayleigh number and the dimensionless activation energy are similar. But the smaller temperature jump $\Delta T$ can yield a higher dimensionless temperature offset $T_{off}^* = T_s/\Delta T$. The value tested here is $T_{off}^* = 100/150 = 0.667$, while it was $T_{off}^* = 0.176$ so far in this study considering Earth's mantle. In the previous Section 4.2, the scaling was shown to correctly predict the effect of varying $T_{off}^*$ for one value of $E_a^*$ (300 kJ mol$^{-1}$, corresponding to $E_a^* = 21.2$ when using $\Delta T = 1700$K for Earth's mantle). The tests done here can be seen as a further verification, with varying $E_a$ for a value of $T_{off}^*$ different from what has been tested so far. $E_a$ will be varied from 30 to 70 kJ mol$^{-1}$ ($E_a^* = 24$ to 56 for $\Delta T = 150$K for a generic icy moon here).

In order to obtain a scaling approach valid over a wide range of parameters, without any interior temperature predictions derived from previous studies, I chose in Section 3 to explore different viscosity laws with parameters that give the same total viscosity contrast $\gamma = \ln(\eta(T_s)/\eta(T_b))$. But when the FK formulation is used as an approximation of the Arrhenius law, it is more interesting to set $\theta^*$ in the FK approximation so that it gives an interior temperature similar to one obtained with the Arrhenius law. Previous authors (e.g. Solomatov & Moresi 1996; Reese et al. 1999; Deschamps & Sotin 2001; King 2009; Deschamps & Lin 2014; Harel et al. 2020) write $\theta^*$ as

$$ \theta_{\text{VTS}}^* = \frac{E_a \Delta T}{R T_i^2} = \frac{E_a^*}{(T_i^* + T_{\text{off}}^*)^2} \quad (46) $$

This corresponds to the 'viscous temperature scale' approach (eq. 5), denoted by the subscript 'VTS' here. Eq. (46) can be combined with $1 - T_i^* = a/\theta_{\text{VTS}}^*$ (eq. 29) to predict $T_i^*$, and it can then be shown that

$$ \theta_{\text{VTS}}^* = a \left[ 1 + T_{\text{off}}^* + \frac{E_a^*}{a} - \left( \frac{E_a^*}{a} \left( 1 + T_{\text{off}}^* + \frac{E_a^*}{4a} \right) \right)^{1/2} \right]^{-1} \quad (47) $$

where $a = 1.051$ was determined numerically from experiments in Section 3.2. I carried out experiments with the real Arrhenius law and with the FK approximation with this corresponding value of $\theta_{\text{VTS}}^*$, with two purposes: (i) check that the scaling approach proposed in Section 3.4 is valid in this context of icy moons and (ii) quantify the differences in terms of interior temperature $T_i^*$ and heat flux Nu between the 'viscous temperature scale' FK approximation and the real Arrhenius law. Results are shown in Fig. 11.

The FK approximation using $\theta_{\text{VTS}}^*$ systematically underpredicts the interior temperature $T_i$ and overpredicts the heat flux compared to the Arrhenius law. This is in agreement with Reese et al. (1999) and Harel et al. (2020), who also found this overprediction of the heat flux with the FK approximation. They indicated a heat flux ~20 per cent higher with the FK approximation compared to the Arrhenius law. The actual excess depends on the activation energy. Here, I find that the heat flux overprediction decreases from 20 to 11 per cent for $E_a$ increasing from 40 to 70 kJ mol$^{-1}$ (lower values of $E_a$ result in a non-stagnant regime with the FK approximation, which is not covered by the present scaling). For any set of parameters, as long as a stagnant-lid regime is observed, the predictive approach proposed in Section 3.4 can be used to estimate the expected difference between the FK approximation and the real Arrhenius law.

### 4.4 New scaling of the FK approximation

The predictive approach of Section 3.4 helps evaluate how the so-called FK approximation differs from the real Arrhenius law, as shown in the previous section. When the FK formulation is used as an approximation of the Arrhenius law, it has been proposed that $\theta$ should be written as $\theta = E_a/(RT_i^2)$. With $T_i^*$ predicted by eq. (29), the corresponding value of $\theta^*$ is then given by eq. (47). The interior temperature $T_i$ then obtained with the FK approximation significantly differs from the one obtained with the real Arrhenius law. A lower temperature and a higher heat flux is always obtained with the FK approximation (see Fig. 11). King (2009) showed that it is possible to get the same Nusselt number and RMS of the velocity within the model when one properly scales the parameters for the FK formulation compared to the Arrhenius law. He chose the parameters for the FK form via trial and error. I hereunder propose an approach relying on the scalings of Section 3.4 to choose those parameters.

As $T_i^*$ can be predicted for both the real Arrhenius law and for the FK approximation with the method given in Section 3.2, the parameter $\theta^*$ for the FK formulation can be chosen directly to give the same temperature $T_i^*$ as the Arrhenius law that it should approximate. It was shown that for the FK approximation, the predicted interior temperature is such that $1 - T_i^* = 1.051/\theta^*$, so that

$$ \theta_{T_i}^* = \frac{1.051}{1 - T_{i\text{Arr}}^*} \quad (48) $$

is the parameter that should be used to obtain the temperature $T_{i\text{Arr}}^*$. But the interior viscosity $\eta_{i\text{FK}}^*(T_{i\text{Arr}}^*)$ (eq. 13) with this approximation is then different from the one for the Arrhenius law $\eta_{i\text{Arr}}^*(T_{i\text{Arr}}^*)$ (eq. 12). With a different viscosity, the velocities within the convective part of the model, and therefore the heat transfer efficiency, are different. One should thus also scale the viscosity reference so that the same interior viscosity is obtained for both formulations. Since the FK approximation implies an interior temperature $T_i^*$ which follows eq. (48), and given the FK formulation (eq. 13), the dimensionless viscosity $\eta_{i\text{FK}}^*(T_i^*)$ in the well-mixed convective part of the model is always $\eta_{i\text{FK}}^*(T_i^*) = \exp(1.051) \approx 2.86$, regardless of the choice made for $\theta^*$. The interior viscosity for the real Arrhenius
<table><caption>**Table 2.** Parameters used for the ice I outer layer of an icy satellite. Data are similar to that used by Deschamps & Lin (2014) and references therein.</caption><thead><tr><th>Parameter</th><th>Symbol</th><th>Value</th></tr></thead><tbody><tr><td>Layer thickness</td><td><i>D</i></td><td>150 km</td></tr><tr><td>Gravity acceleration</td><td><i>g</i></td><td>1.3 m s<sup>-2</sup></td></tr><tr><td>Thermal expansion</td><td><i>α</i></td><td>1.6 × 10<sup>-4</sup> K<sup>-1</sup></td></tr><tr><td>Thermal diffusivity</td><td><i>κ</i></td><td>1.5 × 10<sup>-4</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Reference density</td><td><i>ρ</i><sub>b</sub></td><td>917 kg m<sup>-3</sup></td></tr><tr><td>Reference (bottom) temperature</td><td><i>T</i><sub>ref</sub> = <i>T</i><sub>b</sub></td><td>250 K</td></tr><tr><td>Surface temperature</td><td><i>T</i><sub>s</sub></td><td>100 K</td></tr><tr><td>Reference (bottom) viscosity</td><td><i>η</i><sub>ref</sub> = <i>η</i><sub>b</sub></td><td>5 × 10<sup>13</sup> Pa.s</td></tr><tr><td>Reference (bottom) Rayleigh number</td><td>Ra<sub>ref</sub> = Ra<sub>b</sub></td><td>1.32 × 10<sup>7</sup></td></tr><tr><td>Activation energy</td><td><i>E</i><sub>a</sub></td><td>30–70 kJ mol<sup>-1</sup></td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Observed (symbols) and predicted (lines) values (a) of the interior temperature $T_i$ and (b) of the heat flux. Physical parameters, for a generic icy moon, are given in Table 2. $\theta_{VTS}^*$ is computed as a function of $E_a^*$ with eq. (47).</figcaption></figure>

law $\eta_{Arr}^*(T_{iArr}^*)$ can be obtained from eq. (12), and it can be shown that for all cases tested here, it is lower than $\eta_{FK}^*(T_i^*) = \exp(1.051)$ and it increases with $E_a^*$. The coefficient $C$ by which the reference viscosity $\eta_b$ should be divided so that $\eta_{FK}^*(T_{iArr}^*) = \eta_{Arr}^*(T_{iArr}^*)$, or equivalently by which the Rayleigh number should be multiplied, is then:

$$C = \frac{\eta_{bArr}}{\eta_{bFK}} = \frac{Ra_{bFK}}{Ra_{bArr}} = \exp\left(1.051 - E_a^* \left(\frac{1}{T_{iArr}^* + T_{off}^*} - \frac{1}{1 + T_{off}^*}\right)\right). \quad (49)$$

Additional simulations were done to check this, for $\eta_{bArr} = 10^{21}$ Pa.s for the Arrhenius law with $E_a$ varying from 50 to 300 kJ mol$^{-1}$. Results are shown in Fig. 12 and Table 3, illustrating the validity of the approach. With the same interior temperature and heat flux for the two viscosity laws, the vertical profiles of temperature (not shown here) are also equivalent.

The scaling laws proposed here are valid only for a true stagnant-lid regime. Cases where non negligible deformations occur at the surface of the model do not match the predicted values of either the interior temperature or the Nusselt number (light symbols without a black edge in Fig. 12). The criterion chosen here to exclude non-stagnant cases relies on the RMS of the surface strain rate, that has to be less than 5 per cent of the RMS of the velocity within the model for the case to be considered as truly stagnant. Through a rapid analysis, I find that for the real Arrhenius law, this criterion is met

when the total viscosity contrast $\gamma = \ln(\eta(T_s)/\eta(T_b))$ between the top and bottom of the layer is larger than around 16 [i.e. $\eta(T_s)/\eta(T_b) \gtrsim 9 \times 10^6$]. With the surface and bottom temperatures $T_s = 300$ K and $T_b = 2000$ K used here, this corresponds to values of the activation energy $E_a$ larger than around 32 kJ mol$^{-1}$ (all cases in red in Fig. 12). For the FK approximation, a stagnant-lid regime is obtained if $\gamma > 11$ [$\eta(T_s)/\eta(T_b) \gtrsim 6 \times 10^4$]. If $\theta_T^*$ is used (eq. 48), this corresponds to $E_a \gtrsim 120$ kJ mol$^{-1}$. This indicates that the method proposed here to approximate the Arrhenius law by the FK approximation is not valid for low values of $E_a$, even when they actually correspond to the stagnant-lid regime with the real Arrhenius law (here for instance for $E_a$ comprised between 32 and 120 kJ mol$^{-1}$).

King (2009) indicated that even when the parameters for the FK approximation were properly scaled to match the Arrhenius law in terms of the Nusselt number and RMS velocity, the topography and geoid at the surface of the model still differed significantly. This was not studied here, but since the total viscosity contrast $\gamma$, and therefore the surface viscosities, are very different (see column ‘Viscosity contrast’ in Table 3), I expect the same conclusions. The extent of the differences in geoid and topography as a function of $E_a$ remains to be studied for the Arrhenius law with values of $E_a$ that yield viscosity contrasts $\gamma$ larger than 20 ($E_a > 60$ kJ mol$^{-1}$ for Earth’s mantle parameters in Table 1).
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Observed (symbols) and predicted (lines) values of the interior temperature $T_i^*$ (a) and dimensionless heat flux (Nusselt number Nu) (b) as a function of the activation energy $E_a$, for $\eta_b = 10^{21}$ Pa.s for cases with the Arrhenius law, and with $\eta_b$ given in Table 3 for the FK approximation. Other parameters are given in Table 1.</figcaption></figure>

<table><caption>**Table 3.** Scaled parameters $\theta_{T_i}^*$ (eq. 48) and $C = \eta_{bArr}/\eta_{bFK} = Ra_{bFK}/Ra_{bArr}$ (eq. 49), so that the FK approximation matches the interior temperature and heat flux of the corresponding Arrhenius law for varying $E_a$. Arrhenius law is for $\eta_b = 10^{21}$ Pa.s ($Ra_b = 3.99 \times 10^7$). The total viscosity contrast $\gamma = \ln[\eta(T_s)/\eta(T_b)]$ is given for the real Arrhenius law (Arr subscript) and for the equivalent FK approximation. The predicted interior temperature and Nusselt number (Section 3.4) are given in the columns ‘Scaling predictions’. Results of simulations are given in the four right-hand columns.</caption><thead><tr><th rowspan="2">E<sub>a</sub> (kJ mol<sup>-1</sup>)</th><th colspan="2">Parameters</th><th colspan="2">Viscosity contrast</th><th colspan="3">Scaling predictions</th><th colspan="3">Simulations results</th></tr><tr><th>θ<sup>*</sup><sub>T<sub>i</sub></sub></th><th>C</th><th>γ<sub>Arr</sub></th><th>γ<sub>FK</sub></th><th>T<sup>*</sup><sub>i</sub></th><th>Nu</th><th>T<sup>*</sup><sub>IArr</sub></th><th>T<sup>*</sup><sub>IFK</sub></th><th>Nu<sub>Arr</sub></th><th>Nu<sub>FK</sub></th></tr></thead><tbody><tr><td>50</td><td>7.0082</td><td>1.844</td><td>17.04</td><td>7.01</td><td>0.850</td><td>9.884</td><td>0.847</td><td>0.760<sup>a</sup></td><td>9.676</td><td>18.742<sup>a</sup></td></tr><tr><td>75</td><td>8.5197</td><td>1.687</td><td>25.56</td><td>8.52</td><td>0.877</td><td>7.557</td><td>0.876</td><td>0.840<sup>a</sup></td><td>7.247</td><td>12.059<sup>a</sup></td></tr><tr><td>100</td><td>9.9595</td><td>1.582</td><td>34.08</td><td>9.96</td><td>0.895</td><td>6.132</td><td>0.894</td><td>0.886<sup>a</sup></td><td>6.114</td><td>7.830<sup>a</sup></td></tr><tr><td>125</td><td>11.3582</td><td>1.506</td><td>42.60</td><td>11.36</td><td>0.908</td><td>5.164</td><td>0.906</td><td>0.911</td><td>5.244</td><td>5.624</td></tr><tr><td>150</td><td>12.7306</td><td>1.448</td><td>51.12</td><td>12.73</td><td>0.917</td><td>4.461</td><td>0.920</td><td>0.921</td><td>4.275</td><td>4.496</td></tr><tr><td>200</td><td>15.4260</td><td>1.366</td><td>68.16</td><td>15.43</td><td>0.932</td><td>3.509</td><td>0.936</td><td>0.934</td><td>3.357</td><td>3.512</td></tr><tr><td>250</td><td>18.0810</td><td>1.309</td><td>85.20</td><td>18.08</td><td>0.942</td><td>2.891</td><td>0.945</td><td>0.946</td><td>2.695</td><td>2.747</td></tr><tr><td>300</td><td>20.7113</td><td>1.268</td><td>102.24</td><td>20.71</td><td>0.949</td><td>2.458</td><td>0.952</td><td>0.953</td><td>2.286</td><td>2.333</td></tr></tbody></table>

Other parameters are given in Table 1.

## 5 CONCLUSIONS

A new scaling for the interior temperature in the well-mixed part of a fluid undergoing stagnant-lid convection is proposed. This scaling is equivalent to the one relying on a viscous temperature scale $\Delta T_v$ (Davaille & Jaupart 1993) in the case of the FK approximation, but I show that this latter approach is not valid in the case of an Arrhenius viscosity law. The scalings proposed here for the interior temperature and heat flux, even if not expressed as direct equations, offer a rapid method to check by how much viscosity law approximations will differ from the real Arrhenius law, as shown in Section 4. A Python script is provided in the Supporting Information to compute the predicted interior temperature $T_i^*$ and heat flux Nu, from which an approximated vertical profile of temperature can be deduced.

In studies of the thermal history of planets, other rocky bodies and icy satellites, it is important to keep in mind that a material defined by a constant activation energy $E_a$ undergoes variable viscosity contrasts when the thermal state of the system evolves. This cannot

be simply approached by an FK approximation with a constant viscosity contrast $\gamma = \theta^*$. Indeed, when put in a dimensionless form, the Arrhenius law depends on a dimensionless temperature offset $T_{off}^* = T_s/\Delta T$, with $\Delta T$ evolving through time. If a FK approximation is used as a simplification of the Arrhenius law, the parameter $\theta^*$ should be scaled as a function of $T_i^*$ (see Section 4.4), which, as $\Delta T$, changes through time. To get the same heat loss as with the approximated Arrhenius law, the reference viscosity should also be modified, effectively rendering the scaling of the FK approximation not straightforward.

For further applications to the thermal history of rocky bodies in the stagnant-lid regime, the role of internal heating will be studied next with a real Arrhenius law. The scaling would also need to be extended to 3-D spherical convection. However, the SFS2D convection code used here relies on directly solving for the stream function in eq. (11), which I infer is what allows a fast numerical convergence in models with very high viscosity contrasts. This formulation can be extended to a 2-D cylindrical or 2-D spherical annulus geometry, but not easily to a 3-D Cartesian or spherical geometry. As a first

<footnote><sup>a</sup>Indicates cases that are not in the stagnant-lid regime for the FK approximation.</footnote>
step, a 2-D spherical annulus geometry (Hernlund & Tackley 2008) will be implemented.

In this study, no attempt was made to constrain the range of parameters that generate stagnant-lid convection. The transitional regime between a fully mobile top surface and the stagnant-lid regime was not considered. I only checked that the surface was not deformed in order to exclude results that were not in a fully stagnant-lid state. It will be important in a future study to check how this transition occurs when considering the Arrhenius law, and not the FK approximation. A few preliminary results here point to the fact that stagnant-lid convection is achieved with a lesser value of the viscosity contrast with the real Arrhenius law compared to the FK approximation (see for instance Fig. 12, where some cases with low values of $E_a$ give a stagnant-lid regime with the real Arrhenius law and are non-stagnant with the corresponding scaled FK approximation).

With either a 2-D spherical annulus or 3-D spherical geometry, Yao et al. (2014) and Guerrero et al. (2018) showed that the development of the stagnant-lid regime depended on the curvature of the shell. For a fixed size of a planet or moon, stagnant-lid convection requires a larger viscosity contrast when the core is smaller (i.e. for a more curved layer). The curvature of the layer was also shown to strongly impact the interior temperature, with lower interior temperatures obtained for smaller cores. It will be important to study how these results are affected by the use of the Arrhenius law instead of the FK approximation for the viscosity temperature dependence.

## SUPPORTING INFORMATION

Supplementary data are available at GJI online.

suppl_data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## ACKNOWLEDGMENTS

The author thanks Julian Lowman, an anonymous reviewer and the editor Gaël Choblet for their very useful comments and the time taken to help improve the presentation of this study.

## DATA AVAILABILITY

All data generated for this article was obtained with the convection code SFS2D, which was written by the author. It is available on request to the author (cecile.grigne@univ-brest.fr). A Python script to solve the scaling laws of Section 3.4 is given as Supporting Information.

## REFERENCES

* Christensen, U., 1984. Heat transport by variable viscosity convection and implications for the Earth's thermal evolution, *Phys. Earth planet. Inter.*, **35**, 264–282.

* Davaillé, A. & Jaupart, C., 1993. Transient high-Rayleigh-number thermal convection with large viscosity variations, *J. Fluid. Mech.*, **253**, 141–166.

* Deschamps, F., 2021. Stagnant lid convection with temperature-dependent thermal conductivity and the thermal evolution of icy worlds, *Geophys. J. Int.*, **224**, 1870–1889.

* Deschamps, F. & Lin, J.-R., 2014. Stagnant lid convection in 3D-Cartesian geometry: scaling laws and applications to icy moons and dwarf planets, *Phys. Earth planet. Inter.*, **229**, 40–54.

* Deschamps, F. & Sotin, C., 2000. Inversion of two-dimensional numerical convection experiments for a fluid with a strongly temperature-dependent viscosity, *Geophys. J. Int.*, **143**(1), 204–218.

* Deschamps, F. & Sotin, C., 2001. Thermal convection in the outer shell of large icy satellites, *J. geophys. Res.*, **106**(E3), 5107–5121.

* Deschamps, F. & Vilella, K., 2021. Scaling laws for mixed-heated stagnant-lid convection and application to Europa, *J. geophys. Res.*, **126**(10), e2021JE006963, doi:10.1029/2021JE006963.

* Dumoulin, C., Doin, M.-P. & Fleitout, L., 1999. Heat transport in stagnant lid convection with temperature- and pressure-dependent Newtonian or non-Newtonian rheology, *J. geophys. Res.*, **104**, 12 759–12 778.

* Ferrick, A.L. & Korenaga, J., 2023. Generalizing scaling laws for mantle convection with mixed heating, *J. geophys. Res.*, **128**(5), e2023JB026398, doi:10.1029/2023JB026398.

* Foley, B. & Becker, T., 2009. Generation of plate-like behavior and mantle heterogeneity from a spherical, viscoplastic convection model, *Geochem. Geophys. Geosyst.*, **10**(8), doi:10.1029/2009GC002378.

* Fowler, A., 1985. A simple model of convection in the terrestrial planets, *Geophys. Astrophys. Fluid Dyn.*, **31**, 283–309.

* Gerya, T., 2019. *Introduction to Numerical Geodynamic Modelling*, 2nd edn, Cambridge Univ. Press.

* Grasset, O. & Parmentier, E., 1998. Thermal convection in a volumetrically heated, infinite Prandtl number fluid with strongly temperature-dependent viscosity: implications for planetary thermal evolution, *J. geophys. Res.*, **103**, 18 171–18 181.

* Grigné, C., Labrosse, S. & Tackley, P., 2005. Convective heat transfer as a function of wavelength. Implications for the cooling of the Earth, *J. geophys. Res.*, **110**(B3), doi:10.1029/2004JB003376.

* Guerrero, J.M., Lowman, J.P., Deschamps, F. & Tackley, P.J., 2018. The influence of curvature on convection in a temperature-dependent viscosity fluid: implications for the 2-D and 3-D modeling of moons, *J. geophys. Res.*, **123**(7), 1863–1880.

* Harel, L., Dumoulin, C., Choblet, G., Tobie, G. & Besserer, J., 2020. Scaling of heat transfer in stagnant lid convection for the outer shell of icy moons: influence of rheology, *Icarus*, **338**, doi:10.1016/j.icarus.2019.113448.

* Hernlund, J.W. & Tackley, P.J., 2008. Modeling mantle convection in the spherical annulus, *Phys. Earth planet. Inter.*, **171**(1), 48–54.

* Hewitt, J.M., McKenzie, D.P. & Weiss, N.O., 1980. Large aspect ratio cells in two-dimensional thermal convection, *Earth planet. Sci. Lett.*, **51**(2), 370–380.

* Hussmann, H., Sotin, C. & Lunine, J.I., 2015. 10.18 - Interiors and evolution of Icy Satellites, in *Treatise on Geophysics*, Vol. **10**, 2nd edn, pp. 605–635, ed. Schubert, G., Elsevier.

* Jarvis, G.T., 1984. Time-dependent convection in the Earth's mantle, *Phys. Earth planet. Inter.*, **36**(3), 305–327.

* Karato, S. & Wu, P., 1993. Rheology of the upper mantle: a synthesis, *Science*, **260**, 771–778.

* King, S., 2009. On topography and geoid from 2-D stagnant lid convection calculations, *Geochem. Geophys. Geosyst.*, **11**(3), doi:10.1029/2008GC002250.

* Kirby, S., 1983. Rheology of the lithosphere, *Rev. Geophys. Space Phys.*, **21**, 1458–1487.

* Korenaga, J., 2009. Scaling of stagnant-lid convection with Arrhenius rheology and the effects of mantle melting, *Geophys. J. Int.*, **179**, 154–170.

* Langemeyer, S.M., Lowman, J.P. & Tackley, P.J., 2021. Global mantle convection models produce transform offsets along divergent plate boundaries, *Commun. Earth Environ.*, **2**(1), 69, doi:10.1038/s43247-021-00139-1.

* McKenzie, D., Roberts, J. & Weiss, N., 1973. Numerical models of convection in the Earth's mantle, *Tectonophysics*, **19**(2), 89–103.

* Moresi, L.-N. & Solomatov, V., 1995. Numerical investigation of 2D convection with extremely large viscosity variations, *Phys. Fluids*, **7**(9), 2154–2162.
Morris, S. & Canright, D., 1984. A boundary-layer analysis of Benard convection in a fluid of strongly temperature-dependent viscosity, *Phys. Earth planet. Inter.*, **36**, 355-373.

Nakakuki, T., Hamada, C. & Tagawa, M., 2008. Generation and driving forces of plate-like motion and asymmetric subduction in dynamical models of an integrated mantle–lithosphere system, *Phys. Earth planet. Inter.*, **166**(3), 128–146.

Nataf, H. & Richter, F., 1982. Convection experiments in fluids with highly temperature-dependent viscosity and the thermal evolution of the planets, *Phys. Earth planet. Sci.*, **29**, 320–329.

Ogawa, M., Schubert, G. & Zebib, A., 1991. Numerical simulations of three-dimensional thermal convection in a fluid with strongly temperature-dependent viscosity, *J. Fluid Mech.*, **233**, 299–328.

Olson, P., 1987. A comparison of heat transfer laws for mantle convection at very high Rayleigh numbers, *Phys. Earth planet. Inter.*, **48**(1), 153–160.

Peaceman, D. & Rachford, H., Jr, 1955. The numerical solution of parabolic and elliptic differential equations, *J. Soc. Indust. Appl. Math.*, **3**(1), 28–41.

Reese, C., Solomatov, V. & Moresi, L.-N., 1999. Non-Newtonian stagnant lid convection and magmatic resurfacing on Venus, *Icarus*, **139**(1), 67–80.

Richter, F., Nataf, H. & Daly, S., 1983. Heat transfer and horizontally averaged temperature of convection with large viscosity variations, *J. Fluid Mech.*, **129**, 173–192.

Roe, P., 1986. Characteristic-based schemes for the Euler equations, *Ann. Rev. Fluid Mech.*, **18**, 337–365.

Rolf, T., Capitanio, F. & Tackley, P., 2017. Constraints on mantle viscosity structure from continental drift histories in spherical mantle convection models, *Tectonophysics*, **746**, 339–351.

Schubert, G., Turcotte, D.L. & Olson, P., 2001. Mantle Convection in the Earth and Planets, Cambridge Univ. Press.

Solomatov, V., 1995. Scaling of temperature- and stress-dependent viscosity convection, *Phys. Fluids*, **7**, 266–274.

Solomatov, V. & Moresi, L.-N., 1996. Stagnant lid convection on Venus, *J. geophys. Res.*, **101**(E2), 4737–4753.

Solomatov, V. & Moresi, L.-N., 2000. Scaling of time-dependent stagnant lid convection: application to small-scale convection on Earth and other terrestrial plantes, *J. Geophys. Res.*, **105**(B9), 21795–21817.

Stein, C. & Hansen, U., 2013. Arrhenius rheology versus Frank–Kamenetskii rheology – implications for mantle dynamics, *Geochem. Geophys. Geosyst.*, **14**(8), doi:10.1002/ggge.20158.

Tackley, P., 1996. Effects of strongly variable viscosity on three-dimensional compressible convection in planetary mantles, *J. geophys. Res.*, **101**, 3311–3332.

Tackley, P., 2008. Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid, *Phys. Earth planet. Inter.*, **171**(1), 7–18.

Thiriet, M., Breuer, D., Michaut, C. & Plesa, A.-C., 2019. Scaling laws of convection for cooling planets in a stagnant lid regime, *Phys. Earth planet. Sci.*, **286**, 138–153.

Tosi, N. et al., 2015. A community benchmark for viscoplastic thermal convection in a 2-D square box, *Geochem. Geophys. Geosyst.*, **16**(7), 2175–2196.

Turcotte, D. & Oxburgh, E., 1967. Finite amplitude convection cells and continental drift, *J. Fluid Mech.*, **28**, 29–42.

Turcotte, D. & Oxburgh, E., 1969. Convection in a mantle with variable physical properties, *J. geophys. Res.*, **74**(6), 1458–1474.

Turcotte, D. & Schubert, G., 2002. *Geodynamics*, Cambridge Univ. Press.

Weertman, J., 1983. Creep deformation of ice, *Ann. Rev. Earth planet. Sci.*, **11**(1), 215–240.

Yao, C., Deschamps, F., Lowman, J.P., Sanchez-Valle, C. & Tackley, P.J., 2014. Stagnant lid convection in bottom-heated thin 3-D spherical shells: influence of curvature and implications for dwarf planets and icy moons, *J. geophys. Res.*, **119**(8), 1895–1913.

## APPENDIX: BENCHMARK OF SFS2D

The benchmark by Tosi *et al.* (2015) consists of five cases in 2-D boxes, with a viscosity that depends on temperature (according to the Frank–Kamenetskii approximation) and pressure, for both Newtonian and non-linear rheologies. The full description of the cases can be found in Tosi *et al.* (2015). Parameters for the five cases are given in their Table 1 and the results obtained by 11 numerical codes are given in their Table 2.

For direct comparison, the results obtained with SFS2D are presented here in Table A1 in the same order as in Tosi *et al.* (2015), using the same symbols. For all cases, results of SFS2D are in good agreement with the median results for the 11 codes of the benchmark. A representation of some mean values for the temperature, heat flux and velocity obtained with SFS2D for Cases 1 to 4 of the benchmark, which are all steady-state, compared to the results of the 11 participating codes in Tosi *et al.* (2015) is also given in Fig. A1, showing the good agreement between SFS2D and the majority of the codes.

A bifurcation analysis (Case 5b) was also part of the benchmark, with a non-linear rheology defined by a yield stress whose value is systematically increased. A first transition between a mobile lid and a periodic lid regime, followed by a second transition between the periodic lid and the stagnant-lid regime, are observed. Results for SFS2D are shown in Fig. A2, using the same scales and transition zone (grey area) as fig. 7 of Tosi *et al.* (2015). This bifurcation analysis confirms that SFS2D agrees with the majority of the codes tested in the benchmark.
<table><caption>**Table A1.** Tosi *et al.* (2015) benchmark results, given as minimum, maximum and median values obtained among the 11 participating codes, and results for SFS2D, all for a resolution 128 × 128. Results are shown for three of the tested flux limiter schemes (minmod, koren and superbee) for the advection-diffusion of the temperature (eq. 8).</caption><thead><tr><th></th><th colspan="3">Tosi et al. benchmark</th><th colspan="3">SFS2D</th></tr><tr><th></th><th>min.</th><th>max.</th><th>median</th><th>minmod</th><th>koren</th><th>superbee</th></tr></thead><tbody><tr><td colspan="7" style="text-align: center;">Case 1</td></tr><tr><td>&lt;T&gt;</td><td>0.7758</td><td>0.779</td><td>0.7759</td><td>0.7764</td><td>0.7761</td><td>0.7758</td></tr><tr><td>Nu<sub>top</sub></td><td>3.3129</td><td>3.5889</td><td>3.4214</td><td>3.4223</td><td>3.4185</td><td>3.4152</td></tr><tr><td>Nu<sub>bot</sub></td><td>3.313</td><td>3.4259</td><td>3.4142</td><td>3.4223</td><td>3.4185</td><td>3.4152</td></tr><tr><td>u<sub>RMS</sub></td><td>248.9252</td><td>296.6156</td><td>249.573</td><td>252.1515</td><td>250.0018</td><td>248.0620</td></tr><tr><td>u<sub>RMS</sub><sup>surf</sup></td><td>1.1114</td><td>1.8999</td><td>1.8698</td><td>1.8786</td><td>1.8797</td><td>1.8805</td></tr><tr><td>u<sub>RMS</sub><sup>max</sup></td><td>1.5329</td><td>2.6477</td><td>2.6104</td><td>2.6190</td><td>2.6206</td><td>2.6218</td></tr><tr><td>&lt;W&gt;</td><td>2.369</td><td>2.5548</td><td>2.4246</td><td>2.4162</td><td>2.4149</td><td>2.4146</td></tr><tr><td>&lt;φ&gt;/Ra</td><td>2.31916</td><td>2.4333</td><td>2.4182</td><td>2.4163</td><td>2.4150</td><td>2.4146</td></tr><tr><td>δ</td><td>&lt;0.01 per cent</td><td>9.22 per cent</td><td>&lt;0.003 per cent</td><td>&lt;0.002 per cent</td><td>&lt;0.0007 per cent</td><td></td></tr><tr><td colspan="7" style="text-align: center;">Case 2</td></tr><tr><td>&lt;T&gt;</td><td>0.5951</td><td>0.6071</td><td>0.6034</td><td>0.6033</td><td>0.6030</td><td>0.6023</td></tr><tr><td>Nu<sub>top</sub></td><td>8.5115</td><td>9.2396</td><td>8.556</td><td>8.5529</td><td>8.5423</td><td>8.5314</td></tr><tr><td>Nu<sub>bot</sub></td><td>8.399</td><td>9.2409</td><td>8.5491</td><td>8.5529</td><td>0.5423</td><td>8.5314</td></tr><tr><td>u<sub>RMS</sub></td><td>140.1871</td><td>160.499</td><td>140.7087</td><td>140.5655</td><td>140.2271</td><td>139.9392</td></tr><tr><td>u<sub>RMS</sub><sup>surf</sup></td><td>102.85</td><td>123.8687</td><td>104.585</td><td>104.6307</td><td>104.5176</td><td>104.4226</td></tr><tr><td>u<sub>RMS</sub><sup>max</sup></td><td>119.412</td><td>162.423</td><td>121.766</td><td>121.7703</td><td>121.6821</td><td>121.5005</td></tr><tr><td>η<sub>min</sub></td><td>1.9607 × 10<sup>-5</sup></td><td>2.08 × 10<sup>-5</sup></td><td>1.9854 × 10<sup>-5</sup></td><td>1.9854 × 10<sup>-5</sup></td><td>1.9854 × 10<sup>-5</sup></td><td>1.9854 × 10<sup>-5</sup></td></tr><tr><td>η<sub>max</sub></td><td>0.00477</td><td>1.9421</td><td>1.7887</td><td>1.9898</td><td>1.9721</td><td>1.9503</td></tr><tr><td>&lt;W&gt;</td><td>7.5199</td><td>8.7223</td><td>7.5521</td><td>7.5416</td><td>7.5398</td><td>7.5288</td></tr><tr><td>&lt;φ&gt;/Ra</td><td>7.469</td><td>10.026</td><td>7.55153</td><td>7.5339</td><td>7.5271</td><td>7.5211</td></tr><tr><td>δ</td><td>&lt;0.01 per cent</td><td>13 per cent</td><td>0.104 per cent</td><td>0.169 per cent</td><td>0.102 per cent</td><td></td></tr><tr><td colspan="7" style="text-align: center;">Case 3</td></tr><tr><td>&lt;T&gt;</td><td>0.7241</td><td>0.7305</td><td>0.7275</td><td>0.7276</td><td>0.7276</td><td>0.7274</td></tr><tr><td>Nu<sub>top</sub></td><td>2.9311</td><td>3.0918</td><td>3.0324</td><td>3.0326</td><td>3.0311</td><td>3.0300</td></tr><tr><td>Nu<sub>bot</sub></td><td>2.9311</td><td>3.041</td><td>3.03025</td><td>3.0326</td><td>3.0311</td><td>3.0300</td></tr><tr><td>u<sub>RMS</sub></td><td>99.8701</td><td>111.6121</td><td>100.1208</td><td>100.2847</td><td>100.0313</td><td>99.8945</td></tr><tr><td>u<sub>RMS</sub><sup>surf</sup></td><td>1.356</td><td>2.0916</td><td>2.07299</td><td>2.0799</td><td>2.0812</td><td>2.0810</td></tr><tr><td>u<sub>RMS</sub><sup>max</sup></td><td>1.8806</td><td>2.9207</td><td>2.89501</td><td>2.9044</td><td>2.9049</td><td>2.9061</td></tr><tr><td>η<sub>min</sub></td><td>4.7907 × 10<sup>-5</sup></td><td>0.0001</td><td>4.8047 × 10<sup>-5</sup></td><td>4.8212 × 10<sup>-5</sup></td><td>4.8019 × 10<sup>-5</sup></td><td>4.8002 × 10<sup>-5</sup></td></tr><tr><td>η<sub>max</sub></td><td>0.901</td><td>1.0</td><td>0.9988</td><td>1.0</td><td>1.0</td><td>1.0</td></tr><tr><td>&lt;W&gt;</td><td>2.0028</td><td>2.056</td><td>2.034</td><td>2.0303</td><td>2.0307</td><td>2.0309</td></tr><tr><td>&lt;φ&gt;/Ra</td><td>1.865</td><td>2.037</td><td>2.0302</td><td>2.0324</td><td>2.0325</td><td>2.0330</td></tr><tr><td>δ</td><td>&lt;0.01 per cent</td><td>13.29 per cent</td><td>0.103 per cent</td><td>0.089 per cent</td><td>0.105 per cent</td><td></td></tr><tr><td colspan="7" style="text-align: center;">Case 4</td></tr><tr><td>&lt;T&gt;</td><td>0.527304</td><td>0.5364</td><td>0.5277</td><td>0.5277</td><td>0.5275</td><td>0.5274</td></tr><tr><td>Nu<sub>top</sub></td><td>6.5572</td><td>7.4376</td><td>6.6156</td><td>6.6163</td><td>6.6119</td><td>6.6077</td></tr><tr><td>Nu<sub>bot</sub></td><td>6.5243</td><td>7.4376</td><td>6.6158</td><td>6.6163</td><td>6.6119</td><td>6.6077</td></tr><tr><td>u<sub>RMS</sub></td><td>78.6652</td><td>98.8912</td><td>79.1105</td><td>79.0665</td><td>79.0245</td><td>78.9918</td></tr><tr><td>u<sub>RMS</sub><sup>surf</sup></td><td>74.1719</td><td>93.5057</td><td>75.0827</td><td>75.0843</td><td>75.0598</td><td>75.0387</td></tr><tr><td>u<sub>RMS</sub><sup>max</sup></td><td>87.6118</td><td>123.046</td><td>88.9146</td><td>88.8714</td><td>88.8390</td><td>88.8113</td></tr><tr><td>η<sub>min</sub></td><td>1.9167 × 10<sup>-4</sup></td><td>1.9900 × 10<sup>-4</sup></td><td>1.9220 × 10<sup>-4</sup></td><td>1.9183 × 10<sup>-4</sup></td><td>1.9184 × 10<sup>-4</sup></td><td>1.9184 × 10<sup>-4</sup></td></tr><tr><td>η<sub>max</sub></td><td>0.0057</td><td>1.9834</td><td>1.6773</td><td>1.9250</td><td>1.9161</td><td>1.9218</td></tr><tr><td>&lt;W&gt;</td><td>5.5903</td><td>6.7510</td><td>5.6216</td><td>5.6104</td><td>5.6085</td><td>5.6074</td></tr><tr><td>&lt;φ&gt;/Ra</td><td>5.5434</td><td>7.7860</td><td>5.6174</td><td>5.6098</td><td>5.6079</td><td>5.6068</td></tr><tr><td>δ</td><td>&lt;0.01 per cent</td><td>13.29 per cent</td><td>0.011 per cent</td><td>0.011 per cent</td><td>0.011 per cent</td><td>0.011 per cent</td></tr><tr><td colspan="7" style="text-align: center;">Case 5a</td></tr><tr><td>P</td><td>0.0760</td><td>0.0792</td><td>0.0769</td><td>0.0770</td><td>0.0770</td><td>0.0770</td></tr><tr><td>&lt;T&gt;<sub>min</sub></td><td>0.6491</td><td>0.7119</td><td>0.6510</td><td>0.6513</td><td>0.6511</td><td>0.6509</td></tr><tr><td>&lt;T&gt;<sub>max</sub></td><td>0.6671</td><td>0.7119</td><td>0.6693</td><td>0.6694</td><td>0.6692</td><td>0.6690</td></tr><tr><td>Nu<sub>top</sub><sub>min</sub></td><td>2.6787</td><td>2.9597</td><td>2.6902</td><td>2.6808</td><td>2.6804</td><td>2.6799</td></tr><tr><td>Nu<sub>top</sub><sub>max</sub></td><td>2.9597</td><td>7.4314</td><td>7.4019</td><td>7.4324</td><td>7.4243</td><td>7.4165</td></tr><tr><td>u<sub>RMS</sub><sub>min</sub></td><td>41.2488</td><td>76.5018</td><td>41.4920</td><td>41.2626</td><td>41.2191</td><td>41.1832</td></tr><tr><td>u<sub>RMS</sub><sub>max</sub></td><td>76.5018</td><td>100.5380</td><td>98.8230</td><td>99.5580</td><td>99.9434</td><td>99.6186</td></tr><tr><td>&lt;φ&gt;<sub>min</sub></td><td>1.3307</td><td>1.8644</td><td>1.3447</td><td>1.3310</td><td>1.3324</td><td>1.3340</td></tr><tr><td>&lt;φ&gt;<sub>max</sub></td><td>1.8644</td><td>9.5400</td><td>9.3495</td><td>9.5166</td><td>9.5003</td><td>9.4868</td></tr></tbody></table>
<figure><img src="image_15.png" /><figcaption>**Figure A1.** Results for cases 1 to 4 obtained with SFS2D (red star) compared to the 11 codes of the benchmark (blue ticks), for the mean temperature <T>, the top heat flux Nu<sub>top</sub>, the RMS of the velocity u<sub>RMS</sub> and the maximum of the surface velocity u<sub>max</sub><sup>surf</sup>. SFS2D results are for a resolution of 128 × 128 and the minmod flux limiter.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure A2.** Bifurcation analysis (Case 5b) for SFS2D, as carried out in Tosi *et al.* (2015, see their fig. 7).</figcaption></figure>