<figure><img src="image_126.png" /><figcaption>FIGURE 5. Dimensionless steady-state volume-averaged temperature as a function of Rayleigh–Roberts number: black diamonds represent data from this study (heterogeneous convection experiments) and empty squares data from experiments in homogeneous fluids (Limare et al. 2015). Lines represent fits obtained with a fixed exponent (table 3); heterogeneous (full line) and homogeneous (dashed line).</figcaption></figure>

parameters in the scaling law are equal to those for homogeneous fluids. Moreover, they are almost identical to those obtained numerically in isoviscous and infinite Prandtl fluids encased in tanks with rigid walls (Vilella et al. 2018) (table 3).

In the final steady state, the convective heat flux at the top evacuates all the heat released within the layer. Thus, equation (4.1) for the volume-averaged temperature can be turned into an equation for the heat flux at the top noted $\phi_{\infty}$:

$$ \phi_{\infty} = \left(\frac{k}{Ch}\right)^{1/(1+\beta)} \left(\frac{g\alpha\rho h^4}{k\kappa\mu}\right)^{-\beta/(1+\beta)} T_{vol}^{1/(1+\beta)}. \quad (4.2) $$

This heat flux is transported by conduction through the upper unstable boundary layer. We note that, for $\beta = -1/4$, $\phi_{\infty} \propto T_{vol}^{4/3}$ and the total fluid thickness $h$ gets cancelled in this expression. In this case, the dynamics of the upper boundary layer is controlled locally, independently of the deep fluid region below. As shown by Vilella et al. (2018), the volume-averaged temperature $T_{vol}$ and the temperature difference across the upper boundary layer, noted $\Delta T_{TBL}$, are both scaled to $Ra_H^{-1/4}$, so that $\phi_{\infty} \propto \Delta T_{TBL}^{4/3}$, which conforms to the well-known local heat flux law for actively convecting layers (Townsend 1964).

## 4.2. Transient thermal evolution

The transient thermal evolution of the whole fluid is governed by the heat balance equation:

$$ \rho c_p h \frac{d\langle T \rangle}{dt} = Hh - \phi, \quad (4.3) $$
<table><caption>TABLE 2. List of experiments performed in this study, their characteristic volume-average temperatures and time constants. $T_{vol}$ is the steady-state volumetrically averaged temperature and $\tau_e$ is the thermal relaxation time according to (4.6); $T_{12max}$ is the maximum value of the volume-average temperature contrast between the two layers and $\tau_{12max}$ is the time at which this maximum occurs.</caption><thead><tr><th>Number</th><th><i>T</i><sub>vol</sub> (°C)</th><th><i>T</i><sub>12max</sub> (°C)</th><th><i>τ</i><sub>e</sub> (min)</th><th><i>τ</i><sub>12max</sub> (min)</th><th><i>t</i><sub>0</sub> (min)</th><th><i>t</i><sub>end</sub> (min)</th></tr></thead><tbody><tr><td>1</td><td>9.0</td><td>13.0</td><td>34</td><td>40</td><td>45</td><td>150</td></tr><tr><td>2</td><td>13.3</td><td>14.8</td><td>39</td><td>40</td><td>30</td><td>60</td></tr><tr><td>3</td><td>11.7</td><td>14.0</td><td>34</td><td>70</td><td>41</td><td>174</td></tr><tr><td>4</td><td>12.9</td><td>10.6</td><td>48</td><td>50</td><td>42</td><td>155</td></tr><tr><td>5</td><td>14.3</td><td>10.6</td><td>52</td><td>65</td><td>55</td><td>228</td></tr><tr><td>6</td><td>7.1</td><td>8.2</td><td>25</td><td>70</td><td>54</td><td>166</td></tr><tr><td>7</td><td>7.1</td><td>8.2</td><td>30</td><td>70</td><td>128</td><td>200</td></tr><tr><td>8</td><td>5.5</td><td>5.7</td><td>40</td><td>110</td><td>50</td><td>470</td></tr><tr><td>9</td><td>4.7</td><td>4.0</td><td>72</td><td>170</td><td>70</td><td>2000</td></tr><tr><td>10</td><td>2.9</td><td>4.1</td><td>56</td><td>200</td><td>60</td><td>500</td></tr><tr><td>11</td><td>6.6</td><td>7.5</td><td>45</td><td>100</td><td>40</td><td>206</td></tr><tr><td>11</td><td>5.5</td><td>5.8</td><td>37</td><td>50</td><td>24</td><td>118</td></tr><tr><td>12</td><td>4.5</td><td>4.2</td><td>70</td><td>120</td><td>110</td><td>660</td></tr><tr><td>13</td><td>3.8</td><td>2.7</td><td>64</td><td>150</td><td>100</td><td>620</td></tr><tr><td>14</td><td>6.0</td><td>6.2</td><td>43</td><td>90</td><td>39</td><td>258</td></tr><tr><td>15</td><td>3.0</td><td>2.4</td><td>49</td><td>100</td><td>58</td><td>374</td></tr><tr><td>16</td><td>2.7</td><td>2.0</td><td>41</td><td>90</td><td>100</td><td>853</td></tr><tr><td>18</td><td>6.8</td><td>6.0</td><td>54</td><td>200</td><td>80</td><td>631</td></tr><tr><td>19</td><td>7.4</td><td>4.3</td><td>44</td><td>45</td><td>8</td><td>65</td></tr><tr><td>20</td><td>9.2</td><td>5.2</td><td>35</td><td>40</td><td>12</td><td>40</td></tr><tr><td>21</td><td>9.6</td><td>9.3</td><td>23</td><td>80</td><td>30</td><td>200</td></tr><tr><td>22</td><td>10.1</td><td>9.9</td><td>26</td><td>30</td><td>10</td><td>140</td></tr><tr><td>23</td><td>7.0</td><td>7.4</td><td>35</td><td>40</td><td>10</td><td>100</td></tr><tr><td>24</td><td>8.5</td><td>7.1</td><td>18</td><td>30</td><td>20</td><td>160</td></tr><tr><td>25</td><td>7.9</td><td>5.7</td><td>18</td><td>30</td><td>20</td><td>170</td></tr><tr><td>26</td><td>7.4</td><td>6.9</td><td>19</td><td>40</td><td>30</td><td>190</td></tr><tr><td>27</td><td>9.2</td><td>5.6</td><td>52</td><td>60</td><td>50</td><td>210</td></tr><tr><td>28</td><td>6.3</td><td>5.6</td><td>16</td><td>20</td><td>50</td><td>100</td></tr><tr><td>29</td><td>6.3</td><td>4.6</td><td>17</td><td>20</td><td>90</td><td>130</td></tr><tr><td>30</td><td>5.9</td><td>7.9</td><td>24</td><td>55</td><td>50</td><td>85</td></tr><tr><td>31</td><td>5.9</td><td>9.0</td><td>21</td><td>40</td><td>90</td><td>125</td></tr><tr><td>32</td><td>6.2</td><td>7.9</td><td>16</td><td>25</td><td>92</td><td>115</td></tr><tr><td>33</td><td>5.5</td><td>6.0</td><td>16</td><td>20</td><td>34</td><td>51</td></tr><tr><td>34</td><td>5.5</td><td>8.2</td><td>8</td><td>15</td><td>50</td><td>85</td></tr><tr><td>35</td><td>6.5</td><td>12.3</td><td>9</td><td>20</td><td>65</td><td>75</td></tr><tr><td>36</td><td>6.7</td><td>6.0</td><td>39</td><td>40</td><td>18</td><td>49</td></tr><tr><td>37</td><td>6.5</td><td>3.1</td><td>36</td><td>35</td><td>30</td><td>75</td></tr><tr><td>38</td><td>8.3</td><td>11.1</td><td>11</td><td>20</td><td>65</td><td>80</td></tr></tbody></table>

where $\langle T \rangle$ and $\phi$ are the time-dependent volume-averaged temperature and heat flux at the top of the tank, respectively. This equation can also be rewritten as an equation for the heat flux, which evacuates heat produced internally as well as sensible heat extracted through cooling (appearing as $-\rho c_p (d\langle T \rangle/dt)$). In transient conditions, the unstable boundary layer at the top is much thinner than the whole fluid and rapidly
<table><caption>TABLE 3. Parameters of empirical best-fit power laws for the volume-average ($T_{vol}/\Delta T_{H}$) temperature in steady state at high values of the Rayleigh-Roberts number ($Ra_{H} \ge 10^{5}$) for heterogeneous and homogeneous laboratory experiments. Data for homogeneous fluid layers are taken from the experiments by Limare et al. (2015).</caption><thead><tr><th></th><th></th><th>Exponent left to vary</th><th>Exponent fixed to -1/4</th></tr></thead><tbody><tr><td>T<sub>vol</sub>/&Delta;T<sub>H</sub></td><td>Heterogeneous</td><td>3.27(&plusmn;1.44) Ra<sub>H</sub><sup>-0.245(&plusmn;0.037)</sup></td><td>3.51(&plusmn;0.18) Ra<sub>H</sub><sup>-1/4</sup></td></tr><tr><td>T<sub>vol</sub>/&Delta;T<sub>H</sub></td><td>Homogeneous</td><td>2.07(&plusmn;1.44) Ra<sub>H</sub><sup>-0.216(&plusmn;0.036)</sup></td><td>3.38(&plusmn;0.16) Ra<sub>H</sub><sup>-1/4</sup></td></tr><tr><td>T<sub>vol</sub>/&Delta;T<sub>H</sub></td><td>Numerical, Boussinesq</td><td>/</td><td>3.20(&plusmn;0.11) Ra<sub>H</sub><sup>-1/4</sup></td></tr></tbody></table>

adjusts to the interior thermal structure. One can thus assume that the local heat flux expression (4.2) remains valid at all times as a function of the instantaneous value of the volume-averaged temperature $\langle T \rangle$ (figure 3). This standard and well-tested approximation was validated by the careful laboratory experiments of Katsaros et al. (1977). Scaling temperature by the steady-state value $T_{vol}$ in the heat balance equation leads to the following time scale:

$$ \tau = C \frac{h^2}{\kappa} Ra_H^\beta = C \tau_c Ra_H^\beta, \qquad (4.4) $$

where $\tau_c = h^2/\kappa$ is the diffusive time scale for the whole fluid layer. Using these scales, the dimensionless heat balance equation is:

$$ \frac{d\langle T^* \rangle}{dt} = 1 - \langle T^* \rangle^{1/(1+\beta)}, \qquad (4.5) $$

where $\langle T^* \rangle = \langle T \rangle / T_{vol}$ is the dimensionless volume-averaged temperature. Integrating this equation, we find that calculated values of the heat flux and volume-averaged temperature are very close to the measured ones at all times (figure 4).

The agreement between data and model predictions that is shown in figure 4 is obtained for an experiment with a low $B$ value and a relatively large $Ra_H$ value (experiment 21, table 1). We have checked that this is true for all experiments. In particular, we show results for an experiment with a high $B$ value and a low $Ra_H$ value (experiment 10) in the supplementary material (figure S1). Differences between the experimental data and the theoretical predictions are slightly larger than for experiments at larger values of $Ra_H$. This is expected because the $\phi_\infty - T_{vol}$ scaling law with power-law exponent $\beta = -1/4$ is only valid for $Ra_H > 6 \times 10^5$ (Limare et al. 2015; Vilella et al. 2018), such that the upper boundary layer is very much thinner than the reservoir. Deviations of $T_{vol}$ from this scaling account for most of the scatter in figure 5. In order to verify further that the experimental data are consistent with these scalings, we have determined the duration of the thermal transient, independently of any model calculation or scaling. One procedure would consist of setting a threshold value for either the heat flux or volume-averaged temperature at some fraction of their final steady-state values and determining the times when this is reached. Results depend on the rather arbitrarily chosen threshold and on the accuracy of temperature determinations, therefore we used another method. The evolution of both the heat flux and the volume-averaged temperature seem to be
<figure><img src="image_127.png" /><figcaption>FIGURE 6. Dimensionless thermal relaxation time $\tau_e/\tau_c$ for convection as a function of the bulk Rayleigh–Roberts number $Ra_H$. Line represents the fit obtained by setting the power-law exponent to the value of $-1/4$ (table 4).</figcaption></figure>

<table><caption>TABLE 4. Parameters of empirical best-fit power laws for the characteristic time constant $\tau_e$.</caption><thead><tr><th></th><th>Exponent left to vary</th><th>Exponent fixed to -1/4</th></tr></thead><tbody><tr><td>τ<sub>e</sub>/τ<sub>c</sub></td><td>2.00(±0.46) Ra<sub>H</sub><sup>-0.222(±0.038)</sup></td><td>3.00(±0.23) Ra<sub>H</sub><sup>-1/4</sup></td></tr></tbody></table>

very close to a simple exponential relaxation (figure 4). This allows the determination of a characteristic time $\tau_e$ such that, for example, the time-dependent heat flux is:

$$ \phi(t) = Hh(1 - e^{-t/\tau_e}). \quad (4.6) $$

The exponential approximation of (4.6) allows a very good fit to the experimental data (figure 4a) and we have determined $\tau_e$ values for all experiments through a best-fit procedure. With this method, we can use both the temperature and heat flux data, and we can assess the overall trend towards equilibrium. Further, substituting for this approximate expression in the bulk heat balance equation leads to a very good agreement with the temperature measurements (figure 4b). If the above scaling arguments are correct, $\tau_e$ should scale with $\tau_c Ra_H^\beta$, which is tested successfully in figure 6 and table 4. For convenience, we shall use the value of $\beta = -1/4$ throughout the paper, rather than some empirical value within an uncertainty range.

For planetary studies, the important implication is that the bulk thermal evolution of a planet (which is usually referred to as ‘secular cooling’) is not sensitive to the exact distribution of its internal heat sources, for reasons that will be explained below. As shown below, however, other features of convection, such as planforms and distribution and size of upwellings, depend on the distribution of heat sources as well as on differences in the properties of the two fluids.
**5. Two different convection regimes**

The topography of the interface that separates the two fluids provides a natural marker of flow structure and is of particular interest for comparison with tomographic images of the Earth's mantle. As argued by Foureël *et al.* (2017), deformation is both a consequence and a driver of convective motions due to enhanced heat generation in the lower fluid. The interface is well defined initially because it is characterized by an almost step change of composition and transmitted light intensity. As time progresses, however, mixing of the two fluids is generated at the interface which therefore becomes blurred. We identify an interface when there is a core of lower fluid with a composition that is close to the initial value. Such a core is separated from neighbouring fluid by a sharp composition gradient and we locate the interface using a threshold composition value. Due to the sharp gradient, the position of the interface is weakly sensitive to the exact threshold value that is adopted, as shown in the supplementary material. At large distances from core regions, thin threads of lower fluid can survive for very long times but they are distributed through large volumes and do not affect the fluid motion. In these conditions, they behave as passive tracers (see § 6.1).

In all our experiments, there was some entrainment of upper fluid into the lower fluid, as described later, but it was the upper fluid that ended up ingesting all the lower fluid. In principle, mixing could proceed the other way around with the lower fluid engulfing increasing amounts of upper fluid until the mixture occupies the whole tank. This did not happen in the present set of experiments because convective motions are systematically more intense in the upper layer, due to the higher Rayleigh-Roberts number.

The interface topography is defined as the distribution of interface height above the tank base and develops in two different regimes which can be understood using two end-member cases. For a small density and viscosity contrasts, the interface does not act as a barrier and deforms passively, such that the two fluids undergo convective motions that are determined at the scale of the whole tank ('doming' regime). For large values of the density contrast, the interface remains flat but mixing still occurs, due to the shearing and extraction of thin slivers of lower fluid ('stratified' regime). In principle, one should differentiate between mixing and mingling although they are parts of the same process in miscible fluids. They both describe the entrainment of one fluid into the other, which proceeds in two different ways. One is such that one fluid protrudes into the other and folds over, engulfing parts of the other fluid. In the other process, shear at the interface acts to tear a thin schlieren out of one fluid. These two processes are involved to different degrees in the two regimes of interface deformation. They are both such that entrained fluid parcels get stretched and thinned, up to a point when diffusion becomes effective and eradicates composition differences. Mixing applies to this whole sequence and mingling describes the intermediate process of generating small parcels of one fluid that get carried by the other one. Mingling leads to mixtures that are heterogeneous at a local scale and yet can be homogeneous at large scales.

**5.1. *The stratified regime***

All else being equal, this regime is observed for intrinsically denser lower fluids than in the dome regime, which prevents large-scale protrusions of lower fluid into the upper one (movie 2 in supplementary material). Figure 7 shows three different snapshots of the flow structure, which document the topography of the lower layer
<figure><img src="image_128.png" /><figcaption>FIGURE 7. (Colour online) Interface topography and distribution of the average vertical velocity above and below the interface at $t = 40$ min (a), 100 min (b) and 160 min (c) for experiment 7 in the stratified regime. Space dimensions are in mm.</figcaption></figure>

together with the horizontal distribution of the average vertical velocity field in each layer. The interface develops a morphology of cusp-like ridges encircling basins. The basins grow deeper with time and eventually extend to the base of the tank, which becomes exposed to the upper fluid in several locations. These areas gradually widen until they occupy the whole base of the tank. The cusps that delimit the basins are due to converging flows feeding upwellings in the upper fluid, whereas the basins are associated with downwellings (figure 7).

The cusp and basin morphology of the interface leads to a peaked probability density function (PDF) for the interface height with a long thinning tail at high values. In order to characterize these functions in a quantitative manner, we use gamma distributions:

$$f(x, k, \theta) = \frac{x^{k-1}e^{-(x/\theta)}}{\Gamma(k)\theta^k}, \qquad (5.1)$$

where $\Gamma$ stands for the gamma function, $\theta$ is a scale factor and $k$ is a parameter. As shown in figure 8, these distributions allow a good fit to the experimental data, as shown by the very large values of the correlation coefficient ($R > 0.99$). More importantly, these distributions allow us to track how the distribution changes with time. At $t = 0$ min, one would expect the cumulative distribution (CDF) to be a step function, corresponding to a perfectly flat interface between the two layers. At that initial recording time, the interface has already developed undulations and the height distribution shows up as a narrow peak, close to ‘normal’ distribution, which is the limiting form of the gamma distribution at large values of $k$. With time, as the cusp and basin morphology develops, the distribution changes and the $k$ parameter decreases steadily from an initial value of 150, to a value of 10 at $t = 40$ min and then to 4.5 at $t = 100$ min. At large times, the distribution is close to ‘exponential’, which is the limiting form of the gamma distribution for $k = 1$. At $t = 160$ min, for example, the best-fit $k$ value is 1.2 (figure 8).
<figure><img src="image_129.png" /><figcaption>FIGURE 8. (Colour online) Probability density function (a), cumulative density function (b) of the interface height at different times for experiment 7 in the stratified regime. Dashed line indicate gamma distribution function fits.</figcaption></figure>

## 5.2. *The doming regime*

We illustrate the characteristics of the dome regime using experiment 21 (table 1, movie 1 in supplementary material). In this regime, large domal structures grow steadily into columns that eventually reach the top of the tank. As they rise through the tank, these structures lose heat to the colder surroundings but do not fall back down readily as in Rayleigh-Bénard experiments (Davaille 1999b; Davaille, Le Bars & Carbonne 2003; Le Bars & Davaille 2004a,b). This is because they are made of lower fluid with a higher than average heat production rate, which sustains a temperature difference with the surrounding fluid. In some cases, the domes do deflate and go down, but they do not do so over a large depth range, which sets this behaviour apart from the ‘lava lamp’ oscillatory regime of double-diffusive convection (Turner 1974).

In experiment 21, convection develops first in the upper layer (figure 9a, top map) and acts to deform the lower layer. An internal circulation develops within the large domes, with a central upwelling surrounded by a peripheral downwelling (visible in the bottom part of figure 9b). This contrasts with the flow pattern in a homogeneous internally heated layer, which is characterized by narrow downwellings and a diffuse upward return flow (Vilella et al. 2018). Here, the upwellings of the upper layer get localized at the edges of the domes (figure 9b,c).

Figure 10 shows both the probability density function (PDF) and the cumulative density function (CDF) of the interface height. The interface is almost flat initially, such that the PDF is normal with a narrow peak and a small standard deviation. As convection develops, the PDF becomes increasingly asymmetrical (blue: 30 min and magenta: 60 min) with a peak that widens and eventually disappears (black curve). The gradual change of interface deformation is perhaps better illustrated by the cumulative density function (figure 10b), which starts as a step function and tends towards a linear function spanning the whole interval. These changes reflect the steady upward growth of domes which eventually stretch over the whole tank thickness.

The domes are large-scale coherent structures made of weakly diluted lower fluid which stands in upper fluid and gets sheared out in thin slivers. We shall document later how mixing proceeds and limit ourselves to the dome dimensions. The domes evolve into nearly cylindrical columns with diameters that barely change with height
<figure><img src="image_130.png" /><figcaption>FIGURE 9. (Colour online) Snapshots of the convective flow structure at three different times (a) 40 min (b) 90 min and (c) 130 min for experiment 21 in the dome regime. Three different types of data are shown at each time. The middle figure shows the interface topography, with grey areas indicating where the lower layer has been thinned to zero. The top and lower figure show the horizontal distribution of the averaged vertical velocity in the upper and lower fluids, respectively. Space dimensions are in mm.</figcaption></figure>

<figure><img src="image_131.png" /><figcaption>FIGURE 10. (Colour online) Probability density function (a), cumulative density function (b) of the interface height at different times for experiment 21 in the dome regime. Dashed lines indicate gamma distribution function fits except for $t = 160$ min (uniform distribution).</figcaption></figure>

above the base and remain stable for extended lengths of time (figure 11). Stable dome diameters at mid-height in the reservoir scale approximately with the lower layer thickness (figure 12).

Fits to the gamma distribution are shown with dashed lines in figure 10 for the dome regime. The distribution is again close to ‘normal’ at $t = 0$ min ($k = 100$), but it cannot be fit very well by gamma distributions as the domes grow. Correlation coefficients are smaller than for the stratified regime ($R \sim 0.98$), but still acceptable in principle. At late times, however, the distribution is clearly close to ‘uniform’ (with
<figure><img src="image_132.png" /><figcaption>FIGURE 11. (Colour online) Average diameter of domes scaled to the total fluid thickness *h* as a function of height above base at several times in experiment 21 (in the dome regime).</figcaption></figure>

<figure><img src="image_133.png" /><figcaption>FIGURE 12. Average dimensionless diameter of domes (scaled to the total reservoir thickness) as a function of dimensionless lower layer thickness (equivalent to the volume fraction of lower fluid in the tank). Data are from all experiments in the dome regime.</figcaption></figure>

$R > 0.99$ at $t = 160$ min, figure 10). This contrasts with the stratified regime, for which the height distribution is never close to uniform, as shown by low values of the correlation coefficient ($R < 0.8$).

5.3. Regime diagram

Some experiments are such that convection starts in the stratified regime and then shifts to the dome regime due to mixing, and they are lumped together with those that are in the dome regime at all times. Depending on the regime, the CDF varies with time in very different ways and tends towards two very different functions. The
<figure><img src="image_134.png" /></figure>

FIGURE 13. (Colour online) Regime diagrams in (Ra<sub>H</sub>, B) space. The buoyancy number B is calculated with the conduction temperature scale ΔT from (2.15) in (a) and with the measured temperature difference T<sub>12max</sub> in (b), corresponding to two different values of the buoyancy number noted B<sub>cond</sub> and B<sub>eff</sub> respectively. Triangles and squares stand for experiments in the dome and stratified regimes, respectively. Note that one experiment in the dome regime lies in the stratified domain for B = B<sub>cond</sub> and that it is shifted to the boundary between the two regimes for B = B<sub>eff</sub>.

goodness of fit to a uniform distribution is excellent in the dome regime and poor
in the stratified one, as shown above, which provides an unambiguous discrimination
criterion. All else being equal, the buoyancy ratio must be larger than some threshold
value to maintain a stratified regime over the whole duration of an experiment.

We found no measurable impact of the viscosity contrast on the boundary between
the two regimes, showing that the buoyancy number is the deciding control variable.
Figure 13 summarizes observations from all experiments in (B, RaH) space. The
buoyancy number is calculated using either the conduction scale ΔT from (2.15)
(figure 13a), corresponding to a value of the buoyancy number noted Bcond, or
the maximum temperature difference that develops between the two fluids in the
course of an experiment, noted T12max (see § 6), corresponding to a value noted Beff
(figure 13b). We discuss later the choice of this temperature difference, which scales
with ΔT RaH⁻¹/⁴ (see § 6.3).

The threshold *B* value that separates the two regimes varies by a factor of only 2 over the large range of bulk Rayleigh-Roberts numbers that has been investigated (figure 13). We show results for the two different definitions of the buoyancy number for the sake of discussion, but do not expect that the *B*<sub>cond</sub> value is appropriate since it relies on a static configuration and temperature distribution. When the *B*<sub>cond</sub> value is used, two features of the regime diagram are not satisfactory (figure 13a): one of the dome experiments (experiment 15) shows up in the stratified domain and the threshold buoyancy number for a stratified regime decreases with increasing Rayleigh-Roberts number. One expects instead that it takes increasingly large intrinsic density contrasts to stabilize the lower layer as convection becomes increasingly powerful. With the *B*<sub>eff</sub> values, these anomalies disappear (figure 13b). In this case, the boundary between the two regimes is characterized by *B*<sub>eff</sub> values in a 3–5 range.

In experiments with two fluids layers in a Rayleigh-Bénard set-up, with fixed
temperatures at the top and bottom of the tank, Davaille *et al.* (2003) observed
regimes that are similar to the present ones. They used the total temperature difference
<figure><img src="image_135.png" /><figcaption>FIGURE 14. (Colour online) Distributions of temperature (a) and vertical velocity (b) in a vertical cross-section at a position of 100 mm from the front wall and at time = 60 min in experiment 4 in the ‘dome’ regime.</figcaption></figure>

across the reservoir, which is fixed, to calculate *B* values and obtained threshold values of ≈0.6. Their result cannot be compared directly with the present ones, in part because we are using the volume-averaged temperatures of the two fluids and in part because temperatures are continuously changing during an experiment. Using volume-averaged temperatures, the threshold *B* value for the Rayleigh-Bénard set-up of Davaille *et al.* (2003) would hover around a value of approximately 1.2.

# 6. Time changes of thermal structure and the progression of mixing

## 6.1. *Temperature, velocity and compositional fields*

Figure 14 illustrates the relationship between the flow structure and the temperature field in the dome regime. The hottest regions are all made of lower fluid, where heat production is larger than average, and it would be difficult to distinguish this overall pattern from that of a Rayleigh-Bénard set-up with the same initial layering arrangement. In this particular case, two large domes have almost reached the top of the tank, where fluid is exposed to the cold boundary, and yet they do not recede. The vertical velocity distribution shows that they develop an internal circulation with a central upwelling surrounded by downwellings. The sharpest velocity gradients occur at their edges and are responsible for the tearing out of lower fluid into the upper one. Both the internal circulation inside the domes and filament entrainment of the lower layer are visible in figure 15. The most interesting feature of this form of convection is that the domes remain stable for long time intervals. Occasionally, the domes fall back down and fold over, encapsulating small blobs of denser upper fluid (figure 15e,f). In this regime, mixing proceeds through the two different mechanisms that have been described above, but each one is active in one fluid only. The lower fluid entrains upper fluid through folding and encapsulating, whereas the upper fluid tears thin slivers out of lower fluid. For geochemical applications, it is important to note that the lower fluid maintains a well-defined interface with the upper fluid even when it gets diluted to some extent.
<figure><img src="image_136.png" /><figcaption>FIGURE 15. (Colour online) Mixing in the dome regime (experiment 4) at a position of 100 mm from the front wall and time = 0 min (a), 25 min (b), 40 min (c), 55 min (d), 70 min (e), 85 min (f). Black contours correspond to a given pixel intensity value threshold and delineate the boundary between the two fluids as described in § 5. The large dome on the left develops an overhang on its right-hand side. Note the thin upwellings that emanate from the interface and that are made of diluted lower fluid.</figcaption></figure>

Vertical profiles of the horizontally averaged temperature illustrate the thermal evolution of the two fluids and its relation to the bulk thermal structure. Early profiles are consistent with the presence of two well-separated layers, with the lower
<figure><img src="image_137.png" /><figcaption>FIGURE 16. (Colour online) Averaged vertical temperature profiles at different times during experiment 4, which is the same experiment as in figure 14. Colours indicate time evolution from green to blue, red and black. The presence of a warmer lower layer ($h_1 = 15$ mm) is clearly visible in early stages.</figcaption></figure>

one heating up more rapidly than the upper one due to its larger heat production (figure 16). A temperature gradient develops over a large fluid thickness above the initial position of the lower layer, due to the large interface deformation. Without information on the initial structure of the tank, this could be interpreted as due to an unstable thermal boundary layer above a heated horizontal boundary. As hot domes of lower fluid rise, they raise temperatures in the upper fluid and eradicate the sharp temperature gradient that marked the presence of a distinctive lower layer. Late temperature profiles are characterized by a single well-developed thermal boundary layer at the top of an interior region. The upper part of this boundary layer barely changes whilst the temperature distribution below keeps evolving.

In the stratified regime, the interface between the two fluids does not deform by large amounts and does not deviate markedly from its initial position (figure 17). Most of the activity is in the upper fluid. In marked contrast to the experiments in the dome regime, vertical temperature profiles preserve a trace of the initial two-layer configuration for long times (figure 18). These profiles also show a gradual decrease of the thickness of the lower layer, which is indicated by a break in the profiles. In this regime, mixing proceeds exclusively by the entrainment of stretched filaments at the top of cusps of the lower fluid (figure 19).

## 6.2. *The progression of mixing between the two fluids*

We have found that the thermal structure evolves on a different time scale than that of mixing. As discussed above, all experiments end up with a steady-state thermal structure that is identical to that for a homogeneous reservoir. As discussed above, this does not imply that the two fluids are truly mixed in the chemical sense, such that there are no variations of composition at a local scale. In order to determine how mixing progresses, we have focussed on the lower fluid and measured two different variables. One is the volume of lower fluid that is separated from the upper fluid by a
<figure><img src="image_138.png" /></figure>

FIGURE 17. (Colour online) Distributions of temperature (a) and vertical velocity (b) in a vertical cross-section at time = 90 min in the stratified regime (experiment 7). The lower layer was 12 mm thick initially. This vertical cross-section was obtained at a distance of 100 mm from the front wall.

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

FIGURE 18. (Colour online) Averaged vertical temperature profiles at different times in the stratified regime (experiment 7). Colours indicate time evolution from green to blue, red and black. Temperatures remain high in a thin region at the base of the tank. Note that, at late times, the temperature profile barely changes in the upper part of the tank whilst the interior structure still evolves at a significant rate.

well-defined interface, as discussed above, and the other is the area occupied by lower fluid at the base of the tank. The latter allows us to track the gradual disappearance of a lower layer independently of volume changes, which is appropriate for the Earth because the starting volume of primordial material at the base of the mantle is ill constrained.
<figure><img src="image_140.png" /><figcaption>FIGURE 19. (Colour online) Entrainment in the stratified regime (experiment 7) at time = 0 min (a), 70 min (b), 90 min (c), 110 min (d), 130 min (e), 150 min (f). The interface deforms in a number of cusps separated by large basins in the lower fluid. Thin plumes rise from these cusps. This vertical cross-section was obtained at a distance of 100 mm from the front wall.</figcaption></figure>

The volume and area of lower fluid decrease continuously during an experiment (figure 20a,b). Mixing/mingling does not begin at time $t = 0$ when the first scans are carried out after the microwave emitter is switched on, because it takes a finite
<figure><img src="image_141.png" /><figcaption>FIGURE 20. Dome regime: evolution in time of the lower layer volume scaled by its initial value $V_0$ (a), of the occupied area at the lower boundary scaled by its initial value $A_0$ (b), of the volume average temperature of the lower layer $T_1$ (c), of the top layer $T_2$ (d), volume-average temperature contrast between the two layers $T_{12}$ (e) and of the whole volume $\langle T \rangle$ (f) for experiment 21; $t_0$ is the time corresponding to the onset of mixing and $t_{end}$ marks the end of the well-defined lower fluid reservoir.</figcaption></figure>

temperature difference across the layers to initiate convection. Thus, changes of volume or area of lower fluid become measurable after some time. Mixing/mingling may be considered complete when the lower fluid has been completely swept away from the base of the tank, corresponding to a zero area, and when the volume that is contained in individual parcels, or a deformed layer at the base of the tank, drops below detection level. Defining departure from an initial value or a hit at a target final value requires rather arbitrary thresholds and depends on measurement accuracy as well as on the time interval between two frames. We chose instead to use data in the intermediate range because these focus on the times when mixing is most active, as opposed to early and late phases when the data are either ramping up or winding down. These data are close to a linear trend, indicating that the mixing rate does not vary by large factors in the course of an experiment. Variations of the mixing rate
are due to long-period fluctuations of the interior structure, due for example to the
rise and fall cycle of lower fluid domes when they exist. We extrapolate the trends
in the volume and area data to zero, which provides an estimate of the mixing time,
noted $t_{end}$. As shown below, we found that the volume and area data lead to values
of $t_{end}$ that are nearly identical. The best-fit line also provides an estimate of the time
when mixing starts to have a detectable impact on volume and area, noted $t_o$.

In the dome regime (figure 20), the volume (figure 20a) and area (figure 20b) of
lower fluid begin to deviate from their initial values and go to zero at about the
same times. The onset of mixing at a significant rate is linked to breaks in the time
evolutions of the average temperature of the lower fluid $T_1$ (figure 20c) and of the
temperature difference between the two fluids $T_{12}$ (figure 20e). Both temperatures peak
at approximately the same time as mixing progresses. One should note that, as might
be expected, the temperature difference decreases to zero together with the volume
and area of lower fluid. Due to mixing, the amount of lower fluid diminishes steadily,
which favours thermal equilibrium with the upper fluid.

In the stratified regime (figure 21), there is a marked difference with the dome
regime early in an experiment. Mixing is effective before the cusp and basin
morphology reaches maximum amplitude. Due to the resilience of the lower layer,
the lower fluid sees its volume decrease even though it is spread over the whole base
of the reservoir. The area of lower fluid begins to decrease at a larger time but goes
to zero at about the same time as the volume. At this time the temperature difference
between the two fluids drops to zero rather abruptly because there is no detectable
lower fluid present.

6.3. *Changes of thermal structure*

We have tracked changes of thermal structure as a function of time using four
different temperature values: the average temperatures of the upper and lower fluids
$T_1$ and $T_2$, their difference $T_{12}$ and the volume-averaged temperature for the whole
tank $\langle T \rangle$ (figure 20c-f in the dome regime and figure 21c-f in the stratified regime).
The evolutions of all these temperatures are modulated by low-amplitude long-period
oscillations which reflect the changes of thermal structure and mixing rate that occur
in the tank.

In the dome regime, the average temperature of the lower fluid, whose volume is
decreasing with time as mixing progresses, reaches a peak value and then decreases as
the ever-shrinking parcels of lower fluid equilibrate with the colder surrounding fluid.
As the upper fluid ingests lower fluid, it gradually warms up, but sensible heat plays
a minor role: the slow rate of warming is due to the enhanced heat production. The
progression of mixing is illustrated clearly by a decreasing thermal contrast between
the two fluids. In this regime, mixing is effected by the two different processes that
have already been described. Both act to decrease the thermal contrast between the
two fluids in two opposite ways, the former (filament entrainment) by enhancing
heat production in the upper fluid and the latter (folding of lower layer) by lessening
heat production in the lower fluid. We found that the latter is consistently much less
effective than the former.

The temperature difference between the two layers, noted $T_{12}$, peaks at some
specific time in both regimes. The time of the thermal peak marks a change in the
impact of mixing on the system, such that the influence of the lower fluid on the
bulk thermal structure begins to wane. We have determined the magnitude of the peak
temperature difference $T_{12}$, noted $T_{12max}$, and the time at which it occurs, noted $\tau_{12max}$
<figure><img src="image_142.png" /><figcaption>FIGURE 21. Stratified regime: time evolution of the lower layer volume scaled by its initial value $V_0$ (a), of the occupied area at the lower boundary scaled by its initial value $A_0$ (b), of the volume average temperature of the lower layer $T_1$ (c), of the top layer $T_2$ (d), volume-average temperature contrast between the two layers $T_{12}$ (e) and of the whole volume $\langle T \rangle$ (f) for experiment 34; $t_0$ (vol, area) marks the time when mixing has a measurable effect on the volume and area, respectively, and $t_{end}$ marks the disappearance of a well-defined lower fluid reservoir.</figcaption></figure>

(table 2). The magnitude of the peak is well constrained but the time is not due to two factors. The peak may be quite flat and spread over a finite time interval in some cases. In addition, measurements are repeated at time intervals of 5 min because of the large number of scans that are made, which makes for a poor time resolution. For our present purposes, we do not need accurate determinations of the peak time and the experimental procedure was the result of a compromise between measurement repetition and data coverage. We have scaled $T_{12max}$ and $\tau_{12max}$ values to the reference
<figure><img src="image_143.png" /></figure>

FIGURE 22. Volume-average temperature contrast between the two layers $T_{12max}/\Delta T$ (circles) as a function of the Rayleigh–Roberts number. The line represents a power-law best fit with the exponent fixed at $-1/4$ and calculated with data for $Ra_H > 10^5$ (see also table 5).

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

FIGURE 23. Dependence of the time scale $\tau_{12max}/\tau_c$ of the volume-average temperature contrast on the Rayleigh–Roberts number. The line corresponds to a power-law fit with an exponent fixed at $-1/4$ (table 5).

conduction temperature difference $\Delta T$ and time scale $\tau_c$ (figures 22 and 23). The data are consistent with power-law relationships with the bulk Rayleigh–Roberts number $Ra_H$ and we have determined the corresponding coefficients using a best-fit procedure. The best-fit exponents are remarkably close to values for internally heated homogeneous layers (table 5).

The peak time data in figure 23 obviously do not provide strong evidence for a power-law relationship but what concerns us here is a comparison between this time
<table><caption>TABLE 5. Parameters of empirical best-fit power laws for the volume-average temperature difference between the two layers $T_{12max}$ and the time at which this maximum occurs $\tau_{12max}$ at high values of the Rayleigh–Roberts number ($Ra_H \ge 10^5$).</caption><thead><tr><th>Variable</th><th>Exponent left to vary</th><th>Exponent fixed to -1/4</th></tr></thead><tbody><tr><td>T<sub>12max</sub>/&Delta;T</td><td>4.86(&plusmn;4.14) Ra<sub>H</sub><sup>-0.200(&plusmn;0.043)</sup></td><td>10.05(&plusmn;0.93) Ra<sub>H</sub><sup>-1/4</sup></td></tr><tr><td>&tau;<sub>12max</sub>/&tau;<sub>c</sub></td><td>3.87(&plusmn;2.17) Ra<sub>H</sub><sup>-0.227(&plusmn;0.064)</sup></td><td>5.12(&plusmn;1.15) Ra<sub>H</sub><sup>-1/4</sup></td></tr></tbody></table>

and the thermal relaxation time for the reservoir, $\tau_e$ ($\tau_{12max} \approx 5.1 Ra_H^{-1/4}$ whereas $\tau_e = 3.0 Ra_H^{-1/4}$). Thus, at the time of the peak thermal contrast, $\langle T \rangle$ is within 18% of its final value and the bulk thermal structure of the reservoir is close to the final steady state, but mixing is far from completion. At that time, for example, the volume of lower fluid at the base of the tank is approximately 70% of its initial value in the dome regime and has barely changed in the stratified regime (figures 20 and 21). This shows that thermal steady state is achieved long before complete mingling or mixing.

We are interested in the large-scale flow pattern and mixing processes, which are intimately linked to one another and which are driven by thermally induced density differences. The peak temperature contrast between the two layers, noted $T_{12max}$, is associated with the active mixing phase and is therefore appropriate for calculating the buoyancy number. This buoyancy number, noted $B_{eff}$, is used in the regime diagram of figure 13(b).

### 6.4. Mixing time

Mixing proceeds mostly by the entrainment of thin slivers of one fluid into the other, which occurs at the interface between the two fluids. Thus, we can apply the same reasoning as Davaille (1999b), who derived a relationship for the entrainment rate of a denser layer based on scalings for the shear stress that is imparted on the interface between two convecting fluid layers:

$$t_{mix RB} = h_1 \frac{1}{C_{RB}} \frac{B_l^2}{\kappa} \left(1 + \frac{\gamma}{B_l}\right) \left(\frac{\kappa \mu_2}{\alpha \rho g \theta}\right)^{1/3}, \quad (6.1)$$

where $C_{RB}$ is a proportionality constant and $B_l$ is a buoyancy factor related to the temperature difference between the two layers $\theta$ ($B_l = \Delta\rho/\alpha\rho\theta$). The $(1 + \gamma/B_l)$ was added to account for the viscosity contrast effect in the experiments. The volume of lower fluid is small compared to that of the upper fluid, implying that mixing does not affect significantly the physical properties of the upper fluid. This is why Davaille (1999b) used the physical properties of the upper fluid.

Davaille (1999a) scaled the temperature difference between the layers $\theta$ to the thermal contrast across the tank, which remains constant in the Rayleigh–Bénard configuration of her experiments and determined a best-fit experimental coefficient $C_{RB} = 0.61 \pm 0.05$. We cannot follow the same line of reasoning because the total thermal contrast evolves with time and eventually settles to zero, which is only relevant for the final homogenized fluid. Following Fourel et al. (2017) we use again the maximum temperature difference between the two fluids:

$$t_{mix IH} = h_1 \frac{1}{C_{IH}} \frac{B_{eff}^2}{\kappa} \left(1 + \frac{\gamma}{B_{eff}}\right) \left(\frac{\kappa \mu_2}{\alpha \rho g T_{12max}}\right)^{1/3}. \quad (6.2)$$
<figure><img src="image_145.png" /><figcaption>FIGURE 24. Experimental versus theoretical mixing time in internally heated convection experiments: $t_{experimental}$ is equal to $t_{end}$, which is deduced from both the volume and area data at the lower boundary; $t_{mixIH}$ is calculated from (6.2) with an experimental constant $C_{IH} = 0.08$.</figcaption></figure>

The results are shown in figure 24. This scaling holds remarkably well over 2 orders of magnitude with an experimental constant of $C_{IH} = 0.08 \pm 0.01$. The data exhibit some scatter around the best-fit relationship but this scatter is related neither to the convection regime (i.e. dome or stratified) nor to the viscosity contrast ($\gamma$ value).

Two asymptotic expressions for the mixing time shed light on the viscosity control of mixing phenomena. If $\gamma \gg 1$:

$$t_{mixIH} \sim \gamma T_{12max}^{-(4/3)} \qquad (6.3)$$

If $\gamma \ll 1$, the mixing time does not depend on the viscosity of the lower fluid:

$$t_{mixIH} \sim T_{12max}^{-(7/3)} \qquad (6.4)$$

# 7. Discussion

## 7.1. Internal versus basal heating

In a homogeneous fluid layer, internal heating breaks the symmetry between top and bottom that exists in a Rayleigh-Bénard configuration, where heat is supplied through the base of the tank and extracted at the top. In the latter, convection can be described as sets of equally strong upwellings and downwellings issuing from thermal boundary layers at the top and bottom. In the internal heating configuration, there is only one unstable boundary layer at the top of the reservoir, implying that convection proceeds by strong cold downwellings and a diffuse upward return flow with weak thermal anomalies (Vilella et al. 2018). Adding a denser basal layer with higher than average heat production could be considered as akin to basal heating for the upper fluid, but this may be misleading. In many cases, the large interface deformation that occurs is responsible for vertical temperature variations over a thick region at the base of the tank (figure 3), which could be mistaken for a basal thermal boundary layer due to
heating from below. This lower region is such that heat transport is not dominated by conduction and is a transient feature, as shown by figure 3.

Our experiments show that, in a two-layer stratified configuration, internal heating leads to structures and flow patterns that are very similar to those of Rayleigh-Bénard set-ups. In the dome regime, the lower layer gets distorted and split into patches or piles with a range of shapes that are the same in both types of experimental conditions. Basal patches act to localize upwellings and plumes at their edges independently of the mode of heating. One important difference with a Rayleigh-Bénard set-up is that the large protrusions of lower fluid that grow into the upper fluid in the dome regime are remarkably steady and may persist over large fractions of the mixing time. The domes are made of intrinsically denser fluid, but they are heated internally at a larger rate than the upper fluid, which maintains a positive temperature anomaly. They are also supported by stresses imparted by adjacent upwellings in the upper fluid.

## 7.2. Earth's mantle heterogeneity

Geochemical and geophysical data indicate that the mantle is heterogeneous, but the origins and locations of the different mantle components remain debated. A major stumbling block is that the initial mantle composition, acquired when the Earth began to behave as an isolated system with no external input, is not known precisely (Jackson & Carlson 2012; Gale et al. 2013). In addition, the primordial mantle may well have been stratified. Today, after 4.5 billion years of evolution, the upper mantle is made of different components mixed in variable proportions and two extreme compositions have been determined, called for the sake of convenience 'enriched' and 'depleted' respectively. Here, 'enriched' is only a relative term defined with respect to 'depleted' (in other words it does not necessarily imply material that has effectively been enriched with respect to the starting composition). The enriched components are brought by mantle plumes coming from large depths.

Our experiments show that, starting from a two-layer reservoir, the upper fluid becomes a mixture of the two different fluids with no detectable large-scale changes of concentration with depth, save for residual patches of lower fluids at the base which slowly get eroded by plumes. These plumes are made of two components: upper fluid, which has become a mixture of the initial lower and upper fluids, and pure (or weakly diluted) initial lower fluid. Translated into geochemical parlance, these plumes are therefore enriched with respect to the surrounding fluid, but they are not made of a pure primordial end member. At any given time, they are not made with exactly the same proportions of the two components and hence do not have the same composition. This is consistent with the geochemical data (Gale et al. 2013).

This brief overview of geochemical constraints emphasizes the key contribution of mantle plumes, which sample deep mantle material. As the Earth has been undergoing convection for most of its lifetime, one must address two questions. One is whether or not primordial mantle material can survive and avoid being mixed into the bulk mantle for more than four billions years. This depends on many physical properties that are not known precisely, such as the intrinsic density of primordial mantle material. As the mixing time depends mostly on processes at the interface between the two materials and on the temperature difference between them independently of the mode of heating, the present study cannot add much to earlier analyses (Davaille et al. 2002). Using available constraints on physical properties and temperature differences in the Earth's mantle, Davaille et al. (2002) have concluded that primordial material can survive for many billions of years at the base of the mantle. Our results lead to the same
<figure><img src="image_146.png" /></figure>

FIGURE 25. (Colour online) Time variation of the volume (a) and the volume-averaged temperatures (b) of an enriched basal reservoir for $h_1 = 400$ km, $F = 3$, $Ra_H = 10^8$, $B_{eff} = 1.2$, $\gamma = 10$ and $\Delta\rho/\rho = 1\%$.

conclusion if the viscosity and density contrasts are large enough. Figure 25 shows how the amount of enriched material and its volume-averaged temperature evolve with time for parameter values that are representative of mantle conditions: $h_1 = 400$ km, $a = 0.13$, enrichment factor $F = 3$, $\mu_2 = 10^{21}$ Pa s, viscosity ratio $\gamma = 10$, intrinsic density contrast $\Delta\rho/\rho = 1\%$, $Ra_H = 10^8$ and a starting temperature of $1600\,^\circ$C. With $\alpha = 2 \times 10^{-5} \, \text{K}^{-1}$, the buoyancy number is $B_{eff} = 1.2$, implying a dome regime. In this example, the key difference with the Rayleigh-Bénard configuration is that the time and temperature scales do not depend on an arbitrarily fixed temperature difference between the top and bottom of the reservoir: they are deduced from the total amount and distribution of heat-producing elements. This calculation is just one in a set of many and is only meant to illustrate the time and temperature scales involved in mantle conditions. It shows a lower fluid temperature that initially rises and peaks after almost a billion years, and a lower fluid volume that has decreased by approximately 40% after 4.5 billion years. Interestingly, the average temperature difference between the lower and upper fluids stays quite constant at $\approx 400\,^\circ$C for several billion years.

A second question is whether or not a basal heat flux is required to generate deep plumes. This is associated with a subsidiary question: has the basal heat flux always been as strong as today in the past? Heating the mantle from below depends on the cooling of the core, which in turn depends on the secular evolution of the coupled core-mantle system. Thus, the answer to the subsidiary question requires a full thermal model of the Earth, which is outside the scope of this paper. We have shown, however, that deep mantle plumes with the required characteristics can be generated by enriched lower mantle material with no heat coming from the core.

### 7.3. The cooling of the Earth

The cooling of the Earth over geological time has long been a key issue and at the origin of an ancient and famous controversy between physicists and geologists (England, Molnar & Richter 2007). Today, the controversy has shifted to the mantle convection regime and its relation to the surface heat loss. The present-day cooling rate can be deduced from the global heat budget for the planet (Jaupart et al. 2015) and values over much of Earth's history are derived from the source temperatures of past basalt eruptions (Herzberg, Condie & Korenaga 2010). One expects that a
hotter mantle would be much less viscous and would convect more rapidly than
today, implying large rates of heat loss in the past. Physical models relying on
relationships between the surface heat flux and the interior temperature, such as
equation (4.2), involve the mantle viscosity, which depends strongly on temperature.
They predict that the rate of heat loss has decreased more rapidly than heat production
and underestimate the present-day heat loss (Labrosse & Jaupart 2007). In order to
correct this discrepancy, one must increase the adjustment time of mantle convection
to temperature changes. One solution is to invoke a stratified mantle (McKenzie
& Richter 1981) but, as shown in this study, this is only viable over a long time
interval if the intrinsic density contrast is large enough. The structure of a highly
heterogeneous lowermost mantle that has emerged in the last two decades is consistent
with the remnants of a primordial layer that has been extensively deformed by
convection, which is not likely to add much thermal inertia to the mantle as a whole.

7.4. *Temperature at the core–mantle boundary*

Thermal conditions at the core–mantle boundary, which dictate how the geodynamo operates, are strongly affected by the presence of a basal primordial layer. Numerical simulations of the coupled thermal evolution of Earth’s mantle and core demonstrate this and can only reproduce the present-day size of the inner core if such material is included (Nakagawa & Tackley 2005, 2014). These calculations, however, do not allow for enhanced heat generation in this material and have not been aimed at mixing phenomena. Our experiments have their own shortcomings because they do not include an analogue for the core, implying that they should be regarded as exploratory. They show that an enriched basal layer leads to an early temperature rise followed by a slow decline at the lower boundary. For application to the Earth, one must pay attention to the initial conditions but it is clear that anomalously high heat production in the lowermost mantle is bound to reduce the magnitude of core cooling significantly. Whether it acts to prevent core cooling altogether for some time is an intriguing possibility.

8. Conclusions

We have studied how convection in a two-layer heterogeneous internally heated reservoir evolves with time. The two fluids progressively mix with one another leading to a final homogenized reservoir, such that the heat flux at the top and the average interior temperature depend on the bulk Rayleigh–Roberts number as in a homogeneous fluid. Thermal steady state is reached before mixing has proceeded to completion. The temperature difference between the two fluids reaches a maximum during the mixing phase, as the volume of lower fluid is decreasing. At late times, the reservoir may be in an apparent thermal steady state but may still contain highly deformed patches of lower fluid at its base. These patches are associated with upwellings and plumes emanating from their peripheral regions which look like they are due to heating from below the reservoir.

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 and supported by the Tellus/PNP program of INSU-CNRS.
This is IPGP contribution number 4019.
**Supplementary movies and materials**

Supplementary movies and materials are available at
https://doi.org/10.1017/jfm.2019.243.

REFERENCES

BLOOMFIELD, V. A. & DEWAN, R. K. 1971 Viscosity of liquid mixtures. *J. Phys. Chem.* **75** (20), 3113–3119.

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

DAVAILLE, A. 1999a Simultaneous generation of hotspots and superswells by convection in a heterogeneous planetary mantle. *Nature* **402**, 756–760.

DAVAILLE, A. 1999b Two-layer thermal convection in miscible viscous fluids. *J. Fluid Mech.* **379**, 223–253.

DAVAILLE, A., GIRARD, F. & LE BARS, M. 2002 How to anchor hotspots in a convecting mantle? *Earth Planet. Sci. Lett.* **203** (2), 621–634.

DAVAILLE, A., LE BARS, M. & CARBONNE, C. 2003 Thermal convection in a heterogeneous mantle. *C. R. Geosci.* **335** (1), 141–156.

DAVAILLE, A. & LIMARE, A. 2015 7.03 - Laboratory studies of mantle convection. In *Treatise on Geophysics*, 2nd edn. (ed. G. Schubert), pp. 73–144. Elsevier.

DESCHAMPS, F., COBDEN, L. & TACKLEY, P. J. 2012 The primitive nature of large low shear-wave velocity provinces. *Earth Planet. Sci. Lett.* **349**, 198–208.

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

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

ENGLAND, P., MOLNAR, P. & RICHTER, F. 2007 John Perry’s neglected critique of Kelvin’s age for the Earth: a missed opportunity in geodynamics. *GSA Today* **17** (1), 4–9.

FOUREL, L., LIMARE, A., JAUPART, C., SURDUCAN, E., FARNETANI, C. G., KAMINSKI, E. C., NEAMTU, C. & SURDUCAN, V. 2017 The Earth’s mantle in a microwave oven: thermal convection driven by a heterogeneous distribution of heat sources. *Exp. Fluids* **58** (8), 90.

FRENCH, S. W. & ROMANOWICZ, B. 2015 Broad plumes rooted at the base of the Earth’s mantle beneath major hotspots. *Nature* **525**, 95–113.

GALE, A., DALTON, C. A., LANGMUIR, C. H., SU, Y. & SCHILLING, J.-G. 2013 The mean composition of ocean ridge basalts. *Geochem. Geophys. Geosyst.* **14** (3), 489–518.

GARNERO, E. J. & HELMBERGER, D. V. 1996 Seismic detection of a thin laterally varying boundary layer at the base of the mantle beneath the central-pacific. *Geophys. Res. Lett.* **23** (9), 977–980.

GARNERO, E. J., MCNAMARA, A. K. & SHIM, S.-H. 2016 Continent-sized anomalous zones with low seismic velocity at the base of Earth’s mantle. *Nature Geosci.* **9**, 481–489.

GOLUSKIN, D. 2016 *Internally Heated Convection and Rayleigh-Bénard Convection*. Springer.

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

VAN DER HILST, R. D. & KÁRASON, H. 1999 Compositional heterogeneity in the bottom 1000 kilometers of earth’s mantle: toward a hybrid convection model. *Science* **283** (5409), 1885–1888.

HOFMANN, A. W. 2003 Sampling mantle heterogeneity through oceanic basalts: isotopes and trace elements. *Treatise Geochem.* **2**, 568.

JACKSON, M. G. & CARLSON, R. W. 2012 Homogeneous superchondritic 142Nd/144Nd in the mid-ocean ridge basalt and ocean island basalt mantle. *Geochem. Geophys. Geosyst.* **13** (6), Q06011.
JAUPART, C., LABROSSE, S., LUCAZEAU, F. & MARESCHAL, J.-C. 2015 Temperatures, heat and energy in the mantle of the Earth. In *Treatise of Geophysics*, 2nd edn. (ed. G. Schubert), vol. 7, pp. 223–270. Elsevier.

JAVOY, M. & KAMINSKI, E. 2014 Earth's Uranium and Thorium content and geoneutrinos fluxes based on enstatite chondrites. *Earth Planet. Sci. Lett.* **407**, 1–8.

JELLINEK, A. M. & MANGA, M. 2002 The influence of a chemical boundary layer on the fixity, spacing and lifetime of mantle plumes. *Nature* **418**, 760–763.

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.

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

LABROSSE, S. & JAUPART, C. 2007 Thermal evolution of the earth: secular changes and fluctuations of plate characteristics. *Earth Planet. Sci. Lett.* **260** (3), 465–481.

LE BARS, M. & DAVAILLE, A. 2004a Large interface deformation in two-layer thermal convection of miscible viscous fluids. *J. Fluid Mech.* **499**, 75–110.

LE BARS, M. & DAVAILLE, A. 2004b Whole layer convection in a heterogeneous planetary mantle. *J. Geophys. Res.* **109**, B03403.

LEPOT, S., AUMAÎTRE, S. & GALLET, B. 2018 Radiative heating achieves the ultimate regime of thermal convection. *Proc. Natl Acad. Sci. USA* **115** (36), 8937–8941.

LI, M., MCNAMARA, A. K., GARNERO, E. J. & YU, S. 2017 Compositionally-distinct ultra-low velocity zones on earths core-mantle boundary. *Nature Commun.* **8** (1), 177.

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. & RICHTER, F. M. 1981 Parametrized thermal convection in a layered region and the thermal history of the Earth. *J. Geophys. Res.* **86** (B12), 11667–11680.

MONTELLI, R., NOLET, G., DAHLEN, F. A. & MASTERS, G. 2006 A catalogue of deep mantle plumes: new results from finite-frequency tomography. *Geochem. Geophys. Geosyst.* **7** (11), Q11007.

NAKAGAWA, T. & TACKLEY, P. J. 2005 Deep mantle heat flow and thermal evolution of the Earth's core in thermochemical multiphase models of mantle convection. *Geochem. Geophys. Geosyst.* **6** (8), Q08003.

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.

OLSON, P. 1984 An experimental approach to thermal convection in a two-layered mantle. *J. Geophys. Res.* **89** (B13), 11293–11301.

OLSON, P. & KINCAID, C. 1991 Experiments on the interaction of thermal convection and compositional layering at the base of the mantle. *J. Geophys. Res.* **96** (B3), 4347–4354.

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.

RICHTER, F. M. & MCKENZIE, D. P. 1981 On some consequences and possible causes of layered mantle convection. *J. Geophys. Res.* **86** (B7), 6133–6142.

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.

SCHUBERT, G., TURCOTTE, D. & OLSON, P. 2001 Mantle Convection in the Earth and Planets. Cambridge University Press.
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.

TORSVIK, T. H., SMETHURST, M. A., BURKE, K. & STEINBERGER, B. 2006 Large igneous provinces generated from the margins of the large low-velocity provinces in the deep mantle. *Geophys. J. Intl* **167** (3), 1447–1460.

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

TURCOTTE, D. L., PAUL, D. & WHITE, W. M. 2001 Thorium-uranium systematics require layered mantle convection. *J. Geophys. Res.* **106** (B3), 4265–4276.

TURNER, J. S. 1974 Double-diffusive phenomena. *Annu. Rev. Fluid Mech.* **6** (1), 37–54.

VILELLA, K., LIMARE, A., JAUPART, C., FARNETANI, C., FOUREL, L. & KAMINSKI, E. 2018 Fundamentals of laminar free convection in internally heated fluids at values of the Rayleigh-Roberts number up to $10^9$. *J. Fluid Mech.* **846**, 966–998.