1. Introduction
Some rocky planets and large satellites, such as the Moon (Laneuville et al. Reference Laneuville, Wieczorek, Breuer and Tosi2013) and Mars (Reese & Solomatov Reference Reese and Solomatov2006), maintain a significant hemispherical dichotomy in the thickness of their crust. The cause of this phenomenon remains poorly understood, with proposed causes of the lunar crustal dichotomy, including the effect of radiative heating from Earth (Roy, Wright & Sigurđsson Reference Roy, Wright and Sigurđsson2014), asymmetric radiogenic heating (Laneuville et al. Reference Laneuville, Wieczorek, Breuer and Tosi2013), tidal heating (Garrick-Bethell et al. Reference Garrick-Bethell, Perera, Nimmo and Zuber2014) and degree-1 modes of mantle convection (Morison et al. Reference Morison, Labrosse, Deguen and Alboussière2019). However, each of these approaches raises significant further questions. The time scale of millennia for crustal formation required by Roy et al. (Reference Roy, Wright and Sigurđsson2014) is at odds with the more orthodox scale of tens to hundreds of millions of years for crustal formation by flotation (Elkins-Tanton, Burgess & Yin Reference Elkins-Tanton, Burgess and Yin2011). Asymmetric distribution of radiogenic elements requires a chemical dichotomy in the initial conditions. Tidal heating is, to leading order, an effect at spherical harmonic degree 2. It is difficult to see how this mechanism would induce a major difference between the two hemispheres. The degree-1 mode of mantle convection suggested by Morison et al. (Reference Morison, Labrosse, Deguen and Alboussière2019) does not provide a physical mechanism, but rather posits that effects deep in the mantle can perhaps influence the crust. Additionally, large impacts have been suggested as the cause of Mars's dichotomy (Andrews-Hanna, Zuber & Banerdt Reference Andrews-Hanna, Zuber and Banerdt2008; Marinova, Aharonson & Asphaug Reference Marinova, Aharonson and Asphaug2008; Golabek et al. Reference Golabek, Keller, Gerya, Zhu, Tackley and Connolly2011).
For both Mars (Elkins-Tanton Reference Elkins-Tanton2012) and the Moon (Wieczorek et al. Reference Wieczorek2006), a significant fraction of the crust (if not all) is thought to have originated in the time period when a large proportion of the planet was molten, in what is known as a magma ocean. While ultimately driven by radiative heat loss at the surface, the manner of cooling of both magma oceans (Elkins-Tanton Reference Elkins-Tanton2012; Michaut & Neufeld Reference Michaut and Neufeld2022) and the solid mantle (Labrosse & Jaupart Reference Labrosse and Jaupart2007) is primarily convective. In these settings the temperature dependence of the viscosity of the convecting material is of major importance, given the range in viscosity between a pure melt of 
${\sim }0.1$ Pa s and a pure solid at 
${\sim }10^{21}$ Pa s (Reese & Solomatov Reference Reese and Solomatov2006). Due to the inherent difficulty in driving fluid motions in the cold, very viscous stagnant lid, a thick surface boundary layer commonly occurs in convecting fluids with a temperature-dependent viscosity ratio of several orders of magnitude (Davaille & Jaupart Reference Davaille and Jaupart1994; Solomatov Reference Solomatov1995). Similar stagnant lids have also been observed and studied in laboratory experiments (Richter, Nataf & Daly Reference Richter, Nataf and Daly1983; Davaille & Jaupart Reference Davaille and Jaupart1993) and in numerical models (Stengel, Oliver & Booker Reference Stengel, Oliver and Booker1982; Christensen Reference Christensen1984; Morris & Canright Reference Morris and Canright1984; Ogawa, Schubert & Zebib Reference Ogawa, Schubert and Zebib1991; Moresi & Solomatov Reference Moresi and Solomatov1995), both of which have provided useful insight into the behaviour of convection with temperature-dependent rheologies.
A stagnant lid occurs when the cold boundary layer is too viscous to support convective motions. The interior of the fluid remains vigorously stirred by convection and is therefore well characterised by a near-uniform temperature (Nataf & Richter Reference Nataf and Richter1982; Richter et al. Reference Richter, Nataf and Daly1983), while the stagnant lid transmits heat by conduction. A thinner unstable boundary layer therefore connects the stagnant lid to the convecting interior, and acts in a similar manner to boundary layers found in isoviscous convection, and is defined by having the viscosity across it vary by approximately only one order of magnitude (Morris & Canright Reference Morris and Canright1984; Smith Reference Smith1988; Davaille & Jaupart Reference Davaille and Jaupart1993). The temperature contrast across the thermal boundary layer at the base of the stagnant lid may therefore be characterised by a rheological temperature scale
\begin{equation} T_\mu=\left[-\mu/\frac{{\rm d}\mu}{{\rm d} T}\right]_{\bar{T}}, \end{equation}
where 
$\mu$ and 
$T$ denote viscosity and temperature, respectively. The rheological temperature is evaluated at the temperature of the mixed interior, 
$\bar {T}$.
For large viscosity contrasts, and hence temperature differences much greater than 
$T_\mu$, the convective heat flux associated with the combination of the stagnant lid and the flow of the interior is significantly lower than would be expected for a constant-viscosity fluid at the temperature of the well-mixed interior (Davaille & Jaupart Reference Davaille and Jaupart1993; Solomatov Reference Solomatov1995). Instead, the convective heat flux 
$F_H$ is set by the viscosity of the well-mixed interior and the much smaller rheological temperature scale
\begin{equation} F_H\sim\mu^{{-}1/3}T_\mu^{4/3}=\mu\left(\frac{{\rm d}\mu}{{\rm d} T}\right)^{{-}4/3}. \end{equation}A planet in a thermal steady state exhibits heat fluxes which must balance across the unstable boundary layer, so that in a steady state the convective heat flux from the mixed interior is equal to the conductive heat flux into the stagnant lid. Since the planet evolves slowly, the convective heat flux may be taken to be quasi-steady, and the thermal evolution given by a global heat balance. In such a setting the crust may then form by extraction of buoyant liquid from the partially molten mantle and compaction in the stagnant lid (Michaut & Neufeld Reference Michaut and Neufeld2022).
It is worth noting that the phase diagram as characterised by solidus and liquidus of silicate rocks varies considerably over the large pressure scales found in planetary interiors (Maurice et al. Reference Maurice, Tosi, Samuel, Plesa, Hüttig and Breuer2017), and accordingly there is a large pressure dependence of the viscosity of well-mixed partial melts (Elkins-Tanton Reference Elkins-Tanton2012). These combined effects suggest that the pressure-dependent rheology may have a significant impact on the evolution of planetary surfaces and the global heat budget of their interiors. Pressure- and temperature-dependent viscosities give rise to pressure- and temperature-dependent rheological temperature scales, which are reflected in both a temperature and pressure, or equivalently depth, dependence of the convective heat fluxes. Since the heat flux is independent of the depth of the mixed interior in classical stagnant-lid theory (Davaille & Jaupart Reference Davaille and Jaupart1993), it is the pressure at the top of the mixed interior that is used in the calculations of the viscosity and the corresponding heat flux. Since the pressure within the lid is hydrostatic, the heat flux may be considered to be a function of mixed interior temperature and stagnant-lid thickness. Thicker stagnant lids will lead to a higher pressure and hence higher viscosity, giving rise to the potential for a reduction in convective heat flux. This may in turn lead to further thickening of the stagnant lid. Conversely, for thinner stagnant lids the lower pressures may imply a lower viscosity, and hence higher heat flux reducing stagnant-lid growth.
Here, we investigate the evolution and stability of the growth of a planetary stagnant lid whose rheology is both temperature and pressure dependent. The close association between stagnant lids and crustal formation (Breuer & Spohn Reference Breuer and Spohn2006; Morschhauser, Grott & Breuer Reference Morschhauser, Grott and Breuer2011; Thiriet et al. Reference Thiriet, Michaut, Breuer and Plesa2018; Michaut & Neufeld Reference Michaut and Neufeld2022) suggests that this instability may play an important role in amplifying small initial perturbations to produce the large crustal dichotomy seen in some planetary bodies. In § 2, we begin by describing in detail the physical model, detailing the effect of both the temperature and pressure dependence of the phase diagram on the rheology of a solidifying magma ocean. We use boundary-layer arguments to constrain the dependence of the convective heat flux on temperature and pressure, and hence depth. In addition, we enforce energy conservation across the unstable boundary layer to determine the time-dependent growth of the stagnant lid. Subsequently, in § 3 we present the results of our model for the symmetric growth of a uniform planetary stagnant lid. In order to provide a concrete example, we numerically examine a radially symmetric model of the Moon, and then use this insight to set up a generic symmetric stagnant lid chosen to be representative to use as a base state for our subsequent perturbation analysis. In § 4, we consider the possible growth of perturbations to this symmetric stagnant lid, deriving a condition for their exponential growth throughout the evolution of the stagnant lid. We find that pressure dependence of the melt fraction, and hence rheology and convective heat flux, leads to an instability that grows most quickly at large wavelengths. Finally, in § 5 we discuss the implications of our predictions for observations of crustal dichotomy on planetary bodies, and discuss for which planetary bodies our analysis is relevant.
2. Description of the physical model
2.1. Overview
We consider the cooling and crystallisation of a planetary magma ocean from an initial state with a high melt fraction. Surface radiation or atmospheric convection drives cooling and crystallisation of a surface boundary layer which, because of the sensitive dependence of the mixture viscosity on temperature and solid fraction, results in a stagnant lid of thickness 
$h(t)$, which evolves as a function of time 
$t$ (see figure 1). To simplify the subsequent analysis, and to focus on the stability of the evolving stagnant lid we approximate surface radiative cooling with a fixed, cold surface temperature 
$T_s$. As the lid is stagnant, the thermal evolution of the lid is determined by diffusion and radiogenic heating alone.
Figure 1. Diagram of model set-up, including an asymmetric planetary structure and a close-up of the structure of the stagnant lid and unstable boundary layer.
Rapid convection ensures the interior of the magma ocean is well mixed, as assumed in the work of Michaut & Neufeld (Reference Michaut and Neufeld2022). However, we add the additional complicating factor of a pressure- and hence depth-dependent solidus and liquidus, similarly to Elkins-Tanton (Reference Elkins-Tanton2012). Vigorous convection ensures that the crystals within the magma ocean remain in suspension, and the slurry remains in thermodynamic equilibrium. Thus the solid fraction 
$\phi =\phi (T,p)$ is a function of temperature and pressure – for example, decreasing pressure reduces the solidus and liquidus, causing a decrease in 
$\phi$.
In many studies (Elkins-Tanton Reference Elkins-Tanton2012; Laneuville et al. Reference Laneuville, Wieczorek, Breuer and Tosi2013; Tosi et al. Reference Tosi, Grott, Plesa and Breuer2013; Solomatov Reference Solomatov2015; Maurice et al. Reference Maurice, Tosi, Samuel, Plesa, Hüttig and Breuer2017) the temperature of the convecting interior of rocky planets are assumed to follow an adiabatic, or isentropic, temperature profile. Such a temperature profile may be defined in terms of pressure 
$p$ and potential temperature 
$T_p$. We define an adiabatic temperature profile 
$T_i(T_p,p)$ by
\begin{equation} T_i(T_p,0)=T_p, \end{equation}and
\begin{equation} \frac{\partial T_i}{\partial p}=\frac{\alpha T_i}{\rho c_p}. \end{equation}Hence
\begin{equation} T_i=T_p\exp\left[\int_0^p\frac{\alpha(T_i(T_p,p'),p')} {c_p(T_i(T_p,p'),p')\rho(T_i(T_p,p'),p')}{\rm d} p'\right]. \end{equation}
Here, 
$p$ is the pressure, 
$\alpha (T,p)$ is the thermal expansivity of the fluid, 
$\rho$ is its density and 
$c_p(T,p)$ is its constant-pressure specific heat capacity. The density and thermal expansivity of the fluid are related by
\begin{equation} \alpha=-\frac{1}{\rho_0}\frac{{\rm d}\rho}{{\rm d} T}, \end{equation}
where 
$\rho _0$ is the reference density of the fluid, and we use the Boussinesq approximation
\begin{equation} |\rho-\rho_0|\ll\rho_0. \end{equation}
The mixed interior's potential temperature 
$T_p(t)$ is therefore uniform in space, but may change in time to account for heat conservation. Together with the assumption that pressure is hydrostatic to leading order, the pressure, temperature and solid fraction at any point in the modelled planet may be calculated from the phase diagram, the potential temperature and the radial coordinate. The mixed interior's potential temperature may be increased by mechanisms such as radiogenic heating and heat transfer from the planetary core, and decreased by heat lost to the stagnant lid.
Any difference between the convective heat flux into the unstable boundary layer and the conductive heat flux up into the base of the stagnant lid is accounted for by a Stefan-like growth condition in order to ensure heat conservation at the base of the stagnant lid (see (2.27)). This growth is regulated by the temperature difference 
$\delta T$ across the unstable boundary layer at the interface between the stagnant lid and the interior.
As a further simplification we assume that all our thermodynamic parameters remain constant, with the exception of adjustments to the heat capacity 
$c_p$ and thermal expansivity 
$\alpha$ in the partial melt to account for the release of latent heat and the density difference between solid and liquid phases. In reality all properties of silicates may be temperature and pressure dependent.
2.2. Viscosity of the crystal-magma slurry
The effective viscosity of a fluid with a dense suspension of particles (or crystals) is very complex. Studies such as Einstein (Reference Einstein1906, Reference Einstein1911) suggest a linear relationship between solid fraction and viscosity at low solid fraction, whereas those such as Roscoe (Reference Roscoe1952) suggest a viscosity that diverges at a critical solid fraction. For thermodynamically reactive systems the solid fraction must itself be calculated, and for simplicity here we consider it to be determined by a linear relationship between the solidus and liquidus.
Partial melts in silicate rock can persist over a wide range of temperatures and pressures. For example, Maurice et al. (Reference Maurice, Tosi, Samuel, Plesa, Hüttig and Breuer2017) predict a difference of over 500 K between the solidus and liquidus temperatures at zero pressure in peridotite. In linearised form their equations are
$$\begin{gather} T_{sol}=1400\,\text{K}+149.5\,\text{K}\,\text{GPa}^{{-}1}p, \end{gather}$$
$$\begin{gather}T_{liq}=1977\,\text{K}+64.1\,\text{K}\,\text{GPa}^{{-}1}p, \end{gather}$$
where 
$T_{sol}$ and 
$T_{liq}$ are the solidus and liquidus temperatures, respectively.
In order to take into account the finite viscosity of the pure solid, we use the results of Costa, Caricchi & Bagdassarov (Reference Costa, Caricchi and Bagdassarov2009) for partially solid systems. This slurry viscosity, together with an Arrhenius law for the viscosity of the solid and liquid phases, results in a continuous viscosity function
\begin{equation}
\mu(T,p)=\left\{\begin{array}{@{}ll} \mu_s(T,p) & T\leq
T_{sol}(p),\\ \mu_l(T,p)f(\phi) & T_{sol}(p)< T<
T_{liq}(p),\\ \mu_l(T_p) & T\geq T_{liq}l(p).
\end{array}\right.
\end{equation}
Here, the solidus and liquidus are functions of pressure, 
$T_{sol}(p)$ and 
$T_{liq}(p)$, respectively, as given by (2.6) and (2.7), and the viscosity of the solid matrix and interstitial melt are
$$\begin{gather} \mu_s=\mu_{solidus}\exp\left[\frac{E_\mu}{R_g}\left(\frac{1}{T}-\frac{1}{T_{sol}(p)}\right)\right], \end{gather}$$
$$\begin{gather}\mu_l=\mu_{liquidus}\exp\left[\frac{E_\mu}{R_g}\left(\frac{1}{T}-\frac{1}{T_{liq}(p)}\right)\right], \end{gather}$$
respectively, where 
$E_\mu$ is the activation energy, 
$R_g$ is the ideal gas constant and 
$\mu _{liquidus}$ and 
$\mu _{solidus}$ are set as constants for simplicity. In this paper we assume 
$E_\mu$ is constant, though in reality it would depend on both temperature and pressure.
A more realistic model would also take into account the activation volume for creep to find the viscosity in the solid, this being the value 
$V_a$ in a viscosity equation of the form
\begin{equation} \mu_s\sim\exp\left(\frac{E_\mu+V_ap}{R_gT}\right). \end{equation}However, in silicate rock this volume itself varies significantly with temperature and pressure (Borch & Green Reference Borch and Green II1987; Durham et al. Reference Durham, Mei, Kohlstedt, Wang and Dixon2009), and our approximation using the solidus temperature is not unheard of in the field of geophysics (Weertman Reference Weertman1970, Reference Weertman1978).
At temperatures between the solidus and liquidus, we use a model of suspension rheology by Costa et al. (Reference Costa, Caricchi and Bagdassarov2009), which is described by the empirical function
\begin{equation} f(\phi)=\frac{1+(\phi/\phi_{\ast})^\delta}{\left[\xi+(1-\xi) \text{erfc}\left(\dfrac{\sqrt{\rm \pi}(\phi/\phi_{\ast}) \left(1+(\phi/\phi_{\ast})^\gamma\right)} {2(1-\xi)}\right)\right]^{B_E\phi_{\ast}}}, \end{equation}
where 
$\phi$ is the solid fraction, 
$\phi_{\ast}$ is the critical solid fraction, 
$\delta$, 
$\gamma$ and 
$\xi$ are empirical constants and 
$B_E=\frac {5}{2}$ is the Einstein constant (Einstein Reference Einstein1906, Reference Einstein1911). Costa et al. (Reference Costa, Caricchi and Bagdassarov2009) suggests that 
$\delta +\gamma =13$, and we use 
$\gamma =4$ throughout this work. We define 
$\xi$ by the requirement that the viscosity is continuous at solid fraction 
$\phi =1$, so that
\begin{equation} \mu_l(T_{sol}(p),p)f(1)=\mu_{solidus}. \end{equation}While this viscosity function's complicated nature stands in contrast to the simple nature of the functions we use in the pure solid and pure melt, it holds the advantage of being one of the few differentiable viscosity functions for suspensions that cover the entire range of phases from pure solid to pure fluid.
The viscosity is plotted in figure 2(a) as a function of the potential temperature for three representative pressures, and in figure 2(b) as a function of pressure for three potential temperatures. We note that the viscosity increases by several orders of magnitude around the critical solid fraction, but remains bounded above by the viscosity of the pure solid. Notably, the viscosity function is also strongly pressure dependent around this melt fraction, due to the dependence of the solid fraction on pressure through the phase diagram.
Figure 2. Viscosity, heat flux and rheological temperature scales as functions of potential temperature and pressure.
2.3. Modelling the convective heat flux
The viscosity of a mixture of rock and magma is highly dependent on temperature, pressure and solid fraction, with lower temperatures and higher pressures leading to a higher mixture viscosity. Thus, for a sufficiently steep temperature gradient, the viscosity of the hot interior of a planet is significantly less than that at the surface. Despite the large buoyancy forces associated with surface cooling, in this regime the surface is very cold and viscous, and therefore stagnant. However, a region on the lower edge of the stagnant lid remains relatively warm and low viscosity, and is therefore unstable and can engage with the convection of the fluid below.
In a steady state, the convective heat flux into the bottom of the stagnant lid is balanced by heat conduction upwards through the stagnant lid. Davaille & Jaupart (Reference Davaille and Jaupart1993) found that in this case the convective heat flux obeys
\begin{equation} F_H\approx Ck\left(\frac{g\alpha\rho}{\mu\kappa}\right)^{1/3} T_\mu^{4/3}, \end{equation}
when the viscosity depends solely on temperature. Here, 
$C\approx 0.47$ is a dimensionless constant, 
$k$ is the thermal conductivity of the fluid, 
$\mu$ is its dynamic viscosity and 
$\kappa$ its thermal diffusivity and 
$g$ is the gravitational field strength (which itself may vary with depth), all evaluated at the temperature of the mixed interior. The pressure dependence of the partial melt's viscosity introduces a pressure dependence, and hence a depth dependence, to the heat flux. We incorporate this depth dependence of the heat flux into the stagnant lid by evaluating 
$T_\mu$ at the hydrostatic pressure at the base of the stagnant lid. We denote the temperature at the upper boundary of the mixed interior as
\begin{equation} \bar{T}=T_i(T_p,p({-}h)). \end{equation}Thus the pressure- and temperature-dependent rheological temperature scale may be defined as
\begin{equation} T_\mu=\frac{\mu(\bar{T},p({-}h))}{-\left.\dfrac{\partial\mu}{\partial T}\right|_{\bar{T},p({-}h)}}, \end{equation}
and hence the heat flux is naturally dependent on the stagnant lid thickness 
$h$.
The temperature and pressure vary across the unstable boundary layer, and the viscosity of the slurry may likewise vary. It is this variation in viscosity which sets the depth of the unstable boundary layer, and hence the resultant heat flux. Here, we review the scaling for the boundary-layer structure and convective heat flux for isoviscous convection before extending the analysis to temperature dependent rheologies, covering the experimental results of Davaille & Jaupart (Reference Davaille and Jaupart1993), then examining the pressure-dependent case.
For isoviscous convection, the classical mathematical picture at large Rayleigh number is as follows. At the cooled upper boundary (or equivalently a heated lower one) the behaviour of the flow adjacent to the boundary layer is independent of the thickness of the convection cell. The non-dimensional convective heat flux, or Nusselt number, obeys
\begin{equation} Nu=\frac{F_H d}{k{\rm \Delta} T}=f(Ra,Pr). \end{equation}Here
\begin{equation} Ra=\frac{\rho g \alpha {\rm \Delta} T d^3}{\mu \kappa}, \end{equation}
is the Rayleigh number, 
$d$ is the height of the convection cell, 
${\rm \Delta} T$ is the temperature difference across the convection cell and
\begin{equation} Pr=\frac{\mu}{\rho\kappa}, \end{equation}
is the Prandtl number. If the heat flux 
$F_H$ is to be independent of 
$d$, it follows that
\begin{equation} F_H\sim\frac{k{\rm \Delta} T}{d}Ra^{1/3}. \end{equation}When the viscosity of the fluid is temperature dependent, the extent of the unstable boundary layer is additionally set by the variation in viscosity across the boundary layer. Davaille & Jaupart (Reference Davaille and Jaupart1993) studied the convection of a temperature-dependent fluid theoretically and experimentally, and found that the heat flux in the presence of a stagnant lid is well approximated by
\begin{equation} F_H=C\frac{kT_\mu}{d}Ra_\nu^{1/3}=Ck \left(\frac{\rho g\alpha}{\mu\kappa}\right)^{1/3}T_\mu^{4/3}=C \frac{k{\rm \Delta} T}{d}Ra^{1/3}\left(\frac{{\rm \Delta} T}{T_\mu}\right)^{{-4}/{3}}, \end{equation}
for a constant 
$C$, with 
$Ra_\nu$ denoting the Rayleigh number with 
${\rm \Delta} T$ replaced by 
$T_\mu$.
Davaille & Jaupart (Reference Davaille and Jaupart1993) found the experimental value 
$C\approx 0.47\pm 0.03$. Here, we also consider a temperature drop across the unstable boundary layer given by 
$BT_\mu$. Davaille & Jaupart (Reference Davaille and Jaupart1993) found that 
$B\approx 2$, although the numerical work of Thiriet et al. (Reference Thiriet, Breuer, Michaut and Plesa2019) found a value of 
$B\approx 2.54$.
We thus consider the temperature of the base of the stagnant lid to be the value 
$T_u$ which satisfies
\begin{equation} \mu(T_u,p({-}h))={\rm e}^B\mu(\bar{T},p({-}h)), \end{equation}and set the rheological temperature scale
\begin{equation} T_\mu=\frac{\bar{T}-T_u}{B}. \end{equation}Thus the convective heat flux, including both the temperature and pressure dependence of the rheology, is given by
\begin{equation} F_H=Ck\left(\frac{\rho_0 g\alpha}{\mu(\bar{T},p({-}h))\kappa}\right)^{1/3} \left(\frac{\bar{T}-T_u}{B}\right)^{4/3}. \end{equation}Finally, conduction down the adiabatic temperature gradient must be incorporated, so the overall heat flux is
\begin{equation} F_H=Ck\left(\frac{\rho_0 g\alpha}{\mu(\bar{T},p({-}h))\kappa}\right)^{1/3} \left(\frac{\bar{T}-T_u}{B}\right)^{4/3}+\frac{gk\alpha\bar{T}}{cp}. \end{equation} A heat flux function corresponding to the viscosity function depicted in figures 2(a) and 2(b) is plotted in figures 2(c) and 2(d). Additionally, the corresponding values of 
$T_\mu$ are plotted in figures 2(e) and 2( f). Note that 
$T_\mu$ reaches a local minimum where 
$\mu$ changes most rapidly as a function of 
$T$, and hence so does 
$F_H$.
2.4. Heat conservation at the base of the stagnant lid
At the interface between the stagnant lid and mixed interior, two heat fluxes are involved, the convective heat flux from the planetary interior and the conductive heat flux through the stagnant lid. The convective heat flux 
$F_H$ has been described above and its behaviour is summarised in (2.25). In a steady state, this heat flux from the magma ocean is balanced by a conductive heat flux through the stagnant lid given by
\begin{equation} F_{cond}={-}k\left.\frac{\partial T}{\partial z}\right|_{z={-}h}, \end{equation}
where 
$k$ is the thermal conductivity at the base of the stagnant lid and 
$z$ is the local vertical coordinate. In contrast, during transient growth of the stagnant lid, conservation of energy at the interface between the stagnant lid and the mixed interior suggests that the specific heat capacity of the boundary layer cannot be neglected. In this case, the evolution of the stagnant lid is given by
\begin{equation} {\rm \Delta} H \frac{{\rm d} h}{{\rm d} t} = F_{cond} - F_H, \end{equation}(Schubert, Cassen & Young Reference Schubert, Cassen and Young1979; Breuer, Spohn & Wüllner Reference Breuer, Spohn and Wüllner1993) where the difference in specific enthalpy between the mixed interior and the bottom of the stagnant lid is given by
\begin{equation} {\rm \Delta} H = \int_{T_u}^{\bar{T}} \rho \left[ c_p - L \frac{{\rm d}\phi}{{\rm d} T} \right] {\rm d} T. \end{equation}
Here, 
$L$ is the latent heat and 
$\phi$ is the solid fraction. Given that this latent heat is applicable over a range of temperatures between the solidus and the liquidus, we may absorb it into an effective heat capacity
\begin{equation}
c_p(T,p)=\left\{\begin{array}{@{}ll}
c_p+L/(T_{liq}(p)-T_{sol}(p)) & T_{sol}(p)< T< T_{liq}(p)\\
c_p & \text{otherwise} \end{array}\right.
\end{equation}2.5. Heat conservation in the mixed interior
In addition to a statement of energy conservation at the base of the stagnant lid, heat must be conserved in the mixed interior. The secular variation of the heat in the mixed interior is determined by the heat flux across the core–mantle boundary and from radiogenic heating throughout the volume. These heat sources are balanced by removal by the heat flux into the stagnant lid 
$F_H$. Heat conservation in the mixed interior is therefore given by
\begin{equation} C_m(\bar{T},h) \frac{{\rm d}\bar{T}}{{\rm d} t} = F_{CMB} + Q(t) \rho\hat{V}(h) -F_H \hat{A}(h), \end{equation}where the effective heat capacity of the entire mixed interior is given by
\begin{equation} C_m=\underset{V_{interior}}{\int}\rho c_p(T,p)\frac{\partial T}{\partial\bar{T}}{\rm d} V, \end{equation}
and where 
$F_{CMB}$ is the total heat flux across the core–mantle boundary, 
$Q(t)$ is the (spatially uniform) rate of radiogenic heating per unit mass as a function of time, 
$\hat {V}(h)$ is the volume of the mixed interior and 
$\hat {A}(h)$ is the area of the interface between the mixed interior and the stagnant lid. The temperature 
$T$ at any given point in the calculation of 
$C_m$ is constrained by being on the same adiabat as the temperature 
$\bar {T}$ at radius 
$R-h$.
For the sake of simplicity, and motivated by the small size of the lunar core, in the following analysis we set 
$F_{CMB}=0$. Similarly, we take the rate of radiogenic heating to be of the form
\begin{equation} Q(t)=Q_0\times 2^{{-}t/t_{half}}, \end{equation}
for an initial rate of radiogenic heating 
$Q_0$ per unit mass and a half-life 
$t_{half}$, which we set as 
$t_{half}=1\,\text {Ga}$ (
$10^9$ years) to give a decay rate of the correct order of magnitude for the Earth's moon.
3. Symmetric growth of the stagnant lid
3.1. Background: symmetric computation
Working from (2.27) and (2.30), we numerically calculate the evolution of a spherically symmetric planetary model, in order to illustrate and understand the effect of the pressure-dependent rheology, and hence heat flux, on stagnant-lid growth. We base our calculation on an approximation of parameters relevant to Earth's moon. Table 1 shows the parameters used in our calculations.
Table 1. Parameters used in calculation.
We calculate the evolution of a representative series of stagnant lid growth curves as follows. A high initial temperature is chosen, with a thin stagnant-lid thickness taken that allows for a linear temperature profile with no stagnant-lid growth, so that
\begin{equation} \frac{k(T_u-T_s)}{h}=F_H. \end{equation}
We start from a mixed interior potential temperature of 1800 K, which allows the stagnant lid to start from a negligible thickness of order 1 m. The surface temperature is fixed at 
$T_s=250\,\text {K}$, and the stagnant lid begins with a linear temperature profile so that the conductive heat flux at its lower edge matches the convective heat flux. This state then evolves according to (2.27) and (2.30), with thermal conduction and radiogenic heating in the stagnant lid given by
\begin{equation} \rho c_p\frac{\partial T}{\partial t}=\boldsymbol{\nabla}\boldsymbol{\cdot}(k\boldsymbol{\nabla} T)+Q\rho. \end{equation}
We use fifth-order finite differences to calculate temperature derivatives in the stagnant lid, resulting in a system of ordinary differential equations that we solve using a backwards differentiation method. The results of our simulations of stagnant-lid growth and planetary thermal evolution are shown in figure 3. In figure 3(a) we plot the potential temperature 
$T_p(t)$, which characterises the cooling of the planetary interior, and in figure 3(b) we plot the stagnant-lid thickness, which characterises the surface boundary layer. In these simulations, the value of 
$Q_0$ was changed between calculations, with all other parameters held constant. At early times, 
${\sim }0\unicode{x2013}10\,\text {Ma}$ (millions of years), we find that the high initial interior temperatures lead to a high convective heat flux, which causes rapid cooling of the mixed interior. This causes the convective heat flux to decrease, resulting in rapid growth of the stagnant lid. As the stagnant lid is very thin at this point, the time scale for diffusion across it remains small, so it retains an approximately linear temperature profile with a close balance between convective and conductive heat fluxes. Thus
\begin{equation} h\approx \frac{k(T_u-T_s)}{F_H}. \end{equation}
Figure 3. Symmetric calculation of lunar interior potential temperature (a) and stagnant-lid thickness (b). A clear separation between a low 
$Q_0$, low heat flux regime and a high 
$Q_0$, high heat flux regime is visible.
After 
${\sim }10\,\text {Ma}$, radiogenic heating in the interior balances the rapid surface loss of heat, and so for intermediate times the interior temperature is quasi-static. Beyond 100 Ma, our solutions diverge into two main groups. At higher rates of radiogenic heating, the solutions come close to a balance between radiogenic heating and convective heat flux (see figure 4(b) for an example), with radioactive decay preventing a true steady state from being reached. At lower levels of heating, no such balance is approached (see figure 4(a) for an example), and solutions tend towards having the upper boundary of the mixed interior at the critical solid fraction where viscosity changes most rapidly as a function of temperature. This results in a low rheological temperature scale and hence a low convective heat flux. Stagnant-lid growth in this case is thus primarily controlled by the enthalpy difference across the unstable boundary layer, in an analogous manner to the classical Stefan problem with no convective heat flux.
Figure 4. Overall heat fluxes: 
$Q_0= \ (\textit {a})\ 1.0\times 10^{-11}\,\text {W}\,\text {kg}^{-1}$ and (b) 
$3.0\times 10^{-11}\,\text {W}\,{\rm kg}^{-1}$.
However, for all the values we use, the Moon spends much of its early evolution undergoing quasi-steady changes. This evolution is seen in the growth of the stagnant lid and an increase in the temperature of the mixed interior, with a thinner stagnant lid forming at larger initial radiogenic heating rates. We note that, on a plot of mixed interior potential temperature against stagnant-lid thickness, our solutions tend towards one of two linear regions, approximately following contours of the convective heat flux, as shown in figure 5. For a smaller rate of radiogenic heating (
$Q_0=1.0\times 10^{-11}\,\text {W}\,\text {kg}^{-1}$), this contour lies in a region where the convective heat flux increases with stagnant-lid thickness. However, the other contour for larger rates of radiogenic heating (
$Q_0=3.0\times 10^{-11}\,\text {W}\,\text {kg}^{-1}$) lies in a region where this heat flux decreases with stagnant lid thickness, leading to the possibility of instability and asymmetric stagnant-lid growth. We return to the question of the stability of the stagnant lid in § 4 and instead characterise these modes of symmetric growth as a base state to be assessed.
Figure 5. Symmetric calculation in stagnant-lid thickness–potential temperature space. Ticks on the plots are every 100 Ma. These calculations were run over a time scale of 1 Ga.
Given the important role that the viscosity plays in the heat flux and therefore the thermal evolution of the planetary interior, we also plot the viscosity as a function of time and depth for two simulations (see figure 6). We therefore incorporate the slow growth of the stagnant lid and change in mixed interior temperature into the symmetric base state of our subsequent stability analysis.
Figure 6. Plots of the structure of the viscosity with depth for differing radiogenic heating rates, (a) 
$Q_0=1.0\times 10^{-11}\,\text {W}\,\text {kg}^{-1}$, and (b) 
$Q_0=3.0\times 10^{-11}\,\text {W}\,\text {kg}^{-1}$. In both plots the dashed lines mark the lower boundary of the stagnant lid, and viscosities outside the range of the colour bar are not shown.
3.2. Quasi-steady evolution of a planar symmetric state and the potential for instability
The growth of the stagnant lid is seen to slowly evolve, driven by the conductive cooling from the surface and limited by the convective heat flux from the planetary interior. That this convective heat flux is itself pressure, and hence depth, dependent raises an important possibility. If the convective heat flux decreases more rapidly with depth than the conductive cooling of a thicker stagnant lid then the growth of the stagnant lid may be unstable to perturbations from the symmetric base state. Here, we first describe, in non-dimensional form, the evolution of the symmetric base state before considering perturbations to this initially symmetric evolution in the following section.
The results of the calculations in § 3.1 suggest that we may consider the evolution of the stagnant lid and planetary interior to be approximately quasi-steady. Furthermore, we note that the base state and the following stability analysis are not specific to the Moon, but are generic to any spherical rocky planetary body with pressure- and temperature-dependent viscosity in its convecting interior, provided it exhibits a conductive stagnant lid.
We assume that the stagnant-lid thickness 
$h(t)$ does not vary spatially, but grows at a slow rate 
$dh/dt=V_0$. Similarly, the potential temperature 
$T_p(t)$ of the convecting interior is taken to be a function of time but not space. The mantle mixed interior temperature 
$\bar {T}$ likewise varies slowly at the rate 
$\dot {\bar {T}}={\rm d}\bar {T}/{\rm d} t$. For the purposes of our stability analysis, the calculation of 
$\dot {\bar {T}}$, as detailed in (2.30), is not as important as the resulting value. For simplicity, we assume a quasi-steady stagnant-lid base-state temperature profile 
$T_0(z,t)$ of the form
\begin{equation} T_0=T_s+{\rm \Delta} T \varTheta_0(\eta), \end{equation}
where 
$T_s$ is the fixed surface temperature at 
$z=0$, 
${\rm \Delta} T=\bar {T}-\delta T-T_s$ and 
$\varTheta _0$ is the non-dimensional base-state temperature profile, which is a function of 
$\eta =-z/h(t)$. As in § 2.1, 
$\delta T$ is the temperature difference across the unstable boundary layer. Including the effects of thermal conduction and radiogenic heat production, the temperature in the stagnant lid obeys
\begin{equation} \rho c_p \frac{\partial T_0}{\partial t}=\boldsymbol{\nabla} \boldsymbol{\cdot}(k\boldsymbol{\nabla} T_0)+Q\rho, \end{equation}
where 
$k=\rho c_p\kappa$ is the thermal conductivity, and 
$\kappa$ is the thermal diffusivity. If 
$k$ is constant and we consider first a simpler, planar geometry, then (3.4) and (3.5) may be written non-dimensionally as
\begin{equation}
\frac{{\rm
d}^2\varTheta_0}{{\rm d}\eta^2}+Pe
\eta\frac{{\rm d}\varTheta_0}{{\rm d}\eta}-\tau
\varTheta_0+\tilde{Q}=0,\quad
\varTheta_0(0)=0,\quad\varTheta_0(1)=1,
\end{equation}where
\begin{equation} Pe=hV_0/\kappa, \end{equation}is the Péclet number, which characterises the speed of lid growth compared with thermal diffusion. The non-dimensional rate of change of the mixed interior's potential temperature is given by
\begin{equation} \tau=\dot{\bar{T}}h^2/(\kappa {\rm \Delta} T), \end{equation}and the non-dimensional rate of radiogenic heating is given by
\begin{equation} \tilde{Q}=Qh^2/(c_p\kappa{\rm \Delta} T). \end{equation}
In general, (3.6a–c) cannot be solved analytically, but we plot numerically calculated profiles of 
$\varTheta _0(\eta )$ in figure 7 for representative values of 
${Pe}$, 
$\tau$ and 
$\tilde {Q}$. We find that increasing 
${Pe}$ or, 
$\tilde {Q}$, or decreasing 
$\tau$, gives rise to a negative second derivative of 
$\varTheta _0$, and correspondingly a smaller value of 
$\varTheta _0'(1)$. These curves represent non-dimensional temperature profiles in the stagnant lid under various conditions. Their slope at 
$\eta =1$ is proportional to the conductive heat flux into the base of the stagnant lid, and so is important to understanding its stability.
Figure 7. Example plots of the non-dimensional thermal structure of the stagnant lid in planar geometry, 
$\varTheta _0(\eta )$, for differing Péclet numbers 
${Pe}$, rates of mixed interior temperature change 
$\tau$ and radiogenic heating rates 
$\tilde {Q}$.
We also note that energy conservation in the unstable boundary layer may be written in the form of an evolution equation for the growth of the stagnant lid. In dimensional form this is
\begin{equation} \rho c_p \delta T V_0={-}k\left.\frac{\partial T_0}{\partial z}\right|_{z={-}h}-F_H(\bar{T},h), \end{equation}or non-dimensionally
\begin{equation} Pe \frac{\delta T}{{\rm \Delta} T}=\varTheta_0'(1)-\frac{F_H(\bar{T},h)h_0}{k{\rm \Delta} T}. \end{equation}
If the convecting material is of roughly constant composition, the two main variables affecting 
$F_H$ are temperature and pressure. Given that the pressure is hydrostatic to leading order, we can replace pressure with depth below the surface in our calculation of 
$F_H$.
The solution of (3.6a–c) and (3.10) allows us to determine the thermal structure, and hence the heat flux during the symmetric, quasi-steady growth of the stagnant lid. However, the heat flux is both pressure and temperature dependent, and hence may depend on stagnant-lid thickness. In particular, the convective heat flux from the interior may decrease with increasing stagnant-lid thickness and hence the growth rate of the stagnant lid may increase with depth leading to the potential for instability. In what follows, we consider spatial (azimuthal) perturbations to the evolving stagnant-lid thickness.
4. Growth of variations in the stagnant lid
4.1. Planar perturbations
To examine the possibility of instabilities in the growth of the stagnant lid, we consider a perturbation to the quasi-steady symmetric base state with planar geometry described in § 3. This is a simplification of the more realistic spherical geometry covered in § 4.2. Here, we consider perturbations to a uniform thickness 
$h=h_0$ such that 
$h=h(x,t)$ for horizontal coordinate 
$(x)$ and time 
$(t)$. We assume that any changes to the amplitude of perturbations occur quickly compared with changes in the base state. A non-uniform perturbation may be caused by external forcing, such as spatially variable surface heating or impacts, or internally by chemical differentiation or hot plumes in the convecting interior. However, we neglect any long-term local effects of any such plumes, with the convective heat flux remaining a function of 
$h$ and 
$\bar {T}$ only. With this perturbation, energy conservation at the base of the stagnant lid may now be written as
\begin{align} \rho c_p \delta T (V-V_0) &={-}k\left[\left.\frac{\partial}{\partial z}\right|_{z={-}h}(T-T_0)+ \left.\frac{\partial T_0}{\partial z}\right|_{z={-}h}-\left.\frac{\partial T_0}{\partial z}\right|_{z={-}h_0}\right] \nonumber\\ &\quad -\left[F_H(\bar{T},h)-F_H(\bar{T},h_0)\right], \end{align}
where 
$V={\partial h}/{\partial t}$ denotes the stagnant-lid growth rate in the perturbed state. We consider small-amplitude perturbations to the stagnant-lid thickness of the form
\begin{equation} h=h_0[1+\epsilon(t)\cos(ax)], \end{equation}
for wavenumber 
$a>0$ and dimensionless perturbation amplitude 
$\epsilon (t)\ll 1$. A positive value of 
$a$ is chosen to ensure the overall convective heat flux from the interior is unchanged to leading order, avoiding the more complicated case in which perturbations to the stagnant lid significantly affect the interior temperature. As we study linear perturbations here, we may assume that 
$\epsilon =\epsilon _0\exp (\sigma \kappa t/h_0^2)$ for some non-dimensional growth rate 
$\sigma$ and initial amplitude 
$\epsilon _0$.
To first order in 
$\epsilon$ this suggests a stagnant-lid thermal profile determined by the quasi-steady condition
\begin{equation} T=T_s+[\bar{T}(t)-\delta T-T_s]\varTheta(\eta), \end{equation}together with the thermal diffusion equation
\begin{equation} \left.\frac{\partial T}{\partial t}\right|_z =\kappa\nabla^2T+\frac{Q}{c_p}, \end{equation}and boundary conditions
$$\begin{gather} T(0,t)=T_s, \end{gather}$$
$$\begin{gather}T({-}h_0,t)=\bar{T}-\delta T+[h(x)-h_0]\left. \frac{\partial T_0}{\partial z}\right|_{z={-}h_0}. \end{gather}$$Thus the perturbed temperature profile may be written as
\begin{equation} T=T_0(z)-{\rm \Delta} T \varTheta_0'(1)\epsilon(t)\cos(ax) \varTheta_a(\eta), \end{equation}
where again 
$\eta =-z/h_0$, and the vertical structure of the perturbed thermal field, 
$\varTheta _a(\eta )$, obeys
\begin{equation} \frac{{\rm d}^2\varTheta_a}{{\rm d}\eta^2}+ Pe \ \eta\frac{{\rm d}\varTheta_a}{{\rm d}\eta}- [(ah_0)^2+\tau]\varTheta_a=0, \end{equation}subject to boundary conditions (4.5) and (4.6), which may be written as
$$\begin{gather} \varTheta_a(0)=0, \end{gather}$$
$$\begin{gather}\varTheta_a(1)=1, \end{gather}$$
respectively. Substituting the form of thermal perturbations, (4.7), into a statement of energy conservation, (4.1), we see that the dimensionless growth rate 
$\sigma$ is given by
\begin{equation} \sigma=\frac{h_0^2}{\kappa}\frac{1}{\epsilon}\frac{{\rm d}\epsilon}{{\rm d} t}= \frac{{\rm \Delta} T}{\delta T}\left[\varTheta_0''(1)-\varTheta_0'(1) \varTheta_a'(1) -\frac{h_0^2}{k{\rm \Delta} T}\left.\frac{\partial F_H}{\partial h}\right|_{\bar{T},h_0}\right]. \end{equation}This gives the growth rate of the perturbations as a function of the diffusive thermal profile in the stagnant lid, and the depth dependence of the convective heat flux from the mantle below. Equation (4.11) implies that these modes grow if and only if the heat flux increases with stagnant-lid thickness at a greater rate than the temperature gradient in the stagnant lid,
\begin{equation} \frac{h_0^2}{k{\rm \Delta} T}\left.\frac{\partial F_H}{\partial h}\right|_{\bar{T},h_0}< \varTheta_0''(1)-\varTheta_0'(1) \varTheta_a'(1). \end{equation}4.2. Spherical geometry
To improve the applicability of these results in planetary settings, we recalculate the growth rate of perturbations to the stagnant-lid thickness in a spherical geometry. We let the surface be at radius 
$r=R$, and redefine a spherical 
$\eta =(R-r)/h$.
We again assume a quasi-steady thermal profile given by (3.4) and hence solve for thermal diffusion, given by (3.5) now in spherical coordinates. Along with boundary conditions at the surface and base of the stagnant lid, the thermal perturbations therefore satisfy
$$\begin{gather} \frac{{\rm d}^2\varTheta_0}{{\rm d}\eta^2}+\left[ Pe \eta-\frac{2}{\zeta^{{-}1}-\eta}\right] \frac{{\rm d}\varTheta_0}{{\rm d}\eta}-\tau\varTheta_0+\tilde{Q}=0, \end{gather}$$
$$\begin{gather}\varTheta_0(0)=0, \end{gather}$$
$$\begin{gather}\varTheta_0(1)=1. \end{gather}$$
The primary effect of this coordinate change is to introduce the additional 
$( (2/r) (\partial T / \partial r))$ term that must be added to the Laplacian of the temperature field, in order to account for the curvature inherent in spherical coordinates. Here, 
$\zeta =h_0/R$ is the ratio of the stagnant lid's thickness to the planetary radius, so that the limit 
$\zeta \rightarrow 0$ recovers the planar case.
We consider perturbations to the uniform stagnant-lid thickness 
$h_0$, in the form of a spherical harmonic function of degree 
$l\in \mathbb {N}$,
\begin{equation} h=h_0\left[1+\epsilon_0 \exp({\sigma\kappa t/h_0^2}) Y_{lm}(\theta,\varphi)\right]. \end{equation}
As with the planar case, this spherical wavenumber 
$l$ is non-zero in order to have an unchanged mean stagnant-lid thickness, so that we cover cases in which the interior temperature is not affected to leading order in 
$\epsilon$ by perturbations to the stagnant lid's thickness. Equation (4.8) can thus be written, with appropriate boundary conditions, as
$$\begin{gather} \frac{{\rm d}^2\varTheta_l}{{\rm d}\eta^2}+\left[ Pe \ \eta- \frac{2}{\zeta^{{-}1}-\eta}\right]\frac{{\rm d}\varTheta_l}{{\rm d}\eta}- \left[\tau+\frac{l(l+1)}{\left(\zeta^{{-}1}-\eta\right)^2}\right] \varTheta_l=0, \end{gather}$$
$$\begin{gather}\varTheta_l(0)=0, \end{gather}$$
$$\begin{gather}\varTheta_l(1)=1. \end{gather}$$
Note that, in the limit of small 
$\zeta$ and large 
$l$, we recover (4.8) with 
$ah_0$ replaced by 
$l\zeta$.
In figure 8, we show a comparison of temperature profiles derived from our numerical calculations in § 3.1 with quasi-steady results from a calculation of the base state in spherical geometry. While there remains a difference between these two profiles, the quasi-steady approximation provides a significant improvement on a linear temperature profile in both cases. The instability condition (4.12) and growth rate equation (4.11) hold irrespective of whether the temperature profiles are truly quasi-steady – the quasi-steady calculations of the temperature profiles in (3.6a–c), (4.8), (4.13) and (4.17) are simply used to give values for the derivatives of 
$\varTheta$ at 
$\eta =1$.
Figure 8. Comparison of temperature profiles from numerical calculations with those from the equivalent quasi-steady calculations, at 
$t=250\,\text {Ma}$.
The planar equations, (4.11) and (4.12), describing the evolution of the perturbations and the conditions for growth, may now be written as
\begin{equation} \sigma=\frac{{\rm \Delta} T}{\delta T}\left[\varTheta_0''(1)-\varTheta_0'(1) \varTheta_l'(1) -\frac{h_0^2}{k{\rm \Delta} T}\left.\frac{\partial F_H}{\partial h}\right|_{\bar{T},h_0}\right], \end{equation}
and again we find that the stagnant lid is unstable to the growth of perturbations when 
$(\sigma > 0)$, that is when
\begin{equation} \frac{h_0^2}{k{\rm \Delta} T}\left.\frac{\partial F_H}{\partial h}\right|_{\bar{T},h_0}< \varTheta_0''(1)-\varTheta_0'(1) \varTheta_l'(1). \end{equation}
Again, from (4.21) we see that instabilities grow when the depth dependence of the convective heat flux is less than the depth dependence of the conductive heat loss through the stagnant lid, thus leading to the possibility of enhanced growth. If 
$\varTheta _0''(1)-\varTheta _0'(1)\varTheta _l'(1)$ is increased, then the conductive heat flux depends less strongly on the amplitude of the perturbation, so the instability is possible for a wider range of possible rheologies. The occurrence of this difference in the depth dependence of convective and conductive fluxes having smaller magnitude may be approximated for a quasi-steady temperature profile by rearranging (4.13) and (4.17) at the base of the stagnant lid (
$\eta =1$) from which we find that
$$\begin{gather} \varTheta_0'(1)=\frac{\tau-\tilde{Q}-\varTheta_0''(1)}{ Pe -\dfrac{2}{\zeta^{{-}1}-1}}, \end{gather}$$
$$\begin{gather}\varTheta_l'(1)=\frac{\tau+\dfrac{l(l+1)}{\left(\zeta^{{-}1}-1\right) ^2}-\varTheta_l''(1)}{ Pe -\dfrac{2}{\zeta^{{-}1}-1}}. \end{gather}$$
Changes to 
${Pe}$, 
$\tau$, 
$\tilde {Q}$ and 
$\zeta$ may have complicated effects on 
$\varTheta _0$ and 
$\varTheta _l$, and hence on the growth rate of this instability. However, an increase in the spherical degree 
$l$ will, for a given value of 
$\varTheta _l''(1)$, cause an increase in 
$\varTheta _l'(1)$, resulting in a lower threshold value of 
$({h_0^2}/{k{\rm \Delta} T})({\partial F_H}/{\partial h})$ (see figure 9). Thus small but non-zero spherical wavenumbers promote this instability.
Figure 9. Example plots of 
$\varTheta _0''(1)-\varTheta _0'(1)\varTheta _l'(1)$, for 
$\zeta =0.1$, 
$\tau =0.1$, 
$\tilde {Q}=0.5$. Larger values of the spherical wavenumber 
$l$ give rise to a larger absolute value of 
$\varTheta _0''(1)-\varTheta _0'(1)\varTheta _l'(1)$, and hence result in a greater likelihood of stability via (4.21).
5. Results and discussion
We compare the results of this stability analysis with the series of full numerical calculations described in § 3.1. To make this comparison we must calculate the non-dimensional parameters as functions of values available from the simulations. For given values of the base-state stagnant-lid thickness 
$h_0$, the temperature 
$\bar {T}$ at the top of the mixed interior, and the rate 
$Q$ of radiogenic heating, we calculate the non-dimensional parameters 
$\zeta$ and 
$\tilde {Q}$. We then calculate the parameter 
$\tau$ using (2.30), and solve (3.10) and (4.13) to find the Péclet number 
${Pe}$ and the temperature profile in the stagnant lid, assuming a fixed surface temperature. We thus obtain the value of the Péclet number and interior cooling rate, as plotted in figures 10 and 11, and hence the growth rate of perturbations as plotted in figure 12. A comparison of the trajectories in figure 5 shows that for hundreds of Ma the simulations with 
$Q_0\geq 2.0\times 10^{-11}\,\text {W}\,\text {kg}^{-1}$ sustain positive growth rates for 
$l=1$ and many higher spherical harmonics. Figure 13 shows a selection of such growth rates as a function of time for two numerical simulations, clearly showing the positive growth rates in the case with the higher rate of radiogenic heating. The primary reason for the difference in growth rates between simulations with different rates of radiogenic heating is that those with more heating retain stagnant-lid thicknesses and interior potential temperatures corresponding to the critical solid fraction 
$\phi _{\ast}$, at which suspended crystals begin to lock together. The depth dependence of the convective heat flux, 
$\partial F_H/\partial h$, is correspondingly large in magnitude during this onset. Figure 14 shows a comparison of growth rates for a wide range of wavenumbers. At larger wavenumbers, or smaller wavelengths, growth rates are negative due to the increased effectiveness of conductive cooling through the stagnant lid. Similarly, lower temperatures promote stability by being associated with solid fractions far above the critical value, and hence less strongly temperature- and pressure-dependent viscosities and heat fluxes.
Figure 10. Absolute value of the Péclet number as a function of stagnant-lid thickness and mixed interior potential temperature, calculated using (4.13) and (3.10), for 
$Q=2\times 10^{-11}\,\text {W}\,\text {kg}^{-1}$. The Péclet number is positive to the left of the dashed line and negative to the right.
Figure 11. Absolute value of the rate of change of mixed interior potential temperature 
$\tau$ as a function of stagnant-lid thickness and mixed interior potential temperature, calculated using (4.13), (3.10) and (2.30), for 
$Q=2\times 10^{-11}\,\text {W}\,\text {kg}^{-1}$. Here, 
$\tau$ is positive to the left of the dashed line and negative to the right.
Figure 12. Growth rate for 
$l=1$, for 
$Q=2\times 10^{-11}\,\text {W}\,\text {kg}^{-1}$. The maroon line is the trajectory of the numerical calculation with 
$Q_0=2\times 10^{-11}\,\text {W}\,\text {kg}^{-1}$. The growth rate is positive to the right of the dashed line and negative to the left.
Figure 13. Growth rate as a function of time during the numerical calculation for (a) 
$Q_0=1.0\times 10^{-11}\,\text {W}\,\text {kg}^{-1}$, and (b) 
$Q_0=3.0\times 10^{-11}\,\text {W}\,\text {kg}^{-1}$.
Figure 14. Growth rate as a function of 
$l$, for a range of potential temperatures 
$\bar {T}$, with radiogenic heating rate 
$Q=2\times 10^{-11}\,\text {W}\,\text {kg}^{-1}$, and unperturbed stagnant lid thickness 
$h_0=100\,\text {km}$.
However, these higher wavenumbers do have lower growth rates than 
$l=1$, leading to dominance of a hemispherical dichotomy, at least during the initial growth of the instability.
It is worth noting that we have not calculated the nonlinear growth of these modes, so cannot draw a clear conclusion on the effect of finite amplitude perturbations. However, we do now have a potential mechanism for the growth of a small-scale perturbation to a full hemispheric dichotomy. Additionally, there remains the possibility that the 
$l=1$ mode was initially excited to a much larger amplitude than higher-wavenumber modes by some other mechanism, for example by radiative heat from the Earth, a large impact, or long-term convective modes in the mixed interior. This initial dichotomy would then be amplified by the instability of the stagnant lid.
The most important factor in promoting this instability is having a sufficiently large negative value for the non-dimensional depth (and hence pressure) dependence of the convective heat flux, the left-hand side of the inequality that is the condition for perturbations to grow
\begin{equation} \frac{h_0^2}{k{\rm \Delta} T}\left.\frac{\partial F_H}{\partial h}\right|_{ \bar{T},h_0}< \varTheta_0''(1)-\varTheta_0'(1)\varTheta_l'(1). \end{equation}
This occurs when the solid fraction is around the value that causes the most rapid change in viscosity (see figure 2). Other parameters which may affect the rapidity and likelihood of the growth of the instability are 
$\tau$ corresponding to heating and cooling of the mixed interior, the Péclet number 
${Pe}$ corresponding to growth and shrinkage of the stagnant lid and the monotonic effect of the spherical wavenumber 
$l$ – larger values suppress the instability.
The condition (4.21) may also be expressed in a dimensional form. If we define the quantity 
$T_l$ by the relationship
\begin{equation} T_l(r)\epsilon(t)Y_{lm}(\theta,\varphi)=T(r,\theta,\varphi)-T_0(r), \end{equation}or equivalently
\begin{equation} T_l(r)={\rm \Delta} T\varTheta_0'(1)\varTheta_l\left(\frac{R-r}{h_0}\right), \end{equation}then (4.21) is equivalent to
\begin{equation} \frac{\partial F_H}{\partial h}<{-}\frac{F_{cond}}{h_0}= \frac{k}{h_0}\left[\left.\frac{\partial T_l}{\partial r}\right|_{r=R-h} +\left(\left.\frac{\partial T_0}{\partial r}\right|_{r=R-h} -\left.\frac{\partial T_0}{\partial r}\right|_{r=R-h_0}\right)\right]. \end{equation}While stagnant-lid convection has been studied for decades in the context of solid state mantle convection (e.g. Davaille & Jaupart Reference Davaille and Jaupart1994), our study relies on the application to a stagnant-lid-type model to partial melts around the critical solid fraction. As such, it is particularly applicable to the solidification of a cooling planetary magma ocean. The rheology of the partial melt in such circumstances is highly complex (Costa et al. Reference Costa, Caricchi and Bagdassarov2009; Keller & Suckale Reference Keller and Suckale2019), with considerable uncertainty and likely non-Newtonian behaviour, meaning our model is certainly simplified in this regard. Nevertheless, the instability of stagnant-lid growth for rheologies which are both pressure and temperature dependent is generic, and hence we anticipate that our results will be applicable to similar rheological models.
This model also assumes that solids remain fully entrained during convection. This condition has been considered unlikely in the absence of local entrainment mechanisms (Solomatov & Stevenson Reference Solomatov and Stevenson1993), but these may be provided by turbulent flow. In particular, the authors of Solomatov & Stevenson (Reference Solomatov and Stevenson1993) note in a later paper that crystal entrainment is very likely in magma oceans (Solomatov, Olson & Stevenson Reference Solomatov, Olson and Stevenson1993) due to the high Rayleigh numbers, which may be up to 
$10^{30}$ in the case of Earth's moon (Elkins-Tanton Reference Elkins-Tanton2012).
The model may be further extended with more complicated factors relating to planetary solidification that are not covered here. Radiogenic elements are thought to have been enriched in the last material to solidify in the lunar crust (Wieczorek et al. Reference Wieczorek2006), meaning they were far from a uniform distribution. A move of radioactive material from the interior to the stagnant lid would affect our non-dimensional parameters by reducing the interior heating rate 
$\tau$ and increasing the stagnant-lid radiogenic heating rate 
$\tilde {Q}$, with both effects promoting a stronger concavity of the temperature profile in the stagnant lid. Non-uniform heating within the stagnant lid would lead to a more complicated temperature profile, or in the event of horizontal variations would require an altogether different model, as in Laneuville et al. (Reference Laneuville, Wieczorek, Breuer and Tosi2013). Additionally, strong contrasts in viscosity in the planetary interior may preclude it from being fully mixed, resulting in a separation into multiple mixed shells (Lebrun et al. Reference Lebrun, Massol, Chassefière, Davaille, Marcq, Sarda, Leblanc and Brandeis2013). In this case, the heat flux from below into the outermost mixed shell may be merged into the calculation of 
$\dot {\bar {T}}$, and hence its non-dimensional analogue 
$\tau$. Thus such a scenario does not have a major effect on our stability analysis.
At large absolute values of 
${Pe}$ and 
$\tau$, rapid changes in the stagnant lid's thickness and the interior temperature may lead to the non-dimensional parameters 
${Pe}$, 
$\tau$, 
$\zeta$ and 
$\tilde {Q}$ changing on relatively rapid time scales. In this event the approximation of a quasi-steady temperature profile in the base state will break down. However, the condition (4.21) for the onset of growth of the instability remains the same, provided 
$\varTheta _0$ and 
$\varTheta _l$ are adjusted to non-quasi-steady forms.
An interesting extension is suggested by the recent analysis of Vilella et al. (Reference Vilella, Limare, Jaupart, Farnetani, Fourel and Kaminski2018), who demonstrate a deviation from classical adiabatic temperature profiles in a convecting isoviscous mantle if the heating rate from radiogenic heating is comparable to or greater than the magnitude of the rate of temperature changes due to compression and decompression of a fluid parcel moving vertically. In our study the rate of change of temperature from radiogenic heating is 
$Q/c_p$, and from compression is
\begin{equation} w Tg\alpha/c_p. \end{equation}Here,
\begin{equation} w=\frac{F_H}{T_\mu\rho_0 c_p} \end{equation}is the convective vertical velocity scale. For our calculations the ratio between these rates
\begin{equation} \frac{QT_\mu\rho c_p}{Tg\alpha F_H}, \end{equation}
is below 0.03 at all times. In the cases with 
$Q_0\geq 2.0\times 10^{-11}\,\text {W}\,\text {kg}^{-1}$, namely those in which the instability is found, this ratio has a mean value below 
$3\times 10^{-3}$. In the case of a larger ratio, Vilella et al. (Reference Vilella, Limare, Jaupart, Farnetani, Fourel and Kaminski2018) suggest an increase in temperature at the upper part of the mixed interior, weakening or reversing the temperature gradient found in the adiabatic case. This suggests that rapid stagnant-lid growth would result in the base of the stagnant lid encountering cooler, more viscous material lower in the mixed interior, corresponding to lower convective heat flux and enhancing any instability.
While the pressure-dependent rheology this instability relies on is unlikely to be found on scales smaller than those of planetary mantle/magma ocean convection, it does have analogies in other areas of fluid mechanics, in particular in the solidification of binary alloys. Chemical gradients and their effect on the liquidus play a role similar to that of pressure in our analysis. In what is known as constitutional supercooling, a region adjacent to the solidification front may become locally cooler than the liquidus, due to thermal diffusion outpacing chemical diffusion. A protusion of solid material will extend into this region, resulting in a locally increased solidification rate in a positive feedback loop, with sinusoidal perturbations growing exponentially (Mullins & Sekerka Reference Mullins and Sekerka1964).
The instability covered in this paper relies on the planet or satellite being studied having a stagnant lid and being large enough to be spherical, and having a rocky composition. In terms of the present-day composition of the solar system, this restricts the instability to Earth's moon, Mars and Mercury. Smaller bodies such as Ceres have much smaller lithostatic pressure gradients, thus implying that they may not have a sufficiently strong depth dependence to any convective heat flux. Venus and Io are both very different from the scenario of a conductive stagnant lid, as their high rates of vulcanism result in heat release, in what is known as heat piping (O'Reilly & Davies Reference O'Reilly and Davies1981; Turcotte Reference Turcotte1989; Kankanamge & Moore Reference Kankanamge and Moore2016). In comparison, Earth's surface plate tectonics involve significant surface divergence and recycling, and are hence in a very different fluid dynamical regime from stagnant-lid scenarios.
6. Conclusions
Here we have shown that the pressure dependence of the solidus and liquidus of silicates leads to a strong pressure on the viscosity of a well-mixed partial melt, and hence on the heat flux from convecting planetary interiors. In cases where this convective heat flux changes more rapidly with depth than the conductive heat flux in the lithosphere, we have shown that a spherically symmetric stagnant lid will become unstable to asymmetric perturbations to its thickness. This instability provides a potential mechanism for small-amplitude initial variations in thickness to grow to a much larger scale. Smaller-wavenumber modes grow fastest at infinitesimal amplitude, leading to the promotion of degree-1 perturbations in particular. We suggest that this may lead to the development of hemispherical dichotomies seen in the crusts of Earth's moon and Mars, and may reconcile weak and short-lived forcing with the strong hemispherical differences seen in their lithospheres at present.
Funding
This work was funded by a Royal Society University Research Fellowships extension award to J.N. and C.M. has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No. 101001689).
Declaration of interests
The authors report no conflict of interest.
Open access
