Towards the ultimate regime in Rayleigh–Darcy convection

Numerical simulations are used to probe Rayleigh–Darcy convection in ﬂuid-saturated porous media towards the ultimate regime. The present three-dimensional dataset, up to Rayleigh–Darcy number Ra = 8 × 10 4 , suggests that the appropriate scaling of the Nusselt number is Nu = 0 . 0081 Ra + 0 . 067 Ra 0 . 61 , ﬁtting the computed data for Ra (cid:2) 10 3 . Extrapolation of current predictions to the ultimate linear regime yields the asymptotic law Nu = 0 . 0081 Ra , about 16 % less than indicated in previous studies. Upon examination of the ﬂow structures near the boundaries, we conﬁrm previous indications of small ﬂow 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 ﬂow. The present ﬁndings pave the way for more accurate modelling of geophysical systems, with special reference to geological CO 2 sequestration.


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 Ra, expressing the ratio of buoyancy to dissipative forces, and the resulting heat flux is characterized in terms of the Nusselt number, Nu.The scaling of Nu with Ra in natural convection is traditionally described via power-law dependency, Nu ∼ Ra α , although α 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 α should gradually increase from about 1/3 at moderate 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 Ra, the interior part of the flow is perfectly mixed, and temperature gradients are confined to near-wall boundary layers whose thickness is δ ∝ Ra −1 (whereas in Rayleigh-Bénard convection δ ∝ Ra −1/3 ); hence Nu ∝ δ −1 ∝ 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 Ra 10 5 .Three-dimensional simulations carried out in recent times up to Ra = 2 × 10 4 (Hewitt, Neufeld & Lister 2014) also seem to support the establishment of a shifted linear variation of Nu with Ra.However, given the limited range of 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 Ra (Hewitt 2020).
The quest for the Nu scaling in the high-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 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 Ra up to ∼O(10 5 ÷ 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-Ra range of Rayleigh-Darcy convection, using a database of three-dimensional numerical simulations up to Ra = 8 × 10 4 .Our aim is to extract reliable trends for the Nusselt number in a wide range of Ra, and possibly to infer robust asymptotic estimates for the ultimate linear regime.

Methodology
We consider a fluid-saturated porous medium in a three-dimensional domain (see figure 1) with uniform porosity φ and uniform permeability κ.The flow, which is incompressible and governed by the Darcy law, is characterized by an unstable density difference ( ρ * ) 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 g is directed) and by θ * the fluid temperature, the physical instance we consider corresponds to a fluid-saturated cubic volume with side length h * heated at the bottom, θ * ( y * = 0) = θ * max , and cooled at the top, θ * ( y * = h * ) = θ * min .The evolution of the temperature field is controlled by the advection-diffusion equation where t * is time, 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, ρ * , is a linear function of temperature, with ρ * = ρ * (θ * min ) − ρ * (θ * 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 with μ the fluid viscosity (constant), P * the pressure and 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 ρ * κ/μ, and the domain height, h * , respectively.Using the set of dimensionless variables and introducing the reduced pressure p * , we obtain the dimensionless form of the governing equations (2.1), (2.3a,b): where Ra = g ρ * κh * /(φDμ) = V * h * /(φD) is the Rayleigh-Darcy number.The wall boundary conditions for velocity and temperature then read as v( y = 0) = 0, θ(y = 0) = 1, (2.7a) 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 = φ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 and comparison between simulations at different 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, Ra = x * c /x * d , which may be regarded as the dimensionless height of the domain (Slim 2014).

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, with ∂p/∂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 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 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 Ra 5 × 10 3 , have shown excellent agreement with the results of Hewitt et al. (2014).

Heat transfer
We have carried out numerical simulations in the Rayleigh number range 2.5 × 10 3 < Ra < 8 × 10 4 .Note that, because of the time step restrictions, the simulation at Ra = 8 × 10 4 -using 6144 × 6144 × 2048 grid points -is almost 4096 times more computationally expensive than the simulation at Ra = 1 × 10 4 -using 768 × 768 × 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 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 Nu as a function of Ra (filled circles), compared with data available in the literature (open triangles, Hewitt et al. 2014).Also shown in light grey is the high-Ra portion of the (Ra, 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-Ra data (Ra 2 × 10 4 ), and amounting to a shifted linear variation, Nu = 0.0096Ra + 4.6.As seen from figure 2, this fit overestimates our data by about 9 % at the highest 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 Nu/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 Nu/Ra = 0.0096 should be reached for Ra 10 5 .Our data (filled circles) for Ra > 2 × 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 Nu/Ra = 0.0081 + 0.067Ra Towards the ultimate regime in Rayleigh-Darcy convection correlation yields good prediction of Nu over a much larger range of Ra, starting at Ra ∼ 2.5 × 10 3 , with error no larger than 1.5 %.Extrapolation to higher Ra than reached in the present computational campaign would yield the asymptotic scaling Nu = 0.0081Ra, (3.3) in the ultimate regime.Based on our data, we expect that to have deviations of less than 5 % from the ultimate regime, Ra in excess of 5 × 10 5 is needed.At those huge Rayleigh numbers, which closely represent the case of geological 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 ∼16 %.These findings are contrasted in figure 3 with the case of two-dimensional flow, for which the data generally support Nu ≈ 0.0069Ra + 2.75, as suggested by Hewitt et al. (2012), and attainment of a linear trend at Ra ∼ 3 × 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, δ, which, in Rayleigh-Darcy convection, is δ ∼ 30/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/Ra from the bottom boundary.Indication of the domain size in diffusive-convective units, x * d , for the simulations at lower Ra is given by the white boxes in panel (a).The complex flow pattern near the boundary, which appears to be still developing at Ra < 2 × 10 4 , seems to attain a roughly self-similar structure at larger 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 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. 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 80 000 40 000 20 000 10 000 supercells).The size of the supercells λ 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 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).
In figure 5, we show k r P(k r ) as a function of the normalized radial wavenumber k r /Ra, where k r = k 2 x + k 2 z , 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 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 /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 Ra, thus indicating growth of the supercells with 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 /Ra), which in this case are x + k 2 z , for all 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 /Ra and k i /Ra indicate the peaks of the spectra.The corresponding wavenumber is explicitly reported in the main panel.
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).

Conclusions
We have reported results of numerical simulation of Rayleigh-Darcy convection inside a cubic porous layer up to Rayleigh-Darcy number Ra = 8 × 10 4 .Exploiting the numerical dataset we show that the Nusselt number variation at sufficiently high Ra can be approximately expressed as Nu = 0.0081Ra + 0.067Ra 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 Ra 5 × 10 5 , yields the asymptotic scaling Nu = 0.0081Ra -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 Ra 5 × 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.t) is reported for a portion of the simulation, after the steady state is achieved (the beginning of the steady-state regime is indicated by t = t ss ).We observe that Nu(t) (solid line) oscillates within 1 % to 3 % of the mean value (dashed line).For all simulations, Nu(t) is shown in the range corresponding to ±10 % of the time-averaged value.Since each point of the Nu(t) curve is the result of a space average (over the entire top and bottom wall), the larger Ra is, the larger the number of points over which the space average is evaluated, and so the lower the amplitude of the fluctuations.

Figure 1 .
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, θ * ( y * = 0) = θ * max , and cooled at the top, θ * ( y * = h * ) = θ * min , and boundaries in the x * and z * directions are assumed to be periodic.The gravity acceleration (g) points downwards.The temperature distribution θ * for the case Ra = 1 × 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 * /Ra from the top boundary).

Figure 4 .
Figure 4. Temperature distribution in a plane close to the bottom boundary (x, y = 50/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 Ra is given in panel (a) by the white boxes.The size of the superstructures, estimated as λ i = 2πRa/k i with k i /Ra defined as in figure 5, is also shown in each panel (black bars).

911Figure 5 .
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/Ra, as in figure 4), shown as a function of the radial wavenumber k r = k 2x + k 2 z , for all 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 /Ra and k i /Ra indicate the peaks of the spectra.The corresponding wavenumber is explicitly reported in the main panel.

911Figure 6 .
Figure 6.Temperature distribution from the simulations S 1 (a,b) and S 8 (c,d) in a vertical slice located at x = 1/2.The dimensionless domain size in diffusive-convective units is indicated.A close-up view of the near-wall region, indicated with black squares in panels (a,c), is shown in panels (b,d).

Figure 7 .
Figure 7. Time evolution of instantaneous horizontally-averaged Nusselt number.Nu is shown for simulations S 1 (a), S 2 (b), S 3 (c), S 4 (d) and S 8 (e).The behaviour of Nu(t) is reported for a portion of the simulation, after the steady state is achieved (the beginning of the steady-state regime is indicated by t = t ss ).We observe that Nu(t) (solid line) oscillates within 1 % to 3 % of the mean value (dashed line).For all simulations, Nu(t) is shown in the range corresponding to ±10 % of the time-averaged value.Since each point of the Nu(t) curve is the result of a space average (over the entire top and bottom wall), the larger Ra is, the larger the number of points over which the space average is evaluated, and so the lower the amplitude of the fluctuations.

Table 1 .
Parameters of the main three-dimensional simulations performed in the present study.For each simulation, we report Rayleigh number Ra, grid resolution N x × N z × N y , Nusselt number Nu, and number of samples N s over which Nu is averaged (see appendix B for further details).Additional three-dimensional simulations (not listed here) have been performed at Ra = 2.5 × 10 3 , 5 × 10 3 .
Hewitt et al. (2014)Nu = 0.0081 Ra + 0.067 Ra 0.61 Nu = 0.0096 Ra + 4.6 Hewitt et al. (2014)Figure2.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, Nu = 0.0081Ra + 0.067Ra 0.61 , is shown as a black solid line.Results obtained in previous simulations byHewitt et al. (2014)(open triangles) and the extrapolated shifted linear scaling (blue dash-dotted line) are shown for comparison.The high-Ra portion of the (Ra, Nu) parameter space covered by previous investigations is rendered by the light grey rectangle.
Hewitt et al. (2012)amounts to linear dependence of Nu over Ra augmented with a sublinear corrective term, shown as a black solid line in figure 3. Compared with previous estimates, this 911 R4-5 https://doi.org/10.1017/jfm.2020.11782015)(, anrespectively) and three-dimensional domains (Hewitt et al. 2014) ( ) is also shown, withopen 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 Nu/Ra = 0.096 + 4.6/Ra predicted byHewitt et al. (2014).The scaling law Nu/Ra = 0.0069 + 2.75/Ra proposed byHewitt 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.911R4-6 https://doi.org/10.1017/jfm.2020.1178