# [

###### Abstract

It is a well established result of linear theory that the influence of differing mechanical boundary conditions, i.e., stress-free or no-slip, on the primary instability in rotating convection becomes asymptotically small in the limit of rapid rotation (Chandrasekhar 1961). This is accounted for by the diminishing impact of the viscous stresses exerted within Ekman boundary layers and the associated vertical momentum transport by Ekman pumping (Niiler & Bisshopp 1965; Heard & Veronis 1971). By contrast, in the nonlinear regime recent laboratory experiments and supporting numerical simulations are now providing evidence that the efficiency of heat transport remains strongly influenced by Ekman pumping in the rapidly rotating limit (Stellmach et al. 2014; Cheng et al. 2015). In this paper, a reduced model is developed for the case of low Rossby number convection in a plane layer geometry with no-slip upper and lower boundaries held at fixed temperatures. A complete description of the dynamics requires the existence of three distinct regions within the fluid layer: a geostrophically balanced interior where fluid motions are predominately aligned with the axis of rotation, Ekman boundary layers immediately adjacent to the bounding plates, and thermal wind layers driven by Ekman pumping in between. The reduced model uses a classical Ekman pumping parameterization to alleviate the need for spatially resolving the Ekman boundary layers. Results are presented for both linear stability theory and a special class of nonlinear solutions described by a single horizontal spatial wavenumber. It is shown that Ekman pumping (which correlates positively with interior convection) allows for significant enhancement in the heat transport relative to that observed in simulations with stress-free boundaries. Without the intermediate thermal wind layer the nonlinear feedback from Ekman pumping would be able to generate a heat transport that diverges to infinity. This layer arrests this blowup resulting in finite heat transport at a significantly enhanced value. With increasing buoyancy forcing the heat transport transitions to a more efficient regime, a transition that is always achieved within the regime of asymptotic validity of the theory, suggesting this behavior may be prevalent in geophysical and astrophysical settings. As the rotation rate increases the slope of the heat transport curve below this transition steepens, a result that is in agreement with observations from laboratory experiments and direct numerical simulations.

Rotating convection with Ekman pumping ] A nonlinear model for rotationally constrained convection with Ekman pumping

K. Julien et al.]
Keith Julien^{†}^{†}thanks: Email address for correspondence: julien@colorado.edu,Jonathan M. Aurnou,
Michael A. Calkins,

Edgar Knobloch,
Philippe Marti,
Stephan Stellmach and

Geoffrey M. Vasil

## 1 Introduction

Rotating Rayleigh-Bénard convection (RBC), i.e., a rotating horizontal fluid layer heated from below and cooled from above, provides a canonical framework for the study of fluid phenomena influenced by rotation and thermal forcing. It has proven to be an indispensable framework for understanding fluid motions in geophysical and astrophysical systems including planetary atmospheres and interiors (Aurnou et al. 2015), rapidly rotating stars (Miesch 2005), and open ocean deep convection (Marshall & Schott 1999). In many of these examples the dominant influence of rotation results in geostrophy where the Coriolis force is balanced by the pressure gradient force. Small departures from this dominant force balance, referred to as quasigeostrophy, are known to be capable of producing turbulent fluid motions characterized by anisotropic eddies elongated along the rotation axis (Sakai 1997). The effect of spatial anisotropy on fluid mixing and global transport properties such as heat and energy transport remains an important and largely unanswered question.

In the context of the Boussinesq approximation for incompressible motions the rotating RBC problem is completely specified via three nondimensional parameters, namely, the Rayleigh number , the Ekman number and the Prandtl number , defined by

(1.0) |

Here is the layer depth, is the destabilizing temperature difference, is the rotation rate of the system, is the kinematic viscosity, is the thermal diffusivity, is the gravitational acceleration and is the coefficient of thermal expansion. The Rayleigh number measures the magnitude of the thermal forcing and the Ekman number measures the importance of viscous forces relative to the Coriolis force. The Prandtl number is the ratio of the thermal and viscous diffusion timescales in the system and describes the thermophysical properties of the working fluid. Another dimensionless parameter of importance is the convective Rossby number that measures the relative importance of thermal forcing and the Coriolis force. Specifically, rotationally constrained convection is characterized by . Importantly, it is known that rotation imparts rigidity to the fluid in the RBC problem that delays the onset of convection until a critical Rayleigh number (Chandrasekhar 1961) is reached with associated . At onset, motions are columnar with a horizontal scale . By definition of the convective Rossby number, the rotationally constrained branch is defined as and characterized by , or equivalently, where is the reduced Rayleigh number (Julien et al. 2012a). The range of permissible can thus be vast covering as much as five decades in geophysical and astrophysical settings where . Rich dynamics are observed along this rotationally constrained branch, ranging from coherent laminar to highly turbulent states (Sprague et al. 2006; Julien et al. 2012b; Rubio et al. 2014).

The efficiency of heat transport, as measured by the nondimensional Nusselt number , is perhaps the most common result reported in the literature given that the functional dependence is tied to the underlying dynamics. Here is the heat flux and is the volumetric heat capacity. Attempts have been made to characterize the rotationally constrained regime by a heat transport scaling law of the form with (Rossby 1969; Liu & Ecke 1997; King et al. 2009; Ecke & Niemela 2014). At fixed , combined investigations of laboratory experiments and DNS with no-slip boundaries have reported values in the range with a trend towards the upper bound occurring at the lowest Ekman number (Figure 1) (Cheng et al. 2015). At sufficiently large a transition to the weakly rotating or non-rotating regime occurs, characterized by , where (see Figure 1). Comparative studies in the presence of stress-free boundaries using direct numerical simulations (DNS) have reported similar findings to those seen in Figure 1 but with substantially smaller exponents at (King et al. 2009; Stellmach et al. 2014). Irrespective of the boundary conditions, King et al. (2009) have established that the transition to the nonrotating scaling law occurs at smaller and smaller as the rotation rate increases, i.e., where . This result is attributed to the loss of geostrophic balance in the thermal boundary layer where the local Rossby number reaches unity while the Rossby number in the bulk remains small (Julien et al. 2012a).

The difference in the heat transport scaling exponents and has been associated with Ekman pumping (Stellmach et al. 2014) – the vertical momentum transport that results from the transition from an interior geostrophic balance to a dominant boundary layer force balance between the Coriolis and viscous force in the presence of no-slip boundaries. In the linear regime, Ekman pumping promotes the destabilization of the fluid layer, an effect quantified by the positive difference - in the critical Rayleigh numbers (Niiler & Bisshopp 1965; Heard & Veronis 1971). This difference is asymptotically small in the limit indicating that differing mechanical boundaries become asymptotically indistinguishable. By contrast, laboratory experiments and supporting DNS indicate that the discrepancy due to pumping in the fully nonlinear regime appears to remain finite as (Stellmach et al. 2014; Cheng et al. 2015). Unfortunately, surveying the high –low regime remains a prohibitive challenge for both laboratory experiments and DNS. Laboratory experiments are constrained by their inability to access the rotationally constrained branch at sufficiently low owing to the decreasing accuracy of global heat transport measurements resulting from unknown heat leaks (King et al. 2009; Ecke & Niemela 2014). DNS studies are restricted by spatiotemporal resolution constraints imposed by the presence of Ekman boundary layers and fast inertial waves propagating on an timescale (Nieves et al. 2014; Stellmach et al. 2014).

Although an impediment to DNS, the stiff character of the governing fluid equations in the limit provides a possible path forward for simplifying, or reducing, the governing equations. Indeed, a system of reduced NonHydrostatic Quasigeostrophic Equations (NH-QGE) appropriate for rotating RBC in the presence of stress-free boundaries have been successfully derived and utilized by Julien and collaborators (Julien et al. 1998; Sprague et al. 2006; Julien et al. 2006; Julien & Knobloch 2007; Grooms et al. 2010; Julien et al. 2012a, b; Rubio et al. 2014; Nieves et al. 2014). The NH-QGE, which filter fast inertial waves while retaining inertial waves propagating on slow advective timescales, enable parameter space explorations in the high –low limit. The NH-QGE have been used to identify various flow regimes and heat transfer scaling behavior as a function of , thus providing a valuable roadmap for both DNS and laboratory experiments. Importantly, comparisons with DNS for stress-free boundaries (Stellmach et al. 2014) have established good quantitative agreement in both heat transport and flow morphology at . For and large this approach has established an ultimate exponent , a result that corresponds to a dissipation-free scaling law in the rotationally constrained turbulence regime (Julien et al. 2012a, b).

In the present work we extend the asymptotic theory resulting in the NH-QGE to the case of no-slip boundary conditions in which Ekman boundary layers are present. This is the situation pertinent to the laboratory and of relevance to geophysical scenarios such as convection in the Earth’s liquid iron outer core which is bounded from below by the solid iron inner core and above by a rocky mantle. In agreement with the linear investigation of Heard & Veronis (1971), our nonlinear analysis shows that the presence of no-slip boundaries requires the existence of three distinct regions each characterized by a different dominant physical balance: and hereafter denoted as the outer , middle and inner regions (see Figure 2).

The outer region corresponds to the fluid interior (i.e., the bulk region of depth ). Within this region fluid motions are horizontally non-divergent and in pointwise geostrophic balance. The dynamics are asymptotically described by the NH-QGE investigated by Julien & collaborators (1998, 2006, 2007, 2012a,b, 2014). Similar to the classical quasigeostrophic equations (Charney 1948, 1971; Eady 1949; Pedlosky 1987; Vallis 2006) toroidal or vortical fluid motions aligned with the axis of rotation are forced solely by the vortex stretching associated with the linear Coriolis force. The NH-QGE differ from the classical quasigeostrophic equations in that vertical motions required by vortex stretching are now comparable in magnitude to horizontal motions and thus the nonhydrostatic inertial acceleration force must be retained in the vertical momentum balance. This results in a dynamic, as opposed to a diagnostic, evolution of the vertical velocity field. In the presence of stress-free boundaries the NH-QGE are sufficient to completely describe all the ‘slow’ dynamics occurring within the fluid layer.

The inner regions are the Ekman boundary layers immediately adjacent to the horizontal boundaries. These layers are required to attenuate the interior geostrophic velocity fields to zero. Within them, the geostrophic balance of the interior is relaxed and as a consequence fluid motions become horizontally divergent with cross-isobaric flow. Mass conservation requires that vertical motions are induced, a process referred to as Ekman pumping. The dynamics within the Ekman layers are described by a classical set of reduced linear equations (Greenspan 1969), where as a consequence of the spatial anisotropy and the vortical magnitude observed in rotating RBC (Chandrasekhar 1961), vertical pumping velocities of magnitude are found (Niiler & Bisshopp 1965; Heard & Veronis 1971). As done in large-scale oceanic and atmospheric applications, this linear property of the dynamics can be successfully utilized and the effects of this layer on the fluid interior can be parameterized by the application of simple pumping boundary conditions (Pedlosky 1987). This parameterized approach has been applied successfully in numerous previous works ranging from the investigation of Stewartson layer instabilities (Schaeffer & Cardin 2005), spherical convection (Aubert et al. 2003; Calkins et al. 2012), rotating RBC (Stellmach et al. 2014), and is also employed in the present work.

As shown by Heard & Veronis (1971) Ekman pumping induces thermal fluctuations that cannot be regulated within the Ekman boundary layer to satisfy the thermal boundary conditions. Middle layers of depth arise that separate the Ekman boundary layers from the geostrophically balanced fluid interior. Importantly, the requirement that thermal fluctuations vanish at the fixed temperature bounding plates necessitates the introduction of vertical diffusion for temperature fluctuations. This is an insignificant process in the interior in the rapidly rotating limit (Sprague et al. 2006; Julien et al. 2012a). However, as shown in this paper, the middle region is characterized by a thermal wind balance, i.e., a geostrophic and hydrostatic diagnostic balance, and the dynamics of this region evolve according to an advection-diffusion equation for the temperature fluctuations forced by Ekman pumping. Specifically, the magnitude of the vertical advection of the mean temperature field within the middle layers increases with Rayleigh number until it becomes a dominant source of buoyancy production that is comparable to that produced in the geostrophic interior. This enhancement is due to the intensification of vortical motions, and hence vertical transport, in the vicinity of the boundaries. Convective fluxes driven by Ekman transport are also comparable to that produced in the geostrophic interior. Therefore, changes to the heat transport are to be expected, as observed in laboratory experiments and simulations (Stellmach et al. 2014; Cheng et al. 2015).

We show that the asymptotic results for the different regions may be combined into a single composite system of reduced equations, which we refer to as the CNH-QGE (Composite NonHydrostatic QuasiGeostrophic Equations). On the rotationally constrained branch of RBC, the CNH-QGE are valid in the interval . The linear stability properties of these equations are shown to be consistent with earlier work (Niiler & Bisshopp 1965; Heard & Veronis 1971). Moreover, comparisons with the NH-QGE for single-mode (or single horizontal wavenumber) solutions show that significant departures in heat transport occur in the interval .

The remainder of the paper is organized as follows. In section 2, we present the formulation of the rotating RBC problem in terms of the incompressible Navier-Stokes equations. In section 3, we present the asymptotic development in the presence of no-slip bounding plates using the ‘Method of Composite Expansions’ (see section 4.2, Nayfeh (2008)). In particular, the Ekman and thermal wind boundary layers are identified and analyzed. The significance of Ekman pumping is also deduced by determining the Rayleigh number threshold at which the resulting vertical pumping velocity induces order one changes in the heat transport. It is also established that such a transition always occurs on the rotationally constrained branch of RBC. In section 4, we present the composite reduced model, i.e., the CNH-QGE where the dynamics of each layer are combined into a unified description (see Equations (4.0)-(4.0)). Results from the model are presented in section 5, and establish that this unified model captures the physics associated with no-slip boundary conditions observed in laboratory experiments and DNS. Concluding remarks are given in section 6.

## 2 Governing Equations

We consider thermally driven fluid flows that are characterized by the dimensional scales of length , velocity , time , pressure and destabilizing temperature jump . Assuming a Cartesian coordinate system we adopt the Rayleigh-Bénard configuration of a plane-parallel geometry rotating about the -axis with constant angular velocity in the presence of constant gravity . The nondimensional equations of motion are given by the Boussinesq equations

(2.0) |

(2.0) | |||||

(2.0) |

for the velocity field , temperature and the pressure , where the material derivative . The nondimensional parameters that appear are defined as

(2.0) |

respectively denoting the Rossby, Euler, buoyancy, Reynolds and Péclet numbers. In the present work we are interested in the rotationally constrained regime characterized by and aspect ratio for columnar structures of depth (Julien et al. 1998, 2006; Sprague et al. 2006; Julien et al. 2012b).

For fluid motions in a statistically stationary state, the nondimensional vertical heat transport, i.e., the Nusselt number , is given by

(2.0) |

obtained upon averaging equation (2) over time and the horizontal cross-section. Here , where is the horizontal cross-sectional area. This result indicates that the heat flux through the layer is constant at every vertical level.

## 3 Asymptotic Development

It is well established that in the geophysically and astrophysically relevant regimes
the presence of fast inertial waves, propagating on timescales, poses a severe restriction on DNS.
^{}^{}The discretized equations resulting from the Boussinesq equations will, in general, be coupled through the Coriolis force
. This coupling is routinely eliminated by an explicit treatment in many timestepping algorithms,
i.e., by its relegation and evaluation at previous steps. This favorable numerical situation occurs at the expense of imposing severe timestepping restrictions. Implicit treatment circumvents this issue. However, prohibitive timestepping restrictions persist owing to the Ekman-dependent CFL time constraint
associated with advective nonlinearites.
The evolution of turbulent eddies is insensitive to these waves which can be filtered from the governing equations by asymptotic reduction methods in much the same manner as done in atmospheric and oceanic sciences for stably-stratified layers. Indeed, Julien and collaborators (Julien et al. 2006; Sprague et al. 2006; Julien & Knobloch 2007) have established that an asymptotic reduction of the governing equations (2)-(2) can be deduced upon using as a small parameter and introducing the distinguished limits

(3.0) |

For the appropriate choice of the horizontal diffusive velocity scale, , it follows that

(3.0) |

such that . Here corresponds to the reduced Rayleigh number involving the Prandtl number .

In accord with Figure 2, we anticipate the existence of three distinct regions: the bulk and middle regions, and the Ekman layers, each with their respective nondimensional depths and . We therefore employ a multiple scales expansion in the vertical direction and time,

(3.0) |

where the slow vertical coordinate of the bulk is defined by , the fast coordinate of the Ekman layer defined by , and the slow time by . We find that this setup necessitates the decomposition of the fluid variables into mean (horizontally averaged) and fluctuating (horizontally varying) components, respectively denoted by overbars and primes, e.g.,

(3.0) |

Averaging over the fast time is also required,

(3.0) |

We note a posteriori that, unlike derivations of the reduced dynamics in the presence of stress-free boundaries (Sprague et al. 2006), the consideration of no-slip boundaries and Ekman layers requires the separation of spatial and time-averaging operations as performed here.

To proceed, all fluid variables are decomposed into outer , middle and inner components representing the fluid bulk, middle regions, and Ekman layers. For example,

Here, capitalizations are used to identify the individual contributions of the fluid variables to each region and the notation denotes the upper and lower boundaries, respectively. The boundary layer coordinates are assumed to increase away from the physical boundaries at Each region of the fluid layer may be accessed by the following actions for the outer, middle, and inner limits

Similar expressions hold for the upper inner and middle layers located at upon replacing with . By definition, the middle variables are identically zero in the outer region, while the inner variables are identically zero in both the middle and outer regions. Contributions to regions or involving outer variables, indicated in (3a) and (3a), are obtained by Taylor-expanding and then taking the appropriate limit. Hereinafter, for notational convenience, contractions and are used when referring to both upper and lower regions. Furthermore, contractions such as omit reference to the dependence on other variables (i.e., ) and refer to the vertical coordinate.

Asymptotic series in powers of are now introduced for all fluid variables and substituted into the governing equations (2)-(2). An order by order analysis is then performed. The asymptotic procedure using decompositions of the form (3) with (3)-(3) is referred to as the ‘Method of Composite Expansions’ (see section 4.2, Nayfeh (2008)). Specifically, rather than using the ‘Method of Matched Asymptotic Expansions’ – i.e., first determining the inner and outer expansions analytically, matching them, and then forming the composite expansion (Van Dyke 1975) – here variables of the form (3) are automatically valid everywhere provided they satisfy the physical boundary conditions at the bounding plates.

### 3.1 The Outer Region: Nonhydrostatic Quasigeostrophic Equations

Within the fluid interior, inner and middle variables are identically zero. We introduce expansions in powers of of the form

(3.0) |

where (Sprague et al. 2006; Julien et al. 2006; Julien & Knobloch 2007). The leading order mean component satisfies the motionless hydrostatic balance

(3.0) |

together with . The leading order convective dynamics are found to be incompressible and geostrophically balanced, i.e.,

(3.0) |

where denotes the geostrophic operator. By definition, all outer variables are independent of
and so .^{}^{}
If the small-scale dependence were retained, geostrophy would automatically imply the Proudman-Taylor (PT) constraint
on the small vertical scale (Proudman 1916; Taylor 1923).
The diagnostic balance given by (3.0) is solved by

(3.0) |

for the streamfunction and vertical velocity . Here we adopt the definition with , so that the pressure is now identified as the geostrophic streamfunction. It follows that leading order motions are horizontally nondivergent with . Three-dimensional incompressiblity is captured at the next order, , where

(3.0) |

This results in the production of subdominant ageostrophic motions driven by vertical gradients in . The prognostic evolution of these variables is obtained from balances at next order,

obtained by projecting onto the null space of (Sprague et al. 2006; Calkins et al. 2013). This amounts to performing the projections and on (3.1). On noting that this procedure results in asymptotically reduced equations for the vertical vorticity , vertical velocity , and thermal anomaly :

(3.0) | |||||

(3.0) | |||||

(3.0) |

The evolution of the mean temperature field is deduced at upon averaging over the fast scales :

(3.0) |

In a statistically stationary state, it follows that

(3.0) |

Equations (3.0)-(3.0), (3.0) and (3.0) constitute the asymptotically reduced system referred to as the NonHydrostatic QuasiGeostrophic Equations (NH-QGE). A notable feature in the NH-QGE is the absence of (higher order) vertical advection. This is a hallmark characteristic of quasigeostrophic theory. Equation (3.0) states that vertical vorticity, or toroidal motions, are affected by horizontal advection, vortex stretching arising from the linear Coriolis force, and horizontal diffusion, while equation (3.0) shows that vertical motions are affected by horizontal advection, unbalanced pressure gradient, horizontal diffusion and buoyancy. The buoyancy forces are captured by the fluctuating and mean temperature equations (3.0, 3.0).

The system (3.0)-(3.0) is accompanied by impenetrable boundary conditions

(3.0) |

together with thermal conditions, hereafter taken to be the fixed temperature conditions

(3.0) |

In the presence of impenetrable boundaries, the boundary limits or of equation (3.0) for the temperature fluctuations reduce to the advection-diffusion equation

(3.0) |

The horizontally-averaged variance for such an equation evolves according to

(3.0) |

Therefore decreases monotonically to zero in time implying the implicit thermal boundary condition . Together with the impenetrability condition (3.0), the vertical momentum equation (3.0) implies that motions along the horizontal boundaries are implicitly stress-free with

(3.0) |

If rapidly rotating RBC in the presence of stress-free boundary conditions is the primary objective, a well-posed closed system is obtained from (3.0)-(3.0), (3.0) and (3.0) together with impenetrable boundary conditions (3.0) and fixed mean temperature boundary conditions (3.0).

#### 3.1.1 Validity of the NH-QGE

In the presence of stress-free boundary conditions the NH-QGE remain valid throughout the entire flow domain provided geostrophy, Eq. (3.0), holds. This remains so provided the local Rossby number . Given , one finds^{}^{}The Landau notation little- denotes a function that is of lower order of magnitude than a given function, that is, the function is of a lower order than the function .

(3.0) |

In a detailed investigation of the NH-QGE, Julien et al. (2012a) have shown that this criterion is violated at the transitional value

(3.0) |

In this regime the thermal boundary layers experience a loss of geostrophic balance. This provides an upper bound for comparisons of NH-QGE with DNS with stress-free boundary conditions. Indeed, as illustrated in Figure 3(a) for , within the regime of validity, , good quantitative agreement in the heat transport measurements is observed (Stellmach et al. 2014). Simulations (Julien et al. 2012b; Nieves et al. 2014) of the NH-QGE prior to this transition have revealed four different flow morphologies subsequently confirmed by both DNS (Stellmach et al. 2014) and laboratory experiments (Cheng et al. 2015) as is increased (see Figure 4): a cellular regime, a convective Taylor column (CTC) regime consisting of weakly interacting shielded columns, a plume regime where CTCs have lost stability, and finally a geostrophic turbulence regime that is also associated with an inverse energy cascade that produces a depth-independent large-scale dipole vortex pair. All regimes are identifiable by changes in the heat transport exponent: the CTC regime exhibits a steep heat transport scaling law where , whereas the geostrophic turbulence regime is characterized by a dissipation-free scaling law (Julien et al. 2012a) and an inverse turbulent energy cascade (Julien et al. 2012a, b; Rubio et al. 2014; Favier et al. 2014; Guervilly et al. 2014; Stellmach et al. 2014). A rigorous upper bound heat transport result for the NH-QGE, , where is a constant, has also been reported (Grooms 2015; Grooms & Whitehead 2015).

On the other hand, comparison between the NH-QGE and DNS with no-slip boundaries and laboratory experiments (Stellmach et al. 2014) reveals substantial differences (Figure 3(b)). Specifically, a steep scaling law in the heat transport is observed in the DNS study. Moreover, DNS and laboratory results both suggest that the steep scaling continues as (Figure 1). It is now evident from the implicitly enforced stress-free boundary condition (3.0) that the reduced NH-QGE system and its solutions cannot be uniformly continued to impenetrable no-slip boundaries, where

(3.0) |

In the presence of no-slip boundaries, it is well-known that the viscous boundary layers are Ekman layers of depth (Greenspan 1969). Within these layers the geostrophic velocity field in the bulk must be reduced to zero. In the following, we proceed with an analysis of the rotationally constrained regime with the intent of extending the NH-QGE to the case of no-slip boundaries.

### 3.2 Inner Region: Ekman Boundary Layers

To avoid duplication, we focus on the lower Ekman boundary layer at with nondimensional depth (an identical analysis applies for the upper boundary layer at ). This depth arises as a result of the spatially anisotropic structure of rapidly rotating convection (Heard & Veronis 1971). In dimensional terms the Ekman layer depth since .

The equations that capture the Ekman layer dynamics are obtained by taking the inner limit of the governing equations (2)-(2) about using (3). We pose an asymptotic inner expansion in powers of of the form

(3.0) |

and utilize, a posteriori, knowledge that the middle layer variables have the asymptotic form

(3.0) |

and therefore do not contribute to leading order. Subtracting the contributions of the outer region then yields

(3.0) | |||||

(3.0) | |||||

(3.0) |

Given that the mean components are identically zero, the primed notation is omitted. It follows from equation (3.0) that the pressure within the Ekman boundary layer is the same as that outside it and we thus take . The classical linear Ekman equations are therefore

(3.0) | |||||

(3.0) |

After a simple reformulation and the introduction of no-slip boundary conditions, we have

(3.0) |

where we have utilized in advance that the leading order middle layer variables . Since the flow within the Ekman layer is horizontally divergent, equation (3.0) implies the presence of vertical motions with velocity .

The classical solution (Greenspan 1969) is found within the Ekman layer at and is given by

(3.0) | |||||

(3.0) | |||||

(3.0) |

Figure 5 illustrates sample boundary layer profiles obtained from DNS. The projected horizontal velocities are in perfect agreement with (3.0) and (3.0). These structures are robust throughout the boundary layer, thus providing confirmation of the existence of a linear Ekman layer (Stellmach et al. 2014). Application of mass conservation (3.0) yields the inner and outer vertical velocities

As , we see that and hence that

(3.0) |

This relation, often referred to as the Ekman pumping boundary condition, constitutes an exact parameterization of the linear Ekman layer in rotating RBC and represents Ekman pumping when and Ekman suction when . Hereafter, we do not distinguish between these two cases and refer to this phenomenon generically as Ekman pumping. Use of the Ekman pumping boundary condition alleviates the need to resolve the velocity dynamics that occur within the Ekman layer. A similar analysis of the Ekman layer at gives

(3.0) |

Equations (3.0)-(3.0) exemplify the following important distinction between the asymptotic procedure performed here and the linear analyses of Niiler & Bisshopp (1965) and Heard & Veronis (1971) which implicitly assume that as opposed to (parameterized boundary conditions were not uncovered in these articles). The latter approach yields a perturbative analysis that captures a range of for which the effect of Ekman pumping remains asymptotically close to the stress-free problem. However, for a sufficiently large , determined below in §3.3.1, the asymptotic expansion breaks down and becomes nonuniform. For completeness, this result is summarized in the Appendix. We find that maintaining asymptotic uniformity requires that pumping be elevated to contribute to the leading order vertical velocity, as pursued here. This approach enables a complete exploration of Ekman pumping for all for which the NH-QGE remain valid.

Several studies (Barcilon 1965; Faller & Kaylor 1966; Dudis & Davis 1971) have established that solutions to the classical linear Ekman layer are unstable to small horizontal scale disturbances that evolve on rapid timescales that are filtered from the reduced dynamics. This occurs when the boundary layer Reynolds number . Although this is within the realm of possibility for the rotationally constrained regime according to (3.0), to alter the pumping parameterization such an instability must also reach amplitudes comparable to the bulk vorticity. No evidence of this is presently seen in DNS.

The Ekman layer solutions (3.0)-(3.0) together with the application of the exact parameterizations (3.0)-(3.0) indicate that the outer velocity fields can be continued uniformly to the computational boundaries at . Evidence for this two-layer structure is provided in Figure 6 which illustrates the RMS velocity profiles at various magnifications obtained from DNS at . Specifically, the figure demonstrates that the velocity structure is unaffected by the presence of the thermal boundary layer, defined in terms the maxima of the RMS of (purple line, plot (b)).

These results suggests that Ekman layers in DNS of rotating RBC may be replaced by the parameterized pumping boundary conditions

(3.0) |

Stellmach et al. (2014) have demonstrated the accuracy of this boundary layer parameterization via comparison of the heat transport obtained from DNS with no-slip boundaries and DNS where the Ekman pumping conditions are used. With all other details being identical in the two DNS studies, excellent quantitative agreement is reached. This result establishes that it is the presence of Ekman layers that is responsible for the strong differences in heat transport observed between stress-free and no-slip boundaries (Figure 3).

The inner component of the pumping velocity defined in (3.2) gives rise to an inner temperature fluctuation satisfying

(3.0) |

Utilizing (3.2), the solution to this equation is given by

(3.0) |

The limiting values as a function of and are

(3.0) |

It thus follows that temperature fluctuations within the Ekman layer are of magnitude . This observation yields an estimate of the convective heat transport,

(3.0) |

which is smaller in magnitude by a factor of when compared to that occurring in the convective interior, namely . We can thus conclude that the observed enhancement of heat flux as measured by (Figure 3) cannot occur directly within the Ekman layer given the vorticity bound (3.0).

Inspection of the RMS thermal profiles at increasing magnification obtained at shows no visible boundary layer structure in the vicinity of the Ekman layer (cf. Figure 7, plots (b) and (c)). The RMS profiles reveal that the thermal boundary layer extends much farther into the interior than the Ekman layer observed in Figure 6. Indeed, we recall that the mean temperature is an outer variable independent of the fast spatial variables and (Sprague et al. 2006).

#### 3.2.1 The Significance of Ekman Pumping

It now remains to determine the nature of the thermal response in the immediate vicinity of the Ekman layer. Focusing again on the lower boundary, this requires consideration of how the reduced outer equation for the temperature fluctuations (3.0) can be continued to the physical boundaries. On Taylor-expanding all fluid variables within the Ekman layer, the outer component of the pumping velocity induces outer temperature fluctuations that satisfy

(3.0) |

Here, we have set .

In a statistically stationary state, averaging the equation for the thermal variance obtained from (3.0) gives the balance

(3.0) |

Clearly, in contrast to stress-free boundaries where , pumping induces the enhanced thermal response

(3.0) |

immediately outside the Ekman layer. The associated enhancement in the convective flux is estimated as