cooling down at a bulk rate of approximately 100 K Gyr⁻¹ (Jaupart, Labrosse & Mareschal 2007). Extrapolating this budget back in time is fraught with many uncertainties. Although it is recognized that mantle convection is the main physical process involved in the cooling of the Earth, there is still considerable debate about the relationship between heat loss and internal temperature. The difficulty stems from the complexity of the Earth's mantle system, which loses heat through three different mechanisms: depletion of primordial heat that was accumulated during the accretion process, depletion of internal heat sources through radioactive decay and, last but not least, loss of radioactive elements to continental crust that does not participate in convective overturn. Although this is seldom discussed, the latter plays a major role. According to present estimates, the mantle may have lost as much as half of its radioactive elements to the continental crust (Jaupart et al. 2007). In such conditions, convection has become increasingly sensitive to heat supply from the core. The extraction of radioactive elements from the mantle is effected mostly in subduction zones where downgoing plates lose their radioactive elements to fluids and magmas that get incorporated in continental crust (Taylor & McLennan 1985). Downgoing plates therefore act as local sources of depleted material and are responsible for a heterogeneous distribution of mantle heat sources. Yet another consequence of continental growth is a change of upper boundary conditions. Whereas surface motions occur in the oceanic domain, a situation that can be approximated by a free surface, continents behave as rigid boundaries instead. One might therefore expect a gradual change of behaviour of the upper boundary through time, from that of a free surface in early stages to that of a rigid one later on. Our sister planet Venus, for example, currently does not allow surface motions (Turcotte 1995).

Dealing with the secular evolution of the Earth requires the handling of many different processes. Important goals include the rate of surface renewal and the rate of volcanism and degassing so that one can deduce the mass of the atmosphere from that of the solid planet. These require estimates of both heat loss, which sets the net rate of surface renewal, and internal temperature, which sets the depth range and rate of melting, as a function of time. Because of the variety of conditions that come into play, for want of an all-encompassing physical model that is presently beyond our reach, calculations are best tackled using a so-called 'parametrized' approach, such that the physics of convection is collapsed into a single equation relating the surface heat flux to the temperature difference across the upper boundary layer. This approach has been used in many natural systems, including lava flows (Garel et al. 2012), magma reservoirs (Jaupart & Brandeis 1986; Worster, Huppert & Sparks 1990), silicate planets (McKenzie & Weiss 1975; McNamara & van Keken 2000; Deschamps et al. 2012), the atmosphere (Krishnamurti 1997) and lakes (Neralla & Danard 1975), amongst others. Recently, it has also been used to study the evolution of super-Earths, i.e. massive planets in remote planetary systems (Kite, Manga & Gaidos 2009). The 'parametrized' approach relies on a scaling law, first introduced by Townsend (1964), which states that the surface heat flux $\phi$ is determined locally by the dynamics of the upper boundary layer independently of the total thickness of fluid

$$ \phi = C\lambda \left( \frac{\rho g\alpha}{\kappa\mu} \right)^{1/3} \Delta T_{TBL}^{4/3}, \quad (1.1) $$

where C is a proportionality constant, $\lambda$ is the thermal conductivity, $\rho$ is the density, $g$ is the acceleration of gravity, $\alpha$ is the thermal expansion coefficient, $\kappa$ is the thermal diffusivity, $\mu$ is the dynamic viscosity of the fluid and $\Delta T_{TBL}$ is the
temperature contrast across the upper thermal boundary layer (TBL). A large number of studies, both theoretical and experimental, have been devoted to delineate its domain of validity with a focus on the power-law exponent. However, the effects of the boundary conditions (Grasset & Parmentier 1998; Choblet & Parmentier 2009) and of the mode of heating of the fluid layer (Hansen, Yuen & Malevsky 1992; Sotin & Labrosse 1999; Parmentier & Sotin 2000) have received less attention, but they clearly deserve scrutiny. For example, for free surface boundary conditions, the constant of proportionality C changes by 25% between the two modes of heating (i.e. internal heating due to distributed heat sources or heating from below). According to the small number of studies available, the impact of the mechanical boundary conditions is more important and may be as large as 50% (Jaupart & Mareschal 2011). Such differences may seem small but they cannot be ignored, all the more so because the boundary conditions and dominant heating mode may change with time in a telluric planet, as explained above. Changing the value of C by 50% implies a 36% change of $\Delta T$, or $\approx 360$ K for $\Delta T \approx 10^3$ K. In comparison, the total amount of cooling of the Earth over the last 3 Gyr is only $\approx 200$ K (Herzberg, Condie & Korenaga 2010).

The characteristics of Rayleigh-Bénard convection at the high values of the Rayleigh number relevant to telluric planets have been investigated extensively both in the laboratory and with numerical calculations. In addition, a comprehensive theoretical framework leading to scaling laws for the main variables of interest over wide ranges of Rayleigh and Prandtl numbers is available (Grossmann & Lohse 2000). The vast majority of these studies, however, have dealt with rigid boundaries, which is not appropriate for the Earth. Except for the study by Katsaros et al. (1977), free surface boundary conditions have only been investigated using numerical calculations. As regards the characteristics of convection in fluid layers that are internally heated, our current knowledge is more fragmentary. Numerical calculations are only available for free surface boundary conditions in the limit of infinite Prandtl number (Parmentier, Sotin & Travis 1994; Grasset & Parmentier 1998; Sotin & Labrosse 1999; Parmentier & Sotin 2000) and there have been very few laboratory experiments owing to difficulties in controlling the distribution of the volumetric heating rate in a large volume. Most of these experiments aimed at documenting the planform of convection (Tritton & Zarraga 1967; Schwiderski & Schwab 1971; Kulacki & Goldstein 1972; Tasaka et al. 2005; Takahashi et al. 2010) and only one of them involved temperature and heat flux measurements (Kulacki & Nagle 1975). That latter study was conducted in water, a fluid with a Prandtl number of order 1, with a very small number of temperature probes. The experimental set-up allowed the determination of a scaling law for the heat flux but was ill-suited to the measurement of the vertical temperature profile through the layer, which is of interest in itself. The dearth of studies on internally heated fluids has motivated us to develop a novel experimental technique that achieves a uniform rate of volumetric heating in a fluid layer. The technique can also be used in other configurations, including a heterogeneous distribution of heat sources in controlled conditions. In order to validate the technique, the present work has been limited to a homogeneous distribution.

Our novel experimental design relies on the absorption of microwave (MW) radiation in a fluid layer with rigid boundaries. We are able to achieve high values of the Rayleigh-Roberts number, which is the analogue of the Rayleigh number for an internally heated fluid, at high Prandtl numbers. Non-invasive techniques are used to measure both temperature and velocity fields within the fluid layer and to
determine the horizontal planform of convection. We determine vertical profiles of the horizontally averaged temperature and convective heat flux and we derive scaling laws for the temperature difference and for the thickness of the upper thermal boundary layer. In order to assess the accuracy of the laboratory measurements, we carry out high-resolution numerical calculations in exactly the same three-dimensional (3D) configuration using the convection code Stag3D (Tackley 1993). These calculations account for the exact temperature dependence of the working fluid physical properties. For comparison with previous studies, we also carry out numerical simulations for free boundaries with constant fluid properties. Scaling laws for the main variables of interest are derived for the two types of boundary conditions. We show that, for a given bulk heat flux, the temperature difference across the upper boundary layer is significantly higher below a rigid boundary than below a free one. The importance of this is illustrated by comparing predictions for planets with free boundaries, like the Earth, and planets with rigid ones, like Venus.

## 2. Laboratory experiments and numerical simulations

### 2.1. Relevant dimensionless numbers

Convection generated by bottom heating and top cooling (Rayleigh-Bénard convection) is described by two dimensionless numbers: the Rayleigh number and the Prandtl number. The Rayleigh number, Ra, defines the vigour of convection and represents the ratio of the driving thermal buoyancy forces over the thermal and viscous dissipation,

$$Ra = \frac{\rho g \alpha \Delta T h^3}{\kappa \mu}, \qquad (2.1)$$

where $\Delta T$ is the temperature difference between the top and bottom of the layer and $h$ is the layer thickness. Convection starts when Ra exceeds a critical value (Chandrasekhar 1961), and follows a sequence of transitions toward chaos as Ra increases. The second parameter, the Prandtl number, represents the ratio of momentum diffusivity over heat diffusivity,

$$Pr = \frac{\nu}{\kappa}, \qquad (2.2)$$

where $\nu = \mu/\rho$ is the kinematic viscosity. When $Pr \gg 1$ inertial effects are negligible compared to viscous ones and the fluid motion stops as soon as the heat source is cut off. This is the case for telluric mantles, where $Pr > 10^{23}$.

In the purely internally heated case, the temperature scale for convection is related to the internal heating rate,

$$\Delta T_H = \frac{H h^2}{\lambda}, \qquad (2.3)$$

where $H$ is the heat generated per unit volume. The resulting Rayleigh-Roberts number (Roberts 1967) is

$$Ra_H = \frac{\rho g \alpha H h^5}{\lambda \kappa \mu}. \qquad (2.4)$$

Predicting the behaviour of complex natural systems such as planetary mantles requires the determination of scaling laws derived from fundamental physical principles that have been tested against experimental or numerical results. These scaling laws are applicable to natural systems only if the dynamic similarity is respected, i.e. the same boundary conditions (mechanical, thermal, aspect ratio), same rheology and similar balances between the various physical effects described by the dimensionless numbers.
<figure><img src="image_88.png" /><figcaption>FIGURE 1. (Colour online) Scheme of the microwave-based convection experiment.</figcaption></figure>

## 2.2. Laboratory experiments in the microwave oven

Figure 1 shows the structure of the MW oven, which includes a power generator driven by an embedded control system and an antenna through which the MW radiation propagates towards the working fluid, where it is uniformly absorbed. Our innovative design of the MW circuits guiding the MW radiation into the fluid ensures that a uniform MW field distribution is continuously maintained throughout the heating process (Surducan et al. 2014). The volumetrically heated fluid is cooled from above with an aluminium heat exchanger, which is temperature-controlled by a thermostatic bath. The tank (30 × 30 × 5 cm³) is made of poly(methyl methacrylate), so the bottom and side boundaries are as close as possible to adiabatic. The experimental mechanical boundary conditions are rigid.

The experimental fluids are transparent hydroxyethylcellulose-water mixtures whose viscosity can be varied within a wide range, depending on the polymer concentration. The fluid rheology was characterized with a Thermo Scientific Haake rheometer RS600, its density and thermal expansion were measured with a DMA 5000 Anton Paar densimeter (Limare et al. 2013) and its thermal diffusivity was measured by the photopyroelectric method (Dadarlat & Neamtu 2009). Values of $\Delta T_H$ and $Ra_H$ as well as fluid properties for the 30 experiments are listed in table 1. Experimental values of $Pr$ (3 × 10² < $Pr$ < 3 × 10⁴) are large enough for viscous effects to dominate over inertial ones (Davaille & Limare 2007). The experimental $Ra_H$ is between 5 × 10⁴ and 2 × 10⁷. The fluid properties are temperature-dependent: over the range of experimental conditions, the strongest viscosity reduction is 0.15 with respect to the value at surface temperature $T_0$, thermal expansion increases by a factor of maximum 5.2, whereas the temperature variation of thermal diffusivity is negligible.

During the experiments, the convecting fluid was scanned with a laser sheet over half of the tank size. The scattered light is registered by a charge-coupled device (CCD) E-lite camera (1.4 Mpixels, 17 Hz) from LaVision, allowing the measurement of both temperature and velocity fields without perturbing the flow. More specifically, to measure the temperature field we seeded the fluid with six types of thermochromic liquid crystal (TLC) evenly distributed between 19.3 and 36.6 °C. Each TLC produces one bright contour (i.e. an isotherm) when illuminated by a monochromatic light (Davaille et al. 2011). We calculated the temperature field of each two-dimensional (2D) cross-section by interpolating the isotherms. The surface temperature $T_0$ was
<table><caption>TABLE 1. Laboratory experimental parameters: ΔT<sub>H</sub> is calculated with H determined from the heat flux at the top surface at steady state; Ra<sub>H</sub> is calculated with fluid parameters at mean temperature T<sub>mean</sub> at steady state; μ<sub>0</sub> and α<sub>0</sub> are the viscosity and the thermal expansion at surface temperature T<sub>0</sub>. We also indicate the ratios of viscosity and thermal expansion at T<sub>0</sub> with respect to their values at T<sub>max</sub> to quantify the departure from the Boussinesq approximation.</caption><thead><tr><th>ΔT<sub>H</sub> (K)</th><th>Ra<sub>H</sub></th><th>ΔT<sub>TBL</sub>/ΔT<sub>H</sub></th><th>δ<sub>TBL</sub>/h</th><th>μ<sub>0</sub> (Pa s)</th><th>μ<sub>0</sub>/μ<sub>T<sub>max</sub></sub></th><th>α<sub>0</sub> (K<sup>-1</sup>)</th><th>α<sub>0</sub>/α<sub>T<sub>max</sub></sub></th></tr></thead><tbody><tr><td>94</td><td>5.01 × 10<sup>4</sup></td><td>0.1905</td><td>0.318</td><td>10.6</td><td>2.63</td><td>1.75 × 10<sup>-4</sup></td><td>0.46</td></tr><tr><td>142</td><td>8.22 × 10<sup>4</sup></td><td>0.1604</td><td>0.292</td><td>13.9</td><td>3.43</td><td>1.11 × 10<sup>-4</sup></td><td>0.30</td></tr><tr><td>148</td><td>9.17 × 10<sup>4</sup></td><td>0.1754</td><td>0.332</td><td>15.3</td><td>3.95</td><td>8.77 × 10<sup>-5</sup></td><td>0.23</td></tr><tr><td>72</td><td>1.03 × 10<sup>5</sup></td><td>0.1634</td><td>0.340</td><td>3.04</td><td>2.01</td><td>1.91 × 10<sup>-4</sup></td><td>0.58</td></tr><tr><td>114</td><td>1.48 × 10<sup>5</sup></td><td>0.1557</td><td>0.301</td><td>3.04</td><td>2.10</td><td>1.70 × 10<sup>-4</sup></td><td>0.46</td></tr><tr><td>81</td><td>2.99 × 10<sup>5</sup></td><td>0.1355</td><td>0.288</td><td>0.740</td><td>1.44</td><td>1.75 × 10<sup>-4</sup></td><td>0.58</td></tr><tr><td>154</td><td>3.05 × 10<sup>5</sup></td><td>0.1594</td><td>0.271</td><td>4.71</td><td>3.79</td><td>1.03 × 10<sup>-4</sup></td><td>0.27</td></tr><tr><td>167</td><td>3.47 × 10<sup>5</sup></td><td>0.1643</td><td>0.271</td><td>5.54</td><td>4.59</td><td>6.49 × 10<sup>-5</sup></td><td>0.17</td></tr><tr><td>75</td><td>4.02 × 10<sup>5</sup></td><td>0.1545</td><td>0.308</td><td>0.789</td><td>2.10</td><td>1.63 × 10<sup>-4</sup></td><td>0.55</td></tr><tr><td>83</td><td>4.21 × 10<sup>5</sup></td><td>0.1189</td><td>0.232</td><td>0.761</td><td>1.93</td><td>1.71 × 10<sup>-4</sup></td><td>0.60</td></tr><tr><td>94</td><td>6.58 × 10<sup>5</sup></td><td>0.1183</td><td>0.261</td><td>0.744</td><td>2.04</td><td>1.75 × 10<sup>-4</sup></td><td>0.52</td></tr><tr><td>94</td><td>6.69 × 10<sup>5</sup></td><td>0.1127</td><td>0.208</td><td>0.761</td><td>2.56</td><td>1.71 × 10<sup>-4</sup></td><td>0.59</td></tr><tr><td>144</td><td>1.07 × 10<sup>6</sup></td><td>0.1007</td><td>0.229</td><td>0.769</td><td>2.51</td><td>1.68 × 10<sup>-4</sup></td><td>0.51</td></tr><tr><td>136</td><td>1.11 × 10<sup>6</sup></td><td>0.1139</td><td>0.225</td><td>0.769</td><td>2.66</td><td>1.68 × 10<sup>-4</sup></td><td>0.49</td></tr><tr><td>131</td><td>1.32 × 10<sup>6</sup></td><td>0.1121</td><td>0.214</td><td>0.585</td><td>2.52</td><td>1.82 × 10<sup>-4</sup></td><td>0.53</td></tr><tr><td>57</td><td>1.44 × 10<sup>6</sup></td><td>0.1185</td><td>0.245</td><td>0.108</td><td>1.49</td><td>1.72 × 10<sup>-4</sup></td><td>0.69</td></tr><tr><td>248</td><td>2.47 × 10<sup>6</sup></td><td>0.0800</td><td>0.195</td><td>0.839</td><td>3.31</td><td>1.50 × 10<sup>-4</sup></td><td>0.41</td></tr><tr><td>153</td><td>2.50 × 10<sup>6</sup></td><td>0.0946</td><td>0.220</td><td>0.334</td><td>2.39</td><td>1.66 × 10<sup>-4</sup></td><td>0.51</td></tr><tr><td>145</td><td>3.05 × 10<sup>6</sup></td><td>0.0826</td><td>0.169</td><td>0.650</td><td>2.47</td><td>4.06 × 10<sup>-4</sup></td><td>0.75</td></tr><tr><td>170</td><td>3.13 × 10<sup>6</sup></td><td>0.0927</td><td>0.201</td><td>0.332</td><td>2.57</td><td>1.67 × 10<sup>-4</sup></td><td>0.49</td></tr><tr><td>94</td><td>3.14 × 10<sup>6</sup></td><td>0.0849</td><td>0.182</td><td>0.108</td><td>1.83</td><td>1.72 × 10<sup>-4</sup></td><td>0.66</td></tr><tr><td>100</td><td>3.34 × 10<sup>6</sup></td><td>0.0975</td><td>0.169</td><td>0.108</td><td>1.75</td><td>1.72 × 10<sup>-4</sup></td><td>0.64</td></tr><tr><td>131</td><td>4.01 × 10<sup>6</sup></td><td>0.1067</td><td>0.233</td><td>0.381</td><td>5.74</td><td>1.36 × 10<sup>-4</sup></td><td>0.47</td></tr><tr><td>108</td><td>4.04 × 10<sup>6</sup></td><td>0.0871</td><td>0.177</td><td>0.108</td><td>1.97</td><td>1.72 × 10<sup>-4</sup></td><td>0.62</td></tr><tr><td>141</td><td>4.88 × 10<sup>6</sup></td><td>0.0733</td><td>0.199</td><td>0.106</td><td>1.76</td><td>1.75 × 10<sup>-4</sup></td><td>0.61</td></tr><tr><td>150</td><td>5.62 × 10<sup>6</sup></td><td>0.0639</td><td>0.157</td><td>0.108</td><td>1.95</td><td>1.71 × 10<sup>-4</sup></td><td>0.59</td></tr><tr><td>126</td><td>6.04 × 10<sup>6</sup></td><td>0.0898</td><td>0.200</td><td>0.313</td><td>7.13</td><td>1.80 × 10<sup>-4</sup></td><td>0.59</td></tr><tr><td>138</td><td>7.18 × 10<sup>6</sup></td><td>0.0645</td><td>0.144</td><td>0.198</td><td>1.94</td><td>3.78 × 10<sup>-4</sup></td><td>0.80</td></tr><tr><td>249</td><td>1.04 × 10<sup>7</sup></td><td>0.0702</td><td>0.136</td><td>0.368</td><td>6.79</td><td>1.44 × 10<sup>-4</sup></td><td>0.44</td></tr><tr><td>256</td><td>1.47 × 10<sup>7</sup></td><td>0.0655</td><td>0.148</td><td>0.307</td><td>6.82</td><td>1.83 × 10<sup>-4</sup></td><td>0.50</td></tr></tbody></table>

chosen so that most of the 2D cross-sections contained several isotherms, thereby ensuring a good resolution of the thermal boundary layer. This is important because the temperature structure of the TBL is used to calculate the surface heat flux φ, and hence the actual volumetric heat source H, since at steady state φ = H × h, where h is the tank height.

To measure the velocity field, we seeded the fluid with small hollow glass spheres, behaving as passive tracers. Using particle image velocimetry (package DaVis from LaVision) we calculated the velocity field by cross-correlating successive images. We derived the 3D vertical velocity field by interpolation of a large number of vertical cross-sections obtained at different laser positions during an experiment. We then determined the horizontal planforms by contouring the distribution of vertical velocity in a horizontal plane. Figure 2 shows a vertical cross-section of the convecting
<figure><img src="image_89.png" /><figcaption>FIGURE 2. Vertical cross-section through the convecting fluid: isotherms (a) and the associated velocity field (b) for Ra$_H$ = 3 × 10⁵.</figcaption></figure>

fluid for an experiment at Ra$_H$ = 3.0 × 10⁵. The thermal structure revealed by the isotherms (figure 2a) clearly indicates that convection is characterized by several cold instabilities originating from the top boundary layer. The corresponding velocity field (figure 2b) shows that cold instabilities generate a localized downwelling flow, whereas broad zones of non-buoyant upwellings are associated with the return flow.

## 2.3. Numerical simulations

To validate the laboratory experiments, and to further extend the parameter space, we conduct numerical simulations in 3D Cartesian geometry using the parallel code Stag3D by Tackley (1993). This well-known code solves the equations governing convection in the limit of infinite Prandtl number and has been benchmarked successfully on several occasions: the mass conservation

$$ \nabla \cdot \mathbf{u} = 0, \qquad (2.5) $$

conservation of momentum

$$ 2\nabla \cdot \bar{\mu}\mathbf{e} - \nabla p = Ra_H \bar{\alpha}(T - T_0)\hat{z}, \qquad (2.6) $$

and conservation of energy

$$ \nabla \cdot \nabla T + H = \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T = \frac{\mathrm{D}T}{\mathrm{D}t}, \qquad (2.7) $$

where **u** is the velocity vector, **e** is the strain rate tensor, **p** is the dynamic pressure, **T** is the temperature and $\hat{z}$ is a unit vector in the vertical direction. Fluid properties are assumed constant except for the thermal expansivity and dynamic viscosity, which are normalized to their values at top temperature $T_0$ ($\bar{\alpha}$ and $\bar{\mu}$ are the dimensionless thermal expansivity and viscosity, respectively). The Cartesian domain is divided into 512 × 512 × 64 elements, with the same aspect ratio of 6 as in the experimental tank, corresponding to a horizontal and vertical resolution of 0.6 mm and 0.8 mm, respectively. The top boundary condition is isothermal whereas the bottom one is adiabatic. The vertical sides of the domain are reflecting.

Our first set of simulations is designed to reproduce the experimental conditions, in particular the rigid mechanical boundary conditions at top and bottom and the fluid properties. Thermal expansion and viscosity are temperature-dependent following the measurement made on laboratory fluids (table S1 in the supplementary material
<figure><img src="image_90.png" /><figcaption>FIGURE 3. Convection planforms obtained from the 2D horizontal cross-section of the vertical velocity: downwellings (black) and upwellings (white); (a–c) experimental results; (d–f) numerical simulations; (a,d) $Ra_H = 6 \times 10^4$, (b,e) $Ra_H = 1.5 \times 10^5$, (c,f) $Ra_H = 10^6$. Lateral dimensions are in millimetres. Gray shades are scaled to ± root mean squared velocity.</figcaption></figure>

available at http://dx.doi.org/10.1017/jfm.2015.347), while other fluid properties are constant. Figure 3 shows convection planforms that have been derived from the horizontal distribution of the vertical velocity component in both laboratory experiments (a–c) and numerical calculations (d–f). Three different patterns can be defined using either method. For $Ra_H < 10^5$ a steady spoke-like pattern is found, in agreement with previous studies (Ichikawa et al. 2006; Takahashi et al. 2010). For $Ra_H > 10^6$ the pattern is time-dependent and involves numerous isolated plume-like downwellings. In the intermediate Rayleigh–Roberts number range ($10^5 < Ra_H < 10^6$), the complex pattern is best described as a combination of residual spoke features and sheet-like downwellings.

For the sake of comparison with simple theoretical scalings, we have also carried out a second set of numerical simulations with rigid boundary conditions but with constant physical properties of the fluid. In a third set of calculations with constant fluid properties, we have investigated the effect of free-slip boundary conditions (see tables S2 and S3 in the supplementary material). Finally, we have changed the bottom boundary condition (no slip versus free slip) independently of the upper boundary condition and found no impact on the thermal structure and thickness of the upper boundary layer. This may be expected at large values of the Rayleigh–Roberts number because the dynamics of the upper boundary layer are locally determined. Figure S1 in the supplementary material illustrates this conclusion.

### 3. Scaling laws for the thermal boundary layer

Figure 4(a) shows the vertical profile of the horizontally averaged temperature at steady state. The thermal structure of the convecting layer can be split into an upper boundary layer and a convective interior; there is no basal thermal boundary layer, as no heat is supplied from below. An important feature is that the fluid interior has a slightly negative temperature gradient.
<figure><img src="image_91.png" /><figcaption>FIGURE 4. (Colour online) Numerical simulation using constant fluid parameters and rigid boundary conditions at Ra$_H$ = 10$^6$. (a) Vertical profile of the horizontally averaged temperature. (b) Vertical profiles of the three components of the heat balance (3.3): conductive heat flux (blue), cumulative heat generated (red) and convective heat flux (green). The thickness of the upper thermal boundary layer δ$_{TBL}$ is determined as the vertical position corresponding to the maximum convective heat flux. The temperature difference across the upper thermal boundary layer ΔT$_{TBL}$ corresponds to the temperature at a depth z = δ$_{TBL}$.</figcaption></figure>

If the local temperature is decomposed into its horizontal average $\bar{T}(z, t)$ and a fluctuation $\theta$, then the horizontally averaged heat equation can be written as

$$ \rho C_p \left[ \frac{\partial \bar{T}}{\partial t} + \frac{\partial \overline{w\theta}}{\partial z} \right] = \lambda \frac{\partial^2 \bar{T}}{\partial z^2} + H, \quad (3.1) $$

where $w$ is the vertical velocity and $C_p$ is the heat capacity. At steady state, the horizontally averaged temperature does not depend on time and (3.1) can be reduced to

$$ 0 = -\frac{d}{dz} \left( -\lambda \frac{d\bar{T}}{dz} + \rho C_p \overline{w\theta} \right) + H. \quad (3.2) $$

By integrating (3.2) between $z=h$ (the base of the fluid layer) and $z=0$, we obtain the convective heat flux at any depth $z$

$$ \rho C_p \overline{w\theta} = \lambda \frac{d\bar{T}}{dz} + H(h-z). \quad (3.3) $$

This equation involves three components, the convective and conductive heat fluxes and also the cumulative heat generation, which are illustrated in figure 4(b).

We derive scalings for the characteristics of the upper boundary layer using the viscous dissipation–buoyancy flux balance. For the large values of the Prandtl number of relevance here, the kinetic boundary layer extends over the whole fluid layer. The integral buoyancy flux for the layer is

$$ \int_0^h \rho \alpha g \overline{w\theta} \, dz = \frac{\alpha g}{C_p} \int_0^h \left[ \lambda \frac{d\bar{T}}{dz} + H(h-z) \right] \, dz = \frac{\alpha g}{C_p} \left( -\lambda \Delta T + H \frac{h^2}{2} \right). \quad (3.4) $$
For large values of the Rayleigh number, we may neglect $\lambda\Delta T$ compared to $Hh^2/2$. To derive the kinetic dissipation scale, we use velocity scale $U$ and assume that the relevant length scale for the flow is the fluid layer depth $h$ (Grossmann & Lohse 2000). The kinetic dissipation equation can now be written as

$$ \int_{0}^{h} \rho \alpha g \overline{w \theta} dz \sim \mu \left(\frac{U}{h}\right)^2 h \sim \frac{\alpha g}{C_p} H h^2, \quad (3.5) $$

where dissipation balances buoyancy. At steady state the surface heat flux $\phi = Hh$ and is such that

$$ \phi \sim \lambda \frac{\Delta T_{TBL}}{\delta_{TBL}}, \quad (3.6) $$

where $\Delta T_{TBL}$ and $\delta_{TBL}$ are the temperature difference and thickness of the upper thermal boundary layer, respectively (see figure 4). For a closure equation, we use the balance between horizontal advection and vertical diffusion in the boundary layer and obtain

$$ \delta_{TBL} \sim \left( \frac{\kappa h}{U} \right)^{1/2}. \quad (3.7) $$

Substituting (3.7) in (3.6) we obtain

$$ Hh \sim \lambda \Delta T_{TBL} \left( \frac{U}{\kappa h} \right)^{1/2} \quad (3.8) $$

and

$$ \frac{Hh^2}{\lambda} \sim \Delta T_{TBL} \left(\frac{U}{h}\right)^{1/2} \frac{h}{\kappa^{1/2}}. \quad (3.9) $$

Using the definition of the temperature scale in (2.3), the dissipation equation (3.5) and $\kappa = \lambda/C_p$, we finally obtain

$$ \frac{\Delta T_{TBL}}{\Delta T_H} \sim \left(\frac{h}{U}\right)^{1/2} \frac{\kappa^{1/2}}{h} \sim \left(\frac{\mu C_p}{\alpha g H h}\right)^{1/4} \frac{\kappa^{1/2}}{h} \sim \left(\frac{\mu \lambda \kappa}{\alpha g H h^5}\right)^{1/4}. \quad (3.10) $$

We therefore obtain the explicit dependence of the dimensionless temperature contrast $\Delta T_{TBL}/\Delta T_H$ as a power function of $Ra_H$,

$$ \frac{\Delta T_{TBL}}{\Delta T_H} = C_T Ra_H^\beta, \quad (3.11) $$

where $C_T$ is a constant scaling coefficient and $\beta = -1/4$.

One can determine a scaling law for the boundary layer thickness in two different ways. The first one is to use the heat flux and temperature drop (3.6), which leads to the same power-law scaling as the temperature drop. Alternatively, one can rely on the characteristics of the vertical convective heat flux profile and define the base of the thermal boundary layer where the convective heat flux (green curve) reaches its maximum value (Davaille & Jaupart 1993). This second method has the advantage of allowing us to check that the experimentally determined heat flux profile changes as a function of the Rayleigh number in a self-consistent manner. This leads to the following scaling law for the TBL thickness:

$$ \frac{\delta_{TBL}}{h} = C_{\delta} Ra_{H}^{\beta}, \quad (3.12) $$

where $C_\delta$ is a constant scaling coefficient.
<figure><img src="image_92.png" /><figcaption>FIGURE 5. (a,b) Dimensionless TBL temperature (a) and thickness (b) versus Ra<sub>H</sub> for experiments (full circles) and numerical simulations (empty squares) obtained with rigid boundary conditions and temperature dependence of the fluid parameters. (c,d) Numerical simulations results for the dimensionless TBL temperature (c) and thickness (d) versus Ra<sub>H</sub>. Black lines are power-law best fits with $\beta = -1/4$. Legends indicate the implemented conditions.</figcaption></figure>

Figure 5(a,b) shows the dimensionless TBL temperature and thickness for laboratory experiments (full circles) and corresponding numerical simulations (empty squares) also provided in table 1 and table S1–S3 in the supplementary material. The experimental and numerical results are in good agreement over a large range of Rayleigh number. The larger scatter of the laboratory determinations can be attributed to errors inherent in the interpolation of discrete isotherms to obtain the temperature field. The power-law scalings that have been derived above capture the main trends of the data for both the temperature contrast and the boundary layer thickness. As discussed below, these scalings rely on the Boussinesq approximation, which is not entirely valid for these experiments. For $\beta = -1/4$ we can estimate the best-fitting scaling coefficients: $C_T = 3.56$, $C_\delta = 7.36$ for the experiments, and similar values ($C_T = 3.58$, $C_\delta = 7.50$) are found for the numerical simulations (table 2). The linear fit was calculated using standard linear regression analysis performed in Matlab with input data log(RaH) and log(ΔT<sub>TBL</sub>/ΔT<sub>H</sub>) and log(δ<sub>TBL</sub>/h), respectively.

The next issue we want to address is to what extent these coefficients are affected by the mechanical boundary conditions and by the departure from the Boussinesq approximation. Figure 5(c,d) shows our three sets of numerical simulations: with rigid boundary conditions and temperature-dependent fluid properties (empty squares),
<table><caption>TABLE 2. Scaling law constants for the thermal boundary layer when $\beta = -1/4$.</caption><thead><tr><th>Conditions</th><th>C<sub>T</sub></th><th>C<sub>δ</sub></th></tr></thead><tbody><tr><td>Numerical: free slip, α = α<sub>0</sub>, μ = μ<sub>0</sub></td><td>2.49</td><td>5.90</td></tr><tr><td>Numerical: rigid, α = α<sub>0</sub>, μ = μ<sub>0</sub></td><td>3.41</td><td>7.08</td></tr><tr><td>Numerical: rigid, α = f(T), μ = f(T)</td><td>3.58</td><td>7.50</td></tr><tr><td>Experimental: rigid, α = f(T), μ = f(T)</td><td>3.56</td><td>7.36</td></tr></tbody></table>

<table><caption>TABLE 3. Scaling law constants for the thermal boundary layer when the $\beta$ value is left to vary.</caption><thead><tr><th>Conditions</th><th>β<sub>T</sub></th><th>C<sub>T</sub></th><th>β<sub>δ</sub></th><th>C<sub>δ</sub></th></tr></thead><tbody><tr><td>Numerical: free slip, α = α<sub>0</sub>, μ = μ<sub>0</sub></td><td>-0.246 ± 0.012</td><td>2.35</td><td>-0.232 ± 0.011</td><td>4.51</td></tr><tr><td>Numerical: rigid, α = α<sub>0</sub>, μ = μ<sub>0</sub></td><td>-0.241 ± 0.013</td><td>2.96</td><td>-0.235 ± 0.009</td><td>5.63</td></tr><tr><td>Numerical: rigid, α = f(T), μ = f(T)</td><td>-0.215 ± 0.041</td><td>2.21</td><td>-0.195 ± 0.037</td><td>3.50</td></tr><tr><td>Experimental: rigid, α = f(T), μ = f(T)</td><td>-0.222 ± 0.042</td><td>2.47</td><td>-0.167 ± 0.040</td><td>2.32</td></tr></tbody></table>

with rigid boundary conditions and constant fluid parameters (full squares), and with free-slip boundary conditions and constant fluid parameters (full triangles). For rigid boundary conditions, the use of constant fluid parameters significantly reduces the scatter, and the best-fitting coefficients are $C_T = 3.41$, $C_\delta = 7.08$ (table 2). The most striking result is the difference between rigid and free-slip boundary conditions: at all values of $Ra_H$, the ratio $\Delta T_{TBL}/\Delta T_H$ is always lower for free-slip than for rigid boundary conditions, so that the best-fitting coefficient for free slip is $C_T = 2.49$, a value in agreement with Parmentier & Sotin (2000). In other words, $C_T$ for free slip is 37% lower than for rigid boundary conditions, and $C_\delta = 5.90$ is 22% lower. This result has important implications for the interior temperature of planets as discussed in the following section. Here we note that the scaling laws (3.11) and (3.12) are valid from particularly low values of Rayleigh-Roberts number ($Ra_H \sim 6 \times 10^4$), for which convection is time-independent, to the highest values explored ($Ra_H = 10^9$). Values for the power-law coefficients have also been obtained with $\beta$ left to vary (table 3). We focus on $\Delta T_{TBL}$ because it is determined more precisely than $\delta_{TBL}$ (the peak value can be determined more precisely than the location of the peak). In the example of figure 4(a), the spatial resolution of our set-up (0.6 mm) leads to an uncertainty of approximately 12% for the boundary layer thickness, whereas the measurement error on the temperature difference is only 0.3%. We find that, within the range of uncertainty, values for exponent $\beta$ are consistent with the result of our simple analysis (3.11). As expected, for constant fluid properties, the exponent $\beta$ is closer to -1/4 than for the experiments and numerical simulations with temperature-dependent fluid properties.

## 4. Planetary application: mantle potential temperature for the Earth and Venus

Our scaling laws can be used to estimate mantle potential temperature in telluric planets dominated by internally heated convection. We are aware that our analysis ignores the depth and temperature dependence of mantle properties, yet it is worthwhile to investigate the impact of the boundary conditions on the internal thermal structure of the Earth and Venus, and this can be done in a first step with a simple
<table><caption>TABLE 4. Parameters for the Earth and Venus used in our estimation.</caption><thead><tr><th>Parameters</th><th>Earth</th><th>Venus</th></tr></thead><tbody><tr><td>Thermal expansion coefficient (K<sup>-1</sup>)</td><td>3 &times; 10<sup>-5</sup></td><td>3 &times; 10<sup>-5</sup></td></tr><tr><td>Thermal conductivity (W K<sup>-1</sup> m<sup>-1</sup>)</td><td>3</td><td>3</td></tr><tr><td>Thermal diffusivity (m<sup>2</sup> s<sup>-1</sup>)</td><td>10<sup>-6</sup></td><td>10<sup>-6</sup></td></tr><tr><td>Mean mantle density (kg m<sup>-3</sup>)</td><td>4000</td><td>3500</td></tr><tr><td>Viscosity (Pa s)</td><td>10<sup>21</sup></td><td>10<sup>21</sup></td></tr><tr><td>Gravitational acceleration (m s<sup>-2</sup>)</td><td>9.81</td><td>8.87</td></tr><tr><td>Surface temperature (&deg;C)</td><td>0</td><td>462</td></tr><tr><td>Mantle depth (km)</td><td>2890</td><td>2932</td></tr><tr><td>Bulk silicate radioactive heating (TW)</td><td>19</td><td>16</td></tr><tr><td>Crust radioactive heating (TW)</td><td>8</td><td>(?)</td></tr><tr><td>Mantle radioactive heating (TW)</td><td>9–17</td><td>(?)</td></tr><tr><td>Present-day mantle cooling/heating (TW)</td><td>1–29</td><td>(?)</td></tr><tr><td>Heat flux from the core (TW)</td><td>5–17</td><td>3–14</td></tr><tr><td>Total mantle heat Q (TW)</td><td>38 &plusmn; 2</td><td>3–30</td></tr></tbody></table>

physical system. These two telluric planets share many similarities, but the mechanical top boundary conditions differ, since on Venus there is no active plate tectonics. On both planets, the bottom boundary condition is set by the liquid outer core and hence is closer to free-slip behaviour than to a rigid one. However, as explained in § 2.3, the characteristics of the upper boundary layer are not sensitive to the behaviour of the lower boundary. The same is true for the internal thermal structure, except for a thin lower region, which is not relevant to the present discussion. Additional differences between our experiments and planetary mantles are the presence of sidewalls and geometrical characteristics (i.e. a horizontal layer and a spherical shell). In a large-aspect-ratio tank such as ours, however, the bulk convection characteristics are affected only weakly by local sidewall motions. This is shown by comparing our results with those of Parmentier & Sotin (2000), for example, who carried out high-precision numerical calculations in domains with different values of the aspect ratio. We also expect that the spherical shell geometry does not influence the upper boundary layer characteristics greatly (see Jarvis, Glatzmaier and Vangelov 1995).

For the Earth, the total mantle heat source available for convection is known: $Q = 38 \pm 2$ TW (Jaupart et al. 2007). Using the parameters given in table 4, we obtain $Ra_H = 3 \times 10^9$ and $\Delta T_H = 10^5$ K. We calculate $\Delta T_{TBL}$ according to (3.11) for the range of $C_T$ previously determined, since low $C_T$ values (free-slip condition) are more appropriate to moving oceanic lithosphere whereas high $C_T$ values (rigid condition) better apply to continents. Figure 6(a) illustrates the mantle potential temperature versus the coefficient $C_T$, showing an increasing trend from free slip to rigid. Grey rectangles represent a compilation of mantle potential temperature values (1280–1450°C) from McKenzie & Bickle (1988), Jaupart et al. (2007) and Putirka et al. (2007) and temperatures (1450–1550°C) at the base of the continents (Jaupart & Mareschal 2011). These temperature ranges superpose well on our estimates, showing that a simple temperature scaling law is able to describe the present thermal state of the Earth, as inferred by geological models. Our scaling laws are based on a homogeneous distribution of internal heat sources, and it is not clear whether or not this is also true for planetary mantles. It has been proposed that the Earth's mantle is stratified, with a lower region that is enriched compared to the bulk, but current
<figure><img src="image_93.png" /></figure>

FIGURE 6. (Colour online) (a) Estimates of the Earth's potential temperature using the temperature scaling law, as a function of the coefficient $C_T$. The solid line was obtained for a value of total mantle heat source $Q=38$ TW, and dashed lines for $Q=38 \pm 2$ TW. Rectangles represent estimates of the potential oceanic and continental upper mantle temperatures (McKenzie & Bickle 1988; Jaupart et al. 2007; Putirka, Perfit, Ryerson & Jackson 2007; Jaupart & Mareschal 2011). (b) Estimates of Venus's potential temperature using the temperature scaling law for rigid (solid line) and free-slip (dashed line) boundary conditions, as a function of the total mantle heat source. Red curves represent the solidus potential temperature at a pressure corresponding to the thermal boundary layer thickness: rigid (solid line) and free slip (dashed line).

constraints on the uranium and thorium contents of the Earth are not sufficiently precise for a definitive conclusion (see the discussion of mantle heat production in Jaupart et al. (2007)). Radioactive decay is not the only energy source driving convective motions in the interior of planets, and its contribution is probably not as large as the loss of sensible heat (Jaupart et al. 2007). In all cases, one may reasonably expect that the dynamics of the upper boundary layer are not sensitive to the details of the vertical (or radial) distribution of heat sources. This deserves
an independent study with a degree of stratification that is left to vary because, as mentioned above, current constraints on the Earth's mantle are not robust. In fact, as stated in the introduction, one motivation for our experimental design was to allow us to investigate the consequences of a heterogeneous distribution of heat sources.

Surface boundary conditions on Venus may alternate between rigid, as present day, and free-slip during short periods of resurfacing (Turdotte 1995; Armann & Tackley 2012). Contrary to the Earth, the total mantle heat source is poorly constrained (3–30 TW; see table 4). The bulk silicate radioactive heat source, scaled to the Earth’s mass, is 16 TW (Smrekar & Sotin 2012). The crustal radioactive heating is unknown, although the five measurements by Vega and Venera landers indicate concentrations of heat-producing elements similar to mid-oceanic ridge basalts and oceanic islands basalts (Turdotte 1995). Therefore, Venus’s present-day crust might contain a non-negligible quantity of radioactive elements that do not contribute to mantle convection. Heat flux from the core, estimated by Nimmo (2002) to be between 1.2 and 3.7 TW, could be as high as 3 and 14 TW if recalculated with recent thermal conductivity values (Gomi et al. 2013). Given all these uncertainties, we calculate Venus’s potential temperature over a range of total mantle heat sources to be between 3 and 30 TW (figure 6b).

For rigid boundary condition ($C_T = 3.4$, black, solid line), the potential temperature is higher than for free slip ($C_T = 2.5$, black, dashed line). High potential temperatures may lead to widely spread zones of partial melting beneath the thermal boundary layer. To check this hypothesis, we calculate the dry peridotite solidus temperature (Smrekar & Sotin 2012) at the pressure corresponding to the base of the TBL. Equation (3.12) gives a TBL depth between 75 and 140 km for free slip and 90–170 km for rigid boundary conditions. Solidus temperatures are then converted to potential temperatures (red lines in figure 6b). We find that in the free-slip case, partial melting would occur for mantle heat source values >19 TW, whereas for the rigid case, partial melting would start at values >14 TW. Another implication is that a change of the top boundary condition from free slip to rigid favours partial melting since the increase in potential temperature is larger than the increase in solidus temperature induced by the thermal boundary layer thickening.

## 5. Conclusion

We present a new experimental method to generate internal heating by microwave absorption. This prototype offers the ability to reach high $Ra_H$ and $Pr$ numbers, relevant for planetary convection. Our experimental results are compared with numerical simulations conducted in 3D Cartesian geometry, thereby providing the first cross-validation of experimental and numerical studies of convective viscous systems heated from within. We find that thermal boundary layer temperature and thickness scale with $Ra_H^{-1/4}$ as theoretically predicted by scaling arguments on the dissipation of kinetic energy (Jaupart & Mareschal 2011). We further test numerically the effect of the top mechanical boundary condition. At a given $Ra_H$, the boundary layer temperature is 37% higher for rigid than for free-slip boundary condition and its thickness is 22% greater.

Temperature scaling laws are used to evaluate the mantle temperatures of the Earth and Venus. For the Earth, our estimates match well the potential temperature of the oceanic upper mantle and the temperature at the base of continents inferred from geological data. For Venus, our results indicate the probable presence of partial melting underneath its lithosphere.
Finally, our specific microwave-based method offers the new perspective, unattained up to now experimentally, to selectively heat different zones of a convecting fluid, analogous to heterogeneous convection in the presence of chemical reservoirs with distinct concentration of radioactive isotopes.

## Acknowledgements

We would like to thank the reviewers for their fruitful comments that improved the manuscript and J.-M. Chomaz for his editorial handling. This work was funded by the ANR-11-IS04-0004 project for the French team and by the 1 RO-FR-22-2011 Romanian-French bilateral project for the Romanian team. Numerical computations were performed on the S-CAPAD platform, IPGP, France, and using HPC resources from GENCI-IDRIS (grant 2013-047033). This is IPGP contribution number 3640.

## Supplementary data

Supplementary data are available at http://dx.doi.org/10.1017/jfm.2015.347.

## REFERENCES

* ARMANN, M. & TACKLEY, P. J. 2012 Simulating the thermochemical magmatic and tectonic evolution of Venus's mantle and lithosphere: two-dimensional models. *J. Geophys. Res.* **117**, E12003.

* CHANDRASEKHAR, S. 1961 *Hydrodynamic and Hydromagnetic Stability*. Oxford University Press.

* CHOBLET, G. & PARMENTIER, E. M. 2009 Thermal convection heated both volumetrically and from below: implications for predictions of planetary evolution. *Phys. Earth Planet. Inter.* **173**, 290–296.

* DADARLAT, D. & NEAMTU, C. 2009 High accuracy photopyroelectric calorimetry of liquids. *Acta Chim. Slov.* **56**, 225–236.

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

* DAVAILLE, A. & LIMARE, A. 2007 Laboratory studies on mantle convection. In *Treatise of Geophysics* (ed. D. Bercovici & G. Schubert), pp. 89–165. Elsevier.

* DAVAILLE, A., LIMARE, A., TOUITOU, F., KUMAGAI, I. & VATTEVILLE, J. 2011 Anatomy of a laminar starting thermal plume at high Prandtl number. *Exp. Fluids* **50** (2), 285–300.

* DESCHAMPS, F., YAO, C., TACKLEY, P. J. & SANCHEZ-VALLE, C. 2012 High Rayleigh number thermal convection in volumetrically heated spherical shells. *J. Geophys. Res.* **117**, E09006.

* GAREL, F., KAMINSKI, E., TAIT, S. & LIMARE, A. 2012 An experimental study of the surface thermal signature of hot subaerial isoviscous gravity currents: implications for thermal monitoring of lava flows and domes. *J. Geophys. Res.* **117**, B02205.

* GOMI, H., OHTA, K., HIROSE, K., LABROSSE, S., CARACAS, R., VERSTRAETE, M. J. & HERNLUND, J. W. 2013 The high conductivity of iron and thermal evolution of the Earths core. *Phys. Earth Planet. Inter.* **224**, 88–103.

* GRASSET, O. & PARMENTIER, E. M. 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**, 171–181.

* GROSSMANN, S. & LOHSE, D. 2000 Scaling in thermal convection: a unifying theory. *J. Fluid Mech.* **407**, 27–56.

* HANSEN, U., YUEN, D. A. & MALEVSKY, A. V. 1992 Comparison of steady-state and strongly chaotic thermal convection at high Rayleigh number. *Phys. Rev. A* **46**, 4742–4754.

* HERZBERG, C., CONDIE, K. & KORENAGA, J. 2010 Thermal history of the Earth and its petrological expression. *Earth Planet. Sci. Lett.* **292**, 79–88.

* ICHIKAWA, H., KURITA, K., YAMAGISHI, Y. & YANAGISAWA, T. 2006 Cell pattern of thermal convection induced by internal heating. *Phys. Fluids* **18**, 038101,1–4.
JARVIS, G. T., GLATZMAIERAND, G. A. & VANGELOV, V. I. 1995 Effects of curvature, aspect ratio and plan form in two- and three-dimensional spherical models of thermal convection. *Geophys. Astrophys. Fluid Dyn.* **79** (1-4), 147-171.

JAUPART, C. & BRANDEIS, G. 1986 The stagnant bottom layer of convecting magma chambers. *Earth Planet. Sci. Lett.* **80**, 183-199.

JAUPART, C., LABROSSE, S. & MARESCHAL, J.-C. 2007 Temperatures, heat and energy in the mantle of the Earth. In *Treatise of Geophysics* (ed. D. Bercovici & G. Schubert), pp. 253-303. Elsevier.

JAUPART, C. & MARESCHAL, J.-C. 2011 *Heat Generation and Transport in the Earth*. Cambridge University Press.

KATSAROS, K. B., LIU, W. T., BUSINGER, J. A. & TILLMAN, J. E. 1977 Heat thermal structure in the interfacial boundary layer measured in an open tank of water in turbulent free convection. *J. Fluid Mech.* **83**, 311-335.

KITE, E. S., MANGA, M. & GAIDOS, E. 2009 Geodynamics and rate of volcanism on massive Earth-like planets. *Astrophys. J.* **700**, 1732-1749.

KRISHNAMURTI, R. 1997 Convection induced by selective absorption of radiation: a laboratory model of conditional stability. *Dyn. Atmos. Oceans* **27**, 367-382.

KULACKI, F. A. & GOLDSTEIN, R. J. 1972 Thermal convection in a horizontal fluid layer with volumetric heat sources. *J. Fluid Mech.* **271**, 271-287.

KULACKI, F. A. & NAGLE, M. E. 1975 Natural convection in a horizontal fluid layer with volumetric energy sources. *J. Heat Transfer* **97**, 204-211.

LIMARE, A., SURDUCAN, E., SURDUCAN, V., NEAMTU, C., DI GIUSEPPE, E., VILELLA, K., FARNETANI, C. G., KAMINSKI, E. & JAUPART, C. 2013 Microwave-based laboratory experiments for internally-heated mantle convection. In *Processes in Isotopes and Molecules* (PIM 2013): AIP Conference Proceedings, vol. 1565, pp. 14-16. AIP.

MCKENZIE, D. & BICKLE, M. 1988 The volume and composition of melt generated by extension of lithosphere. *J. Petrol.* **29**, 625-679.

MCKENZIE, D. & WEISS, N. O. 1975 Speculations on the thermal and tectonic history of the Earth. *J. R. Astron. Soc.* **42**, 131-174.

MCNAMARA, A. K. & VAN KEKEN, P. E. 2000 Cooling of the earth: a parameterized convection study of whole versus layered models. *Geochem. Geophys. Geosyst.* **1**, 2000GC000045.

NERALLA, V. R. & DANARD, M. B. 1975 Incorporation of parameterized convection in the synoptic study of large scale effects of the great lakes. *Mon. Weath. Rev.* **103**, 388-405.

NIMMO, F. 2002 Why does Venus lack a magnetic field? *Geology* **30**, 387-990.

PARMENTIER, E. M. & SOTIN, C. 2000 Three-dimensional numerical experiments on thermal convection in a very viscous fluid: implications for the dynamics of the thermal boundary layer at high Rayleigh number. *Phys. Fluids* **12**, 609-617.

PARMENTIER, E. M., SOTIN, C. & TRAVIS, B. J. 1994 Turbulent 3-D thermal convection in an infinite Prandtl number, volumetrically heated fluid: implication for mantle dynamics. *Geophys. J. Intl* **116**, 241-251.

PUTIRKA, K. D., PERFIT, M., RYERSON, F. J. & JACKSON, M. G. 2007 Ambient and excess mantle temperatures, olivine thermometry, and active vs passive upwelling. *Chem. Geol.* **241**, 177-206.

ROBERTS, P. H. 1967 Convection in horizontal layers with internal heat generation. Theory. *J. Fluid Mech.* **30**, 33-49.

SCHWIDERSKI, E. W. & SCHWAB, H. J. A. 1971 Convection experiments with electrolytically heated fluid layers. *J. Fluid Mech.* **48**, 703-717.

SMREKAR, S. E. & SOTIN, C. 2012 Constraints on mantle plumes on Venus: implications for volatile history. *Icarus* **217**, 510-523.

SOTIN, C. & LABROSSE, S. 1999 Three-dimensional thermal convection in an iso-viscous, infinite Prandtl number fluid heated from within and below: applications to the transfer of heat through planetary mantles. *Phys. Earth Planet. Inter.* **112**, 171-190.

SURDUCAN, E., SURDUCAN, V., LIMARE, A., NEAMTU, C. & DI GIUSEPPE, E. 2014 Microwave heating device performing internal-heating convection experiments, applied to Earth's mantle dynamics. *Rev. Sci. Instrum.* **85**, 124702.
TACKLEY, P. J. 1993 Effects of strongly temperature-dependent viscosity on time-dependent, 3-dimensional models of mantle convection. *Geophys. Res. Lett.* **20**, 2187–2190.

TAKAHASHI, J., TASAKA, Y., MURAI, Y., TAKEDA, Y. & YANAGISAWA, T. 2010 Experimental study of cell pattern formation induced by internal heat sources in a horizontal fluid layer. *Int. J. Heat Mass Transfer* **53**, 1483–1490.

TASAKA, Y., KUDOH, Y., TAKEDA, Y. & YANAGISAWA, T. 2005 Experimental investigation of natural convection induced by internal heat generation. *J. Phys.: Conf. Series* **14**, 168–179.

TAYLOR, S. R. & MCLENNAN, S. M. 1985 *The Continental Crust: Its Composition and Evolution*. Blackwell Scientific Publications.

TOWNSEND, A. A. 1964 Natural convection in water over an ice surface. *Q. J. R. Meteorol. Soc.* **90**, 248–259.

TRITTON, D. J. & ZARRAGA, M. M. 1967 Convection in horizontal fluid layers with internal heat generation experiments. *J. Fluid Mech.* **30**, 21–31.

TURCOTTE, D. L. 1995 How does Venus lose heat? *J. Geophys. Res.* **100**, 16931–16940.

WORSTER, M. G., HUPPERT, H. E. & SPARKS, R. S. J. 1990 Convection and crystallization in magma cooled from above. *Earth Planet. Sci. Lett.* **101** (1), 78–89.
# Fundamentals of laminar free convection in internally heated fluids at values of the Rayleigh–Roberts number up to $10^9$

Kenny Vilella¹,²,†, Angela Limare¹, Claude Jaupart¹, Cinzia G. Farnetani¹, Loic Fourel¹,³ and Edouard Kaminski¹

¹Institut de Physique du Globe, Univ. Paris Diderot, Sorbonne Paris Cité, CNRS, F-75005 Paris, France

²Institute of Earth Sciences, Academia Sinica, Taipei 11529, Taiwan

³Geological Survey of Norway, 7040 Trondheim, Norway

(Received 29 August 2017; revised 15 February 2018; accepted 10 April 2018)

Motions in the solid mantle of silicate planets are predominantly driven by internal heat sources and occur in laminar regimes that have not been systematically investigated. Using high-resolution numerical simulations conducted in three dimensions for a large range of Rayleigh–Roberts numbers ($5 \times 10^3 \le Ra_H \le 10^9$), we have determined the characteristics of flow in internally heated fluid layers with both rigid and free slip boundaries. Superficial planforms evolve with increasing $Ra_H$ from a steady-state tessellation of hexagonal cells with axial downwellings to time-dependent clusters of thin linear downwellings within large areas of nearly isothermal fluid. The transition between the two types of planforms occurs as a remarkable flow polarity reversal over a small $Ra_H$ range, such that downwellings go from isolated cylindrical structures encircled by upwellings to thin interconnected linear segments outlining polygonal cells. In time-dependent regimes at large values of $Ra_H$, linear downwellings dominate the flow field at shallow depth but split and merge at intermediate depths into nearly cylindrical plume-like structures that go through the whole layer. With increasing $Ra_H$, the number of plumes per unit area and their velocities increase whilst the amplitude of thermal anomalies decreases. Scaling laws for the main flow characteristics are derived for $Ra_H$ values in a $10^6–10^9$ range. For given $Ra_H$, plumes are significantly colder, narrower and wider apart beneath free boundaries than beneath rigid ones. From the perspective of planetary studies, these results alert to the dramatic changes of convective planform that can occur along secular cooling.

**Key words:** convection, pattern formation, plumes/thermals

## 1. Introduction

Geological activity is a consequence of the convective motions that develop in the Earth's mantle. In order to understand how it is generated, forward models of convection of increasing complexity have multiplied over the years (e.g. McKenzie,

†Email address for correspondence: vilella@earth.sinica.edu.tw
Roberts & Weiss 1973; Christensen & Hofmann 1994; Samuel & Farnetani 2003; Nakagawa & Tackley 2014). These models are, however, inherently questionable because of fundamental uncertainties on the sources and amount of energy involved. According to recent studies (e.g. Jaupart et al. 2015), approximately half of the Earth's heat loss today is due to secular cooling at a rate of approximately 100 K per billion years and the other half comes from the radioactive decay of uranium, thorium and potassium in silicate rocks of the crust and mantle. Large uncertainties arise when one attempts to focus on the mantle, since the amount of heat received from the core, induced by its secular cooling, remains poorly constrained (Labrosse, Poirier & Le Mouël 2001; Lay et al. 2006; Cottaar & Buffett 2012; Gomi et al. 2013). This has led to models that can be described as semi-empirical because they rely heavily on current observations such as the motions of surface plates and the distribution of internal density anomalies deduced from gravity and seismological studies (Mitrovica & Forte 2004; Deschamps & Tackley 2008, 2009). These models are focused on the current convection regime and incorporate most of the available constraints, implying that they are difficult to test and ill-suited for studies of planetary evolution back in time (Bello et al. 2014). An alternative strategy has been to try to infer the controls on convective motions from the observations themselves. For instance, Deschamps, Rogister & Tackley (2018) have proposed to constrain the presence of chemical heterogeneities at the base of the Earth's mantle with the topography of the core–mantle boundary, while Mitrovica & Forte (2004) have combined models of Earth's mantle convection with glacial isostatic adjustment data to infer the mantle viscosity profile. On Earth, however, one runs the risk of being mistaken because of the many complicating factors involved, which include the presence of continents or multiple mantle phase changes (e.g. Gurnis 1988; Guillou & Jaupart 1995; Tackley 1998b; Stixrude & Lithgow-Bertelloni 2011; Rolf, Coltice & Tackley 2012). Thus, whether all the available observations form a self-consistent set, such that they can all be accounted for by the workings of a single convective system, is an important issue. In silicate planets other than Earth, one has to work with fewer observations but, with a sound understanding of the characteristics of convection, one can still make useful inferences. For example, the main features of the gravity field and surface topography on Venus have been explained by a small number of mantle plumes, nine to be precise, over a uniform background (Smrekar & Sotin 2012). If this inference is valid, it leads to powerful constraints on the energy budget and interior viscosity of the planet.

In studies of silicate planets and planetary bodies, one cannot draw only from the comprehensive body of knowledge that is available for Rayleigh–Bénard convection because interior motions are largely driven by internal heat sources instead of heat supplied to the base. Heat sources include long-period radionuclides as well as contributions from secular cooling and tidal heating, although the latter is only dominant in some planetary bodies such as the icy satellites of Saturn and Jupiter (Hussmann et al. 2010; Schubert et al. 2010). In spite of their relevance to planetary evolution, studies of internally heated convective systems have lagged behind those for Rayleigh–Bénard set-ups (Goluskin 2015). Moreover, such systems are particularly interesting from theoretical and physical perspectives because they involve a single unstable thermal boundary layer that is not affected by motions generated at another boundary (Goluskin & Spiegel 2012; Seis 2013; Goluskin 2015). For instance, this may enable us to better capture the dynamics of the thermal boundary layer, especially the mechanisms generating downwellings and upwellings.

For a fluid layer with internal heating rate *H* above an adiabatic base, dimensional analysis of the mass, momentum and energy conservation equations in the Boussinesq
(1903) approximation leads to two dimensionless numbers, the Rayleigh–Roberts number $Ra_H$ and the Prandtl number $Pr$ (Roberts 1967; Schubert, Turcotte & Olson 2001):

$$Ra_H = \frac{\rho g \alpha H d^5}{\eta \kappa k}, \qquad (1.1)$$

$$Pr = \frac{\eta}{\rho \kappa}, \qquad (1.2)$$

where $\rho$ is the density, $g$ the acceleration of gravity, $\alpha$ the thermal expansion coefficient, $d$ the thickness of the fluid layer, $\eta$ viscosity, $\kappa$ thermal diffusivity and $k$ thermal conductivity. Note that all fluid properties are taken in the reference state. The influence of the Prandtl number has only been investigated in a cursory fashion, despite the very large gap that exists between values for water ($Pr \approx 7$), the most common experimental fluid, and for the solid interior of silicate planets ($Pr > 10^{21}$). In the latter, convection is in a laminar regime even at very large values of $Ra_H$ and it is appropriate to work in the limit of infinite $Pr$ (Schubert et al. 2001).

For values of $Ra_H$ up to $10^5$, steady-state motions are organized in periodic planforms that become increasingly complex as the Rayleigh–Roberts number is increased. For $Pr \approx 7$, these planforms have been explored systematically in laboratory experiments (Tritton & Zarraga 1967; Kulacki & Nagle 1975; Kulacki & Emara 1977; Takahashi et al. 2010) and numerical simulations (Ichikawa et al. 2006; Glover & Generalis 2009). Ordered by $Ra_H$ value, they take the form of squares, hexagons, both with an axial downwelling, and finally spokes, i.e. hexagons with an axial downwelling that is star-shaped instead of cylindrical. Their horizontal wavelength increases with $Ra_H$ in a phenomenon called ‘cell dilatation’. The same planforms are observed at larger $Pr$ values, as shown by both laboratory experiments ($3 \times 10^2 \le Pr \le 3 \times 10^4$) and numerical simulations in the infinite $Pr$ limit (Carrigan 1982, 1985; Houseman 1988; Parmentier & Sotin 2000; Limare et al. 2015), but it seems that cell dilatation does not occur. Moreover, their stability domains appear to be different from those for $Pr \approx 7$, but have not been determined precisely yet (Davaille & Limare 2015).

Very few studies have been made for intermediate values of $Ra_H$ in the $10^5-10^6$ range (Limare et al. 2015), despite the dramatic change of convection planform that is observed. Most of studies focus on cases with $Ra_H \gtrsim 10^6$ where convective motions become time-dependent. In this regime, scaling laws for the thickness of the upper thermal boundary layer and for the temperature difference across it have been derived from simple dimensional considerations and have been compared successfully to numerical and laboratory results (Cheung 1977; Kulacki & Emara 1977; Parmentier & Sotin 2000; Limare et al. 2015; Vilella & Kaminski 2017). In addition, scaling arguments and direct numerical simulations in the infinite $Pr$ limit indicate that the number of downwellings per unit area is proportional to $Ra_H^{1/4}$ (Parmentier & Sotin 2000).

Here, we aim at filling gaps in our understanding of laminar convection in internally heated horizontal fluid layers above an adiabatic base. This work builds on two previous studies (Limare et al. 2015; Vilella & Kaminski 2017). Limare et al. (2015) have compared the characteristics of convective flows in the laboratory at large $Pr$ values ($3 \times 10^2 < Pr < 3 \times 10^4$) and in direct numerical simulations in the infinite $Pr$ limit with exactly the same tank dimensions and temperature-dependent fluid properties. They have observed the same planforms in both types of studies and
have shown that the thickness and temperature difference across the upper boundary layer scale as $Ra_H^{-1/4}$ down to surprisingly small values of $Ra_H$ ($\approx 8 \times 10^5$). Vilella & Kaminski (2017) have focused on the upper thermal boundary layer and have shown that it hosts large lateral temperature variations even at high $Ra_H$ values, with important implications for magma generation in silicate planets. Using high-precision three-dimensional (3-D) direct numerical simulations in the infinite Pr limit, they have established scalings for both the horizontally averaged structure and the hottest regions of the upper thermal boundary layer. In this paper, we use the same numerical runs and add seven new ones to extract a comprehensive and consistent data set of thermal and flow characteristics, focusing on cold downwellings for the latter, over a very large range of Rayleigh-Roberts numbers ($5 \times 10^3-10^9$). We determine precisely the various planforms that occur, their domains of existence, their geometrical shapes and average dimensions. We document how the thermal structure of the fluid layer evolves as the Rayleigh-Roberts number increases. We also derive scaling laws for several variables of great geophysical interest, such as the average plume size and vertical velocity and the amplitude of thermal anomalies in the fluid interior, and compare them successfully to numerical results. We consider two sets of mechanical boundary conditions (free slip or rigid) in order to investigate how the interior convective motions depend on the dynamics of the unstable boundary layer at the top.

The paper is organized as follows. We begin by describing the various planforms of convection that are generated at the top of the layer and determine their respective domains of existence as a function of $Ra_H$. We investigate how motions in the fluid interior are related to the shallow planform. We then determine the main characteristics of the cold downwellings that are generated at the upper boundary and that go through the fluid layer and propose simple scaling laws for them. In a final section, we discuss some implications of our results for studies of the Earth and other planetary bodies.

## 2. Numerical simulations

We consider an incompressible fluid layer heated internally and cooled from above in the Boussinesq (1903) approximation and in the infinite Pr limit. The top boundary is kept at a constant temperature set to zero and the bottom one is adiabatic. We seek scaling laws that encapsulate the main physical aspects of the flow in compact form, which is best achieved for fluids with constant physical properties. Identical boundary conditions, which are either free slip or rigid, are imposed at the top and bottom. This is done to investigate the impact of boundary conditions on the dynamics of the thermal boundary layer (TBL) at the top and on the convective motions in the fluid interior. Using the following scales, $d$ for spatial coordinates, $d^2/\kappa$ for time, $\kappa/d$ for velocity, $\eta\kappa/d^2$ for pressure and finally $Hd^2/k$ for temperature, the conservation equations are written in dimensionless form as follows:

$$0 = \nabla \cdot \mathbf{v}, \qquad (2.1)$$

$$0 = -\nabla P + \nabla^2 \mathbf{v} + Ra_H T \mathbf{n}_z, \qquad (2.2)$$

$$\frac{DT}{Dt} = \frac{\partial T}{\partial t} + \mathbf{v} \cdot \nabla T = \nabla^2 T + 1, \qquad (2.3)$$

where $\mathbf{v}$, $P$ and $T$ are velocity, pressure and temperature, respectively, and where $t$ is time and $\mathbf{n}_z$ the unit vector in the vertical direction. Calculations were
<table><caption>TABLE 1. Input parameters for the 45 numerical simulations analysed in this paper: the Rayleigh–Roberts number ($Ra_H$), the number of grid elements used in X:Y:Z directions, the horizontal dimensions X:Y (scaled to the fluid layer thickness). $Nd^2$ is the dimensionless number of downwellings per unit area at mid-depth (see text). ‘NA’ indicates that it was not possible to obtain a reliable count of downwellings; ‘–’ indicates that no calculation was carried out.</caption><thead><tr><th>Ra<sub>H</sub></th><th>Grid size</th><th>Horizontal dimensions</th><th>Nd<sup>2</sup><br>rigid</th><th>Nd<sup>2</sup><br>free slip</th></tr></thead><tbody><tr><td>5 × 10<sup>3</sup></td><td>1024 × 1024 × 64</td><td>16:16</td><td>0.180</td><td>0.148</td></tr><tr><td>10<sup>4</sup></td><td>1024 × 1024 × 64</td><td>16:16</td><td>0.210</td><td>0.152</td></tr><tr><td>1.8 × 10<sup>4</sup></td><td>1024 × 1024 × 64</td><td>16:16</td><td>0.223</td><td>—</td></tr><tr><td>2.2 × 10<sup>4</sup></td><td>1024 × 1024 × 64</td><td>16:16</td><td>0.230</td><td>—</td></tr><tr><td>2.8 × 10<sup>4</sup></td><td>1024 × 1024 × 64</td><td>16:16</td><td>0.219</td><td>0.137</td></tr><tr><td>3.2 × 10<sup>4</sup></td><td>1024 × 1024 × 64</td><td>16:16</td><td>0.230</td><td>0.121</td></tr><tr><td>3.5 × 10<sup>4</sup></td><td>1024 × 1024 × 64</td><td>16:16</td><td>0.215</td><td>0.121</td></tr><tr><td>4.5 × 10<sup>4</sup></td><td>1024 × 1024 × 64</td><td>16:16</td><td>—</td><td>0.117</td></tr><tr><td>5 × 10<sup>4</sup></td><td>1024 × 1024 × 64</td><td>16:16</td><td>0.221</td><td>0.117</td></tr><tr><td>10<sup>5</sup></td><td>1024 × 1024 × 64</td><td>16:16</td><td>NA</td><td>NA</td></tr><tr><td>1.4 × 10<sup>5</sup></td><td>1024 × 1024 × 64</td><td>16:16</td><td>NA</td><td>NA</td></tr><tr><td>2.2 × 10<sup>5</sup></td><td>512 × 512 × 64</td><td>6:6</td><td>NA</td><td>NA</td></tr><tr><td>3.2 × 10<sup>5</sup></td><td>512 × 512 × 64</td><td>6:6</td><td>NA</td><td>NA</td></tr><tr><td>6.8 × 10<sup>5</sup></td><td>512 × 512 × 64</td><td>6:6</td><td>1.14</td><td>1.12</td></tr><tr><td>10<sup>6</sup></td><td>512 × 512 × 64</td><td>6:6</td><td>1.61</td><td>1.19</td></tr><tr><td>1.6 × 10<sup>6</sup></td><td>512 × 512 × 64</td><td>6:6</td><td>1.87</td><td>1.37</td></tr><tr><td>3.2 × 10<sup>6</sup></td><td>512 × 512 × 64</td><td>6:6</td><td>2.38</td><td>1.61</td></tr><tr><td>5 × 10<sup>6</sup></td><td>512 × 512 × 64</td><td>6:6</td><td>2.70</td><td>1.82</td></tr><tr><td>6.3 × 10<sup>6</sup></td><td>512 × 512 × 64</td><td>6:6</td><td>2.84</td><td>1.85</td></tr><tr><td>7.9 × 10<sup>6</sup></td><td>512 × 512 × 64</td><td>6:6</td><td>3.08</td><td>1.97</td></tr><tr><td>1.4 × 10<sup>7</sup></td><td>512 × 512 × 64</td><td>6:6</td><td>3.46</td><td>2.25</td></tr><tr><td>10<sup>8</sup></td><td>384 × 384 × 128</td><td>6:6</td><td>6.49</td><td>3.71</td></tr><tr><td>3.2 × 10<sup>8</sup></td><td>384 × 384 × 128</td><td>6:6</td><td>9.17</td><td>4.68</td></tr><tr><td>10<sup>9</sup></td><td>512 × 512 × 256</td><td>4:4</td><td>14.5</td><td>6.87</td></tr></tbody></table>

performed in 3-D Cartesian geometry using Stag3D (Tackley 1998a, 2008), a robust code that has been tested thoroughly against several others. Time-stepping is implemented for the heat equation using a finite volume formulation and the mass and momentum conservation equations are solved using a finite difference multigrid technique. The accuracy of this code was assessed by a detailed comparison with laboratory experiments (Limare et al. 2015) and precision tests are described in the supplementary material available at https://doi.org/10.1017/jfm.2018.316. We have run a large number of numerical simulations in order to determine the characteristics of convection over a wide $Ra_H$ range (table 1). We made sure that the calculations allow excellent resolution of the temperature and velocity fields by decreasing the grid spacing as we increased the Rayleigh–Roberts number. At $Ra_H = 10^9$, the largest value of the Rayleigh–Roberts number studied, there are nine grid points through the thermal boundary layer. We also made sure that the horizontal resolution was optimal, such that the intricacies of the convection planform were adequately accounted for. We show in the supplementary material that results do not depend on grid size. In order to minimize lateral edge effects, computational domains had high aspect
ratios (width over thickness) and reflecting vertical boundaries. The smallest aspect ratio was equal to 4 and was only used for calculations at the largest Ra<sub>H</sub> value (Ra<sub>H</sub> = 10<sup>9</sup>). In this case, convective motions develop over small horizontal spatial scales, such that the calculations allowed a representative statistical sampling of the complex flow field and such that the impact of the lateral boundaries was limited to narrow peripheral regions. In all other cases, aspect ratios were set at 6 or more (table 1). In a few cases, we checked that calculations for different aspect ratios lead to essentially identical results (see the supplementary material for a detailed comparison).

Calculations are time-dependent and are started from different initial conditions depending on the Ra<sub>H</sub> value. At low Ra<sub>H</sub> associated with steady-state regimes, the whole fluid layer is initially at a constant temperature with random, small amplitude, perturbations. This ensures that the convective planform of the steady-state solution does not depend on an arbitrary initial spatial pattern. At high Ra<sub>H</sub>, truly steady-state regimes cannot be achieved and time series of the local values of temperature and velocity exhibit fluctuations at all times. We therefore consider ‘statistical steady-state’ conditions, such that temperature and heat flux values remain within the same ranges over several ‘overturn times’. The overturn time is defined to be the time taken by a fluid parcel to go full circle, from the upper boundary layer to the base of the domain and then back up. For these calculations (high Ra<sub>H</sub>), the initial condition is the temperature field for a solution at a lower Ra<sub>H</sub> value, a strategy that significantly reduces the computation time. We checked that the final result does not depend on this initial condition in a few cases.

### 3. Convective planforms

#### 3.1. Preliminary note

Before we describe the convective planforms in more detail, it is useful to recall how the characteristics of the upper thermal boundary layer depend on the model parameters, especially on the Rayleigh–Roberts number and on the boundary conditions. For 10<sup>5</sup> ≤ Ra<sub>H</sub> ≤ 10<sup>9</sup>, the temperature difference and thickness of the boundary layer, noted ΔT<sub>TBL</sub> and δ, respectively, are such that

$$ \frac{\Delta T_{TBL}}{Hd^2/k} = C_{\Delta T} Ra_H^{-1/4} \quad \text{and} \quad \frac{\delta}{d} = C_{\delta} Ra_H^{-1/4}, \qquad (3.1a,b) $$

where C<sub>ΔT</sub> and C<sub>δ</sub> are proportionality constants that depend on the boundary conditions (Limare et al. 2015, and see table 2). These relationships can be derived from local scaling arguments (Vilella & Kaminski 2017), indicating that the dynamics of the upper boundary layer is locally controlled, independently of the total fluid layer thickness. Note that thickness δ corresponds to the base of the boundary layer, which is defined with the horizontally averaged temperature at the point where the convective flux is maximum. An alternative definition relies on the heat flux through the upper boundary, such that q = ΔT<sub>TBL</sub>/δ* in dimensionless variables. In steady state, heat conservation requires that q = 1, implying that δ* = ΔT<sub>TBL</sub>. Limare et al. (2015) chose the former definition in order to verify that the vertical temperature profile in the boundary layer conforms to local scalings. By contrast, for Ra<sub>H</sub> ≤ 10<sup>5</sup> the relationship (3.1) is no longer valid, indicating an interdependence between the upper boundary layer and the fluid interior. Another key result is that regardless of the Ra<sub>H</sub> value, beneath a rigid top, the boundary layer is significantly thicker and hosts a larger temperature difference than beneath a free slip surface. Variations of behaviour
<table><caption>TABLE 2. Proportionality constants in the scaling laws for the temperature difference and thickness of the upper thermal boundary layer (3.1) from Limare et al. (2015).</caption><thead><tr><th>Boundary conditions</th><th>C<sub>ΔT</sub></th><th>C<sub>δ</sub></th></tr></thead><tbody><tr><td>Free slip</td><td>2.49</td><td>5.90</td></tr><tr><td>Rigid</td><td>3.41</td><td>7.08</td></tr></tbody></table>

observed when changing the boundary conditions as well as the Ra$_H$ value are likely to affect the convective planforms. These differences will be identified and quantified in the following sections.

## 3.2. Steady-state regimes

By definition, steady-state regimes are such that the geometrical pattern of convection does not change over time. It is difficult, however, to prove that true steady state has been achieved in numerical simulations, since it would require extremely long computation times. We consider that the solution is steady-state when the pattern of convection does not change significantly over several overturn times (which has been defined above). This requires a dimensionless computational time 5–50 with a typical value of ≈10.

The analysis is carried out in a horizontal plane located at dimensionless height z=0.9 above the base. The flow field can be split into a number of ‘convective cells’, which are defined by one downwelling and its ‘capture’ area, i.e. the area in the upper boundary layer where fluid flows laterally towards the downwelling. For visualization purposes, we also use a Voronoi procedure in which each cell is built around one downwelling and includes all the points that are closer to it than to any other one.

We first discuss results obtained for free slip boundaries. Near the threshold value of Ra$_H$ for the onset of convection (Ra$_{H,cr}$ = 868), the theoretical analysis by Roberts (1967) indicates that three planforms are possible, described as 2-D rolls and hexagons with either an axial upwelling or an axial downwelling. By analogy with Rayleigh-Bénard systems, we may expect that hexagons with an axial downwelling are preferred over other ones due to the different thermal boundary conditions at the top and bottom (Cross & Hohenberg 1993). Our calculations confirm this (figures 1a and 2). The preferential pattern is made of hexagons with a cylindrical downwelling at their centre, which have been called ‘down’ hexagons, but there are also a few pentagons and rare squares and heptagons.

As Ra$_H$ increases, the planform grades into a spoke pattern (figure 1b–d). This pattern is still made of hexagons, but features an axial downwelling that is star-shaped instead of cylindrical, with six spokes perpendicular to the six sides of the encasing hexagon. The planform evolves progressively as Ra$_H$ increases, such that spokes extend to increasing distances from the hexagon centres. This makes it difficult to determine a threshold value for the Rayleigh-Roberts number that marks the appearance of this planform unambiguously. We have set the transition at the stage when the downwellings outlines become noticeably spiky (Ra$_H$ ≈ 2.8 × 10$^4$, figure 1c). This is obviously an approximate procedure, but the range of Ra$_H$ values for the transition is less than 0.2 log units, as may be seen from table 1 and figure 1(b,c).

Increasing Ra$_H$ further, the spokes become increasingly prominent and their tips eventually reach the sides of the encasing hexagons. This leads to a spectacular change of flow pattern that marks the onset of time-dependent convection (Ra$_H$ ≈ 10$^5$ figure 1e,f).
<figure><img src="image_94.png" /></figure>

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

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

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

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

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

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

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

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

FIGURE 1. Convective planform as a function of the Rayleigh–Roberts number ($Ra_H$) for free slip boundaries at the top ($z=1$) and bottom ($z=0$). The temperature field is shown for a horizontal plane located at $z=0.9$. Blue and red colours correspond to fluid that is colder and hotter than the horizontal average at that depth, respectively. The colour scale is not the same for all the panels in order to enhance visibility of the convective structure. The resolution and the domain aspect ratio are changed for panels (g–i) (table 1) due to the small dimensions of the convective cells. The black bar beneath each panel corresponds to the thickness of the fluid layer, i.e. the height of the computation domain.

### 3.3. The ‘sheet’ regime

For $Ra_H \gtrsim 10^5$, the flow adopts a new configuration that is observed over a range of $Ra_H$ spanning approximately 0.3 log units for free slip boundaries and about 0.7 log units for rigid ones. Instead of continuous upwelling regions encircling individual downwellings, the flow gets organized in arrays of thin line downwellings that encapsulate broad upwellings (figure 1f). These downwellings are no longer separated from one another by continuous upwelling areas and will be referred to as ‘sheets’. They converge to three-way junction nodes and outline a pattern suggestive
<figure><img src="image_103.png" /></figure>

FIGURE 2. Convective planform defined from the temperature field in figure 1(a) for Ra$_H$ = 5 × 10$^3$ (a) and in figure 1(e) for Ra$_H$ = 5 × 10$^4$ (b). Temperatures are shown as deviations from the local horizontal average, with blue and red standing for temperature that is colder and hotter than average, respectively. Black lines in both panels are the Voronoi tessellation, which indicates a dominantly hexagonal pattern. Thick blue lines in (b) illustrate how hexagons may form out of a spoke pattern, in that the case the centre downwelling would no longer be present.

of hexagons, but with a polarity that is the opposite of that for smaller values of
the Rayleigh-Roberts number: they now occur along the sides of polygonal cells
instead of at their centre, which was described as a switch from ‘down’ hexagons to
‘up’ hexagons by Carrigan (1985). This polarity reversal is associated with a marked
increase of cell dimensions. This regime, which will be called the sheet regime, is
of particular interest for reasons that will be explained later. Sheets are similar to
‘line plumes’ emanating from elongated regions in an unstable boundary layer, which
are the predominant flow structure in turbulent Rayleigh-Bénard convection (Adrian,
Ferreira & Boberg 1986; Theerthan & Arakeri 1994; Shishkina & Wagner 2008;
Zhou & Xia 2010).

Two features make the sheet regime distinctive. One is that the length and spacing of sheets exceed the depth of the fluid layer by large amounts. It is unfortunately difficult to determine a representative length scale due to the irregular arrangement of the sheets, their variable lengths and their small number even in a domain that extends laterally over sixteen times the layer thickness. The other feature is that it is a time-dependent pattern, albeit one that evolves very slowly, such that the sheets move horizontally across the domain. The ‘drift’ velocity is smaller than the average vertical velocity by at least one order of magnitude. The sheet regime is therefore neither truly steady-state nor truly time-dependent. Nevertheless, for the sake of simplicity, we consider in the following the sheet regime as a steady-state regime.

3.4. *Time-dependent flow regimes*

For free slip boundaries, the sheet regime eventually segues into a time-dependent
one with no obvious geometrical pattern when Ra<sub>H</sub> ≈ 3.2 × 10<sup>5</sup> (figure 1h). The
planform is shown at a dimensionless depth of 0.1 below the top, in the upper part
of the boundary layer (δ = 0.25 in this particular case), corresponding to what would
be seen by an external observer. It is ill-defined at that level but we shall see later that it changes as one probes deeper in the fluid and one can determine a meaningful spacing between adjacent downwellings. Individual sheet segments are unstable and develop into a number of discrete focused downwellings (figure 1g). The lifetime of an individual sheet is smaller than the time scale of its lateral drift motion by approximately two orders of magnitude. Sheets are rarely isolated and belong to clusters made of nodal points and branches. With three sheets per nodal point on average, the pattern can be interpreted as a set of ‘truncated hexagons’ (figure 1i for Ra$_H$ = $10^8$).

### 3.5. Rigid boundaries

We use the same procedure and descriptive terms to characterize convective planforms between rigid boundaries. One notable difference is that downwellings are in larger numbers than in cases with free slip boundaries (figures 1 and 3). They are also wider. This second difference is, however, not apparent in figures 1 and 3, which illustrate the geometrical planforms at shallow depth, and will be shown by measurements of the cross-sectional area of downwellings in horizontal planes at mid-depth in the fluid interior. Downwellings are fed by laterally converging flows that extend through the upper boundary layer, and their size reflects the thickness of that boundary layer, which is thicker beneath a rigid boundary (table 2).

The planforms are essentially identical to those for free slip boundaries but the critical Ra$_H$ values for transitions are not the same (figure 4). This is not surprising as the critical value of Ra$_H$ for the onset of convection is larger for rigid boundaries (Ra$_{H,cr}$ = 2772 following Roberts 1967) than for free slip ones (Ra$_{H,cr}$ = 868 following Kulacki & Goldstein 1975). The spoke and sheet regimes are observed over a wider Ra$_H$ range than with free slip boundaries. The sheet regime is heralded by a gradual change of planform illustrated in figure 3(d,e), which is such that gaps between neighbouring spokes get bridged progressively, leading to a polygonal array of line downwellings. As in cases with free boundaries, downwellings go from isolated plumes at cell centres to a continuous network of thin sheets encircling broad upwellings (figures 3f,g). This can be described again as a switch from ‘down’ to ‘up’ hexagons with a marked increase of cell dimensions. This transition occurs at Ra$_H$ ≈ 40 Ra$_{H,cr}$ ≈ $10^5$, which is consistent with the careful laboratory observations of Carrigan (1985).

For Ra$_H$ ≥ $10^6$ (figure 3i), the planform can be described again as a set of truncated hexagons, perhaps more convincing than in the free slip calculations. The vast majority of the downwelling sheets belong to networks that can be traced over large horizontal distances, with nodal points where three sheets converge. Such three-way junctions are consistent with a hexagonal pattern, but there are also a very small number of four-way junctions. Individual sheets typically undergo the following cycle. A new downwelling forms in the upper boundary layer and gets stretched towards an already existing sheet due to the horizontal flow that feeds it. Once several (generally three) downwelling sheets are joined in this manner, downward flow is amplified in the nodal area and draws fluid from the sheets, which coalesce and eventually disappear. This can be described as a hexagon-building process that does not go to completion. An important feature is that sheets are not stable structures as they go down in the fluid, as illustrated in the following section.

### 3.6. Vertical structure of downwellings

We have so far chosen to document the convective planforms near the top of the fluid layer because this is what an observer can see from ‘outside’ (i.e. with
<figure><img src="image_104.png" /></figure>

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

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

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

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

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

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

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

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

FIGURE 3. Convective planform as a function of the Rayleigh–Roberts number ($Ra_H$) for rigid boundaries at the top ($z = 1$) and bottom ($z = 0$). The temperature field is shown for a horizontal plane located at $z = 0.9$. Blue and red colours correspond to fluid that is colder and hotter than the horizontal average at that depth, respectively. The colour scale is not the same for all the panels in order to enhance visibility of the convective structure. The resolution and the domain aspect ratio are changed for panels (g–i) (table 1) due to the small size of convective cells. The black bar beneath each panel corresponds to the thickness of the fluid layer, i.e. the box height. Note that the typical dimensions of the convective cells that are shown in (f) and (g) are very close to one another when scaled to layer thickness.

observations of the planet surface), but this may not be sufficient for studies of deep planetary motions. In steady-state regimes, one expects that the same basic planform is maintained through the layer interior, save for the intricate small-scale fabric of the spoke pattern. In time-dependent regimes at high values of $Ra_H$, we may expect that the local arrays of thin downgoing sheets that prevail at shallow depth are not stable over a large vertical extent. To address these issues, we have studied the convective planforms at several depths in the fluid layer.
<figure><img src="image_113.png" /><figcaption>FIGURE 4. Diagram showing the change of convective planform as a function of the Rayleigh–Roberts number ($Ra_H$) for both rigid (black bar) and free slip (red bar) boundaries. Threshold values of $Ra_H$ for the onset of convection come from stability analysis (Roberts 1967; Kulacki & Goldstein 1975) and CR stands for conductive regime, whereas the other transitions are based on the present numerical simulations (figures 1 and 3).</figcaption></figure>

Figure 5 shows how downwellings evolve with increasing depth for cases spanning the whole range of Rayleigh–Roberts number investigated. The same changes are observed for both types of boundaries and are illustrated only for free slip ones. In the steady-state hexagon (figure 5a) and sheet regimes (figure 5c), downwellings keep the same spacing and shape at all depths, save for the very base of the fluid layer where they spread laterally. In the spoke regime (figure 5b), however, the star-shaped structures that prevail at shallow depths gradually coarsen into nearly cylindrical ones. In time-dependent regimes at high values of $Ra_H$ (figure 5d,e), sheets split and merge into a few individual downwellings with nearly equidimensional horizontal cross-sections.

In steady-state regimes, the planform is essentially the same at all depths, indicating that the flow structure is stable over the whole thickness of fluid. This contrasts with the downward change that is observed in time-dependent regimes at high $Ra_H$ values. In these cases, the boundary layer characteristics can be derived from local arguments (Parmentier & Sotin 2000; Limare et al. 2015; Vilella & Kaminski 2017), suggesting that the shallow convective planform is determined locally, rather than by dynamics at the scale of the whole layer. The transition from a global to a local control of the shallow convective planform is correlated with a change in the behaviour of sheets, which extend to the base of the fluid layer in the regime of the same name (figure 5c) but become unstable when the Rayleigh–Roberts number exceeds a threshold value. As sheets get destabilized, they feed nearly equidimensional and isolated descending plumes that go through the whole fluid layer (figure 5d,e).

## 4. Thermal structure of convection in the fluid interior

In the following, it will be useful to break down temperature as the sum of a horizontal average, $T_{avg}$, and a fluctuation, $\theta$, also called an ‘anomaly’ to follow common geophysical terminology. For clarity purposes, temperatures will be made dimensionless using either the bulk temperature scale $\Delta T_H = Hd^2/k$ or the local horizontal average $T_{avg}$.

### 4.1. Temperature distribution

The different regimes that have been observed are associated with different temperature distributions that reflect changes in the balance between downwellings and upwellings.
<figure><img src="image_114.png" /><figcaption>FIGURE 5. Temperature field (T) within and around downwellings at different heights above base ($z = 1$ is at the top) for different values of the Rayleigh–Roberts number ($Ra_H$) in a fluid layer with free slip boundaries. All the panels extend over three times the fluid layer thickness in both horizontal directions. For a given Rayleigh–Roberts number, all panels are centred at the same horizontal coordinates in order to track the downward change of downwelling structure.</figcaption></figure>

Figure 6 illustrates how the horizontal temperature distribution varies with depth for four representative cases. In all cases, independently of the $Ra_H$ value, the spread of temperatures increases from zero at the top to a maximum at the base of the upper boundary layer, and then decreases towards the bottom of the fluid layer. This is due to internal heating, as shown by the steady increase of the minimum temperature with increasing depth (the cold profile in figure 6). In the steady-state hexagon regime at $Ra_H = 5 \times 10^3$, temperatures are spread across the range due to the nearly equal areas occupied by upwellings and downwellings. In contrast, the temperature distribution
<figure><img src="image_115.png" /><figcaption>FIGURE 6. Horizontal distribution of temperature as a function of depth for different values of the Rayleigh–Roberts number ($Ra_H$) in a layer with free slip boundaries. Temperatures have been made dimensionless with the bulk scale $\Delta T_H = Hd^2/k$. Solid black lines correspond to the minimum and maximum temperatures recorded at each depth ('cold' and 'hot' temperature profiles), while solid light grey lines correspond to vertical profiles of the horizontally averaged temperature. Temperature values at each depth are shown on a white to black colour scale depending on their probability of occurrence. At $Ra_H = 10^9$, due to the overwhelming areal dominance of upwellings, the vertical profile of the most probable temperature is essentially identical to that of the horizontally averaged temperature.</figcaption></figure>
<figure><img src="image_116.png" /><figcaption>FIGURE 7. Horizontally averaged temperature profile obtained for two different values of the Rayleigh–Roberts number ($Ra_H$) with either free slip (solid red lines) or rigid (solid black lines) boundary conditions.</figcaption></figure>

in the spoke regime at $Ra_H = 5 \times 10^4$ is markedly peaked due to the large areas of upward return flow where lateral gradients are small (figure 5b). ‘Peaked’ distributions are also a feature of time-dependent regimes, with a temperature range that decreases with increasing $Ra_H$. In these cases, the flow can be described as a set of narrow downwellings going through nearly isothermal fluid. The gradual change from spread out to increasingly peaked distributions is interrupted in the sheet regime, as shown for $Ra_H = 1.4 \times 10^5$ in figure 6, which is yet another peculiarity of this regime.

## 4.2. Vertical distribution of the horizontally averaged temperature

One significant feature of convection in internally heated fluid layers is that the horizontally averaged temperature $T_{avg}$ decreases significantly with increasing depth below the upper boundary layer (figures 6 and 7). This is in marked contrast to Rayleigh–Bénard convection where the decrease of $T_{avg}$ only occurs for a narrow vertical extent at the base of the boundary layer. We find that the magnitude of the vertical gradient of $T_{avg}$ in the fluid interior decreases with increasing $Ra_H$, as already noted by Parmentier & Sotin (2000) and Moore (2008).

Vertical profiles of the horizontally averaged temperature are shown in figure 7 for two values of $Ra_H$ ($10^6$ and $10^8$) and for both types of boundary conditions. In all cases, the temperature at the bottom of the fluid is significantly smaller than at the base of the upper boundary layer, implying that the fluid interior is stably stratified. As $Ra_H$ increases, the vertical gradient of $T_{avg}$ at mid-depth decreases in parallel with the temperature contrast across the upper boundary layer. One notes again the important influence of the mechanical boundary conditions, which affect the magnitudes of both the temperature difference across the layer and the interior
<table><caption>TABLE 3. Some characteristics of convection for the 45 numerical simulations carried out for this study. $T_{vol}k/Hd^2$ is the volume-averaged temperature and $T_{bot}k/Hd^2$ is the basal temperature. The other variables correspond to the average properties of downwellings at mid-depth: $A_i/d^2$ is the average horizontal cross-section, $\Delta T_i k/Hd^2$ is the average temperature contrast and $W_i d/\kappa$ is the average vertical velocity. ‘NA’ indicates that it was not possible to obtain a reliable count of downwellings; ‘—’ indicates that no calculation was carried out.</caption><thead><tr><th rowspan="2">Ra<sub>H</sub></th><th colspan="2">T<sub>vol</sub>k/Hd<sup>2</sup></th><th colspan="2">T<sub>bot</sub>k/Hd<sup>2</sup></th><th colspan="2">A<sub>i</sub>/d<sup>2</sup></th><th colspan="2">ΔT<sub>i</sub>k/Hd<sup>2</sup></th><th colspan="2">W<sub>i</sub>d/κ</th></tr><tr><th>Rigid</th><th>Free slip</th><th>Rigid</th><th>Free slip</th><th>Rigid</th><th>Free slip</th><th>Rigid</th><th>Free slip</th><th>Rigid</th><th>Free slip</th></tr></thead><tbody><tr><td>5 × 10<sup>3</sup></td><td>0.269</td><td>0.188</td><td>0.368</td><td>0.209</td><td>2.33</td><td>1.89</td><td>0.558</td><td>0.372</td><td>4.26</td><td>10.6</td></tr><tr><td>10<sup>4</sup></td><td>0.221</td><td>0.163</td><td>0.275</td><td>0.163</td><td>1.80</td><td>1.36</td><td>0.461</td><td>0.324</td><td>8.31</td><td>19.4</td></tr><tr><td>1.8 × 10<sup>4</sup></td><td>0.198</td><td>—</td><td>0.228</td><td>—</td><td>1.48</td><td>—</td><td>0.422</td><td>—</td><td>12.4</td><td>—</td></tr><tr><td>2.2 × 10<sup>4</sup></td><td>0.191</td><td>—</td><td>0.215</td><td>—</td><td>1.36</td><td>—</td><td>0.410</td><td>—</td><td>14.9</td><td>—</td></tr><tr><td>2.8 × 10<sup>4</sup></td><td>0.185</td><td>0.143</td><td>0.203</td><td>0.126</td><td>1.37</td><td>0.773</td><td>0.399</td><td>0.284</td><td>17.4</td><td>50.7</td></tr><tr><td>3.2 × 10<sup>4</sup></td><td>0.181</td><td>0.142</td><td>0.199</td><td>0.124</td><td>1.25</td><td>0.790</td><td>0.392</td><td>0.286</td><td>19.5</td><td>54.6</td></tr><tr><td>3.5 × 10<sup>4</sup></td><td>0.180</td><td>0.141</td><td>0.193</td><td>0.121</td><td>1.28</td><td>0.709</td><td>0.392</td><td>0.281</td><td>20.4</td><td>62.8</td></tr><tr><td>4.5 × 10<sup>4</sup></td><td>—</td><td>0.138</td><td>—</td><td>0.116</td><td>—</td><td>0.615</td><td>—</td><td>0.275</td><td>—</td><td>78.5</td></tr><tr><td>5 × 10<sup>4</sup></td><td>0.172</td><td>0.137</td><td>0.181</td><td>0.114</td><td>1.07</td><td>0.570</td><td>0.377</td><td>0.272</td><td>27.9</td><td>87.0</td></tr><tr><td>10<sup>5</sup></td><td>0.157</td><td>0.118</td><td>0.161</td><td>0.106</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>1.4 × 10<sup>5</sup></td><td>0.146</td><td>0.107</td><td>0.152</td><td>0.0984</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>2.2 × 10<sup>5</sup></td><td>0.132</td><td>0.0977</td><td>0.140</td><td>0.0876</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>3.2 × 10<sup>5</sup></td><td>0.124</td><td>0.0867</td><td>0.131</td><td>0.0759</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td><td>NA</td></tr><tr><td>6.8 × 10<sup>5</sup></td><td>0.103</td><td>0.0722</td><td>0.107</td><td>0.0628</td><td>0.133</td><td>0.0730</td><td>0.209</td><td>0.134</td><td>127</td><td>191</td></tr><tr><td>10<sup>6</sup></td><td>0.0944</td><td>0.0657</td><td>0.0966</td><td>0.0565</td><td>0.0913</td><td>0.0647</td><td>0.189</td><td>0.121</td><td>148</td><td>221</td></tr><tr><td>1.6 × 10<sup>6</sup></td><td>0.0854</td><td>0.0588</td><td>0.0869</td><td>0.0503</td><td>0.0804</td><td>0.0550</td><td>0.169</td><td>0.108</td><td>170</td><td>260</td></tr><tr><td>3.2 × 10<sup>6</sup></td><td>0.0734</td><td>0.0488</td><td>0.0742</td><td>0.0412</td><td>0.0621</td><td>0.0446</td><td>0.145</td><td>0.0890</td><td>213</td><td>343</td></tr><tr><td>5 × 10<sup>6</sup></td><td>0.0661</td><td>0.0446</td><td>0.0666</td><td>0.0373</td><td>0.0557</td><td>0.0379</td><td>0.130</td><td>0.0808</td><td>250</td><td>394</td></tr><tr><td>6.3 × 10<sup>6</sup></td><td>0.0628</td><td>0.0422</td><td>0.0632</td><td>0.0352</td><td>0.0526</td><td>0.0359</td><td>0.124</td><td>0.0765</td><td>271</td><td>432</td></tr><tr><td>7.9 × 10<sup>6</sup></td><td>0.0597</td><td>0.0398</td><td>0.0599</td><td>0.0331</td><td>0.0484</td><td>0.0327</td><td>0.117</td><td>0.0721</td><td>291</td><td>469</td></tr><tr><td>1.4 × 10<sup>7</sup></td><td>0.0527</td><td>0.0349</td><td>0.0527</td><td>0.0287</td><td>0.0428</td><td>0.0283</td><td>0.103</td><td>0.0630</td><td>349</td><td>575</td></tr><tr><td>10<sup>8</sup></td><td>0.0336</td><td>0.0218</td><td>0.0333</td><td>0.0176</td><td>0.0231</td><td>0.0143</td><td>0.0660</td><td>0.0393</td><td>656</td><td>1184</td></tr><tr><td>3.2 × 10<sup>8</sup></td><td>0.0252</td><td>0.0165</td><td>0.0248</td><td>0.0132</td><td>0.0154</td><td>0.00958</td><td>0.0495</td><td>0.0297</td><td>981</td><td>1860</td></tr><tr><td>10<sup>9</sup></td><td>0.0195</td><td>0.0126</td><td>0.0191</td><td>0.00990</td><td>0.00856</td><td>0.00560</td><td>0.0383</td><td>0.0229</td><td>1480</td><td>2856</td></tr></tbody></table>

thermal gradient. They also imply different vertical temperature profiles in the lower parts of the fluid layer, where a sub-layer connects the stably stratified interior to the adiabatic lower boundary. There, rigid behaviour impedes the lateral spreading of cold fluid coming from downwellings, leading to a thicker sub-layer than in the free slip case. In order to characterize how the bulk thermal structure varies with the Rayleigh–Roberts number, we have determined the volume-averaged and basal temperatures, noted $T_{vol}$ and $T_{bot}$, respectively, in all our numerical simulations. These results are reported in table 3 together with others that will be discussed in another section below.

Values of both $T_{vol}$ and $T_{bot}$ are shown as functions of the Rayleigh–Roberts number in figure 8 for all simulations. Overall, both temperatures decrease steadily with increasing $Ra_H$ with a marked change of trend for $Ra_H \approx 10^5$, in conjunction with the occurrence of the sheet regime. For $Ra_H \ge 6.8 \times 10^5$, the data appear to conform to power-law relationships with $Ra_H$ and best-fit parameter values are reported in table 4. The power-law exponents take values that are very close to $-1/4$, which is the same value as that for the temperature contrast across the upper boundary
<figure><img src="image_117.png" /></figure>

FIGURE 8. Dependence of (a) the volume-averaged temperature ($T_{vol}/k/Hd^2$) and (b) the horizontally averaged basal temperature ($T_{bot}/k/Hd^2$) as a function of the Rayleigh-Roberts number ($Ra_H$). Results are shown for both free slip (red) and rigid (black) boundaries. Symbols indicate the convection planform that is achieved: a hexagonal pattern (square), a spoke pattern (star), a sheet regime (diamond) and the truncated hexagonal pattern (circle). Error bars, often smaller than the symbol size, indicate variations that are observed in time-dependent regimes. Solid lines are best fits to the data with power-law exponents set equal to -1/4. Values for the various constants in the power-law relationships are listed in tables 4 and 5. Data with open symbols have not been included in the best-fit procedure.

layer. We thus surmise that $T_{vol}$ and $T_{bot}$ also scale as $Ra_H^{-1/4}$ and have implemented another best-fit procedure with the value of the power-law exponent set to -1/4. As shown in table 4, this allows excellent fits to the data over a $Ra_H$ range spanning more than three orders of magnitude, with deviations that are typically less than ±1%. Results for the volume-averaged temperature are slightly less satisfactory for rigid boundaries but remain within ±4% of the simple -1/4 scaling law. For the sake of completeness, we have also determined best-fit power-law exponents for the steady-state regimes (table 5) but they do not account for the data very well as can be seen from the large errors.

### 4.3. *Temperature fluctuations*

Temperature fluctuations at mid-depth in the layer are calculated as departures from the horizontal average at the same depth scaled to that average. Their distributions are strongly asymmetric, as shown in figure 9, in marked contrast to those of Rayleigh-Bénard convection, which are very close to being symmetrical (and Gaussian in the 'soft turbulence' regime, e.g. Castaing *et al.* 1989). Their spread of values decreases as $Ra_H$ increases. In time-dependent regimes at large $Ra_H$ values (figure 9b), they are markedly peaked with a slightly positive peak value that is due to the broad and rather homogeneous passive return flows that surround downgoing plumes. The distribution for the sheet regime stands out of all the others with its two peaks, corresponding to two well-separated populations for upwellings and downwellings (figure 9a).
<table><caption>TABLE 4. Parameters of best-fit power laws for the volume-averaged ($T_{vol}k/Hd^2$) and basal ($T_{bot}k/Hd^2$) temperatures for time-dependent regimes at large values of the Rayleigh-Roberts number ($Ra_H$). Results in the right-hand column have been obtained by setting the power-law exponents to $-1/4$.</caption><caption>Ra<sub>H</sub> &ge; 6.8 &times; 10<sup>5</sup></caption><thead><tr><th>Variable</th><th>Boundary conditions</th><th>Exponent left to vary</th><th>Exponent fixed</th></tr></thead><tbody><tr><td rowspan="2">T<sub>vol</sub>k/Hd<sup>2</sup></td><td>free slip</td><td>1.796(&plusmn;0.016) Ra<sub>H</sub><sup>-0.240(&plusmn;0.001)</sup></td><td>2.127(&plusmn;0.035) Ra<sub>H</sub><sup>-1/4</sup></td></tr><tr><td>rigid</td><td>2.267(&plusmn;0.125) Ra<sub>H</sub><sup>-0.229(&plusmn;0.003)</sup></td><td>3.191(&plusmn;0.111) Ra<sub>H</sub><sup>-1/4</sup></td></tr><tr><td rowspan="2">T<sub>bot</sub>k/Hd<sup>2</sup></td><td>free slip</td><td>1.839(&plusmn;0.053) Ra<sub>H</sub><sup>-0.252(&plusmn;0.002)</sup></td><td>1.770(&plusmn;0.009) Ra<sub>H</sub><sup>-1/4</sup></td></tr><tr><td>rigid</td><td>2.503(&plusmn;0.108) Ra<sub>H</sub><sup>-0.235(&plusmn;0.003)</sup></td><td>3.201(&plusmn;0.080) Ra<sub>H</sub><sup>-1/4</sup></td></tr></tbody></table>

<table><caption>TABLE 5. Parameters of empirical best-fit power laws for the volume-averaged ($T_{vol}k/Hd^2$) and basal ($T_{bot}k/Hd^2$) temperatures in steady-state regimes at low values of the Rayleigh-Roberts number ($Ra_H$).</caption><caption>Ra<sub>H</sub> &le; 5 &times; 10<sup>4</sup></caption><thead><tr><th>Variable</th><th>Boundary conditions</th><th></th></tr></thead><tbody><tr><td rowspan="2">T<sub>vol</sub>k/Hd<sup>2</sup></td><td>free slip</td><td>0.5785(&plusmn;0.290) Ra<sub>H</sub><sup>-0.135(&plusmn;0.025)</sup></td></tr><tr><td>rigid</td><td>1.329(&plusmn;0.438) Ra<sub>H</sub><sup>-0.192(&plusmn;0.033)</sup></td></tr><tr><td rowspan="2">T<sub>bot</sub>k/Hd<sup>2</sup></td><td>free slip</td><td>1.848(&plusmn;0.704) Ra<sub>H</sub><sup>-0.260(&plusmn;0.037)</sup></td></tr><tr><td>rigid</td><td>4.891(&plusmn;2.262) Ra<sub>H</sub><sup>-0.309(&plusmn;0.090)</sup></td></tr></tbody></table>

## 5. The characteristics of downwellings

In an internally heated fluid layer, downwellings occur as discrete entities associated with large deviations from the horizontally averaged temperature. We therefore focus on them to describe convection and determine some of their characteristics, such as the average number per unit area, the horizontal cross-section, the average temperature and the vertical velocity. As discussed above, convective planforms may change with increasing depth below the boundary layer and we shall restrict our analysis to structures at mid-depth in the fluid interior, away from the upper thermal boundary layer. Results for the two types of boundary conditions are qualitatively similar and we restrict the discussion to free slip cases for the sake of simplicity.

### 5.1. Detection method

In order to deal with the very large number of results, we have sought for an automatic procedure to detect downwellings and to determine their lateral dimensions. The number of downwellings and their average characteristics vary slightly with time in time-dependent regimes, suggesting that our finite computational domains do not allow a comprehensive statistical sampling of all the possible flow configurations. Steady average properties could be obtained by enlarging the computational domain, which would thus host a larger number of downwellings. An alternative strategy
<figure><img src="image_118.png" /><figcaption>FIGURE 9. Distribution of the dimensionless time-averaged temperature fluctuation $\theta/T_{\text{avg}} = (T - T_{\text{avg}})/T_{\text{avg}}$ in the fluid layer at mid-depth ($z=0.5$) for (a) moderate values and (b) large values of the Rayleigh–Roberts number ($Ra_H$) with free slip boundaries at the top and bottom. Distributions for time-dependent regimes are obtained by stacking a large number of instantaneous ones.</figcaption></figure>

is to average results obtained at different times. Because our calculations are time-dependent, the latter is easier to implement and was therefore adopted.

Several methods have been used to detect downwellings and upwellings in convective systems. Most of these have relied on the temperature field, including the definition of a threshold value that separates an individual plume from the background (Labrosse 2002; Zhong 2005; Galsa & Lenkey 2007; Zhou & Xia 2010). Others have used the convective heat flux, the local dissipation rate or the conditional average of the velocity on the temperature (Ching et al. 2004; Shishkina & Wagner 2008). The latter, more elaborate, methods were devised to identify coherent structures in the complex temperature fields of turbulent Rayleigh–Bénard convection. In our study, flow is in a laminar regime and downwellings clearly stand out of the background temperature field so that the former, and simpler, method was sufficient.

We expect that downwellings, which are driven by a downward buoyancy force, are associated with negative temperature anomalies. In time-dependent regimes, however, one observes tiny areas with very small negative anomalies that are not part of coherent structures and that are highly transient. In order to exclude these small anomalies, we introduce a threshold temperature $T_{\text{ref}}$, related to a ‘detection parameter’ $p_c$,

$$ \frac{T_{\text{ref}} - T_{\text{avg}}}{T_{\text{avg}}} = -p_c \qquad (5.1) $$

so that fluctuations at any depth are properly scaled to the local temperature field. Obviously, the temperature threshold decreases with increasing $p_c$. As described above, at some distance of the unstable boundary layer, downwellings take the shapes of laminar plumes with nearly cylindrical outlines, where temperature decreases to a minimum value at the downwelling centre. Viewed in a horizontal cross-section, a plume thus appears as an area where temperature varies radially in a systematic
<figure><img src="image_119.png" /><figcaption>FIGURE 10. Number of structures with negative thermal anomalies (N) detected as a function of detection parameter $p_c$ for a horizontal plane located at mid-depth ($z=0.5$) in a layer with free slip boundaries. Downwellings are identified using the algorithm described in the text. In this figure, N is not scaled to the horizontal area in order to show the total number of hits and its variation as a function of $p_c$. Results are shown for different values of the Rayleigh–Roberts number $Ra_H$. Note the difference between the steady-state regime at $Ra_H = 5 \times 10^3$ and the two time-dependent ones at larger values of $Ra_H$.</figcaption></figure>

manner. For each threshold temperature value, we have determined contour lines that close upon themselves. With increasing $p_c$, the procedure selects parts of the plumes that are increasingly colder, such that it hits the same plume over increasingly smaller areas. Eventually, of course, the temperature threshold decreases below the smallest downwelling temperature and the procedure returns a count of zero. This method is not efficient for the sheet regime due to the contorted structures and their complex temperature fields, and hence was not used in this regime.

Figure 10 shows the number of structures that are hit as a function of the detection parameter. In the steady-state regime at $Ra_H = 5 \times 10^3$, temperature varies smoothly in the horizontal plane in a well-defined periodic pattern and negative temperature fluctuations are always found in coherent structures. In this case, the number of ‘hits’ is constant for $0 \le p_c \le 0.6$ and drops down sharply for $p_c \gtrsim 0.7$, which corresponds to the coldest temperature. Results are more complicated for time-dependent regimes (figure 10). The number of hits decreases rapidly for small values of the detection parameter and stabilizes to an almost constant value over a large range of $p_c$ values, typically from close to 0 to 0.4–0.5. The initial decrease is due to the tiny areas with small negative anomalies that have been mentioned above, because a large proportion of these small anomalies are detected for very small values of the detection parameter. In contrast, the $p_c$ ‘plateau’ is due to coherent structures that, as explained above, are expected to be detected over a finite $p_c$ interval. The number of hits starts to decrease above a certain $p_c$ value more gradually than in steady-state regimes, due to the more heterogeneous temperature fields.

The performance of the detection procedure is illustrated in figure 11, which shows the temperature field at mid-depth for two time-dependent regimes and for
<figure><img src="image_120.png" /><figcaption>FIGURE 11. Temperature field for rigid (a,b) and free slip (c,d) boundaries in a plane located at mid-depth for Ra<sub>H</sub> = 10<sup>6</sup> (a,c) and 10<sup>9</sup> (b,d). Results for Ra<sub>H</sub> = 10<sup>9</sup> are only shown in part of the horizontal plane, over a dimensionless 2:2 square, in order to facilitate visualization. Note that the dimension of the horizontal plane is different for Ra<sub>H</sub> = 10<sup>6</sup> (6:6) and 10<sup>9</sup> (2:2). Blue and red colours correspond to material that is colder and hotter than the local horizontal average, respectively. The colour scale is not the same in all panels and was chosen to enhance convective structures. Yellow contours illustrate the downwelling outlines that have been determined using the detection procedure described in the text. Values of the detection parameter p<sub>c</sub> were 0.12 and 0.03 for cases with free slip and rigid boundaries, respectively (see text). The numerical resolution and domain aspect ratio are not the same for all these calculations (see table 1) and were changed to deal with the decreasing size of convective cells.</figcaption></figure>

both free slip and rigid boundaries. Downwellings appear as coherent structures with well-defined radial distributions of negative temperature values. They have predominantly circular outlines at Ra<sub>H</sub> = 10<sup>6</sup> and more complex and contorted ones at Ra<sub>H</sub> = 10<sup>9</sup>. One notes that, away from the downwellings, the background temperature field is smooth for Ra<sub>H</sub> = 10<sup>6</sup> and rugged for Ra<sub>H</sub> = 10<sup>9</sup>. In the latter case, the
temperature field is heterogeneous and includes downwellings of different strengths, i.e. characterized by different temperature anomalies. A few downwellings have detached from the upper boundary layer and leave faint ‘ghost’ thermal traces that are picked up by the detection procedure at small $p_c$ values (figure 10). These must not be included in the tally of downwellings, which is easily achieved by raising the detection threshold. A snapshot at a randomly chosen time may not capture a fully representative temperature field and we have repeated the tally several times. Deviations about the mean are small and will be shown.

Our goal is to obtain robust results for all numerical simulations such that they can be compared to one another. We therefore need to follow the same procedure for each calculation in order to avoid potential artefacts due to changing criteria. In other words, we have to consider the same detection parameter $p_c$ for each set of numerical simulations, and to do so we proceed as follows. Steady-state regimes are straightforward with stable counts even for very small values of the detection parameter. For time-dependent regimes, close inspection of the number of hits shows that there is a local minimum at the lower end of the $p_c$ plateau (figure 10). This can be understood from figure 11, where it may be seen that the temperature has more than a single local minimum in a few downwellings with contorted outlines. Such structures would be counted as one at very low $p_c$ values and as two, and sometimes three, at slightly higher $p_c$ values. We conclude that the most appropriate choice is the first local minimum in the number of hits. The uncertainty is very small because the number of counts in the intermediate $p_c$ range does not change by more than two or three units out of a few tens. A constant value of $p_c = 0.12$ is appropriate for free slip boundaries at all $Ra_H$ values. Detection levels are different for cases with rigid boundaries owing to the different downwelling pattern and an appropriate choice is $p_c = 0.03$. Figure 11 shows the downwellings that have been identified by this procedure as well as their outlines (which appear as yellow contours).

## 5.2. Data

With the method that has just been described, we were able to analyse a large number of calculations with minimum effort. We have determined the values of four different variables at mid-depth in the layer: the number of downwellings per unit area, which allows calculation of the average area of a convection cell made of a downwelling and its associated return flow, the average areal extent of the negative thermal anomalies $A_i$, the average thermal anomaly, $\Delta T_i$, and the average vertical velocity, $W_i$. Values are made dimensionless using $d^2$ as a scale for area, $Hd^2/k$ as temperature scale and $\kappa/d$ as velocity scale.

Figure 12 shows results for all calculations save those for the sheet regime in the intermediate $Ra_H$ range. Data for steady-state regimes at $Ra_H \le 5 \times 10^4$ and time-dependent ones at $Ra_H \ge 6.8 \times 10^5$ clearly define two different sets of relationships as a function of $Ra_H$. The most significant difference lies in the number of downwellings per unit area, which remains approximately constant in steady-state cases in contrast to time-dependent ones. Thus, the characteristic spacing between downwellings scales with the total layer thickness for the former and not for the latter. This fundamental change of convection characteristics does not occur as a bifurcation at some threshold $Ra_H$ value but is achieved in the transitional sheet regime over a range of Rayleigh–Roberts numbers. All other variables vary monotonically as $Ra_H$ increases in both steady-state and time-dependent cases. The size of downwellings,
<figure><img src="image_121.png" /></figure>

FIGURE 12. Dependence of (a) the number of downwellings per unit area (Nd²), (b) their average horizontal cross-section (Aᵢ/d²), (c) their average temperature contrast (ΔTᵢk/Hd²), (d) their average vertical velocity (Wᵢd/κ) at mid-depth in the fluid layer as a function of the Rayleigh–Roberts number (RaH). Results are shown for both free slip (red) and rigid (black) boundaries. Downwellings are identified and contoured using an automatic procedure described in the text. Different symbols indicate the convection pattern stable at a given RaH: a hexagonal pattern (square), a spoke pattern (star) and the truncated hexagonal pattern (circle). Error bars, often smaller than the symbol size, indicate variations that are observed in time-dependent regimes. Solid lines are best fits to the data with power laws predicted by the scaling analysis of § 5.3. Values for the various constants in the power-law relationships are listed in tables 6 and 7. Data with open symbols have not been included in the best-fit procedure.

as measured by the areal extent of the thermal anomaly, decreases markedly as RaH increases whilst velocity values increase. Data for the two types of boundary conditions exhibit parallel trends but differ by significant factors. For example, the magnitude of thermal anomalies is 35% larger for rigid boundaries than for free slip ones. The greatest difference is in the number of downwellings in time-dependent cases, which is larger by 60% in rigid cases.
<table><caption>TABLE 6. Parameters of empirical best-fit power laws for the characteristics of downwellings in a planform located at mid-depth for steady-state regimes at low values of the Rayleigh–Roberts number (Ra<sub>H</sub>).</caption><thead><tr><th colspan="2">Ra<sub>H</sub> &le; 5 &times; 10<sup>4</sup></th></tr><tr><th>Property</th><th>Boundary conditions</th></tr></thead><tbody><tr><td rowspan="2">Nd<sup>2</sup></td><td>free slip</td><td>0.433(&plusmn;0.300) Ra<sub>H</sub><sup>-0.120(&plusmn;0.064)</sup></td></tr><tr><td>rigid</td><td>0.0933(&plusmn;0.0650) Ra<sub>H</sub><sup>0.084(&plusmn;0.065)</sup></td></tr><tr><td rowspan="2">A<sub>i</sub>/d<sup>2</sup></td><td>free slip</td><td>152.4(&plusmn;63.73) Ra<sub>H</sub><sup>-0.514(&plusmn;0.040)</sup></td></tr><tr><td>rigid</td><td>34.75(&plusmn;15.43) Ra<sub>H</sub><sup>-0.320(&plusmn;0.043)</sup></td></tr><tr><td rowspan="2">&Delta;T<sub>i</sub>/k/Hd<sup>2</sup></td><td>free slip</td><td>1.117(&plusmn;0.252) Ra<sub>H</sub><sup>-0.132(&plusmn;0.022)</sup></td></tr><tr><td>rigid</td><td>2.218(&plusmn;0.814) Ra<sub>H</sub><sup>-0.167(&plusmn;0.036)</sup></td></tr><tr><td rowspan="2">W<sub>i</sub>/&kappa;</td><td>free slip</td><td>0.00428(&plusmn;0.00095) Ra<sub>H</sub><sup>0.916(&plusmn;0.022)</sup></td></tr><tr><td>rigid</td><td>0.00521(&plusmn;0.00263) Ra<sub>H</sub><sup>0.793(&plusmn;0.049)</sup></td></tr></tbody></table>

In order to facilitate comparisons and to evaluate the sensitivity of the variables to Ra<sub>H</sub>, we have determined power-law relationships using a best-fit procedure, such that:

$$N = C_N(1/d^2)Ra_H^{\beta_N}, \qquad (5.2)$$

$$A_i = C_A d^2 Ra_H^{\beta_A}, \qquad (5.3)$$

$$\Delta T_i = C_T(Hd^2/k)Ra_H^{\beta_T}, \qquad (5.4)$$

$$W_i = C_W(\kappa/d)Ra_H^{\beta_W}, \qquad (5.5)$$

where C<sub>N</sub>, C<sub>A</sub>, C<sub>T</sub> and C<sub>W</sub> are dimensionless proportionality constants and where β<sub>N</sub>, β<sub>A</sub>, β<sub>T</sub> and β<sub>W</sub> are power-law exponents. Since variables do not follow the same trends for Ra<sub>H</sub> &le; 5 &times; 10<sup>4</sup> and for Ra<sub>H</sub> &ge; 6.8 &times; 10<sup>5</sup> (figure 12), we have determined two sets of parameters which are listed in tables 6 and 7. Overall, power-law relationships allow very good fits to the numerical results for steady-state regimes and an excellent fit for time-dependent regimes. The contrasting exponents for steady-state and time-dependent cases confirm the fundamental difference of convection characteristics between the two. Limitations of a best-fit procedure where both the proportionality constant and the exponent are left to vary are well known. Variations of the exponent can be balanced by changes of the proportionality coefficient, which is likely to induce large uncertainties on both parameters. In the next section, we use simple scaling arguments to derive values for the various exponents and compare them to the numerical results.

## 5.3. Scaling analysis

We restrict the analysis to time-dependent regimes which are such that the dynamics of the upper boundary layer is independent of the total fluid thickness (Vilella & Kaminski 2017). In this case, the boundary layer thickness and the temperature contrast across it can be derived from local scaling arguments (table 2). We argue that the interior temperature anomaly ΔT<sub>i</sub> scales with ΔT<sub>TBL</sub> the temperature contrast
<table><caption>TABLE 7. Parameters of best-fit power laws for the characteristics of downwellings in a planform located at mid-depth for time-dependent regimes at large values of the Rayleigh–Roberts number (Ra$_H$). Results in the right-hand column have been obtained by setting exponents to values derived from scaling arguments (see text).</caption><caption>Ra<sub>H</sub> &ge; 6.8 &times; 10<sup>5</sup></caption><thead><tr><th>Property</th><th>Boundary conditions</th><th>Exponent left to vary</th><th>Exponent fixed</th></tr></thead><tbody><tr><td rowspan="2">Nd<sup>2</sup></td><td>free slip</td><td>0.0410(&plusmn;0.0060) Ra<sub>H</sub><sup>0.244(&plusmn;0.009)</sup></td><td>0.0375(&plusmn;0.0008) Ra<sub>H</sub><sup>1/4</sup></td></tr><tr><td>rigid</td><td>0.0226(&plusmn;0.0047) Ra<sub>H</sub><sup>0.309(&plusmn;0.012)</sup></td><td>0.0598(&plusmn;0.0060) Ra<sub>H</sub><sup>1/4</sup></td></tr><tr><td rowspan="2">A<sub>i</sub>/d<sup>2</sup></td><td>free slip</td><td>7.594(&plusmn;1.864) Ra<sub>H</sub><sup>-0.343(&plusmn;0.015)</sup></td><td>12.73(&plusmn;0.763) Ra<sub>H</sub><sup>-3/8</sup></td></tr><tr><td>rigid</td><td>8.772(&plusmn;3.719) Ra<sub>H</sub><sup>-0.328(&plusmn;0.025)</sup></td><td>19.19(&plusmn;1.770) Ra<sub>H</sub><sup>-3/8</sup></td></tr><tr><td rowspan="2">&Delta;T<sub>i</sub>/Hd<sup>2</sup></td><td>free slip</td><td>3.407(&plusmn;0.137) Ra<sub>H</sub><sup>-0.242(&plusmn;0.002)</sup></td><td>3.875(&plusmn;0.053) Ra<sub>H</sub><sup>-1/4</sup></td></tr><tr><td>rigid</td><td>4.604(&plusmn;0.154) Ra<sub>H</sub><sup>-0.231(&plusmn;0.002)</sup></td><td>6.294(&plusmn;0.199) Ra<sub>H</sub><sup>-1/4</sup></td></tr><tr><td rowspan="2">W<sub>i</sub>/&kappa;</td><td>free slip</td><td>1.319(&plusmn;0.073) Ra<sub>H</sub><sup>0.370(&plusmn;0.003)</sup></td><td>1.216(&plusmn;0.013) Ra<sub>H</sub><sup>3/8</sup></td></tr><tr><td>rigid</td><td>1.504(&plusmn;0.116) Ra<sub>H</sub><sup>0.331(&plusmn;0.005)</sup></td><td>0.732(&plusmn;0.053) Ra<sub>H</sub><sup>3/8</sup></td></tr></tbody></table>

across the upper boundary layer, which is such that $\Delta T_{TBL} \sim (Hd^2/k)Ra_H^{-1/4}$ (3.1). This implies that $\beta_T = -1/4$, which is close to the empirically derived values (table 7, figure 12c). Setting $\beta_T = -1/4$, we find that the numerical data for both free slip and rigid boundaries are within a few per cent of the best-fit relationships. Interestingly, the proportionality coefficients in the scaling relationships for $\Delta T_i$ and $\Delta T_{TBL}$ increase by almost the same factor when one moves from free slip boundaries to rigid ones, which confirms that it is indeed the temperature contrast across the boundary layer that sets the amplitude of thermal anomalies in the fluid interior.

In order to derive values for the other power-law exponents, we take into account the structure of the flow field. For $10^6 \le Ra_H \lesssim 10^8$, the flow field may be described as a set of downwellings that go through the whole fluid layer and show up as isolated negative thermal anomalies in an almost uniform background. Such a simple description may not capture all the complex structures that are observed at the largest values of Ra$_H$ and we shall evaluate its relevance a posteriori. We thus consider that each downwelling draws fluid from the upper boundary layer over an area $\Delta S$ and, on average, transports heat at a rate which is equal to

$$Q = q\Delta S = Hd\Delta S, \qquad (5.6)$$

where $q = Hd$ is the heat flux through the top surface. $\Delta S$ is inversely proportional to the number of downwellings per unit area, therefore

$$Q \sim \frac{Hd}{N}. \qquad (5.7)$$

Downwellings are defined by their average thermal contrast with respect to the average background temperature, $\Delta T_i$. They develop circular cross-sections $A_i$ below the upper boundary layer, with radius $R$, such that $R^2 \sim A_i$. The associated velocity scale is

$$W_i \sim \frac{\rho g \alpha \Delta T_i R^2}{\eta}. \qquad (5.8)$$
Conservation of the energy flux implies that

$$Q \sim \rho C_p W_i \Delta T_i R^2, \qquad (5.9)$$

where $C_p$ is heat capacity. Combining (5.8) and (5.9), one obtains

$$W_i \sim \left( \frac{\kappa \rho g \alpha Q}{\eta k} \right)^{1/2}, \qquad (5.10)$$

the Stokes velocity scale for a laminar plume (Batchelor 1954), which is an appropriate scale for the transient plumes generated in time-dependent regimes (Kaminski & Jaupart 2003). We note that this scale depends explicitly on neither the temperature anomaly nor the plume cross-section. Substituting (5.7) in the velocity scale (5.10) together with the Rayleigh–Roberts number (1.1) and power-law expressions for $N$ (5.2) and $W_i$ (5.5), we obtain

$$2\beta_W = 1 - \beta_N. \qquad (5.11)$$

Within uncertainty, the empirical exponents of table 7 are fully consistent with this relationship. For free slip cases, for example, $\beta_W \approx 0.370$ and $\beta_N \approx 0.244$. Similarly, substituting for (5.7) in (5.9), we obtain

$$0 = \beta_N + \beta_A + \beta_T + \beta_W, \qquad (5.12)$$

which is also consistent with the empirically derived exponents and their uncertainties (table 7).

In a laminar regime, the radius of the downwelling is set by a balance between vertical heat advection and horizontal diffusion, such that

$$R \sim \left( \frac{\kappa z}{W_i} \right)^{1/2}, \qquad (5.13)$$

where $z$ is the distance from source, taken here as the distance to the upper boundary. We determine the velocity at mid-depth, i.e. at a fixed distance equal to $d/2$, implying that $R \sim (\kappa d/W_i)^{1/2}$ and hence that $A_i \sim R^2 \sim \kappa d/W_i$. This implies in turn that

$$\beta_A = -\beta_W. \qquad (5.14)$$

This is again consistent with the data in table 7 and their uncertainties. Using (5.12), one deduces further that $\beta_T = -\beta_N$. This is also verified, albeit less satisfactorily for rigid boundaries than for free slip ones. We shall come back to this point at the end of this section.

The above relationships between the four power-law exponents are all consistent with the empirical best-fit values, which supports the validity of the physical arguments. Substituting for the value of $\beta_T$ (-1/4), we find that $\beta_W = 3/8$, $\beta_A = -3/8$ and $\beta_N = 1/4$. These exponents allow excellent fits to the numerical data with few exceptions (figure 12). All the data are typically within less than ±2% of the best-fit power laws for free slip boundaries (except for $A_i$ where the maximum deviation rises to ±6%), but results are slightly less satisfactory for rigid boundaries (table 7). Departures from the proposed scalings may be attributed to two factors. One is due to the assumption that the downwellings issue from point sources, whereas they emanate
from a boundary layer of finite thickness. This problem is more severe for rigid boundary cases due to their thicker boundary layers. Another factor is that the shapes of the downwellings are not as regular in cross-section as we have assumed, as shown by figure 5. This is most pronounced for flows between rigid boundaries at large Ra$_H$ values (figure 11) but, even in these cases, this only degrades significantly the quality of the fits to velocity values and to the number of plumes (table 7). Nevertheless, the power laws for these two variables remain accurate to within ±10% over the whole range of Rayleigh–Roberts number investigated.

### 5.4. Comparison with previous studies

The above results and scaling laws have few equivalents in previous studies. The number of downwellings was determined by Parmentier & Sotin (2000) for free slip boundaries. Their values are consistent with the same Ra$^{1/4}$ scaling law but are somewhat larger than ours (by approximately 30%). This difference may be attributed to the smaller aspect ratio of their computation domain, which was equal to either 2 or 1. We have found indeed that the number of downwellings increases slightly but systematically as the aspect ratio is decreased and argue that this is due to the preferential generation of downwellings at the walls and corners of the domain (supplementary material). In fact, there can be little doubt that the results of Parmentier & Sotin (2000) bear the influence of the vertical bounding walls: for the same Rayleigh–Roberts number ($Ra_H = 10^8$), the value of $Nd^2$ changes when lateral boundary conditions are switched from reflecting to periodic.

At large values of $Ra_H$, the downwellings are generated in a thin upper boundary layer, much as those of Rayleigh–Bénard (RB) convection, and it is worthwhile to compare the two cases. The latter type of convection is controlled by the Rayleigh number, which is defined as

$$Ra = \frac{\rho g \alpha \Delta T^* d^3}{\eta \kappa}, \qquad (5.15)$$

where $\Delta T^*$ is the temperature difference across the layer and where all the other variables are identical to those of this paper. In high-resolution numerical simulations of RB convection, Zhong (2005) found that the number of downwellings increases approximately with $Ra^{1/3}$ over a limited range of Rayleigh numbers and seems to level off for Rayleigh numbers larger than approximately $10^8$. The former behaviour is reminiscent of that of this paper but the latter is not: as shown by figure 12, the number of downwellings increases steadily over the whole $Ra_H$ range investigated. We note that the same behaviour was observed by Parmentier & Sotin (2000) over the same $Ra_H$ range with a completely different numerical method. One might expect that the number of downwellings cannot increase indefinitely due to a ‘crowding’ effect, but, in this aspect, there is a major difference between basal heated and internally heated convection. The former has both active upwellings and downwellings competing for space in contrast to the latter, which only involves active downwellings. Thus, it may well be that crowding limitations in internally heated convective fluids become important at values of the Rayleigh–Roberts number that are larger than those of this paper. One might remark, however, that the fractional area of downwellings scales as $NA_i \sim Ra_H^{-1/8}$ (table 7), which decreases with increasing $Ra_H$.
# 6. Discussion and implications for convection in planetary bodies

## 6.1. General features of internally heated convection

Our numerical simulations illustrate the important effect of boundary conditions on the downwelling characteristics: for a given regime and hence a given Rayleigh-Roberts number, downwellings are significantly colder, narrower and wider apart beneath free slip boundaries than beneath rigid ones. With the four scaling laws that have been derived here and adding those for the thickness and temperature difference across the upper thermal boundary layer (3.1) as well as those for the volume-averaged and bottom temperatures, one can determine all the key characteristics of time-dependent convection from knowledge of the $Ra_H$ value. For geophysical/planetological applications, however, this may not be manageable, partly because calculation of the Rayleigh-Roberts number is fraught with many uncertainties and partly because other factors come into play, such as large variations of physical properties as a function of temperature and pressure and the presence of continents and oceans at the Earth's surface. It may be more useful, therefore, to adopt a different perspective; all the variables of interest here form a self-consistent set, such that knowledge of any one of them leads to estimates for all the others.

In rocky planets with a solid mantle, it is easier to obtain information on the upper boundary layer structure than on the deeper flow characteristics. With respect to a Rayleigh-Bénard set-up, internal heat sources imply a different relationship between downwellings in the fluid interior and the unstable boundary layer at the top. In Rayleigh-Bénard convection, a fixed temperature difference is maintained across the fluid layer, which leads to the formation of boundary layers at both the top and bottom and to upwellings and downwellings of similar strength. We have briefly commented on the different temperature distributions that are generated by the two types of convection. In laminar Rayleigh-Bénard convection at high values of $Ra$ ($\gtrsim 10^5$), the average downwelling velocity $W_{RB}$ and the boundary layer thickness $\delta_{RB}$ scale as $(\kappa/d)Ra^{2/3}$ and $dRa^{-1/3}$, respectively (Turcotte & Oxburgh 1967; Galsa & Lenkey 2007). Note that these relationships may not be valid for $Ra$ values that are larger than about $10^7$ (e.g. Zhong 2005). One deduces that

$$W_{RB} \sim \frac{\kappa}{d} \left( \frac{\delta_{RB}}{d} \right)^{-2} . \qquad (6.1)$$

In an internally heated convecting layer, combining the scaling laws derived above with (3.1) for the boundary layer thickness, we obtain

$$W_i \sim \frac{\kappa}{d} \left( \frac{\delta}{d} \right)^{-3/2} . \qquad (6.2)$$

Alternatively, one could derive relationships between velocity and the temperature difference across the boundary layer because it may be easier to obtain constraints on temperature than on thickness in some cases. The key point, however, is the same: the relationship between downwellings and boundary layer structure is fundamentally different for Rayleigh-Bénard and internally heated systems.

## 6.2. Subduction zones on Earth

The Earth's mantle is clearly more complex than the internally heated fluid layer studied here. Nevertheless, it is worth evaluating our results in this context because
they illustrate the peculiarities of internally heated systems and the range of dynamical regimes that can be achieved. Our point is not to claim that our calculations provide realistic representations of the Earth's mantle but to argue that one should be careful when attributing Earth's characteristic features to a particular mechanism.

The Earth's present-day convective regime is characterized by a small number of 'plates' created at mid-ocean ridges and subducted along trenches that stretch over large horizontal distances. Surprising features are the large range of plate dimensions and velocities, which typically span approximately one order of magnitude, and the complex network of subduction zones which do not conform to any simple geometrical pattern (Mora et al. 2013). Seismic images of the mantle interior show that mid-ocean ridges, which are the surface expressions of local upwellings, are not underlain by deep-seated thermal anomalies, in marked contrast to subduction zones (e.g. Ritsema et al. 2011). Such characteristics are hallmarks of an internally heated convective system. Seismological data also indicate that subducting plates remain whole as they go through the upper mantle. In this regard, the sheet regime is of special interest because of its large-scale irregular network of narrow linear downwellings that extend to the base of the layer. The most challenging fact is probably that Earth seems to be the only rocky planet of the solar system hosting plate tectonics. Planet Venus, which is comparable in size and composition, has a convective mantle, as shown by the many volcanoes that dot its surface, but no plates and no subduction zones (Schubert et al. 2010). Another conundrum deals with when plate tectonics started on Earth, which remains hotly debated today. Some authors have proposed that plate tectonics is essentially as old as the planet and others have argued that its onset coincided with a major compositional change that occurred 3 Gyr ago in the subcontinental lithospheric mantle (Shirey & Richardson 2011). Answers to both questions may be sought in transient evolutionary models, such that plate tectonics is only active for a finite length of time. As a planet cools down and its heat sources get depleted, the Rayleigh-Roberts decreases and the convection regime and planform are bound to change. The sheet regime serves to emphasize the dramatic planform transformations that can occur.

Subducting slabs can be linked to the downgoing sheets that have been documented here. Subduction zones occasionally form triple junctions, such as in the Boso subduction margin, central Japan (Yujiro et al. 1989), which are strongly reminiscent of the planforms of this study. Our calculations show that sheets are not the only form of convection that can occur and that they require values of the Rayleigh-Roberts number above a certain threshold. Planets with small Ra<sub>H</sub> values, due either to their small sizes or to their small rates of internal heat production, are not expected to be in the same regime and would likely involve hexagons with an axial upwelling. By the same token, the fact that the downgoing sheets of Earth remain whole over the total vertical extent of the upper mantle provides another constraint: sheets break up into nearly cylindrical downgoing plumes if the Rayleigh-Roberts number exceeds a critical value. Thus, plate tectonics, which is the current convection regime on Earth, may not be the standard regime for all silicate planets and, furthermore, may not have prevailed early in Earth's history, as suggested by Jellinek & Jackson (2015) amongst others. One should also note that the sheet regime exists over a much larger range of Ra<sub>H</sub> values for rigid boundaries than for free ones, which illustrates the sensitivity of the superficial convection planform and the bulk flow characteristics to the mechanical boundary condition at the top.
## 7. Conclusion

Three-dimensional numerical simulations of laminar convection in an internally heated fluid layer cooled from above, illustrate the different superficial planforms that may be generated. As the Rayleigh-Roberts number increases, we have found four different planforms described as hexagons with an axial cylindrical downwelling, hexagons with a spoke-like axial downwelling, networks of linear downwellings or sheets and finally truncated hexagons. Rigid and free slip boundaries allow the same planforms but are associated with marked differences in the spacing and other characteristics of convective motions. In time-dependent regimes at large Ra<sub>H</sub> values (≥10<sup>6</sup>), complex superficial planforms morph into arrays of nearly equally spaced cylindrical plume-like downwellings with characteristics that depend on the Ra<sub>H</sub> value in systematic ways. Temperature distributions in both the horizontal and vertical directions are markedly different from those of Rayleigh-Bénard convection. Scaling laws for the main variables of interest, including the thermal anomalies, velocities, dimensions and spacing of descending plumes, have been developed and are fully consistent with numerical results over a Ra<sub>H</sub> range spanning three orders of magnitude. They illustrate important differences of boundary layer dynamics with respect to Rayleigh-Bénard convection.

The sheet regime deserves more attention because of its peculiar characteristics and also because it illustrates that dramatic reorganizations of flow structure in rocky planets can be caused by small changes of internal heat release and physical properties. Much is likely to be gained by investigating how such reorganizations proceed and how sensitive they are to temperature-dependent fluid properties.

## Acknowledgements

We thank three anonymous reviewers for their thorough and constructive comments that helped us to produce a much improved paper. We also thank J. Wettlaufer for his editorial handling of the manuscript. This work was funded by the ANR-11-IS04-0004 project and supported by the Tellus/PNP program of INSU-CNRS. Numerical computations were performed on the S-CAPAD platform, IPGP, France and using HPC resources from GENCI-IDRIS (Grant 2013-047033).

## Supplementary material

Supplementary material is available at https://doi.org/10.1017/jfm.2018.316.

## REFERENCES

* ADRIAN, R. J., FERREIRA, R. T. D. S. & BOBERG, T. 1986 Turbulent thermal convection in wide horizontal fluid layers. *Exp. Fluids* **4** (3), 121–141.

* BATCHELOR, G. K. 1954 Heat convection and buoyancy effects in fluids. *Q. J. R. Meteorol. Soc.* **80**, 339–358.

* BELLO, L., COLTICE, N., ROLF, T. & TACKLEY, P. J. 2014 On the predictability limit of convection models of the Earth’s mantle. *Geochem. Geophys. Geosyst.* **15**, 2319–2328.

* BOUSSINESQ, J. 1903 *Théorie analytique de la chaleur mise en harmonie avec la thermodynamique et avec la théorie mécanique de la lumière, Tome II.* Gauthier-Villars.

* CARRIGAN, C. R. 1982 Multiple-scale convection in the Earth’s mantle: a three-dimensional study. *Science* **215**, 965–967.

* CARRIGAN, C. R. 1985 Convection in an internally heated, high Prandtl number fluid: a laboratory study. *Geophys. Astrophys. Fluid Dyn.* **32**, 1–21.
CASTAING, B., GUNARATNE, G., HESLOT, F., KADANOFF, L., LIBCHABER, A., THOMAE, S., WU, X. Z., ZALESKI, S. & ZANETTI, G. 1989 Scaling of hard thermal turbulence in Rayleigh-Bénard convection. *J. Fluid Mech.* **204**, 1–30.

CHEUNG, F. B. 1977 Natural convection in a volumetrically heated fluid layer at high Rayleigh numbers. *Intl J. Heat Mass Transfer* **20**, 499–506.

CHING, E. S. C., GUO, H., SHANG, X. D., TONG, P. & XIA, K. Q. 2004 Extraction of plumes in turbulent thermal convection. *Phys. Rev. Lett.* **93** (12), 124501.

CHRISTENSEN, U. R. & HOFMANN, A. W. 1994 Segregation of subducted oceanic crust in the convecting mantle. *J. Geophys. Res.* **99**, 19867–19884.

COTTAAR, S. & BUFFETT, B. 2012 Convection in the Earth's inner core. *Phys. Earth Planet. Inter.* **198**, 67–78.

CROSS, M. C. & HOHENBERG, P. C. 1993 Pattern formation outside of equilibrium. *Rev. Mod. Phys.* **65** (3), 851–1112.

DAVAILLE, A. & LIMARE, A. 2015 Laboratory studies of mantle convection. In *Treatise on Geophysics*, vol. 7, pp. 73–144. Elsevier.

DESCHAMPS, F., ROGISTER, Y. & TACKLEY, P. J. 2018 Constraints on core-mantle boundary topography from models of thermal and thermo-chemical convection. *Geophys. J. Intl* **212**, 164–188.

DESCHAMPS, F. & TACKLEY, P. J. 2008 Searching for models of thermo-chemical convection that explain probabilistic tomography. Part I. Principles and influence of rheological parameters. *Phys. Earth Planet. Inter.* **171**, 357–373.

DESCHAMPS, F. & TACKLEY, P. J. 2009 Searching for models of thermo-chemical convection that explain probabilistic tomography. Part II. Influence of physical and compositional parameters. *Phys. Earth Planet. Inter.* **176**, 1–18.

GALSA, A. & LENKEY, L. 2007 Quantitative investigation of physical properties of mantle plumes in three-dimensional numerical models. *Phys. Fluids* **19**, 116601.

GLOVER, G. M. C. & GENERALIS, S. C. 2009 Pattern competition in homogeneously heated fluid layers. *Engng Appl. Comput. Fluid Mech.* **3**, 164–174.

GOLUSKIN, D. 2015 *Internally Heated Convection and Rayleigh-Bénard Convection*. Springer International Publishing.

GOLUSKIN, D. & SPIEGEL, E. A. 2012 Convection driven by internal heating. *Phys. Lett. A* **377**, 83–92.

GOMI, H., OHTA, K., HIROSE, K., LABROSSE, S., CARACAS, R., VERSTRAETE, M. J. & HERNLUND, J. W. 2013 The high conductivity of iron and thermal evolution of the Earth's core. *Phys. Earth Planet. Inter.* **224**, 88–103.

GUILLOU, L. & JAUPART, C. 1995 On the effect of continents on mantle convection. *J. Geophys. Res.* **100** (B12), 24217–24238.

GURNIS, M. 1988 Large-scale mantle convection and the aggregation and dispersal of supercontinents. *Nature* **332** (6166), 695–699.

HOUSEMAN, G. 1988 The dependence of convection planform on mode of heating. *Nature* **332**, 346–349.

HUSSMANN, H., CHOBLET, G., LAINEY, V., MATSON, D. L., SOTIN, C., TOBIE, G. & VAN HOOLST, T. 2010 Implications of rotation, orbital states, energy sources, and heat transport for internal processes in icy satellites. *Space Sci. Rev.* **153**, 314–348.

ICHIKAWA, H., KURITA, K., YAMAGISHI, Y. & YANAGISAWA, T. 2006 Cell pattern of thermal convection induced by internal heating. *Phys. Fluids* **18**, 038101.

JAUPART, C., LABROSSE, S., LUCAZEAU, F. & MARESCHAL, J. C. 2015 Temperatures, heat and energy in the mantle of the Earth. In *Treatise on Geophysics*, vol. 7, pp. 223–270. Elsevier.

JELLINEK, A. M. & JACKSON, M. G. 2015 Connections between the bulk composition, geodynamics and habitability of Earth. *Nat. Geosci.* **8**, 587–593.

KAMINSKI, E. & JAUPART, C. 2003 Laminar starting plumes in high-Prandtl-number fluids. *J. Fluid Mech.* **478**, 287–298.
KULACKI, F. A. & EMARA, A. A. 1977 Steady and transient thermal convection in a fluid layer with uniform volumetric energy sources. *J. Fluid Mech.* **83**, 375–395.

KULACKI, F. A. & GOLDSTEIN, R. J. 1975 Hydrodynamic instability in fluid layers with uniform volumetric energy sources. *Appl. Sci. Res.* **31**, 81–109.

KULACKI, F. A. & NAGLE, M. E. 1975 Natural convection in a horizontal fluid layer with volumetric energy sources. *Trans. ASME J. Heat Transfer* **97**, 204–211.

LABROSSE, S. 2002 Hotspots, mantle plumes and core heat loss. *Earth Planet. Sci. Lett.* **199**, 147–156.

LABROSSE, S., POIRIER, J. P. & LE MOUËL, J. L. 2001 The age of the inner core. *Earth Planet. Sci. Lett.* **190**, 111–123.

LAY, T., HERNLUND, J., GARNERO, E. J. & THORNE, M. S. 2006 A post-Perovskite lens and D'' heat flux beneath the central Pacific. *Science* **314** (5803), 1272–1276.

LIMARE, A., VILELLA, K., DI GIUSEPPE, E., FARNETANI, C., KAMINSKI, E., SURDUCAN, E., SURDUCAN, V., NEAMTU, C., FOUREL, L. & JAUPART, C. 2015 Microwave-heating laboratory experiments for planetary mantle convection. *J. Fluid Mech.* **1565**, 14–18.

MCKENZIE, D. P., ROBERTS, J. M. & WEISS, N. O. 1973 Numerical models of convection in the Earth's mantle. *Tectonophysics* **19**, 89–103.

MITROVICA, J. X. & FORTE, A. M. 2004 A new inference of mantle viscosity based upon joint inversion of convection and glacial isostatic adjustment data. *Earth Planet. Sci. Lett.* **225**, 177–189.

MOORE, W. B. 2008 Heat transport in a convecting layer heated from within and below. *J. Geophys. Res.* **113**, B11407.

MORA, G., SETON, M., QUEVEDO, L. & MÜLLER, R. D. 2013 Organization of the tectonic plates in the last 200 Myr. *Earth Planet. Sci. Lett.* **373**, 93–101.

NAKAGAWA, T. & TACKLEY, P. J. 2014 Influence of combined primordial layering and recycled MORB on the coupled thermal evolution of Earth's mantle and core. *Geochem. Geophys. Geosyst.* **15**, 619–633.

PARMENTIER, E. M. & SOTIN, C. 2000 Three-dimensional numerical experiments on thermal convection in a very viscous fluid: implications for the dynamics of a thermal boundary layer at high Rayleigh number. *Phys. Fluids* **12**, 609–617.

RITSEMA, J., DEUSS, A., VAN HEIJST, H. J. & WOODHOUSE, J. H. 2011 S40rts: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements. *Geophys. J. Intl* **184** (3), 1223–1236.

ROBERTS, P. H. 1967 Convection in horizontal layers with internal heat generation: theory. *J. Fluid Mech.* **30**, 33–49.

ROLF, T., COLTICE, N. & TACKLEY, P. J. 2012 Linking continental drift, plate tectonics and the thermal state of the Earth's mantle. *Earth Planet. Sci. Lett.* **351**, 134–146.

SAMUEL, H. & FARNETANI, C. G. 2003 Thermochemical convection and helium concentrations in mantle plumes. *Earth Planet. Sci. Lett.* **207**, 39–56.

SCHUBERT, G., HUSSMANN, H., LAINEY, V., MATSON, D. L., MCKINNON, W. B., SOHL, F., SOTIN, C., TOBIE, G., TURRINI, D. & VAN HOOLST, T. 2010 Evolution of icy satellites. *Space Sci. Rev.* **153**, 447–484.

SCHUBERT, G., TURCOTTE, D. L. & OLSON, P. 2001 *Mantle Convection in the Earth and Planets*. Cambridge University Press.

SEIS, C. 2013 Laminar boundary layers in convective heat transport. *Commun. Math. Phys.* **324**, 995–1031.

SHIREY, S. B. & RICHARDSON, S. H. 2011 Start of the Wilson cycle at 3 Ga shown by diamonds from the subcontinental mantle. *Science* **333**, 434–436.

SHISHKINA, O. & WAGNER, C. 2008 Analysis of sheet-like thermal plumes in turbulent Rayleigh-Bénard convection. *J. Fluid Mech.* **599**, 383–404.

SMREKAR, S. E. & SOTIN, C. 2012 Constraints on mantle plumes on Venus: implications for volatile history. *Icarus* **217**, 510–523.
STIXRUDE, L. & LITHGOW-BERTELLONI, C. 2011 Thermodynamics of mantle minerals? Part II. Phase equilibria. *Geophys. J.* **184**, 1180–1213.

TACKLEY, P. J. 1998a Self-consistent generation of tectonic plates in three-dimensional mantle convection. *Earth Planet. Sci. Lett.* **157**, 9–22.

TACKLEY, P. J. 1998b Three-dimensional simulations of mantle convection with a thermochemical CMB boundary layer: *D''*? In *The Core–Mantle Boundary Region* (ed. M. Gurnis, M. E. Wysession, E. Knittle & B. A. Buffett), pp. 231–253. AGU.

TACKLEY, P. J. 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**, 7–18.

TAKAHASHI, J., TASAKA, Y., MURAI, Y., TAKEDA, Y. & YANAGISAWA, T. 2010 Experimental study of cell pattern formation induced by internal heat sources in a horizontal fluid layer. *Intl J. Heat Mass Transfer* **53**, 1483–1490.

THEERTHAN, S. A. & ARAKERI, J. H. 1994 Planform structure of turbulent Rayleigh–Bénard convection. *Intl Commun. Heat Mass Transfer* **21** (4), 561–572.

TRITTON, D. J. & ZARRAGA, M. N. 1967 Convection in horizontal layers with internal heat generation. Experiments. *J. Fluid Mech.* **30**, 21–31.

TURCOTTE, D. L. & OXBURGH, E. R. 1967 Finite amplitude convective cells and continental drift. *J. Fluid Mech.* **28** (1), 29–42.

VILELLA, K. & KAMINSKI, E. 2017 Fully determined scaling laws for volumetrically heated convective systems, a tool for assessing habitability of exoplanets. *Phys. Earth Planet. Inter.* **266**, 18–28.

YUJIRO, O., TETSUZO, S., HINAKO, A., HIDEKAZU, T., KANTARO, F. & HIDETSUGU, T. 1989 Structure and development of the Sagami trough and the Boso triple junction. *Tectonophysics* **160**, 135–150.

ZHONG, S. 2005 Dynamics of thermal plumes in three-dimensional isoviscous thermal convection. *Geophys. J. Intl* **162**, 289–300.

ZHOU, Q. & XIA, K. Q. 2010 Physical and geometrical properties of thermal plumes in turbulent Rayleigh–Bénard convection. *New J. Phys.* **12** (7), 075006.
# Convection in an internally heated stratified heterogeneous reservoir

Angela Limare<fnref n="1" />,†, Claude Jaupart<fnref n="1" />, Edouard Kaminski<fnref n="1" />, Loic Fourel<fnref n="2" /> and Cinzia G. Farnetani<fnref n="1" />

(Received 16 November 2018; revised 22 February 2019; accepted 22 March 2019)

The Earth's mantle is chemically heterogeneous and probably includes primordial material that has not been affected by melting and attendant depletion of heat-producing radioactive elements. One consequence is that mantle internal heat sources are not distributed uniformly. Convection induces mixing, such that the flow pattern, the heat source distribution and the thermal structure are continuously evolving. These phenomena are studied in the laboratory using a novel microwave-based experimental set-up for convection in internally heated systems. We follow the development of convection and mixing in an initially stratified fluid made of two layers with different physical properties and heat source concentrations lying above an adiabatic base. For relevance to the Earth's mantle, the upper layer is thicker and depleted in heat sources compared to the lower one. The thermal structure tends towards that of a homogeneous fluid with a well-defined time constant that scales with $Ra_H^{-1/4}$, where $Ra_H$ is the Rayleigh–Roberts number for the homogenized fluid. We identified two convection regimes. In the dome regime, large domes of lower fluid protrude into the upper layer and remain stable for long time intervals. In the stratified regime, cusp-like upwellings develop at the edges of large basins in the lower layer. Due to mixing, the volume of lower fluid decreases to zero over a finite time. Empirical scaling laws for the duration of mixing and for the peak temperature difference between the two fluids are derived and allow extrapolation to planetary mantles.

**Key words:** mantle convection, mixing and dispersion

## 1. Introduction

Convection in the Earth's mantle is largely driven by internal energy sources, involving heat released by the radioactive decay of long-lived isotopes of uranium, thorium and potassium and sensible heat extracted through secular cooling (Schubert, Turcotte & Olson 2001; Jaupart et al. 2015). The Earth's mantle is also heated from below by the Earth's core but the magnitude of the basal heat flux is not known precisely. Today's mantle motions are organized in two radically different planforms and spatial scales. The dominant planform is a set of thin subduction zones stretching

<footnote n="1">¹Institut de Physique du Globe, Université de Paris, CNRS, F-75005 Paris, France</footnote>

<footnote n="2">²Solid Earth Geology Team, Geological Survey of Norway (NGU), Trondheim NO-7491, Norway</footnote>

<footnote>† Email address for correspondence: limare@ipgp.fr</footnote>
over large distances at the edges of oceanic plates, where cold lithosphere goes down. These downwellings are associated with positive seismic velocity anomalies that can be traced to great depths (van der Hilst & Kárason 1999). The return flow proceeds through mid-ocean ridges, but these are not underlain by deep-seated seismic anomalies, indicating that flow is essentially of a passive nature (Ritsema et al. 2011). These are hallmarks of convection in an internally heated system that is cooled from above, where motions are driven by density variations generated in an unstable thermal boundary layer at the top (Roberts 1967; Kulacki & Goldstein 1972; Parmentier, Sotin & Travis 1994; Goluskin 2016; Vilella et al. 2018). In the Earth, convection also involves a number of narrow upwellings feeding intraplate volcanoes called ‘hotspots’, such as beneath Hawaii and Reunion islands. Many prominent hotspots stand above nearly vertical negative seismic velocity anomalies that extend down to the top of the core (Montelli et al. 2006; French & Romanowicz 2015). These observations appear to support the traditional mantle plume model that calls for basal heating by the Earth’s core. These plumes are not distributed at random, however, and are linked to large anomalous structures in the lowermost mantle. High-resolution seismic studies have revealed that the ≈200–300 km thick so-called D'' basal layer is in fact not laterally continuous. It is more appropriate to refer to a basal region of several hundred kilometres thickness where material is anomalous and much more heterogeneous than the rest of the lower mantle (Garnero & Helmberger 1996; Garnero, McNamara & Shim 2016). Seismic wave speeds are anomalously low in two broad regions called large low-shear-velocity provinces (LLSVPs), which extend over a thickness of as much as 1000 km above the core–mantle boundary and which are difficult to reconcile with convection in a homogeneous mantle. Of particular interest is that hotspot seismic anomalies preferentially lie at the edges of these LLSVPs (Torsvik et al. 2006; Deschamps, Cobden & Tackley 2012; Garnero et al. 2016; Li et al. 2017).

Geochemical data confirm that the Earth's mantle is heterogeneous and indicate the presence of different materials, primordial mantle, mantle that has been depleted through melting and possibly a second type of primordial mantle inherited from early planetary formation and differentiation processes (Hofmann 2003; Gale et al. 2013). These contain different amounts of radioactive elements, implying that heat sources are not uniformly distributed (Turcotte, Paul & White 2001; Javoy & Kaminski 2014). One type of material is the consequence of melt extraction at shallow depths beneath ocean ridges and trenches, which generates a depleted oceanic lithosphere and enriched continental crust, processes which have been active over most of Earth's history. A second type may be inherited from early planetary processes (Javoy & Kaminski 2014) and may have survived for several billions of years. How it was generated remains unclear. Contributing mechanisms include the settling of iron phases and crystallization in an early magma ocean phase. Iron-silicate and melt-solid equilibria depend on pressure, implying that the primordial mantle composition was likely to depend on depth. The locations of the different mantle materials remain debated. What is known with certainty is that continental crust has been extracted from the mantle, leaving a residue that has been qualified as 'depleted'. Because continental crust is a concentrate of heat-producing elements, 'depletion' adequately reflects the impact of crust formation on the mantle with respect to internal heating (Turcotte et al. 2001). The upper mantle is made of different materials mixed in variable proportions, including primordial mantle brought by deep plumes (Gale et al. 2013). Mass balance constraints are met with a three reservoir structure, made of a primordial basal region, a depleted mid-mantle and an only slightly depleted upper
mantle (Jackson & Carlson 2012; Gale et al. 2013). The difference between the upper and mid-mantle reservoirs is due to subduction, which goes through the former and continuously injects newly depleted material into the latter.

The presence of undepleted primordial material with higher than average heat production at the bottom of the mantle makes for a deep supply of heat able to generate upwelling activity. Thus, mantle plumes may be due to heat coming from the core as well as heat produced in a basal reservoir. Discriminating between these two contributions is important because they involve different physical controls and are not likely to follow the same time evolution. There is little doubt that the Earth's core cannot host significant amounts of radioactive elements, implying that it can only heat the mantle if it is cooling down. The heat flux out of the core depends on thermal coupling with the highly viscous mantle, and has probably changed with time. In an initially thermally well-mixed planet, for example, this heat flux would have been zero. In similar fashion, the amount of primordial material at the base of the mantle cannot be taken as constant because convective motions are bound to induce mixing with the overlying material.

Starting from two layers of different materials, mixing may proceed by two different processes, the folding of one fluid over the other and the tearing out of thin schlieren at the interface (Olson & Kincaid 1991). The latter process operates at very small scales and presents a difficult challenge for direct numerical simulations (Deschamps & Tackley 2008, 2009). Quantitative laboratory studies have been carried out in Rayleigh-Bénard set-ups with fixed temperatures at the top and bottom or a fixed heat flux (Richter & McKenzie 1981; Olson 1984; Olson & Kincaid 1991; Davaille 1999a,b). These experiments reveal a wealth of phenomena such as the oscillatory motions of plumes initiated in a dense lower layer and the overturn of the layers. They also show that protrusions of the lower layer into the upper one act to anchor upwellings at the same locations for long time intervals, with important implications for mantle convection (Davaille, Girard & Le Bars 2002; Jellinek & Manga 2002). A recent study by Lepot, Aumaître & Gallet (2018) sheds light on heat transport in a fluid layer where internal heat sources are concentrated in a basal region. To the best of our knowledge, however, no experiments are available for internally heated compositionally stratified reservoirs, which leaves an important gap in studies of mantle convection.

In the Earth's mantle, the distribution of heat sources affects the pattern and amplitude of convective motions and in turn gets modified by the flow. As explained above, one cannot consider that either the heat flux or temperature at the core-mantle boundary have remained constant through time. It is thus worthwhile to focus on a stand alone system powered by its heat sources with no heat supplied by a separate system. In order to study mixing phenomena over long time intervals quantitatively, we rely on laboratory experiments. We investigate the behaviour of a stratified reservoir that is cooled from above and that has a larger concentration of heat sources in a lower layer above an adiabatic boundary. Due to mixing, the internal structure becomes increasingly removed from the starting one. We track the time-dependent spatial pattern and temperature differences that drive convection and investigate how they depend on the control variables. In the Earth's mantle, complete mixing, such that discriminating between individual material components is no longer possible at the smallest scale of relevance, involves solid-state diffusion and is unlikely due to the very small values of the diffusion coefficient, as indicated by geochemical data (Hofmann 2003; Gale et al. 2013). We determine the time for the pervasive mingling of the two fluids, such that the initial two-layer configuration has been
eradicated completely, and work out how and when the reservoir starts behaving as a homogeneous one.

The paper is organized as follows. We describe the experimental set-up and protocol as well as the relevant control variables and dimensionless numbers. We determine the bulk thermal evolution and surface heat flux and show that, after a small lead time, both depend weakly on the distribution of internal heat sources. We identify and characterize two convection regimes based on the statistical distribution of the interface depth. These two regimes depend on a buoyancy number which scales density differences due to temperature to the intrinsic density contrast between the two fluids and on the intensity of convective motions as measured by a Rayleigh number. We describe in detail the changes of thermal structure that are induced by mixing and derive empirical scaling laws for the temperature excess and the lifetime of the enriched lower layer. We discuss some implications of our results for the Earth's mantle in a final section.

## 2. Dynamical regimes and governing parameters

### 2.1. Dimensionless numbers for homogeneous internally heated convection

In a homogeneous internally heated fluid layer that is cooled from above and that has an adiabatic base, the relevant temperature scale is:

$$ \Delta T_H = \frac{Hh^2}{\lambda}, \qquad (2.1) $$

where $h$ is the reservoir thickness, $\lambda$ is thermal conductivity and $H$ is the rate of heat generation per unit volume. Using standard scales for time, velocity and length, the governing Boussinesq equations lead to two dimensionless numbers, the Rayleigh–Roberts number $Ra_H$ and the Prandtl number $Pr$ (Roberts 1967):

$$ Ra_H = \frac{\rho g \alpha H h^5}{\lambda \kappa \mu}, \qquad (2.2) $$

and

$$ Pr = \frac{\nu}{\kappa}, \qquad (2.3) $$

where $\rho$ is the density at some reference temperature, $g$ is the acceleration of gravity, $\alpha$ is the thermal expansion coefficient, $\kappa$ is the thermal diffusivity, $\mu$ is the dynamic viscosity and $\nu = \mu/\rho$ is the kinematic viscosity. For sufficiently large values of $Pr$, inertial effects are negligible compared to viscous effects and one can work in the infinite Prandtl limit. In these conditions, which are certainly those of telluric planets with $Pr > 10^{23}$, the characteristics of convection depend only on $Ra_H$.

### 2.2. Dimensionless numbers for a two-layer system

We investigate convection in an initially two-layer reservoir involving two miscible fluids with different physical properties and thicknesses (figure 1). The upper surface is kept at a constant temperature and the lower one is adiabatic. The two fluids have identical coefficients of thermal expansion, thermal conductivities and heat capacities, as in the Earth's mantle, where the small variations of composition that exist have no significant impact on these properties. In contrast, viscosity values and rates of
<figure><img src="image_122.png" /><figcaption>FIGURE 1. (Colour online) Conduction temperature profile (a) and the corresponding density profile (b) in a two-layer reservoir. (a) Temperature elevation in layer 2 due to internal heating only (light blue), temperature elevation in layer 2 due to the heat flux at its base (green), total temperature profile in layer 2 (red) and temperature profile in layer 1 due to internal heating only (blue). (b) Density profile in layer 2 (red) and in layer 1 (blue).</figcaption></figure>

internal heat generation in mantle material are very sensitive to trace amounts of water and radioactive elements, respectively. The rheological properties of mantle rocks also depend on pressure, stress, grain size and mineral composition, but these effects will be ignored here for simplicity. We therefore consider two fluids which differ by their intrinsic densities, viscosities and heat generation rates. Our working fluids are dilute aqueous solutions of sodium chloride salt and hydroxyethylcellulose, a compound that induces large viscosity changes for very small concentrations.

In the following, indices 1 and 2 will refer to the lower and upper layers, respectively. Equations of state are as follows:

$$ \rho_i = \rho_{io}[1 - \alpha(T - T_o)], \qquad (2.4) $$

where $i = 1, 2$ refers to each fluid. Reference temperature $T_o$ is that of the upper surface and can be taken as equal to zero. The two fluids occupy a total thickness $h$ initially split into two layers with thicknesses $h_1$ and $h_2$, such that $h = h_1 + h_2$. Thickness ratio $a = h_1/h$ is equal to the volume fraction of the lower fluid and remains a key control variable even when the layered structure has been destroyed. With heat production rates $H_1$ and $H_2$ in the lower and upper fluids, respectively, the total rate of heat released in the reservoir is:

$$ Hh = H_1h_1 + H_2h_2. \qquad (2.5) $$

We define an enrichment factor for the lower fluid as follows:

$$ F = \frac{H_1}{H}, \qquad (2.6) $$

and a viscosity ratio:

$$ \gamma = \frac{\mu_1}{\mu_2}. \qquad (2.7) $$
One last dimensionless number is needed to characterize the two-layer system, which is the buoyancy number (i.e. the ratio of the stabilizing density anomaly to the destabilizing thermal anomaly):

$$B = \frac{\rho_{10} - \rho_{20}}{\rho \alpha \Delta T} = \frac{\Delta \rho}{\rho \alpha \Delta T}, \qquad (2.8)$$

where $\rho_{10} - \rho_{20} = \Delta\rho$ is the intrinsic density difference, that is the density difference due to composition only, as opposed to the actual density difference which depends on composition and temperature, and $\Delta T$ is an appropriate scale for the temperature contrast between the two fluids. In a major difference with Rayleigh-Bénard experiments, this temperature contrast is not fixed and changes with time as mixing progresses. Temperature scales for the two layers can be determined using a conduction equilibrium reference state (figure 1a). The appropriate boundary conditions are zero heat flux at the reservoir base (at $z=0$) and zero temperature at the top (at $z=h$). Integrating the diffusion heat equation and applying the continuity of temperature and heat flux at the interface (at $z=h_1$), one obtains:

$$\Delta T_1(z) = H_1 h_1 \frac{h_2}{\lambda} + \frac{H_2 h_2^2}{2\lambda} + \frac{H_1 h_1^2}{2\lambda} \left(1 - \frac{z^2}{h_1^2}\right), \qquad (2.9)$$

$$\Delta T_2(z) = H_1 h_1 \frac{h-z}{\lambda} + \frac{H_2 h_2^2}{2\lambda} \left(1 - \frac{(z-h_1)^2}{h_2^2}\right). \qquad (2.10)$$

Figure 1(b) shows the corresponding density distribution. The temperature differences across the two layers are, respectively:

$$\Delta T_2 = \frac{H_2 h_2^2}{2\lambda} + \frac{H_1 h_1 h_2}{2\lambda} = \frac{1}{2} \Delta T_H (1-a)(1+aF), \qquad (2.11)$$

$$\Delta T_1 = \frac{H_1 h_1^2}{2\lambda} = \frac{1}{2} \Delta T_H [1 + a(F - 1)], \qquad (2.12)$$

where $\Delta T_H$ is the temperature scale for a homogeneous layer of thickness $h$ and heat generation $H$ (2.1). The average temperature in each layer is obtained by integrating (2.9) and (2.10):

$$\overline{\Delta T_1} = \frac{1}{2} \frac{H_2 h_2^2}{\lambda} + \frac{H_1 h_1 h_2}{\lambda} + \frac{1}{3} \frac{H_1 h_1^2}{\lambda}, \qquad (2.13)$$

$$\overline{\Delta T_2} = \frac{1}{3} \frac{H_2 h_2^2}{\lambda} + \frac{1}{2} \frac{H_1 h_1 h_2}{\lambda}. \qquad (2.14)$$

The global temperature difference between the two layers is thus:

$$\Delta T = \overline{\Delta T_1} - \overline{\Delta T_2} = \frac{1}{6} \Delta T_H (1 - a + 2aF). \qquad (2.15)$$

This temperature difference corresponds to a reference thermal structure and illustrates the influence of several dimensionless numbers. It is likely to be very different from that of the fully convecting system, however, and we shall also consider an 'effective' temperature difference determined in the course of an experiment. In the following, we shall use two different temperature scales. To deal with the coupling between the two fluids, we use both the conduction thermal contrast $\Delta T$ (2.15) and the 'effective'
temperature difference, and we shall show that the latter is related to the former. For the bulk evolution of the whole reservoir, it is more appropriate to use the bulk scale $\Delta T_H$.

Convection in the lower layer is driven by internal heat sources only, and hence can be characterized by its local Rayleigh–Roberts number $Ra_{H1}$:

$$Ra_{H1} = \frac{\rho_{1o}g\alpha H_1 h_1^5}{\lambda\kappa\mu_1}. \qquad (2.16)$$

In contrast, the upper layer is heated by both internal heat sources and by the lower fluid. Thus, convection depends not only on the Rayleigh–Roberts number $Ra_{H2}$ for the layer, defined as above with the relevant properties, but also on the relative importance of the two types of heat input. We introduce two dimensionless numbers. The first one, noted $Ra_2$, relies on the total heat flux at the top of the upper fluid in steady-state conditions and can be written as a function of $Ra_{H2}$ and a second dimensionless number noted $E$:

$$Ra_2 = \frac{\rho_{2o}g\alpha \left( H_1 \frac{h_1}{h_2} + H_2 \right) h_2^5}{\lambda\kappa\mu_2} = Ra_{H2}(E+1), \qquad (2.17)$$

where $E$ is the ratio between the two different contributions to the heat flux at the top of the reservoir.

$$E = \frac{H_1 h_1}{H_2 h_2}. \qquad (2.18)$$

All the dimensionless numbers rely on the initial configuration of two layers separated by a horizontal interface, which may not be relevant in practice as the interface gets distorted and as the two fluids progressively mingle with one another. The viscosity of the final homogenized fluid can be estimated as follows (Bloomfield & Dewan 1971):

$$\mu = \mu_1 a \mu_2^{1-a}. \qquad (2.19)$$

Using the bulk heat production for the reservoir, $H$ (2.5), this allows calculation of a ‘bulk’ Rayleigh–Roberts number for the homogenized reservoir, denoted by $Ra_H$ as above.

It is useful to refer to the individual Rayleigh number for the two layers because it specifies their convection regimes in initial stages, but there are only five independent dimensionless numbers, which can be chosen to be $Ra_H$, $\gamma$, $F$, $a$ and $B$, to which one should add $Pr$ for completeness. All the other numbers can be derived from this list. For example:

$$E = \frac{aF}{1 - aF}. \qquad (2.20)$$

If $F=1$, the lower layer has the same amount of heat sources as the upper one and the basal temperature is that of a homogeneous fluid layer, $(1/2)\Delta T_H$. If $F = 1/a$, there are no heat sources in the upper layer and the bottom temperature is $\Delta T_H(1 - a/2)$, which is larger and illustrates the importance of heat sources at the base of the reservoir. Maintaining the constraint that $F = 1/a$, decreasing the value of $a$ ($a \to 0$) implies an increasing amount of heat sources in an increasingly thinner lower layer and produces a situation similar to Rayleigh–Bénard convection in the upper layer with a basal temperature set to $\Delta T_H$.
**3. Laboratory set-up and experimental protocol**

We rely on our new experimental set-up which was designed specifically to study convection due to internal heat sources. A complete description and discussion of measurement precision and accuracy may be found in Fourel *et al.* (2017).

**3.1. *Experimental techniques***

Working fluids were prepared in our laboratory in order to investigate a large range of dimensionless parameters. All were aqueous solutions of salt and hydroxyethylcellulose, which are fully miscible and whose intrinsic densities and viscosities can be varied within large ranges at small cost. Sodium chloride salt increases both density and microwave absorption. New fluids were used for each experiment and physical properties were determined over the relevant temperature range. Viscosity was measured to better than a few per cent uncertainty with a Thermo Scientific Haake rheometer RS600. Prandtl numbers were always larger than $10^2$, which ensures the dominance of viscous stresses over inertia, as in the Earth's mantle (Davaille & Limare 2015). Over the typical temperature range of an experiment (10°C), viscosity varies by a factor of about 0.6. A careful comparison between high-precision numerical calculations with or without such variations shows that viscosity changes of this magnitude have no significant effect on the variables of interest and on scaling laws (Limare *et al.* 2015). For density and thermal expansion coefficient, we used a DMA 5000 Anton Paar densimeter with an accuracy of one part per million. Thermal diffusivity and conductivity were determined by the photopyroelectric method (Dadarlat & Neamtu 2009), and dielectric properties with an Agilent N5230A vector network analyser. The salt diffusion coefficient in polymerized fluids such as ours was estimated for different polymer concentration and decreases with increasing polymer concentration (and hence with increasing fluid viscosity) (Davaille 1999b). At the interface between two layers of different viscosities, chemical diffusion is limited by the fluid with the smallest diffusion coefficient (in the most viscous one). In this series of experiments, relevant values of the diffusion coefficient lie in a $10^{-12}-10^{-13} \text{ m}^2 \text{ s}^{-1}$ range. For a reference time interval of 3 h (the typical duration of an experiment), these lead to diffusion lengths $\sqrt{Dt}$ in a 0.03–0.1 mm range. Thus, diffusion is not important at the scale of these experiments.

The experimental set-up is shown in figure 2. The 30 cm wide and 5 cm high tank is placed inside a modified microwave oven that achieves a laterally uniform microwave absorption and internal heat release (Surducan *et al.* 2014). The top is made of an aluminium heat exchanger connected to a thermostated bath which allows a constant temperature boundary condition. The bottom boundary is made of a thick poly(methyl-methacrylate) plate which can be considered as adiabatic. The performance of this laboratory set-up was assessed through a thorough comparison with numerical calculations reproducing the exact same configuration, including the tank dimensions and boundary conditions, for a homogeneous fluid layer over a large range of Rayleigh–Roberts numbers (Limare *et al.* 2015).

*In situ* temperature determinations were achieved via laser induced fluorescence (LIF) by adding a combination of two fluorescent dyes in order to separate between composition and temperature contributions to fluorescence. Spatial resolution is set by the 0.2 mm pixel size of the digital camera, enabling excellent resolution of the temperature profile in thermal boundary layers at the top of the tank, which were always thicker than 2 mm. A laser sheet scans half of the tank whilst two CCD cameras acquire images in different spectral ranges, so that temperature and
<figure><img src="image_123.png" /><figcaption>FIGURE 2. (Colour online) Experimental set-up.</figcaption></figure>

composition can both be measured simultaneously. In addition, a particle image velocimetry system allows determination of the velocity field. Three-dimensional distributions are constructed by interpolation of the two-dimensional data sets. Scans were obtained at a spacing of 1 cm over a half-tank width (15 cm), which allows a representative sampling of the different fields. Dye concentrations in the lower layer are three times larger than in the upper one, which makes for a sharp change of fluorescence between the two fluids and allows us to track the interface that separates them.

Each experiment followed the same protocol. The fluids were left for at least 12 h in the temperature-controlled laboratory so that they were initially at room temperature. The tank was first completely filled with the upper fluid and air bubbles were carefully removed. Next, a known volume of denser lower layer fluid was injected at the bottom whilst the excess upper fluid was removed. The dense fluid was left to spread across the whole tank and settle in a layer of uniform thickness. This generated an initially stably stratified system. At the start of an experiment, thermostated water at room temperature was made to circulate through the heat exchanger at the top and the microwave source was turned on at some prescribed power. In these conditions, both layers were heating simultaneously at rates that depended on their respective microwave absorption intensities. We did not start by cooling an initially hot initial layer because we sought to avoid an early phase of convection driven exclusively by boundary layer instabilities at the top of the upper fluid independently of internal heat production.

## 3.2. *Experiments*

Additional information on fluids properties and other experimental conditions are given in the supplementary material (table S1), available at https://doi.org/10.1017/jfm.2019.243. We performed 38 experiments investigating large ranges of dimensionless numbers (table 1). For the sake of simplicity, we restricted our attention to lower layers that were thinner and with a higher rate of heat generation than the upper one, which is relevant to the undepleted reservoir that is likely to exist at the base
<table><caption>TABLE 1. List of experiments performed in this study, their dimensionless numbers, their temperature scales and convection regime: domes (D) and stratified (S); $B_{cond}$ was obtained using (2.8) with $\Delta T$ the steady-state conduction temperature scale (2.15). Rayleigh numbers are calculated for fluid properties at the volume-averaged temperature, whereas the viscosity ratio is given at a reference temperature (temperature of the top surface).</caption><thead><tr><th>Number</th><th>γ</th><th>a</th><th>F</th><th>B<sub>cond</sub></th><th>Ra<sub>H1</sub></th><th>Ra<sub>H2</sub></th><th>E</th><th>Ra<sub>H</sub></th><th>ΔT<sub>H</sub> (°C)</th><th>ΔT (°C)</th><th>Regime</th></tr></thead><tbody><tr><td>1</td><td>1.19</td><td>0.28</td><td>2.59</td><td>0.68</td><td>751</td><td>1.4 × 10<sup>4</sup></td><td>2.6</td><td>1.8 × 10<sup>5</sup></td><td>64.7</td><td>23.4</td><td>D</td></tr><tr><td>2</td><td>0.16</td><td>0.26</td><td>2.68</td><td>0.54</td><td>1.5 × 10<sup>3</sup></td><td>6.7 × 10<sup>3</sup></td><td>2.3</td><td>1.2 × 10<sup>5</sup></td><td>64.7</td><td>23.0</td><td>D</td></tr><tr><td>3</td><td>3.18</td><td>0.32</td><td>2.41</td><td>0.67</td><td>492</td><td>9.1 × 10<sup>3</sup></td><td>3.4</td><td>1.3 × 10<sup>5</sup></td><td>64.7</td><td>24.0</td><td>D</td></tr><tr><td>4</td><td>1.77</td><td>0.30</td><td>2.57</td><td>0.76</td><td>396</td><td>5.7 × 10<sup>3</sup></td><td>3.4</td><td>8.9 × 10<sup>4</sup></td><td>64.7</td><td>24.2</td><td>D</td></tr><tr><td>5</td><td>2.11</td><td>0.30</td><td>2.77</td><td>1.01</td><td>490</td><td>5.7 × 10<sup>3</sup></td><td>4.9</td><td>1.1 × 10<sup>5</sup></td><td>64.7</td><td>25.5</td><td>D</td></tr><tr><td>6</td><td>2.58</td><td>0.26</td><td>2.78</td><td>0.79</td><td>937</td><td>5.9 × 10<sup>4</sup></td><td>2.6</td><td>5.5 × 10<sup>5</sup></td><td>64.7</td><td>23.6</td><td>D</td></tr><tr><td>7</td><td>1.24</td><td>0.24</td><td>3.13</td><td>1.10</td><td>1.5 × 10<sup>3</sup></td><td>4.8 × 10<sup>4</sup></td><td>3.2</td><td>5.9 × 10<sup>5</sup></td><td>64.7</td><td>24.7</td><td>S</td></tr><tr><td>8</td><td>4.61</td><td>0.18</td><td>3.10</td><td>1.55</td><td>25</td><td>3.2 × 10<sup>4</sup></td><td>1.3</td><td>1.3 × 10<sup>5</sup></td><td>32.4</td><td>10.5</td><td>S</td></tr><tr><td>9</td><td>4.61</td><td>0.30</td><td>2.48</td><td>3.03</td><td>114</td><td>4.5 × 10<sup>4</sup></td><td>3.0</td><td>4.9 × 10<sup>4</sup></td><td>14.8</td><td>5.4</td><td>S</td></tr><tr><td>10</td><td>0.07</td><td>0.20</td><td>3.13</td><td>2.74</td><td>286</td><td>3.0 × 10<sup>3</sup></td><td>1.7</td><td>3.4 × 10<sup>4</sup></td><td>14.8</td><td>5.1</td><td>S</td></tr><tr><td>11</td><td>0.07</td><td>0.20</td><td>3.13</td><td>1.25</td><td>626</td><td>6.6 × 10<sup>3</sup></td><td>1.7</td><td>7.4 × 10<sup>4</sup></td><td>32.4</td><td>11.1</td><td>D</td></tr><tr><td>12</td><td>0.20</td><td>0.12</td><td>3.52</td><td>1.32</td><td>34</td><td>2.5 × 10<sup>4</sup></td><td>0.73</td><td>8.9 × 10<sup>4</sup></td><td>32.4</td><td>9.3</td><td>D</td></tr><tr><td>13</td><td>0.20</td><td>0.24</td><td>2.78</td><td>2.40</td><td>393</td><td>3.7 × 10<sup>3</sup></td><td>2.0</td><td>6.7 × 10<sup>4</sup></td><td>14.8</td><td>5.2</td><td>S</td></tr><tr><td>14</td><td>0.30</td><td>0.22</td><td>2.89</td><td>2.49</td><td>267</td><td>6.9 × 10<sup>3</sup></td><td>1.7</td><td>5.2 × 10<sup>4</sup></td><td>14.8</td><td>5.1</td><td>S</td></tr><tr><td>15</td><td>0.30</td><td>0.16</td><td>3.55</td><td>1.68</td><td>148</td><td>2.4 × 10<sup>4</sup></td><td>1.3</td><td>1.4 × 10<sup>5</sup></td><td>32.4</td><td>10.7</td><td>D</td></tr><tr><td>16</td><td>0.69</td><td>0.20</td><td>2.79</td><td>2.09</td><td>108</td><td>1.5 × 10<sup>4</sup></td><td>1.3</td><td>8.8 × 10<sup>4</sup></td><td>14.8</td><td>4.7</td><td>S</td></tr><tr><td>17</td><td>2.09</td><td>0.17</td><td>3.11</td><td>3.15</td><td>30</td><td>3.1 × 10<sup>4</sup></td><td>1.1</td><td>1.2 × 10<sup>5</sup></td><td>14.8</td><td>4.6</td><td>S</td></tr><tr><td>18</td><td>10.4</td><td>0.20</td><td>2.87</td><td>1.67</td><td>9</td><td>1.8 × 10<sup>4</sup></td><td>1.3</td><td>6.5 × 10<sup>4</sup></td><td>32.4</td><td>10.5</td><td>S</td></tr><tr><td>19</td><td>1.05</td><td>0.14</td><td>3.02</td><td>0.67</td><td>37</td><td>7.5 × 10<sup>4</sup></td><td>0.73</td><td>2.4 × 10<sup>5</sup></td><td>48.0</td><td>13.6</td><td>D</td></tr><tr><td>20</td><td>1.05</td><td>0.14</td><td>3.02</td><td>0.50</td><td>50</td><td>1.0 × 10<sup>5</sup></td><td>0.73</td><td>3.2 × 10<sup>5</sup></td><td>64.7</td><td>18.4</td><td>D</td></tr><tr><td>21</td><td>23.3</td><td>0.14</td><td>4.02</td><td>0.95</td><td>21</td><td>4.8 × 10<sup>5</sup></td><td>1.3</td><td>1.3 × 10<sup>6</sup></td><td>92.5</td><td>30.6</td><td>D</td></tr><tr><td>22</td><td>12.2</td><td>0.12</td><td>4.23</td><td>1.02</td><td>9</td><td>3.0 × 10<sup>5</sup></td><td>1.0</td><td>7.5 × 10<sup>5</sup></td><td>87.8</td><td>27.7</td><td>D</td></tr><tr><td>23</td><td>0.12</td><td>0.10</td><td>4.46</td><td>1.67</td><td>227</td><td>2.0 × 10<sup>5</sup></td><td>0.80</td><td>6.8 × 10<sup>5</sup></td><td>48.1</td><td>14.4</td><td>S</td></tr><tr><td>24</td><td>5.0</td><td>0.14</td><td>4.29</td><td>1.06</td><td>64</td><td>2.7 × 10<sup>5</sup></td><td>1.5</td><td>1.0 × 10<sup>6</sup></td><td>92.5</td><td>31.5</td><td>S</td></tr><tr><td>25</td><td>9.0</td><td>0.10</td><td>4.81</td><td>1.17</td><td>13</td><td>7.8 × 10<sup>5</sup></td><td>0.92</td><td>1.8 × 10<sup>6</sup></td><td>92.5</td><td>28.6</td><td>S</td></tr><tr><td>26</td><td>38.5</td><td>0.08</td><td>4.28</td><td>1.00</td><td>2</td><td>2.6 × 10<sup>6</sup></td><td>0.52</td><td>4.1 × 10<sup>6</sup></td><td>92.4</td><td>24.7</td><td>S</td></tr><tr><td>27</td><td>3.1</td><td>0.16</td><td>3.08</td><td>1.24</td><td>3.8</td><td>9.1 × 10<sup>3</sup></td><td>0.97</td><td>3.0 × 10<sup>4</sup></td><td>32.3</td><td>9.85</td><td>D</td></tr><tr><td>28</td><td>0.96</td><td>0.14</td><td>4.55</td><td>1.07</td><td>2.6 × 10<sup>3</sup></td><td>1.9 × 10<sup>6</sup></td><td>1.7</td><td>9.7 × 10<sup>6</sup></td><td>92.5</td><td>32.9</td><td>S</td></tr><tr><td>29</td><td>0.96</td><td>0.18</td><td>4.04</td><td>0.92</td><td>1.8 × 10<sup>4</sup></td><td>2.6 × 10<sup>6</sup></td><td>2.7</td><td>2.1 × 10<sup>7</sup></td><td>102</td><td>38.6</td><td>S</td></tr><tr><td>30</td><td>0.82</td><td>0.15</td><td>4.05</td><td>0.90</td><td>2.2 × 10<sup>3</sup></td><td>1.1 × 10<sup>6</sup></td><td>1.5</td><td>5.7 × 10<sup>6</sup></td><td>92.5</td><td>31.8</td><td>S</td></tr><tr><td>31</td><td>0.68</td><td>0.26</td><td>3.09</td><td>0.80</td><td>3.5 × 10<sup>4</sup></td><td>3.5 × 10<sup>5</sup></td><td>4.1</td><td>6.6 × 10<sup>6</sup></td><td>92.5</td><td>36.2</td><td>D</td></tr><tr><td>32</td><td>1.05</td><td>0.30</td><td>2.87</td><td>0.77</td><td>1.8 × 10<sup>5</sup></td><td>8.5 × 10<sup>5</sup></td><td>6.2</td><td>2.5 × 10<sup>7</sup></td><td>92.5</td><td>37.4</td><td>D</td></tr><tr><td>33</td><td>1.05</td><td>0.10</td><td>4.81</td><td>0.89</td><td>1.4 × 10<sup>3</sup></td><td>9.5 × 10<sup>6</sup></td><td>0.92</td><td>2.8 × 10<sup>7</sup></td><td>102</td><td>31.6</td><td>S</td></tr><tr><td>34</td><td>1.19</td><td>0.14</td><td>5.26</td><td>0.74</td><td>1.4 × 10<sup>4</sup></td><td>7.5 × 10<sup>6</sup></td><td>2.8</td><td>5.1 × 10<sup>7</sup></td><td>162</td><td>63.0</td><td>S</td></tr><tr><td>35</td><td>1.18</td><td>0.28</td><td>3.32</td><td>0.68</td><td>3.0 × 10<sup>5</sup></td><td>1.0 × 10<sup>6</sup></td><td>14</td><td>5.2 × 10<sup>7</sup></td><td>162</td><td>69.7</td><td>D</td></tr><tr><td>36</td><td>0.11</td><td>0.24</td><td>2.84</td><td>0.55</td><td>5.1 × 10<sup>3</sup></td><td>2.6 × 10<sup>4</sup></td><td>2.1</td><td>4.1 × 10<sup>5</sup></td><td>48.0</td><td>17.0</td><td>D</td></tr><tr><td>37</td><td>0.06</td><td>0.21</td><td>2.88</td><td>0.87</td><td>3.5 × 10<sup>3</sup></td><td>2.8 × 10<sup>4</sup></td><td>1.5</td><td>3.2 × 10<sup>5</sup></td><td>32.3</td><td>10.8</td><td>D</td></tr><tr><td>38</td><td>0.66</td><td>0.36</td><td>2.63</td><td>0.56</td><td>4.7 × 10<sup>5</sup></td><td>1.5 × 10<sup>5</sup></td><td>19</td><td>2.1 × 10<sup>7</sup></td><td>162</td><td>68.5</td><td>D</td></tr></tbody></table>

of the Earth's mantle. The lower fluid was allowed to be more and less viscous than the upper one and viscosity contrasts between the two fluids were varied over more than two orders of magnitude ($0.06 \le \gamma \le 38$). We were interested in large
Rayleigh-Roberts values appropriate for the Earth's mantle and did not investigate conditions near the stability threshold. In all experiments Ra₂ > RaH₁, due to the large impact of the layer thickness on the Rayleigh-Roberts number, a condition that is appropriate for the Earth's mantle. Based on values of Rayleigh-Roberts number of the lower layer, we investigated cases such that the lower layer would have been able or unable to undergo convection on its own (table 1). Values for the Rayleigh-Roberts for the final homogenized reservoir were in a 3 × 10⁴ – 5 × 10⁷ range, which straddles the threshold between steady and time-dependent convection regimes (Vilella et al. 2018).

# 4. Thermal structure and heat flux through the upper boundary

## 4.1. *Time evolution*

All experiments followed the same basic evolution, with an interface that started to deform when convection set in. The volume of lower fluid separated from the upper fluid by a well-defined interface decreased steadily until there was no longer evidence for the existence of two different fluid regions in the tank. We describe later how the interface deformed and how mingling proceeded and discuss here the evolution of convection from a purely thermal perspective.

We derived vertical profiles of the horizontally averaged temperature and calculated the heat flux at the top of the tank using a least-squares linear fit to the three uppermost temperature values. Figure 3 shows how the vertical profile of the horizontally averaged temperature changes with time for experiment 21. At small times, the profiles bear the influence of two different fluid layers, with an upper region that appears well mixed beneath an upper boundary layer and a lower region with a large negative temperature gradient. Later in the experiment, there is little evidence for a basal fluid layer, which is due to the large deformations of the interface, such that the lower fluid is not confined below a well-defined horizon (as later shown in § 5.2). With increasing time, upper regions warm up whilst lower ones cool down until a steady-state profile with the hallmarks of an internally heated fluid layer is achieved, with a small but well-defined stable temperature gradient in the fluid interior (Vilella et al. 2018). This late evolution proceeds by changes of internal thermal structure that maintain the volume-average temperature and the heat flux at the top almost constant. At these late times, there is no longer a lower layer at the base but one can still observe thin slivers of lower fluid in the tank interior. There is no detectable heterogeneity at the scale of the fluid motions and the reservoir can be considered as both homogeneous from a thermal standpoint and heterogeneous from a chemical standpoint.

We have tracked the time evolution of convection using the heat flux $\phi$ at the top and the volume-averaged temperature in the tank $\langle T \rangle$ (figure 4). Both variables tended towards steady-state values, noted $\phi_\infty$ and $T_{vol}$, over a well-defined time. There were small and simultaneous fluctuations in both the heat flux and temperature records, which will be discussed later. At steady state, the heat flux at the top was equal to the total amount of heat released in the tank interior, as required in equilibrium conditions.

To demonstrate that the apparent thermal steady state corresponds to that of a homogeneous fluid, we checked that the two thermal structures conform to the same scaling laws. This is evaluated in figure 5, where the steady-state volume-averaged temperature $T_{vol}$ is shown as a function of the ‘bulk’ Rayleigh-Roberts number RaH (4.1). Also shown in this figure are data from experiments in homogeneous fluids published in a previous study (Limare et al. 2015) and listed in the supplementary
<figure><img src="image_124.png" /><figcaption>FIGURE 3. (Colour online) Vertical profiles of the horizontally averaged temperature as a function of time for experiment 21. Colours indicate time evolution from green to blue, red and black. Profiles are taken every 10 min from green (0 min) to blue, red and black (180 min). The profiles tend to that for a homogeneous fluid layer in equilibrium with its heat sources, with a small stable interior bulk thermal stratification.</figcaption></figure>

<figure><img src="image_125.png" /><figcaption>FIGURE 4. (Colour online) Dimensionless heat flux (a) and volume-averaged temperature (b) as a function of time normalized to the conduction time for the whole layer, $τ_c = h^2/κ$, for experiment 21. Black diamonds represent experimental data, black lines represent exponential fits and red, dashed lines represent the results of transient calculations from (4.3) and (4.5).</figcaption></figure>

material in table S2. At large values of $Ra_H$ (>10⁵), Vilella et al. (2018) have shown that the two variables are related to one another by the following scaling law:

$$ \frac{T_{vol}}{\Delta T_H} = C R a_H^\beta, \qquad (4.1) $$

with an exponent $β$ close to the value of -1/4 for Boussinesq fluids with constant physical properties (table 3); C is a constant depending on the mechanical boundary conditions, rigid in our case. Within their error ranges, best-fit values of the two