Rapidly rotating radiatively driven convection: experimental and numerical validation of the ‘geostrophic turbulence’ scaling predictions

Abstract

We experimentally and numerically characterize rapidly rotating radiatively driven thermal convection, beyond the sole heat transport measurements reported by Bouillaut et al. (Proc. Natl Acad. Sci., vol. 118, 2021, e2105015118). Based on a suite of direct numerical simulations (DNS) and additional processing of the experimental data collected by Bouillaut et al. (Proc. Natl Acad. Sci., vol. 118, 2021, e2105015118), we report the simultaneous validation of the scaling predictions of the ‘geostrophic turbulence’ regime – the diffusivity-free or ‘ultimate’ regime of rapidly rotating convection – for the heat transport and the temperature fluctuations. Following such cross-validation between DNS and laboratory experiments, we further process the numerical data to validate the ‘geostrophic turbulence’ scaling predictions for the flow velocity and horizontal scale. Radiatively driven convection thus appears as a versatile set-up for the laboratory observation of the diffusivity-free regimes of various convective flows of geophysical and/or astrophysical interest.

1. Introduction

Thermal convection is a key process driving natural turbulent flows, be it in stellar and planetary interiors or in the open ocean (Stevenson 1979; Marshall & Schott 1999; Aurnou et al. 2015; de Jong & de Steur 2016; Hindman, Featherstone & Julien 2020; Vasil, Julien & Featherstone 2021). In these various contexts, global rotation strongly affects the resulting flow and its transport properties through the action of the Coriolis force (Ecke & Shishkina 2023). The Reynolds number takes prohibitively large values in natural settings, beyond direct reach of laboratory and numerical models. The standard approach thus consists in inferring scaling laws based on dimensional analysis, asymptotic theory, experimental studies and numerical studies, with the goal of extending these scaling laws to the extreme parameter values that characterize natural flows.

The scaling relation between the convective heat flux and the overall temperature drop has received most attention, as it characterizes both the heat transport properties of the system and the energetics of the flow. When it comes to the modelling of natural flows, a general belief is that the tiny molecular values of the thermal and momentum diffusivities of the system should not enter this scaling relation (Spiegel 1971; Marshall & Schott 1999). Such a ‘diffusivity-free’ argument can be traced back to the zeroth law of turbulence, understood in a broad sense: in a fully turbulent system, the large-scale quantities should be related in a way that does not involve the tiny molecular diffusivities.

For non-rotating convection, the diffusivity-free argument alone allows one to conclude on the scaling relation between the dimensionless heat flux (the Nusselt number $Nu$) and the dimensionless temperature drop (the temperature-based Rayleigh number $Ra^{(\Delta T)}$). This leads to the so-called ‘ultimate’ scaling regime of thermal convection (Kraichnan 1962; Spiegel 1963). By contrast, rotating convection involves an additional dimensionless parameter (the Ekman number $E$) and the diffusivity-free argument alone does not suffice to predict the scaling behaviour of the system. Stevenson (1979) therefore introduced additional scaling arguments, many of which were inferred (perhaps surprisingly) from the linear stability analysis of rapidly rotating convection. Combining these additional arguments with the diffusivity-free assumption, Stevenson (1979) obtained definite predictions for the scaling behaviour of the main quantities of interest: heat transport, temperature fluctuations, flow speed and characteristic scales. A more mathematically rooted way to arrive at the same scaling predictions stems from the work of Julien, Knobloch & Werne (1998), who carried out an asymptotic expansion of the governing equations in the rapidly rotating regime. This procedure results in a reduced set of equations that involves fewer dimensionless control parameters: instead of the Rayleigh number $Ra^{(\Delta T)}$ and the Ekman number $E$ arising independently in the equations, the reduced set of equations involves only a reduced Rayleigh number $\widetilde {Ra}=Ra^{(\Delta T)} E^{4/3}$. Provided one starts from the reduced set of equations of Julien et al. (1998), the sole diffusivity-free assumption yields all the scaling predictions of Stevenson (1979), confirming Stevenson's impressive physical insight. The corresponding fully turbulent scaling regime of rapidly rotating convection is referred to as the ‘Geostrophic Turbulence’ (GT) scaling regime in the more recent literature (Cheng et al. 2018).

An ongoing line of research consists in trying to reach the GT regime using Direct Numerical Simulation (DNS). This approach has proven successful in a range of idealized models of increasing numerical complexity. DNS of the reduced equations (Julien et al. 2012; Plumley et al. 2016) point to the GT regime at large reduced Rayleigh number. DNS of the full Boussinesq equations for the Rayleigh–Bénard set-up with free-slip boundary conditions also point to the GT scaling regime for the heat transport (Stellmach et al. 2014). Departing from the standard Rayleigh–Bénard set-up, Barker, Dempsey & Lithwick (2014) introduced internal heating and cooling in the vicinity of the lower and upper boundaries, respectively, and validated the various scaling predictions of the GT regime (see also Currie et al. 2020). More recently, Song, Shishkina & Zhu (2024) have reported DNS of Rayleigh–Bénard convection with no-slip boundaries; the heat transport scaling law appears to be compatible with the GT prediction for the lowest value of the Ekman number considered in their study, while the GT scaling predictions for the flow speed and flow structure appear to be satisfied over a more extended region of parameter space.

At the experimental level, there are challenging constraints to observe the GT regime in the laboratory. Achieving high Rayleigh number and low Ekman number suggests building tall Rayleigh–Bénard cells while keeping the radius of the tank small enough to avoid centrifugal effects (Ecke & Niemela 2014; Cheng et al. 2015, 2018, 2020; Zhang et al. 2020; Wedi et al. 2021). An issue with such cigar-shaped containers is the emergence of wall modes and boundary zonal flows that can contribute significantly to the overall heat transport (Favier & Knobloch 2020; Zhang, Ecke & Shishkina 2021; Wedi et al. 2022). It is therefore unclear whether the GT regime can be achieved in the current generation of rotating Rayleigh–Bénard experiments, and determining the optimal laboratory design that could potentially lead to the observation of the GT regime in rotating Rayleigh–Bénard convection remains an intense topic of ongoing research (Cheng et al. 2018; Kunnen 2021; Terrien, Favier & Knobloch 2023).

Departing from the standard Rayleigh–Bénard set-up, we have introduced an experimental apparatus where convection is driven radiatively, through the absorption of visible light by a dyed fluid. Through a combination of radiative heating and effective internal cooling, one can bypass the top and bottom boundary layers and observe regimes of thermal convection that are controlled by the bulk turbulent flow, such as the ‘ultimate’ regime of (non-rotating) turbulent convection (Lepot, Aumaître & Gallet 2018; Bouillaut et al. 2019; Miquel et al. 2019, 2020). Global rotation was recently added to the laboratory set-up, providing the first experimental observation of the GT heat transport scaling relation (Bouillaut et al. 2021). This experimental configuration shares many similarities with the numerical set-up introduced by Barker et al. (2014) (see also Currie et al. 2020), differing primarily in the distribution of internal heat sources and sinks, in the boundary conditions, and in the precise definitions of the diagnostic variables. The goal of the present study is to characterize the turbulent state of radiatively driven rotating convection beyond the sole heat transport, and assess the validity of the GT scaling predictions for the temperature fluctuations, the flow speed and the characteristic scale of the flow. One issue at the experimental level is that the fluid is opaque, making global velocity measurements challenging (see Bouillaut et al. 2022 for velocity estimates in the non-rotating experiment). We thus carried out a hybrid experimental–numerical study of radiatively driven convection combining further processing of the temperature measurements in Bouillaut et al. (2021) with in-depth numerical diagnostics of the speed and structure of the flow. We validate the numerical approach as a good model of the system by intercomparison of the experimental and numerical data for temperature, before leveraging the three-dimensional (3-D) DNS to characterize the velocity field.

2. Radiatively driven rotating convection

2.1. Theoretical set-up

We consider a horizontal fluid layer of height $H$ radiatively heated from below and thermally insulated at all boundaries (see figure 1 for a schematic of the laboratory implementation). The fluid absorbs light at the same rate for all incoming wavelengths, which leads to a typical absorption length $\ell$. Following Beer–Lambert's law and denoting as $P$ the radiative flux (per unit surface) impinging on the bottom boundary $z=0$ of the fluid layer, the absorption of light induces an internal heat source that decreases exponentially with height $z$:

(2.1)\begin{equation} Q_H(z)=\frac{P}{\ell} \,{\rm e}^{{-}z/\ell} . \end{equation}

Together with such radiative heating, we consider a uniform internal cooling term $Q_C=-(1-{\rm e}^{-H/\ell })P/H < 0$. The volume integral of this cooling term over the fluid domain is opposite to that of the heating term. That is, the overall heat input by the radiative heat source is exactly balanced by the overall heat removed by the uniform heat sink, which ensures a statistically stationary state for the temperature and velocity fields (see § 2.2 for the experimental implementation of the effective cooling term).

Figure 1. (a) Experimental and (b) numerical implementations of radiatively driven rotating convection. Absorption of an upward flux of light by the dyed fluid induces an internal heat source that decreases exponentially with height measured from the bottom of the fluid domain, over an absorption length $\ell$. Secular heating of the fluid induces uniform effective internal cooling compensating the radiative heat source on vertical average.

The fluid layer rotates at a rate $\varOmega \boldsymbol {e}_z$ with respect to an inertial frame and is subject to gravity $-g\boldsymbol {e}_z$. We restrict attention to the range of parameters where the centrifugal acceleration is negligible (Horn & Aurnou 2018, 2019). Following the Boussinesq approximation (Spiegel & Veronis 1960), the fluid properties are assumed to be constant and uniform, with the exception of the density $\rho$, whose variations are retained in the buoyancy force only and are assumed to vary linearly with temperature. The velocity field is thus divergence free. Denoting the reduced pressure field as $p({\boldsymbol x},t)$, the temperature variable as $\theta ({\boldsymbol x},t)$ and the velocity field in the rotating frame as $\boldsymbol {u}({\boldsymbol x},t)$, the governing equations read:

(2.2)\begin{gather} \partial_t\boldsymbol{u} + (\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla})\boldsymbol{u} +2\varOmega \boldsymbol{e}_z \times \boldsymbol{u} ={-}\boldsymbol{\nabla}p + \alpha g \theta \boldsymbol{e}_z + \nu \nabla^2\boldsymbol{u} , \end{gather}

(2.3)\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0 , \end{gather}

(2.4)\begin{gather}\partial_t \theta + \boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}\theta = \kappa \nabla^2 \theta + \frac{Q_H(z)+Q_C}{\rho C} , \end{gather}

where $C$ denotes the specific heat capacity of the fluid, $\alpha$ the coefficient of thermal expansion, $\rho$ the reference density, $\kappa$ the thermal diffusivity and $\nu$ the kinematic viscosity.

The set of equations above involves four dimensionless control parameters:

(2.5a–d)\begin{equation} Ra^{(P)}=\frac{\alpha g P H^4}{\rho C \kappa^2 \nu},\quad E=\frac{\nu}{{2}\varOmega H^2},\quad Pr=\frac{\nu}{\kappa}\quad \text{and} \quad \tilde{\ell}=\frac{\ell}{H} , \end{equation}

where the flux-based Rayleigh number $Ra^{(P)}$ characterizes the strength of the thermal forcing (strength of the imposed radiative heat flux), the Ekman number $E$ characterizes the strength of the global rotation, with $E \ll 1$ for rapid rotation, the Prandtl number $Pr$ characterizes the relative magnitudes of momentum and thermal diffusivities, and the dimensionless absorption length $\tilde {\ell }$ characterizes the spatial structure of the radiative heat source.

We wish to characterize the statistically steady state of the system. Most studies focus primarily on the overall heat transport properties of the system: how large is the overall temperature drop $\Delta T$ that emerges across the fluid layer as a result of the radiative heating? That is, for a given value of the flux-based Rayleigh number $Ra^{(P)}$, one would like to determine the emergent temperature-based Rayleigh number $Ra^{(\Delta T)}$, or, equivalently, the Nusselt number $Nu$, defined as

(2.6a,b)\begin{equation} Ra^{(\Delta T)}=\frac{\alpha g \Delta T H^3}{\nu \kappa}, \quad Nu =\frac{Ra^{(P)}}{Ra^{(\Delta T)}} = \frac{PH}{\rho C \kappa \Delta T} . \end{equation}

The goal of the present study is to go beyond the sole quantification of the overall heat transport, with a more in-depth characterization of the temperature and velocity fields: what is the typical flow speed? The typical horizontal scale $\ell ^\perp$ of the flow? How do local temperature fluctuations compare to the mean temperature drop? With these questions in mind we introduce the Reynolds number $Re$, the fluctuation-based Rayleigh number $Ra^{(\theta )}$ and the dimensionless horizontal length scale $\ell _*^\perp$, defined as

(2.7a–c)\begin{equation} Re = \frac{\sqrt{\langle \overline{\boldsymbol{u}^2} \rangle}H}{\nu} , \quad Ra^{(\theta)}=\frac{\alpha g \theta_{std} H^3}{\nu \kappa} , \quad \ell_*^\perp =\frac{\ell^\perp}{H} , \end{equation}

where $\langle {\cdot }\rangle$ denotes space average, $\bar {\cdot }$ denotes time average and $\theta _{std}=\sqrt {\overline {(\theta -\bar {\theta })^2}}$ denotes the standard deviation of $\theta$ at a given $z$ (see below).

2.2. Laboratory implementation

The experimental set-up, introduced by Bouillaut et al. (2021), is a rotating version of the radiatively driven convection set-up designed by Lepot et al. (2018). A schematic is provided in figure 1. A cylindrical tank rotating at a rate $\varOmega$ around the vertical axis contains a mixture of water and dye. A powerful spotlight shines at the tank from below. The light passes through a layer of cool water filtering out infrared radiation before reaching the transparent bottom plate of the tank. Absorption of visible light by the dye induces an internal heat source of the form (2.1), where the absorption length $\ell$ is directly set by the (uniform) concentration of the dye. The dimensionless absorption length is kept constant in the present study. Specifically, we adopt $\ell =0.048H$ which proves sufficient to bypass the throttling bottom boundary layers and induce diffusivity-free regimes of thermal convection (Lepot et al. 2018; Bouillaut et al. 2019, 2021).

While there is no cooling mechanism in the experimental apparatus, secular heating of the body of fluid amounts to an effective uniform heat sink denoted as $Q_C$ above. Indeed, consider the radiatively heated set-up in the absence of cooling. That is, we consider (2.4) with $Q_C$ set to zero, and we denote the temperature field as $T({\boldsymbol x},t)$ instead of $\theta ({\boldsymbol x},t)$. Space integration over the entire body of fluid, using insulated boundary conditions at all boundaries, yields

(2.8)\begin{equation} \frac{{\rm d} \langle T \rangle}{{\rm d}t}=\frac{P}{\rho CH}(1-{\rm e}^{{-}H/\ell}), \end{equation}

which indicates that the space-averaged temperature $\langle T \rangle$ increases linearly with time. This result holds as long as heat losses through the boundaries of the tank are negligible, a valid approximation when $\langle T \rangle$ is within a few degrees of room temperature. Now, introduce the variable $\theta ({\boldsymbol x},t)=T({\boldsymbol x},t)-\langle T \rangle (t)$. One can check that the evolution equation for $\theta$ is precisely (2.4), with both the heating term $Q_H$ and the effective cooling term $Q_C$ included. Additionally, the buoyancy force in (2.2) can be equivalently cast as $\alpha g T {\boldsymbol e}_z$ or as $\alpha g \theta {\boldsymbol e}_z$, the difference between the two being absorbed by the pressure gradient. We conclude that (2.2)–(2.4) indeed model the convective dynamics arising in the laboratory set-up.

We leave a free surface at the top of the dyed water, which allows us to easily vary the fluid height $H$ between $10$ and $25$ cm. The second dimensional control parameter is the rotation rate $\varOmega$, which we vary between $10$ and $85$ rpm. The corresponding region of the dimensionless parameter space is shown in figure 2, and we checked in Bouillaut et al. (2021) that centrifugal effects do not impact the measurements.

Figure 2. Parameter space spanned by the present dataset. Circles are experimental data ($Pr \approx 7$) while triangles are DNS ($Pr=7$).

We characterize the temperature field using three thermocouples located at height $z=0$, $z=0.25H$ and $z=0.75H$ along the axis of the cylinder. An experimental run consists in filling the tank up to a height $H$ with dyed water at approximately 10 $^\circ$C, spinning the tank at a given rate $\varOmega$ for at least ten minutes to achieve solid-body rotation, and then turning on the spotlight. After some transient, the timeseries of all three probes exhibit a common linear drift at a rate given by the right-hand side of (2.8) (examples of timeseries are provided in Bouillaut et al. 2021). Following the analysis above, any temperature difference between two probes exhibits a statistically steady signal provided (i) the initial transient phase has decayed and (ii) the fluid temperature is reasonably close to room temperature (typically ${\pm }5\,^\circ {\rm C}$). Similarly, temperature fluctuations around the mean drift are statistically steady when conditions (i) and (ii) are met. The experimental data consist of the temperature drop $\Delta T$ measured by Bouillaut et al. (2021), which we complement with an estimate of the temperature fluctuations. To wit, we quantify the temperature fluctuations at $z=0.25H$ by subtracting the linear drift from the timeseries before computing the root-mean-square (r.m.s.) fluctuations of the resulting statistically steady signal. This leads to the standard deviation $\theta _{std}$ entering the definition of $Ra^{(\theta )}$.

2.3. Numerical implementation

To characterize the system beyond the sole quantities accessible in the laboratory experiment, we performed DNS of the governing equations (2.2)–(2.4) in a horizontally periodic domain using the pseudo-spectral solver Coral (Miquel 2021), validated against both analytical results (Miquel et al. 2020) and solutions computed with the Dedalus software (Burns et al. 2020). Coral employs a Chebyshev–Fourier spatial decomposition and an implicit–explicit time-stepping scheme. In the present study, the upper and lower boundaries are insulated ($\partial _z\theta =0$) and impenetrable ($w=0$). To replicate as closely as possible the experiment, the kinematic boundary conditions are no-slip at the bottom ($u=v=0$) and free-slip at the top ($\partial _z u=\partial _z v = 0$).

The suite of numerical simulations is focused on the region of parameter space where the experimental data from Bouillaut et al. (2021) point to the GT scaling regime, see figure 2. The Ekman number varies between $5\times 10^{-7}$ and $10^{-4}$, while the flux-based Rayleigh number varies between $10^{11}$ and $10^{12}$. The Prandtl number is $Pr=7$ for all the DNS. The aspect ratio of the numerical domain is in the range $[0.35,1]$, ensuring that, for each DNS run, at least four wavelengths of the most unstable mode fit along a horizontal direction of the numerical domain. The typical resolution ranges from $(N_x,N_y,N_z)=(150,150,256)$ for $Ra^{(P)}=10^{11}$ to $(N_x,N_y,N_z)=(300,300,512)$ for $Ra^{(P)}=10^{12}$, where $N_x$ and $N_y$ denote the number of Fourier modes in the horizontal directions, while $N_z$ denotes the number of Chebyschev polynomials in the vertical directions. The non-uniform Chebyschev grid is particularly welcome for the present problem, allowing us to have typically eight grid points inside the Ekman layer after de-aliasing.

Initial conditions are chosen as either small-amplitude noise, or a solution computed in a previous run. After the initial transient has subsided, we compute the emergent dimensionless parameters of interest by averaging over the statistically steady regime. Because of the periodic boundary conditions, the numerical system is invariant to translations in the horizontal directions. Based on this invariance, when performing a time average $\bar {q}$ of some quantity $q$ extracted from the DNS, we also include a horizontal area average to speed up convergence. Whenever possible, we extract the exact numerical counterpart of the experimental measurements: the temperature difference $\Delta T$ is computed by temporally and horizontally averaging $\theta$ at $z=0$ and $z=0.75H$, before subtracting the two. The fluctuation-based Rayleigh number $Ra^{(\theta )}$ is based on $\theta _{std}$ evaluated at $z=0.25H$. Beyond the experimentally measured quantities, DNS gives access to the full statistics of the temperature and velocity fields. We thus also extract the r.m.s. velocity $\langle \overline {{\boldsymbol u}^2} \rangle ^{1/2}$ over the entire fluid domain to compute the emergent Reynolds number $Re$. Finally, we extract a characteristic horizontal length scale $\ell ^\perp$ of the flow as

(2.9)\begin{equation} \ell^\perp=\left( \frac{\langle \overline{\psi^2} \rangle}{\langle \overline{ \boldsymbol{u}^2} \rangle} \right)^{1/2} , \end{equation}

where $\psi (x,y,z,t)$ denotes the so-called toroidal streamfunction, defined as minus the inverse horizontal Laplacian of the vertical vorticity (for a purely horizontal flow, this definition leads to the standard streamfunction of the two-dimensional (2-D) flow in a given constant-$z$ plane, see e.g. Julien & Knobloch (2007) for details on the toroidal/poloidal decomposition).

In figure 3, we show vertical slices of the temperature field for fixed $Ra^{(P)}=10^{12}$ and decreasing values of the Ekman number (increasing rotation rate). One observes thinner structures with strong vertical coherence as the rotation rate increases. As compared with standard rotating Rayleigh–Bénard convection, the top-down asymmetry of the present system is visible in the temperature snapshots. This asymmetry arises predominantly as a result of the asymmetry of the distribution of heat sources and sinks, and to a lesser extent as a result of the asymmetry between the (free-slip) top and (no-slip) bottom boundary conditions. One observes stronger vertical temperature variations in the vicinity of the heating region, with a quieter region in the upper half of the domain. The temperature fluctuations measured by the probe located at $z=0.75H$ are thus weaker than the fluctuations measured at $z=0.25H$, although we checked that they display the same scaling behaviour with the control parameters of the system. More generally, and anticipating the results in the next sections, the good agreement between the scaling laws measured in the present set-up and in the idealized top-down-symmetric set-up of Barker et al. (2014) indicates that the top-down asymmetry does not impact the scaling behaviour of the various quantities of interest in the GT regime.

Figure 3. Vertical temperature slices extracted from DNS with $Ra^{(P)}=10^{12}$. From left to right, the diffusivity-free flux-based Rayleigh number (3.12) is ${\mathcal {R}}=\{1.31\times 10^{-6}, 1.63\times 10^{-7}, 2.04\times 10^{-8},2.55\times 10^{-9}\}$, and the convective Rossby number (5.1) is $Ro=\{5.05\times 10^{-2}, 2.66\times 10^{-2}, 1.77\times 10^{-2}, 1.16\times 10^{-2}\}$. The flow develops thinner columnar structures as the rotation rate increases (left to right).

3. The geostrophic turbulence regime: theoretical background

3.1. Asymptotically fast rotation, $E\ll 1$

As discussed at the outset, deriving scaling predictions for the various emergent quantities of interest proves more challenging for rotating convection than for non-rotating convection, because of the additional dimensionless parameter $E$ in the former case. Progress can be made in the rapidly rotating regime $E \ll 1$, as initially proposed by Stevenson (1979) and put on firm analytical footing by Julien et al. (1998) through an asymptotic expansion of the equations in powers of $E^{1/3}$ (Aurnou, Horn & Julien 2020). As for any asymptotic expansion, a challenging part of the analysis consists in inferring the correct scaling of the various quantities of interest with the small parameter $E^{1/3}$. We propose here a recipe to determine these scalings based on the following observation: the reduced set of equations for rapid rotation must hold near the threshold of instability. That is, the sought scalings can be inferred from the structure of the most unstable eigenmode arising near the threshold for instability. This approach provides a shortcut to deriving the predictions of the geostrophic turbulence scaling regime. The first step consists in writing the scaling in $E$ of the various quantities characterizing the most unstable eigenmode computed through linear stability analysis (see Chandrasekhar 1961 for rotating Rayleigh–Bénard convection and Bouillaut (2022) for the rotating radiatively driven set-up). In the rapidly rotating regime $E \ll 1$, the marginal state is characterized by the following asymptotic relations:

(3.1a)\begin{gather} \begin{array}{r} \text{Threshold temperature-based }\\ \text{Rayleigh number:} \end{array}\quad Ra^{(\Delta T)}_c \sim E^{{-}4/3}, \end{gather}

(3.1b)\begin{gather}{\text{Growthrate (or angular frequency):}} \quad {\sigma \sim \frac{\kappa}{H^2} E^{{-}2/3}}, \end{gather}

(3.1c) \begin{gather}\text{Horizontal length scale:} \quad \ell^{{\perp}} \sim H E^{1/3}, \end{gather}

(3.1d) \begin{gather}\text{Vertical length scale:} \quad \ell_z \sim H, \end{gather}

(3.1e)\begin{gather}\begin{array}{r} \text{Ratio of temperature to velocity}\\ \text{fluctuations:} \end{array}\quad {\theta_{std}} \sim \frac{\Delta T H}{\kappa} \langle \overline{w^2} \rangle^{1/2} E^{2/3}, \end{gather}

(3.1f)\begin{gather}\begin{array}{r} \text{Relation between the various velocity}\\ \text{components:} \end{array}\quad \langle \overline{u^2} \rangle \sim \langle \overline{v^2} \rangle \sim \langle \overline{w^2} \rangle . \end{gather}

Based on the asymptotic behaviour of the temperature-based Rayleigh number, one introduces the reduced Rayleigh number $\widetilde {Ra}=Ra^{(\Delta T)} E^{4/3}$. We focus on the rapidly rotating near-threshold regime corresponding to the distinguished limit $E \to 0$ with fixed $\widetilde {Ra}={{O}}(E^0)$. Equation (3.1b) above indicates the behaviour of the growthrate of the most unstable mode in this distinguished limit, or simply the behaviour of the frequency of oscillation at threshold when convection arises through a Hopf bifurcation (Chandrasekhar 1961). More generally, in the distinguished limit of interest here, all the $\sim$ symbols in (3.1a)–(3.1f) can be replaced by an equals sign, at the expense of multiplying the right-hand side by a generic function ${\mathcal {F}}(\widetilde {Ra},Pr)$ (in the following, the symbol ${\mathcal {F}}$ denotes a generic functional dependence that a priori differs between successive equations). Indeed, at the instability threshold, $\widetilde {Ra}$ is constant and one must recover the scalings (3.1a)–(3.1f) in Ekman number.

The saturation level of the various fields is determined based on the dominant nonlinearity of the equations (Stevenson 1979). Here the nonlinearities are all of advective type, entering the expression of the total derivative as

(3.2)\begin{equation} \partial_t + (\boldsymbol{u}\boldsymbol{\cdot} \boldsymbol{\nabla}) = \partial_t + \boldsymbol{u}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla}_\perp+ w \partial_z . \end{equation}

In the linear regime of the instability (exponential growth), the term $\partial _t$ is of the order of the growth rate $\sigma$, while the advective nonlinearities are negligible. Saturation arises when these advective nonlinearities become comparable to $\sigma$. One easily checks using (3.1a)–(3.1f) that vertical advection is negligible as compared with horizontal advection:

(3.3)\begin{equation} w \partial_z \sim \frac{w}{H} \sim \frac{u_\perp}{H} \sim \frac{\ell^\perp}{H} \times \frac{u_\perp}{\ell^\perp} \sim E^{1/3} \boldsymbol{u}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla}_\perp\ll \boldsymbol{u}_\perp \boldsymbol{\cdot} \boldsymbol{\nabla}_\perp . \end{equation}

Saturation thus arises when $\boldsymbol {u}_\perp \boldsymbol {\cdot } \boldsymbol {\nabla }_\perp \sim \langle \overline {{\boldsymbol u}^2} \rangle ^{1/2}/\ell ^\perp$ becomes comparable to $\sigma$. Substitution of the scalings (3.1a)–(3.1f) yields the scaling for the velocity at saturation:

(3.4)\begin{equation} \langle \overline{{\boldsymbol u}^2} \rangle^{1/2} = \frac{\kappa}{H} E^{{-}1/3} {\mathcal{F}}(\widetilde{Ra},Pr) , \end{equation}

which we recast as a scaling for the Reynolds number:

(3.5)\begin{equation} \frac{\langle \overline{w^2} \rangle^{1/2} H}{\nu} \sim \frac{\langle \overline{{\boldsymbol u}^2} \rangle^{1/2} H}{\nu} = E^{{-}1/3} {\mathcal{F}}(\widetilde{Ra},Pr). \end{equation}

With the heat flux $P/(\rho C)$ scaling as the convective flux $\langle \overline {w \theta } \rangle \sim \langle \overline {w^2} \rangle ^{1/2} \theta _{std}$, one can deduce the scaling in Ekman number of all the quantities of interest listed above at finite distance from threshold based on the combination of (3.4) with (3.1):

(3.6a)\begin{gather} Nu = {\mathcal{F}}(\widetilde{Ra},Pr) , \end{gather}

(3.6b)\begin{gather}\frac{\ell^\perp}{H} = E^{1/3} {\mathcal{F}}(\widetilde{Ra},Pr), \end{gather}

(3.6c)\begin{gather}\frac{\ell_z}{H} = {\mathcal{F}}(\widetilde{Ra},Pr), \end{gather}

(3.6d)\begin{gather}\frac{\theta_{std}}{\Delta T} = E^{1/3} {\mathcal{F}}(\widetilde{Ra},Pr) , \end{gather}

where, again, the symbol ${\mathcal {F}}$ denotes a generic function that differs between successive equations. The equations above provide the dominant scalings in Ekman number behind the derivation of the reduced model by Julien et al. (1998).

3.2. Diffusivity-free regime

In practical terms, the distinguished limit considered above allows one to consider increasingly large values of $\widetilde {Ra}$ as the Ekman number $E$ decreases. For extremely low values of $E$, one can even consider values of $\widetilde {Ra}$ that are much greater than one. The leap of faith is then to assume that there is a regime of low-enough Ekman number for (3.5)–(3.6) to hold, but far enough from threshold ($\widetilde {Ra} \gg 1$) to reach a diffusivity-free scaling regime (numerical integration of the associated reduced model does provide evidence for this, see e.g. Julien et al. 2012). The diffusivity-free scaling argument then amounts to demanding that the functions ${\mathcal {F}}$ in (3.5)–(3.6d) be such that $\kappa$ and $\nu$ can be crossed out from both sides of the equations, which finally leads to the GT scaling predictions:

(3.7)\begin{gather} Nu \sim \frac{\widetilde{Ra}^{3/2}}{Pr^{1/2}} \sim (Ra^{(\Delta T)})^{3/2} E^2 Pr^{{-}1/2}, \end{gather}

(3.8)\begin{gather}\frac{\ell^\perp}{H} \sim E^{1/3} \sqrt{\frac{\widetilde{Ra}}{Pr}} \sim E \sqrt{\frac{Ra^{(\Delta T)}}{Pr}}, \end{gather}

(3.9)\begin{gather}\frac{\ell_z}{H} \sim 1, \end{gather}

(3.10)\begin{gather}\frac{\langle {\boldsymbol u}^2 \rangle^{1/2} H}{\nu} \sim E^{{-}1/3} \frac{\widetilde{Ra}}{Pr} \sim \frac{Ra^{(\Delta T)} \, E}{Pr}, \end{gather}

(3.11)\begin{gather}\frac{\theta_{std}}{\Delta T} \sim E^{1/3} \sqrt{\frac{\widetilde{Ra}}{Pr}} \sim E \sqrt{\frac{Ra^{(\Delta T)}}{Pr}} . \end{gather}

To highlight the diffusivity-free form of these predictions, we combine the emergent parameters of interest with $E$ and $Pr$ to form dimensionless combinations that do not involve the molecular diffusivities $\kappa$ and $\nu$. For radiatively driven convection, the central control parameter is then the diffusivity-free flux-based Rayleigh number ${\mathcal {R}}$, defined as

(3.12)\begin{equation} {\mathcal{R}}= \frac{Ra^{(P)} E^3}{Pr^2} . \end{equation}

Denoting with a star the diffusivity-free form of the emergent dimensionless parameters, we introduce

(3.13a–d)\begin{equation} Nu_* = \frac{Nu E}{Pr} , \quad Ra^{(\theta)}_*=\frac{Ra^{(\theta)} E^2}{Pr} , \quad Re_*=E\, Re , \quad \ell_*^\perp=\frac{\ell^\perp}{H} . \end{equation}

Such diffusivity-free dimensionless numbers allow us to recast the GT scaling predictions in a particularly compact form reported in table 1. Namely, each diffusivity-free emergent parameter evolves as some power-law in ${\mathcal {R}}$ with a specific exponent.

Table 1. Scaling predictions for the various emergent quantities of interest in the geostrophic turbulence (GT) scaling regime, expressed in terms of the diffusivity-free flux-based Rayleigh number ${\mathcal {R}}= {Ra^{(P)} E^3}/{Pr^2}$.

4. Experimental and numerical assessment of the scaling predictions

4.1. Heat transport

As discussed at the outset, a central question in turbulent convection is the scaling relation between the heat flux and the overall temperature drop. The GT scaling prediction for the heat transport has been previously validated in Bouillaut et al. (2021) by plotting the experimental data for $Nu_*$ as a function of ${\mathcal {R}}$. We reproduce this plot in figure 4 for completeness, to which we add the numerical data. We confirm the excellent agreement between the experimental data and the GT scaling prediction $Nu_* \sim {\mathcal {R}}^{3/5}$ for ${\mathcal {R}} \leq 3\times 10^{-7}$. The numerical data are also in excellent agreement with the prediction. The compensated plot further highlights this agreement and indicates that the prefactors of the scaling laws inferred from the experimental and numerical data are fully compatible, the numerical prefactor being greater by approximately $15$ %.

Figure 4. Diffusivity-free Nusselt number $Nu_*$ as a function of the diffusivity-free flux-based Rayleigh number ${\mathcal {R}}$. Same symbols as in figure 2. The dotted line indicates the GT scaling exponent $3/5$. The dashed line is the prediction of the non-rotating ‘ultimate’ regime, characterized by an exponent of $1/3$ in this representation. The inset shows $Nu_*/{\mathcal {R}}^{3/5}$ versus ${\mathcal {R}}$.

4.2. Temperature fluctuations

We now turn to the temperature fluctuations, which we also measure both experimentally and numerically. To assess the validity of the GT scaling prediction, in figure 5, we plot the diffusivity-free fluctuation-based Rayleigh number $Ra^{(\theta )}_*$ as a function of ${\mathcal {R}}$. Again, this representation leads to a very good collapse of the data, indicating diffusivity-free (or ‘ultimate’) dynamics. There is arguably more scatter in the experimental data than in the numerical data, because the latter benefit from an area horizontal average together with the time average when computing the r.m.s. temperature fluctuations (see § 2.3). Still, the theoretical prediction ${\mathcal {R}}^{3/5}$ is validated over four decades in ${\mathcal {R}}$, with a prefactor that is greater by approximately 20 % for the experimental data than for the numerical data. Somewhat surprisingly, the range of validity of the GT prediction seems to extend beyond ${\mathcal {R}} = 3\times 10^{-7}$, that is, the prediction for $Ra^{(\theta )}_*$ is validated over a more extended range of ${\mathcal {R}}$ than the prediction for $Nu_*$.

Figure 5. Diffusivity-free fluctuation-based Rayleigh number $Ra^{(\theta )}_*$ as a function of ${\mathcal {R}}$. Same symbols as in figure 2. The dotted line indicates the GT scaling prediction ${\mathcal {R}}^{3/5}$, while the dashed line indicates the scaling prediction ${\mathcal {R}}^{2/3}$ of the non-rotating ‘ultimate’ scaling regime. Inset: $Ra^{(\theta )}_*/{\mathcal {R}}^{3/5}$ versus ${\mathcal {R}}$.

4.3. Convective flow speed

As compared with the laboratory experiment, DNS allows us to readily characterize the velocity field. We thus turn to the r.m.s. velocity in the fluid domain inferred from the numerical simulations. In figure 6, we plot the diffusivity-free Reynolds number $Re_*$ as a function of ${\mathcal {R}}$. Once again, we obtain an excellent collapse of the numerical data onto a single master curve, which indicates diffusivity-free dynamics. The master curve is in excellent agreement with the GT prediction ${\mathcal {R}}^{2/5}$, as shown by the eye-guide in the main figure and by the compensated plot in the inset. Also shown in the figure is the prediction ${\mathcal {R}}^{1/3}$ associated with the non-rotating ultimate regime of thermal convection (Bouillaut et al. 2022), which is incompatible with the present rapidly rotating data points. We also note that the GT scaling prediction for $Re_*$ appears to be valid over a more extended range of ${\mathcal {R}}$ than the prediction for $Nu_*$. The reason may be that the exponents $2/5$ and $1/3$ of the rotating and non-rotating scaling predictions are rather close, possibly inducing a very smooth and extended cross-over region between the two power-law behaviours.

Figure 6. Diffusivity-free Reynolds number $Re_*$ extracted from the DNS as a function of ${\mathcal {R}}$. Same symbols as in figure 2. The dotted line indicates the GT scaling prediction ${\mathcal {R}}^{2/5}$, while the dashed line indicates the scaling prediction ${\mathcal {R}}^{1/3}$ of the non-rotating ‘ultimate’ scaling regime. Inset: $Re_*/{\mathcal {R}}^{2/5}$ versus ${\mathcal {R}}$.

4.4. Horizontal length scale

To further characterize the convective flow, we now turn to its characteristic horizontal length scale, extracted from the DNS based on the definition in (2.9). In figure 7, we plot $\ell _*^\perp$ as a function of ${\mathcal {R}}$ for the numerical dataset. Once again, this representation leads to a good collapse of the data, although the scatter appears greater than for the previous quantities of interest. This is likely a consequence of the shallower GT scaling prediction ${\mathcal {R}}^{1/5}$ for this quantity of interest. The eyeguide in the main figure indicates good agreement with this prediction, which is further confirmed by the compensated plot in the inset.

Figure 7. Dimensionless horizontal scale $\ell _*^\perp$ as a function of ${\mathcal {R}}$. The dotted line indicates the GT scaling prediction. Inset: $\ell _*^\perp /{\mathcal {R}}^{1/5}$ versus ${\mathcal {R}}$.

Comments are in order regarding the definition (2.9) of the horizontal length scale of the flow. Indeed, previous studies have reported different scaling behaviours using different proxies for the horizontal scale of the flow, only a subset of which agree with the GT prediction (Oliver et al. 2023; de Vries, Barker & Hollerbach 2023). One reason is that many of these proxies directly involve the dissipative scales and thus the viscosity. For instance, one should refrain from using the Taylor scale or the Kolmogorov scale of the flow when investigating GT scaling, as these scales are (partly) controlled by the viscous scale. Similarly, one should refrain from using the r.m.s. vorticity of the flow to diagnose the horizontal flow scale, as the r.m.s. vorticity is typically controlled by the small viscous scales. In other words, one expects to observe the zeroth law of turbulence (and therefore the GT predictions) for large-scale quantities, such as the horizontal integral scale of the flow defined in (2.9). As discussed in the following section, the no-slip bottom boundary condition of the present set-up prevents the emergence of domain-scale vortices, such that $\ell ^\perp$ in (2.9) is always smaller than the horizontal extent of the domain. This seems to be a simplification as compared with studies employing stress-free boundary conditions.

5. Discussion

Based on a hybrid experimental–numerical study, we have validated the scaling predictions of the GT regime for rapidly rotating radiatively driven convection. The scaling predictions for the heat transport and temperature fluctuations are validated using both experimental and numerical data. Following such cross-validation, the scaling predictions for the flow speed and horizontal scale are validated through further processing of the numerical data.

The GT scaling predictions for the temperature fluctuations, the flow speed and the horizontal scale of the flow are observed over many decades in diffusivity-free flux-based Rayleigh number ${\mathcal {R}}$, and they seem to arise even before the Nusselt number displays GT scaling. The same observation seems to hold for the rotating Rayleigh–Bénard system, where numerical studies indicate that the observation of the GT prediction for the heat transport requires more extreme parameter values than the GT predictions for the flow speed and flow structures (Julien et al. 2012; Song et al. 2024). At the experimental level, Vogt, Horn & Aurnou (2021) report velocity measurements that are compatible with the GT scaling prediction, even though the GT scaling regime for the heat transport remains elusive in rotating Rayleigh–Bénard experiments. Similarly, Guervilly, Cardin & Schaeffer (2019) report diffusivity-free behaviour for the Reynolds number and the azimuthal scale of the flow in a numerical study of rotating convection in a spherical shell, while the Nusselt number has still to reach the GT regime. These observations confirm that GT scaling for the Nusselt number arises further into the asymptotic regime of control parameters than GT scaling for the flow scale and flow speed. Although the diffusivity-free flux-based Rayleigh number ${\mathcal {R}}$ is a natural control parameter of our experiment, regime transitions in the Rayleigh–Bénard set-up are more often characterized in terms of the convective Rossby number:

(5.1)\begin{equation} Ro = \sqrt{\frac{\alpha g \Delta T}{4 H \varOmega^2}} = \sqrt{\frac{\mathcal{R}}{Nu_*}} , \end{equation}

the rapidly rotating regime corresponding typically to $Ro \lesssim 0.1$. We provide the values of $Ro$ for each panel of figure 3 in the figure caption to illustrate the transition to rapidly rotating convection. The GT asymptote in figure 4 corresponds to $Nu_*\simeq 2.8 \, {\mathcal {R}}^{3/5}$, which after substitution in (5.1) yields $Ro=0.6 \, {\mathcal {R}}^{1/5}$. As discussed above, clear GT scaling for the Nusselt number is observed for ${\mathcal {R}} \lesssim 3\times 10^{-7}$. The corresponding transition convective Rossby number is $Ro=0.6 \, (3\times 10^{-7})^{1/5} = 0.03$, within the expected range.

As mentioned in the introduction, the idealized numerical study that is closest to the present set-up is that of Barker et al. (2014). The authors employ a combination of local internal heating and cooling localized near the boundaries, together with an imposed background temperature gradient that is relaxed to a target value, with the goal of ensuring no average flux through the fixed-temperature top and bottom boundaries (the background temperature gradient is fixed in some simulations, while it obeys a relaxation equation when fluctuations are too strong). With this approach, the authors successfully validate the scaling predictions of Stevenson (1979), provided the heat transport is characterized using the temperature gradient at mid-depth and the Reynolds number is based on the vertical velocity component. Our set-up can be seen as an experimental realization of such internally heated and cooled convection, further confirming the numerical results of Barker et al. (2014) while avoiding the subtleties of their set-up and diagnostics. For instance, the boundaries of the laboratory experiment are thermally insulated to a good approximation, without the need for the adjusted background temperature gradient of Barker et al. (2014). Similarly, our DNS employs Chebyshev expansions to readily account for such insulated boundaries (contrasting with the Fourier decomposition of Barker et al. (2014)). Additionally, we follow the standard practice of defining the Nusselt number based on a temperature drop, showing consistent behaviour with the temperature gradient considered by Barker et al. (2014). Finally, we believe that the solid boundary at the bottom of the experimental tank – and associated no-slip boundary condition in the DNS – also contributes to ‘simplifying’ the set-up, in the sense that it prevents (or delays) any accumulation of horizontal kinetic energy into large-scale vortices (Favier, Silvers & Proctor 2014; Guervilly, Hughes & Jones 2014; Stellmach et al. 2014; Kunnen et al. 2016; Favier et al. 2019; Aguirre Guzmán et al. 2020; Maffei et al. 2021). (Indeed, running a few numerical test cases with a free-slip bottom boundary condition did lead to large-scale-vortex formation sufficiently far from threshold, while large-scale vortices were not observed for the realistic no-slip bottom boundary condition. When they develop, these large scale vortices modify the scaling for the horizontal scale and the horizontal speed of the flow.) Overall, our study shows excellent agreement with the pioneering study of Barker et al. (2014) while offering a ‘realistic’ set-up, in the sense that it can be realized in the laboratory.

To conclude, radiatively driven convection offers a robust way of achieving the diffusivity-free (or ultimate) regimes of thermal convection in the laboratory, with and without global rotation. Additional physical ingredients could be included in the experimental and numerical set-ups, with the goal of investigating the diffusivity-free regimes of a broader range of natural convective flows. Magnetic field crucially influences convection in planetary and stellar interiors and could be the natural next step.