Towards the ultimate regime in Rayleigh–Darcy convection

Abstract

Numerical simulations are used to probe Rayleigh–Darcy convection in fluid-saturated porous media towards the ultimate regime. The present three-dimensional dataset, up to Rayleigh–Darcy number $\textit {Ra}=8\times 10^4$, suggests that the appropriate scaling of the Nusselt number is $\textit {Nu}=0.0081\textit {Ra}+0.067\textit {Ra}^{0.61}$, fitting the computed data for $\textit {Ra}\gtrsim 10^3$. Extrapolation of current predictions to the ultimate linear regime yields the asymptotic law $\textit {Nu}=0.0081 \textit {Ra}$, about $16\,\%$ less than indicated in previous studies. Upon examination of the flow structures near the boundaries, we confirm previous indications of small flow cells hierarchically nesting into supercells, and we show evidence that the supercells at the boundary are the footprints of the megaplumes that dominate the interior part of the flow. The present findings pave the way for more accurate modelling of geophysical systems, with special reference to geological $\textrm {CO}_2$ sequestration.

1. Introduction

When a fluid-saturated porous layer undergoes bottom-up heating, it is prone to exhibit Rayleigh–Darcy convection (Horton & Rogers 1945; Lapwood 1948; Graham & Steen 1994). 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. 2012; Lepot, Aumaître & Gallet 2018; Zhu et al. 2018; Wang, Zhou & Sun 2020), 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 1962).

Classical theory suggests that the appropriate scaling in the ultimate regime of Rayleigh–Darcy convection should be linear (Malkus 1954; Howard 1964). 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. 2004; Wen et al. 2012; Hassanzadeh, Chini & Doering 2014). This prediction is supported by existing numerical studies, mainly based on two-dimensional computations (Otero et al. 2004; Hewitt, Neufeld & Lister 2012; Wen et al. 2012; Wen, Corson & Chini 2015; De Paoli, Zonta & Soldati 2016) 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 2014) 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 2020).

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. 2012; Huppert & Neufeld 2014; Riaz & Cinar 2014; Emami-Meybodi et al. 2015; De Paoli 2021). 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 2014; De Paoli et al. 2016; De Paoli, Zonta & Soldati 2017). 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

(2.1)\begin{equation} \phi\frac{\partial \theta^{*}}{\partial t^{*}}+\boldsymbol{\nabla} \boldsymbol{\cdot}(\boldsymbol{u}^{*}\theta^{*}-\phi D \boldsymbol{\nabla}\theta^{*})=0 { ,} \end{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,

(2.2)\begin{equation} \rho^{*}(\theta^{*})=\rho^{*}(\theta^{*}_{min})-\Delta\rho^{*}\frac{\theta^{*}-\theta^{*}_{min}}{\theta^{*}_{max}- \theta^{*}_{min}}{ ,} \end{equation}

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

(2.3a,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.

Figure 1. Sketch of the computational domain – a cube with side length $h^*$ – used to study Rayleigh–Darcy convection. The flow is heated at the bottom, $\theta ^*(y^*=0)=\theta ^*_{{max}}$, and cooled at the top, $\theta ^*(y^*=h^*)=\theta ^*_{{min}}$, and boundaries in the $x^*$ and $z^*$ directions are assumed to be periodic. The gravity acceleration ($\boldsymbol {g}$) points downwards. The temperature distribution $\theta ^*$ for the case $\textit {Ra}=1\times 10^4$ is also shown for illustrative purposes on the side boundaries and in a plane very close to the top boundary (i.e. at a distance of $50h^*/\textit {Ra}$ from the top boundary).

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

(2.4a–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.3a,b):

(2.5)\begin{gather} \frac{\partial \theta}{\partial t}+\boldsymbol{\nabla} \boldsymbol{\cdot} \left( {\boldsymbol{u}} \theta -\frac{1}{\textit{Ra}} \boldsymbol{\nabla} \theta\right)=0 ,\end{gather}

(2.6a,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

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

(2.7b)\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 2013; Wen et al. 2018). 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 2014), 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 2014).

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 2000), which has been extensively used for DNS of wall-bounded neutrally-buoyant and unstably-stratified turbulent flows (Pirozzoli 2014; Pirozzoli et al. 2017). 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,

(2.8)\begin{equation} \nabla^2 p = \frac{\partial \theta}{\partial y}, \end{equation}

with $\partial p / \partial y = 0$ at walls to satisfy the impermeability condition, which is then fed into (2.6a,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 1985; Orlandi 2000).

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. (2014).

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,

(3.1)\begin{equation} \textit{Nu}=-\left \langle \frac{1}{2}\int_{0}^{1}\int_{0}^{1}\left.\frac{\partial \theta}{\partial y}\right\vert_{y=0}+\left.\frac{\partial \theta}{\partial y}\right\vert_{y=1} \,\text{d}x\,\text{d}z\right\rangle, \end{equation}

with angle brackets indicating time averaging.

Table 1. Parameters of the main three-dimensional simulations performed in the present study. For each simulation, we report Rayleigh number $\textit {Ra}$, grid resolution $N_{x}\times N_{z}\times N_{y}$, Nusselt number $\textit {Nu}$, and number of samples $N_{s}$ over which $\textit {Nu}$ is averaged (see appendix B for further details). Additional three-dimensional simulations (not listed here) have been performed at $\textit {Ra}=2.5\times 10^3, 5\times 10^3$.

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. 2014). 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. (2014) 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.

Figure 2. Nusselt number as a function of the Rayleigh–Darcy number. Results obtained from the present numerical simulations are shown with filled circles, and the fitting curve, $\textit {Nu}=0.0081\textit {Ra}+0.067 \textit {Ra}^{0.61}$, is shown as a black solid line. Results obtained in previous simulations by Hewitt et al. (2014) (open triangles) and the extrapolated shifted linear scaling (blue dash-dotted line) are shown for comparison. The high-$\textit {Ra}$ portion of the $(\textit {Ra},\textit {Nu})$ parameter space covered by previous investigations is rendered by the light grey rectangle.

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. 2014) and two-dimensional (Hewitt et al. 2012; Wen et al. 2015; De Paoli et al. 2016) domains. According to the prediction of Hewitt et al. (2014), 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 2020). Data fitting then results in a functional representation of the type

(3.2)\begin{equation} \textit{Nu}/\textit{Ra}=0.0081+0.067\textit{Ra}^{-0.39} , \end{equation}

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

(3.3)\begin{equation} \textit{Nu} = 0.0081 \textit{Ra} , \end{equation}

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. (2014) 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. (2012), 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.

Figure 3. Compensated Nusselt number as a function of the Rayleigh–Darcy number, $\textit {Ra}$. Results obtained from the present numerical simulations, which we have run in three-dimensional and two-dimensional domains, are shown by filled circles ($\bullet$) and diamonds ($\blacklozenge$), respectively. For comparison purposes, a collection of previous data obtained in both two-dimensional domains (Hewitt et al. 2012; De Paoli et al. 2016; Wen et al. 2015) ($\square$, $\triangledown$ and $\triangleright$, respectively) and three-dimensional domains (Hewitt et al. 2014) ($\triangle$) is also shown, with open symbols. The black solid line denotes the best fit of our three-dimensional data, as from (3.2), and the blue solid line the shifted linear scaling $\textit {Nu}/\textit {Ra}=0.096+4.6/\textit {Ra}$ predicted by Hewitt et al. (2014). The scaling law $\textit {Nu}/\textit {Ra}=0.0069+2.75/\textit {Ra}$ proposed by Hewitt et al. (2012) for the two-dimensional case is shown with a solid red line. The asymptotic predictions for each scaling law are reported as dashed lines.

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. (2004), Hewitt et al. (2012), 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}$.

Figure 4. Temperature distribution in a plane close to the bottom boundary $(x,y=50/\textit {Ra},z)$, for simulations $S_8$ (a), $S_4$ (b), $S_2$ (c) and $S_1$ (d). The domain size is also explicitly indicated in diffusive–convective units in each panel. Indication of the domain size – rescaled using the diffusive–convective length scale $x^*_d$ – for the simulations at lower $\textit {Ra}$ is given in panel (a) by the white boxes. The size of the superstructures, estimated as $\lambda _i=2{\rm \pi} \textit {Ra}/k_i$ with $k_i/\textit {Ra}$ defined as in figure 5, is also shown in each panel (black bars).

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. 2013; Amooie, Soltanian & Moortgat 2018). 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. 2012; Wen et al. 2015; De Paoli et al. 2016), 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 4a) – 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. 2018; Krug, Lohse & Stevens 2020) or those of solar granules in the Sun's photosphere (Bahng & Schwarzschild 1961; Miesch 2005; Rieutord & Rincon 2010; Witze 2020) (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 4a).

Figure 5. Power spectra $k_{r}P(k_{r})$ (solid lines and symbols) of temperature distribution near the boundary (sampled at the location $y=50/\textit {Ra}$, as in figure 4), shown as a function of the radial wavenumber $k_{r}=\sqrt {k_{x}^{2}+k_{z}^{2}}$, for all $\textit {Ra}$ considered here. Results of the $S_3$ simulation are omitted for ease of reading. In the inset of the figure we plot the power spectra of the temperature distribution at the centre (midplane) of the domain. Labels $k_i/\textit {Ra}$ and $k'_i/\textit {Ra}$ indicate the peaks of the spectra. The corresponding wavenumber is explicitly reported in the main panel.

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. 2012; Hewitt, Neufeld & Lister 2013; De Paoli et al. 2016; Hewitt & Lister 2017).

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.