1. Introduction
Convection plays a decisive role in transport for geophysical and astrophysical environments. As such, a fundamental understanding of convection, as affected by rotation, will improve our knowledge of planetary and stellar dynamics. However, in comparison with boundary heating, convection driven by internal sources remains less studied. In planetary interiors, internal heating primarily due to the radioactive decay of isotopes is crucial to the evolution of plate tectonics and the generation of a dynamo (Schubert, Turcotte & Olson Reference Schubert, Turcotte and Olson2001; Ricard Reference Ricard2015). In the cores of planets where dynamos typically occur, before inner core nucleation, internal heating by secular cooling drives the dynamo (Roberts Reference Roberts2015; Dormy Reference Dormy2025). Crucially, for planetary and astrophysical dynamos, the convection that creates large-scale magnetic fields is constrained by rotation (Tilgner Reference Tilgner2015), so an insight into internal heating with rotation is required.
Before considering rotation, the dynamics of internally heated convection (IHC) on its own warrants an introduction. In contrast to Rayleigh–Bénard convection (RBC), where fluid in a periodic plane layer is heated from below and cooled from above, with internal heating, the choice of heating and boundary conditions results in different flows (Goluskin Reference Goluskin2016). For example, if net internal heating is zero and the boundaries are thermally insulating, a purely internally heated flow is generated devoid of the effects of boundary layers (Bouillaut et al. Reference Bouillaut, Lepot, Aumaître and Gallet2019; Miquel et al. Reference Miquel, Lepot, Bouillaut and Gallet2019; Bouillaut et al. Reference Bouillaut, Flesselles, Miquel, Aumaître and Gallet2022). This set-up allows for diffusivity-free heat transport (Lepot, Aumaître & Gallet Reference Lepot, Aumaître and Gallet2018). However, if the net heating is unity, two distinct flows occur based on the choice of boundary condition at the bottom. If the top is isothermal and the bottom an insulator (denoted as IH3; Goluskin Reference Goluskin2016), all internal heating leaves via the upper boundary and the flow shares similarities with RBC. A Nusselt number, quantifying the ratio of the averaged total vertical heat transport to conduction, can be defined. However, if both the top and bottom are isothermal (denoted IH1 in previous literature), heat can pass through both boundaries. In the latter scenario, convection quantifies the asymmetry in the heat flux out of the boundaries. As the mean conductive heat transport is zero, an exact Nusselt number does not exist, and the flow described contains notable differences from RBC.
A fundamental task in the study of fluid dynamics involves characterising the statistical properties of the flow as a function of the control parameters of the system (Doering Reference Doering2019; Lohse & Shishkina Reference Lohse and Shishkina2023). For IHC, this process is less thoroughly studied, especially for low Prandtl numbers, due to computational and experimental difficulties (see Goluskin Reference Goluskin2016 for a review of experimental results). The non-dimensional control parameters are the Rayleigh number
$R$
, quantifying the ratio of destabilisation by internal heating to diffusion, the Prandtl number
$\textit{Pr}$
, the ratio of viscous to thermal diffusivity, and the aspect ratio
$\varGamma$
, comparing the horizontal with the vertical scale of the system. Two properties of interest in IHC that quantify the effect of convective heat transport when the flow is laminar, chaotic and fully turbulent, are
$\overline {\langle wT \rangle }$
and
$\overline {\langle T \rangle }$
, the mean vertical heat transport by convection, and the mean temperature. For notation in this paper,
$\langle \boldsymbol{\cdot }\rangle$
denotes horizontal and time averaging, and
$\overline {\boldsymbol{\cdot }}$
denotes vertical averaging. Unlike RBC,
$\overline {\langle wT \rangle }$
and
$\overline {\langle T \rangle }$
in IHC cannot be related a priori and describe two separate emergent quantities of the system. Furthermore, the differences between IHC and RBC are observable by comparison of the horizontally averaged temperature and root mean square (r.m.s.) velocity distributions. The IHC temperature profiles show the existence of a stably stratified boundary layer at the bottom with a thickness larger than the unstable thermal boundary layer at the top (Goluskin & Spiegel Reference Goluskin and Spiegel2012; Goluskin & Van der Poel Reference Goluskin and Van der Poel2016). The r.m.s. velocity highlights the plumes from the upper boundary and asymmetric convective winds (relative to the midplane) within the domain.
The effect of rotation alters turbulent convection in several ways. Irrespective of buoyancy, the Coriolis force creates an asymmetry in velocity gradients and introduces an Ekman boundary layer (Greenspan Reference Greenspan1968; Pedlosky Reference Pedlosky2013). The non-dimensional Ekman number
$E$
, the ratio of viscous diffusion to rotation, gauges the effect of rotation and the creation of anisotropy in the flow. In the limit of zero Ekman number, the flow is subject to the Taylor–Proudman constraint, where the velocity is constant parallel to the axis of rotation and the dynamics is two-dimensional (Gallet Reference Gallet2015).
Rotating convection has been thoroughly studied in the plane layer when the heating source is the boundaries (Zhong, Ecke & Steinberg Reference Zhong, Ecke and Steinberg1993; Julien et al. Reference Julien, Legg, McWilliams and Werne1996; Knobloch Reference Knobloch1998; Vorobieff & Ecke Reference Vorobieff and Ecke2002; Stevens, Clercx & Lohse Reference Stevens, Clercx and Lohse2013; Aguirre Guzmán et al. Reference Aguirre Guzmán, Madonia, Cheng, Ostilla-Mónico, Clercx and Kunnen2021; Song, Shishkina & Zhu Reference Song, Shishkina and Zhu2024), notably also in the asymptotic limit of rapid rotation (Julien et al. Reference Julien, Knobloch and Werne1998, Reference Julien, Knobloch, Rubio and Vasil2012a
,
Reference Julien, Rubio, Grooms and Knoblochb
, Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016; Sprague et al. Reference Sprague, Julien, Knobloch and Werne2006). The first studies on the onset of convection, summarised in Chandrasekhar (Reference Chandrasekhar1961), highlighted the exact role of rotation in inhibiting the onset of convection and in the dominant length scales of the flow. Specifically, the critical Rayleigh number for linear instability scales as
$E^{-4/3}$
, while the dominant non-dimensional length scales of the system go as
$E^{1/3}$
. These scalings were analytically demonstrated for internal heating in Arslan (Reference Arslan2024). Studies of rotating RBC indicate that the dominant length scale remains significant as the Rayleigh and Ekman numbers are varied such that the flow becomes turbulent. Furthermore, numerical experiments have shown that the
$E {-} R$
parameter space is rich with a variety of flow states (see Kunnen Reference Kunnen2021 and Ecke & Shishkina Reference Ecke and Shishkina2023 for reviews). For sufficiently small
$E$
and large enough
$R$
for convection, the flow is geostrophically turbulent with characteristics different to buoyancy-driven turbulence (Song et al. Reference Song, Shishkina and Zhu2024). Recent investigations demonstrate geostrophic turbulence with net zero heating in an insulated domain, confirming that the regime can be achieved with internal heating, albeit one that differs from that in this work, by studying the scaling laws for heat transport in the flow (Bouillaut et al. Reference Bouillaut, Miquel, Julien, Aumaître and Gallet2021; Hadjerci et al. Reference Hadjerci, Bouillaut, Miquel and Gallet2024).
It is worth highlighting that rotating IHC in alternative geometries is also relevant to geophysics. For example, for an internally heated fluid in a sphere, linear stability analysis and simulations have documented the flow state for a small range of Rayleigh numbers above the onset of convection (Jones, Soward & Mussa Reference Jones, Soward and Mussa2000). An important discovery has been that, in spherical geometries, subcritical convection occurs for low Prandtl numbers (
$\textit{Pr} \leqslant 0.1$
) (Guervilly & Cardin Reference Guervilly and Cardin2016; Kaplan et al. Reference Kaplan, Schaeffer, Vidal and Cardin2017). In addition to the sphere, studies of rotating spherical shells combine internal and boundary heating to drive convection (Gastine, Wicht & Aubert Reference Gastine, Wicht and Aubert2016). Both geometries model the core of the Earth at different periods in its history, albeit without the effects of magnetism. Returning to the plane layer, the energy method also predicts subcritical convection for non-rotating IHC (Goluskin Reference Goluskin2016). However, the effect of rotation on nonlinear stability is unknown as rotating IHC with isothermal boundaries remains largely unexplored.
This work is motivated by the lack of studies of rotating IHC in the plane layer and, in part, by the analytical results obtained in Arslan (Reference Arslan2024). In Arslan (Reference Arslan2024), the author established the onset of convection for large
$\textit{Pr}$
(
$\textit{Pr} \gt 1$
) and rotation (small
$E$
) and proved lower bounds on the heat flux out of the bottom boundary
$\mathcal{F}_b$
and the mean temperature
$\overline {\langle T \rangle }$
in different areas of the
$E {-} R$
space. That work uses the auxiliary functional method (Goulart & Chernyshenko Reference Goulart and Chernyshenko2012; Chernyshenko Reference Chernyshenko2022) to construct a variational problem for
$\mathcal{F}_B$
and
$\overline {\langle T \rangle }$
, the solution of which gives a lower bound in terms of
$R$
and
$E$
. The method can prove sharp bounds on the long-time behaviour of ordinary and partial differential equations (Fantuzzi et al. Reference Fantuzzi, Goluskin, Huang and Chernyshenko2016; Tobasco, Goluskin & Doering Reference Tobasco, Goluskin and Doering2018; Rosa & Temam Reference Rosa and Temam2022). In the case of the Navier–Stokes equations, quadratic auxiliary functionals have been used to prove bounds on turbulent flows (see Fantuzzi, Arslan & Wynn Reference Fantuzzi, Arslan and Wynn2022 for a detailed review). Although the method has proven sharp bounds for specific flows, in other cases, the bounds appear not to capture the dynamics.
Bounds on
$\overline {\langle wT \rangle }$
, for IHC, are an example where the sharpness of bounds remain unknown, even though the proofs require delicate mathematical constructions and estimates (Arslan et al. Reference Arslan, Fantuzzi, Craske and Wynn2021a
,
Reference Arslan, Fantuzzi, Craske and Wynnb
; Kumar et al. Reference Kumar, Arslan, Fantuzzi, Craske and Wynn2022; Arslan & Rojas Reference Arslan and Rojas2025). The bounds on long-time averages indicate that IHC has distinct features from RBC. Recent work shows that, even when the internal heating is non-uniform, open questions remain, and that the limit of heating at the boundary remains unavailable (Arslan et al. Reference Arslan, Fantuzzi, Craske and Wynn2024). One means of testing the sharpness of bounds is via direct numerical simulations of a representative range of parameter values. In general, proving results on rotating convection is challenging with a standard application of bounding methods because the Coriolis force does not appear in energy identities. The bounds proven for RBC or IHC, either consider the infinite Prandtl limit (Constantin, Hallstrom & Putkaradze Reference Constantin, Hallstrom and Putkaradze1999; Doering & Constantin Reference Doering and Constantin2001; Yan Reference Yan2004; Arslan Reference Arslan2024) or a reduced set of equations where time and length are rescaled by factors of
$E^{1/3}$
(Grooms & Whitehead Reference Grooms and Whitehead2014; Pachev et al. Reference Pachev, Whitehead, Fantuzzi and Grooms2020). While we cannot span the entire range of parameters, this work provides the first results for the variations of
$\overline {\langle wT \rangle }$
and
$\overline {\langle T \rangle }$
with
$R$
and
$E$
by direct numerical simulation in the rotationally affected regime.
The paper is organised as follows. Section 2 describes the simulation method. In § 3.1, we analyse the onset of convection in rotating IHC, while § 3.2 describes the resulting flow morphology. Global quantities are analysed in § 3.3, while their relationship to upper bounds is discussed in § 3.4. The remaining results sections discuss the behaviour of local temperature and velocity statistics (§ 3.5), boundary layer behaviour (§ 3.6) and dissipation rates (§ 3.7). Finally, § 4 briefly summarises the results obtained and Appendix A summarises the simulation data.
2. Numerical method

Figure 1. A non-dimensional schematic diagram for rotating uniform IHC. The upper and lower plates are at the same temperature, and the domain is periodic in the
$x$
and
$y$
directions and rotates about the
$z$
axis. Here,
$\mathcal{F}_B$
and
$\mathcal{F}_T$
are the mean heat fluxes out the bottom and top plates,
$\overline {\langle T \rangle }$
the mean temperature and
$g$
is the acceleration due to gravity.
We consider a fluid layer confined by two horizontal walls separated a distance
$d$
, and periodic in both horizontal directions with equal periodicity lengths
$\varGamma d$
. The fluid has kinematic viscosity
$\nu$
, thermal diffusivity
$\kappa$
, density
$\rho$
, specific heat capacity
$c_p$
and thermal expansion coefficient
$\alpha$
. Gravity acts downwards with strength
$g$
, while the fluid rotates at rate
$\varOmega$
about the vertical axis and is uniformly heated internally at a volumetric rate
$H$
. Density variations are modelled with the Boussinesq approximation where they only contribute to buoyancy forces. Figure 1 shows a schematic of the problem.
To non-dimensionalise the system, we use
$d$
as the characteristic length scale,
$d^2/\kappa$
as the time scale and
$d^2H/(\kappa \rho c_p)$
as the temperature scale (Roberts Reference Roberts1967). The velocity of the fluid
$\boldsymbol {u}(\boldsymbol {x},t)=u(\boldsymbol {x},t)\boldsymbol {e}_1+v(\boldsymbol {x},t)\boldsymbol {e}_2+w(\boldsymbol {x},t)\boldsymbol {e}_3$
, pressure
$p((x),t)$
, and temperature
$T(\boldsymbol {x},t)$
then satisfy the following non-dimensional equations:



subject to no-slip (
$\boldsymbol {u}=0$
) and isothermal (
$T=0$
) boundary conditions at both walls.
The system’s response is governed by three non-dimensional parameters: the Rayleigh, Ekman and Prandtl numbers, defined as

The convective Rossby number, a parameter often used to directly compare the relative importance of buoyancy and rotation is obtained as
$Ro = E\sqrt {R/\textit{Pr}}$
.
This choice for non-dimensionalisation follows convention (Goluskin & Van der Poel Reference Goluskin and Van der Poel2016; Kazemi, Ostilla-Mónico & Goluskin Reference Kazemi, Ostilla-Mónico and Goluskin2022; Arslan Reference Arslan2024), ensuring consistency with prior work, but it is not the only reasonable choice. Unlike RBC, where the conventional non-dimensionalisation results in velocity and temperature scales are
$\mathcal{O}(1)$
for
$\textit{Pr}=1$
, IHC exhibits velocity and temperature scales that vary with
$R$
. Our choice of non-dimensionalisation results in increasing velocity and decreasing temperature scales as
$R$
increase, as discussed later. We note that this is not the only possible choice of non-dimensionalising variables: an alternative choice is to use a free-fall-based velocity scale similar to Rayleigh–Bénard, with
$d^2 H/\kappa$
as a temperature scale and
$\sqrt {\alpha gd^3H/\kappa }$
as a velocity scale. This would instead result in decreasing velocity and temperature scales with
$R$
, but for
$\textit{Pr}\sim \mathcal{O}(1)$
, these would be of the same order of magnitude, which may benefit the numerical stability of the code. However, it would have the drawback of making comparisons with existing work less immediate.
We solve (2.1)–(2.3) using the open-source, second-order centred finite difference code AFiD (Van Der Poel et al. Reference Van Der Poel, Ostilla-Mónico, Donners and Verzicco2015), which has been extensively validated for IHC (Goluskin & Van der Poel Reference Goluskin and Van der Poel2016; Kazemi et al. Reference Kazemi, Ostilla-Mónico and Goluskin2022). The non-dimensional periodicity length
$\varGamma$
is selected such that it is at least the value in Goluskin & Van der Poel (Reference Goluskin and Van der Poel2016) which ensures that
$\overline {\langle T \rangle }$
and
$\overline {\langle wT \rangle }$
remain independent of domain size.
Following Goluskin & Van der Poel (Reference Goluskin and Van der Poel2016), the resolution adequacy is measured by the exact relationships


and ensuring that they do not diverge numerically by less than 1 %. This is a stricter resolution requirement than looking at the averaged Kolmogorov or Batchelor scales, and has been shown to be adequate for thermal convection (Stevens, Verzicco & Lohse Reference Stevens, Verzicco and Lohse2010). For the non-rotating cases, we find that the resolutions in Goluskin & Van der Poel (Reference Goluskin and Van der Poel2016) can be reduced without loss of accuracy. In the most adverse cases, the discretisation used in this manuscript has a maximum point distance of
$1.5$
Kolmogorov length scales. In the non-rotating case, the upper thermal boundary layer, whose extent can be estimated by
$\langle \overline {T}\rangle /(({1}/{2})+\overline {\langle wT \rangle } )$
, is resolved by at least six points. As rotation increases, the number of points in the
$z$
direction must be increased to resolve the increasingly thinner Ekman (velocity) boundary layers.
In this study, we will set
$\textit{Pr}=1$
, and vary
$R-E$
in the ranges
$R\in [3.16\times 10^5,10^{10}]$
and
$E\in [10^{-6},\infty )$
. The full parameter space explored is discussed in the next section, with a summary of results and resolutions provided in the Appendix. Except for the initial run, simulations are initialised from adjacent points in the
$R$
$E$
parameter space. After transients, statistics are collected for a time interval of length
$t_{av}$
which ensured that the time-averaged values of
$\overline {\langle T \rangle }$
and
$\overline {\langle wT \rangle }$
were less than 1 % different from those corresponding to half the run time. The values of
$t_{av}$
are available in the Appendix. In general
$t_{av}$
is dependent on
$R$
, as the flow becomes faster with increasing
$R$
, and hence the necessary averaging time decreases.
3. Results
3.1. Onset of convection
We start our exploration by probing the onset of convection. As in RBC, rotation stabilises the system against vertical motion, delaying the onset of convection. However, unlike RBC, where a precise relationship between the Rayleigh and Ekman numbers determines the onset under rotation (Chandrasekhar Reference Chandrasekhar1961), in IHC only an asymptotic scaling has been determined, with an analytic exact relation for all Ekman numbers not known, due to the mathematical nature of the governing equations (Arslan Reference Arslan2024). This scaling follows
$R_L\sim E^{-4/3}$
, where
$R_L$
is the critical Rayleigh number above which the flow is linearly unstable (Arslan Reference Arslan2024). The exponent
$-4/3$
matches that of rotating RBC, reflecting a similarity between the two systems.
Using numerical simulations, we mapped the
$\textit{RE}$
parameter space to estimate the prefactor in this relationship. Figure 2 presents the simulations conducted, displayed both in
$\textit{RE}$
space and in
$E\!{\tilde {R}}$
space, where
$\tilde {R}=\textit{RE}^{4/3}$
measures the degree of supercriticality in rotating IHC. The panels also show
$Ro=1$
, marking the region where rotation significantly influences the dynamics. When rotation is significant, independent of
$R$
we find convection for all cases with
$\tilde {R}\geqslant 68.3$
while the flow remains quiescent when
$\tilde {R}\leqslant 46.8$
. This means that the critical
$\tilde {R}$
for the onset of convection when the system is rotating must lie somewhere between those values.
When assessing system stability, sufficient simulation time is crucial for the flow to adjust to parameter changes. Increasing rotation initially suppresses turbulence, leading to a gradual rise in mean temperature, which destabilises the system. If given enough time (which may be rather long), turbulence may re-emerge. Thus, we define a purely conductive state as one where the mean temperature deviates by less than 1 % from its conductive value and turbulence decay persists.

Figure 2. Explored parameter space in
$\textit{RE}$
(a) and
$E\!\tilde {R}$
(b) variables. Simulations conducted are marked as green squares, except when they produce purely conductive states and are marked as black stars. On both graphs, the dashed black line marks
$Ro=1$
and the dash-dot line marks
$\tilde {R}=50$
.
3.2. Flow morphology
Figure 3 illustrates how the flow evolves with increasing rotation (decreasing
$E$
) at
$R=10^{10}$
. Across all cases, plumes or other convective structures originate primarily in the top boundary layer, which is unstably stratified. This mixes the fluid, and produces turbulence which then mixes the stably stratified lower region.
The flow structure changes significantly with rotation. Without rotation, convection is dominated by downwards plumes. As rotation increases, plumes elongate vertically and become more organised, eventually forming Taylor columns similar to those in rotating RBC. At
$E=10^{-6}$
, flow structures are highly constrained horizontally, and turbulence is nearly absent due to the low supercriticality (
$\tilde {R}=100$
). This means that our simulations are mainly in the parameter space corresponding to intermediate, or ‘rotation-affected’ convection. While the smallest
$E$
values considered correspond to
$Ro=0.1$
, the flow remains insufficiently supercritical to enter the geostrophic turbulence regime, instead exhibiting quasi-laminar structures.

Figure 3. Volumetric visualisation of the instantaneous temperature field for
$R=10^{10}$
and from a–f no rotation,
$E=10^{-4}$
(
$Ro=10$
),
$E=3.16\times 10^{-4}$
(
$Ro=3.16$
),
$E=10^{-5}$
(
$Ro=1$
),
$E=3.16\times 10^{-5}$
(
$Ro=0.31$
) and
$E=10^{-6}$
(
$Ro=0.1$
). The flow can be seen to change structures from plume dominated to column dominated.
3.3. Global quantities
In this section we analyse the system response throughout the parameter space by examining three global parameters:
$\mathcal{F}_B$
,
$\overline {\langle T \rangle }$
and the wind Reynolds number
${\textit{Re}}_w=U_wd/\nu$
, where
$U_w$
is a characteristic velocity for the wind defined as

The left panels of figure 4 shows the behaviour of
$\mathcal{F}_B$
, i.e. the fraction of heat transported through the bottom plate, through the
$R$
$E$
parameter space. This is the quantity conventionally used when representing mathematical upper bounds. In our set up of IHC, all heat must leave through the boundaries. Therefore, the fractions corresponding to heat leaving from the top and bottom boundaries must be equal to unity. Mean vertical convective heat flux,
$\overline {\langle wT \rangle }$
, causes more heat to leave through the top boundary than the bottom. The exact relationships which link the variables are

Therefore, representing
$\overline {\langle wT \rangle }$
is equivalent to representing
$\mathcal{F}_B$
.
Results for non-rotating IHC (
$E=\infty$
) from Goluskin & Van der Poel (Reference Goluskin and Van der Poel2016) are included for comparison, and an excellent agreement between both sets of simulations can be appreciated, despite the lower resolutions used in this manuscript. Two dependencies are apparent from the graph: first, that as
$R$
increases, the fraction of heat transported through the bottom plate
$\mathcal{F}_B$
decreases regardless of the value of
$E$
. While this decrease is slow, going from
$40\,\%$
at the lowest
$R$
to the lowest value seen in the graph of
$15\,\%$
for rotating IHC, the monotonic trend of
$\mathcal{F}_B$
is apparent.
While it may seem intuitive that increased driving leads to higher convective vertical transport and more heat leaving from the top, this is not the case for two-dimensional (2-D) IHC (Goluskin & Spiegel Reference Goluskin and Spiegel2012). Goluskin & Van der Poel (Reference Goluskin and Van der Poel2016) postulated that this was due to the existence of a competing effect: increased turbulence leads to shear-driven mixing of the bottom boundary layer, which allows for more heat to escape. In 2-D simulations, the geometrical constraints on the flow result in organisation. A large-scale circulation develops, which heavily increases the shear-driven mixing. In three dimensions, this circulation does not form for the non-rotating case. In our simulations, we do not observe non-monotonic behaviour of
$\mathcal{F}_B(R)$
, despite rotation changing the flow morphology, further confirming that this non-monotonicity is a product of the 2-D simulations.
The second behaviour observed is that there are regions of parameter space where rotation results in lower values of
$\mathcal{F}_B$
. This is especially pronounced for the larger values of
$R$
seen in the bottom left panel of figure 4. This means that the overall level of vertical convective transport increases. Furthermore, the value of
$E$
where
$\mathcal{F}_B$
is minimum can be observed to depend on
$R$
. We will return to this phenomena, which is akin to increases in the heat transport seen for rotating RB at comparable values of the Rayleigh number.

Figure 4. Global responses. (a,b,c)
$\mathcal{F}_B$
(a),
$\overline {\langle T \rangle }$
, (b) and
${\textit{Re}}_w$
(c) against
$R$
for several values of
$E$
. (d,e, f) same quantities against
$E$
for several values of
$R$
. Theoretical bounds are not shown as they lie significantly outside the data range.
The middle panels of figure 4 show the behaviour of the volumetrically averaged temperature
$\overline {\langle T \rangle }$
in the parameter space. Both patterns of behaviour seen in
$\mathcal{F}_B$
are also present here: the mean temperature decreases with
$R$
, and certain finite values of
$E$
lead to decreases of the mean temperature when compared with those for the non-rotating case. We also note that the mean temperature scales roughly like
$R^{-0.2}$
. This scaling is only an approximation and disappears when the values are examined closely. However, what is clear is that the temperature scale of the system diminishes as
$R$
increases, due to the higher level of turbulent convection.
Finally, the right panels of figure 4 show the behaviour of the wind Reynolds number. For this quantity, the dependence on driving parameters is simpler:
$R$
increases velocity fluctuations, while decreasing
$E$
damps these. No ‘optimal’ rotation rate which increases
${\textit{Re}}_w$
is observed, demonstrating that the variations of
$\mathcal{F}_B$
and
$\overline {\langle T \rangle }$
are due to flow organisation and not to increased wind. We also note that
${\textit{Re}}_w$
very approximately scales as
${\textit{Re}}_w\sim R^{0.4}$
. While the exact value of this exponent is not important, what is clear is that the scaling is less steep than one which would correspond to a free-fall velocity, i.e.
$R^{0.5}$
. This is because the velocity scale is a response of the system, and because the temperature scale is decreasing with
$R$
(as seen from the behaviour of
$\overline {\langle T \rangle }$
), the forcing and hence the scaling exponent of
${\textit{Re}}_w$
is reduced as a consequence. We note that our results are different from the scalings observed for
${\textit{Re}}_w$
in 2-D simulations by Wang, Lohse & Shishkina (Reference Wang, Lohse and Shishkina2020), which see an approximate free-fall scaling. This can again be attributed to the large differences in velocity statistics observed in IHC when comparing 2-D and 3-D simulations (Goluskin & Van der Poel Reference Goluskin and Van der Poel2016).

Figure 5. Optimal transport: normalised
$\langle wT \rangle$
(a) and
$\overline {\langle T \rangle }$
(b) against
$1/Ro$
for several values of
$R$
.
To further analyse the phenomena of optimal transport, in figure 5 we show the behaviour of
$\overline {\langle wT \rangle }$
and
$1/{\overline {\langle T \rangle } }$
normalised by their values for the non-rotating case (
$E=\infty$
) as a function of the inverse Rossby number. This representation allows for a better collapse in the independent variable, while displaying the relative increases more clearly. Three regimes can be observed: (i) for small values of
$1/Ro$
, rotation does not affect the flow, (ii) for intermediate values of
$1/Ro$
, rotating the flow organises the motions and increases convective transport resulting in larger values of
$\overline {\langle wT \rangle }$
and smaller values of
$\overline {\langle T \rangle }$
, (iii) for large values of
$1/Ro$
, rotation suppresses convection resulting in smaller
$\overline {\langle wT \rangle }$
and larger
$\overline {\langle T \rangle }$
. We note that the collapses are not perfect: the changing velocity scale is reflected in slight shifts in the curves.
We also note that the relative increases in
$\overline {\langle wT \rangle }$
are much larger than those of
$1/{\overline {\langle T \rangle } }$
, reaching up to
$25\,\%$
in the former, while remaining below
$10\,\%$
for the latter. Furthermore, the increases in
$1/{\overline {\langle T \rangle } }$
are very
$R$
dependent, and are almost eliminated for the highest value of
$R$
. This is in line with what is seen for the Nusselt number in rotating RBC: a ‘rotation-affected’ regime, where the Ekman boundary layers formed at the no-slip plates transfer heat effectively and cause small increases in the heat transfer that vanish at higher thermal drivings due to changes in flow organisation (Kunnen Reference Kunnen2021). This change in flow organisation can be also appreciated in figure 3. Furthermore, for
$R=10^{10}$
in our simulations, the heat transfer increase has practically disappeared from the right panel of figure 5 in line with what is expected from rotating RBC.
However, for
$\overline {\langle wT \rangle }$
the increases remain for all
$R$
. As mentioned earlier, the asymmetry between heat escaping the top and bottom boundary layers is also driven by the shear mixing of the stably stratified bottom. The convective structures present in rotating turbulence result in smaller
${\textit{Re}}_w$
, and mix the bottom less well. This second mechanism is probably behind the persistent increases in
$\overline {\langle wT \rangle }$
: mean vertical convective heat flux compensates the decreased amount of heat escaping the bottom plate, and increases the asymmetry between boundary layers.
3.4. Comparison with bounds
The global quantities can also be compared with the rigorous bounds proven with the methods mentioned in the introduction. We therefore include a brief comparison with the results proven in Arslan (Reference Arslan2024). An exact evaluation of the scaling laws of
$\overline {\langle wT \rangle }$
and
$\overline {\langle T \rangle }$
was not made from the results, nevertheless, a qualitative comparison with the rigorous bounds are possible. However, we keep in mind that the bounds in Arslan (Reference Arslan2024) are proven for (2.2) in the limit of infinite
$\textit{Pr}$
, while the numerical investigation was carried out at
$\textit{Pr}=1$
. An extension of infinite
$\textit{Pr}$
bounds for rotating convection to finite
$\textit{Pr}$
was made in Tilgner (Reference Tilgner2022). The author shows that, for fixed
$R$
and
$E$
in rotating RBC, the bounds at finite
$\textit{Pr}$
carry the correction
$(\sqrt {1 + A } + A)^2$
, where
$A = (0.088 + 0.077 c_0 Ra^{1/2})\textit{Pr}^{-1}Ra^{1/2}\,E^{1/3}$
and
$c_0={\rm max}(0,L/2)$
depends on the domain length
$L$
in the horizontal directions. When
$\textit{Pr}\rightarrow \infty$
the correction term tends to 1 and the bound on heat transport is given by that obtained in the infinite
$\textit{Pr}$
limit. A similar, partial, result for arbitrary
$\textit{Pr}$
in IHC is unknown and would be required for a direct assessment of the sharpness of the bounds to the direct numerical simulations.
The range of parameter space explored in this work, as shown in figure 2, falls within regimes III and IIIa of Arslan (Reference Arslan2024) (see table 1 and figure 8). Both regimes correspond to rotationally affected convection, where III is the area in
$R{-}E$
space where
$ R^{-2} \lesssim E \lesssim 5$
while is where to
$R^{-2/3} \lesssim E \lesssim R^{-5/9}$
. For
$\mathcal{F}_B$
, the bounds have been proven for regime III, the rotationally affected regime, and are that
$\mathcal{F}_B \gtrsim R^{-2/3}E^{2/3} + R^{-1/2}E^{1/2} \ln {|1 - R^{-1/3}E^{1/3}|}$
. Comparing with the first row of figure 4, for fixed
$E$
,
$\mathcal{F}_B$
decreases with
$R$
, as in the bound. For fixed
$R$
the lower bound matches the behaviour only at higher
$E$
(slower rotation). On the other hand, the bounds for
$\overline {\langle T \rangle }$
, in the region where we have carried out simulations, are split into regimes III where
${\overline {\langle T \rangle } } \gtrsim R^{-2/7}E^{2/7}$
and IIIa where
${\overline {\langle T \rangle } } \gtrsim R^{-1}E^{-1}$
. Regime IIIa captures the increased influence of rotation on convection. The middle column of figure 4 matches the behaviour of the bounds. In particular, fixing
$R$
there are, at the very least, two distinct behaviours of
$\overline {\langle T \rangle }$
as
$E$
varies, and this correlates with the bounds in III and IIIa. Therefore, the bounds on
$\overline {\langle T \rangle }$
qualitatively match the simulation results, while those for
$\mathcal{F}_B$
can be improved for smaller
$E$
.
3.5. Temperature and velocity statistics
To further understand the results obtained, we turn to an examination of the local behaviour of the fluid quantities. Starting with temperature, figure 6 shows the mean
$\langle T \rangle$
and r.m.s.
$\langle T^\prime \rangle = \sqrt {\langle T^2\rangle - \langle T\rangle ^2 }$
temperature profiles for several values of
$R$
and
$E$
. In the top panels, both mean temperature and fluctuations generally decrease with increasing
$R$
, consistent with the decline in
$\overline {\langle T \rangle }$
. Without rotation, the mean profiles exhibit a well-mixed bulk with a nearly constant temperature, a sharp gradient in the top boundary layer and a more gradual gradient in the stably stratified bottom boundary layer. As
$R$
increases, the well-mixed region expands, and the temperature profile flattens.
In contrast, the fluctuations show two peaks: a more pronounced one near the top boundary layer and a weaker one near the bottom. The minimum fluctuation region is located between
$z=0.2$
and
$z=0.4$
, and shifts downward as
$R$
increases. We identify this minimum with the interface between stably and unstably stratified layers. This is a better diagnosis than the position of the mean temperature maximum, as the latter is located closer to the top plate than the bottom plate due to the temperature gradient in the bulk, and does not correspond to the physical reality of the system.
With rotation, the shape of the mean temperature profiles remain largely unchanged except for the emergence of a bulk temperature gradient, which becomes steeper as
$E$
decreases. This gradient, opposite in sign to the weak positive gradient seen without rotation, resembles a phenomenon observed in rotating RBC (Kunnen Reference Kunnen2021) and is quantified in figure 7. Initially, the bulk remains well mixed with a slight positive gradient at low
$R$
, but as rotation strengthens, the gradient reverses. Once
$1/Ro\gt 1$
, the gradient grows rapidly in absolute terms.
Temperature fluctuations also respond to rotation: fluctuations in the bottom boundary layer decrease, while those in the top increase. This suggests that rotation reduces the mixing in the bottom region, and supports the idea mentioned earlier that increases in
$\overline {\langle wT \rangle }$
are driven by less thorough mixing of the stably stratified layer.

Figure 6. Mean
$\langle T \rangle$
(a,c) and fluctuation
$\langle T^\prime \rangle$
(b,d) profiles for several values of
$R$
in the non-rotating case (a,b) and for several values of
$E$
at
$R=10^{10}$
(c,d).

Figure 7. Temperature gradient in the mid-gap as a function of
$E$
(a) and
$1/Ro$
(b).

Figure 8. Velocity fluctuations in the vertical (a,c) and horizontal (b,d) directions for several values of
$R$
in the non-rotating case (a,b) and for several values of
$E$
at
$R=10^{10}$
(c,d).
We now examine the velocity statistics. To quantify the fluctuations in the horizontal direction we define the horizontal velocity
$u_h=\sqrt {u_x^2+u_y^2}$
. Figure 8 presents both the horizontal
$\langle u_h^\prime \rangle$
and vertical
$\langle u_z^\prime \rangle$
velocity fluctuations. The vertical velocity follows a straightforward trend aligning with expectations: fluctuations (and transport) increase with
$z$
before decreasing near the velocity boundary layer at the top plate. When rotation is introduced, vertical velocity is increasingly suppressed as
$E$
decreases. However, the overall shape of the profile remains relatively unchanged, aside from the fact that at the smallest
$E$
, the stably stratified layer becomes increasingly quiet.
In contrast, horizontal motions exhibit more complex behaviour. In the absence of rotation, the upper half of the system behaves similarly to RBC: fluctuations peak near the wall, marking the edge of the boundary layer, while the bulk region maintains an almost uniform fluctuation level. Introducing rotation preserves this overall structure but thins the boundary layer and enhances fluctuations as it transitions into an Ekman layer.
The lower half of the domain features two distinct regions in
$\langle u_h^\prime \rangle$
fluctuations, with a local maximum (for small
$E$
) or an inflection point (for large
$E$
) close to the wall, and a second local maximum a bit further away. The peak closest to the wall corresponds to the velocity boundary layer, while the second peak is caused by the interface between stable and unstably stratified regions. This peak roughly aligns with the position of the minimum in
$\langle T^\prime \rangle$
discussed earlier. As rotation increases, these effects become more pronounced: the first peak shifts closer to the wall as the velocity boundary layer becomes a thin Ekman layer, while the second peak moves outwards as reduced mixing allows the stably stratified layer to extend further into the flow. This showcases the presence of two distinct regions near the bottom plate, the velocity boundary layer and the stably stratified region, and that their interaction is intricate.
3.6. Boundary layer behaviour
Before examining the exact relationships for the dissipation, the extent and type of the boundary layers present in the system must be analysed. We begin with the thermal boundary layers whose thickness
$\delta _\kappa$
is determined using the peak in the
$\langle T^\prime \rangle$
profiles. Figure 9 illustrates the dependence of
$\delta _\kappa$
on
$R$
and
$E$
. In general, increasing
$R$
leads to thinner thermal boundary layers due to enhanced convection. The bottom boundary layer extends further into the flow compared with the top one, as stable stratification inhibits heat transport. While rotation generally increases the boundary layer thickness, some non-monotonic behaviour is observed, likely resulting from the competing effects of rotational suppression of turbulence and improved flow organisation, which enhances heat transport.

Figure 9. Extent of the thermal boundary layers
$\delta _\kappa$
at the bottom (a) and top (b) plates, as well as the ratio between their sizes (c). Symbols as in figure 7.
Turning to the velocity boundary layers, the top boundary layer behaves in a straightforward manner. Its thickness can be defined from the peak in
$\langle u_h^\prime \rangle$
, with its dependence on
$E$
and
$R$
shown in the left panel of figure 10. At high rotation rates, the data show reasonable collapse when plotted as a function of
$E$
, following the expected scaling
$\delta _\nu |_{z=1}\sim \mathcal{O}(E^{1/2})$
. For sufficiently small
$E$
, the dependence on
$R$
drops out as expected
The bottom boundary layer, however, is more difficult to quantify, as the fluctuation peak only becomes distinct at sufficiently small
$E$
. As a result, its thickness can only be precisely determined once the boundary layer has transitioned into an Ekman layer. In the right panel of figure 10, we present measurements of the bottom boundary layer extent for cases where the peak is clearly identifiable. The available data points are limited in quantity, but we can again observe a scaling similar to
$\delta _\nu |_{z=0}\sim \mathcal{O}(E^{1/2})$
, with minimal
$R$
dependence at small
$E$
.

Figure 10. Extent of the (a) and (b) velocity boundary layer
$\delta _\nu |$
as a function of
$E$
.

Figure 11. (a) numerically obtained velocity profile for
$R=10^{10}$
and
$E=10^{-6}$
compared with the theoretical Ekman solution fitted manually and using a least-squares fit. (b) size ratio between thermal and velocity boundary layers at the top plate.
It is important to note that these boundary layers are not sufficiently thin to exhibit a ‘pure’ Ekman dynamics, as seen in rotating RBC at small
$E$
. In the left panel of figure 11, we follow Aguirre Guzmán et al. (Reference Aguirre Guzmán, Madonia, Cheng, Ostilla-Mónico, Clercx and Kunnen2022) and attempt to fit the theoretical solution for the Ekman boundary layer to the profiles obtained. In general, we cannot reproduce the peak, neither using both a least-squares fit to the full profile which fails to capture the size of the peak, nor by choosing the parameters manually to capture the size of the peak as this fails to capture the remaining behaviour.
Finally, the right panel of figure 11 presents the ratio of thermal to velocity boundary layer thicknesses at the top plate. This ratio has been proposed as a diagnostic for the transition to rotation-dominated convection or even the onset of geostrophic turbulence (King et al. Reference King, Stellmach, Noir, Hansen and Aurnou2009; King & Aurnou Reference King and Aurnou2012). However, in our simulations, this cross-over occurs at relatively large values of
$E$
(i.e. weak rotation) compared with RBC. This supports the argument by Kunnen et al. (Reference Kunnen, Ostilla-Mónico, Van Der Poel, Verzicco and Lohse2016) that this ratio alone does not govern the transition to rotation-dominated convection.

Figure 12. Profiles of
$\langle \epsilon _\nu \rangle$
(a,c) and
$\langle \epsilon _\theta \rangle$
(b,d) dissipation rates for several values of
$R$
for the non-rotating case (a,b) and several values of
$E$
for
$R=10^{10}$
(c,d).
3.7. Dissipation rates
Having described the phenomenology of the flow, we conclude our analysis by examining the dissipation rates
$\epsilon _\nu$
and
$\epsilon _\theta$
. These quantities are related to the global quantities
$\overline {\langle wT \rangle }$
and
$\overline {\langle T \rangle }$
via (2.5)–(2.6). Investigating their local behaviour and the relative contributions from the bulk and boundary layers has been a key approach in RBC to estimate the global response of the system and identify the factors limiting heat transport.
The left column of figure 12 shows the profiles of
$\langle \epsilon _\nu \rangle$
for various values of
$R$
and
$E$
. Generally, these profiles feature a bulk region with nearly constant values of
$\epsilon _\nu$
and near-wall regions with local maxima and minima. At the top wall, kinetic energy dissipation shows a peak, consistent with expectations. However, the bottom wall displays a more complex structure: a peak at the wall, and a local minimum just beyond the velocity boundary layer. As rotation increases, this minimum becomes more pronounced, while the dissipation peak at the top boundary layer strengthens due to the intensified kinetic energy dissipation within the Ekman boundary layer.
The right column of figure 12 shows the profiles of
$\langle \epsilon _\theta \rangle$
for the same values of
$R$
and
$E$
. Unlike
$\epsilon _\nu$
,
$\epsilon _\theta$
decreases with increasing
$R$
, reflecting the overall reduction in temperature fluctuations. Local maxima appear across all cases at the plates, and a distinct minimum can be observed farther from the bottom wall. The location of this minimum roughly corresponds to that of
$\langle T^\prime \rangle$
seen in figure 6, which we previously associated with the edge of the stably stratified layer. Increasing rotation does not significantly alter the qualitative shape of the
$\epsilon _\theta$
profiles or the magnitude of the peak at the top plate, although the position of the local minimum shifts slightly.

Figure 13. Absolute (a,c) and relative (b,d) contributions to the viscous dissipation rate
$\overline {\langle \epsilon _\nu \rangle }$
for non-rotating IHC (a,b) and rotating IHC at
$R=10^{10}$
(c,d).
Following Grossmann & Lohse (Reference Grossmann and Lohse2000) and Wang et al. (Reference Wang, Lohse and Shishkina2020), we decompose the dissipation rates
$\overline {\langle \epsilon _\theta \rangle }$
and
$\overline {\langle \epsilon _\nu \rangle }$
into contributions from the bulk and boundary layers to better understand the scaling laws governing
$\overline {\langle T \rangle }$
and
$\overline {\langle wT \rangle }$
. Defining the bottom velocity boundary layer extent is challenging, preventing us from isolating its contribution to
$\overline {\langle \epsilon _\nu \rangle }$
. However, inspection of the
$\langle \epsilon _\nu \rangle$
profiles in figure 12 suggests that dissipation in the bottom boundary layer is negligible compared with the bulk contribution, justifying the assumption that a combined bulk – bottom boundary layer contribution would primarily originate from the bulk region.
Figure 13 shows the absolute and relative contributions of
$\overline {\langle \epsilon _\nu \rangle }$
from the top boundary and from the bulk (including the bottom boundary layer). In non-rotating cases,
$\overline {\langle \epsilon _\nu \rangle }$
increases in the bulk with increasing
$R$
, but decreases at the top BL. The relative contributions indicate that over
$80\,\%$
of
$\overline {\langle \epsilon _\nu \rangle }$
comes from the bulk, suggesting that the dependence of
$\overline {\langle wT \rangle }$
on
$R$
is primarily governed by bulk viscous dissipation.

Figure 14. Absolute (a,c) and relative (b,d) contributions to the thermal dissipation rate
$\overline {\langle \epsilon _\theta \rangle }$
for non-rotating IHC (a,b) and rotating IHC at
$R=10^{10}$
(c,d).
With increasing rotation, the behaviour becomes more complex. While
$\overline {\langle \epsilon _\nu \rangle }$
in the bulk remains nearly constant before decreasing rapidly at small values of
$E$
, dissipation in the top boundary layer increases significantly, reaching relative contributions of over
$20\,\%$
. The maxima in
$\overline {\langle wT \rangle } (E)$
can be directly linked to increases in
$\overline {\langle \epsilon _\nu \rangle }$
within the top boundary layer, preceding the decline of
$\overline {\langle \epsilon _\nu \rangle }$
in the bulk.
A similar analysis for
$\overline {\langle \epsilon _\theta \rangle }$
is shown in figure 14. In the absence of rotation, the dominant contribution comes from the top boundary layer. A rough scaling of
$R^{-0.2}$
is observed for thermal dissipation in this region, indicating that the scaling laws governing
${\overline {\langle T \rangle } }(R)$
are constrained by the flow dynamics in this region. As
$R$
increases, bulk contributions become increasingly significant, suggesting that, for much larger values of
$R$
, the bulk may eventually dominate the behaviour of
$\overline {\langle T \rangle }$
.
The behaviour of
$\overline {\langle \epsilon _\theta \rangle }$
under rotation differs significantly from that of
$\overline {\langle \epsilon _\nu \rangle }$
. While dissipation in the top thermal boundary layers decreases slightly at high rotation rates, this effect is quickly overshadowed by a surge in bulk dissipation, which surpasses boundary layer contributions and limits the decrease of
$\overline {\langle T \rangle }$
. This is another explanation on why
$\overline {\langle T \rangle }$
and
$\overline {\langle wT \rangle }$
exhibit different trends around ‘optimal transport’ from a dissipation perspective, and why
$\overline {\langle T \rangle }$
is more strongly dependent on
$R$
, as seen from figure 5.
Our findings indicate that, for
$\textit{Pr}=1$
,
$\overline {\langle \epsilon _\nu \rangle }$
is primarily controlled by bulk dissipation, whereas
$\overline {\langle \epsilon _\theta \rangle }$
is dominated by dissipation in the top boundary layer for low rotation rates and by bulk dissipation for high rotation rates. This contrasts with the 2-D simulations of Wang et al. (Reference Wang, Lohse and Shishkina2020), which suggest that
$\overline {\langle \epsilon _\nu \rangle }$
is primarily governed by boundary layer contributions. The discrepancy likely arises from the behaviour of IHC in two dimensions, where the presence of a large-scale circulation significantly biases velocity statistics, causing the velocities to be higher and hence resulting in much larger dissipations (Goluskin & Van der Poel Reference Goluskin and Van der Poel2016).
4. Conclusion and outlook
In this study, we have conducted the first 3-D simulations of rotating IHC with isothermal boundaries, where the fluid is cooled from both plates. We demonstrated that the asymptotic stability results derived in Arslan (Reference Arslan2024) accurately describe the behaviour observed in our simulations at
$\textit{Pr}=1$
. We then explored the effects of rotation on IHC over a parameter range of
$R \in [3.16\times 10^5,10^{10}]$
and
$E \in [10^{-6},\infty )$
. Weak rotation organises the plumes originating from the top boundary layer, enhancing heat transport even as the velocity fluctuations are reduced. This effect is reflected on both global quantities
$\overline {\langle T \rangle }$
and
$\overline {\langle wT \rangle }$
, with a more pronounced impact on the latter. We also analysed local flow behaviour, revealing that an Ekman boundary layer forms near the top plate, while near the bottom plate, a complex interaction occurs between the boundary layers and the stably stratified region. This effect is unique to IHC in the IH1 configuration, and is affected by rotation. Finally, we examined the dissipation rates, concluding that the system was dominated by boundary layer dissipation for low rotation rates and by bulk dissipation for high rotation rates for
$\overline {\langle \epsilon _\theta \rangle }$
and by bulk dissipation for
$\overline {\langle \epsilon _\nu \rangle }$
, a significant difference from the2-Dsimulations.
This study has two main limitations. First, we have fixed
$\textit{Pr}=1$
, and further exploration of different Prandtl numbers would provide deeper insight into the system’s complex dynamics, especially into the interplay between boundary layer(s) and stably stratified regions. System asymmetry has been shown to substantially depend on
$\textit{Pr}$
(Goluskin & Van der Poel Reference Goluskin and Van der Poel2016; Wang et al. Reference Wang, Lohse and Shishkina2020). Second, our thermal driving is not strong enough to access the fully rotation-dominated regime, where
$E$
can be varied without significantly suppressing convection. Future simulations extending into this regime will be essential for a more comprehensive understanding of rotating IHC. In an extended parameter range the sharpness of rigorous bounds can be further tested and possibly new flow states can be probed for rotating IHC.
Acknowledgements
R.O.M. acknowledges support from the Emergia Program of the Junta de Andalucía (Spain). A.A. acknowledges funding from the ERC (agreement no. 833848-UEMHP) under the Horizon 2020 program and the SNSF (grant number 219247) under the MINT 2023 call. The authors thank R. Kunnen for helpful comments. We also thank the Systems Unit of the Information Systems Area of the University of Cadiz for computer resources and technical support.
Declaration of interests
The authors report no conflict of interest.
Appendix. Summary of results
Table 1 presents a summary of all numerical results.
Table 1. For caption see next page.
