3.1. Region-I (Interior)

Given the scaling of (2.4)–(2.7), we can read off the leading-order interior dynamical equations by inspection. That is, for the bulk interior, all variables retain their original amplitude scalings, and the Coriolis and hydrostatic terms are the only remaining in (2.4)–(2.7) as ${\textit{E}} \to 0$ . The bulk dynamics then only contains ${\mathcal{O}}( \varepsilon ^{0} )$ amplitudes and the only time derivative occurs in the temperature equation (2.8). It is possible to derive this explicitly by considering the leading-order balances of the full asymptotic expansion. From the momentum equations and continuity,

(3.13) \begin{eqnarray} U_{0} = -\partial _{Y}P_{0}, \quad V_{0} = \partial _{X} P_{0}, \quad W_{0} = 0, \quad \Theta _{0} = \partial _{Z} P_{0}. \end{eqnarray}

The pressure gives the horizontal velocities and temperature in a thermal-wind balance (i.e. a joint hydrostatic and geostrophic balance).

The temperature equation retains its original nonlinear form (except for vertical advection),

(3.14) \begin{eqnarray} (\partial _{t} + U_{0} \partial _{X} + V_{0}\partial _{Y}) \, \Theta _{0} = \left( \partial _{X}^{2} + \partial _{Y}^{2} + \partial _{Z}^{2} \right)\, \Theta _{0}. \end{eqnarray}

The combined (3.13) and (3.14) almost form a closed system, with two important caveats.

The first major caveat is that the temperature evolution equation contains no obvious sources. Defining the horizontal velocity vector, $U_{\bot } = (U_{0}, V_{0})$ , the available thermal energy evolution,

(3.15) \begin{align} \frac {\textrm d}{{\textrm d}t} \int _{0}^{1} \!\int _{\mathcal{A}} \frac {\Theta _{0}^{2}}{2} \,\text{d}{X} \,\text{d}{Y} \,\text{d}{Z} & + \int _{0}^{1} \!\int _{\partial \mathcal{A}} \left [ \frac {\Theta _{0}^{2}}{2} U_{\bot } - \Theta _{0} \boldsymbol{\nabla} _{\bot } \Theta _{0} \right ] \boldsymbol{\cdot} \hat {n} \,\text{d}{\ell } \,\text{d}{Z}\nonumber \\[5pt]& =-\int _{0}^{1} \!\int _{\mathcal{A}} \left [ |\boldsymbol{\nabla} _{\bot }\Theta _{0}|^{2} + |\partial _{Z}\Theta _{0}|^{2}\right ] \text{d}{X} \,\text{d}{Y} \,\text{d}{Z} \leqslant 0, \end{align}

where $\,\text{d}{\ell }$ is the integration measure along the boundary, $\partial {\mathcal{A}}$ , and $\hat {n}$ is the outward unit normal vector. Equation (3.15) shows that the boundary terms must force the system strongly enough to overcome persistent dissipation. Without the boundary terms (like in conventional bulk convection), the available thermal energy would monotonically decrease. Driving (3.14) requires forcing from the sidewall boundary layers, which we explain in § 3.2.

The second caveat is that (3.14) determines the evolution only for the temperature. Obtaining the velocities requires the pressure and obtaining the pressure requires solving hydrostatic balance (3.13). The difficulty comes from not knowing the depth-independent, or barotropic, component of the pressure,

(3.16) \begin{eqnarray} \langle P_{0} \rangle \equiv \int _{0}^{1}P_{0} \, \text{d}Z. \end{eqnarray}

We derive an equation for $\!\left \langle P_{0} \right \rangle$ by considering the top and bottom boundaries together with their velocity boundary conditions in § 3.3.

3.2. Region-II (thermal Stewartson layer)

In the Stewartson sidewall layers, focusing on $X=\Gamma$ and $- \infty \lt x \le 0$ , the physical container boundary corresponds to $x=0$ and $X=\Gamma$ , while the boundary of the interior corresponds to $X=\Gamma$ and $x \to - \infty$ . The leading-order momentum equation (2.4) in the $X$ -direction gives

(3.17) \begin{eqnarray} && u_{-1} = u_{-1/2} = 0. \end{eqnarray}

The leading-order diffusion of temperature from (2.8) gives

(3.18) \begin{eqnarray} \partial _{x}^{2}\Theta _{0} = 0, \end{eqnarray}

which implies a constant leading-order temperature across the boundary layer. Strictly speaking, (3.18) contains a linearly growing solution in $x$ . The possibility of unbounded amplitudes as $|x| \to \infty$ preclude this solution. Therefore, the leading-order temperature exists as a pure interior dynamical variable,

(3.19) \begin{eqnarray} \Theta _{0} = \Theta _{0}(X,Y,Z). \end{eqnarray}

While (3.19) implies that $\Theta _{0}$ does not vary across the boundary layer, there exist temperature fluctuations in the directions along the wall both in $Y$ and $Z$ . The interior thermally forces the boundary layer via tangential gradients.

Given a temperature field in the interior, the leading-order momentum and continuity equations follow:

(3.20) \begin{align} \partial _{x}p_{0} - v_{-1} &= 0, \end{align}

(3.21) \begin{align} \partial _{Y} p_{0} + u_{0} - \partial _{x}^{2}v_{-1} &= 0, \end{align}

(3.22) \begin{align} \partial _{Z} p_{0} - \partial _{x}^{2} w_{-1} &= \Theta _{0}, \end{align}

(3.23) \begin{align} \partial _{x}u_{0} + \partial _{Y}v_{-1} + \partial _{Z}w_{-1} &= 0. \end{align}

The system is similar to the Stewartson layer balances derived in other work (Stewartson Reference Stewartson1957; van Heijst Reference van Heijst1983; Kunnen et al. Reference Kunnen, Clercx and van Heijst2013). The main difference is that (3.20)–(3.23) have enough derivatives to satisfy all velocity boundary conditions at the physical wall and thermal anomalies force the whole boundary layer from the interior. There is no need for the classical ${\textit{E}}^{\, 1/4}$ exterior layer typically required to satisfy all velocity boundary conditions.

The interior forces the vertical velocity via (3.22). Leading-order balances also imply $\partial _{x}^{2} \Theta _{1/2} = 0$ . To see how the boundary layer feeds back to the interior, the next-order temperature equation for $\Theta _{1}$ ,

(3.24) \begin{eqnarray} v_{-1} \, \partial _{Y} \Theta _{0} + w_{-1} \, ( \partial _{Z}\Theta _{0} - R) = \partial _{x}^{2} \Theta _{1}. \end{eqnarray}

The crucial connection between the interior and boundary layer comes from the full form of the insulating thermal boundary condition at the physical sidewall,

(3.25) \begin{eqnarray} \partial _{X}\Theta _{0} + \partial _{x}\Theta _{1} = 0 \quad \text{at} \quad x=0,\ X=\Gamma . \end{eqnarray}

This couples $\Theta _{0}$ , from the interior, with $\Theta _{1}$ , from (3.24). Integrating (3.24) over boundary layer coordinate,

(3.26) \begin{eqnarray} x\in (-\infty , 0\, ] = {\mathbb{R}}^{-}, \end{eqnarray}

gives

(3.27) \begin{eqnarray} q_{Y} \, \partial _{Y} \Theta _{0} + q_{Z} \, (\partial _{Z}\Theta _{0} - R) = \partial _{x}\Theta _{1}\big |_{x=0} = -\partial _{X}\Theta _{0}\big |_{X=\Gamma } ,\end{eqnarray}

where the $\partial _{x}\Theta _{1}\big |_{x\to -\infty }$ term vanishes because $v_{-1}$ and $w_{-1}$ decay exponentially away from $x=0$ (see later), and there is no linearly growing homogeneous solution in $x$ with bounded amplitude as $x\to -\infty$ . Equation (3.27) expresses energy conservation between the boundary layer and the interior.

The total momentum fluxes within the sidewall layers are

(3.28) \begin{eqnarray} q_{Y} = \int _{\mathbb{R}^{-}} \! \! v_{-1} \,\text{d}{x}, \quad q_{Z} = \int _{\mathbb{R}^{-}} \! \! w_{-1} \,\text{d}{x}. \end{eqnarray}

These purely surface quantities depend on the tangent directions, $(Y, Z)$ . They couple to the ${\mathcal{O}}( \varepsilon ^{0} )$ dynamics because while $v_{-1}, \, w_{-1} \sim {\mathcal{O}}( \varepsilon ^{-1} )$ , the boundary layer width, $\,\text{d}{x} \sim {\mathcal{O}}( \varepsilon )$ . The tangential velocities act exactly like Dirac- $\delta$ distributions at the boundary of the bulk, $X=\Gamma$ .

To further constrain $q_{Y}$ and $q_{Z}$ , we integrate the continuity equation

(3.29) \begin{eqnarray} \partial _{Y}q_{Y} + \partial _{Z}q_{Z} = u_{0}\big |_{x=-\infty } = U_{0}\big |_{X=\Gamma }. \end{eqnarray}

The integration limit $u_{0}\big |_{x=0} = 0$ because of impenetrability at the physical container wall. Equation (3.29) expresses mass conservation between the boundary layer and the interior.

To close the system, we need another relation between $q_{Y}$ and $q_{Z}$ , which comes from directly solving (3.20)–(3.23). We represent the horizontal velocities strictly in terms of the pressure

(3.30) \begin{align} u_{0} &= -\partial _{Y}p_{0} + \partial _{x}^{3}p_{0}, \end{align}

(3.31) \begin{align} v_{-1} &= \partial _{x}p_{0}. \end{align}

The following coupled system relates the pressure and vertical velocity:

(3.32) \begin{align} \partial _{Z}w_{-1} + \partial _{x}^{4}p_{0} &= 0, \end{align}

(3.33) \begin{align} - \partial _{x}^{2}w_{-1} + \partial _{Z}p_{0} &= \Theta _{0}, \end{align}

which collapses into a single system for the vertical velocity,

(3.34) \begin{eqnarray} \left (\partial _{Z}^{2} + \partial _{x}^{6}\right ) w_{-1} = 0. \end{eqnarray}

3.2.1. Bulk onset comparison

Equation (3.34) should look familiar; it is half of the traditional asymptotic stability condition in the rapidly rotating regime when ${\textit{Ra}} \sim {\textit{E}}^{\,4/3}$ . In terms of the $R$ and $\varepsilon$ parameters, the full bulk leading-order balance,

(3.35) \begin{eqnarray} \left (\partial _{Z}^{2} + \partial _{x}^{6}\right ) w_{-1} = \varepsilon \, R \, \partial _{x}^{2} w_{-1} . \end{eqnarray}

If $\varepsilon \, R \sim {\mathcal{O}}( 1 )$ , the left-hand side balances with $\partial _{x}^{2} \to -k^{2}$ , $\partial _{Z}^{2} \to -\pi ^{2}$ and we have the traditional stability condition for rotating convection (Chandrasekhar Reference Chandrasekhar1961),

(3.36) \begin{eqnarray} \text{Bulk onset:} \quad \varepsilon \, R_{c} \sim \frac {\pi ^{2} + k^{6}}{k^{2}} \ge \frac {3 \pi ^{4/3}}{2^{2/3}} \quad \text{most unstable at} \quad k_{c} = \frac {\pi ^{1/3}}{2^{1/6}}. \end{eqnarray}

Here, $R_{c}$ is the critical reduced Rayleigh number and $k_{c}$ is the critical wavenumber for instability.

Equation (3.36) represents a triple balance between vortex stretching (i.e. $\partial _{Z}^{2}$ ) and diffusion (i.e. $\partial _{x}^{6}$ ) on the left, with buoyancy on the right. However, if $R \sim {\mathcal{O}}( 1 )$ , the right-hand side of (3.35) drops out and the bulk no longer supports harmonic $x$ dependence. In that case, the only possible solutions happen with exponential $x$ dependence, which can only remain bounded near a wall.

3.2.2. Boundary layer solutions

We solve the system for no-slip physical sidewall boundary conditions,

(3.37) \begin{eqnarray} && u_{0} = v_{-1} = w_{-1} = 0 \quad \text{at} \; x=0. \end{eqnarray}

We could alternatively use stress-free conditions; however, doing so gives exactly the same final relationship between $q_{Y}$ and $q_{Z}$ . We also impose $w_{-1} = 0$ at $Z=0,1$ .

Solving (3.32) and (3.33) requires Fourier decomposing the velocities, pressure and temperature in the $Z$ direction,

(3.38) \begin{align} \Theta _{0} &= \sum _{n\ge 1}\widehat {\Theta }_{n}(X,Y)\,\sin (n \pi Z), \end{align}

(3.39) \begin{align} p_{0} &= \sum _{n\ge 0}\hat {p}_{n}(x,X,Y)\, \cos (n \pi Z), \end{align}

(3.40) \begin{align} u_{0} &= \sum _{n\ge 0}\hat {u}_{n}(x,X,Y)\cos (n \pi Z), \end{align}

(3.41) \begin{align} v_{-1} &= \sum _{n\ge 0}\hat {v}_{n}(x,X,Y)\, \cos (n \pi Z), \end{align}

(3.42) \begin{align} w_{-1} &= \sum _{n\ge 1}\hat {w}_{n}(x,X,Y)\, \sin (n \pi Z). \end{align}

Solving for the Fourier coefficients, which decay as $x\to -\infty$ and satisfy either no-slip or stress-free boundary conditions at $x=0$ , gives

(3.43) \begin{align} \hat {p}_{n} &= \widehat {P}_{n} + \widehat {Q}_{n} \, \mathcal{X}(\alpha _{n} x) , \end{align}

(3.44) \begin{align} \hat {u}_{n} &= \widehat {U}_{n}\left ( 1 - \mathcal{X}(\alpha _{n} x) \right ), \end{align}

(3.45) \begin{align} \hat {v}_{n} &= \alpha _{n} \, \widehat {Q}_{n} \,\mathcal{X}'(\alpha _{n} x), \end{align}

(3.46) \begin{align} \hat {w}_{n} &= \alpha _{n} \, \widehat {Q}_{n} \,\mathcal{X}'(\alpha _{n} x), \end{align}

where $\alpha _{n} = (n\pi )^{1/3}$ . The boundary layer profile function,

(3.47) \begin{eqnarray} \mathcal{X}(x) = \frac {e^{ \frac {x}{2} } \cos \left( \dfrac {\sqrt {3}\,x}{2} + \delta \right) } {\cos (\delta )}, \end{eqnarray}

where the phase $\delta = \pm \pi /6$ depends on the tangential boundary conditions. In both cases,

(3.48) \begin{eqnarray} \mathcal{X}''(x) - \mathcal{X}'(x) + \mathcal{X}(x) = 0, \quad \mathcal{X}(0) = 1, \end{eqnarray}

which satisfies the normal velocity, $\hat {u}_{n} = 0$ at $x=0$ . For the tangential velocities,

(3.49) \begin{align} & \,\,\quad\text{No-slip:} \quad \hat {v}_{n} = \hat {w}_{n} = 0 \quad \iff \quad \delta = +\frac {\pi }{6}, \quad \mathcal{X}'(0) = 0, \quad \end{align}

(3.50) \begin{align} &\text{Stress-free:} \quad \partial _{x}\hat {v}_{n} = \partial _{x}\hat {w}_{n} = 0 \quad \iff \quad \delta = -\frac {\pi }{6}, \quad \mathcal{X}''(0) = 0. \end{align}

Figure 3 shows the velocity amplitudes for the no-slip case.

The capital-letter variables, $\widehat {U}_{n}$ , $\widehat {P}_{n}$ and $\widehat {\Theta }_{n}$ are the Fourier transforms of the interior bulk quantities that satisfy thermal wind balance at $X=\Gamma$ ,

(3.51) \begin{eqnarray} \widehat {U}_{n} = - \partial _{Y} \widehat {P}_{n}, \quad - n \pi \widehat {P}_{n} = \widehat {\Theta }_{n}. \end{eqnarray}

The overall amplitude satisfies

(3.52) \begin{eqnarray} (\partial _{Y} + n \pi )\, \widehat {Q}_{n} = \widehat {U}_{n}, \end{eqnarray}

which is the Fourier-space version of mass conservation at the bulk boundary (3.29).

Figure 3. Thermal Stewartson boundary layer solutions for no-slip tangential velocities: the blue line shows $1-\mathcal{X}(x)$ . The orange line shows $\mathcal{X}'(x)$ . The orange line compares well with figure 2(b) from Buell & Catton (Reference Buell and Catton1983) for ${\textit{Ta}} \approx 1400^{2}$ .

There are a couple of notable comments regarding (3.43)–(3.46). The first is that both $\hat {v}_{n}$ and $\hat {w}_{n}$ exist strictly within the boundary layer and do not continue into the bulk interior. The second comment is that $\hat {p}_{n}$ and $\hat {u}_{n}$ both contain components that are independent of $x$ . These terms take precisely the same form as the interior geostrophic and hydrostatic balance found in § 3.1.

The most important comment regarding (3.45) and (3.46) is that the coefficients of $\hat {v}_{k,n}$ and $\hat {w}_{k,n}$ take identical form,

(3.53) \begin{eqnarray} \hat {v}_{n} = \hat {w}_{n} \end{eqnarray}

and

(3.54) \begin{eqnarray} \widehat {Q}_{n} = \int _{\mathbb{R}^{-}} \! \hat {v}_{n}\,\text{d}{x} = \int _{\mathbb{R}^{-}} \! \hat {w}_{n}\,\text{d}{x}. \end{eqnarray}

That is, the $\widehat {Q}_{n}$ amplitude gives the Fourier coefficients for both $q_{Y}$ and $q_{Z}$ , closing the system of equations. However, we still need to collect all the relevant information and pose the system in a local abstract form, not in Fourier space.

Finally, regarding the bulk, all relevant quantities do not depend on the value of $\delta = \pm \pi /6$ . This observation accords with the findings from Liao et al. (Reference Liao, Zhang and Chang2006) that the leading-order stability is independent of the sidewall velocity boundary conditions.

3.2.3. Hilbert transform

While the Fourier amplitudes of $v_{-1}$ , $q_{Y}$ , $w_{-1}$ and $q_{Z}$ are all identical, the local functions are not equal because $v_{-1}$ and $q_{Y}$ are Fourier cosine series in $Z$ , but $w_{-1}$ and $q_{Z}$ are Fourier sine series. We need an operator to transform between the different parities, which naturally leads to the Hilbert transform, $\mathcal{H}$ , in the $Z$ direction. For the present analysis, the utility of the Hilbert transform lies in its action on a Fourier sine and cosine basis. For $n \ne 0$ ,

(3.55) \begin{eqnarray} \mathcal{H}\!\left [ \,\sin (n \pi Z)\, \right ] = -\cos (n \pi Z), \quad \mathcal{H}\!\left [ \cos (n \pi Z) \right ] = + \sin (n \pi Z), \end{eqnarray}

which alters the parity of the basis elements in the same manner as a derivative, but without the multiplication by the wavenumber, $n \pi$ . Additionally, the Hilbert transform possesses the following useful properties

(3.56) \begin{eqnarray} \mathcal{H}^{2}[f] = -f + \left \langle f\right \rangle , \quad \partial _{Z}\mathcal{H}\!\left [ f \right ] = \mathcal{H}\!\left [ \partial _{Z}f \right ]. \end{eqnarray}

The Benjamin–Ono equation (Benjamin Reference Benjamin1967; Ono Reference Ono1975) for small-amplitude surface gravity waves over an infinitely deep interior employs the Hilbert transform for reasons similar to the current context. In particular, (3.56) implies

(3.57) \begin{eqnarray} \mathcal{H}\partial _{Z} \equiv |\partial _{Z}|, \end{eqnarray}

for trigonometric functions. Therefore, we may factor the Laplacian-like operator in (3.34)

(3.58) \begin{eqnarray} \left (\partial _{Z}^{2} + \partial _{x}^{6}\right ) = \left (|\partial _{Z}| + \partial _{x}^{3}\right )\!\left (-|\partial _{Z}| + \partial _{x}^{3}\right ), \end{eqnarray}

which separates modes that decay as $x \to \infty$ versus $x \to -\infty$ .

Our final closure for the boundary momentum fluxes is

(3.59) \begin{eqnarray} q_{Y} = - \mathcal{H}\!\left [ q_{Z} \right ], \end{eqnarray}

which means they share the same Fourier coefficients but with different vertical parity.

3.2.4. Barotropic mode

Finally, we derive properties of the barotropic mode in the Stewartson layer. We first take the vertical average of the horizontal velocities

(3.60) \begin{align} \left \langle v_{-1}\right \rangle &= \partial _{x}\!\left \langle p_{0}\right \rangle\!, \end{align}

(3.61) \begin{align} \left \langle u_{0}\right \rangle &= -\partial _{Y} \!\left \langle p_{0}\right \rangle + \partial _{x}^{3}\!\left \langle p_{0}\right \rangle\!, \end{align}

(3.62) \begin{align} \partial _{x}\!\left \langle u_{0}\right \rangle &+ \partial _{Y}\!\left \langle v_{-1}\right \rangle = \partial _{x}^{4}\!\left \langle p_{0}\right \rangle = 0. \end{align}

Therefore, $\!\left \langle p_{0}\right \rangle = \!\left \langle P_{0}\right \rangle$ within the sidewall layer, with $\!\left \langle u_{0}\right \rangle = -\partial _{Y} \!\left \langle P_{0}\right \rangle$ . Impenetrability at the physical container boundary, along with fixing the overall pressure gauge condition, implies

(3.63) \begin{eqnarray} \!\left \langle P_{0}\right \rangle = 0 \quad \text{at}\quad \; x=0, \; X=\Gamma , \end{eqnarray}

which is the only sidewall boundary condition needed for no-slip top and bottom boundary conditions.

Free-slip top and bottom boundaries also require the next-order barotropic momentum balances and continuity in the sidewall layer,

(3.64) \begin{align} \!\left \langle f\right \rangle + \partial _{Y} \!\left \langle p_{1}\right \rangle + \!\left \langle u_{1}\right \rangle &= \partial _{x}^{2} \!\left \langle v_{0}\right \rangle\!, \end{align}

(3.65) \begin{align} \partial _{x} \!\left \langle p_{1}\right \rangle + \partial _{X} \!\left \langle P_{0}\right \rangle - \!\left \langle v_{0}\right \rangle &= 0, \end{align}

(3.66) \begin{align} \partial _{x}\!\left \langle u_{1}\right \rangle + \partial _{X}\!\left \langle U_{0}\right \rangle + \partial _{Y} \!\left \langle v_{0}\right \rangle &= 0, \end{align}

where

(3.67) \begin{eqnarray} \!\left \langle f\right \rangle = \frac {\partial _{x} \!\left \langle u_{0}\, v_{-1}\right \rangle + \partial _{Y}\!\left \langle v_{-1}\, v_{-1}\right \rangle }{\sigma }. \end{eqnarray}

The tangential solution in the sidewall layer,

(3.68) \begin{eqnarray} \!\left \langle v_{0}\right \rangle = \partial _{X}\!\left \langle P_{0}\right \rangle - \int _{-\infty }^{x} (x'-x) \, \langle \,f(x') \rangle \,\text{d}{x'}. \end{eqnarray}

A no-slip boundary condition, $\!\left \langle v_{0}\right \rangle = 0$ , at $x=0$ , $X=\Gamma$ implies

(3.69) \begin{align} \sigma \, \partial _{X}\!\left \langle P_{0}\right \rangle &= \int _{\mathbb{R}^{-}} \!\! x \, \left [\, \partial _{x} \!\left \langle u_{0}\, v_{-1}\right \rangle + \partial _{Y}\!\left \langle v_{-1}\, v_{-1}\right \rangle \, \right ] \,\text{d}{x} \end{align}

(3.70) \begin{align} \nonumber &= \ \!\left \langle \, U_{0} \, q_{Y} \right \rangle \int _{\mathbb{R}^{-}} \!\! x \, \partial _{x} \left [ (1-\mathcal{X}(x))\, \mathcal{X}'(x) \right ] \,\text{d}{x}\, + \\ & \quad \ 2\,\!\left \langle \, q_{Y} \, \partial _{Y} q_{Y} \right \rangle \int _{\mathbb{R}^{-}} \!\! x \, \mathcal{X}'(x)^{2} \,\text{d}{x}. \end{align}

Computing the coefficients,

(3.71) \begin{eqnarray} \int _{\mathbb{R}^{-}}\!\! x \, \partial _{x} \left [ (1-\mathcal{X}(x))\, \mathcal{X}'(x) \right ] \,\text{d}{x} = - \frac {1}{2}, \quad \int _{\mathbb{R}^{-}} \!\! x \, \mathcal{X}'(x)^{2} \,\text{d}{x} = - \frac {3}{4}. \end{eqnarray}

Both coefficients are independent of $\delta = \pm \pi /6$ . Putting everything together,

(3.72) \begin{eqnarray} \partial _{X}\!\left \langle P_{0}\right \rangle = - \frac { 3\!\left \langle \, q_{Y}\, \partial _{Y} q_{Y} \right \rangle-\left \langle \, q_{Y} \,\partial _{Y} P_{0} \right \rangle }{2\sigma } \quad \text{at}\quad \; X=\Gamma . \end{eqnarray}

3.3. Region-III (Ekman layers)

We still need to fix the barotropic pressure, $\!\left \langle P_{0}\right \rangle$ , in the interior, which requires considering the vertical vorticity along with the top and bottom boundary conditions. Geostrophic balance relates the vertical vorticity and the pressure,

(3.73) \begin{eqnarray} \zeta _{0} \equiv \partial _{X}V_{0} - \partial _{Y} U_{0} = \left (\partial _{X}^{2} + \partial _{Y}^{2}\right ) P_{0}. \end{eqnarray}

Taking the horizontal curl of the momentum equations and using the continuity equation,

(3.74) \begin{eqnarray} \frac {1}{\sigma }\left (\partial _{t} + U_{0} \partial _{X} + V_{0} \partial _{Y} \right )\zeta _{0} - \left (\partial _{X}^{2} + \partial _{Y}^{2} \right ) \zeta _{0} = \partial _{Z}W_{3} + \partial _{Z}^{2} \zeta _{0}. \end{eqnarray}

Notably, the leading-order rotational terms enforce

(3.75) \begin{eqnarray} \partial _{Z}W_{i} = 0 \quad \text{for}\quad \; i = 0,\, \frac {1}{2}, \,1, \, \frac {3}{2}, \, 2,\, \frac {5}{2}. \end{eqnarray}

The first non-constant vertical velocity is $W_{3}$ .

3.3.1. Stress-free boundaries

For stress-free boundaries (2.11), applying the averaging operator to (3.74) gives

(3.76) \begin{eqnarray} {({z-SF}): \,}\quad \partial _{t}\!\left \langle \zeta _{0}\right \rangle + \partial _{X}\!\left \langle U_{0} \,\zeta _{0}\right \rangle + \partial _{Y}\!\left \langle V_{0} \, \zeta _{0} \right \rangle = \sigma \left (\partial _{X}^{2} + \partial _{Y}^{2}\right )\!\left \langle \zeta _{0} \right \rangle\!. \end{eqnarray}

The $\partial _{Z}$ terms on the right-hand side vanish because of impenetrability and vanishing stress,

(3.77) \begin{align} W_{3}(Z=0) &= W_{3}(Z=1) = 0, \end{align}

(3.78) \begin{align} \partial _{Z}\zeta _{0}(Z=0) &= \partial _{Z}\zeta _{0}(Z=1) = 0. \end{align}

Equation (3.76) determines the missing depth-independent pressure component needed to obtain the velocities in (3.14). The depth-dependent pressure (hence temperature) forces (3.76) through the nonlinear baroclinic interactions, e.g. $\!\left \langle U_{0} \zeta _{0}\right \rangle - \!\left \langle U_{0}\right \rangle \!\left \langle \zeta _{0}\right \rangle$ , and the two systems are coupled. We solve (3.76) together with boundary conditions $\!\left \langle P_{0}\right \rangle =0$ and (3.72).

3.3.2. No-slip boundaries

For no-slip top and bottom boundaries, one expects an Ekman-pumping effect of the order (Zhang & Roberts Reference Zhang and Roberts1998)

(3.79) \begin{eqnarray} W_{{Ekman}} \sim {\textit{E}}^{\,1/2}\, \zeta _{0} = \varepsilon ^{3/2} \, \zeta _{0}. \end{eqnarray}

Together, (3.75) and (3.79) imply a depth-independent Ekman velocity,

(3.80) \begin{eqnarray} \partial _{Z}W_{3/2} = 0. \end{eqnarray}

The Ekman layers determine the barotropic pressure.

Given an interior flow, the Ekman boundary layer response follows as a well-known exercise (Greenspan Reference Greenspan1969; Pedlosky Reference Pedlosky1987). Region-III lies beyond the sidewall Stewartson layer, and hence $\partial _{x} \to 0$ . Within the top and bottom boundary layers,

(3.81) \begin{align} \partial _{z}P_{0} &= 0, \end{align}

(3.82) \begin{align} - v_{0} + \partial _{X}P_{0} &= \partial _{z}^{2}u_{0}, \end{align}

(3.83) \begin{align} u_{0} + \partial _{Y}P_{0} &= \partial _{z}^{2}v_{0}, \end{align}

(3.84) \begin{align} \partial _{X} u_{0} + \partial _{Y}v_{0} + \partial _{z} w_{3/2} &= 0. \end{align}

The solutions to (3.81)–(3.84) are well known for $Z=0$ and the boundary layer coordinate $0 \le z \lt \infty$ , and for $Z=1$ and the boundary layer coordinate $-\infty \lt z \le 0$ . In both cases,

(3.85) \begin{eqnarray} \left [ \begin{array}{c} u_{0} \\ v_{0} \\ \end{array} \right ] = \left [ \begin{array}{rr} \text{se}(z) & \text{ce}(z) \\ -\text{ce}(z) & \text{se}(z) \\ \end{array} \right ]\boldsymbol{\cdot} \left [ \begin{array}{c} \partial _{X}P_{0} \\ \partial _{Y}P_{0} \\ \end{array} \right ], \end{eqnarray}

where $P_{0}$ is evaluated at $Z=0$ or $Z=1$ depending on the boundary layer and the functions,

(3.86) \begin{eqnarray} \text{se}(z) = e^{- \frac {|z|}{\sqrt {2}}} \sin \left (\dfrac {|z|}{\sqrt {2}}\right ), \quad \text{ce}(z) = 1- e^{- \frac {|z|}{\sqrt {2}}} \cos \left (\dfrac {|z|}{\sqrt {2}}\right ). \end{eqnarray}

Both $\text{se}(0)=\text{ce}(0)=0$ .

For the vertical velocity at either boundary,

(3.87) \begin{eqnarray} w_{3/2} = \text{sign}(z)\, \frac {\zeta _{0}}{\sqrt {2}} \left (\text{ce}(z) - \text{se}(z) \right ). \end{eqnarray}

The Ekman pumping velocity must connect to the interior as $|z| \to \infty$ ,

(3.88) \begin{eqnarray} \lim _{|z| \to \infty } \! w_{3/2} = \text{sign}(z) \frac {\zeta _{0}}{\sqrt {2}}. \end{eqnarray}

However, in the interior, $\partial _{Z}W_{3/2} = 0$ . Therefore,

(3.89) \begin{eqnarray} \zeta _{0}(X,Y,Z=0) + \zeta _{0}(X,Y,Z=1) = 0, \end{eqnarray}

which indirectly fixes the barotropic pressure. In terms of the local pressure,

(3.90) \begin{eqnarray} P_{0}|_{Z=0} + P_{0}|_{Z=1} = 2\,\Phi (X,Y), \end{eqnarray}

where $\Phi (X, Y)$ represents a two-dimensional harmonic function (the factor of two provides a later simplification). The Stewartson layer supports no barotropic pressure, and $\Phi$ ensures $\left \langle P_{0}\right \rangle$ vanishes at the sidewall.

Given hydrostatic balance, $\partial _{Z}P_{0}(Z) = \Theta _{0}(Z)$ , the pressure decomposes into barotropic and baroclinic components,

(3.91) \begin{eqnarray} P_{0}(Z) = \left \langle P_{0}\right \rangle + \int _{0}^{Z} \! Z_{0}\, \Theta _{0}(Z_{0})\,\text{d}{Z_{0}} - \int _{Z}^{1} \! (1-Z_{0})\, \Theta _{0}(Z_{0})\,\text{d}{Z_{0}}. \end{eqnarray}

Applying (3.90), the barotropic pressure is

(3.92) \begin{eqnarray} \!\left \langle P_{0}\right \rangle = \Phi (X,Y) - \int _{0}^{1} \! \left (Z-\frac {1}{2}\right )\Theta _{0}(X,Y,Z) \,\text{d}{Z}. \end{eqnarray}

The Stewartson layer equations imply $\!\left \langle P_{0}\right \rangle = 0$ within the side wall layer, thus fixing the harmonic function, $\Phi (X, Y)$ , when applying (3.92) at $X=\Gamma$ .

While the Ekman layers in the bulk dictate leading-order dynamics of the barotropic mode, the Ekman layers within the sidewall Stewartson layers generate a lower-order passive response. That is, it is possible to compute the effects of Ekman pumping after the fact (e.g § 5.2), but we can press ahead with calculating the leading-order dynamics if we are not interested in the correction terms.

There is, however, one notable aspect. The balances in the sidewall naturally produce velocities $u_{0}, v_{-1}, w_{-1}$ . That is, the tangential velocities are ${\mathcal{O}}( {\textit{E}}^{\,-1/3} )$ compared with ${\mathcal{O}}( 1 )$ for the normal component. However, the normal component achieves its largest amplitudes in the form of an Ekman response.

Like before, $\partial _{z}p_{0} = 0$ . Now, the horizontal velocities are isotropic,

(3.93) \begin{align} - v_{-1} + \partial _{x}p_{0} &= \partial _{z}^{2} \, u_{-1}, \end{align}

(3.94) \begin{align} u_{-1} &= \partial _{z}^{2}\, v_{-1} . \end{align}

The Ekman flux is now larger than within the bulk, but is still subdominant compared with the $w_{-1}$ convection; i.e. $\partial _{x} u_{-1} + \partial _{z} w_{-1/2} =0$ ,

(3.95) \begin{eqnarray} w_{-1/2} = \text{sign}(z)\, \frac {\partial _{x}^{2}p_{0}}{\sqrt {2}} \left (\text{ce}(z) - \text{se}(z) \right ). \end{eqnarray}

The $w_{-1/2}$ does not automatically satisfy the no-slip condition at $x=0$ , which requires an ${\mathcal{O}}( {\textit{E}}^{\,1/2} )$ corner layer (Kerswell & Barenghi Reference Kerswell and Barenghi1995) that remains passive with respect to the leading-order dynamics. Solving the corner layer requires a somewhat technical direct numerical calculation (Burns et al. Reference Burns, Fortunato, Julien and Vasil2022). While the sidewall Ekman response remains passive in the fully asymptotic regime, we should expect noticeable corrections to heat transport as Rayleigh number transitions out of a strongly rotationally dominated (asymptotic) regime (Stellmach et al. Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014; Julien et al. Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016).

4.1. Geometry and operators

The container is now a general upright direct product of a $z$ direction with a smooth but otherwise arbitrary horizontal area,

(4.1) \begin{eqnarray} \mathcal{V} = \mathcal{A} \times [\, 0,H \, ], \end{eqnarray}

with e.g. Cartesian horizontal coordinates $(x,y) \in \mathcal{A}$ , and vertical coordinate $z \in [0,H]$ , with local unit vector, $\hat {z}$ . The boundary of the horizontal area is $\partial \mathcal{A}$ with the local outward normal vector, $\hat {n}$ . In the case of a general smooth $\partial \mathcal{A}$ , the analysis in § 3.2 would follow the same reasoning with $x$ representing a local signed-distance coordinate (Hester & Vasil Reference Hester and Vasil2023). We represent the coordinate along the tangential direction to $\partial \mathcal{A}$ as $\ell$ , with local unit vector, $\hat {\ell } = \hat {z} \times \hat {n}$ , assuming a right-handed coordinate system. For short-hand, we denote derivatives, $\partial _{z} = \hat {z} \boldsymbol{\cdot} \boldsymbol{\nabla}$ , $\partial _{n} = \hat {n} \boldsymbol{\cdot} \boldsymbol{\nabla}$ and $\partial _{\ell } = \hat {\ell } \boldsymbol{\cdot} \boldsymbol{\nabla}$ . In the two most relevant cases,

(4.2) \begin{align} \text{Cartesian:} \quad & \partial _{n} = \partial _{x}, \quad \partial _{\ell } = \partial _{y}, \end{align}

(4.3) \begin{align} \text{Cylindrical:} \quad & \partial _{n} = \partial _{r}, \quad \partial _{\ell } = \frac {1}{r}\,\partial _{\phi }. \end{align}

We define the respective barotropic average and Hilbert transform,

(4.4) \begin{eqnarray} \!\left \langle p\right \rangle = \frac {1}{H}\int _{0}^{H} \! \! p(z') \,\text{d}{z'} , \quad \mathcal{H}\!\left [ q \right ](z) = \text{p.v.} \!\left \langle q(z')\, \cot \left( \dfrac {\pi (z'-z)}{H}\right) \right \rangle , \end{eqnarray}

where p.v. denotes the principle value in the Hilbert transform integral.

4.2. Dynamics

For $(x,y) \in \mathcal{A}$ ,

(4.5) \begin{align} 2 \Omega \, \partial _{z}u = g \alpha \, \hat {z} \times \boldsymbol{\nabla} T, \quad & \partial _{t} T + u \boldsymbol{\cdot} \boldsymbol{\nabla} T = \kappa \nabla ^{2} T, \end{align}

(4.6) \begin{align} T\big |_{z=0} = T_{{bot}}, \quad & T\big |_{z=H} = T_{{top}}. \end{align}

For $(x,y) \in \partial \mathcal{A}$ ,

(4.7) \begin{align} q_{\ell } = - \mathcal{H}\!\left [ q_{z} \right ], \quad \partial _{\ell } q_{\ell } + \partial _{z} q_{z} = u_{n}, \quad& (q_{\ell } \,\partial _{\ell } + q_{z}\, \partial _{z}) \, T = - \kappa \, \partial _{n} T , \end{align}

(4.8) \begin{align} q_{z} = 0 \quad \text{for}\quad\& \quad z = 0,H. \end{align}

For vertical no-slip,

(4.9) \begin{eqnarray} \zeta \big |_{z=0} + \zeta \big |_{z=H} = 0. \end{eqnarray}

For vertical stress-free,

(4.10) \begin{eqnarray} \partial _{t}\!\left \langle \zeta \right \rangle + \boldsymbol{\nabla} \boldsymbol{\cdot} \, \!\left \langle u \,\zeta \right \rangle = \nu \, \nabla^{2} \!\left \langle \zeta \right \rangle , \quad -2\nu \, \!\left \langle u_{\ell }\right \rangle \! \big |_{\partial \mathcal{A}} = \!\left \langle \, q_{\ell } \left ( 4\,\partial _{\ell }q_{\ell } + \partial _{z}q_{z} \right )\,\right \rangle . \end{eqnarray}

4.3. Heat transport

Heat transport is an essential diagnostic in convection-dominated systems, with much work for laboratory and natural systems still ongoing (Doering, Toppaladoddi & Wettlaufer Reference Doering, Toppaladoddi and Wettlaufer2019; Schumacher & Sreenivasan Reference Schumacher and Sreenivasan2020; Vasil, Julien & Featherstone Reference Vasil, Julien and Featherstone2021; Lohse & Shishkina Reference Lohse and Shishkina2024). Rapid rotation adds many additional nuances (Julien et al. Reference Julien, Knobloch, Rubio and Vasil2012; Stellmach et al. Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014; Julien et al. Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016). In studies of bulk convection, heat transport from wall modes tends to confound measurements, and efforts have attempted to quantify and/or limit their influence (Ecke et al. Reference Ecke, Zhang and Shishkina2022; Terrien, Favier & Knobloch Reference Terrien, Favier and Knobloch2023; Zhang et al. Reference Zhang, Reiter, Shishkina and Ecke2024).

At first glance, the reduced wall-mode equations do not appear to contain any means for vertical heat transport; there is no bulk vertical velocity. The resolution occurs when taking into account the sidewall heat fluxes. All the vertical heat transport happens in the sidewall layer, which manifests within the bulk via $q_{z}$ .

Defining the horizontal average, $\overline {T}$ , the mean thermal equation,

(4.11) \begin{eqnarray} \partial _{t}\overline {T} + \partial _{z} \left [ \frac {1}{|\mathcal{A}|}\oint _{\partial \mathcal{A}} \! q_{z} T \,\text{d}{\ell } - \kappa \, \partial _{z} \overline {T}\right ] = 0. \end{eqnarray}

The integrated boundary term follows straightforwardly from integrating the bulk equation and applying (4.7). When the system is in a statistically steady state, the Nusselt number is a global constant with

(4.12) \begin{eqnarray} {\textit{Nu}} = \frac {1}{|\mathcal{A}|}\oint _{\partial \mathcal{A}} \! q_{z} \Theta \,\text{d}{\ell } - \kappa \, \partial _{z} \overline {T}, \end{eqnarray}

where in this case, $\Theta = T - \overline {T}$ .

When considering the boundary terms in (3.15), we also now find the non-trivial wall-modes ‘power integral’ (Howard Reference Howard1963) that regulates overall input and dissipation,

(4.13) \begin{eqnarray} \frac {1}{2} \frac {d}{{\textrm d}t} \!\big\langle \overline {\Theta ^{2}} \big\rangle = \frac { 1 }{|\mathcal{A}|} \oint _{\partial \mathcal{A}} \!\left \langle q_{z} \,\partial _{z}\overline {T}\, \Theta \right \rangle \,\text{d}{\ell } - \kappa \!\left \langle \overline {|\boldsymbol{\nabla} \Theta |^{2}}\right \rangle . \end{eqnarray}

In both cases, it is clear the vertical velocity behaves precisely as a Dirac- $\delta$ distribution supported on the boundary; e.g. in local Cartesian coordinates with the wall at $x=0$ ,

(4.14) \begin{eqnarray} w(x,y,z) = q_{z}(y,z) \, \delta (x), \end{eqnarray}

where the $\delta$ -function replaces the detailed behaviour within the sidewall layer.

5.1. Semi-infinite channel

This section reproduces the final results of Herrmann & Busse (Reference Herrmann and Busse1993) with the advantage that now, all aspects are physically apparent. We assume a semi-infinite domain with the wall-normal coordinate, $- \infty \lt x \le 0$ , and harmonic dependence along the wall, $\sim e^{i k y}$ , with wavenumber $k$ .

The following exponential profile satisfies the boundary condition at $x=0$ :

(5.6) \begin{eqnarray} \theta=e^{\,\beta \, x + i\, (k y +\omega t) } \sin (\pi z) + \text{c.c.} \quad \text{where}\quad \; \beta = \frac {k (k+i \pi ) R }{\pi \left (k^{2}+\pi ^{2}\right )}. \end{eqnarray}

Assuming $R\gt 0$ , the temperature perturbations decays as $x\to -\infty$ . The bulk evolution equation gives the dispersion relation for the complex-valued frequency

(5.7) \begin{eqnarray} \omega &=& \frac {2k^{3}R^{2} }{\pi (\pi ^{2}+k^{2})^{2} } - i\, \frac {k^{2} (k^{2}-\pi ^{2}) R^{2} - \pi ^{2} (\pi ^{2}+k^{2})^{3}}{\pi ^{2} (\pi ^{2}+k^{2})^{2}}. \end{eqnarray}

The real-valued growth rate,

(5.8) \begin{eqnarray} \gamma = -\textrm{Im}(\omega ) = \frac {k^{2} (k^{2}-\pi ^{2}) R^{2} - \pi ^{2} (\pi ^{2}+k^{2})^{3}}{\pi ^{2} (\pi ^{2}+k^{2})^{2}}. \end{eqnarray}

The growth rate is positive for

(5.9) \begin{eqnarray} R \ge R_{{c}}(k) = \frac {\pi (k^{2}+\pi ^{2})^{3/2}}{k (k^{2}-\pi ^{2})^{1/2} } \ge R_{{c}} (k_{{c}}) = \pi ^{2}\sqrt {6\sqrt {3}} \ \approx \ 31.8167. \end{eqnarray}

Zhang et al. (Reference Zhang, Reiter, Shishkina and Ecke2024) fit a quadratic curve through their stability data near the minimum of the form

(5.10) \begin{align} \frac {R_{{c}}(k)}{R_{{c}}(k_{{c}})} = 1 + \xi ^{2} ( k - k_{{c}})^{2} + {\mathcal{O}}((k - k_{{c}})^{3} ). \end{align}

They find $\xi \approx 0.18$ for ${\textit{E}} = 10^{-6}$ . Equation (5.9) implies $\xi = \sqrt {2 - \sqrt {3}}/\pi \approx 0.164769$ .

The respective critical wavenumber and frequency are

(5.11) \begin{eqnarray} k_{{c}}=\pi \sqrt {2+\sqrt {3}} \approx 6.0690, \quad \omega _{{c}} = 2 \pi ^{2}\sqrt {3 (2 + \sqrt {3})} \approx 66.0487. \end{eqnarray}

The sign convention for frequency means that $\omega \gt 0$ moves retrograde to the background rotation, $\Omega$ . The respective phase and group speeds,

(5.12) \begin{eqnarray} \frac {\omega _{{c}}}{k_{{c}}} = 2 \sqrt {3} \, \pi \approx 10.8828, \quad \frac {\partial \omega }{\partial k}\Big |_{k_{{c}},R_{{c}}} = 2 (\sqrt {3} - 2)\,\pi \approx -1.68357, \end{eqnarray}

within $\approx \,$ 20 % of the group velocity of $-2.1$ reported in the laboratory experiments of Ning & Ecke (Reference Ning and Ecke1993) for ${\textit{E}}=10^{-3}$ . Notably, the group speed is always prograde above onset.

Even though we remove the explicit Stewartson boundary layer, at onset, the temperature perturbations remain quite localised near the wall (on the bulk scale, $x$ ), where

(5.13) \begin{eqnarray} \beta _{{c}} = \pi \sqrt {3+2 \sqrt {3}} + i \, \pi \,3^{1/4} \ \approx \ 7.9873 + 4.1345 \, i. \end{eqnarray}

In physical units, the critical wavelength along the wall is approximately $\lambda _{y} \approx 1.04\, H$ and the $e$ -folding length away from the wall is $\lambda _{x} \approx 0.13\, H$ .

For large $R$ , the growth rate,

(5.14) \begin{eqnarray} \gamma \sim \frac {R^{2}}{\pi ^{2}} -\frac {k^{4}+3 R^{2}}{k^{2}} + {\mathcal{O}}( 1 ) \quad \text{as} \;\quad R \to \infty . \end{eqnarray}

The fastest-growing wavenumber, $k$ maximises the second term with $k \sim 3^{1/4} \sqrt {R}$ leaving

(5.15) \begin{eqnarray} \gamma \sim \frac {R^{2}}{\pi ^{2}} - 2 \sqrt {3} R + {\mathcal{O}}( 1 ) \quad \text{as} \;\quad R \to \infty . \end{eqnarray}

The severe dependence of $\gamma$ on $R$ makes nonlinear simulations very stiff at high values of criticality. For example, when $R = 5 \, R_{{c}}(k_{{c}})$ , $\gamma \approx 2500$ . Fast time scales exist for the propagation as well. As $R\to \infty$ , the phase and group speeds are exactly opposite,

(5.16) \begin{eqnarray} c_{p} \sim \frac {2 R}{\sqrt {3} \pi } - \frac {4\pi }{3}, \quad c_{g} \sim -\frac {2 R}{\sqrt {3} \pi } - 4\pi . \end{eqnarray}

Reintroducing physical parameters implies the dimensional growth rate in the large- $R$ limit,

(5.17) \begin{eqnarray} \gamma ^{*} \sim \frac {(g\alpha \Delta T)^{2}}{(2\Omega )^{2} \kappa }, \end{eqnarray}

which is notably independent of the container size.

5.2. Ekman-pumping correction

In § 1, we report the critical Rayleigh numbers, including the leading-order corrections due to Ekman pumping effects in (1.3) and (1.4). The ${\mathcal{O}}({\textit{Ta}}^{7/12})$ correction to the bulk ${\textit{Ra}}$ has an interesting history. The ${\mathcal{O}}({\textit{Ta}}^{5/12})$ correction to the wall ${\textit{Ra}}$ has not, to our knowledge, been calculated before.

In his famous book, Chandrasekhar (Reference Chandrasekhar1961) calculated the leading-order contribution for bulk modes (unaware of wall modes at the time) using sound asymptotic methods,

(5.18) \begin{eqnarray} {\textit{Ra}} \sim (27 \pi ^{4}/4)^{1/3} \, {\textit{Ta}}^{\,2/3} , \quad \text{as} \;\quad {\textit{Ta}} \to \infty . \end{eqnarray}

However, when comparing to numerical calculations up to ${\textit{Ta}} \approx 10^{12}$ , he noticed a discrepancy for no-slip boundary conditions that appeared like an offset in the numerical proportionality factor, commenting ‘…the discussion is not “fine” enough to account for the differences in the constants of proportionality…’ (Chandrasekhar Reference Chandrasekhar1961 § 27d ‘The origin of the $T^{{2}/{3}}$ -law’). The resolution of the issue is that the relative leading-order correction is $\propto {\textit{Ta}}^{-1/12}$ , which produces a noticeable change for realistic numerical and laboratory parameter values. Ekman pumping manifests especially clearly in heat-transport enhancement (King et al. Reference King, Stellmach, Noir, Hansen and Aurnou2009; Stellmach et al. Reference Stellmach, Lischper, Julien, Vasil, Cheng, Ribeiro, King and Aurnou2014; Julien et al. Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016). The upshot is that Chandrasekhar’s discussion was ‘fine’ enough; it simply did not go far enough.

The issue became fully resolved when Zhang & Roberts (Reference Zhang and Roberts1998) showed the leading-order correction results from considering the Ekman pumping boundary conditions of the form we use in § 3.3 (3.88), together with a leading-order interior equation of the form as (3.35). The critical bulk Rayleigh number and wavenumber,

(5.19) \begin{align} {\textit{Ra}} &\sim (27 \pi ^{4}/4)^{1/3} \, {\textit{Ta}}^{\,2/3} - (8192\,\pi ^{4})^{1/6} {\textit{Ta}}^{\,7/12}, \end{align}

(5.20) \begin{align} k &\sim (\pi ^{2}/2)^{1/6} {\textit{Ta}}^{1/6} - (27 \pi ^{4}/4)^{-1/3} {\textit{Ta}}^{1/12}, \end{align}

as ${\textit{Ta}} \to \infty$ . We use analogous techniques to find the equivalent correction for wall modes.

The Ekman flux in the bulk is much too small to produce a significant correction. However, in the sidewall,

(5.21) \begin{align} \partial _{x}p_{1/2} - v_{-1/2} &= 0, \end{align}

(5.22) \begin{align} \partial _{Y} p_{1/2} + u_{1/2} - \partial _{x}^{2}v_{-1/2} &= 0, \end{align}

(5.23) \begin{align} \partial _{Z} p_{1/2} - \partial _{x}^{2} w_{-1/2} &= \Theta _{1/2}, \end{align}

(5.24) \begin{align} \partial _{x}u_{1/2} + \partial _{Y}v_{-1/2} + \partial _{Z}w_{-1/2} &= 0. \end{align}

The boundary condition is the key difference from the leading order,

(5.25) \begin{eqnarray} w_{-1/2} = \mp \frac {\partial _{x}^{2}p_{0}}{\sqrt {2}} \quad \text{at}\quad \; Z = 0, 1. \end{eqnarray}

Projecting the system on to the relevant $2\cos (n\pi Z)$ and $2\sin (n\pi Z)$ functions,

(5.26) \begin{align} \partial _{x}\hat {p}_{n,1/2} - \hat {v}_{n,-1/2} &= 0, \end{align}

(5.27) \begin{align} \partial _{Y} \hat {p}_{n,1/2} + \hat {u}_{n,1/2} - \partial _{x}^{2}\hat {v}_{n,-1/2} &= 0, \end{align}

(5.28) \begin{align} - \alpha _{n}^{3} \, \hat {p}_{n,1/2} - \partial _{x}^{2} \hat {w}_{n,-1/2} &= \widehat {\Theta }_{n,1/2}, \end{align}

(5.29) \begin{align} \partial _{x}\hat {u}_{n,1/2} + \partial _{Y}\hat {v}_{n,-1/2} + \alpha _{n}^{3} \hat {w}_{n,-1/2} &= \! \! \! \! \sum _{m+n = \text{even}} \! \! \! \! \sqrt {2}\,\alpha _{m}^{2} \widehat {Q}_{m} \mathcal{X}''(\alpha _{m} x), \end{align}

where recalling $\alpha _{n} = (n\pi )^{1/3}$ . The right-hand side results from integrating by parts in $Z$ and using the Ekman boundary conditions,

(5.30) \begin{eqnarray} 2 \int _{0}^{1}\cos (n \pi Z) \, \partial _{Z} w(Z) \,\text{d}{Z} = n \pi \, \hat {w}_{n} + 2\left [(-1)^{n} w(1) - w(0) \right ], \end{eqnarray}

where by definition,

(5.31) \begin{eqnarray} \hat {w}_{n} = 2 \int _{0}^{1}\sin (n\pi Z)\, w(Z) \,\text{d}{Z}. \end{eqnarray}

The solution to the system is straightforward but messy. However, we only need the relationship between $q_{Y}$ and $q_{Z}$ to close the system.

The final result in general form (dropping the $X, Y, Z$ bulk notation),

(5.32) \begin{eqnarray} q_{\ell } + \mathcal{H}\!\left [ q_{z} \right ] = -\sqrt {\varepsilon } \, \mathcal{E} [q_{\ell }] \quad \text{where}\quad \; \sqrt {\varepsilon } = {\textit{E}}^{\,1/6}. \end{eqnarray}

The right-hand side contains the highly non-local operator

(5.33) \begin{eqnarray} \widehat {\mathcal{E} [q]}_{n} = \frac {2\sqrt {2}}{\pi ^{2/3}}\sum _{m+n=\text{even}} \frac {m}{(n m)^{1/3} (m+n)} \, \widehat {q}_{m}. \end{eqnarray}

We can now continue to find the ${\mathcal{O}}( \varepsilon ^{1/2} )$ perturbation to $R, \omega ,k$ from the semi-infinite stability problem:

(5.34) \begin{align} \partial _{t}\theta _{1/2} + \partial _{\tau } \theta _{0} &= \nabla^{2}\theta _{1/2} , \end{align}

(5.35) \begin{align} \theta _{1/2} &= \partial _{z}p_{1/2}, \end{align}

(5.36) \begin{align} R_{0}\, q_{z,1/2} + R_{1/2}\, q_{z,0} &= \partial _{x} \theta _{1/2} , \end{align}

(5.37) \begin{align} \partial _{y}q_{y,1/2} + \partial _{z}q_{z,1/2} &= - \partial _{y}p_{1/2} , \end{align}

(5.38) \begin{align} q_{y,1/2} + \mathcal{H}\!\left [ q_{z,1/2} \right ] &= -\frac {\sqrt {2}}{\pi ^{2/3}} q_{y,0}. \end{align}

The system has an extra $R_{1/2}$ term along with an ${\mathcal{O}}( \varepsilon ^{1/2} )$ time scale, $\partial _{\tau }$ . For a given wavenumber along the wall, $k$ ,

(5.39) \begin{eqnarray} R_{1/2} = -\frac {\sqrt {2} \pi ^{1/3 } k \left (k^{2}-3 \pi ^{2}\right ) \sqrt {k^{2}+\pi ^{2}}}{ \left (k^{2}-\pi ^{2}\right )^{3/2}}, \quad \omega _{1/2} = \frac {2\sqrt {2} \,\pi ^{1/3} k \left (k^{2}+\pi ^{2}\right )^{2}}{\left (k^{2} - \pi ^{2}\right )^{2}} . \end{eqnarray}

Optimising the combined reduced Rayleigh number, $R \sim R_{0}(k) + \sqrt {\varepsilon }\, R_{1/2}(k)$ , over $k$ yields

(5.40) \begin{align} R_{c} \sim \sqrt {6\sqrt {3}} \, \pi ^{2} - \sqrt {\varepsilon }\, \sqrt {4\sqrt {3}-6}\, \pi ^{4/3} & \approx 31.8167 - 4.4329 \sqrt {\varepsilon }, \end{align}

(5.41) \begin{align} k_{c} \sim \sqrt {2 + \sqrt {3}}\, \pi + \sqrt {\varepsilon }\,\frac {5\,\left (1+\sqrt {3}\right )}{6} \,\pi ^{1/3} & \approx 6.06909 + 3.33445 \sqrt {\varepsilon }, \end{align}

(5.42) \begin{align} \omega _{c} \sim 2 \sqrt {3 (2+\sqrt {3})} \, \pi ^{2} + \sqrt {\varepsilon }\, \frac {\left (23+13 \sqrt {3}\right )}{3} \pi ^{4/3} & \approx 66.0487 + 69.8097 \sqrt {\varepsilon }, \end{align}

as $\varepsilon \to 0$ . The next section computes the critical onset parameters directly.

5.3. Direct numerical calculations

We validate the asymptotic correction from § 5.2 with high-precision direct numerical calculations of the critical Rayleigh number, wave number and precession frequency over a range of very small but finite Ekman numbers. We solve a linearised version of (2.4)–(2.8) in a finite Cartesian channel assuming $y$ -direction Fourier dependence with a complex-valued growth rate, i.e.

(5.43) \begin{align} \frac {1}{\sigma } (\gamma + i\, \omega )\, u + \frac {\partial _{x} p - v}{{\textit{E}}} &= (\partial _{x}^{2} + \partial _{z}^{2} - k^{2})\, u, \end{align}

(5.44) \begin{align} \frac {1}{\sigma } (\gamma + i\, \omega ) \, v + \frac {ik\, p + u}{{\textit{E}}} &= (\partial _{x}^{2} + \partial _{z}^{2} - k^{2})\, v, \end{align}

(5.45) \begin{align} \frac {1}{\sigma } (\gamma + i\, \omega )\, w + \frac {\partial _{z} p - \theta }{{\textit{E}}} &= (\partial _{x}^{2} + \partial _{z}^{2} - k^{2})\, w, \end{align}

(5.46) \begin{align} \partial _{x} u + ik\, v + \partial _{z} w &= 0, \end{align}

(5.47) \begin{align} (\gamma + i\, \omega )\, \theta - R \, w &= (\partial _{x}^{2} + \partial _{z}^{2} - k^{2})\, \theta , \end{align}

Nominally, the system has no-slip boundary conditions on all walls. For the top and sides of the domain,

(5.48) \begin{eqnarray} u = v = w = 0 \quad \text{at}\quad \; z = 1, \quad \text{and} \quad x = 0,\, \Gamma . \end{eqnarray}

However, vertical reflection symmetry, $z \to 1-z$ , implies we can save computational cost by restricting to only the upper-half domain with mid-plane conditions,

(5.49) \begin{eqnarray} u = v = p = \partial _{z} \theta = 0 \quad \text{at}\quad \; z=\frac {1}{2}. \end{eqnarray}

The perfectly insulating thermal boundary condition at $x=\Gamma$ implies $\partial _{x} \theta = 0$ . We impose perfectly conducting conditions at all other walls, $\theta = 0$ at $x=0$ and $z=1$ . The conducting condition at $x=0$ isolates the dynamics only to the $x=\Gamma$ boundary with quasi-exponential decay into the interior.

We use Dedalus to implement a continuous Galerkin spectral-element solver for (5.43)–(5.49) with a $3 \times 3$ array of elements scaled to capture the extreme boundary-layer behaviour for ${\textit{E}} \ll 1$ . Each subdomain uses a tensor-product discretisation with between 20 and 120 polynomial modes in each direction, chosen so that the solutions in each element are resolved to a truncation error of approximately $ 10^{-6}$ .

We find the minimum onset value for $R$ as a function of $k$ at fixed $E$ by solving the simultaneous equations for the growth rate,

(5.50) \begin{eqnarray} \gamma (k,R) = \partial _{k} \gamma (k,R) = 0. \end{eqnarray}

We solve these equations numerically by computing $\gamma (k,R)$ on a $3\times 3$ grid of $(k,R)$ input values in the vicinity of the leading-order solution and interpolating with a local polynomial approximation of the form

(5.51) \begin{eqnarray} \gamma (k,R) \approx \gamma _{0,0} + \gamma _{1,0}\, k + \gamma _{2,0}\, k^{2} + \gamma _{0,1}\, R. \end{eqnarray}

We update the approximate solution to (5.50) from the polynomial fit via

(5.52) \begin{eqnarray} k \approx - \frac {\gamma _{1,0}}{2 \gamma _{2,0}}, \quad R \approx \frac {\gamma _{1,0}^{2}}{4 \gamma _{0,1} \gamma _{2,0}} -\frac {\gamma _{0,0}}{\gamma _{0,1}}. \end{eqnarray}

The true solution need not lie exactly within the original $3\times 3$ grid; the polynomial approximation can extrapolate slightly outside the fitting range. We perform the procedure with both Chebyshev and Legendre polynomial discretisations to estimate the numerical uncertainty in the critical parameters, which are in the fourth digit or smaller. The asymptotic results are independent of $\sigma$ . We set the Prandtl number $\sigma = 1$ throughout. We use both $\Gamma = 2, 4$ over the whole range of ${\textit{E}}$ . The results agree to within the other uncertainties of the calculation.

Table 1. Direct spectral element calculations of the stability threshold at various Ekman numbers. Values of $R$ , $k$ and $\omega$ have estimated respective uncertainties of $4\times10^{-4}$ , $3\times10^{-4}$ and $3\times10^{-3}$ , with asymptotic analytical results reported to the same number of digits.

We scan over a range of powers of Ekman number with $- \log _{10} {\textit{E}} \ge 6$ . These go to 11. Table 1 reports the onset critical parameters, $(R, k, \omega )$ , from the scans over ${\textit{E}}$ . We perform five-term fits for the critical parameters in powers of $\varepsilon = {\textit{E}}^{\,1/3}$ as

(5.53) \begin{align} R(\varepsilon ) & \approx R_{0} + R_{1/2} \, \varepsilon ^{1/2} + R_{1}\, \varepsilon + R_{3/2} \, \varepsilon ^{3/2} + R_{2}\, \varepsilon ^{2}, \end{align}

(5.54) \begin{align} k(\varepsilon ) & \approx k_{0} + k_{1/2} \, \varepsilon ^{1/2} + k_{1}\, \varepsilon + k_{3/2} \, \varepsilon ^{3/2} + k_{2}\, \varepsilon ^{2}, \end{align}

(5.55) \begin{align} \omega (\varepsilon ) & \approx \omega _{0} + \omega _{1/2} \, \varepsilon ^{1/2} + \omega _{1}\, \varepsilon + \omega _{3/2} \, \varepsilon ^{3/2} + \omega _{2}\, \varepsilon ^{2}. \end{align}

Figure 4. Plots of critical parameters versus Ekman number for $6 \le -\log _{10} {\textit{E}} \le 11$ from the spectral-element calculations of (5.43)–(5.49). (a) Scaled Rayleigh number, $R$ , (b) tangential wave number, $k$ , and (c) precession frequency, $\omega$ . Table 2 gives the coefficients of the best-fit expansions in powers of $E^{1/6}$ . The orange dashed lines show the asymptotic predictions in (5.40)–(5.42).

Figure 4 shows plots over the values covered along with the respective best-fit curves. Table 2 shows the best-fit parameters. The $R_{0}, \,R_{1/2}$ , $k_{0}, \,k_{1/2}$ and $\omega _{0}, \,\omega _{1/2}$ all agree with the predictions in § 5.2 to the number of reported digits.

Table 2. Best-fit coefficients for critical parameters as a function of powers of $\varepsilon = {\textit{E}}^{\,1/3}$ , computed using spectral-element calculations of (5.43)–(5.49) and plotted in figure 4.

5.4. Finite cylinder

In the finite cylinder case, we use polar coordinates, $(r,\phi ,z)$ , where the radius $0 \le r \le \Gamma / 2$ and azimuth angle $0 \le \phi \lt 2 \pi$ , where $\Gamma$ denotes the non-dimensional diameter or aspect ratio.

In this case, we solve the bulk equation first and substitute the result into the boundary condition. The radial dependence uses modified Bessel functions of the first kind, $I_{m}$ , with integer wavenumber, $m \in \mathbb{Z}$ , for $\beta = \sqrt {\pi ^{2} + i\, \omega }$ ,

(5.56) \begin{eqnarray} \theta (r,\phi ,z) = I_{m}(\beta \, r) \, e^{i (m \phi + \omega t)} \sin (\pi z) + \text{c.c.}. \end{eqnarray}

The dispersion relation follows from the boundary condition applied at the outer radius,

(5.57) \begin{eqnarray} R = \pi \beta \left (1 + \frac {\pi \Gamma }{2 i m} \right ) \frac {I_{m}^{\prime}(\beta \, \Gamma /2)}{I_{m}(\beta \, \Gamma /2)} \quad \text{where}\quad \; \textrm{Im}(R) = 0. \end{eqnarray}

Given $m$ and $\Gamma$ , we solve $\textrm{Im}(R) = 0$ using Newton’s method for real-valued $\omega$ . We define $R_{{c}}$ and $\omega _{{c}}$ as their respective values when minimising $R$ over integer $m$ values.

Figure 5. Critical parameters for a finite cylinder with aspect ratio, $\Gamma$ . In each case, the dashed lines show the value for the semi-infinite channel. The critical wavenumber is typically within one of $m \approx \lfloor \Gamma \,k_{{c}}/2 \rfloor$ , where $k_{{c}}$ is the critical wavenumber from the semi-infinite case. Starting from $m=1$ , for low aspect ratio, each cusp feature in $R_{{c}}$ indicates a unit increase in $m$ . The overall minimum happens for $\Gamma \approx 0.366196$ and $m=1$ with $R \approx 26.3789$ and $\omega \approx 82.9854$ .

Figure 6. Mid-plane linear temperature perturbations from two critical modes at aspect ratios $1$ and $3$ , respectively. The $\Gamma =1$ case prefers $m=3$ angular wavenumber, while $\Gamma = 3$ prefers $m=9$ .

Figure 5 shows the dependence of $R_{{c}}$ and $\omega _{{c}}$ on aspect ratio, $\Gamma$ . The cusp-like features in both plots correspond to discrete changes in the optimal $m$ value. At the lowest aspect ratios, the optimal $m=1$ and scales approximately as $m \sim \lfloor \Gamma \,k_{{c}}/2 \rfloor$ . Compared with the various Ekman-number corrections, the finite cylinder has a much stronger, ${\mathcal{O}}( 1 )$ , influence on the onset Rayleigh number.

Figure 6 shows two critical mid-plane temperature perturbation profiles for respective aspect ratios, $\Gamma = 1, 3$ . We report the solutions from applying Newton’s method to (5.56):

for $\Gamma = 1$ , $R_{{c}} \approx 29.26899820255448$ , $\omega _{{c}} \approx 68.93822823201539$ with $m_{{c}} = 3$ ;

for $\Gamma = 3$ , $R_{{c}} \approx 30.84620681276749$ , $\omega _{{c}} \approx 66.29068835812032$ with $m_{{c}} = 9$ .

We also compare to finite ${\textit{E}}=10^{-6}$ calculations of Zhang et al. (Reference Zhang, Reiter, Shishkina and Ecke2024) with $\Gamma =1/2$ , which find $R_{c}\approx 28$ and $\omega _{c} \approx 73.5$ . For asymptotically low ${\textit{E}}$ , with $\Gamma = 1/2$ , we find $R_{{c}} \approx 27.100728659192075$ , $\omega _{{c}} \approx 73.83743415332742$ with $m_{{c}} = 1$ .

5.4.1. Tall aspect ratio

In the past few years, multiple research groups (e.g. Julien et al. Reference Julien, Aurnou, Calkins, Knobloch, Marti, Stellmach and Vasil2016; de Wit et al. Reference de Wit, Aguirre Guzmán, Madonia, Cheng, Clercx and Kunnen2020; Pandey & Sreenivasan Reference Pandey and Sreenivasan2024) have started using tall-aspect-ratio cylinders for rapidly rotating convection experiments; the reason being rotation’s tendency to produce elongated vertical structures. The linear stability result in (5.57) simplifies considerably as $\Gamma \to 0$ . Assuming $\omega \,\Gamma \sim {\mathcal{O}}(1)$ ,

(5.58) \begin{align} R = \frac {2\pi m}{\Gamma } - i \pi \left (\pi -\frac {\omega \,\Gamma }{4(m+1)}\right ) + {\mathcal{O}}(\Gamma ). \end{align}

Imposing $\textrm{Im}(R) = 0$ and minimising over $m$ implies

(5.59) \begin{align} R_{c}\, \Gamma \, \sim \, 2\pi , \quad \omega _{c} \, \Gamma \, \sim \, 8\pi \quad \text{as} \;\quad \Gamma \, \to \, 0, \end{align}

with $m_{c} = 1$ . Notably, the aspect ratio and scaled Rayleigh number exactly cancel the cylinder depth dependence in the onset condition,

(5.60) \begin{align} \frac {g \alpha \Delta T D}{2\Omega \, \kappa } \, = \, 2\pi , \end{align}

where $D$ is the cylindrical diameter. Furthermore, as $ \Gamma \sim {\textit{E}}^{\,1/3}$ , the critical Rayleigh number scales the same as for bulk onset, also with identical horizontal scales, ${\mathcal{O}}({\textit{E}}^{\,1/3})$ .

5.5. Local baroclinic instability

When assembled in one place, it becomes apparent that the wall-mode system is closely akin to the well-known planetary geostrophy (PG) equations, which typically model global-scale stratified ocean dynamics away from coastlines (Vallis Reference Vallis2017, Ch. 5). The PG system has rich dynamical structure. However, a common puzzle is the lack of self-consistent forcing and eventual rundown (Schonbek & Vallis Reference Schonbek and Vallis1999). Several works have considered coupling to smaller-scale quasi-geostrophic (QG) dynamics as an intriguing alternative (see Grooms, Julien & Fox-Kemper Reference Grooms, Julien and Fox-Kemper2011).

This context clarifies that the wall-mode convection equations are a sideways forced PG-QG system that relies on non-hydrostatic QG (similar to Sprague et al. Reference Sprague, Julien, Knobloch and Werne2006; Julien & Knobloch Reference Julien and Knobloch2007) near lateral boundaries and interacts with a barotropic QG system in the stress-free case. However, we can also study the wall-mode equations without the lateral boundary interactions.

In a highly reduced setting, the combined baroclinic temperature and barotropic vorticity equations have a rudimentary form of baroclinic instability. To demonstrate, we neglect all diffusion and boundary conditions. Hence,

(5.61) \begin{eqnarray} \partial _{t} \partial _{z} \psi + J( \psi , \partial _{z}\psi ) = 0, \quad \partial _{t} \!\left \langle \nabla_{\!\bot }^{2} \psi \right \rangle + \left \langle J(\psi ,\nabla_{\!\bot }^{2} \psi )\right \rangle = 0. \end{eqnarray}

The baroclinic equation has a family of exact nonlinear solutions,

(5.62) \begin{eqnarray} \psi = P_{0}(z) + x\, V_{0}(z) + (c + V_{0}(z))\, A(y + c\, t) + \text{c.c.} , \end{eqnarray}

where $V_{0}(z)$ is an arbitrary $y$ -direction velocity profile (without loss of generality, $\!\left \langle V_{0}\right \rangle = 0$ ), $P_{0}(z)$ is and arbitrary hydrostatic pressure, $A(y+c\,t)$ is an arbitrary complex-valued amplitude function and $c$ is a complex-valued phase speed.

The barotropic equation requires the solvability condition,

(5.63) \begin{eqnarray} \!\big\langle (c + V_{0})^{2} \big\rangle = 0 \quad \implies \quad c = \pm \,i\, \| V_{0} \|. \end{eqnarray}

For example, $V_{0} = z-1/2$ , $c = \pm i / \sqrt {12}$ , or $V_{0} = \cos (\pi z)$ , $c = \pm i/\sqrt {2}$ . When $A(y) = A_{0}\,e^{iky}$ for wavenumber $k$ , then the growth rate is $k \| V_{0} \|$ .

6.1. Set-up for numerical implementation

Our Dedalus implementation corresponds to the following description. We use the same non-dimensional form as the asymptotic analysis with control parameters $R$ , $\sigma$ and $\Gamma$ . We use $p(r,\phi ,z)$ as the primary dynamical variable for the bulk. The respective thermal wind relations follow:

(6.1) \begin{eqnarray} \Theta = \partial _{z} p, \quad u_{\bot } = \mathcal{J} \boldsymbol{\cdot} \boldsymbol{\nabla} _{\! \bot \, } p \quad \text{where}\quad \; \mathcal{J} = \hat {\phi } \otimes \hat {r} - \hat {r} \otimes \hat {\phi }. \end{eqnarray}

We represent $p(r,\phi ,z)$ in a vertical cosine series, $\cos (n \pi z)$ . The temperature perturbations, $\Theta (r,\phi ,z)$ , are naturally a sine series $\sin (n \pi z)$ . For $n \ge 1$ , the baroclinic evolution equation

(6.2) \begin{eqnarray} \partial _{t}\Theta - ( \nabla _{\! \bot \, }^{2} + \partial _{z}^{2}) \Theta = - u_{\bot } \boldsymbol{\cdot} \boldsymbol{\nabla} _{\! \bot \, } \Theta . \end{eqnarray}

We use a second-order accurate multi-stage, mixed implicit-explicit time-stepping scheme (SBDF scheme from Wang & Ruuth Reference Wang and Ruuth2008), with linear terms on the left-hand side treated implicitly and the nonlinear right-hand side treated explicitly. The time-step size follows a Courant–Friedrichs–Lewy criterion with a 0.2 safety factor. Because of the fine radial grid spacing near the walls, normal flow through the outer boundary principally determines the time-step size.

For the outer boundary conditions at $r=\Gamma /2$ , we use $q_{z}(\phi ,z)$ as an additional variable and define $q_{\phi } = -\mathcal{H} [ q_{z} ]$ . The boundary conditions are

(6.3) \begin{align} \boldsymbol{\nabla} _{\!\phi \,} q_{\phi } + \partial _{z} q_{z} - \hat {r}\boldsymbol{\cdot} u_{\bot } &= 0, \end{align}

(6.4) \begin{align} R\, q_{z} - \hat {r}\boldsymbol{\cdot} \boldsymbol{\nabla} _{\! \bot \, } \Theta &= (q_{\phi } \boldsymbol{\nabla} _{\!\phi \,} + q_{z} \partial _{z} ) \Theta , \end{align}

where $\boldsymbol{\nabla} _{\!\phi \,} = \hat {\phi } \boldsymbol{\cdot} \boldsymbol{\nabla} _{\!\bot \,} = r^{-1} \partial _{\phi }$ . As before, all linear terms are time-implicit and nonlinear terms are explicit.

The $n=0$ barotropic pressure mode has different dynamics depending on the velocity boundary conditions on the top and bottom. In both cases, we define the vorticity, $\zeta = \nabla _{\!\bot \,}^{2} p$ .

For stress-free, we evolve the $n=0$ vorticity equation,

(6.5) \begin{eqnarray} \partial _{t} \!\left \langle \zeta \right \rangle - \sigma \nabla_{\!\bot \,}^{2} \!\left \langle \zeta \right \rangle = - \boldsymbol{\nabla} _{\!\bot \,} \boldsymbol{\cdot} \left \langle u_{\bot } \zeta \right \rangle , \end{eqnarray}

with, again, the usual linear–nonlinear implicit–explicit time splitting. For boundary conditions at $r =\Gamma /2$ ,

(6.6) \begin{eqnarray} \!\left \langle p\right \rangle = 0, \quad - 2\sigma \, \hat {r} \boldsymbol{\cdot} \boldsymbol{\nabla} _{\!\bot \,} \!\left \langle p\right \rangle = 3 \!\left \langle q_{\phi } \boldsymbol{\nabla} _{\!\phi \,} q_{\phi } \right \rangle - \left \langle q_{\phi } \boldsymbol{\nabla} _{\!\phi \,} p \right \rangle . \end{eqnarray}

For no-slip, the $n=0$ mode involves a bit of subtlety. The Ekman boundary condition directly mixes barotropic and baroclinic variables, and parity mixes sine and cosine series. Vasil, Brummell & Julien (Reference Vasil, Brummell and Julien2008a , Reference Vasil, Brummell and Julien b ) discusses issues like this detail. The problem is that we need to represent expressions like (3.92) spectrally using sine/cosine series. We fix the barotropic pressure via

(6.7) \begin{eqnarray} \nabla_{\!\bot \,}^{2} \!\left \langle p\right \rangle = \nabla_{\!\bot \,}^{2} \!\left \langle K_{\!N} \Theta \right \rangle , \end{eqnarray}

with the boundary condition $\!\langle p\rangle =0$ at $r = \Gamma /2$ . We use the finite spectral representation of $1/2 -z$ ,

(6.8) \begin{eqnarray} K_{\!N}(z) = \sum _{n=1}^{N/2} \frac {\sin (2\pi n z)}{n\pi } \quad \to \quad \frac {1}{2} -z \quad \text{as} \;\quad N \to \infty . \end{eqnarray}

A sufficiently large resolution parameter, $N$ , captures the integral exactly. That is,

(6.9) \begin{eqnarray} \int _{0}^{1} \! K_{N}(z) \, \sin (n \pi z) \,\text{d}{z} = \int _{0}^{1}\! \left (\frac {1}{2}-z \right )\,\sin (n \pi z) \,\text{d}{z} \quad \textrm {(exactly)} \end{eqnarray}

for all $n \le N$ . The advantage is that $K_{N}$ is a natural sine series having a natural inner product with $\Theta$ .

6.2. Weakly nonlinear theory

Weakly nonlinear theory (see, e.g. Hoyle Reference Hoyle2006) is a useful diagnostic for validating calculations slightly above onset. We define the perturbation parameter,

(6.10) \begin{align} \delta = \frac {R}{R_{c}} - 1. \end{align}

As $\delta \to 0$ , we expect pressure perturbations of the form,

(6.11) \begin{align} p &= A(t)\, P_{1,1}(r) \,e^{i(m\phi + \omega _{c} t)} \cos (\pi z) \, \nonumber \\ &\quad + \left ( |A(t)|^{2} \,P_{2,0}(r) + A(t)^{2}\, P_{2,2}(r) \,e^{2i(m\phi + \omega _{c} t)} \right )(\cos (2\pi z) - 1) \, \nonumber \\ &\quad + \left ( |A(t)|^{2} P_{2,0}(a) + A(t)^{2}\,P_{2,0}(a) \left(\frac {r}{a}\right)^{2m} e^{2i(m\phi + \omega _{c} t)} \right ) \, + \text{c.c.} + {\mathcal{O}}\left(\delta ^{3/2}\right), \end{align}

with a complex-valued amplitude varying on a slow time scale

(6.12) \begin{align} A(t) = {\mathcal{O}}(\delta ^{1/2}), \quad A'(t) = {\mathcal{O}}(\delta ^{3/2}). \end{align}

The pressure ansatz contains the leading-order convective modes derived from linear theory and the next-order convective feedback. At a given aspect ratio, we pick the most unstable wavenumber, $m$ , defining the non-dimensional radius, $a=\Gamma /2$ . The radial functions are normalised modified Bessel functions,

(6.13) \begin{align} P_{n,j}(r) \, = \, \frac {I_{m}(r\beta _{n,j})}{ C_{n,j} }, \quad \beta _{n,j} = \sqrt {(n\pi )^{2} + j\, \sqrt {-1} \, \omega }, \end{align}

with the coefficients,

(6.14) \begin{align} C_{1,1} &\, = \, I_{m}(a\beta _{1,1}), \end{align}

(6.15) \begin{align} C_{2,0} &\, = \, \frac {2 \left (\pi ^{2} a^{2}+m^{2}\right )\beta _{2,0} I_{m}^{\prime}(a\beta _{2,0})}{\pi m^{2}}, \end{align}

(6.16) \begin{align} C_{2,2} &\, = \, \frac {4 \pi a (m - i \pi a) \beta _{2,2} I_{m}^{\prime}(a\beta _{2,2}) -2 a m R_{c}I_{m}(a\beta _{2,2}) }{\pi m (\pi a+i m)}, \end{align}

with $P_{n,-j}(r)= P^{*}_{n,j}(r)$ . The normalisations ensure the solution satisfies all horizontal boundary conditions, along with the no-slip barotropic condition (6.7); the stress-free barotropic condition would be considerably more complicated.

The weakly nonlinear amplitude satisfies a complex Ginzburg–Landau equation,

(6.17) \begin{align} A'(t) = (c_{0}\,\delta - c_{1}\, |A(t)|^{2})\, A(t), \end{align}

where

(6.18) \begin{align} c_{0} = & \frac {2 i \pi m \left (\pi ^{2}+i \omega \right ) (\pi a+i m) R_{c}}{a^{2} m^{2} R_{c}^{2}+\pi ^{2} (\pi a+i m)^{2} \left (m^{2}+a^{2} \left (\pi ^{2}+i \omega \right )\right )}, \end{align}

(6.19) \begin{align} c_{1} = & \frac {\frac {4 \pi m a}{m-i \pi a}P_{2,0}(a)\!-\!2i m P_{2,2}(a) \!+\! m \int _{0}^{a} ( (4 \Phi _{2,2}\! -\! 10 P_{2,2}) \textrm{Im}(P_{1,1} P_{1,-1}^{\prime}) \!- \!5i P_{1,1}^{2}P_{1,1}^{\prime} )\,\text{d}{r} }{2\int _{0}^{a} P_{1,1}^{2} r\,\text{d}{r} } \end{align}

for $\Phi _{2,2} = P_{2,2}(a) (r/a)^{2m}$ . Linear theory is sufficient to determine $c_{0}$ . The nonlinear coefficient, $c_{1}$ , requires a messy perturbation series that closes at ${\mathcal{O}}(\delta ^{3/2})$ . There is no spatial modulation because of the discrete critical wavenumber, $m$ .

The amplitude in (6.17) converges to a travelling wave,

(6.20) \begin{align} A(t) \, \to \, \mathcal{A} \,e^{ i \Omega \, t}, \end{align}

with time-independent parameters,

(6.21) \begin{align} |\mathcal{A}|^{2} \, = \, \frac {\textrm{Re}(c_{0}) \,}{\textrm{Re} (c_{1}) } \, \delta ,\quad \Omega \, = \, \frac {\textrm{Im}(c_{0})\, \textrm{Re} (c_{1}) - \textrm{Im}(c_{1})\, \textrm{Re}(c_{0}) }{\textrm{Re} (c_{1}) } \,\delta . \end{align}

The weakly nonlinear solution also gives the Nusselt number near onset. In particular,

(6.22) \begin{align} {\textit{Nu}} - 1 \, = \, - \frac {\partial _{z} \Theta |_{z=0}} {R} \, = \, K\, \delta , \end{align}

where we define the convective conductivity,

(6.23) \begin{align} K = \left . \frac {\partial \log {\textit{Nu}}}{\partial \log R}\right |_{R=R_{c}} \, = \, \frac {2 \pi m^2}{a \left (\pi ^2 a^2+m^2\right )} \frac {\textrm{Re}(c_{0}) \,}{\textrm{Re} (c_{1}) } . \end{align}

As the aspect ratio $\Gamma \to 0$ , the calculation simplifies considerably and $K \to 12/7$ . As the aspect ratio, $\Gamma \to \infty$ , we can use a derivation along a flat Cartesian wall such that

(6.24) \begin{align} K \, \sim \, \frac {0.7018222861264536}{\Gamma } \quad \text{as} \;\quad \Gamma \, \to \, \infty . \end{align}

Like the conductivity, we find similar behaviour for the nonlinear frequency, $\omega = \omega _{c} + \Omega$ . We define the convective Doppler shift,

(6.25) \begin{align} \Upsilon = \left . \frac {\partial \log \omega }{\partial \log R}\right |_{R=R_{c}}. \end{align}

As $\Gamma \to 0$ , $\Upsilon \to 173/252$ . Likewise,

(6.26) \begin{align} \Upsilon \to 1.5336344674049933 \quad \text{as} \;\quad \Gamma \to \infty . \end{align}

Figure 7 shows plots of $K$ and $\Upsilon$ as a function of aspect ratio.

Figure 7. Nusselt number and relative frequency slopes as a function of aspect ratio in a finite cylinder. For the Nusselt number, the horizontal dashed grey line indicates the analytical value of $12/7$ as $\Gamma \to 0$ and the dashed grey curve shows the asymptotic scaling $\sim\! 0.701822/\Gamma$ as $\Gamma \to \infty$ derived from a flat Cartesian wall bounding an infinite bulk interior. For the relative frequency slope, the dashed grey curves indicate the analytical value $173/252$ as $\Gamma \to 0$ and the Cartesian result $\sim \!1.533634$ as $\Gamma \to \infty$ . The jagged behaviour in both plots results from jumps in the critical mode number, $m_{c}$ , as a function of aspect ratio.

6.3. Fully nonlinear results

We first use the methods described previously to simulate the reduced wall-mode equations with no-slip boundary conditions, $R=2R_{{c}}$ , and $\Gamma =4/5$ . For this choice of $\Gamma$ , the most unstable mode at $R=R_{{c}}$ has an azimuthal wavenumber $m=2$ (figure 5). We initialise the simulation with random low-amplitude noise, which triggers the wall-mode instability. At $R=2R_{{c}}$ , the modes $m=2,3,4,5,6$ all have positive growth rates, and $m=3,4$ have the largest growth rates with $\gamma /\omega _{{c}} \approx 2.3$ and $\omega /\omega _{{c}} \approx 3.1, 2.7$ , respectively.

The instability saturates nonlinearly, producing characteristic patterns which propagate retrograde. In figure 8, we plot the temperature near the wall and at mid-height as a function of $\phi$ and $t$ . The downward propagation of the constant temperature stripes shows the retrograde motion. Initially, the dominant mode has $m=3$ , but a coarsening event at $t\approx 0.075$ results in a steady $m=2$ pattern that persists for the rest of the simulation. Others have also found coarsening dynamics in similar systems (Ecke et al. Reference Ecke, Zhong and Knobloch1992; Plaut Reference Plaut2003; Scheel et al. Reference Scheel, Paul, Cross and Fischer2003; Choi et al. Reference Choi, Prasad, Camassa and Ecke2004; Lopez et al. Reference Lopez, Marques, Mercader and Batiste2007; Favier & Knobloch Reference Favier and Knobloch2020). In particular, Ning & Ecke (Reference Ning and Ecke1993), Liu & Ecke (Reference Liu and Ecke1997, Reference Liu and Ecke1999) study the coarsening dynamics for wall modes in great detail. They find that increasing the Rayleigh number slowly (as we do later) maintains $m$ , but a rapid change in Rayleigh (e.g. starting a simulation with $R/R_{{c}}=2$ from noise) usually results in a decrease in $m$ . We chose $\Gamma =4/5$ because the simulations yield an $m=2$ pattern, implying more azimuthal resolution elements per wavelength than simulations with the same resolution but higher $m$ .

Figure 8. Temperature space–time propagation diagram as a function of azimuth, $\phi$ , and time, $t$ . The pattern shows the $R/R_{{c}} = 2$ , no-slip case. The simulation starts with low-amplitude noise. The modes $m=2,3,4,5,6$ are all unstable, with $m=3, 4$ having the largest growth rates. The pattern starts with approximately three full wavelengths, undergoes a coarsening defect and settles into a consistent two-wave pattern.

Having found a steady propagating solution for no-slip boundary conditions and $R=2R_{{c}}$ , we bootstrap simulations to progressively higher and lower $R$ , and simulations with stress-free boundary conditions. We run each simulation until it reaches a statistically steady state. Starting from the saturated state with a similar $R$ means that our subsequent simulations run much shorter than figure 8. In each case, we find the $m=2$ structure persists. Table 3 lists all simulations described in this paper and their spatial resolution.

Table 3. Summary of all simulation parameters. In every case, for $\Gamma = 4/5$ , the critical scaled Rayleigh number, $R_{{c}} = 28.237851887421087$ , and the critical frequency, $\omega _{{c}} = 69.02735982770973$ . For the weakly nonlinear parameters, $K = 0.8028404111646181$ , $\Upsilon = 1.0933189545586768$ . By convention, $\omega \gt 0$ corresponds to retrograde propagation, consistent with the linear instability calculation. For top/bottom boundary conditions, NS is no-slip, and SF is stress-free. While $R,\Gamma ,\sigma$ are input control parameters, figure 12 shows the output parameters, $\omega$ and ${\textit{Nu}}$ . Note that the Prandtl number, $\sigma$ , drops out of the system in the no-slip case.

Figure 9 illustrates the sidewall temperature pattern for different boundary conditions and $R$ . In all cases, alternating hot/cold spots exist near the top/bottom of the domain. These patterns extend towards the mid-plane of the cylinder as $\phi$ increases past the location of the hot/cold spot. The height of these patterns is smaller for no-slip boundary conditions than for stress-free boundary conditions. As $R$ increases, the edges of the hot/cold spots sharpen, producing front-like features that require increasingly higher azimuthal resolution for increased $R$ . The stress-free simulations have sharper, more vertically aligned fronts. For this reason, we could only run simulations up to $R=4R_{{c}}$ for stress-free boundary conditions, while we were able to reach higher values $R=6R_{{c}}$ in the no-slip simulations.

Figure 9. Temperature patterns at the sidewall: (a) the no-slip case and (b) the stress-free case. The rows show $R/R_{{c}} = 2,3,4$ , respectively.

To show the horizontal structure of the modes, we plot the pressure at the top of the cylinder in figure 10. The pressure perturbations extend significantly into the cylinder, although the features become increasingly localised to the walls at $R$ increases. Recall that these visualisations do not directly show the Stewartson layers, which are infinitely thin in this analysis. For stress-free boundary conditions, the pressure fields do not change substantially as $R$ increases, other than becoming more localised near the boundary. However, there are more interesting differences between different $R$ for the no-slip boundary conditions. For low $R$ , the contours of constant pressure slope away from the walls in the counter-clockwise direction, as for the stress-free simulations. However, for higher $R$ , the no-slip simulations have contours of constant pressure that slope away from the walls in the clockwise direction. This change in morphology may be related to the differences in propagation direction between the no-slip and stress-free simulations at higher $R$ .

Figure 10. Pressure patterns at the top: (a) the no-slip case and (b) the stress-free case. The rows show $R/R_{{c}} = 2,3,4$ , respectively.

Figure 11 visualises the momentum fluxes within the Stewartson layers by plotting $\boldsymbol{\nabla} \boldsymbol{\cdot} q$ and $\boldsymbol{\nabla} \times q$ . For the divergence, $\boldsymbol{\nabla} \boldsymbol{\cdot} q = u_{n}$ , the normal velocity at the sidewall. For the curl, $\boldsymbol{\nabla} \times q = \partial _{\ell } q_{z} - \partial _{z}q_{\ell }$ . The temperature at the sidewalls shows sharper azimuthal features in the stress-free simulations than in the no-slip simulations. In particular, $\boldsymbol{\nabla} \times q$ is very sharp and vertically aligned in the stress-free simulation. In the no-slip simulation, $\boldsymbol{\nabla} \boldsymbol{\cdot} q$ and $\boldsymbol{\nabla} \times q$ are much wider and show the same types of sloped features near the top and bottom of the cylinder as the temperature field.

Figure 11. Boundary-layer divergence and curl patterns at the sidewall: (a) the no-slip case and (b) the stress-free case. Both cases correspond to $R/R_{{c}} =4$ .

Next, we focus on heat transport. In the reduced wall-mode equations, the bulk has no vertical velocity. Within the asymptotic approximation, the diffusive heat transport happens entirely in the bulk, while the convective heat transport happens entirely within the Stewartson layer. In polar coordinates,

(6.27) \begin{eqnarray} F_{c} = \frac {1}{|\mathcal{A}|} \int _{0}^{2\pi } \! q_{z} \Theta \,\text{d}{\phi }, \quad F_{d} = -\frac {1}{|\mathcal{A}|}\int _{0}^{2\pi }\!\!\int _{0}^{\frac {\Gamma }{2}} \! \partial _{z} T \,r\,\text{d}{r}\,\text{d}{\phi }, \end{eqnarray}

where $|\mathcal{A}| = \pi \Gamma ^{2}/4 = 4\pi /25 \approx 0.502$ and $\Theta = T - \overline {T}$ . These fluxes sum to the total $F_{t} = F_{c}+F_{d}$ . The dynamical equations (6.2)–(6.4) have the mean linear conduction temperature profile subtracted out, i.e. $- \partial _{z} T_{0} = R$ ; therefore, the Nusselt number ${\textit{Nu}} = F_{t}/R$ .

Figure 12(a) shows ${\textit{Nu}}$ as a function of $R/R_{{c}}$ for both no-slip and stress-free simulations. As $R$ increases from $R_{{c}}$ , we find that ${\textit{Nu}}$ also increases linearly, as ${\textit{Nu}} - 1 \approx 0.807\,(R/R_{c}-1)$ , shown as a dashed line in the plot. Weakly nonlinear theory predicts a coefficient $\approx 0.8028404111646181$ ; see § 6.2 and figure 7. As $R$ increases further, ${\textit{Nu}}$ appears to level off, although we have not been able to run at sufficiently high $R/R_{{c}}$ to probe the asymptotic behaviour as $R$ becomes large.

Figure 12(b) shows the different depth-dependent heat fluxes for the no-slip and stress-free simulations with $R=4R_{{c}}$ . The convective heat flux is zero at the top and bottom boundaries ( $q_z=0$ ), but becomes significant and roughly constant away from the horizontal boundaries. The diffusive heat flux is large and positive in the top and bottom thermal boundary layers, but intriguingly becomes negative near the mid-plane. The (vertical) middle of the cylinder becomes stably stratified due to the wall modes. All of our simulations are thermally equilibrated, where the total heat flux is constant with height.

Figure 12. Diagrams of heat transport. (a) Nusselt number, ${\textit{Nu}}$ , versus the supercriticality, $R/R_{{c}}$ ; by definition ${\textit{Nu}}_{\,{c}} = 1$ at onset. The dashed line shows the approximate linear dependence ${\textit{Nu}}-1 \approx 0.807\, (R/R_{c}-1)$ near onset; weakly nonlinear predicts a slope of $\approx \,$ 0.803 (an $\approx \,$ 0.5 % difference). (b) Depth-dependent diffusive ( $F_{d}$ , blue) and convective fluxes ( $F_{c}$ , red), along with their total ( $F_{t} = F_{c} + F_{d}$ , gold) for simulations with $R=4R_{{c}}$ . The solid lines correspond to the no-slip case, and the dashed lines correspond to stress-free. The independence of $F_{t}$ on depth, $z$ , implies the systems are in statistically steady state with a well-defined ${\textit{Nu}}$ .

Figure 8 shows the instability of the wall modes saturates into nonlinear travelling waves. At onset, we calculate the propagation rate via linear theory; for $\Gamma =4/5$ , $\omega _{{c}} \approx 69$ in the retrograde direction. Previous work has shown the propagation rate varies with $R$ as expected for a supercritical Hopf bifurcation (e.g. Zhong et al. Reference Zhong, Ecke and Steinberg1991; Liu & Ecke Reference Liu and Ecke1999; Favier & Knobloch Reference Favier and Knobloch2020). Figure 13 shows the propagation rate as a function of $R/R_{{c}}$ for no-slip and stress-free simulations. As $R$ increases from $R_{{c}}$ , the propagation rate $\omega$ also increases approximately linearly. At $R=2R_{{c}}$ , the propagation rates are similar for no-slip and stress-free simulations. However, for larger $R$ , the propagation rate increases monotonically for stress-free simulations, whereas the propagation rate decreases for simulations with no-slip boundary conditions. The experiments of Zhong et al. (Reference Zhong, Ecke and Steinberg1993) at ${\textit{E}} \approx 10^{-4}$ similarly found a maximum propagation rate for $R/R_{c} \approx \text{2-3}$ . At $R\approx 5R_{{c}}$ , there is a steady nonlinear solution with no propagation and for higher $R$ , the solution propagates in the prograde direction. This switch in propagation direction may relate to the change in the azimuthal structure of the pressure perturbations shown in figure 10. Horn & Schmid (Reference Horn and Schmid2017), Ravichandran & Wettlaufer (Reference Ravichandran and Wettlaufer2023) and Zhang et al. (Reference Zhang, Reiter, Shishkina and Ecke2024) all find mixtures of both prograde and retrograde propagation, depending on parameter values. Note that the structure of the temperature perturbations at the boundary (e.g. shown up to $R=4R_{{c}}$ in figure 9) does not change significantly as $R$ increases, other than becoming increasingly sharp.

Figure 13. Propagation rate, $\omega / \omega _{{c}}$ versus the supercriticality, $R/R_{{c}}$ . The dashed line shows the approximate linear dependence $\omega /\omega _c-1 \approx 1.05\, (R/R_{c}-1)$ near onset; compared with a predicted slope of approximately 1.09. The rates continue to increase up to $R/R_{{c}} \approx 2$ . At that point, the stress-free case continues the trend, while the no-slip rate declines and eventually becomes prograde beyond $R/R_{{c}} \approx 5$ .

Finally, figure 14 shows the mean zonal flow for no-slip and stress-free simulations with $R=4R_{{c}}$ . At this $R$ , the modes propagate retrograde in both simulations. However, the mean zonal flow at the boundary is positive in the no-slip simulation, while it is (as expected) negative in the stress-free simulation. The change in propagation direction for simulations with no-slip boundary conditions may be related to the development of a prograde mean zonal flow.

Figure 14. Mean zonal flow in meridional planes: (a) no-slip case and (b) stress-free case. Both cases correspond to $R/R_{{c}} =4$ .