# The effects of Ekman pumping on quasi-geostrophic Rayleigh-Bénard convection

###### Abstract

Numerical simulations of three-dimensional, rapidly rotating Rayleigh-Bénard convection are performed using an asymptotic quasi-geostrophic model that incorporates the effects of no-slip boundaries through (i) parameterized Ekman pumping boundary conditions, and (ii) a thermal wind boundary layer that regularizes the enhanced thermal fluctuations induced by pumping. The fidelity of the model, obtained by an asymptotic reduction of the Navier-Stokes equations that implicitly enforces a pointwise geostrophic balance, is explored for the first time by comparisons of simulations against the findings of direct numerical simulations (DNS) and laboratory experiments. Results from these methods have established Ekman pumping as the mechanism responsible for significantly enhancing the vertical heat transport. This asymptotic model demonstrates excellent agreement over a range of thermal forcing for when compared with results from experiments and DNS at maximal values of their attainable rotation rates, as measured by the Ekman number (); good qualitative agreement is achieved for . Similar to studies with stress-free boundaries, four spatially distinct flow morphologies exists. Despite the presence of frictional drag at the upper and/or lower boundaries, a strong non-local inverse cascade of barotropic (i.e., depth-independent) kinetic energy persists in the final regime of geostrophic turbulence and is dominant at large scales. For mixed no-slip/stress-free and no-slip/no-slip boundaries, Ekman friction is found to attenuate the efficiency of the upscale energy transport and, unlike the case of stress-free boundaries, rapidly saturates the barotropic kinetic energy. For no-slip/no-slip boundaries, Ekman friction is strong enough to prevent the development of a coherent dipole vortex condensate. Instead vortex pairs are found to be intermittent, varying in both time and strength. For all combinations of boundary conditions, a Nastrom-Gage type spectrum of kinetic energy is found where the power law exponent changes from to , i.e., from steep to shallow, as the spectral wavenumber increases.

Key words: Bénard convection, geostrophic turbulence, quasi-geostrophic flows

## 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 2005; Gastine et al. 2014), planetary atmospheres and interiors (Aurnou et al. 2015; Heimpel et al. 2016), and deep ocean convection on Earth (Marshall & Schott 1999). 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 . Flows are characterized by the nondimensional Rayleigh (), convective Rossby () and Ekman () numbers defined as

(1.0) |

These, respectively, measure the strength of thermal forcing, the strength of system rotation relative to thermal forcing, and the importance of viscous diffusion relative to . Here denotes gravity, the layer depth, and the thermal expansion coefficient. The Prandtl number, Pr = , 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 and low Rossby numbers , implying small Ekman numbers . For example, estimates for these parameters in Earth’s outer core are , , and (Schubert & Soderlund 2011). 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 (Pope 2000) or in the case of the outer core. The low 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 timescales 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 DNS and the thermal and mechanical capabilities of laboratory experiments.

While the extreme parameter values present a challenge for DNS (Favier et al. 2014; Guervilly et al. 2014; Stellmach et al. 2014; Horn & Shishkina 2015), they allow for the use of perturbation methods that rely on the small size of as an asymptotic expansion parameter. The reduced Non-Hydrostatic Quasi-Geostrophic Model (NH-QGM) (Julien et al. 1998, 2006; Sprague et al. 2006) was developed as a simplification, or reduction, of the Boussinesq equations at low . The model implicitly enforces geostrophy as a pointwise constraint and is asymptotically exact for RRBC in the presence of stress-free boundary conditions as . Of particular interest is the dependence of the heat transport, as measured by the nondimensional Nusselt number , on the external control parameters and the molecular fluid property ; here is the heat flux and is the volumetric heat capacity. Importantly, 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 (Stellmach et al. 2014). Figure 1a illustrates the high level of agreement in the - relation for stress-free boundaries with the asymptotic scaling of .

For the case of no-slip boundaries (see figure 1b) substantial disagreement is found in comparison with the NH-QGM results for the - relation. DNS and laboratory data in water () show markedly increased values of for all values of with a steep scaling exponent of (Cheng et al. 2015; Stellmach et al. 2014). 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 , whereas experiments are performed in cylindrical containers with ) More specifically, it is shown that the divergence in between the NH-QGM, where Ekman boundary layers are presumed passive, and DNS data grows as the Ekman number is lowered in the latter. A recent study by Cheng et al. (2015) (see their figure 6) indicates that the heat transport exponent in the - relation increases monotonically from to as decreases to . By replacing the no-slip boundaries with parameterized Ekman pumping boundary conditions and establishing quantitative agreement, Stellmach et al. (2014) 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 and , have diminishing values as (The latter formulae for represents an exact parameterization of the laminar Ekman boundary layer (Greenspan 1968)). Indeed, the asymptotic NH-QGM assumes that Ekman boundary layers become increasingly insignificant as and are filtered from the dynamics, thus implicitly assuming impenetrable stress-free boundaries (Julien et al. 1998; Sprague et al. 2006).

An asymptotic analysis of the boundary layer regions of RRBC using the Boussinesq equations was recently reported in Julien et al. (2016). Specifically, it was established that the effects of Ekman pumping become significant within the rotationally constrained regime when , 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 parameterized 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 . Thus, the model, called the Composite NH-QGM (CNH-QGM), now includes the Ekman number as an explicit control parameter. The capabilities of the composite model were explored in Julien et al. (2016) 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 influence the path of the realized solution in phase space and therefore capture and characterize the qualitative behavior of the enhanced heat transport in the fully nonlinear regime.

In this paper, we extend the work of Julien et al. (2016) 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 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 Section 2, the CNH-QGM derived in Julien et al. (2016) is documented and the multi-layer structure unique to the presence of no-slip boundaries identified. In Section 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 Section 4.

## 2 Composite Model

The full derivation and discussion of the CNH-QGM is presented in Julien et al. (2016); here we simply provide an overview of their asymptotic model. The leading order velocity field outside the Ekman layer, , is in geostrophic balance where . It follows that the horizontal velocity field solved by is non-divergent where the pressure becomes the geostrophic streamfunction. On defining , the reduced CNH-QGM system governing the evolution of the fluid is given by

(2.0) |

(2.0) |

(2.0) |

(2.0) |

capturing, respectively, the evolution of vertical vorticity , vertical velocity , and temperature at reduced Rayleigh number for a given Prandtl number . The temperature is decomposed into a mean (horizontally-averaged) component evolving on a slow timescale and a small fluctuating component characteristic of RRBC. The overbar indicates the horizontal and fast time averaging operator given by

(2.0) |

where is the area of the box. Here denotes the material derivative with horizontal geostrophic advection and in (2) denotes the ageostrophic velocity field determined through the three-dimensional incompressibility condition

(2.0) |

The system is solved for the RRBC configuration with fixed temperature boundaries

(2.0) |

Velocity boundary conditions are either parameterized Ekman pumping for no-slip (NS)

(2.0) |

or implicit impenetrable stress-free boundaries ( SF)

(2.0) |

with or mixed (NS-SF) combinations thereof (e.g., NS at and SF at ). Equations (2)-(2) with one of the above choices of velocity boundary conditions form a closed reduced system that filters fast inertial waves occurring on timescales and spatial scales. The imposition of strong spatial anisotropy by nondimensionalizing vertical and horizontal scales with and is a key step in the filtering that gives rise to the reduced dynamics (Julien & Knobloch 2007). Velocity and time are nondimensionalized using and and temperature is nondimensionalized by .

The evolution of the barotropic vertical vorticity is a key quantity in the understanding of the inverse energy cascade recently observed in reduced simulations of RRBC (Rubio et al. 2014) as well as DNS (Favier et al. 2014; Guervilly et al. 2014; Stellmach et al. 2014). Depth averaging the vertical vorticity equation (2) according to , with the vortical variables decomposed as , gives

(2.0) |

indicating that the growth of the barotropic vertical vorticity is dependent on the net balance between right-hand-side terms respectively for the barotropic self-interaction, baroclinic (i.e. convective) forcing, barotropic and baroclinic Ekman friction, and barotropic viscous dissipation. Here denotes the Jacobian for the horizontal advection of , when is either or .

We spectrally analyze equation (2) where barotropic and baroclinic fields are decomposed according to the discrete Fourier transform

(2.0) |

with . Specifically, we consider the evolution of the spectral barotropic kinetic energy at each wave number

(2.0) |

obtained by multiplication of the Fourier decomposition of (2) by , followed by integration over physical space. Here and describe how power is converted to barotropic Fourier amplitude through the barotropic and baroclinic triadic interactions involving modes and such that (Rubio et al. 2014; Favier et al. 2014), where

Box Dimensions | CNH-QGM Resolution | CNH-QGM | DNS Resolution | CNH-QGM | DNS | CNH-QGM | DNS | ||
---|---|---|---|---|---|---|---|---|---|

x x | N x N x N | x/L | N x N x N | ||||||

15 | 1/2 | 10 x 10 x 1 | 192 x 192 x 192 | 0.868 | 384 x 384 x 385 | 2x | 1.66x | 4.445 | 4.526 |

15 | 1 | 10 x 10 x 1 | 192 x 192 x 192 | 0.633 | 288 x 288 x 385 | 2x | 1.64x | 4.398 | 4.414 |

15 | 3 | 10 x 10 x 1 | 192 x 192 x 192 | 0.811 | 384 x 384 x 385 | 2x | 1.86x | 5.020 | 4.897 |

15 | 7 | 10 x 10 x 1* | 192 x 192 x 192 | 0.547 | 288 x 288 x 385 | 2x | 1.71x | 6.162 | 5.825 |

20 | 1/2 | 10 x 10 x 1 | 288 x 288 x 192 | 1.04 | 384 x 384 x 385 | 1x | 9.78x | 8.623 | 8.491 |

20 | 1 | 10 x 10 x 1 | 192 x 192 x 192 | 0.903 | 384 x 384 x 385 | 3x | 6.04x | 10.18 | 10.11 |

20 | 3 | 10 x 10 x 1 | 192 x 192 x 192 | 1.22 | 384 x 384 x 385 | 3x | 5.54x | 22.91 | 14.83 |

20 | 7 | 10 x 10 x 1* | 192 x 192 x 192 | 0.826 | 384 x 384 x 385 | 2x | 4.82x | 28.92 | 16.87 |

30 | 1/2 | 10 x 10 x 1 | 384 x 384 x 192 | 0.729 | 384 x 384 x 385 | 5x | 1.07x | 14.60 | 15.00 |

30 | 1 | 10 x 10 x 1 | 384 x 384 x 192 | 0.601 | 288 x 288 x 385 | 1x | 5.01x | 20.44 | 19.61 |

30 | 3 | 10 x 10 x 1 | 288 x 288 x 192 | 1.09 | 384 x 384 x 385 | 1x | 2.14x | 50.56 | 37.29 |

30 | 7 | 10 x 10 x 1* | 288 x 288 x 216 | 0.850 | 288 x 288 x 385 | 2x | 2.36x | 112.2 | 72.86 |

50 | 1/2 | 10 x 10 x 1 | 384 x 384 x 288 | 1.00 | 384 x 384 x 385 | 3x | 5.74x | 29.82 | 30.73 |

50 | 1 | 10 x 10 x 1 | 384 x 384 x 288 | 0.739 | 384 x 384 x 385 | 7x | 2.87x | 31.02 | 30.15 |

50 | 3 | 10 x 10 x 1 | 288 x 288 x 288 | 1.39 | 384 x 384 x 385 | 1x | 1.51x | 74.70 | 62.87 |

50 | 7 | 10 x 10 x 1* | 462 x 462 x 384 | 0.436 | 384 x 384 x 385 | 1x | 1.10x | 157.2 | 89.17 |

70 | 1/2 | 10 x 10 x 1 | 520 x 520 x 384 | 0.898 | 768 x 768 x 769 | 1x | 1.78x | 46.16 | 48.37 |

70 | 1 | 10 x 10 x 1 | 384 x 384 x 384 | 0.687 | 384 x 384 x 385 | 1x | 2.12x | 42.96 | 42.63 |

70 | 3 | 10 x 10 x 1 | 288 x 288 x 192 | 1.61 | 384 x 384 x 385 | 1x | 1.40x | 84.21 | 72.89 |

70 | 7 | 10 x 10 x 1* | 384 x 384 x 384 | 0.923 | 384 x 384 x 385 | 1x | 9.04x | 179.4 | 122.0 |

90 | 1/2 | 10 x 10 x 1 | 520 x 520 x 384 | 1.15 | 768 x 768 x 1025 | 1x | 1.12x | 68.34 | 67.36 |

90 | 1 | 10 x 10 x 1 | 384 x 384 x 384 | 1.02 | 576 x 576 x 513 | 1x | 1.11x | 59.22 | 58.04 |

100 | 3 | 10 x 10 x 1 | 288 x 288 x 192 | 1.76 | 384 x 384 x 385 | 1x | 1.28x | 95.71 | 84.81 |

100 | 7 | 10 x 10 x 1* | 384 x 384 x 384 | 1.03 | 384 x 384 x 513 | 5x | 7.47x | 182.7 | 146.4 |

(2.0) | |||||

(2.0) | |||||

(2.0) |

and denotes the real part and is the Kronecker delta function. The latter terms in (2) are defined as

(2.0) | |||||

(2.0) | |||||

(2.0) |

where

(2.0) |

Here the asterisk denotes the conjugate. We find that baroclinic Ekman friction is negligible compared to barotropic Ekman friction, i.e., .

From (2), one can consider the evolution of the 1D energy spectra, defined as

(2.0) |

where, in the polar representation, and . Also of particular interest are the transfer maps

(2.0) | |||

(2.0) |

detailing the transfer of energy from all wavenumbers of magnitude to magnitude .

## 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 and Prandtl numbers . 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 2009), and a mixed semi-implicit explicit (IMEX) timestepping scheme. The DNS code utilizes the second-order, IMEX Adams Bashforth/backward-differentiation scheme (AB2/BDF2) (Stellmach & Hansen 2008) and the CNH-QGM employs the third-order IMEX Runge-Kutta acheme (RK3) (Spalart et al. 1991). 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 thermal wind layer for the CNH-QGM and the Ekman layer in the presence of no-slip boundaries for the DNS. Specifically, the resolutions selected ensure a minimum of grid points in the thermal wind layer and a minimum of grid points in the Ekman boundary layer, which is calculated using the distance to the first maxima in the vertical spiral profile (Stellmach et al. 2014). The resolution in the horizontal direction is controlled by the viscous dissipation, which is estimated by the Kolmogorov length scale

(3.0) |

The energy dissipation rate 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 (see column 5, Table 1). This criteria was adequate for convergence in the Nusselt number. Moreover, under this criteria, simulations with sufficiently large box sizes (i.e., greater than 10 x 10 x 1) resulted in the same values as those reported in the table.

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 CFL constraint imposed by the explicit treatment of the nonlinear advection terms. In the case of DNS, geostrophy and the separation of timescales 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 magnitude from the timestep (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.

Due to spatiotemporal constraints, DNS with still remains the lower bound for investigations of rotating convection (Stellmach et al. 2014). Similar bounds are also present in laboratory experiments because of thermal and mechanical constraints (Ecke & Niemela 2014; Cheng et al. 2015). The CNH-QGM is capable of explorations at much lower as a result of the asymptotic design. Specifically, we note that represents the case of purely stress-free boundaries considered elsewhere (Sprague et al. 2006; Julien et al. 2012b, a; Rubio et al. 2014). 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 . In subsection 3.4, results are reported from simulations with the Ekman number varying below .

### 3.1 Flow Morphologies

The four convective flow morphologies and their domains of existence as a function of and has been well documented using the NH-QGM for RRBC in the presence of top and bottom stress-free (SF-SF) boundaries (c.f. figure 2 in Julien et al. (2012b)). Figure 2 illustrates that simulations of the CNH-QGM using the parameterized no-slip (NS-NS) conditions yield the same four geostrophic regimes. Ordered by increasing 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 and yield comparable visualization results (c.f. figure 2 in Stellmach et al. (2014)), thereby establishing a first marker for the fidelity of the CNH-QGM.

### 3.2 Heat Transport

A comparison of the - heat transport relations at for DNS and CNH-QGM is shown in figure 3 for (plot a), (plot b) and (plot c). Comparisons with are included in each plot for reference. Confirming the findings of Stellmach et al. (2014) and Cheng et al. (2015) (see figure 1), Ekman pumping results in a significantly enhanced for all compared to simulations with the SF-SF boundaries. Moreover, as with the stress-free simulations (Julien et al. 2012b; Nieves et al. 2014), transitions between regimes are observable in figure 3 by abrupt changes in the local heat transport exponent (reference slopes are provided to guide the eye).

For , the heat transport law (figure 3a,b) begins with for cellular convection, followed by a significant steepening of the exponent for the CTC regime, and finally by an abrupt slope reduction to . This result appears to be consistent with the findings of Julien et al. (2016) 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 . Transitional values from cells to CTCs and CTCs to plumes occur at and respectively. By contrast the much delayed transitional values of and are obtained for the SF-SF case (Nieves et al. 2014). Explanation for this difference stems from the fact that CTCs form as consequence of a thermal boundary layer instability (Julien et al. 2012b), and the presence of Ekman pumping enables the threshold value for instability to be achieved at lower . A further inspection of figure 3a,b reveals that quantitative agreement between CNH-QGM and DNS is only achieved in the cellular regime for . For larger the asymptotic model yields larger values for in the CTC and plume regimes before showing signs of convergence again as the GT regime is approached (figure 3d). 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 for and for are required to enter the geostrophic regime, simulations of which would demand resolutions beginning at 500 x 500 x 500. (page 4, column 2 of Julien et al. 2012a).

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 , 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 to be a threshold value for geostrophically balanced flows (Ecke & Niemela 2014; Cheng et al. 2015). Advancement will require reducing in both laboratory experiments and DNS.

Figure 4a illustrates vertical RMS profiles for the fluctuating temperature fields in the CTC regime at , . 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. 2012a, b). Hence, the observation of increased quantitative agreement in the cellular and GT regimes implies good agreement of bulk dynamics.

It is well-known that for linear onset is oscillatory. Fully nonlinear oscillatory regimes known to be prevalent at (e.g., King & Aurnou 2013) are not the focus of the present study, and we consider only the range . For , we find that the cellular convection regime rapidly gives way to geostrophic turbulence. This regime is characterized by the dissipation free scaling exponent (Julien et al. 2012a). Excellent quantitative agreement is achieved between the CNH-QGM and DNS for heat transport (figure 3c) and flow morphology (figure 4b,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.

### 3.3 Non-local inverse cascade

It has recently been demonstrated that geostrophic turbulence in RRBC with stress-free boundaries yields a nonlocal inverse energy cascade that drives the formation of large scale dipole condensate consisting of a cyclonic and anticyclonic vortex pair (figure 5a). This was established using the asymptotic NH-QGM (Julien et al. 2012b; Rubio et al. 2014) and confirmed by DNS (Favier et al. 2014; Guervilly et al. 2014; Stellmach et al. 2014). We note that in the Rayleigh-Bénard 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. (2016)). 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. 2014). 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 5b). In the presence of NS-NS boundaries, visualizations (figure 5c) 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.

Figure 8 illustrates the 1D spectra of horizontal kinetic energy decomposed into the barotropic and baroclinic components, and , as defined in (2). 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 and the Taylor micro scale. The Taylor micro scale estimates the viscous distortion of eddies and, by estimation from isotropic turbulence theory, is given as , where is defined as in (3). At the thermal forcing of , this implies a ratio of : (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 . Figure 8 therefore illustrates the direct and inverse cascades, each with a wavenumber range that is less than decadal.

Despite having inertial ranges that are less than decadal, we find a compelling picture emerges upon consideration of the properties of energy transfer and the instantaneous power law exponents for the energy spectra. Hereinafter, wavenumbers are normalized to the box scale , such that with or (see Table 1) and (Chandrasekhar 1961). The barotropic component dominates the energy spectra at low wavenumbers for all combinations of boundary conditions. Moreover, we find that the baroclinic kinetic energy saturates in all cases. For SF-SF simulations (Julien et al. 2012b; Rubio et al. 2014; Stellmach et al. 2014), the barotropic kinetic energy 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 . The horizontal baroclinic spectra gains dominance at higher wavenumbers , it also exhibit shallower instantaneous power law exponents . The total kinetic energy spectra exhibits a steep to shallow transition in the power law exponents, a result reminiscent of the Nastrom-Gage spectra of atmospheric and oceanic measurements with a to power law transition (Nastrom & Gage 1985). For both plots in figure 8, it can be seen that the spectra is very similar in comparison with the SF-SF simulation results (displayed in gray) 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 8b (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 and the baroclinic forcing of the barotropic mode obtained from (2) are displayed in figure 9 for all combinations of boundary conditions. Each plot indicates how spectral power is transferred from wavenumber mode to barotropic mode via triadic interaction with catalyst wavenumber mode such that . Several common characteristics are evident in all plots.

The spectral power in is negligible in all maps for all (Figure 9, upper panels). The power signature in the super and sub off-diagonal lines in all barotropic self-interaction maps indicate the existence of a forward or direct cascade, i.e., spectral power being directly transferred from low to high -wavenumbers at constant wavenumber . For , we also find strong evidence for the nonlocal inverse cascade in that power is transferred from the high wavenumbers directly to . Similarly, for , we find energy is being extracted from high wavenumbers and transferred to .

The presence of large scale barotopic motions strongly influences the source of its origin, i.e. the baroclinic forcing in (2). This term may be interpreted as barotropic collinearity between and . The transfer maps (Figure 9, lower panels) clearly indicate a direct forcing of the low wavenumber baroclinic modes by high wavenumber modes. Particularly, one sees that the line plots for are positive and non-negligible for when friction is present and non-negligible when for SF-SF boundaries. The increased wavenumber range in for NS-NS and NS-SF cases is consistent with the fact that it must be balanced by the energy extracted by barotropic friction which is directly proportional to and therefore dominates at low wavenumbers (see (2)). One can also see that the total nonlinear transfer is zero for (line plots in Figure 9) indicating the beginning of an inertial regime from where the observed forcing , dissipation , and are negligible and the energy flux is constant. For smaller wavenumber values the barotropic cascade first removes energy before adding energy at the smallest . Energy in the 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 baroclininc dynamics. If the forcing term is removed from (2), thereby suppressing the barotropic cascade, the alterations to the baroclinic kinetic energy spectra are minimal. This result, first established by Rubio et al. (2014), 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 Simulations

In this section, we extend our simulations to Ekman numbers below those presently attainable by DNS and experiments () to explore the heat transport as . We recall that the asymptotic model CNH-QGM is valid within . Moreover in Julien et al. (2016), it was shown that the transition from a regime adequately described by stress-free boundaries (i.e., ) to a regime dominated by Ekman pumping occurred at for the special class of single-mode solutions. These results may be of particular utility in physical systems such as Earth’s core where (Schubert & Soderlund 2011). Specifically, we consider the cases , and for Prandtl numbers and . To account for the most restrictive resolution requirement in the vertical (i.e., the thermal wind layer), a resolution of x x with 6 points within the thermal wind layer was chosen for and x x with 5 points in the thermal wind layer was chosen for .

Figure 10 shows the Nusselt number as a function of for (a) and (b) . No-slip boundaries are considered with Ekman numbers and and compared to stress-free results with (gray curves). For increasing , it is observed that the curves associated with no-slip boundaries initially coincide with the stress-free result before departure to enhanced values of . The departure locations are consistent with the prediction by Julien et al. (2016), which we approximate simply with . Specifically, with , their theory predicts a departure from the stress-free curves at for (not shown), for (black arrow), and for (red arrow). The transitional value at is below critical onset and thus departure occurs immediately.

In Figure 10a, the departures from the stress-free results are apparent for both and . 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. 2016). These results show a strong departure to higher heat transfer, and the curves for and seem to be crossing near . The results with in Figure 10b show a smaller departure away from the stress-free curve. We note that the vertically coherent, single mode solutions are less relevant for low . Here the flow is dominated by geostrophic turbulence and the small departures from the stress-free may indicate that the stress-free curve is a good leading-order approximation for the heat transfer of low , geostrophically turbulent flows.

## 4 Conclusions

Understanding the details and implications of rotating thermal convection in the geophysically and astrophysically relevant, but extreme, high Rayleigh - low Rossby number 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 presently reside at (Ecke & Niemela 2014; Cheng et al. 2015) compared to planetary and stellar values of (Schubert & Soderlund 2011). 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 quasigeostrophic equations of Charney (1948, 1971). His model, which captures the asymptotic behavior 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. nonhydrostatic) columnar turbulence have also played a vital role when the directions of gravity and rotation are orthogonal (Roberts 1968; Busse 1970, 1986; Calkins et al. 2013). Orthogonality aids the Proudman-Taylor constraint in this configuration, which can co-exist with strong thermal forcing without conflict. Only recently has this constraint of orthogonality been relaxed to produce asymptotic (nonhydrostatic quasigeostrophic) models valid as for stress-free bounding plates (Julien et al. 1998, 2006; Sprague et al. 2006; Julien et al. 2012b; Rubio et al. 2014). 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 .

In this regard, the composite nonhydrostatic quasigeostrophic model was derived by Julien et al. (2016) 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 () and pumping magnitudes () as , as a consequence of the increasing magnitude of thermal fluctuations in an 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 (Rubio et al. 2014; Guervilly et al. 2014; Favier et al. 2014; Stellmach et al. 2014). 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 the Ekman boundary layers (see Table 1). Parameterized 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 non-linear advection terms due to the large velocity amplitudes at the computational boundary.

For , 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 are twofold. First, the value 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 asymptotical model but inherently retained in the DNS. Resolution to 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 manuscript.

For , 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 nonlocal 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. to , 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 both in time and coherence. We note that this finding, from both the CNH-QGM and the DNS, is contrary to Kunnen et al. (2016), whose simulations argue that the inverse cascade vanishes for no-slip plates.

The investigation using low Ekman numbers () demonstrates the effect of decreasing on the heat transfer enhancements. The departure from the stress-free curve remains large for simulations with , even at . However, for , the results suggest that the stress-free curve is a fair leading-order approximation at low . The capacity of this model to probe smaller provides a new avenue for explorations of RRBC, one that should be of utility to DNS and laboratory experiments, both presently limited to , going forward. Moreover, future studies at lower and higher 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 et al. 2015; Yadav et al. 2015) and without magnetic fields (Favier et al. 2014; Guervilly et al. 2014; Stellmach et al. 2014). 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 quasigeostrophic equations, the CNH-QGE represents a promising fundamental model that will complement emerging laboratory experiments and DNS for future investigations at more extreme Ekman numbers.

## Acknowledgments

This work was supported by the National Science Foundation under EAR grants #1320991 & CSEDI #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 number CNS-0821794) and the University of Colorado Boulder.

- Aurnou et al. (2015) Aurnou, J. M., Calkins, M. A., Cheng, J. S., Julien, K., King, E. M., Nieves, D., Soderlund, K. M. & Stellmach, S. 2015 Rotating convective turbulence in earth and planetary cores. Physics of the Earth and Planetary Interiors 246, 52–71.
- Busse (1970) Busse, F. H. 1970 Thermal instabilities in rapidly rotating systems. J. Fluid Mech. 44, 441–460.
- Busse (1986) Busse, F. H. 1986 Asymptotic theory of convection in a rotating, cylindrical annulus. J. Fluid Mech. 173, 545–556.
- Calkins et al. (2013) Calkins, M. A., Julien, K. & Marti, P. 2013 Three-dimensional quasi-geostrophic convection in the rotating cylindrical annulus with steeply sloping endwalls. Journal of Fluid Mechanics 732, 214–244.
- Chandrasekhar (1961) Chandrasekhar, S. 1961 Hydrodynamic and Hydromagnetic Stability. Oxford: Oxford University Press.
- Charney (1971) Charney, J. 1971 Geostrophic turbulence. Journal of Atmoshperic Science 28, 1087–1095.
- Charney (1948) Charney, J. G. 1948 On the scale of atmospheric motions. Geofys. Publ. 17, 3–17.
- Cheng et al. (2015) Cheng, J. S., Stellmach, S., Ribeiro, A., Grannan, A., King., E. M. & Aurnou, J. M. 2015 Laboratory-numerical models of rapidly rotating convection in planetary cores. Geophysics Journal International 201, 1–17.
- Ecke & Niemela (2014) Ecke, R. E. & Niemela, J. J. 2014 Heat transport in the geostrophic regime of rotating Rayleigh-Bénard convection. Phys. Rev. Lett. 113, 114301.
- Favier et al. (2014) Favier, B., Silvers, L. J. & Proctor, M. R. E. 2014 Inverse cascade and symmetry breaking in rapidly rotating Boussinesq convection. Phys. Fluids 26 (096605).
- Gastine et al. (2014) Gastine, T., Heimpel, M. & Wicht, J. 2014 Zonal flow scaling in rapidly-rotating compressible convection. Physics of the Earth and Planetary Interiors 232, 36–50.
- Greenspan (1968) Greenspan, H. P. 1968 The Theory of Rotating Fluids. Cambridge, United Kingdom: Cambridge University Press.
- Guervilly et al. (2014) Guervilly, C., Hughes, D. W. & Jones, C. A. 2014 Large-scale vortices in rapidly rotating Rayleigh–Bénard convection. Journal of Fluid Mechanics 758, 407–435.
- Guervilly et al. (2015) Guervilly, C., Hughes, D. W. & Jones, C. A. 2015 Generation of magnetic fields by large-scale vortices in rotating convection. Physical Review E 91 (4), 041001.
- Heimpel et al. (2016) Heimpel, M., Gastine, T. & Wicht, J. 2016 Simulation of deep-seated zonal jets and shallow vortices in gas giant atmospheres. Nature Geoscience 9, 19â–23.
- Horn & Shishkina (2015) Horn, S. & Shishkina, O. 2015 Toroidal and poloidal energy in rotating Rayleigh-Bénard convection. J. Fluid Mech. 762, 232–255.
- Julien et al. (2016) Julien, K., Aurnou, J., Calkins, M., Knobloch, E., Marti, P., Stellmach, S. & Vasil, G. 2016 A nonlinear model for rotationally constrained convection with ekman pumping. Journal of Fluid Mechanics 798, 50–87.
- Julien & Knobloch (2007) Julien, K. & Knobloch, E. 2007 Reduced models for fluid flows with strong constraints. Journal of Mathematical Physics 48 (6).
- Julien et al. (2006) Julien, K., Knobloch, E., Milliff, R. & Werne, J. 2006 Generalized quasi-geostrophy for spatially anistropic rotationally constrained flows. J. Fluid Mech. 555, 233–274.
- Julien et al. (2012a) Julien, K., Knobloch, E., Rubio, A. M. & Vasil, G. M. 2012a Heat transport in Low-Rossby-number Rayleigh-Bénard Convection. Phys. Rev. Lett. 109 (254503).
- Julien et al. (1998) Julien, K., Knobloch, E. & Werne, J. 1998 A new class of equations for rotationally constrained flows. Theoret. Comput. Fluid Dyn. 11, 251–261.
- Julien et al. (2012b) Julien, K., Rubio, A. M., Grooms, I. & Knobloch, E. 2012b Statistical and physical balances in low Rossby number Rayleigh-Bénard convection. Geophysical & Astrophysical Fluid Dynamics 106 (4-5), 392–428.
- Julien & Watson (2009) Julien, K. & Watson, M. 2009 Efficient multi-dimensional solution of PDEs using Chebyshev spectral methods. Journal of Computational Physics 228, 1480–1503.
- King & Aurnou (2013) King, E. M. & Aurnou, J. M. 2013 Turbulent convection in liquid metal with and without rotation. Proceedings of the National Academy of Sciences 110 (17), 6688–6693, arXiv: http://www.pnas.org/content/110/17/6688.full.pdf.
- Kunnen et al. (2016) Kunnen, R. P. J., Ostilla-MÃ³nico, R., van der Poel, E. P., Verzicco, R. & Lohse, D. 2016 Transition to geostrophic convection: the role of the boundary conditions. Journal of Fluid Mechanics 799, 413–432.
- Marshall & Schott (1999) Marshall, J. & Schott, F. 1999 Open-ocean convection: Observations, theory, and models. Reviews of Geophysics 37 (1), 1–64.
- Miesch (2005) Miesch, M. S. 2005 Large-scale dynamics of the convection zone and tachocline. Living Reviews in Solar Physics 2 (1).
- Nastrom & Gage (1985) Nastrom, G. D. & Gage, K. S. 1985 A climatology of atmospheric wavenumber spectra of wind and temperature observed by commercial aircraft. Journal of the atmospheric sciences 42 (9), 950–960.
- Nieves et al. (2014) Nieves, D., Rubio, A. M. & Julien, K. 2014 Statistical classification of flow morphology in rapidly rotating Rayleigh-Bénard convection. Physics of Fluids 26 (8).
- Pope (2000) Pope, S. B. 2000 Turbulent Flows. Cambridge: Cambridge University Press.
- Roberts (1968) Roberts, P. H. 1968 On the thermal instability of a rotating-fluid sphere containing heat sources. Phil. Trans. R. Soc. A 263, 93–117.
- Rubio et al. (2014) Rubio, A. M., Julien, K., Knobloch, E. & Weiss, J. B. 2014 Upscale energy transfer in three-dimensional rapidly rotating turbulent convection. Phys. Rev. Lett. 112, 144501.
- Schubert & Soderlund (2011) Schubert, G. & Soderlund, K. M. 2011 Planetary magnetic fields: Observations and models. Physics of the Earth and Planetary Interiors 187 (3), 92–108.
- Spalart et al. (1991) Spalart, P. R., Moser, R. D. & Rogers, M. M. 1991 Spectral methods for the Navier-Stokes equations with one infinite and two periodic directions. J. Comp. Phys. 96, 297–324.
- Sprague et al. (2006) Sprague, M., Julien, K., Knobloch, E. & Werne, J. 2006 Numerical simulation of an asymptotically reduced system for rotationally constrained convection. Journal of Fluid Mechanics 551, 141–174.
- Stellmach & Hansen (2008) Stellmach, S. & Hansen, U. 2008 An efficient spectral method for the simulation of dynamos in Cartesian geometry and its implementation on massively parallel computers. Geochem. Geophys. Geosys. 9 (5).
- Stellmach et al. (2014) Stellmach, S., Lischper, M., Julien, K., Vasil, G., Cheng, J. S., Ribeiro, A., King, E. M. & Aurnou, J. M. 2014 Approaching the asymptotic regime of rapidly rotating convection: Boundary layers versus interior dynamics. Phys. Rev. Lett. 113, 254501.
- Yadav et al. (2015) Yadav, R. K., Gastine, T., Christensen, U. R. & Reiners, A. 2015 Formation of starspots in self-consistent global dynamo models: Polar spots on cool stars. Astronomy & Astrophysics 573, A68.