# Asymptotic ultimate regime of homogeneous Rayleigh-Bénard convection on logarithmic lattices

Amaury Barral, Berengere Dubrulle

► To cite this version:

Amaury Barral, Berengere Dubrulle. Asymptotic ultimate regime of homogeneous Rayleigh-Bénard convection on logarithmic lattices. Journal of Fluid Mechanics, 2023, 962, pp.A2. 10.1017/jfm.2023.204. hal-04141097

HAL Id: hal-04141097

https://hal.science/hal-04141097

Submitted on 26 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
Banner appropriate to article type will appear here in typeset article

# Asymptotic Ultimate Regime of Homogeneous Rayleigh-Bénard Convection on Logarithmic Lattices

Amaury Barral¹ and Berengere Dubrulle¹†

¹SPEC/IRAMIS/DSM, CEA, CNRS, University Paris-Saclay, CEA Saclay, 91191 Gif-sur-Yvette, France

(Received xx; revised xx; accepted xx)

We investigate how the heat flux Nu scales with the imposed temperature gradient Ra in Homogeneous Rayleigh-Bénard convection using 1, 2 and 3D simulations on logarithmic lattices. Logarithmic lattices are a new spectral decimation framework which enables us to span an unprecedented range of parameters (Ra, Re, Pr) and test existing theories using little computational power. We first show that known diverging solutions can be suppressed with a large-scale friction. In the turbulent regime, for Pr ≈ 1, the heat flux becomes independent of viscous processes (“asymptotic ultimate regime”, Nu ~ Ra¹/² with no logarithmic correction). We recover scalings coherent with the theory developed by Grossmann & Lohse, for all situations where the large-scale frictions keep a constant magnitude with respect to viscous and diffusive dissipation. We also identify another turbulent friction dominated regime at Pr ≪ 1, where deviations from GL prediction are observed. These two friction dominated regimes may be relevant in some geophysical or astrophysical situations, where large scale friction arises due to rotation, stratification or magnetic field.

## 1. Introduction

Convection is a dynamical process that governs heat transport and mixing in a variety of systems ranging from planetary and astrophysical flows to industrial devices. In that respect, a crucial question is how the heat flux in the system is connected with the temperature gradient. Near equilibrium, where both quantities are small, Fourier laws apply, and the heat flux is simply proportional to the temperature gradient. For larger values, the system enters a non-linear then turbulent regime, where thermal energy in converted into mechanical energy, and the relation becomes nonlinear. The deviations from linearity are quantified by the relation between the Nusselt number, Nu, the ratio between the heat flux and its laminar value, and the Rayleigh number Ra, the non-dimensional temperature gradient.

In fluid mechanics, the paradigmatic system describing convection is a fluid enclosed in a volume, in which thermal energy is injected at the bottom via imposed heat flux or temperature gradient. Its dynamics is described by the Rayleigh-Bénard (RB) equations. Despite decades of theoretical, experimental and numerical developments, the scaling of the heat transfer in RB remains a subject of discussion and active research. In bounded domains at low Ra, a simple argument by Malkus & Chandrasekhar (1954) based on the criticality

<footnote>† Email address for correspondence: berengere.dubrulle@cea.fr</footnote>
<table><caption>**Table 1:** Scaling predictions for HRB observables in the turbulent regime with and without friction. The observables are given by Table (2). DNS stands for Direct Numerical Simulation using regular Fourier modes (Calzavarini et al. 2005) while LL refers to simulations using Fourier modes on a LogLattices (this paper). $U_{ls}^2$ and $\Theta_{ls}^2$ are large scale kinetic and thermal energy. Exponents are computed by fitting over $Ra > 10^7$ (resp. $1 < Pr < 50$) for varying $Ra$ (resp. $Pr$). Errors represent std of fit parameters.</caption><thead><tr><th></th><th>GL Theory f = 0</th><th>DNS f = 0</th><th>LL f = 1</th></tr></thead><tbody><tr><td>Nu ~ &radic;<span style="text-decoration: overline;">Ra</span>Pr</td><td>Pr<sup>0.43</sup> Ra<sup>0.50</sup></td><td>Pr<sup>0.51&plusmn;0.01</sup> Ra<sup>0.53&plusmn;0.03</sup></td><td></td></tr><tr><td>Re ~ &radic;<span style="text-decoration: overline;">Ra</span>/Pr</td><td>Pr<sup>-0.55</sup> Ra<sup>0.5</sup></td><td>Pr<sup>-0.54&plusmn;0.01</sup> Ra<sup>0.54&plusmn;0.01</sup></td><td></td></tr><tr><td>&epsilon;<sub>&theta;</sub> ~ c<sub>1</sub>&radic;<span style="text-decoration: overline;">Re/Ra</span> + c<sub>2</sub> Re &radic;<span style="text-decoration: overline;">Pr/Ra</span></td><td>(Re Pr)<sup>-0.17</sup></td><td>Re<sup>x</sup> Pr<sup>x-0.5</sup>/&radic;<span style="text-decoration: overline;">Ra</span>, 1 &le; x &le; 1.2</td><td></td></tr><tr><td>&epsilon;<sub>u</sub> ~ Re<sup>3</sup>(Pr/Ra)<sup>3/2</sup></td><td>Re<sup>2.77</sup>(Pr/Ra)<sup>3/2</sup></td><td>Re<sup>2.88&plusmn;0.03</sup>(Pr<sup>0.95&plusmn;0.01</sup>/Ra)<sup>3/2</sup></td><td></td></tr></tbody></table>

of the thermal boundary layer gives $Nu \sim Ra^{1/3}$, observed in many experiments (see Ahlers et al. (2009) for review). As we increase $Ra \to \infty$, viscous processes (and their associated boundary layers) are believed to become irrelevant, resulting in an “ultimate regime of convection”, where $Nu \sim Ra^{1/2}$ (hereafter called “asymptotic ultimate regime”) (Spiegel 1963; Grossmann & Lohse 2000), with possible logarithmic corrections (Kraichnan 1962; Grossmann & Lohse 2011) (hereafter called “ultimate regime”). Experimental or numerical observations of the (asymptotic) ultimate regime prove to be very difficult, and no final consensus has been reached so far about its existence in a pure RB setting (Chavanne et al. 1997; Urban et al. 2019; Doering & Constantin 1996; Zhu et al. 2018, 2019a; Roche 2020) (see Ahlers et al. (2009) for a less recent but more synthetic review). When the gravity is artificially increased using centrifugal force, one can indeed observe hints of an ultimate regime (Jiang et al. 2022). On the other hand, various modifications of the RB geometry aiming at modifying the influence of the boundary layers result in experimental observation of a regime where $Nu \sim Ra^{1/2}$: using highly elongated cells (Castaing et al. 2017; Pawar & Arakeri 2016), using rough (Ciliberto & Laroche 1999; Rusaouën et al. 2018; Zhu et al. 2019b; Kawano et al. 2021) or porous (Zou & Yang 2021; Motoki et al. 2022) boundaries, or radiatively heating the flow (Lepot et al. 2018; Bouillaut et al. 2019).

From a numerical point of view, a simple way to remove boundary layers is to consider a triply periodic geometry, and heat the flow via an applied temperature gradient. This setting was first explored by Borue & Orszag (1997); Lohse & Toschi (2003); Calzavarini et al. (2005, 2007) and called Homogeneous Rayleigh-Bénard (HRB) convection. The corresponding scalings and predictions are summarized in Table 1. Although the results of those simulations are consistent with the predictions of Grossmann & Lohse (2000) (hereafter called "GL theory"), they are undermined by several drawbacks: statistics polluted by the growth of uncontrolled exponential instabilities (Calzavarini et al. 2006) of unclear physical relevance, a small $Ra$ and $Pr$ range, sparse data points due to difficulties in running numerically challenging simulations. Indeed, pushing the Rayleigh number to large values increases the numerical burden beyond the capacity of present computers, as the number of grid points needed to describe the flow usually scales like $Re^3$ with $Re \sim Ra^{1/2}$. In an attempt to reduce the number of degrees of freedom, models based on sparse interacting Fourier modes have been recently devised (Campolina & Mailybaev 2018, 2021). Those modes are evenly spaced points in log space (thus thereafter called “log-lattice”) and are interacting via nonlinear equations that are derived from the fluid equations by substituting for the
convolution product a new operator, that can be seen as a convolution on the log-lattice, and preserves all the main symmetries and conservation laws of the original equations. As such, log-lattices are likely to preserve properties of the original equations that are directly linked to these symmetries and conservation laws. This was indeed checked for the Burgers and Navier-Stokes equation in the Fourier space (energy spectrum, energy transfers), over an unprecedented wide range of scales (Campolina & Mailybaev 2021). Another interesting feature of log-lattices is that in 1D, they encompass classical shell models of turbulence for special values of the log-lattice spacing (Campolina & Mailybaev 2021), such as the Sabra shell model of turbulence (Gloaguen et al. 1985; Biferale 2003).

1D shell models of turbulence were used previously in the context of HRB (Ching & Ko 2008) in an effort to increase the Ra and Pr range of results. They successfully display the asymptotic ultimate regime of convection, at the price of tuning several parameters of the model to get rid of the uncontrolled exponential instabilities. This, combined with the 1D nature of the model, renders the informative and conclusive nature of the observations questionable. The goal of the present letter is therefore to re-explore the HRB equation using the log-lattice framework, that allows both the exploration of a wide range of parameters on a large array of wavenumbers, and a flexibility of dimensionality from 1D to 3D, at low numerical cost, and without additional empirical parameters. Given that they preserve all main conservation laws and symmetry of the original HRB equation, many features of the original equation are still valid, like the exact conservation laws of Table 1. Whether the GL theory still applies, and what are the modifications of the asymptotic ultimate regime implied by the log-lattice geometry are interesting open questions that we investigate here. In that respect, the present paper offers an exploration of the analogy and differences between log-lattices and classical fluid dynamics in a more complex case (HRB) than previous examples (Campolina & Mailybaev 2018, 2021).

## 2. Numerical simulations

### 2.1. Generalities

The dynamics of a homogeneous fluid, with coefficient of thermal dilation $\alpha$, viscosity $\nu$ and diffusivity $\kappa$, subject to a temperature gradient $\Delta T$ over a length $H$ and vertical gravity $g$ is given by the HRB set of equations (Lohse & Toschi 2003; Calzavarini et al. 2005, 2006, 2007),

$$ 
\begin{aligned}
\partial_t u + u \cdot \nabla u + \frac{1}{\rho_0} \nabla p = \nu \nabla^2 u + \alpha g \theta z, \\
\partial_t \theta + u \cdot \nabla \theta = \kappa \nabla^2 \theta + u_z \frac{\Delta T}{H}, \\
\nabla \cdot u = 0,
\end{aligned} 
\quad (2.1) $$

where $u$ is the velocity, $\theta$ the temperature fluctuation, $\rho_0$ is the (constant) reference density and $p$ is the pressure. Here, the mean temperature gradient $\Delta T$ acts as a forcing term. This gradient is non-dimensionalized into the Rayleigh number $Ra = \alpha g H^3 \Delta T / (\nu \kappa)$. The Prandtl number $Pr = \nu / \kappa$ is the ratio of the fluid viscosity to its thermal diffusivity. The mean total heat flux is the $z$ direction is $J = \langle u_z \theta \rangle - \kappa \Delta T$ which is adimensionalized into $Nu = J H / \kappa \Delta T$.

Taking global space and time average of the equation (2.1), one can derive (Lohse & Toschi 2003; Calzavarini et al. 2005) two exact relation for the volume averaged kinetic and
<figure><img src="image_2.png" /><figcaption>**Figure 1:** Absolute value of the rate of growth of instability $\sigma = d \log X/d \log t$ where $X = \langle u\theta \rangle$ without large-scale friction ($f = 0$), vs. Rayleigh number. The green dashed line is the theoretical growth rate for $k = k_c = 2\pi\sqrt{3}$, corresponding to Eq. (2.4). The interval $k < k_c$ corresponds to negative values of $\sigma$.</figcaption></figure>

thermal dissipation, which respectively scale as

$$v \langle (\partial_t u_j)^2 \rangle_V = v^3 H^{-4} \text{ Nu Ra Pr}^{-2}, \quad (2.2)$$

$$\kappa \langle (\partial_t \theta)^2 \rangle_V = \kappa H^{-2} (\Delta T)^2 \text{ Nu.} \quad (2.3)$$

Additionally, to get rid of the pressure term, we take the rotational of the above equation ($\omega = \text{rot}u = ik \times u$).

## 2.2. Quantities of interest

### 2.3. Adaptation on log-lattices: HRB with friction

#### 2.3.1. Exponential instabilities in HRB

As first shown by Calzavarini et al. (2006), HRB equations are prone to exponential instabilities, due to the conservation of the total energy. In the absence of large-scale friction, we also observe those instabilities in our log-lattice simulations (Figure 2a). As shown in Figure (1), the growth rate of the instability in the log-lattice simulations matches the theoretical growth rate given by Calzavarini et al. (2006); Schmidt et al. (2012):

$$\sigma \sqrt{\text{Ra Pr}} = \frac{1}{2} \left[ \sqrt{((\text{Pr}+1)k^2)^2 + 4 \text{Pr}(\text{Ra}-k^4)} - (\text{Pr}+1)k^2 \right] \sim \sqrt{\text{Ra}}, \quad (2.4)$$

for $\theta, u \sim e^{\sigma t + i\vec{k} \cdot \vec{x}}$. This expression yields unstable solutions for $\text{Ra} > \text{Ra}_c = k_{\min}^4$ where $k_{\min}$ is the modulus of the smallest mode on the grid, which is $2\pi\sqrt{3}$ in our case.

However, the non-linear behavior of the instability in the log-lattice case is quite different from the one reported by Calzavarini: instabilities tend to extend significantly further and for longer times. Our interpretation is that in our log-lattice model, the modes are not coupled enough to develop the nonlinear saturation. The instabilities widely interfere with the statistical stability of observables and need to be removed for a meaningful analysis. Physically, these exponential ramps originate for a lack of energy sink to absorb the constant energy injection in the bulk by the (fixed) temperature gradient. Previous works on 1D simulations (Ching & Ko 2008) have shown that without a large-scale sink to counteract this source, energy diverges at large scales and scaling laws become incorrect. Therefore, to get rid of the exponential instabilities, we include a large-scale friction $f$ on both $u$ and $\theta$. By doing so, the instability saturates, and we achieve a statistically stationary state for
<figure><img src="image_3.png" /><figcaption>**Figure 2:** Influence of the large-scale friction on the time behavior of the Nusselt number Nu in 3D HRB. (2a) Without friction: we observe the growth of an exponential instability. (2b) With friction: the instability saturates and the dynamics become statistically convergent. Parameters: Ra = $10^6$, Pr = 1, N = 13.</figcaption></figure>

<caption>**Table 2:** Physical quantities expressed as a function of the non-dimensional variables of Eq. (2.5). $\langle \cdot \rangle$ denotes the temporal and spatial average.</caption>

$$ 
\begin{aligned}
\text{Nu} &= \frac{JH}{\kappa\Delta T} - 1 \rightarrow \sqrt{\text{RaPr}} \cdot \langle u_z \theta \rangle - 1 \\
\text{Re} &= \frac{\sqrt{\langle U_i U_i \rangle H}}{\nu} \rightarrow \sqrt{\frac{\text{Ra}}{\text{Pr}}} \cdot \sqrt{\langle u_i u_i \rangle} \\
\epsilon_{\theta} &= \kappa \langle (\partial_i \Theta)^2 \rangle \rightarrow \frac{\langle (\partial_i \theta)^2 \rangle}{\sqrt{\text{RaPr}}} \\
\epsilon_u &= \nu \langle (\partial_i U_j)^2 \rangle \rightarrow \sqrt{\frac{\text{Pr}}{\text{Ra}}} \cdot \langle (\partial_i u_j)^2 \rangle
\end{aligned}
$$

the heat transfer, as displayed in figure 2. Note however that the fluctuations of Nu around the stationary value are very broad, and extend over one or two orders of magnitudes. The same phenomenon was observed in the DNS of HRB (Calzavarini et al. 2005, 2006) and mentioned to be a source of difficulty to achieve reliable results (Borue & Orszag 1997). For this reason, very long simulations are necessary to get steady averages (Pumir & Shraiman 1995; Calzavarini et al. 2006). In DNS, this cannot be achieved without cutting down the resolution, which may impact the reliability of dissipation estimates (Yeung et al. 2018). In the log-lattice framework, we do not have this problem, and we performed high resolution very long time averages on the log of Nu, and represent all quantities in log-log variables.

### 2.3.2. Equations

To investigate the ultimate regime, it is natural to adimensionalize the equation in terms of “inertial quantities”, i.e. using the vertical width $H$ as a unit of length, the free fall velocity $U_{ff} = \alpha g \Delta T H$ as a unit of velocity, and $\Delta T$ as a unit of temperature. Table 2 indicates the form taken by observables after rescaling as indicated. The equations including the temperature gradient and the friction can then be written in terms of velocity as (with the Einstein
convention on summed repeated indices):

155

$$ \partial_t u_i = \mathbb{P} \left[ -u_j \partial_j u_i + \theta \delta_{i=z} + \sqrt{\frac{\mathrm{Pr}}{\mathrm{Ra}}} \nabla^2 u_i - f u_i \delta_{k \approx k_{\mathrm{min}}} \right]_i, \quad (2.5) $$

$$ \partial_t \theta = -u_i \partial_i \theta + u_z + \frac{\nabla^2 \theta}{\sqrt{\mathrm{Ra}} \mathrm{Pr}} - f \theta \delta_{k \approx k_{\mathrm{min}}}, $$

where the Dirac $\delta_{k \approx k_{\mathrm{min}}}$ filters out the small scales, and the projector, given in the Fourier

space by $\mathbb{P}(\mathbf{A}) = \mathbf{A} - \frac{k_i}{k^2} k_j A_j$, accounts for the pressure term under the divergence-free condition. We also looked at those equations expressed in terms of the vorticity $\omega = \nabla \times \mathbf{u}$:

160

$$ \partial_t \omega_i = -\omega_j \partial_j u_i - u_j \partial_j \omega_i + \theta[\nabla \times \mathbf{z}]_i + \sqrt{\frac{\mathrm{Pr}}{\mathrm{Ra}}} \nabla^2 \omega_i - f \omega_i \delta_{k \approx k_{\mathrm{min}}}, \quad (2.6) $$

$$ \partial_t \theta = -u_i \partial_i \theta + u_z + \frac{\nabla^2 \theta}{\sqrt{\mathrm{Ra}} \mathrm{Pr}} - f \theta \delta_{k \approx k_{\mathrm{min}}}, $$

Adding a large-scale friction to damp the inverse cascade is a classical trick-it is e.g.

162 routinely used numerical simulations of 2D turbulence to avoid Bose condensation at $k = 0$

163 and enable stationarity (Sukoriansky et al. 1999). The present case is 3D, but we interpret

164 the formation of exponential ramps as a signature of back-scattering of energy, a feature that

165 was already mentioned previously in shell models of Rayleigh-Bénard convection (Ching

166 & Ko 2008). The addition of the friction is therefore a convenient way to damp the large-

167 scale modes that are generated by the large-scale instability. Such friction is also added in

168 many models of climate, as a subgrid model to account for the friction at the boundary

169 layer that cannot be resolved in the stratified case. The hand waving argument is that, within

170 boundary layers, a shear profile develops, with extraction of energy at the boundaries, which

171 is proportional to the square of the shear. Assuming the shear to be constant in the boundary

172 layer, we can then estimate it by the difference between the velocity at the top of the layer,

173 minus the velocity at the boundary which is zero. In total, the energy pumped by friction is

174 proportional to the square of the velocity, which is exactly the law we have implemented.

175 Such friction is termed Rayleigh friction in the climate community (Stevens et al. 2002) and

176 can actually be seen as a way to take into account the boundary conditions that we have

177 removed in the HRB setting.

### 2.3.3. *Conservations laws for HRB with and without friction*

In the absence of friction, the conservation laws for HRB are given by Eqs. (2.2) and by (2.3).

179 The presence of the friction just adds a supplementary term proportional to $f$ in each equation.

180 The result can be made non-dimensional using $U_{\mathrm{ff}}$, $H$ and $\Delta T$ as units of velocity, length and

181 temperature, resulting in:

183

$$ f \langle u^2 \delta_{k \approx k_{\mathrm{min}}} \rangle + \epsilon_u = \frac{\mathrm{Nu} + 1}{\sqrt{\mathrm{Ra}} \mathrm{Pr}}, \quad (2.7) $$

184

$$ f \langle \theta^2 \delta_{k \approx k_{\mathrm{min}}} \rangle + \epsilon_{\theta} = \frac{\mathrm{Nu} + 1}{\sqrt{\mathrm{Ra}} \mathrm{Pr}}, \quad (2.8) $$

185

From now on, we define $U_{\mathrm{ls}}^2 = \langle u^2 \delta_{k \approx k_{\mathrm{min}}} \rangle$ and $\Theta_{\mathrm{ls}}^2 = \langle \theta^2 \delta_{k \approx k_{\mathrm{min}}} \rangle$.

186
## 2.4. Log-lattices

Log-lattice models fit into the more general framework of REduced Wavenumber set Approximation (REWA) (Grossmann et al. 1994)) or fractal decimated models (Frisch et al. 2012; Lanotte et al. 2015). The spirit of these methods is to use a reduced subset of modes obeying a well-defined hierarchy, so as to stick closer to the observed organized nature of turbulence. In the original REWA models (Grossmann et al. 1994), non-linear-interactions are projectively decreased either in a random manner or such that they are distributed over a fractal set (Frisch et al. 2012; Lanotte et al. 2015). In log-lattice models, the modes reduction is achieved by keeping modes following a geometric progression, thereby allowing to reach very small scales with a very small number of modes. The construction is detailed in Campolina & Mailybaev (2021), where it is shown that fluid equations on log-lattices respect all symmetries of the Euler equations, and retain classical and basic properties of the Navier-Stokes equation, such as constancy of energy flux in the inertial range.

There are several key differences compared to shell models (Brandenburg 1992; Ching & Ko 2008) or the original REWA model. Like in a shell model, simulations are carried out in Fourier space on an logarithmically-decimated grid. Unlike shell models, log-lattices are truly multidimensional, and unlike the original REWA model, the decimation does not have a fixed number of points per shell: $k(n_1, \dots, n_d) = \sum_i \lambda^{n_i} e_i$, $n_i \in \mathbb{Z}$ with $d$ the spatial dimension and $e_i = x, y, z, \dots$. Log-lattices are endowed with a scalar product:

$$ (f, g) = \Re \left( \sum_k f(\mathbf{k}) \overline{g(\mathbf{k})} \right), \qquad (2.9) $$

and a convolution operator:

$$ (f * g)(\mathbf{k}) = \sum_{\substack{p,q \\ p+q=k}} f(p)g(q), \qquad (2.10) $$

that naturally extend the corresponding operators on regular Fourier grids. This ensures that the log-lattice operators respect the symmetries of the Navier-Stokes equation, which ensures the conservation of energy, helicity (3D) or enstrophy (2D), etc. provided that they are conserved in the original equation. The constrain on the interacting triads on log-lattices $\exists p, q \in \lambda^Z : p+q \in \lambda^Z$ restricts the acceptable values of $\lambda$ to three main families: $\lambda = 2$, the plastic number $\lambda = \rho \approx 1.324$, and $\lambda^b - \lambda^a = 1$, $(a, b) \in \mathbb{N}^2$, whose biggest solution is the golden number $\lambda = \phi \approx 1.618$. From a numerical point of view, $\lambda = 2$ is the “fastest” option, as it has both a maximal span for a given number of points, and the least interactions per point. However, as outlined in the next part, we believe that $\lambda = 2$ should be avoided for incompressible simulations. We thereafter perform all our simulations with $\lambda = \phi$, which is the second biggest value of $\lambda$, and has the second least number of interactions per grid point.

## 2.5. Numerical details

### 2.5.1. Configuration

The minimum wave vector of the grid is set to $k_{min} = 2\pi$ to match a simulation on a box of size $L = 1$. The grid size $N$ is then set so as to reach the dissipative scale both for velocity and temperature. We alternate between several initial condition (IC) choices for our simulations: large-scale initialization, Kolmogorov spectrum, flat-spectrum. All those choices are modulated by a weak multiplicative complex noise. We find no significant influence of those initial conditions on the scaling laws. As Ra or Pr increase, the simulations become slower and slower. This sets the upper bound on the range of parameters we can
<figure><img src="image_4.png" /><figcaption>**Figure 3:** Exact conservation laws for $\epsilon_{\theta}$ in 3D results. Black points correspond to varying Ra, gray points correspond to varying Pr. (3a) $\epsilon_{\theta} + f\Theta_{ls}^2$ vs $(Nu+1)/\sqrt{RaPr}$. (3b) Compensated plot $(\epsilon_{\theta} + f\Theta_{ls}^2)/((Nu+1)/\sqrt{RaPr})$ vs $(Nu+1)/\sqrt{RaPr}$.</figcaption></figure>

229 integrate while retaining statistically relevant observables in a reasonable simulation time (one CPU days at most). In 3D, this yields $Ra_{max} \approx 10^{10}$ for Pr = 1 and $Pr_{max} \approx 5 \cdot 10^4$ for Ra = $10^8$. The lower bound is set by the value of the Nusselt number, which must obey

230
231 Nu ≫ 1, the value Nu ≈ 1 corresponding to the laminar regime with trivial scaling laws.

232
233 Finally, integrating equations on log-lattices yields interesting and new numerical challenges.

234 We built our own ODE integrator to solve them, as detailed in Supplementary Materials.

235 Once we have ran a simulation for a long enough time, we compute Nu, $\epsilon_{\theta}$, $\epsilon_{u}$ by taking long time and space averages (with $\langle ab \rangle = \frac{1}{T} \int_{t}^{T} dt(a, b)$) according to table 2. The accuracy

237 of our results is controlled by checking that we recover the exact laws of HRB convection

238 Eqs (2.8) and (2.7). This is shown in Fig. 3 and 4, for all 3D data sets used in the present

239 paper (see Table 3). Furthermore, the ratio between the friction term and the dissipation is

240 shown in Fig. 5.

241 2.5.2. *Simulation sets*

242 The results we obtained come from seven types of simulation that are described in the Table 3.

243 For comparison, we also included in some graphs the results by (Calzavarini et al. 2005),
244 obtained using DNS of the same equations, but at $f=0$.

245 Historically, we performed first vorticity simulations, then velocity simulations, improving
246 the integrator scheme in between to be able to better handle various numerical challenges
247 raised by simulating wavenumbers as high as $k \sim 10^5$ in 3D. For transparency reasons,
248 we decided to include all datasets we had at our disposal, but we believe that the velocity
249 simulations are the more faithful ones, in the sense that they deal better with the small
250 scales at large Rayleigh or Reynolds number. This sensitivity to small scale modeling (and
251 resolution) is also a well-known feature of direct numerical simulations, especially when it
252 comes to statistics of gradients or energy dissipation (Yeung et al. 2018).

253 We have verified that the size of the grid for 3D simulations ($N = 13$) does not affect
254 the mean value of the observables Nu, Re, ..., which is already converged for grids of size
255 $N \ge 6$. However the tail of the pdfs does depend on $N$. Another 3D simulation set at $N = 20$
256 (not shown here, both vs Ra and Pr) displays the same scaling laws as the $N = 13$ case,
257 confirming this analysis.
<figure><img src="image_5.png" /><figcaption>**Figure 4:** Exact conservation laws for $\epsilon_u$ in 3D results. Black points correspond to varying Ra, gray points correspond to varying Pr. (4a) $\epsilon_u + fU_{ls}^2$ vs $(Nu+1)/\sqrt{RaPr}$. (4b) Compensated plot $(\epsilon_u + fU_{ls}^2)/((Nu+1)/\sqrt{RaPr})$ vs $(Nu+1)/\sqrt{RaPr}$.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5:** Ratio between friction $f_x = fX_{ls}^2$ and dissipation $\epsilon_x$ for $x = u, \theta$ (5a) versus Ra at Pr = 1 and (5b) versus Pr at Ra = $10^8$.</figcaption></figure>

### 2.5.3. A case against $\lambda = 2$

This section explains why the log-lattice parameter $\lambda = 2$ is ill-suited to simulating divergence-free equations. It is not specific to HRB simulations, however we believe this issue has not been reported in a publication before.

$\lambda = 2$ is the biggest grid parameter that can be accommodated on a log-lattice. For a fixed grid size $N$ in dimension $D$, it is therefore very tempting to use $\lambda = 2$, since among all the $\lambda$s it spans the greatest range of wavenumbers (the convolution's complexity rises as $O(N^D)$). However, $\lambda = 2$ misrepresents the convection term $u_j \partial_j u_i$.

The heart of the problem is easily understood through a simple 2D example. Consider the convection term $u_x \partial_x \omega + u_y \partial_y \omega$ of a divergence-free flow, with a large-scale initialization $u(k > k_0) = \omega(k > k_0) = 0$ for some $k_0$. From a physical point of view, we expect convection to populate the $k \ge k_0$ region as time advances. However, with $\lambda = 2$, this does not happen, as is demonstrated below.
<table><caption>**Table 3:** Parameters of the data sets used in the present paper. *D* is the dimension. The "velocity" datasets are obtained by integration of Eq. (2.5), while the "vorticity" datasets are obtained by integration of Eq. (2.6). DNS refers to direct simulations of (Calzavarini et al. 2005), using a classical spectral Fourier code (on a regular grid). The ++ label refers to an integration using an improved integrator, using a reshuffling of variable matrices that allows faster simulations. The Ra and Pr column provides the Rayleigh and Prandtl number range of the simulations. *f* is the large scale friction, *N* = 1 + log *k*<sub>max</sub>/log(*φ*), where *k*<sub>max</sub> is the maximal wavenumber of the simulation and *φ*, the golden mean, is a measure of the spatial resolution. For log-lattices, its corresponds to the number of modes in each direction. *N*<sub>av</sub> is the length of the simulation, divided by the large eddy turnover time. It provides the number of decorrelated frames that can be used to estimate statistical averages. The tolerance refers to the absolute and relative tolerances that are fixed equal in all the simulations.</caption><thead><tr><th>Name</th><th>D</th><th>Type</th><th>Ra</th><th>Pr</th><th>f</th><th>N</th><th>N<sub>av</sub></th><th>Tolerance</th><th>Symbol</th></tr></thead><tbody><tr><td>(I)</td><td>1D</td><td>Velocity</td><td>[10<sup>5</sup>, 10<sup>50</sup>]</td><td>1</td><td>1</td><td>120</td><td>-</td><td>10<sup>-3</sup></td><td>●</td></tr><tr><td>(II)</td><td>2D</td><td>Vorticity</td><td>[10<sup>5</sup>, 10<sup>50</sup>]</td><td>1</td><td>1</td><td>20</td><td>-</td><td>10<sup>-3</sup></td><td>◆</td></tr><tr><td>(III)</td><td>3D</td><td>Velocity</td><td>[1, 10<sup>10</sup>]</td><td>1</td><td>1</td><td>13</td><td>&gt; 480</td><td>10<sup>-6</sup></td><td>●</td></tr><tr><td>(IV)</td><td>3D</td><td>Vorticity</td><td>[1, 10<sup>10</sup>]</td><td>1</td><td>1</td><td>13</td><td>&gt; 480</td><td>10<sup>-6</sup></td><td>●</td></tr><tr><td>(V)</td><td>3D</td><td>Velocity</td><td>10<sup>8</sup></td><td>[5 &middot; 10<sup>-4</sup>, 10<sup>2</sup>]</td><td>1</td><td>13</td><td>&gt; 50</td><td>10<sup>-6</sup></td><td>●</td></tr><tr><td>(VI)</td><td>3D</td><td>Vorticity</td><td>10<sup>8</sup></td><td>[5 &middot; 10<sup>-4</sup>, 10<sup>2</sup>]</td><td>1</td><td>13</td><td>&gt; 50</td><td>10<sup>-6</sup></td><td>●</td></tr><tr><td>(VII)</td><td>3D</td><td>++Velocity</td><td>{10<sup>9</sup>, 10<sup>10</sup>, 10<sup>11</sup>}</td><td>[5 &middot; 10<sup>-4</sup>, 10<sup>2</sup>]</td><td>1</td><td>13</td><td>&gt; 80</td><td>10<sup>-6</sup></td><td>-</td></tr><tr><td>(VIII)</td><td>3D</td><td>Velocity</td><td>[10<sup>6</sup>, 10<sup>10</sup>]</td><td>1</td><td>1</td><td>13</td><td>&gt; 50</td><td>10<sup>-6</sup></td><td>●</td></tr><tr><td>Calzavarini</td><td>3D</td><td>DNS</td><td>[10<sup>5</sup>, 10<sup>8</sup>]</td><td>[10<sup>-1</sup>, 10]</td><td>0</td><td>-</td><td>&gt; 64</td><td>-</td><td>◆</td></tr></tbody></table>

In a divergence-free flow, $u_x * \partial_x \omega = -i \left( \frac{\omega k_y}{k^2} * k_x \omega \right)$, $u_y * \partial_y \omega = i \left( \frac{\omega k_x}{k^2} * k_y \omega \right)$ where * denotes a convolution. In a $\lambda = 2$ log-lattice, convolutions are defined as† $f * g(\lambda^n, \lambda^m) = f(\lambda^{n-1}, \lambda^{m-1}) \cdot g(\lambda^{n-1}, \lambda^{m-1})$. Due to the initial conditions, this yields:

$$ (u_x * \partial_x \omega + u_y * \partial_y \omega) (k \approx k_0) = 0 $$

There is no forward convection at all, therefore there can be no forward cascade in such case.

This does not happen for other values of $\lambda$, for which the convolution is evaluated at asymmetric positions. We therefore advise against using $\lambda = 2$ in divergence-free fluids, and suggest to rather use $\lambda = \phi$ (the second-biggest grid parameter).

### 2.5.4. Zero-divergence problem in 1D

In the 1D case, we cannot impose the zero-divergence condition, so that quantities like $u_x \partial_x \theta$ and $\partial_x(u_x \theta)$ are not equivalent. Here, we have followed the same choice than Ching & Ko (2008), and wrote the equation as:

$$ \begin{aligned} \partial_t u &= -u \partial_x u + \theta + \sqrt{\frac{\text{Pr}}{\text{Ra}}} \nabla^2 u - f u \delta_{k \approx k_{\text{min}}}, \\ \partial_t \theta &= -u \partial_x \theta + u + \frac{\nabla^2 \theta}{\sqrt{\text{Ra Pr}}} - f \theta \delta_{k \approx k_{\text{min}}}. \end{aligned} \quad (2.11) $$

† excluding the $k = 0$ mode, which is not used in this paper
<figure><img src="image_7.png" /><figcaption>**Figure 6:** Non-dimensional heat transfer Nu vs Rayleigh number Ra in 1D and 2D. Correspondence between symbols and datasets are given in Table 3. (6a) Nu vs Ra. The gray dashed line corresponds to $Nu \sim \sqrt{Ra}$, corresponding to ultimate regime scaling. (6b) Compensated plot $A \text{Nu}/\sqrt{Ra}$ vs Ra., where A is adjusted to collapse the 1D and 2D data in the ultimate regime.</figcaption></figure>

# 3. Results and Discussion

## 3.1. One and two-dimensional cases

Figure 6 presents the Nu vs Ra scaling in 1D and 2D. The 1D Nu scaling law extends over 50 orders of magnitude in Ra (Fig. 6a), and follows closely the law $Nu \sim \text{Ra}^{1/2}$, as can be checked by the compensated plot in Fig. (6b), in agreement with Ching & Ko (2008). In 2D, the scaling also extends approximately over 30 orders of magnitudes for $Ra > 10^{23}$. Moreover, the compensated plot highlights small fluctuations around this law, see Fig. (6b), due to statistical noise.

## 3.2. In 3D

In 3D, the simulations get significantly more turbulent and results are subject to more statistical fluctuations. Another source of fluctuations comes from a physical phenomenon, associated with the existence of friction. To showcase this effect, we plot in Fig. 5a and Fig. 5b the ratio between the energy dissipated by friction and the energy dissipated by viscosity or diffusivity for both the kinetic energy and the thermal energy.

Fixing $Pr = 1$ and varying Ra between $10^3$ and $10^8$, we observe in Fig. 5a that both $f_u = f U_{ls}^2 / \epsilon_u$ or $f_\theta = f \Theta_{ls}^2 / \epsilon_\theta$ behave in the same way as a function of Ra at low Ra, the dissipation due to friction is small, and gradually increases towards reaching a plateau around $Ra \sim 10^7$, where energy dissipated by frictions reach about 90% of the energy dissipated by viscosity or diffusivity. We can thus define a “non-universal” regime where $f/\epsilon$ depends on Ra, Pr and a “universal” regime where $f/\epsilon$ does not depend on Ra, Pr.

The critical Rayleigh number where the plateau occurs is likely to depend on the Prandtl number. To check this; we now fix $Ra = 10^8$ and vary $Pr$ from several order of magnitude. In Fig. 5b, we then observe an interesting symmetrical behaviour, with respect to $Pr = 1$: decreasing $Pr$, we observe that the energy dissipated by the velocity friction remains of the same order of magnitude than the dissipation by viscosity, while the energy dissipated by thermal friction strongly decays and become negligible. As Pr shifts away from 1, we observe the symmetrical behavior, with velocity friction becoming negligible, while thermal friction remains of the same order of magnitude than the thermal energy dissipation. As we will see,
<figure><img src="image_8.png" /><figcaption>**Figure 7:** Non-dimensional heat transfer Nu vs Rayleigh number Ra in 3D for Pr = 1. Correspondence between symbols and datasets are given in Table 3. The gray dashed line separates the non-universal (left) and the universal (right) friction dominated regimes for data corresponding to Fig. 5. (7a) Nu vs Ra. The black dashed line corresponds to Nu ~ $\sqrt{Ra}$, corresponding to asymptotic ultimate regime scaling. (7b) Compensated plot $Nu/\sqrt{Ra}$ vs Ra.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8:** Scaling of non-dimensional heat transfer Nu as a function of Prandtl number Pr in 3D for Ra = $10^8$. Correspondence between symbols and datasets are given in Table 3. The gray dashed line separates the non-universal (left) and the universal (right) friction dominated regimes for data corresponding to Fig. 5. (8a) Nu vs Pr. The black dashed line corresponds to Nu ~ $\sqrt{Pr}$, corresponding to asymptotic ultimate regime scaling. (8b) Compensated plot $Nu/\sqrt{Pr}$ vs Pr.</figcaption></figure>

this will have an impact on the thermal transport. Note that at small (resp. large) Pr, all the thermal (resp. velocity) modes become concentrated at large scale, where the friction occurs.

Therefore, in the large Pr regime, the kinetic friction and viscous dissipation compete, while at small Pr the same remark holds for the thermal friction and diffusive dissipation. This may then explain the vanishing of the friction in those regime.

We now focus on the regimes where the ratio of friction to dissipation is approximately constant. These regimes are friction dominated, but, as we will see, are characterized by interesting universal scaling regimes.
<figure><img src="image_10.png" /><figcaption>**Figure 9:** Scaling of Reynolds number Re as a function of Rayleigh number Ra in 3D for Pr = 1. Correspondence between symbols and datasets are given in Table 3. The gray dashed line separates the non-universal (left) and the universal (right) friction dominated regimes for data corresponding to Fig. 5. (9a) Re vs Ra. The black dashed line corresponds to Re ~ $\sqrt{Ra}$, corresponding to asymptotic ultimate regime scaling. (9b) Compensated plot Re/√Ra vs Ra.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 10:** Scaling of Reynolds number Re as a function of Prandtl number Pr in 3D for Ra = $10^8$. Correspondence between symbols and datasets are given in Table 3. The gray dashed line separates the non-universal (left) and the universal (right) friction dominated regimes for data corresponding to Fig. 5. (10a) Re vs Pr. The black dashed line corresponds to Re ~ $1/\sqrt{Pr}$, corresponding to ultimate regime scaling. (10b) Compensated plot Re/(1/√Pr) vs Pr.</figcaption></figure>

Figures 7 and 8 presents the 3D Nu vs Ra, Pr scalings. Figures 9 and 10 presents the 3D Re vs Ra, Pr scalings. Scaling are always displayed both directly and in compensated form.

At low Ra, we first observe a transition from a laminar regime, where Nu = 1 up to a turbulent regime starting around Ra ~ $10^7$ at Pr = 1. In this transition regime, the Nusselt number varies approximately like Nu ~ $Ra^{2/3}$, while the Reynolds number remains less than $10^4$, but follows approximates laws Re ~ $Ra^{1/2}$. In this regime, the friction is negligible, as we saw, so that it corresponds to a laminar, frictionless regime.

After this laminar regime, we obtain a turbulent regime around $10^7 < Ra$ for Pr = 1 in
<figure><img src="image_12.png" /><figcaption>**Figure 11:** Scaling of heat transfer Nu as a function of Prandtl number Pr in 3D results at fixed Ra, dataset VII (Table 3) (11a) $Nu/\sqrt{Pr}$ vs Pr for various Ra. (11b) $Re/(1/\sqrt{Pr})$ vs Pr for various Ra.</figcaption></figure>

which Nu ~ Ra¹/² and Re ~ Ra¹/², like GL theory. The exact value of the exponent is provided in Table 1. In this regime, the friction are non-negligible, so that it is a “turbulent friction dominated regime” However, as both ratio $f_u = fU_{ls}^2/\epsilon_u$ or $f_\theta = f\Theta_{ls}^2/\epsilon_\theta$ remain independent of Ra, they do not change the scaling of the total kinetic and thermal energy dissipation. Therefore, the argument developed by GL theory should still apply in this situation, as is indeed observed, with minor corrections due to the small variations of the ratios.

In that respect, it is not surprising that the extent of this regime varies with Pr, as is shown in Fig. 11 for various Ra. At Ra = 10⁸, the “universal GL” regime stops for Pr < ~ 10⁻¹. In this range of parameters, Re is still large, so that the flow is turbulent. However, Nu drops quicker with decreasing Pr than in the universal GL regime, as can be seen from the filled data points in Fig. 8, in parallel with a similar drop for the thermal friction observed in Fig. 5b. This regime seems therefore dependent of the variation on the friction, and is non-universal. In this regime, the Reynolds number variation with Pr is milder than in the universal regime, as can be seen in Fig. 10.

As the Rayleigh number increases, we nevertheless observe in Fig. 11 that the extent of the universal turbulent regime extends towards smaller and smaller values of Pr, so that the universal scaling regime corresponds to an “asymptotic scaling regime” at low value of Pr < 1, valid in the limit of infinite Ra.

Figures 12 and 13 plot the kinetic and thermal dissipation rates $\epsilon_u, \epsilon_\theta$ against GL predictions. In agreement with what has been observed previously, we observe agreement with GL theory in the range of parameters where the friction ratios are approximately constant with the parameters, i.e. at large value of Re Pr. Overall, it is interesting to note that even when the friction is dominant, we can recover the ultimate regime scaling, as long as the velocity friction ratio remain relatively constant as a function of the parameters and neither there is not too big an asymmetry between the two frictions. In regimes where the asymmetry prevails, there are no clear scaling laws that emerge, meaning that the scaling are probably not universal in Ra and Pr only, and that friction depending corrections need to be implemented.
<figure><img src="image_13.png" /><figcaption>**Figure 12:** Scaling of thermal dissipation rate $\epsilon_{\theta}$ compared to the GL prediction Re$\sqrt{Pr/Ra}$ in 3D results. Correspondence between symbols and datasets are given in Table 3. The gray dashed line separates the non-universal (left) and the universal (right) friction dominated regimes for data corresponding to Fig. 5. (12a) $\epsilon_{\theta}\sqrt{RaPr}$ vs RePr. The black dashed line corresponds to the GL prediction $\epsilon_{\theta} \sim \text{Re}(Pr/Ra)^{1/2}$. (12b) Compensated plot $\epsilon_{\theta}\sqrt{RaPr}/\sqrt{RePr}$ vs $\sqrt{RePr}$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13:** Scaling of kinetic dissipation rate $\epsilon_{u}$ compared to the GL prediction Re$^{3}$(Pr/Ra)$^{3/2}$ in 3D results. Correspondence between symbols and datasets are given in Table 3). The gray dashed line separates the non-universal (left) and the universal (right) friction dominated regimes for data corresponding to Fig. 5. (13a) $\epsilon_{u}\sqrt{Ra^{3}Pr}$ vs Re$^{3}Pr^{2}$. The black dashed line corresponds to the GL prediction $\epsilon_{u} \sim \text{Re}^{3}(Pr/Ra)^{3/2}$. (13b) Compensated plot $\epsilon_{u}\sqrt{Ra^{3}Pr}/\sqrt{Re^{3}Pr^{2}}$ vs Re$^{3}Pr^{2}$.</figcaption></figure>

# 4. Conclusion

In this letter, we investigated scaling laws in the Homogeneous Rayleigh Bénard (HRB) equations through a new mathematical framework (“log-lattice”). Using a modified DOPRI solver, we are able to explore a range of parameters and wave-numbers way beyond what is accessible in direct numerical simulations of the equations. By adding a large-scale friction to the HRB equations, we are able to solve the issue of exponentially diverging solutions. This large scale friction becomes negligible when the fluid becomes turbulent enough, so that total energy balance depart from the energy balance considered in GL theory, where no
366 friction is present. Despite this, we still observe scaling law for Nu and Pr that are very close
to the universal turbulent predictions of Grossmann-Lohse (GL) theory: Nu ~ Ra¹/² Pr¹/²,
367
368 Re ~ Ra¹/² Pr⁻¹/², εθ ~ Re(Pr/Ra)¹/², εu ~ Re³(Pr/Ra)³/² for an important range of
parameters, corresponding to situations where the thermal friction is non-negligible and the
369
370 kinetic friction does not vary significantly as a function of the parameters. This is obtained
371 at large enough Ra and for Pr depending on the value of Ra.

372 In addition to this regime, we also observe another turbulent friction dominated regime
373 at Pr ≪ 1. This regime has no simple and universal dependence with the parameter, and
374 depends on the variations of the kinetic friction with the parameters.

375 Our observation show that the inclusion of friction, which is necessary to obtain stationary
376 regimes in the HRB framework, complexifies the phase space but nevertheless allows for
377 the existence of a universal turbulent regime, where scaling laws are very close to the GL
378 friction-less theoretical laws. In some geophysical or astrophysical situations, large scale
379 friction arises due to rotation (Ekman friction), stratification (Rayleigh friction) or magnetic
380 field (Hartman friction), and the two scaling regimes we find (one universal, and one non-
381 universal) may be relevant and could be explored within the log-lattice framework.

382 More generally, we believe that log-lattices, with their unique performances in terms of
383 numerical complexity, due to a spectrally sparse representation and strong mathematical
384 qualities, have a great potential in numerical simulations of geophysical or astrophysical
385 flows. However, as they are still in their infancy, many different paths would benefit from
386 being explored to better understand their strengths and weaknesses. This in particular includes
387 a better understanding of the influence of the numerical scheme which, as discussed in the
388 supplementary materials, may misrepresent the viscosity at high wavenumbers. We believe
389 that methods such as Whalen et al. (2015) could prove most useful in that regard. Other topics
390 of interest include the behavior of observables when $\lambda \to 1$ and the addition of the $k=0$
391 mode would prove very interesting to study. Likewise, in a similar spirit as was done for
392 the REWA model in Grossmann et al. (1996), a detailed comparison of DNS and log-lattice
393 results (which is far from trivial, as there is room for interpretation as to the mathematical
394 meaning of the fields simulated on a log-lattice) would be highly useful.

395 **Acknowledgements.** We thank A. Mailybaev and C. Campolina for initiating us to the log-lattice computations, and for stimulating discussions.

397 **Funding.** This work received funding from the Ecole Polytechnique, from ANR EXPLOIT, grant agreement no. ANR-16-CE06-0006-01 and ANR TILT grant agreement no. ANR-20-CE30-0035.

399 **Declaration of interests.** The authors report no conflict of interest.

REFERENCES

400 AHLERS, GUENTER, GROSSMANN, SIEGFRIED & LOHSE, DETLEF 2009 Heat transfer and large scale dynamics in turbulent Rayleigh-Bénard convection. *Rev. Mod. Phys.* **81**, 503–537.

401 BIFERALE, LUCA 2003 Shell models of energy cascade in turbulence. *Annual Review of Fluid Mechanics* **35**, 441–468.

404 BORUE, VADIM & ORSZAG, STEVEN A 1997 Turbulent convection driven by a constant temperature gradient. *Journal of Scientific computing* **12** (3), 305–351.

406 BOUILLAUT, VINCENT, LEPOT, SIMON, AUMÂTRE, SÉBASTIEN & GALLET, BASILE 2019 Transition to the ultimate regime in a radiatively driven convection experiment. *Journal of Fluid Mechanics* **861**, R5.

408 BRANDENBURG, AXEL 1992 Energy spectra in a model for convective turbulence. *Phys. Rev. Lett.* **69**, 605–608.

410 CALZAVARINI, E., DOERING, C. R., GIBBON, J. D., LOHSE, D., TANABE, A. & TOSCHI, F. 2006 Exponentially growing solutions in homogeneous Rayleigh-Bénard convection. *Phys. Rev. E* **73**, 035301.

412 CALZAVARINI, E., LOHSE, D. & TOSCHI, F. 2007 Homogeneous Rayleigh-Bénard Convection. In *Progress in*
413 *Turbulence II* (ed. Martin Oberlack, George Khujadze, Silke Günther, Tanja Weller, Michael Frewer, Joachim Peinke & Stephan Barth), pp. 181–184. Berlin, Heidelberg: Springer Berlin Heidelberg.

415 CALZAVARINI, ENRICO, LOHSE, DETLEF, TOSCHI, FEDERICO & TRIPICCIONE, RAFFAELE 2005 Rayleigh and Prandtl number scaling in the bulk of Rayleigh-Bénard turbulence. *Physics of Fluids* **17** (5), 055107, arXiv: https://doi.org/10.1063/1.1884165.

418 CAMPOLINA, CIRO S. & MAILYBAEV, ALEXEI A. 2018 Chaotic Blowup in the 3D Incompressible Euler Equations on a Logarithmic Lattice. *Phys. Rev. Lett.* **121**, 064501.

420 CAMPOLINA, CIRO S & MAILYBAEV, ALEXEI A 2021 Fluid dynamics on logarithmic lattices. *Nonlinearity* **34** (7), 4684–4715.

422 CASTAING, B., RUSAOUEN, E., SALORT, J. & CHILLA, F. 2017 Turbulent heat transport regimes in a channel. *Phys. Rev. Fluids* **2**, 062801.

424 CHAVANNE, X., CHILLÀ, F., CASTAING, B., HÉBRAL, B., CHABAUD, B. & CHAUSSY, J. 1997 Observation of the Ultimate Regime in Rayleigh-Bénard Convection. *Phys. Rev. Lett.* **79**, 3648–3651.

426 CHING, EMILY S. C. & KO, T. C. 2008 Ultimate-state scaling in a shell model for homogeneous turbulent convection. *Phys. Rev. E* **78**, 036309.

428 CILIBERTO, S. & LAROCHE, C. 1999 Random Roughness of Boundary Increases the Turbulent Convection Scaling Exponent. *Phys. Rev. Lett.* **82**, 3998–4001.

430 DOERING, CHARLES R. & CONSTANTIN, PETER 1996 Variational bounds on energy dissipation in incompressible flows. III. Convection. *Phys. Rev. E* **53**, 5957–5981.

432 FRISCH, URIEL, POMYALOV, ANNA, PROCACCIA, ITAMAR & RAY, SAMRIDDHI SANKAR 2012 Turbulence in Noninteger Dimensions by Fractal Fourier Decimation. *Phys. Rev. Lett.* **108**, 074501.

434 GLOAGUEN, C., LÉORAT, J., POUQUET, A. & GRAPPIN, R. 1985 A scalar model for MHD turbulence. *Physica D: Nonlinear Phenomena* **17** (2), 154–182.

436 GROSSMANN, SIEGFRIED & LOHSE, DETLEF 2000 Scaling in thermal convection: a unifying theory. *Journal of Fluid Mechanics* **407**, 27–56.

438 GROSSMANN, SIEGFRIED & LOHSE, DETLEF 2011 Multiple scaling in the ultimate regime of thermal convection. *Physics of Fluids* **23** (4), 045108, arXiv: https://doi.org/10.1063/1.3582362.

440 GROSSMANN, SIEGFRIED, LOHSE, DETLEF, L'VOV, VICTOR & PROCACCIA, ITAMAR 1994 Finite size corrections to scaling in high Reynolds number turbulence. *Phys. Rev. Lett.* **73**, 432–435.

442 GROSSMANN, SIEGFRIED, LOHSE, DETLEF & REEH, ACHIM 1996 Developed Turbulence: From Full Simulations to Full Mode Reductions. *Phys. Rev. Lett.* **77**, 5369–5372.

444 JIANG, HECHUAN, WANG, DONGPU, LIU, SHUANG & SUN, CHAO 2022 Experimental Evidence for the Existence of the Ultimate Regime in Rapidly Rotating Turbulent Thermal Convection. *Phys. Rev. Lett.* **129**, 204502.

447 KAWANO, KOKI, MOTOKI, SHINGO, SHIMIZU, MASAKI & KAWAHARA, GENTA 2021 Ultimate heat transfer in 'wall-bounded' convective turbulence. *Journal of Fluid Mechanics* **914**, A13.

449 KRAICHNAN, ROBERT H. 1962 Turbulent Thermal Convection at Arbitrary Prandtl Number. *The Physics of Fluids* **5** (11), 1374–1389, arXiv: https://aip.scitation.org/doi/pdf/10.1063/1.1706533.

451 LANOTTE, ALESSANDRA S., BENZI, ROBERTO, MALAPAKA, SHIVA K., TOSCHI, FEDERICO & BIFERALE, LUCA 2015 Turbulence on a Fractal Fourier Set. *Phys. Rev. Lett.* **115**, 264502.

453 LEPOT, SIMON, AUMAÎTRE, SÉBASTIEN & GALLET, BASILE 2018 Radiative heating achieves the ultimate regime of thermal convection. *Proceedings of the National Academy of Sciences* **115** (36), 8937–8941, arXiv: https://www.pnas.org/doi/pdf/10.1073/pnas.1806823115.

456 LOHSE, DETLEF & TOSCHI, FEDERICO 2003 Ultimate State of Thermal Convection. *Phys. Rev. Lett.* **90**, 034502.

458 MALKUS, W. V. R. & CHANDRASEKHAR, SUBRAHMANYAN 1954 The heat transport and spectrum of thermal turbulence. *Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences* **225** (1161), 196–212, arXiv: https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.1954.0197.

462 MOTOKI, SHINGO, KAWAHARA, GENTA & SHIMIZU, MASAKI 2022 Steady thermal convection representing the ultimate scaling. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences* **380** (2225), 20210037, arXiv: https://royalsocietypublishing.org/doi/pdf/10.1098/rsta.2021.0037.

466 PAWAR, SHASHIKANT S. & ARAKERI, JAYWANT H. 2016 Two regimes of flux scaling in axially homogeneous turbulent convection in vertical tube. *Phys. Rev. Fluids* **1**, 042401.

468 PUMIR, ALAIN & SHRAIMAN, BORIS I. 1995 Persistent Small Scale Anisotropy in Homogeneous Shear Flows. *Phys. Rev. Lett.* **75**, 3114–3117.
470 ROCHE, PHILIPPE-E 2020 The ultimate state of convection: a unifying picture of very high Rayleigh numbers experiments. *New Journal of Physics* **22** (7), 073056.

471 RUSAOUËN, E., LIOT, O., CASTAING, B., SALORT, J. & CHILLA, F. 2018 Thermal transfer in Rayleigh-Bénard cell with smooth or rough boundaries. *Journal of Fluid Mechanics* **837**, 443–460.

474 SCHMIDT, LAURA E., CALZAVARINI, ENRICO, LOHSE, DETLEF, TOSCHI, FEDERICO & VERZICCO, ROBERTO 2012 Axially homogeneous Rayleigh-Bénard convection in a cylindrical cell. *Journal of Fluid Mechanics* **691**, 52–68.

477 SPIEGEL, EDWARD A 1963 A Generalization of the Mixing-Length Theory of Turbulent Convection. *The Astrophysical Journal* **138**, 216.

479 STEVENS, BJORN, DUAN, JIANJUN, MCWILLIAMS, JAMES C., MÜNNICH, MATTHIAS & NEELIN, J. DAVID 2002 Entrainment, Rayleigh Friction, and Boundary Layer Winds over the Tropical Pacific. *Journal of Climate* **15** (1), 30–44.

482 SUKORIANSKY, SEMION, GALPERIN, BORIS & CHEKHLOV, ALEXEI 1999 Large scale drag representation in simulations of two-dimensional turbulence. *Physics of Fluids* **11** (10), 3043–3053, arXiv: https://doi.org/10.1063/1.870163.

485 URBAN, P., HANZELKA, P., KRÁLÍK, T., MACEK, M., MUSILOVÁ, V. & SKRBEK, L. 2019 Elusive transition to the ultimate regime of turbulent Rayleigh-Bénard convection. *Phys. Rev. E* **99**, 011101.

487 WHALEN, PATRICK, BRIO, MOYSEY & MOLONEY, JEROME V 2015 Exponential time-differencing with embedded Runge-Kutta adaptive step control. *Journal of Computational Physics* **280**, 579–601.

488 YEUNG, P. K., SREENIVASAN, K. R. & POPE, S. B. 2018 Effects of finite spatial and temporal resolution in direct numerical simulations of incompressible isotropic turbulence. *Phys. Rev. Fluids* **3**, 064603.

491 ZHU, XIAOJUE, MATHAI, VARGHESE, STEVENS, RICHARD JAM, VERZICCO, ROBERTO & LOHSE, DETLEF 2019a Reply to "Absence of Evidence for the Ultimate Regime in Two-Dimensional Rayleigh-Benard Convection". *Physical Review Letters* **123** (25), 259402.

494 ZHU, XIAOJUE, MATHAI, VARGHESE, STEVENS, RICHARD J. A. M., VERZICCO, ROBERTO & LOHSE, DETLEF 2018 Transition to the Ultimate Regime in Two-Dimensional Rayleigh-Bénard Convection. *Phys. Rev. Lett.* **120**, 144502.

497 ZHU, XIAOJUE, STEVENS, RICHARD J. A. M., SHISHKINA, OLGA, VERZICCO, ROBERTO & LOHSE, DETLEF 2019b $Nu \sim Ra^{1/2}$ scaling enabled by multiscale wall roughness in Rayleigh-Bénard turbulence. *Journal of Fluid Mechanics* **869**, R4.

500 ZOU, SHUFAN & YANG, YANTAO 2021 Realizing the ultimate scaling in convection turbulence by spatially decoupling the thermal and viscous boundary layers. *Journal of Fluid Mechanics* **919**, R3.