1 Introduction
Rotation and thermal forcing are essential ingredients that influence the dynamics of stellar and planetary interiors, planetary atmospheres, stars and terrestrial oceans. Examples include solar and stellar envelopes (Miesch Reference Miesch2005; Gastine, Heimpel & Wicht Reference Gastine, Heimpel and Wicht2014), planetary atmospheres and interiors (Aurnou et al.
Reference Aurnou, Calkins, Cheng, Julien, King, Nieves, Soderlund and Stellmach2015; Heimpel, Gastine & Wicht Reference Heimpel, Gastine and Wicht2016) and deep ocean convection on Earth (Marshall & Schott Reference Marshall and Schott1999). The quintessential paradigm for investigating the fundamentals of rotating thermally forced flows is rotating Rayleigh–Bénard convection (RRBC) about a vertical axis, i.e. convection in a layer of Boussinesq fluid confined between flat, horizontal, rigidly rotating upper and lower boundaries that sustain a destabilizing temperature jump
$\unicode[STIX]{x0394}T>0$
. Flows are characterized by the non-dimensional Rayleigh (
$Ra$
), convective Rossby (
$Ro$
) and Ekman (
$E$
) numbers defined as
$$\begin{eqnarray}Ra=\frac{g\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}TH^{3}}{\unicode[STIX]{x1D708}\unicode[STIX]{x1D705}},\quad Ro=\sqrt{\frac{g\unicode[STIX]{x1D6FC}\unicode[STIX]{x0394}T}{4\unicode[STIX]{x1D6FA}^{2}H}},\quad E=Ro\sqrt{\frac{Pr}{Ra}}.\end{eqnarray}$$
These, respectively, measure the strength of thermal forcing, the strength of system rotation
$\unicode[STIX]{x1D6FA}$
relative to thermal forcing, and the importance of viscous diffusion relative to
$\unicode[STIX]{x1D6FA}$
. Here
$g$
denotes gravity,
$H$
the layer depth and
$\unicode[STIX]{x1D6FC}$
the thermal expansion coefficient. The Prandtl number,
$\mathit{Pr}=\unicode[STIX]{x1D708}/\unicode[STIX]{x1D705}$
, is the ratio of viscous to thermal diffusion rates.
Geophysical and astrophysical flows are often turbulent and rotationally constrained. This is characterized by high Rayleigh numbers
$Ra\gg 1$
and low Rossby numbers
$Ro\ll 1$
, implying small Ekman numbers
$E\ll 1$
. For example, estimates for these parameters in Earth’s outer core are
$Ra\approx 10^{20}$
,
$Ro\approx 10^{-6}$
and
$E\approx 10^{-15}$
(Schubert & Soderlund Reference Schubert and Soderlund2011). Utilization of isotropic turbulence theory for an order-of-magnitude estimate predicts that the ratio between the largest and smallest length scales in a turbulent flow scales as
${\sim}Ra^{3/8}$
(Pope Reference Pope2000) or
$O(10^{7})$
in the case of the outer core. The low
$Ro$
indicates geostrophy, a primary hydrodynamic balance between the pressure gradient force and the Coriolis force. This balance further increases the complexity, with the presence of fast inertial waves, thin Ekman boundary layers adjacent to no-slip surfaces, and disparate time scales characterizing the evolution of mean and fluctuating temperatures. The extreme parameters and their associated properties described above are substantially beyond the spatiotemporal resolution capabilities for direct numerical simulations (DNS) and the thermal and mechanical capabilities of laboratory experiments.
While the extreme parameter values present a challenge for DNS (Favier, Silvers & Proctor Reference Favier, Silvers and Proctor2014; Guervilly, Hughes & Jones Reference Guervilly, Hughes and Jones2014; Stellmach et al.
Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014; Horn & Shishkina Reference Horn and Shishkina2015), they allow for the use of perturbation methods that rely on the small size of
$Ro$
as an asymptotic expansion parameter. The reduced non-hydrostatic quasi-geostrophic model (NH-QGM) (Julien, Knobloch & Werne Reference Julien, Knobloch and Werne1998; Julien et al.
Reference Julien, Knobloch, Milliff and Werne2006; Sprague et al.
Reference Sprague, Julien, Knobloch and Werne2006) was developed as a simplification, or reduction, of the Boussinesq equations at low
$Ro$
. The model implicitly enforces geostrophy as a pointwise constraint and is asymptotically exact for RRBC in the presence of stress-free boundary conditions as
$Ro\sim E^{1/3}\rightarrow 0$
. Of particular interest is the dependence of the heat transport, as measured by the non-dimensional Nusselt number
$Nu=qH/(\unicode[STIX]{x1D70C}_{0}c_{p}\unicode[STIX]{x1D705}\unicode[STIX]{x0394}T)$
, on the external control parameters
$(Ra,E,Ro)$
and the molecular fluid property
$Pr$
; here
$q$
is the heat flux and
$c_{p}$
is the volumetric heat capacity. Importantly,
$Nu$
characterizes the fluid state and transitions between fluid states, and is commonly reported due to the relative ease with which it can be measured. To validate results obtained from the NH-QGM, comparisons were made with DNS of the Boussinesq equations at
$E=10^{-7}$
(Stellmach et al.
Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014). Figure 1(a) illustrates the high level of agreement in the
$Nu$
–
$Ra$
relation for stress-free boundaries with the asymptotic scaling of
${\sim}3/2$
.
Comparison of DNS for RRBC (filled symbols) at
$E=10^{-7}$
and the NH-QGM (open symbols) for (a) stress-free and (b) no-slip boundary conditions. The plots show the Nusselt number
$Nu$
as a function of reduced Rayleigh number
$\widetilde{Ra}=Ra\,E^{4/3}$
for Prandtl numbers
$Pr$
of 1 (squares), 3 (circles) and 7 (diamonds). The two methods show good quantitative agreement in the stress-free case (a) but a large discrepancy for the no-slip case (b). Laboratory data for
$5.79\leqslant Pr\leqslant 6.2$
are also included (asterisks). Plots are reproduced from Stellmach et al. (Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014) and Cheng et al. (Reference Cheng, Stellmach, Ribeiro, Grannan, King and Aurnou2015), with an updated plot (b) including higher-
$\widetilde{Ra}$
runs.
For the case of no-slip boundaries (see figure 1
b), substantial disagreement is found in comparison with the NH-QGM results for the
$Nu$
–
$Ra$
relation. DNS and laboratory data in water (
$Pr\sim 7$
) show markedly increased values of
$Nu$
for all values of
$Ra$
, with a steep scaling exponent of
${\gtrsim}3$
(Stellmach et al.
Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014; Cheng et al.
Reference Cheng, Stellmach, Ribeiro, Grannan, King and Aurnou2015). While there is qualitative agreement between DNS and experiments, differences in the Nusselt-number results can be attributed to differing Prandtl numbers for water and different geometries (DNS are performed in a horizontally periodic box with
$Pr=7$
, whereas experiments are performed in cylindrical containers with
$5.79\leqslant Pr\leqslant 6.2$
). More specifically, it is shown that the divergence in
$Nu$
between the NH-QGM, where Ekman boundary layers are presumed passive, and DNS data grows as the Ekman number
$E$
is lowered in the latter. A recent study by Cheng et al. (Reference Cheng, Stellmach, Ribeiro, Grannan, King and Aurnou2015) (see their figure 6) indicates that the heat transport exponent
$\unicode[STIX]{x1D6FC}$
in the
$Nu$
–
$Ra$
relation
$Nu\propto Ra^{\unicode[STIX]{x1D6FC}}$
increases monotonically from 1 to 3.6 as
$E$
decreases to
${\sim}10^{-7}$
. By replacing the no-slip boundaries with parametrized Ekman pumping boundary conditions and establishing quantitative agreement, Stellmach et al. (Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014) clearly demonstrate that the root cause of this divergence is the presence of Ekman boundary layers. This was an important finding given that the boundary layer width and the associated vertical pumping (or suction) velocities, which respectively scale as
$E^{1/2}H$
and
$w_{E}=E^{1/2}\,\widehat{\boldsymbol{z}}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\times \boldsymbol{u}/\sqrt{2}$
, have diminishing values as
$E\rightarrow 0$
. (The latter formula for
$w_{E}$
represents an exact parametrization of the laminar Ekman boundary layer (Greenspan Reference Greenspan1968).) Indeed, the asymptotic NH-QGM assumes that Ekman boundary layers become increasingly insignificant as
$E\rightarrow 0$
and are filtered from the dynamics, thus implicitly assuming impenetrable stress-free boundaries (Julien et al.
Reference Julien, Knobloch and Werne1998; Sprague et al.
Reference Sprague, Julien, Knobloch and Werne2006).
An asymptotic analysis of the boundary layer regions of RRBC using the Boussinesq equations was recently reported in Julien et al. (Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016). Specifically, it was established that the effects of Ekman pumping become significant within the rotationally constrained regime when
$E^{-1/9}\lesssim Ra\,E^{4/3}\lesssim E^{-1/3}$
, i.e. for sufficiently strong thermal forcing. Consequently, the NH-QGM was updated to include the effect of no-slip boundaries and the resultant Ekman pumping through the use of the parametrized pumping boundary conditions and the introduction of a thermal wind layer. The thermal wind layer was shown to be capable of supporting the enhanced thermal fluctuations necessary for the amplification of heat transport. Both effects are explicitly controlled by the magnitude of the Ekman number
$E$
. Thus, the model, called the composite non-hydrostatic quasi-geostrophic model (CNH-QGM), now includes the Ekman number
$E$
as an explicit control parameter. The capabilities of the composite model were explored in Julien et al. (Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016) by locating and tracking exact single-mode (i.e. single horizontal wavenumber) solutions as a function of reduced Rayleigh number. Single-mode solutions represent an unstable solution branch of the CNH-QGM; despite this, their presence strongly influences the path of the realized solution in phase space and therefore they capture and characterize the qualitative behaviour of the enhanced heat transport in the fully nonlinear regime.
In this paper, we extend the work of Julien et al. (Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016) and explore simulations of the asymptotically reduced CNH-QGM that, similar to DNS, capture all dynamically active spatial scales and thus the realized nonlinear solution. The explicit appearance of the Ekman number
$E$
permits a direct comparison with DNS investigations of the Boussinesq equation for RRBC and provides the capability to assess both the fidelity of the reduced model and the degree to which the DNS has entered the asymptotic regime. In § 2, the CNH-QGM derived in Julien et al. (Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016) is documented and the multi-layer structure unique to the presence of no-slip boundaries identified. In § 3 the numerical method used to solve the CNH-QGM is discussed and results obtained from a suite of simulations presented. Concluding remarks are given in § 4.
2 Composite model
The full derivation and discussion of the CNH-QGM is presented in Julien et al. (Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016); here we simply provide an overview of their asymptotic model. The leading-order velocity field outside the Ekman layer,
$\boldsymbol{u}=(\boldsymbol{u}_{\bot },W)$
, is in geostrophic balance, where
$\widehat{\boldsymbol{z}}\times \boldsymbol{u}_{\bot }=-\unicode[STIX]{x1D735}_{\bot }p$
. It follows that the horizontal velocity field solved by
$\boldsymbol{u}_{\bot }=(-\unicode[STIX]{x2202}_{y},\unicode[STIX]{x2202}_{x},0)\unicode[STIX]{x1D6F9}$
is non-divergent, where the pressure
$p=\unicode[STIX]{x1D6F9}$
becomes the geostrophic streamfunction. On defining
$\unicode[STIX]{x1D716}=E^{1/3}=Ro$
, the reduced CNH-QGM system governing the evolution of the fluid is given by
$$\begin{eqnarray}\displaystyle & \displaystyle \text{D}_{t}^{\bot }\unicode[STIX]{x1D701}-\unicode[STIX]{x2202}_{Z}W=\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D701}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{D}_{t}^{\bot }W+\unicode[STIX]{x2202}_{Z}\unicode[STIX]{x1D6F9}=\frac{\widetilde{Ra}}{Pr}\unicode[STIX]{x1D703}+\unicode[STIX]{x1D6FB}_{\bot }^{2}W, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \text{D}_{t}^{\bot }\unicode[STIX]{x1D703}+W\unicode[STIX]{x2202}_{Z}\overline{\unicode[STIX]{x1D6E9}}+\unicode[STIX]{x1D716}(\unicode[STIX]{x1D735}_{\bot }\boldsymbol{\cdot }(\boldsymbol{u}_{1\bot }^{ag}\unicode[STIX]{x1D703})+\unicode[STIX]{x2202}_{Z}(W\unicode[STIX]{x1D703}-\overline{W\unicode[STIX]{x1D703}}))=\frac{1}{Pr}(\unicode[STIX]{x1D6FB}_{\bot }^{2}+\unicode[STIX]{x1D716}^{2}\unicode[STIX]{x2202}_{ZZ})\unicode[STIX]{x1D703}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle \unicode[STIX]{x2202}_{\unicode[STIX]{x1D70F}}\overline{\unicode[STIX]{x1D6E9}}+\unicode[STIX]{x2202}_{Z}(\overline{W\unicode[STIX]{x1D703}})=\frac{1}{Pr}\unicode[STIX]{x2202}_{ZZ}\overline{\unicode[STIX]{x1D6E9}}, & \displaystyle\end{eqnarray}$$
capturing, respectively, the evolution of vertical vorticity
$\unicode[STIX]{x1D701}=-\unicode[STIX]{x1D6FB}_{\bot }^{2}\unicode[STIX]{x1D6F9}$
, vertical velocity
$W$
and temperature
$\unicode[STIX]{x1D6E9}=\overline{\unicode[STIX]{x1D6E9}}+\unicode[STIX]{x1D716}\unicode[STIX]{x1D703}$
at reduced Rayleigh number
$\widetilde{Ra}=Ra\,\unicode[STIX]{x1D716}^{4}$
for a given Prandtl number
$Pr$
. The temperature is decomposed into a mean (horizontally averaged) component
$\overline{\unicode[STIX]{x1D6E9}}$
evolving on a slow time scale
$\unicode[STIX]{x1D70F}=\unicode[STIX]{x1D716}^{2}t$
and a small fluctuating component
$\unicode[STIX]{x1D703}$
characteristic of RRBC. The overbar indicates the horizontal and fast time averaging operator given by
$$\begin{eqnarray}\overline{f}(Z,T)=\lim _{\unicode[STIX]{x1D70F},A\rightarrow \infty }\frac{1}{\unicode[STIX]{x1D70F}A}\int _{A}f(\boldsymbol{x},Z,t,T)\,\text{d}\boldsymbol{x}\,\text{d}t,\end{eqnarray}$$
where
$A$
is the area of the box. Here
$\text{D}_{t}^{\bot }=\unicode[STIX]{x2202}_{t}+\boldsymbol{u}_{\bot }\boldsymbol{\cdot }\unicode[STIX]{x1D735}_{\bot }$
denotes the material derivative with horizontal geostrophic advection and
$\boldsymbol{u}_{1\bot }^{ag}$
in (2.3) denotes the ageostrophic velocity field determined through the three-dimensional (3D) incompressibility condition
$$\begin{eqnarray}\unicode[STIX]{x1D735}_{\bot }\boldsymbol{\cdot }\boldsymbol{u}_{1\bot }^{ag}+\unicode[STIX]{x2202}_{Z}W=0.\end{eqnarray}$$
The system is solved for the RRBC configuration with fixed temperature boundaries
$$\begin{eqnarray}\overline{\unicode[STIX]{x1D6E9}}=1,\quad \unicode[STIX]{x1D703}=0\quad \text{at}~Z=0\quad \text{and}\quad \overline{\unicode[STIX]{x1D6E9}}=0,\quad \unicode[STIX]{x1D703}=0\quad \text{at}~Z=1.\end{eqnarray}$$
Velocity boundary conditions are either parametrized Ekman pumping for no-slip (NS) boundaries,
$$\begin{eqnarray}\text{NS{-}NS:}\quad W=\unicode[STIX]{x1D716}^{1/2}\frac{\unicode[STIX]{x1D701}}{\sqrt{2}}\quad \text{at}~Z=0\quad \text{and}\quad W=-\unicode[STIX]{x1D716}^{1/2}\frac{\unicode[STIX]{x1D701}}{\sqrt{2}}\quad \text{at}~Z=1,\end{eqnarray}$$
or implicit impenetrable stress-free (SF) boundaries
$$\begin{eqnarray}\text{SF{-}SF:}\quad W=0\quad \text{at}~Z=0,1,\end{eqnarray}$$
with
$\unicode[STIX]{x1D716}=0$
or mixed (NS–SF) combinations thereof (e.g. NS at
$Z=0$
and SF at
$Z=1$
). Equations (2.1)–(2.7) with one of the above choices of velocity boundary conditions form a closed reduced system that filters fast inertial waves occurring on
$O(\unicode[STIX]{x1D6FA}^{-1})$
time scales and
$O(L)$
spatial scales. The imposition of strong spatial anisotropy by non-dimensionalizing vertical and horizontal scales with
$H$
and
$L=\unicode[STIX]{x1D716}H$
is a key step in the filtering that gives rise to the reduced dynamics (Julien & Knobloch Reference Julien and Knobloch2007). Velocity and time are non-dimensionalized using
$\unicode[STIX]{x1D708}/L$
and
$L^{2}/\unicode[STIX]{x1D708}$
and temperature is non-dimensionalized by
$\unicode[STIX]{x0394}T$
.
The evolution of the barotropic vertical vorticity
$\langle \unicode[STIX]{x1D701}\rangle$
is a key quantity in the understanding of the inverse energy cascade recently observed in reduced simulations of RRBC (Rubio et al.
Reference Rubio, Julien, Knobloch and Weiss2014) as well as DNS (Favier et al.
Reference Favier, Silvers and Proctor2014; Guervilly et al.
Reference Guervilly, Hughes and Jones2014; Stellmach et al.
Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014). Depth averaging the vertical vorticity equation (2.1) according to
$\langle ~\rangle =\int _{0}^{1}\text{d}Z$
, with the vortical variables decomposed as
$f=\langle f\rangle +f^{\prime }$
, gives
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}\langle \unicode[STIX]{x1D701}\rangle =-\text{J}[\langle \unicode[STIX]{x1D6F9}\rangle ,\langle \unicode[STIX]{x1D701}\rangle ]-\langle \text{J}[\unicode[STIX]{x1D6F9}^{\prime },\unicode[STIX]{x1D701}^{\prime }]\rangle -\frac{\unicode[STIX]{x1D716}^{1/2}}{\sqrt{2}}(\unicode[STIX]{x1D701}^{\prime }(0)+\unicode[STIX]{x1D701}^{\prime }(1))-\frac{2\unicode[STIX]{x1D716}^{1/2}}{\sqrt{2}}\langle \unicode[STIX]{x1D701}\rangle +\unicode[STIX]{x1D6FB}_{\bot }^{2}\langle \unicode[STIX]{x1D701}\rangle ,\end{eqnarray}$$
indicating that the growth of the barotropic vertical vorticity
$\langle \unicode[STIX]{x1D701}\rangle$
is dependent on the net balance between terms on the right-hand side for, respectively, the barotropic self-interaction, baroclinic (i.e. convective) forcing, barotropic and baroclinic Ekman friction, and barotropic viscous dissipation. Here
$\text{J}[f,g]=\unicode[STIX]{x2202}_{x}\,f\unicode[STIX]{x2202}_{y}g-\unicode[STIX]{x2202}_{y}\,f\unicode[STIX]{x2202}_{x}g$
denotes the Jacobian for the horizontal advection of
$g$
, when
$f$
is either
$\langle \unicode[STIX]{x1D6F9}\rangle$
or
$\unicode[STIX]{x1D6F9}^{\prime }$
.
We spectrally analyse (2.10) where barotropic and baroclinic fields are decomposed according to the discrete Fourier transform,
$$\begin{eqnarray}\langle f\rangle =\mathop{\sum }_{\boldsymbol{k}}\langle \,\hat{f}_{\boldsymbol{k}}\rangle \text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}_{\bot }},\quad f^{\prime }=\mathop{\sum }_{\boldsymbol{k}}\hat{f}_{\boldsymbol{k}}^{\prime }(Z)\text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}_{\bot }},\end{eqnarray}$$
with
$\boldsymbol{k}=(k_{x},k_{y})$
. Specifically, we consider the evolution of the spectral barotropic kinetic energy
$E_{bt}(\boldsymbol{k})=(|\boldsymbol{k}|^{2}|\langle \hat{\unicode[STIX]{x1D6F9}}_{\boldsymbol{k}}\rangle |^{2})/2$
at each wavenumber
$\boldsymbol{k}$
,
$$\begin{eqnarray}\unicode[STIX]{x2202}_{t}E_{bt}(\boldsymbol{k})=T_{\boldsymbol{k}}+F_{\boldsymbol{k}}+Ef_{\boldsymbol{k}}^{bc}+Ef_{\boldsymbol{ k}}^{bt}+D_{\boldsymbol{ k}},\end{eqnarray}$$
obtained by multiplication of the Fourier decomposition of (2.10) by
$-\langle \hat{\unicode[STIX]{x1D6F9}}_{\boldsymbol{k}}\rangle \text{e}^{\text{i}\boldsymbol{k}\boldsymbol{\cdot }\boldsymbol{x}_{\bot }}$
, followed by integration over physical space. Here
$T_{\boldsymbol{k}}=\sum _{\boldsymbol{p},\boldsymbol{q}}T_{\boldsymbol{k}\boldsymbol{p}\boldsymbol{q}}$
and
$F_{\boldsymbol{k}}=\sum _{\boldsymbol{p},\boldsymbol{q}}F_{\boldsymbol{k}\boldsymbol{p}\boldsymbol{q}}$
describe how power is converted to barotropic Fourier amplitude
$\langle \hat{\unicode[STIX]{x1D6F9}}_{\boldsymbol{k}}\rangle$
through the barotropic and baroclinic triadic interactions involving modes
$\boldsymbol{p}$
and
$\boldsymbol{q}$
such that
$\boldsymbol{k}+\boldsymbol{q}+\boldsymbol{p}=\mathbf{0}$
(Favier et al.
Reference Favier, Silvers and Proctor2014; Rubio et al.
Reference Rubio, Julien, Knobloch and Weiss2014), where
$$\begin{eqnarray}\displaystyle & \displaystyle T_{\boldsymbol{k}\boldsymbol{p}\boldsymbol{q}}=b_{\boldsymbol{p}\boldsymbol{q}}\,\text{Re}[\langle \hat{\unicode[STIX]{x1D6F9}}_{\boldsymbol{k}}\rangle \langle \hat{\unicode[STIX]{x1D6F9}}_{\boldsymbol{p}}\rangle \langle \hat{\unicode[STIX]{x1D6F9}}_{\boldsymbol{q}}\rangle ]\unicode[STIX]{x1D6FF}_{\boldsymbol{k}+\boldsymbol{p}+\boldsymbol{q},\mathbf{0}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle F_{\boldsymbol{k}\boldsymbol{p}\boldsymbol{q}}=b_{\boldsymbol{p}\boldsymbol{q}}\,\text{Re}[\langle \hat{\unicode[STIX]{x1D6F9}}_{\boldsymbol{k}}\rangle \langle \hat{\unicode[STIX]{x1D6F9}}_{\boldsymbol{p}}^{\prime }\hat{\unicode[STIX]{x1D6F9}}_{\boldsymbol{q}}^{\prime }\rangle ]\unicode[STIX]{x1D6FF}_{\boldsymbol{k}+\boldsymbol{p}+\boldsymbol{q},\mathbf{0}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle b_{\boldsymbol{p}\boldsymbol{q}}=b_{\boldsymbol{q}\boldsymbol{p}}={\textstyle \frac{1}{2}}(p^{2}-q^{2})(p_{x}q_{y}-p_{y}q_{x}), & \displaystyle\end{eqnarray}$$
$\text{Re}[\ldots \,]$
denotes the real part and
$\unicode[STIX]{x1D6FF}_{\boldsymbol{k}+\boldsymbol{p}+\boldsymbol{q},\mathbf{0}}$
is the Kronecker delta function. The latter terms in (2.12) are defined as
$$\begin{eqnarray}\displaystyle & \displaystyle D_{\boldsymbol{k}}=-|\boldsymbol{k}|^{2}E_{bt}(\boldsymbol{k}), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle Ef_{\boldsymbol{k}}^{bc}=-\frac{1}{\sqrt{2}}\unicode[STIX]{x1D700}^{1/2}(E_{bc}(\boldsymbol{k},Z=0)+E_{bc}(\boldsymbol{k},Z=1)), & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle Ef_{\boldsymbol{k}}^{bt}=-\frac{2}{\sqrt{2}}\unicode[STIX]{x1D716}^{1/2}E_{bt}(\boldsymbol{k}), & \displaystyle\end{eqnarray}$$
where
$$\begin{eqnarray}\displaystyle E_{bc}(\boldsymbol{k},Z)=\langle \hat{\unicode[STIX]{x1D6F9}}_{\boldsymbol{k}}\rangle \hat{\unicode[STIX]{x1D6F9}}_{\boldsymbol{k}}^{^{\prime }\ast }(Z)+\text{c.c.} & & \displaystyle\end{eqnarray}$$
Here the asterisk denotes the conjugate. We find that baroclinic Ekman friction is negligible compared to barotropic Ekman friction, i.e.
$|Ef_{\boldsymbol{k}}^{bc}|\ll |Ef_{\boldsymbol{k}}^{bt}|$
.
From (2.12), one can consider the evolution of the one-dimensional (1D) energy spectrum, defined as
$$\begin{eqnarray}\displaystyle E_{bt}(k)=\int _{0}^{2\unicode[STIX]{x03C0}}\!E_{bt}(\boldsymbol{k})k\,\text{d}\unicode[STIX]{x1D719}_{\boldsymbol{k}}, & & \displaystyle\end{eqnarray}$$
where, in the polar representation,
$k=|\boldsymbol{k}|$
and
$\tan \unicode[STIX]{x1D719}_{\boldsymbol{k}}=k_{y}/k_{x}$
. Also of particular interest are the transfer maps
$$\begin{eqnarray}\displaystyle & \displaystyle T_{kp}=\int k\,\text{d}\unicode[STIX]{x1D719}_{\boldsymbol{k}}\int p\,\text{d}\unicode[STIX]{x1D719}_{\boldsymbol{p}}\mathop{\sum }_{\boldsymbol{q}}T_{\boldsymbol{k}\boldsymbol{p}\boldsymbol{q}}, & \displaystyle\end{eqnarray}$$
$$\begin{eqnarray}\displaystyle & \displaystyle F_{kp}=\int k\,\text{d}\unicode[STIX]{x1D719}_{\boldsymbol{k}}\int p\,\text{d}\unicode[STIX]{x1D719}_{\boldsymbol{p}}\mathop{\sum }_{\boldsymbol{q}}F_{\boldsymbol{k}\boldsymbol{p}\boldsymbol{q}}, & \displaystyle\end{eqnarray}$$
detailing the transfer of energy from all wavenumbers of magnitude
$p$
to magnitude
$k$
.
3 Results
In this section, we present our findings obtained from the suite of numerical simulations documented in table 1. These span a range of reduced Rayleigh numbers
$Ra\,E^{4/3}\leqslant 100$
and Prandtl numbers
$1/2\leqslant Pr\leqslant 7$
. All simulations of the CNH-QGM and the DNS employ a spectral Fourier discretization in the periodic horizontal directions, a sparse Chebyshev-tau discretization in the vertical direction (Julien & Watson Reference Julien and Watson2009), and a mixed semi-implicit explicit (IMEX) time-stepping scheme. The DNS code utilizes the second-order IMEX Adams–Bashforth/backward-differentiation scheme (AB2/BDF2) (Stellmach & Hansen Reference Stellmach and Hansen2008) and the CNH-QGM employs the third-order IMEX Runge–Kutta scheme (RK3) (Spalart, Moser & Rogers Reference Spalart, Moser and Rogers1991). The aspect ratio for all simulations is selected to permit a sufficiently large sampling of linearly unstable convection modes to ensure convergence in the averaged statistics and to permit adequate room to capture the inverse energy cascade. The grid resolutions in all simulations are selected to adequately resolve dynamics on the most restrictive scales. In the vertical direction, this is controlled by the
$O(E^{1/3}H)$
thermal wind layer for the CNH-QGM and the
$O(E^{1/2}H)$
Ekman layer in the presence of no-slip boundaries for the DNS. Specifically, the resolutions selected ensure a minimum of eight grid points in the
$O(E^{1/3}H)$
thermal wind layer and a minimum of eight grid points in the Ekman boundary layer, which is calculated using the distance
$z_{E}=(3/4)\unicode[STIX]{x03C0}\sqrt{2E}$
to the first maxima in the vertical spiral profile (Stellmach et al.
Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014). The resolution in the horizontal direction is controlled by the viscous dissipation, which is estimated by the Kolmogorov length scale
$$\begin{eqnarray}L_{kol}=(\unicode[STIX]{x1D708}^{3}/\unicode[STIX]{x1D700})^{1/4}.\end{eqnarray}$$
The energy dissipation rate
$\unicode[STIX]{x1D700}$
is estimated as the minimum of the horizontally averaged square of the viscous strain rate. The horizontal resolution ensures that the ratio of the smallest resolved length scale to the Kolmogorov length scale is below 2 (see column 5, table 1). This criteria was adequate for convergence in the Nusselt number. Moreover, under this criterion, simulations with sufficiently large box sizes (i.e. greater than
$10L_{c}\times 10L_{c}\times 1$
) resulted in the same
$Nu$
values as those reported in table 1.
Model parameters used in NS–NS simulations for the CNH-QGM and DNS at
$E=10^{-7}$
. Horizontal box dimensions are a function of the critical length scale for stress-free convection,
$L_{c}=4.8154$
for
$Pr\geqslant 1$
and
$L_{c}=6.0292$
for
$Pr<1$
(Chandrasekhar Reference Chandrasekhar1961). Box dimensions marked by an asterisk were simulated with box size
$20L_{c}\times 20L_{c}\times 1$
for CNH-QGM. Simulation resolutions were verified a posteriori. The vertical resolutions ensure at least eight points in both the
$O(E^{1/3}H)$
thermal wind layer for the CNH-QGM and the
$O(E^{1/2}H)$
Ekman layer for DNS. For the horizontal direction, the ratio of the Kolmogorov length scale to grid resolution was bounded below 2. Time averaging for the
$Nu$
results used intervals that ranged from
$3.6\times 10^{-2}$
to
$8.5\times 10^{-4}$
for DNS and from
$2\times 10^{-1}$
to
$1\times 10^{-3}$
for CNH-QGM.
The time resolutions utilized in the simulations of the CNH-QGM and the DNS are displayed in table 1 (columns 8 and 9). Each is limited by the Courant–Friedrichs–Lewy (CFL) constraint imposed by the explicit treatment of the nonlinear advection terms. In the case of DNS, geostrophy and the separation of time scales between the mean and fluctuating temperature fields are not implicitly imposed, which therefore increases the degree of stiffness of any IMEX time integration method. As is evident from the comparison, the asymptotic approach permits a temporal rescaling that automatically factors an
$O(E^{-2/3}Pr)$
magnitude from the time step (see column 7 of table 1). Hence, there is no explicit Ekman-number dependence on the time step and the increase in allowable non-dimensional time steps demonstrates a numerical advantage over DNS.
Owing to spatiotemporal constraints, DNS with
$E=10^{-7}$
still remains the lower bound for investigations of rotating convection (Stellmach et al.
Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014). Similar bounds are also present in laboratory experiments because of thermal and mechanical constraints (Ecke & Niemela Reference Ecke and Niemela2014; Cheng et al.
Reference Cheng, Stellmach, Ribeiro, Grannan, King and Aurnou2015). The CNH-QGM is capable of explorations at much lower
$E$
as a result of the asymptotic design. Specifically, we note that
$E=0$
represents the case of purely stress-free boundaries considered elsewhere (Sprague et al.
Reference Sprague, Julien, Knobloch and Werne2006; Julien et al.
Reference Julien, Knobloch, Rubio and Vasil2012a
,Reference Julien, Rubio, Grooms and Knobloch
b
; Rubio et al.
Reference Rubio, Julien, Knobloch and Weiss2014). In the following three subsections, the fidelity of the asymptotic model is explored via qualitative and quantitative comparisons with DNS imposing no-slip boundaries at
$E=10^{-7}$
. In § 3.4, results are reported from simulations with the Ekman number varying below
$E=10^{-7}$
.
3.1 Flow morphologies
The four convective flow morphologies and their domains of existence as a function of
$\widetilde{Ra}$
and
$Pr$
has been well documented using the NH-QGM for RRBC in the presence of top and bottom stress-free (SF–SF) boundaries (cf. figure 2 in Julien et al. (Reference Julien, Rubio, Grooms and Knobloch2012b
)). Figure 2 illustrates that simulations of the CNH-QGM using the parametrized no-slip (NS–NS) conditions yield the same four geostrophic regimes. Ordered by increasing
$\widetilde{Ra}$
, these include a cellular regime, a regime of shielded convective Taylor columns (CTCs), a plume regime and ultimately the geostrophic turbulence (GT) regime. The DNS with spatially resolved Ekman boundary layers at similar values of
$\widetilde{Ra}$
,
$\mathit{Pr}$
and
$E$
yield comparable visualization results (cf. figure 2 in Stellmach et al. (Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014)), thereby establishing a first marker for the fidelity of the CNH-QGM.
Perspective view of the fluctuating temperature field
$\unicode[STIX]{x1D703}$
, illustrating the four geostrophically balanced flow morphologies obtained from simulations of the CNH-QGE with NS–NS boundaries at
$E=10^{-7}$
(or
$\unicode[STIX]{x1D716}=E^{1/3}=0.0046$
). The plots show: (a) cellular convection at
$\widetilde{Ra}=10$
,
$\mathit{Pr}=7$
; (b) convective Taylor columns at
$\widetilde{Ra}=30$
,
$\mathit{Pr}=7$
; (c) plumes at
$\widetilde{Ra}=70$
,
$\mathit{Pr}=7$
; and (d) geostrophic turbulence at
$\widetilde{Ra}=90$
,
$\mathit{Pr}=1$
. These results should be compared with figure 2 in Stellmach et al. (Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014).
3.2 Heat transport
A comparison of the
$Nu$
–
$\widetilde{Ra}$
heat transport relations at
$E=10^{-7}$
for DNS and CNH-QGM is shown in figure 3 for
$Pr=7$
(a),
$Pr=3$
(b) and
$Pr=1/2$
(c). Comparisons with
$Pr=1$
are included in each plot for reference. Confirming the findings of Stellmach et al. (Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014) and Cheng et al. (Reference Cheng, Stellmach, Ribeiro, Grannan, King and Aurnou2015) (see figure 1), Ekman pumping results in a significantly enhanced
$Nu$
for all
$\widetilde{Ra}$
compared to simulations with the SF–SF boundaries. Moreover, as with the stress-free simulations (Julien et al.
Reference Julien, Rubio, Grooms and Knobloch2012b
; Nieves, Rubio & Julien Reference Nieves, Rubio and Julien2014), transitions between regimes are observable in figure 3 by abrupt changes in the local heat transport exponent
$\unicode[STIX]{x1D6FD}=\text{d}\ln Nu/\text{d}\ln \widetilde{Ra}$
(reference slopes are provided to guide the eye).
(a–c) Comparison of heat transport
$Nu$
as a function of
$\widetilde{Ra}$
between DNS (filled symbols) and the CNH-QGM (open symbols) with NS–NS boundaries at
$E=10^{-7}$
. The plots show the Nusselt number as a function of (a)
$\mathit{Pr}=7$
(diamonds), (b)
$\mathit{Pr}=3$
(circles), (c)
$\mathit{Pr}=1/2$
(triangles) and (a–c)
$\mathit{Pr}=1$
(squares). Data from experiments by Cheng et al. (Reference Cheng, Stellmach, Ribeiro, Grannan, King and Aurnou2015) for
$Pr\approx 7$
are marked by asterisks in (a). Note the increased level of heat transport compared with the stress-free results presented in figure 1. Reference slopes are included as a guide. (d) The relative percentage errors of the Nusselt number, calculated as
$100\,\%\times (|Nu_{CNH}-Nu_{DNS}|/Nu_{CNH})$
, for all
$Pr$
values as a function of
$\widetilde{Ra}$
.
For
$Pr>1$
, the heat transport law (figure 3
a,b) begins with
$\unicode[STIX]{x1D6FD}\approx 3$
for cellular convection, followed by a significant steepening of the exponent for the CTC regime, and finally by an abrupt slope reduction to
$\unicode[STIX]{x1D6FD}\approx 1$
. This result appears to be consistent with the findings of Julien et al. (Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016) using single-mode theory (see their figure 12). To first order, it can be seen that the topology of the curves for the CNH-QGM agree qualitatively with those obtained for DNS for
$Pr>1$
. Transitional values from cells to CTCs and CTCs to plumes occur at
$\widetilde{Ra}\approx 15$
and
$\widetilde{Ra}\approx 30$
, respectively. By contrast, the much delayed transitional values of
$\widetilde{Ra}\approx 22$
and
$\widetilde{Ra}\approx 55$
are obtained for the SF–SF case (Nieves et al.
Reference Nieves, Rubio and Julien2014). The explanation for this difference stems from the fact that CTCs form as a consequence of a thermal boundary layer instability (Julien et al.
Reference Julien, Rubio, Grooms and Knobloch2012b
), and the presence of Ekman pumping enables the threshold value for instability to be achieved at lower
$\widetilde{Ra}$
. A further inspection of figure 3(a,b) reveals that quantitative agreement between CNH-QGM and DNS is only achieved in the cellular regime for
$Pr>1$
. For larger
$\widetilde{Ra}$
, the asymptotic model yields larger values for
$Nu$
in the CTC and plume regimes before showing signs of convergence again as the GT regime is approached (figure 3
d). Estimates for exploring the geostrophic regime indicate a more computationally expensive undertaking that, while possible, is not the focus of the present work. For instance, in the stress-free case, estimated values of
$\widetilde{Ra}=220$
for
$Pr=3$
and
$\widetilde{Ra}=290$
for
$Pr=7$
are required to enter the geostrophic regime, simulations of which would demand resolutions beginning at
$500\times 500\times 500$
(see p. 4, column 2, of Julien et al.
Reference Julien, Knobloch, Rubio and Vasil2012a
).
Comparisons with DNS reveal that coherent structures, such as CTCs and plumes, have an increased vertical spatial coherence in the asymptotic model. These structures are stable conduits that result in an increased heat transport efficiency. We speculate that the difference at this lower-bound threshold value of
$E=10^{-7}$
, which is only just entering the geostrophically balanced regime, is due to subdominant ageostrophic dynamics present in the DNS but absent from the asymptotic model. Subdominant ageostrophic effects include nonlinear vertical and horizontal advection of momenta and heat. These terms destabilize the degree of vertical coherence of CTCs and permit an earlier transition to the geostrophic turbulence regime. Such effects are more prominent in the boundary layer region where velocity and thermal amplitudes are greatest. These findings are consistent with recent reports in the literature that indicate
$E=10^{-7}$
to be a threshold value for geostrophically balanced flows (Ecke & Niemela Reference Ecke and Niemela2014; Cheng et al.
Reference Cheng, Stellmach, Ribeiro, Grannan, King and Aurnou2015). Advancement will require reducing
$E$
in both laboratory experiments and DNS.
Figure 4(a) illustrates vertical root-mean-square (r.m.s.) profiles for the fluctuating temperature fields in the CTC regime at
$\widetilde{Ra}=20$
,
$\mathit{Pr}=7$
. The enhanced thermal fluctuations within the thermal wind layer are evident, with the CNH-QGM overestimating the magnitude of the DNS data, but good topological agreement is found. As a notable reminder, consistent with this observation, both the cellular and GT regimes are dominated by interior dynamics; no thermal boundary exists in the cellular regime, while the bulk in the GT regime has been shown to be the thermal bottleneck for heat transport and thus insensitive to the properties of the more thermally efficient boundary layers (Julien et al.
Reference Julien, Knobloch, Rubio and Vasil2012a
,Reference Julien, Rubio, Grooms and Knobloch
b
). Hence, the observation of increased quantitative agreement in the cellular and GT regimes implies good agreement of bulk dynamics.
R.m.s. fluctuating temperature
$\unicode[STIX]{x1D703}$
profiles at
$E=10^{-7}$
for (a) CTCs at
$\widetilde{Ra}=20$
,
$\mathit{Pr}=7$
, (b) GT at
$\widetilde{Ra}=90$
,
$\mathit{Pr}=1$
and (c) enlargement of the GT boundary layer. Results of CNH-QGM SF–SF (grey, dashed), CNH-QGM NS–NS (black, solid) and DNS NS–NS (red, dot-dashed) cases are displayed. The squares denote the locations of the maximum r.m.s.
It is well known that for
$Pr\lesssim 0.7$
linear onset is oscillatory. Fully nonlinear oscillatory regimes known to be prevalent at
$Pr\ll 1$
(e.g. King & Aurnou Reference King and Aurnou2013) are not the focus of the present study, and we consider only the range
$1/2\leqslant Pr\leqslant 1$
. For
$Pr\leqslant 1$
, we find that the cellular convection regime rapidly gives way to geostrophic turbulence. This regime is characterized by the dissipation-free scaling exponent
$\unicode[STIX]{x1D6FD}\approx 3/2$
(Julien et al.
Reference Julien, Knobloch, Rubio and Vasil2012a
). Excellent quantitative agreement is achieved between the CNH-QGM and DNS for heat transport (figure 3
c) and flow morphology (figure 4
b,c). Specifically, we observe Nusselt numbers of similar magnitudes in the geostrophic regime. The heat transport results thus establish a second marker for the fidelity of the CNH-QGM.
The upper row displays top-view visualizations of the total vertical vorticity
$\unicode[STIX]{x1D701}$
in the GT regime at
$\widetilde{Ra}=90$
,
$\mathit{Pr}=1$
for CNH-QGM: (a) SF–SF
$(E=0)$
, (b) mixed NS–SF
$(E=10^{-7})$
and (c) NS–NS
$(E=10^{-7})$
boundaries. The lower row displays the barotropic (depth-averaged) vorticity
$\langle \unicode[STIX]{x1D701}\rangle$
, illustrating the persistence of large-scale cyclonic and anticyclonic vortices in the stress-free and mixed cases. Large-scale structures form in the presence of no-slip boundaries, but they are transient and less coherent. Colour in the images is scaled between 60 (red) and
$-60$
(blue).
3.3 Non-local inverse cascade
It has recently been demonstrated that geostrophic turbulence in RRBC with stress-free boundaries yields a non-local inverse energy cascade that drives the formation of a large-scale dipole condensate consisting of a cyclonic and anticyclonic vortex pair (figure 5 a). This was established using the asymptotic NH-QGM (Julien et al. Reference Julien, Rubio, Grooms and Knobloch2012b ; Rubio et al. Reference Rubio, Julien, Knobloch and Weiss2014) and confirmed by DNS (Favier et al. Reference Favier, Silvers and Proctor2014; Guervilly et al. Reference Guervilly, Hughes and Jones2014; Stellmach et al. Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014). We note that, in the Rayleigh–Bénard (RB) configuration, the inverse cascade is unaffected by the box size, and it is found that energy always cascades to the largest scales permitted, and the associated coherent structures continually grow to box scale. However, it presently remains an open question whether (i) such a cascade persists in the presence of boundary layer friction, and (ii) if present, how its effects on the fluid dynamics may be quantified (cf. Kunnen et al. Reference Kunnen, Ostilla-Mónico, van der Poel, Verzicco and Lohse2016). Figure 6 illustrates the growth in barotropic kinetic energy occurs for all combinations of boundary conditions used in the NH-QGM. With only bulk dissipation present, slow saturation is observed for SF–SF simulations (Rubio et al. Reference Rubio, Julien, Knobloch and Weiss2014). However, in the presence of NS–NS or NS–SF boundaries, it can be seen that an inverse energy cascade still persists but is both diminished and saturated by Ekman friction (figure 6). For mixed NS–SF boundaries, the dipole condensate persists in a diminished stationary state (figure 5 b). In the presence of NS–NS boundaries, visualizations (figure 5 c) show that no persistent large-scale condensate is established. Only structures that are intermittent in time are observed (figure 7). Here, despite the perpetual flow of energy to large scales, we find that Ekman friction continually causes the large vortical structures to break up. The time series in figure 7 indicate that no large dipole forms and that the depth-independent structures vary with time in intensity and coherence.
Growth and saturation of the barotropic kinetic energy
$E_{bt}(k)$
using the CNH-QGM in the presence of Ekman friction for mixed NS–SF (solid blue) and NS–NS (solid black) boundaries at
$E=10^{-7}$
for
$\widetilde{Ra}=90$
and
$Pr=1$
. The barotropic energy for the SF–SF,
$E=0$
case (solid grey) saturates on a much longer time scale (Rubio et al.
Reference Rubio, Julien, Knobloch and Weiss2014). The total baroclinic kinetic energy
$E_{3D}(k)$
(dashed lines) saturates in all cases.
The top row shows time series from the CNH-QGM of the total vertical vorticity
$\unicode[STIX]{x1D701}$
in the GT regime at
$\widetilde{Ra}=90$
,
$\mathit{Pr}=1$
with NS–NS boundaries at
$E=10^{-7}$
. The bottom row shows the barotropic (depth-averaged) vorticity
$\langle \unicode[STIX]{x1D701}\rangle$
scaled from
$-60$
(blue) to 60 (red). The structures that form for two no-slip boundaries vary with time in both strength and coherence.
Figure 8 illustrates the 1D spectrum of horizontal kinetic energy (KE) decomposed into the barotropic and baroclinic components,
$E_{bt}(k)$
and
$E_{bc}(k)$
, as defined in (2.20). We note that the characteristics of the vertical and horizontal baroclinic spectra are similar. We also note that the ability to build an extensive inertial range by a forward cascade in rotating convection is challenging. Specifically, the baroclinic inertial range is bounded between the injection scale comparable to the linear onset scale
$L/H\sim 4.82E^{1/3}$
and the Taylor microscale. The Taylor microscale estimates the viscous distortion of eddies and, by estimation from isotropic turbulence theory, is given as
$\sqrt{10/\unicode[STIX]{x1D700}}\sim 2.71E^{1/3}$
, where
$\unicode[STIX]{x1D700}$
is defined as in (3.1). At the thermal forcing of
$\widetilde{Ra}=90$
, this implies a ratio of 1 : 1.77 (see figure 8), thus the convective forcing scale still operates very close to the diffusive scales. Furthermore, the inertial range for the inverse cascade is bounded from above by the box size
$L_{box}$
. Figure 8 therefore illustrates the direct and inverse cascades, each with a wavenumber range that is less than decadal.
The time-averaged KE spectrum for (a) mixed NS–SF and (b) NS–NS simulations with
$\mathit{Pr}=1$
,
$\widetilde{Ra}=90$
and
$E=10^{-7}$
. Wavenumbers
$k$
are normalized to the box scale
$L=10L_{c}$
, so
$k=10\tilde{k}/\tilde{k}_{c}$
, where
$\tilde{k}_{c}=1.3048$
. The black line shows the barotropic KE spectrum
$E_{bt}(k)$
with a reference slope of
$-3$
for guidance. The red line shows the horizontal baroclinic KE spectrum
$E_{bc}(k)$
taken at
$Z=1/2$
with the slope of
$-5/3$
plotted for reference. For comparison, the grey lines indicate the baroclinic (dashed) and barotropic (solid) components for the case of SF–SF boundaries. DNS NS–NS results (green symbols) are in agreement (b). The effects of Ekman friction on
$E_{bt}(k)$
and Ekman pumping on
$E_{bc}(k)$
are maximal for the NS–NS case. The spectra exhibit scaling transitions reminiscent of Nastrom & Gage (Reference Nastrom and Gage1985). The arrows mark the energy injection at box wavenumber
$k_{f}$
and the Taylor microscale at
$k_{tm}$
.
Despite having inertial ranges that are less than decadal, we find that a compelling picture emerges upon consideration of the properties of energy transfer and the instantaneous power-law exponents for the energy spectra. Hereinafter, wavenumbers
$\tilde{k}$
are normalized to the box scale
$L_{box}=\unicode[STIX]{x1D706}L_{c}$
, such that
$k=\tilde{k}/\tilde{k}_{box}=\unicode[STIX]{x1D706}\tilde{k}/\tilde{k}_{c}$
, with
$\unicode[STIX]{x1D706}=10$
or 20 (see table 1) and
$\tilde{k}_{c}=1.3048$
(Chandrasekhar Reference Chandrasekhar1961). The barotropic component dominates the energy spectra at low wavenumbers for all combinations of boundary conditions. Moreover, we find that the baroclinic kinetic energy
$E_{bc}(k)$
saturates in all cases. For SF–SF simulations (Julien et al.
Reference Julien, Rubio, Grooms and Knobloch2012b
; Rubio et al.
Reference Rubio, Julien, Knobloch and Weiss2014; Stellmach et al.
Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014), the barotropic kinetic energy
$E_{bt}(k)$
saturates on a longer time scale and exhibits the steepest power law. In the presence of Ekman friction, the barotropic kinetic energy saturates and the power-law exponent is reduced, with
$E_{bt}(k)\sim k^{-3}$
. The horizontal baroclinic spectrum
$E_{bc}(k)$
gains dominance at higher wavenumbers
$(k\approx 3)$
; it also exhibits shallower instantaneous power-law exponents
$E_{bc}(k)\sim k^{-5/3}$
. The total kinetic energy spectrum
$E_{bt}+E_{bc}$
exhibits a steep to shallow transition in the power-law exponents, a result reminiscent of the Nastrom–Gage spectrum of atmospheric and oceanic measurements with a
$-3$
to
$-5/3$
power-law transition (Nastrom & Gage Reference Nastrom and Gage1985). For both plots in figure 8, it can be seen that the spectrum is very similar in comparison with the SF–SF simulation results (displayed in grey), but with a few important distinctions due to the presence of Ekman pumping and Ekman friction. The effect of Ekman friction is evident in the decreased power in the barotropic kinetic energy at low wavenumbers, which is maximal for the NS–NS case. The effect of Ekman pumping is also evident by the increased power levels of the baroclinic kinetic energy at high wavenumbers. Excellent agreement is achieved in comparison with NS–NS DNS results in figure 8(b) (open green squares and circles). The quantitative agreement in the kinetic energy spectra thus establishes the validity of the inverse cascade in the presence of frictional boundaries and is a third marker for the fidelity of the CNH-QGM.
The transfer maps for the nonlinear barotropic self-interaction
$T_{kp}$
and the baroclinic forcing of the barotropic mode
$F_{kp}$
obtained from (2.21) are displayed in figure 9 for all combinations of boundary conditions. Each plot indicates how spectral power is transferred from wavenumber mode
$p$
to barotropic mode
$k$
via triadic interaction with catalyst wavenumber mode
$q$
such that
$\boldsymbol{k}+\boldsymbol{p}+\boldsymbol{q}=\mathbf{0}$
. Several common characteristics are evident in all plots.
The transfer maps for barotropic self-interaction
$T_{kp}$
on the top row and baroclinic-to-barotropic forcing
$F_{kp}$
along the bottom row, each averaged over several times for a case with
$\mathit{Pr}=1$
,
$\widetilde{Ra}=90$
. Plots show (a) SF–SF
$E=0$
, (b) mixed NS–SF boundaries
$E=10^{-7}$
and (c) NS–NS boundaries
$E=10^{-7}$
. Colour bars indicate how power moves from mode
$p$
to mode
$k$
. The vertical lines show the integrated values (summed over
$p$
) for each
$k$
.
The spectral power in
$T_{kp}$
is negligible in all maps for all
$(k,p)\geqslant 15$
(figure 9, upper panels). The power signature in the super- and sub-off-diagonal lines in all barotropic self-interaction
$T_{kp}$
maps indicate the existence of a forward or direct cascade, i.e. spectral power being directly transferred from low to high
$k$
wavenumbers at constant wavenumber
$p$
. For
$p\gg k$
, we also find strong evidence for the non-local inverse cascade, in that power is transferred from the high
$p$
wavenumbers directly to
$k\approx 1$
. Similarly, for
$p\ll k$
, we find energy is being extracted from high
$k$
wavenumbers and transferred to
$p\approx 1$
.
The presence of large-scale barotropic motions strongly influences the source of its origin, i.e. the baroclinic forcing
$-\langle \text{J}[\unicode[STIX]{x1D6F9}^{\prime },\unicode[STIX]{x1D701}^{\prime }]\rangle$
in (2.10). This term may be interpreted as barotropic collinearity between
$\boldsymbol{u}_{\bot }$
and
$\unicode[STIX]{x1D735}_{\bot }\unicode[STIX]{x1D701}^{\prime }$
. The transfer maps (figure 9, lower panels) clearly indicate a direct forcing of the low
$k$
wavenumber baroclinic modes by high
$p$
wavenumber modes. Particularly, one sees that the line plots for
$F(k)=\sum _{p}F_{kp}$
are positive and non-negligible for
$k\lesssim 4$
when friction is present and non-negligible when
$k=1$
for SF–SF boundaries. The increased wavenumber range in
$F(k)$
for NS–NS and NS–SF cases is consistent with the fact that it must be balanced by the energy extracted by barotropic friction
$Ef_{k}^{bt}$
, which is directly proportional to
$E_{bt}(\boldsymbol{k})$
and therefore dominates at low wavenumbers (see (2.18)). One can also see that the total nonlinear transfer
$T(k)=\sum _{p}T_{kp}$
is zero for
$k\gtrsim 8$
(line plots in figure 9) indicating the beginning of an inertial regime from
$8\lesssim k\lesssim 15$
where the observed forcing
$F_{\boldsymbol{k}}$
, dissipation
$D_{\boldsymbol{k}}$
and
$Ef_{\boldsymbol{k}}$
are negligible and the energy flux is constant. For smaller wavenumber values
$k\lesssim 8$
, the barotropic
$T(k)$
cascade first removes energy before adding energy at the smallest
$k$
. Energy in the
$k=1$
mode is maximal and minimal in the SF–SF and NS–NS cases, respectively.
Here, we remark on the efficiency by which the barotropic mode extracts energy from the baroclinic dynamics. If the forcing term
$-\langle \text{J}[\unicode[STIX]{x1D6F9}^{\prime },\unicode[STIX]{x1D701}^{\prime }]\rangle$
is removed from (2.10), thereby suppressing the barotropic cascade, the alterations to the baroclinic kinetic energy spectra are minimal. This result, first established by Rubio et al. (Reference Rubio, Julien, Knobloch and Weiss2014), indicates a symbiotic relationship between baroclinic convective motions and the large-scale barotropic dynamics, i.e. convective fluctuations adjust in the presence of barotropic dynamics without amplitude suppression. This adjustment can be contrasted to the more extreme cases of predator–prey competitive dynamics, in which relaxation oscillations result in large-amplitude variations in both baroclinic and barotropic modes.
3.4 Low-Ekman-number simulations
In this section, we extend our simulations to Ekman numbers below those presently attainable by DNS and experiments (
$E\approx 10^{-7}$
) to explore heat transport as
$E\rightarrow 0$
. We recall that the asymptotic model CNH-QGM is valid within
$\widetilde{Ra}=Ra\,E^{4/3}\lesssim E^{-1/3}$
. Moreover, in Julien et al. (Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016), it was shown that the transition from a regime adequately described by stress-free boundaries (i.e.
$E=0$
) to a regime dominated by Ekman pumping occurred at
$Ra_{t}E^{4/3}\sim E^{-1/9}$
for the special class of single-mode solutions. These results may be of particular utility in physical systems such as the Earth’s core, where
$E\approx 10^{-15}$
(Schubert & Soderlund Reference Schubert and Soderlund2011). Specifically, we consider the cases
$E=10^{-7}$
,
$10^{-9}$
and
$10^{-11}$
for Prandtl numbers 7 and 1. To account for the most restrictive resolution requirement in the vertical (i.e. the
$O(E^{1/3})$
thermal wind layer), a resolution of
$288\times 288\times 288$
with six points within the thermal wind layer was chosen for
$E=10^{-9}$
and
$384\times 384\times 576$
with five points in the thermal wind layer was chosen for
$E=10^{-11}$
.
Figure 10 shows the Nusselt number
$Nu$
as a function of
$\widetilde{Ra}$
for (a)
$Pr=7$
and (b)
$Pr=1$
. No-slip boundaries are considered with Ekman numbers
$E=10^{-7}$
,
$10^{-9}$
and
$10^{-11}$
and compared to stress-free results with
$E=0$
(grey curves). For increasing
$\widetilde{Ra}$
, it is observed that the curves associated with no-slip boundaries initially coincide with the stress-free result before departure to enhanced values of
$Nu$
. The departure locations are consistent with the
$\widetilde{Ra_{t}}=cE^{-1/9}$
prediction by Julien et al. (Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016), which we approximate simply with
$c=1$
. Specifically, with
$c=1$
, their theory predicts a departure from the stress-free curves at
$\widetilde{Ra_{t}}\approx 5.99$
for
$E=10^{-7}$
(not shown),
$\widetilde{Ra_{t}}\approx 10$
for
$E=10^{-9}$
(black arrow) and
$\widetilde{Ra_{t}}\approx 16.7$
for
$E=10^{-11}$
(red arrow). The transitional value at
$E=10^{-7}$
is below critical onset and thus departure occurs immediately.
Heat transfer
$Nu$
as a function of
$\widetilde{Ra}$
for CNH-QGM no-slip simulations with (a)
$Pr=7$
and (b)
$Pr=1$
. Comparison with the stress-free (
$E=0$
) results (grey circles) is made for Ekman numbers
$E=10^{-7}$
(blue diamonds),
$E=10^{-9}$
(black squares) and
$E=10^{-11}$
(red triangles). The arrows indicate the locations of the predicted departure at
$\widetilde{Ra_{t}}=E^{-1/9}$
(Julien et al.
Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016) away from the stress-free curve for
$E=10^{-9}$
and
$E=10^{-11}$
. The departure for
$E=10^{-7}$
is below the onset
$\widetilde{Ra}$
and is not shown.
In figure 10(a), the departures from the stress-free results are apparent for both
$E=10^{-9}$
and
$E=10^{-11}$
. The shape of these curves is largely similar to those predicted by the single-mode solutions, which provide an upper-bound estimate of the heat transfer (cf. figure 10 in Julien et al.
Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016). These results show a strong departure to higher heat transfer, and the curves for
$E=10^{-7}$
and
$E=10^{-9}$
seem to be crossing near
$\widetilde{Ra}=100$
. The results with
$Pr=1$
in figure 10(b) show a smaller departure away from the stress-free curve. We note that the vertically coherent, single-mode solutions are less relevant for low
$Pr$
. Here the flow is dominated by geostrophic turbulence, and the small departures from the stress-free case may indicate that the stress-free curve is a good leading-order approximation for the heat transfer of low-
$Pr$
, geostrophically turbulent flows.
4 Conclusions
Understanding the details and implications of rotating thermal convection in the geophysically and astrophysically relevant, but extreme, high-Rayleigh-number (
$Ra$
), low-Rossby-number (
$Ro$
) regime remains one of the outstanding challenges for both laboratory experiments and direct numerical simulations. Attainable bounds in Ekman number for an adequate survey across a sufficient range of
$Ra$
presently reside at
$E\gtrsim 10^{-7}$
(Ecke & Niemela Reference Ecke and Niemela2014; Cheng et al.
Reference Cheng, Stellmach, Ribeiro, Grannan, King and Aurnou2015) compared to planetary and stellar values of
${\lesssim}10^{-15}$
(Schubert & Soderlund Reference Schubert and Soderlund2011). The utilization of asymptotically reduced models derived from the governing Navier–Stokes equations to gain access to such regimes has a long tradition in geophysical and astrophysical fluid dynamics, beginning with the seminal introduction of the classical hydrostatic quasi-geostrophic equations of Charney (Reference Charney1948, Reference Charney1971). His model, which captures the asymptotic behaviour of geostrophically balanced large-scale eddies in a stably stratified layer while filtering fast inertio-gravity waves, remains a vital theoretical tool in investigating geophysical turbulence. Asymptotically reduced models for unstably stratified (i.e. non-hydrostatic) columnar turbulence have also played a vital role when the directions of gravity and rotation are orthogonal (Roberts Reference Roberts1968; Busse Reference Busse1970, Reference Busse1986; Calkins, Julien & Marti Reference Calkins, Julien and Marti2013). Orthogonality aids the Proudman–Taylor constraint in this configuration, which can coexist with strong thermal forcing without conflict. Only recently has this constraint of orthogonality been relaxed to produce asymptotic (non-hydrostatic quasi-geostrophic) models valid as
$E\rightarrow 0$
for stress-free bounding plates (Julien et al.
Reference Julien, Knobloch and Werne1998, Reference Julien, Knobloch, Milliff and Werne2006; Sprague et al.
Reference Sprague, Julien, Knobloch and Werne2006; Julien et al.
Reference Julien, Rubio, Grooms and Knobloch2012b
; Rubio et al.
Reference Rubio, Julien, Knobloch and Weiss2014). However, comparisons with laboratory investigations and geophysical settings, where no-slip bounding surfaces are inevitable, require an understanding of the roles of Ekman pumping and Ekman friction in the extrapolated limit
$E\rightarrow 0$
.
In this regard, the composite non-hydrostatic quasi-geostrophic model was derived by Julien et al. (Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016) to include the effect of no-slip boundaries. Through an investigation into single-mode solutions, they verified that the presence of Ekman layers formed along a no-slip boundary have a dominant positive effect on the convective heat flux and transport in the RRBC system within the geostrophic regime. This occurs, despite the increasingly narrow widths (
$\propto E^{1/2}H$
) and pumping magnitudes (
$\propto E^{1/2}\unicode[STIX]{x1D708}/H$
) as
$E\rightarrow 0$
, as a consequence of the increasing magnitude of thermal fluctuations in an
$O(E^{1/3}H)$
thermal wind boundary layer. The increased thermal fluctuations are produced by the vertical advection of the mean temperature gradient and effectively compensate the decreasing magnitude of Ekman pumping to yield
$O(1)$
convective fluxes.
In this present study, the fidelity of the reduced composite model is explored through full numerical simulations capturing all dynamically unstable modes, and results are compared against DNS. To surmise the findings, excellent qualitative and quantitative comparisons are found for three metrics: (i) flow visualizations and transitions between all geostrophic regimes, (ii) enhanced heat transport due to Ekman pumping in the Nusselt–Rayleigh relation, and (iii) the inverse barotropic cascade of kinetic energy. The latter in particular has only recently been established and quantified in the presence of stress-free boundaries (Favier et al. Reference Favier, Silvers and Proctor2014; Guervilly et al. Reference Guervilly, Hughes and Jones2014; Rubio et al. Reference Rubio, Julien, Knobloch and Weiss2014; Stellmach et al. Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014). The advantage of the model over DNS is the relaxation of the stiffness imposed by (i) the multi-time-scale properties of the thermal field, and often the explicit time representation of the Coriolis force, and (ii) spatially resolving Ekman boundary layers (see table 1). Parametrized Ekman pumping boundary layers may also be employed in the DNS; however, any gain is ultimately restricted by the explicit time-stepper constraints on the nonlinear advection terms due to the large velocity amplitudes at the computational boundary.
For
$Pr>1$
, good qualitative agreement is achieved in flow visualizations (figure 2) and good quantitative agreement in heat transport is obtained in the cellular regime and as the geostrophic turbulence regime is approached (figure 3). Plausible explanations for discrepancies in the heat transport for the intermediate convective Taylor column and plume regimes at
$E=10^{-7}$
are twofold. First, the value
$E=10^{-7}$
remains asymptotically too large, a fact that motivates future efforts in the Rayleigh–Bénard community at more extreme Ekman numbers (e.g. the Eindhoven TROCONVEX, the Gottingen U-boot and the UCLA Spinlab rotating convection experiments). Secondly, discrepancies may occur as a result of the subdominant effects of nonlinear advection of momenta and temperature omitted in the asymptotic model but inherently retained in the DNS. Resolution of each conjecture requires systematic studies of either lowering of the Ekman number in DNS or inclusion of subdominant terms in the asymptotic equations, both of which are beyond the scope and intent of the present paper.
For
$Pr\leqslant 1$
, where a direct transition from the cellular to the geostrophic turbulence regime occurs, excellent agreement with DNS (figure 3) is found. In the geostrophic turbulence regime, a non-local inverse barotropic cascade persists for all types of boundaries but is both diminished and rapidly saturated by Ekman friction. A Nastrom–Gage kinetic spectrum is found that exhibits a shallow to steep, i.e.
${\approx}-3$
to
${\approx}-5/3$
, transition in the power-law exponent. In the presence of NS–NS boundaries, large-scale vortices are no longer statistically stationary but are found to vary in both time and coherence. We note that this finding, from both the CNH-QGM and the DNS, is contrary to Kunnen et al. (Reference Kunnen, Ostilla-Mónico, van der Poel, Verzicco and Lohse2016), whose simulations argue that the inverse cascade vanishes for no-slip plates.
The investigation using low Ekman numbers (
$E\leqslant 10^{-7}$
) demonstrates the effect of decreasing
$E$
on the heat transfer enhancements. The departure from the stress-free curve remains large for simulations with
$Pr>1$
, even at
$E=10^{-11}$
. However, for
$Pr=1$
, the results suggest that the stress-free curve is a fair leading-order approximation at low
$E$
. The capacity of this model to probe smaller
$E$
provides a new avenue for explorations of RRBC, one that should be of utility to DNS and laboratory experiments, both presently limited to
$E\geqslant 10^{-7}$
, going forward. Moreover, future studies at lower
$E$
and higher
$Pr$
should provide insight into geophysical systems.
Large-scale structures are persistent features of rapidly rotating stars and planets that remain largely ill-understood. The non-local inverse cascade represents a new mechanism for generating persistent large-scale structures now captured in DNS both with (Guervilly, Hughes & Jones Reference Guervilly, Hughes and Jones2015; Yadav et al. Reference Yadav, Gastine, Christensen and Reiners2015) and without magnetic fields (Favier et al. Reference Favier, Silvers and Proctor2014; Guervilly et al. Reference Guervilly, Hughes and Jones2014; Stellmach et al. Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014). We contend that a systematic exploration of such a phenomenon will require a collective approach involving laboratory experiments, DNS and asymptotic simulations. In this regard and in a similar spirit to investigations using the classical quasi-geostrophic equations, the CNH-QGM represents a promising fundamental model that will complement emerging laboratory experiments and DNS for future investigations at more extreme Ekman numbers.
Acknowledgements
This work was supported by the National Science Foundation under EAR grants no. 1320991 and CSEDI no. 1067944 (K.J., P.M.) and by NASA Headquarters under the NASA Earth and Space Science Fellowship Program, grant 15-PLANET15F-0105 (M.P.). We wish to thank the anonymous referees for their useful comments on the manuscript. S.S. gratefully acknowledges the Gauss Centre for Supercomputing (GCS) for providing computing time through the John von Neumann Institute for Computing (NIC) on the GCS share of the supercomputer JUQUEEN at Jülich Supercomputing Centre (JSC) in Germany. This work utilized the Janus supercomputer, which is supported by the National Science Foundation (award no. CNS-0821794) and the University of Colorado Boulder.
