## 1. Introduction

When a fluid-saturated porous layer undergoes bottom-up heating, it is prone to exhibit Rayleigh–Darcy convection (Horton & Rogers Reference Horton and Rogers1945; Lapwood Reference Lapwood1948; Graham & Steen Reference Graham and Steen1994). The main controlling parameter is the Rayleigh–Darcy number $\textit {Ra}$, expressing the ratio of buoyancy to dissipative forces, and the resulting heat flux is characterized in terms of the Nusselt number, $\textit {Nu}$. The scaling of $\textit {Nu}$ with $\textit {Ra}$ in natural convection is traditionally described via power-law dependency, $\textit {Nu}\sim \textit {Ra}^{\alpha }$, although $\alpha$ may change as a result of the change in the structure of the convective flow. For instance, in Rayleigh–Bénard convection (He *et al.* Reference He, Funfschilling, Nobach, Bodenschatz and Ahlers2012; Lepot, Aumaître & Gallet Reference Lepot, Aumaître and Gallet2018; Zhu *et al.* Reference Zhu, Mathai, Stevens, Verzicco and Lohse2018; Wang, Zhou & Sun Reference Wang, Zhou and Sun2020), the general agreement is that $\alpha$ should gradually increase from about $1/3$ at moderate $\textit {Ra}$, to asymptotically approach a value of $1/2$, corresponding to the attainment of an ‘ultimate’ convection regime (Kraichnan Reference Kraichnan1962).

Classical theory suggests that the appropriate scaling in the ultimate regime of Rayleigh–Darcy convection should be linear (Malkus Reference Malkus1954; Howard Reference Howard1964). The main argument is that at high $\textit {Ra}$, the interior part of the flow is perfectly mixed, and temperature gradients are confined to near-wall boundary layers whose thickness is $\delta \propto \textit {Ra}^{-1}$ (whereas in Rayleigh–Bénard convection $\delta \propto \textit {Ra}^{-1/3}$); hence $\textit {Nu} \propto \delta ^{-1} \propto \textit {Ra}$. Other more rigorous arguments have also been provided based on variational calculus, which show that linear scaling yields maximum heat transfer among all suitably energy-constrained flows (Otero *et al.* Reference Otero, Dontcheva, Johnston, Worthing, Kurganov, Petrova and Doering2004; Wen *et al.* Reference Wen, Dianati, Lunasin, Chini and Doering2012; Hassanzadeh, Chini & Doering Reference Hassanzadeh, Chini and Doering2014). This prediction is supported by existing numerical studies, mainly based on two-dimensional computations (Otero *et al.* Reference Otero, Dontcheva, Johnston, Worthing, Kurganov, Petrova and Doering2004; Hewitt, Neufeld & Lister Reference Hewitt, Neufeld and Lister2012; Wen *et al.* Reference Wen, Dianati, Lunasin, Chini and Doering2012; Wen, Corson & Chini Reference Wen, Corson and Chini2015; De Paoli, Zonta & Soldati Reference De Paoli, Zonta and Soldati2016) carried out in the range $\textit {Ra} \leqslant 10^5$. Three-dimensional simulations carried out in recent times up to $\textit {Ra}=2\times 10^4$ (Hewitt, Neufeld & Lister Reference Hewitt, Neufeld and Lister2014) also seem to support the establishment of a shifted linear variation of $\textit {Nu}$ with $\textit {Ra}$. However, given the limited range of $\textit {Ra}$ in which data are available, uncertainties remain regarding the actual establishment of the expected ultimate convection regime, and especially how it is reached starting from finite $\textit {Ra}$ (Hewitt Reference Hewitt2020).

The quest for the $\textit {Nu}$ scaling in the high-$\textit {Ra}$ regime is crucial from a fundamental and also from an applied viewpoint, as the Rayleigh–Darcy model closely mimics flows in a number of geophysical applications, with special reference to geological $\textrm {CO}_2$ sequestration in deep saline aquifers (Hidalgo *et al.* Reference Hidalgo, Fe, Cueto-Felgueroso and Juanes2012; Huppert & Neufeld Reference Huppert and Neufeld2014; Riaz & Cinar Reference Riaz and Cinar2014; Emami-Meybodi *et al.* Reference Emami-Meybodi, Hassanzadeh, Green and Ennis-King2015; De Paoli Reference De Paoli2021). Real-life instances may exhibit $\textit {Ra}$ up to ${\sim }{ {O}}(10^5 \div 10^6)$, and small deviations from the expected ultimate linear scaling can produce large differences in the overall predicted transfer flux and in the corresponding cumulative time integral (Slim Reference Slim2014; De Paoli *et al.* Reference De Paoli, Zonta and Soldati2016; De Paoli, Zonta & Soldati Reference De Paoli, Zonta and Soldati2017). The goal of the present work is to investigate the high-$\textit {Ra}$ range of Rayleigh–Darcy convection, using a database of three-dimensional numerical simulations up to $\textit {Ra}=8\times 10^4$. Our aim is to extract reliable trends for the Nusselt number in a wide range of $\textit {Ra}$, and possibly to infer robust asymptotic estimates for the ultimate linear regime.

## 2. Methodology

We consider a fluid-saturated porous medium in a three-dimensional domain (see figure 1) with uniform porosity $\phi$ and uniform permeability $\kappa$. The flow, which is incompressible and governed by the Darcy law, is characterized by an unstable density difference ($\Delta \rho ^{*}$) induced by heating the flow from the bottom and cooling it from the top. Indicating by $x^*,z^*$ the horizontal directions, by $y^*$ the vertical direction (along which gravity acceleration $\boldsymbol {g}$ is directed) and by $\theta ^{*}$ the fluid temperature, the physical instance we consider corresponds to a fluid-saturated cubic volume with side length $h^{*}$ heated at the bottom, $\theta ^*(y^{*}=0)=\theta ^{*}_{max}$, and cooled at the top, $\theta ^*(y^{*}=h^{*})=\theta ^{*}_{min}$. The evolution of the temperature field is controlled by the advection–diffusion equation

where $t^{*}$ is time, $\boldsymbol {u}^{*}=(u^{*},v^{*},w^{*})$ is the volume-averaged velocity field and $D$ is the thermal diffusivity, which is considered constant here. The superscript $^{*}$ is used to indicate dimensional variables. We assume that fluid density, $\rho ^{*}$, is a linear function of temperature,

with $\Delta \rho ^{*}=\rho ^{*}(\theta ^{*}_{min})-\rho ^{*}(\theta ^{*}_{max})$. Assuming validity of the Boussinesq approximation (Landman & Schotting Reference Landman and Schotting2007; Zonta & Soldati Reference Zonta and Soldati2018), the flow field is fully described by the continuity and Darcy equations

*a*,

*b*)\begin{equation} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}^{*}=0,\quad \boldsymbol{u}^{*}=-\frac{\kappa}{\mu}\left(\boldsymbol{\nabla}P^{*}+\rho^{*}g \boldsymbol{j}\right) { ,} \end{equation}

with $\mu$ the fluid viscosity (constant), $P^{*}$ the pressure and $\boldsymbol {j}$ the vertical unit vector. The walls are assumed to be impermeable and isothermal, and periodicity is assumed in the wall-parallel directions.

### 2.1. Dimensionless equations

Natural velocity and length scales for the system are the buoyancy velocity, $V^{*}=g \Delta \rho ^{*} \kappa / \mu$, and the domain height, $h^{*}$, respectively. Using the set of dimensionless variables

*a*–

*c*)\begin{equation} \theta=\frac{\theta^{*}-\theta^{*}_{min}}{\theta^{*}_{max}-\theta^{*}_{min}} ,\quad t=\frac{t^{*}}{\phi h^{*}/V^{*}},\quad P=\frac{P^{*}}{\Delta \rho^{*}gh^{*}} , \end{equation}

and introducing the reduced pressure $p^{*}$, we obtain the dimensionless form of the governing equations (2.1), (2.3*a*,*b*):

*a*,

*b*)\begin{gather} \boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0,\quad {\boldsymbol{u}}=-\left( \boldsymbol{\nabla} p -\theta \boldsymbol{j}\right) ,\end{gather}

where $\textit {Ra}=g \Delta \rho ^{*} \kappa h^{*} / (\phi D \mu )=V^{*} h^{*} / (\phi D)$ is the Rayleigh–Darcy number. The wall boundary conditions for velocity and temperature then read as

*a*)\begin{gather} v(y=0)=0,\quad \theta(y=0)=1, \end{gather}

*b*)\begin{gather}v(y=1)=0,\quad \theta(y=1)=0. \end{gather}

As previously mentioned, for the physical system investigated here, the buoyancy velocity $V^{*}$ is a natural reference velocity scale (Fu, Cueto-Felgueroso & Juanes Reference Fu, Cueto-Felgueroso and Juanes2013; Wen *et al.* Reference Wen, Akhbari, Zhang and Hesse2018). At the same time, a possible reference length scale is the thickness of the porous layer, ${x}_{c}^{*}=h^*$ (convective scaling). However, an alternative choice for the reference length scale, which is perhaps more related to the physics of the phenomena under investigation, is ${x}_{d}^{*}=\phi D/V^{*}$ (diffusive–convective scaling). Note that ${x}^{*}_{d}$ represents the length over which advection and diffusion balance (Slim Reference Slim2014), and is independent of the physical domain thickness. When rescaled in the latter way, dimensions are bounded within the range ${x}^{*}/{x}_{d}^{*}\in [0,\textit {Ra}]$, and comparison between simulations at different $\textit {Ra}$ is easier. For this reason, lengths in this paper are rescaled with respect to $x_{d}^{*}$. Furthermore, introduction of this length scale also yields another interpretation of the Rayleigh–Darcy number, $\textit {Ra}={x}_{c}^{*}/{x}_{d}^{*}$, which may be regarded as the dimensionless height of the domain (Slim Reference Slim2014).

### 2.2. Computational details

The numerical simulations rely on the modified version of a second-order finite-difference incompressible flow solver, based on staggered arrangement of the flow variables (Orlandi Reference Orlandi2000), which has been extensively used for DNS of wall-bounded neutrally-buoyant and unstably-stratified turbulent flows (Pirozzoli Reference Pirozzoli2014; Pirozzoli *et al.* Reference Pirozzoli, Bernardini, Verzicco and Orlandi2017). The temperature transport equation is advanced in time by means of a hybrid third-order low-storage Runge–Kutta algorithm, whereby the convective terms are handled explicitly and the diffusive terms are handled implicitly, limited to the wall-normal direction. This approach guarantees that the total temperature variance is discretely preserved in the limit of inviscid flow. The pressure field in the forced Darcy equation is determined by solving the Poisson equation resulting from the divergence-free constraint,

with $\partial p / \partial y = 0$ at walls to satisfy the impermeability condition, which is then fed into (2.6*a*,*b*). Fourier expansion along the horizontal periodic directions yields a system of tridiagonal equations in the wall-normal direction for each Fourier mode, which is then solved with standard highly efficient numerical techniques (Kim & Moin Reference Kim and Moin1985; Orlandi Reference Orlandi2000).

The mesh spacing in the wall-parallel directions was decided based on preliminary grid-resolution studies at low $\textit {Ra}$ and inspection of the temperature spectra, to prevent any energy pile-up at the smallest resolved flow scales. Regarding the resolution in the wall-normal direction, we have followed the criterion that twenty points should be placed within the thermal boundary layer edge, identified through the peak location of the temperature variance, and grid points are clustered towards the walls using a hyperbolic tangent stretching function. Given the expected linear growth of the temperature gradients, the number of points in each coordinate direction was increased proportionally to $\textit {Ra}$. The time step is selected so that the CFL number is approximately unity for all the simulations herein reported. Preliminary calculations, carried out at $\textit {Ra} \leqslant 5\times 10^3$, have shown excellent agreement with the results of Hewitt *et al.* (Reference Hewitt, Neufeld and Lister2014).

## 3. Results

### 3.1. Heat transfer

We have carried out numerical simulations in the Rayleigh number range $2.5\times 10^3<\textit {Ra}<8\times 10^4$. Note that, because of the time step restrictions, the simulation at $\textit {Ra}=8\times 10^4$ – using $6144\times 6144\times 2048$ grid points – is almost 4096 times more computationally expensive than the simulation at $\textit {Ra}=1\times 10^4$ – using $768\times 768\times 256$ grid points – hence the computational challenge is substantial. A summary of the main computational parameters is provided in table 1. The heat flux is quantified in terms of the Nusselt number $\textit {Nu}$, evaluated as the mean temperature gradient at the top and bottom boundaries,

with angle brackets indicating time averaging.

In figure 2 we show the computed values of $\textit {Nu}$ as a function of $\textit {Ra}$ (filled circles), compared with data available in the literature (open triangles, Hewitt *et al.* Reference Hewitt, Neufeld and Lister2014). Also shown in light grey is the high-$\textit {Ra}$ portion of the $(\textit {Ra},\textit {Nu})$ parameter space covered in previous three-dimensional numerical simulations. The dash-dot line in figure 2 corresponds to the prediction of Hewitt *et al.* (Reference Hewitt, Neufeld and Lister2014) based on fitting of lower-$\textit {Ra}$ data ($\textit {Ra} \leqslant 2\times 10^4$), and amounting to a shifted linear variation, $\textit {Nu}=0.0096\textit {Ra} + 4.6$. As seen from figure 2, this fit overestimates our data by about $9\,\%$ at the highest $\textit {Ra}$; hence we believe that a different expression may be appropriate to more accurately characterize the trend towards the ultimate linear scaling.

In figure 3 we thus show the compensated Nusselt number $\textit {Nu}/\textit {Ra}$, for a collection of current and previous numerical data obtained in three-dimensional (Hewitt *et al.* Reference Hewitt, Neufeld and Lister2014) and two-dimensional (Hewitt *et al.* Reference Hewitt, Neufeld and Lister2012; Wen *et al.* Reference Wen, Corson and Chini2015; De Paoli *et al.* Reference De Paoli, Zonta and Soldati2016) domains. According to the prediction of Hewitt *et al.* (Reference Hewitt, Neufeld and Lister2014), the asymptotic value $\textit {Nu}/\textit {Ra}=0.0096$ should be reached for $\textit {Ra} \gtrsim 10^5$. Our data (filled circles) for $\textit {Ra}>2\times 10^4$ fall well below that prediction, showing no clear evidence for inception of an ultimate regime. An improved prediction for the pre-asymptotic regime can be constructed by fitting the numerical data and relying on the compelling evidence that the asymptotic scaling should be linear (Hewitt Reference Hewitt2020). Data fitting then results in a functional representation of the type

which amounts to linear dependence of $\textit {Nu}$ over $\textit {Ra}$ augmented with a sublinear corrective term, shown as a black solid line in figure 3. Compared with previous estimates, this correlation yields good prediction of $\textit {Nu}$ over a much larger range of $\textit {Ra}$, starting at $\textit {Ra} \sim 2.5\times 10^3$, with error no larger than $1.5\,\%$. Extrapolation to higher $\textit {Ra}$ than reached in the present computational campaign would yield the asymptotic scaling

in the ultimate regime. Based on our data, we expect that to have deviations of less than $5\,\%$ from the ultimate regime, $\textit {Ra}$ in excess of $5\times 10^5$ is needed. At those huge Rayleigh numbers, which closely represent the case of geological $\textrm {CO}_{2}$ sequestration in deep saline aquifers, differences between the current prediction and previous predictions of Hewitt *et al.* (Reference Hewitt, Neufeld and Lister2014) would be as large as ${\sim }16\,\%$. These findings are contrasted in figure 3 with the case of two-dimensional flow, for which the data generally support $\textit {Nu} \approx 0.0069 \textit {Ra} + 2.75$, as suggested by Hewitt *et al.* (Reference Hewitt, Neufeld and Lister2012), and attainment of a linear trend at $\textit {Ra} \sim 3\times 10^4$. Hence, extrapolation of our results yields the prediction that in the ultimate regime the Nusselt number should be larger by about $18\,\%$ in the three-dimensional than in the two-dimensional case.

### 3.2. Flow structure

It is reasonable to expect that the heat transfer behaviour just discussed is reflected in the behaviour of the flow structure. We characterize the flow structure focusing on its signature in the near-boundary region where the contribution of both diffusion-driven and buoyancy-driven mechanisms is important. An estimate of the location at which diffusion and buoyancy are both important is given by the temperature boundary layer thickness, $\delta$, which, in Rayleigh–Darcy convection, is $\delta \sim 30/\textit {Ra}$ (see Otero *et al.* (Reference Otero, Dontcheva, Johnston, Worthing, Kurganov, Petrova and Doering2004), Hewitt *et al.* (Reference Hewitt, Neufeld and Lister2012), and also appendix A of the present manuscript). Accordingly, in figure 4 we show the temperature distributions in a horizontal plane at $y=50/\textit {Ra}$ from the bottom boundary. Indication of the domain size in diffusive–convective units, $x^*_d$, for the simulations at lower $\textit {Ra}$ is given by the white boxes in panel (*a*). The complex flow pattern near the boundary, which appears to be still developing at $\textit {Ra}<2\times 10^4$, seems to attain a roughly self-similar structure at larger $\textit {Ra}$. Hence, in the following we will mostly refer to the $S_8$ simulation, since the flow features we wish to discuss are emphasized at that value of $\textit {Ra}$.

In figure 4 (and in particular in panel (*a*)) we clearly observe the presence of thin bright filaments which appear to be connected in a quasi-regular pattern of cells (Fu *et al.* Reference Fu, Cueto-Felgueroso and Juanes2013; Amooie, Soltanian & Moortgat Reference Amooie, Soltanian and Moortgat2018). These filaments correspond to regions where high-temperature fluid protrudes from the boundary, and are produced by tiny instabilities of the thermal boundary layer. Those are the three-dimensional counterpart of the protoplumes observed in previous two-dimensional studies (Hewitt *et al.* Reference Hewitt, Neufeld and Lister2012; Wen *et al.* Reference Wen, Corson and Chini2015; De Paoli *et al.* Reference De Paoli, Zonta and Soldati2016), and define the perimeter of rather homogeneous (dark) regions produced by downwelling motions (downwashes) of cold fluid towards the bottom boundary (see also the vertical slices shown in appendix A and comments therein). After impinging on the boundary, these downwashes spread horizontally and interact with the rising protoplumes, determining a dynamic pattern, in which some of the protoplumes eventually cluster into specific regions – thicker bright ridges (see figure 4*a*) – that define the boundaries of seemingly regular superstructures. These ‘supercells’, which encapsulate larger ensembles of smaller cells, have shape and dynamics recalling those of classical Rayleigh–Bénard turbulence (Stevens *et al.* Reference Stevens, Blass, Zhu, Verzicco and Lohse2018; Krug, Lohse & Stevens Reference Krug, Lohse and Stevens2020) or those of solar granules in the Sun's photosphere (Bahng & Schwarzschild Reference Bahng and Schwarzschild1961; Miesch Reference Miesch2005; Rieutord & Rincon Reference Rieutord and Rincon2010; Witze Reference Witze2020) (see also animations in the supplementary movies available at https://doi.org/10.1017/jfm.2020.1178 for a time-resolved rendering of the supercells). The size of the supercells $\lambda _i$, as inferred from the analysis of the wavenumber spectra (see figure 5 and the corresponding discussion), can be appreciated in figure 4. Since at large $\textit {Ra}$ the size of the supercells, which is approximately one-tenth of the box side length for simulation $S_8$ – see figure 4(*a*) – does not seem to evolve once the simulation has reached the steady state, we are keen to attribute to them an important role in the overall heat transfer mechanisms. At the same time, an important role is played by the smaller cells enclosed by supercells. Previous qualitative observations show a complex network of structures composed of small cells hierarchically nesting into supercells. To determine the size of such supercells, we use the time-averaged wavenumber power spectrum of the temperature signal, $P(k)$, which we obtain by taking the two-dimensional Fourier transform of the temperature field in a horizontal plane located near the boundary (such as the one sketched in figure 4*a*).

In figure 5, we show $k_rP(k_r)$ as a function of the normalized radial wavenumber $k_r/\textit {Ra}$, where $k_r=\sqrt {k_x^2+k_z^2}$, and $k_x$ and $k_z$ are the wavenumbers along $x$ and $z$, respectively. Results of the $S_3$ simulation are not shown, for better readability. For all considered values of $\textit {Ra}$, a sharp peak occurs at low wavenumbers, which clearly characterizes the size of the supercells. The wavenumber at which this peak occurs, labelled with $k_i/\textit {Ra}$ in figure 5 (where $i$ refers to the simulation number as in table 1, i.e. $i=1, 2, 4, 8$), decreases for increasing $\textit {Ra}$, thus indicating growth of the supercells with $\textit {Ra}$. Measurements of the temperature power spectrum at the midplane of the domain (inset of figure 5) is also characterized by the presence of sharp peaks (labelled as $k'_i/\textit {Ra}$), which in this case are due to large megaplumes developing in the interior part of the domain (Hewitt *et al.* Reference Hewitt, Neufeld and Lister2012; Hewitt, Neufeld & Lister Reference Hewitt, Neufeld and Lister2013; De Paoli *et al.* Reference De Paoli, Zonta and Soldati2016; Hewitt & Lister Reference Hewitt and Lister2017).

Direct comparison with the spectra near the boundary shows close coincidence of the spectral peaks, thus clearly indicating that supercells at the boundary are the footprint of the megaplumes that dominate the interior part of the flow. At larger $k_{r}$, the spectrum near the boundary features a local minimum, followed by rapid increase. Such increase, which is not present in the spectrum at the domain centre, is naturally associated with the small-scale protoplumes populating the region near the boundary (see the small filaments in figure 4).

## 4. Conclusions

We have reported results of numerical simulation of Rayleigh–Darcy convection inside a cubic porous layer up to Rayleigh–Darcy number $\textit {Ra}=8\times 10^4$. Exploiting the numerical dataset we show that the Nusselt number variation at sufficiently high $\textit {Ra}$ can be approximately expressed as $\textit {Nu}=0.0081\textit {Ra}+0.067\textit {Ra}^{0.61}$, with error bar no larger than $1.5\,\%$. Extrapolation of this prediction to the ultimate regime, which we expect to set in at $\textit {Ra} \gtrsim 5\times 10^5$, yields the asymptotic scaling $\textit {Nu}=0.0081\textit {Ra}$ – almost $16\,\%$ less than inferred in previous studies, and exceeding by about $18\,\%$ the Nusselt number found in two-dimensional simulations. We have further characterized the complex flow structure near the boundaries, which consists of small flow cells hierarchically nesting into supercells, and we have shown evidence that the supercells at the boundary are the footprint of megaplumes dominating the interior part of the flow. We believe that it would be possible to conduct direct investigations of the ultimate convective regime to check our extrapolations, by running a simulation at $\textit {Ra} \simeq 5\times 10^5$. However, we estimate a cost for this simulation of approximately 2 billion CPU hours, well beyond present-day capabilities.

## Supplementary movies

Supplementary movies are available at https://doi.org/10.1017/jfm.2020.1178.

## Acknowledgements

We are grateful to the CINECA supercomputing centre (under the grants ISCRA IsB20–3DSIMCON and PRACE Pra19-4938) and the Vienna Scientific Cluster (VSC) for generous allowance of computational resources. We also acknowledge Dr D. Hewitt and Dr B. Wen for providing data from Hewitt *et al.* (Reference Hewitt, Neufeld and Lister2012, Reference Hewitt, Neufeld and Lister2014) and Wen *et al.* (Reference Wen, Corson and Chini2015), respectively.

## Funding

S.P. and F.Z. gratefully acknowledge financial support from the ERASMUS+ project 29415-EPP–1–2014–1–IT–EPPKA3–ECHE between TU Wien and La Sapienza University. The authors acknowledge financial support from the TU Wien University Library through its Open Access Funding Program.

## Declaration of interests

The authors report no conflict of interest.

## Appendix A. Structure of the flow in a vertical plane

In figure 6 we show the temperature distribution in a vertical $(z,y)$ plane located at $x=1/2$. A close-up view of the temperature field near the boundary is also offered to exhibit the thickness of the thermal boundary layer, which, scaling as $\delta \sim 1/\textit {Ra}$, would otherwise be too small to be observed. Consistently with previous two- and three-dimensional studies (Hewitt *et al.* Reference Hewitt, Neufeld and Lister2012, Reference Hewitt, Neufeld and Lister2014; Wen *et al.* Reference Wen, Corson and Chini2015; De Paoli *et al.* Reference De Paoli, Zonta and Soldati2016), we observe that small fingers of light fluid emerge from the bottom boundary and move upwards, and correspondingly small fingers of heavy fluid emerge from the top boundary and move downwards. The small fingers then merge to form megaplumes (large columnar structures that dominate the core region of the flow) which, under the vigorous action of buoyancy, increase their vertical velocity and reach the opposite boundary. Upon impact with the boundary, the megaplumes are deflected and create the complex flow patterns analysed in detail in the main body of the present manuscript. As expected, the strength and persistence of the megaplumes increase with increasing $\textit {Ra}$.

## Appendix B. Evaluation of the Nusselt number

To estimate the averaged Nusselt number, we first run each simulation (starting from a still fluid, plus random perturbations) until the initial transient is finished. From that time, labelled as $t_{ss}$, the Nusselt number oscillates around a statistically steady value. In figure 7, we show the time evolution of the Nusselt number, $\textit {Nu}(t)$, over a specific time window after $t_{ss}$, for the main simulations herein performed. $\textit {Nu}(t)$ is observed to oscillate within $1\,\%$ to $3\,\%$ of the mean value (indicated by the dashed line). In particular, larger values of $\textit {Ra}$ correspond to lower relative amplitude of the fluctuations. The time-averaged Nusselt number (reported in table 1) is computed over a time window corresponding to 128, 64, 16, 8 and 4 dimensionless time units for $S_1$, $S_2$, $S_3$, $S_4$ and $S_8$, respectively. Naturally, the duration of the time window, $\Delta t=t-t_{ss}$, over which statistics are evaluated is less at the largest $\textit {Ra}$. However, it is important to point out that the Nusselt number is not only averaged in time, but also in space (along $x$ and $z$) over the entire wall surface, which increases the size of the ensemble over which statistics are evaluated. The effective number of samples, $N_{s}$, used to estimate the Nusselt number (and also reported in table 1) is here defined as $N_{s}=N_x\times N_z\times N_T$, where $N_T$ is the time interval (expressed in dimensionless convective units, i.e. turnover times) over which statistics are performed (starting from $t_{ss}$). Once the averaged value of the Nusselt number is estimated, we determine the best-fit scaling as in (3.2), yielding prediction error of no more than $1.5\,\%$, for $\textit {Ra} \geqslant 2.5 \times 10^3$.