Vorticity skewness of finite-amplitude rapidly rotating Rayleigh–Bénard convection

Abstract

Rotating Rayleigh–Bénard convection denotes the convection between a warm plate and a cold plate in a rotating environment. It is a classic model for understanding convective vortices in the atmosphere and ocean. The influence of background rotation on fluid inertia breaks the symmetry between cyclones and anticyclones. Such a symmetry breaking could be represented by vorticity skewness, which still lacks a systematic theory. Rapidly rotating convection with stress-free boundaries and unit Prandtl number is a convenient starting point. The investigation starts from the convective onset stage, where the vortices grow stationarily. Asymptotic analysis shows that the volumetric vorticity skewness $S$ is produced by the interaction between the $n=0,1$ and $n=1,2$ vertical eigenmodes. The $n=0$ (barotropic) mode contributes positively to $S$ mainly by stretching the vertical relative vorticity, an ageostrophic effect. The $n=2$ mode makes a minor negative contribution to $S$ by preferentially intensifying the outflow over the inflow, a non-hydrostatic effect. The theory predicts $S$ to be proportional to the global Rossby number defined with the volumetric standard deviation of vorticity, ${Ro_g}$. The proportional factor does not depend on the Rayleigh and Ekman numbers, agreeing with direct numerical simulations. Then the system enters the equilibrium stage. The stretching of vertical vorticity still contributes to $S$ dominantly. At ${Ro_g}\gtrsim 0.5$, the emergent unsteady flow significantly suppresses the asymmetry between the inflow and outflow strength, and weakens its influence on $S$.

1. Introduction

Fluid motion in a rotating background tends to organize into vortices that spin along the rotating axis. The fluid area with the same or opposite sign of vertical vorticity to the solid-body vorticity is defined as a cyclone or an anticyclone, respectively. The strengths of cyclones and anticyclones are not statistically identical due to the symmetry breaking brought by the influence of background rotation on the fluid inertia (e.g. Polvani et al. 1994; Julien et al. 1996b; Yavneh et al. 1997; Hakim, Snyder & Muraki 2002; Morize, Moisy & Rabaud 2005; Sreenivasan & Davidson 2008; Naso 2015). Background rotation does not induce vorticity asymmetry when it dominates the fluid inertia, i.e. in the quasi-geostrophic limit. This is because cyclones and anticyclones are dominantly produced by stretching and squashing the solid-body vorticity, which are statistically identical processes. Background rotation also does not induce vorticity asymmetry when it is too weak to influence the flow. Thus the asymmetry vanishes at these two ends (Vorobieff & Ecke 2002).

It is well known that the vertical vorticity is stretched more efficiently at the cyclonic region than squashed at the anticyclonic region due to the influence of relative vorticity on the absolute vorticity. This ageostrophic mechanism, briefed as the stretching of vertical relative vorticity, is key for producing vorticity asymmetry in rotating fluids (e.g. Morize et al. 2005); the challenge is quantifying how the asymmetry depends on the solid-body rotation rate. In addition, it remains unclear whether the non-hydrostatic effect, which generates the asymmetry between convergent and divergent flow, could indirectly influence the vorticity asymmetry. This paper explores the vorticity asymmetry of rapidly rotating Rayleigh–Bénard convection, a canonical non-hydrostatic flow near the quasi-geostrophic end.

Rotating Rayleigh–Bénard convection (RRBC) is a prototype model of rotating convection. It is the free convection between a warm lower plate and a cold upper plate in a rotating background (e.g. Bénard 1901; Rayleigh 1916; Chandrasekhar 1953; Nakagawa & Frenzen 1955; Veronis 1959; Boubnov & Golitsyn 1986; Julien et al. 1996b; Ding et al. 2023; Ecke & Shishkina 2023; Anas & Joshi 2024), and has implications for the dynamics of tropical cyclones, tornadoes, open ocean convection, Earth's dynamo, and the interior circulation of gas giants, etc. (Marshall & Schott 1999; Hendricks, Montgomery & Davis 2004; Vasavada & Showman 2005; Aurnou et al. 2015; Horn & Aurnou 2021; Vélez-Pardo & Cronin 2023). Without considering centrifugal acceleration, RRBC could be governed by three independent non-dimensional parameters:

(1.1a–c)\begin{equation} {Ra=}\frac{\beta g\,\Delta T\,H^3}{\nu \kappa }, \quad {E=}\frac{\nu }{fH^2}, \quad{Pr=}\frac{\nu }{\kappa }, \end{equation}

where $\beta$ is the thermal expansion coefficient, $g$ is the gravitational acceleration, $\Delta T$ is the temperature difference between the bottom and top, $H$ is the depth of the fluid, $f$ is the Coriolis parameter (solid-body vorticity) that equals twice the background rotation rate, and $\nu$ and $\kappa$ are the kinematic viscosity and thermal diffusivity, respectively. The Rayleigh number ${Ra}$ depicts the relative strength of the destabilizing effect of buoyancy to the damping effect of viscosity and thermal diffusion. The Ekman number ${E}$ represents the relative strength of viscosity and the rotational effect. The Prandtl number ${Pr}$ is the ratio of kinematic viscosity to thermal diffusivity.

The equilibrium state of RRBC is further classified with the reduced Rayleigh number $\widetilde {{Ra}} \equiv {Ra}\,{{E}}^{{4/3}}$. It represents the extent to which the convective system is supercritical to neutral stability (Chandrasekhar 1961). According to Stellmach et al. (2014), when $\widetilde {{Ra}}\lesssim 25$, the flow is dominated by quasi-steady densely packed columnar vortices that barely interact with each other, called the cellular regime. When $25\lesssim \widetilde {{Ra}}\lesssim 70$, the system is in the convective Taylor column regime. The vortices are packed less densely and more significantly shielded with opposite-sign vorticity, and vortex interaction remains weak. For $\widetilde {{Ra}} \gtrsim 70$, organized vortices break into unsteady plumes and could form barotropic large-scale vortices.

In RRBC, vertical vorticity has been found to skew towards positive generally (Julien et al. 1996a,b, 2012; Vorobieff & Ecke 2002; Kunnen, Geurts & Clercx 2009, 2010b; Favier, Silvers & Proctor 2014; Guervilly, Hughes & Jones 2014; Guzmán et al. 2020; Shi et al. 2020). We call it a ‘cyclonic bias’. Previous studies on the cyclonic bias of RRBC focus on the turbulent regime and have presented three qualitative explanations. First, the convergent flow makes a cyclone more compact, and the divergent flow makes an anticyclone more diluted (Guervilly et al. 2014). This is essentially the stretching of vertical vorticity. Second, the turbulent mixing due to vortex–vortex interaction weakens the convective cell's outflow to produce a diluted anticyclone (Julien et al. 1996a,b; Vorobieff & Ecke 2002; Kunnen, Clercx & Geurts 2010a; Kunnen et al. 2010b; Shi et al. 2020). This effect could be enhanced by Ekman pumping, which intensifies a convective cell's inflow to produce a compact cyclone (e.g. Kunnen, Clercx & Geurts 2006). Third, an anticyclone with vertical relative vorticity lower than $-f$ is susceptible to centrifugal instability (Kunnen et al. 2010b; Favier et al. 2014), but there is no such constraint for cyclonic vorticity. The eddies produced by centrifugal instability might further dilute the anticyclones, so they can also enhance the second factor. These arguments have covered most of the ageostrophic effects and the asymmetry between convergent and divergent flow. Still, a theoretical framework is needed to clarify which mechanism dominates at which stage or regime.

Little progress in theoretical modelling has been made on the vorticity skewness in the weakly nonlinear regime of RRBC, which is the basis for understanding vorticity skewness in the turbulent regime. In the weakly nonlinear regime, the instability is finite-amplitude. Previous studies have diagnosed the vorticity skewness at each horizontal slice (e.g. Julien et al. 1996b; Vorobieff & Ecke 2002; Kunnen et al. 2009), but how it depends quantitatively on the solid-body rotation rate has not been discovered. The pioneering work of Veronis (1959) has presented a comprehensive asymptotic analysis of weakly nonlinear RRBC in the steady state. The vorticity skewness, a nonlinear effect, can be calculated with his second-order special solutions. However, possibly due to the low expectation for a simple result, the analysis seems to have not been done. Our interest in this problem was ignited when plotting the simulated volumetric vorticity skewness $S$ against the global (volumetric) Rossby number ${Ro_g}$, using a series of direct numerical simulations (DNS) with different ${Ra}$ and ${E}$, and fixed ${Pr}=1$. Here, $S$ and ${Ro_g}$ are defined as

(1.2a,b)\begin{equation} S \equiv \frac{ \langle \overline{ {\omega^*_z}^3 } \rangle }{\langle \overline{ {\omega^*_z}^2 } \rangle^{3/2} }, \quad {Ro_{g}} \equiv \frac{\langle \overline{{\omega^*_z}^2} \rangle^{1/2}}{f}, \end{equation}

where $\omega ^*_z$ denotes the (dimensional) vertical vorticity, $\overline {\omega ^*_z}$ denotes its vertical average, and $\langle \omega ^*_z \rangle$ denotes its horizontal average. The volumetric skewness and Rossby number have hardly been used in previous studies of RRBC. The volumetric averaging eliminates many freedoms, showing compactly the vorticity asymmetry of the system. At the convective onset stage, we found $S\propto {Ro_g}$ with a proportional factor approximately 2.5 that hardly depends on ${Ra}$ and ${E}$ in the ${Ro_g}\ll 1$ regime. This inspires us to revisit the finite-amplitude RRBC with asymptotic analysis and explain why $S \propto {Ro_g}$. We start the investigation from the convective onset stage, where the vortices are erect and highly organized. We find that the skewness is contributed positively by the stretching of $\omega ^*_z$ and contributed negatively by vorticity tilting and outflow intensification. Then we extend the theory to the equilibrium stage by accounting for the eddy-induced vertical shear, a stochastic factor that breaks the vertical coherency and terminates the outflow intensification. The flow fields at the two stages are demonstrated in figure 1, as well as in movies 1 and 2 in the supplementary material (available at https://doi.org/10.1017/jfm.2024.571).

Figure 1. Cross-sections of vertical vorticity normalized by $f$ at (a) the convective onset stage and (b) the equilibrium stage. The data of the Ra3 experiment (see table 1) with stress-free boundaries are used ($\widetilde {{Ra}}=23.2$ and ${E}=10^{-4}$). The plots are of $\omega _z^*/f$ at the non-dimensional times (a) $t=8$ and (b) $t=100$. The horizontal plane is at $z=0.2$, and the vertical planes are at $x=1.25$ and $y=1.25$. The simulation domain is doubly periodic, with $0< x<2.5$, $0< y<2.5$ and $0< z<1$. See § 2 for the experimental setting and the non-dimensionalization procedure.

The paper is organized as follows. Section 2 introduces the governing equation and DNS set-up. Section 3 presents the $S$ and ${Ro_g}$ diagnosed from DNS. Section 4 uses vertical mode decomposition and an asymptotic equation set to study $S$ at the convective onset stage. Section 5 extends the theory to the equilibrium stage. Section 6 concludes the research.

2. The governing equations and DNS set-up

This section introduces the variables, the governing equations and the DNS set-up. The dimensional variables (not including constant quantities) are denoted with *. Here, $\boldsymbol {i}, \boldsymbol {j}, \boldsymbol {k}$ are the unit vectors of the Cartesian coordinates, $\boldsymbol {x}^*=(x^*,y^*,z^*)$ is the position, $\boldsymbol {u}^*=(u^*,v^*,w^*)$ is the velocity, $p^*$ is the pressure potential, $T^{*}$ is the perturbation temperature that has subtracted a diffusive-equilibrium linear temperature profile, and $\boldsymbol {\omega }^*=({\omega }_x^*,{\omega }_y^*,{\omega }_z^*)$ is the vorticity, with $\boldsymbol {\omega }^*=\boldsymbol {{\nabla }}^* \times \boldsymbol {u}^*$, where $\boldsymbol {\nabla }^* \equiv \boldsymbol {i}\,\partial /\partial x^* +\boldsymbol {j}\,\partial /\partial y^* + \boldsymbol {k}\, \partial /\partial z^*$ is the gradient operator. We non-dimensionalize the variables in the formulation of Portegies et al. (2008). The convective overturning time scale is $H/W$, where the characteristic vertical velocity $W$ uses the free-fall scaling $W=\sqrt {g\beta \,{\Delta } T\,H}$. The length scale is the domain height $H$. The temperature scale is the temperature difference between the lower and upper plates, ${\Delta }T$. So

(2.1)\begin{equation} \left.\begin{gathered} t^*=t H/W,\quad {(x^*,y^*,z^*)}={(x,y,z)}\,H, \quad \boldsymbol{u}^*=\boldsymbol{u} W, \\ \boldsymbol{\omega}^*=\boldsymbol{\omega}W/H, \quad T^{*}= T\,\Delta T, \quad p^*=p W^2. \end{gathered}\right\} \end{equation}

The convection is between a warm plate at $z=0$ and a cold plate at $z=1$, with a doubly periodic lateral boundary. The flow obeys the incompressible Boussinesq equation:

(2.2)$$\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t}+\left(\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla} \right)\boldsymbol{u}+\frac{1}{{Ro}}\,\boldsymbol{k} \times \boldsymbol{u} =- \boldsymbol{\nabla} p + T \boldsymbol{k}+{\left(\frac{{Pr}}{{Ra}}\right)}^{1/2}{\nabla}^2 \boldsymbol{u}, \end{gather}$$

(2.3) $$\begin{gather}\frac{\partial T }{\partial t}+\left(\boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla}\right) T -w= {\left({Pr\,Ra}\right)}^{-1/2}\,{\nabla}^2 T, \end{gather}$$

(2.4)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}$$

where $\boldsymbol {\nabla }\equiv \boldsymbol {i}\,\partial /\partial x +\boldsymbol {j}\, \partial /\partial y + \boldsymbol {k}\,\partial /\partial z$ is the non-dimensional gradient operator, and ${Ro} \equiv {E}({Ra/Pr})^{1/2}= W/(\,fH)$ is the convective Rossby number that measures the relative strength of thermal forcing to the rotational effect (Julien et al. 1996b). Note that both ${Ro}$ and $\widetilde {{Ra}}\equiv {Ra}\,{E}^{4/3}$ are combinations of ${Ra}$ and ${E}$, but their physical interpretations are different: ${Ro}$ measures the deviation from geostrophic balance at the equilibrium state, and $\widetilde {{Ra}}$ measures the deviation from neutral stability ($\widetilde {{Ra}}\approx 8.7$; Chandrasekhar 1953). For the weakly nonlinear and rapidly rotating RRBC, there is ${Ro} \ll 1$ and $\widetilde {{Ra}} \lesssim O(10^1)$.

The temperature boundary condition is Dirichlet:

(2.5a,b)\begin{equation} T |_{z=0}=0,\quad T |_{z=1}=0. \end{equation}

The impermeable velocity boundary condition is

(2.6a,b)\begin{equation} \boldsymbol{k} \boldsymbol{\cdot}\boldsymbol{u}|_{z=0}=0,\quad \boldsymbol{k} \boldsymbol{\cdot} \boldsymbol{u}|_{z=1} =0. \end{equation}

For the tangential velocity, we study only the stress-free boundary condition:

(2.7a,b)\begin{equation} \boldsymbol{k} \times \left.\frac{\partial \boldsymbol{u}}{\partial z}\right|_{z=0}=0, \quad \boldsymbol{k} \times \left.\frac{\partial \boldsymbol{u}}{\partial z}\right|_{z=1}=0. \end{equation}

The DNS are performed with the Boussinesq solver of Cloud Model 1, version 19.8 (CM1; Bryan & Fritsch 2002). See Appendix A for the detailed numerical setting. The simulation is run in a $[0,2.5]\times [0,2.5]\times [0,1]$ horizontally doubly periodic domain. The initial condition is a spatially uncorrelated random noise on $T$ in the whole domain, obeying a uniform distribution between $-0.25$ and $0.25$. This noise generation method is the default setting of the configured ‘Rayleigh–Bénard convection’ set-up in CM1. The perturbation amplitude uses the default value, which is relatively large but does not cause numerical instability in our simulations. We fix ${Pr}=1$, which is close to the ${Pr}\approx 0.7$ value of air and is mathematically simple (e.g. Vélez-Pardo & Cronin 2023). Because our $Pr$ is above 0.68, the primary instability is stationary, with vortices growing at the same location (Chandrasekhar 1953).

We define the convective onset stage as the stage where the vortices have not been sufficiently deformed and tilted by mutual advection or by any secondary instability (e.g. Küppers & Lortz 1969; Carton 1992). Practically, the end of the convective onset stage is set as the overshooting peak of the time series of ${Ro_g}$, which marks the initiation of vortex interaction that breaks down the stationary pattern.

We perform two groups of experiments, with eight simulations in each group.

1. (i) Changing ${Ra}$ and using a fixed ${E}=10^{-4}$.

2. (ii) Changing ${E}$ and using a fixed ${Ra}=2.5\times 10^{6}$.

The experimental parameters are listed in table 1 and plotted in figure 2. We explore the parameter space around ${Ra}=2.5\times 10^6$ and ${E}=10^{-4}$, because this point has been investigated carefully by Portegies et al. (2008), and detailed background information is available. The critical Rayleigh number ${Ra_c}$ for the stress-free boundary case has an asymptotic expression: ${Ra_c} = 8.6956 (1-1.108{E}^{1/6} + 0.1533 {E}^{1/3}) {E^{-4/3}}$, derived by Homsy & Hudson (1971). The ratio ${Ra}/{Ra_c}$, which is a more accurate measure of supercriticality than $\widetilde {{Ra}}$, ranges from 1.52 to 16.30 in the experiments. The equilibrium states of our experiments are mainly in the cellular regime ($\widetilde {{Ra}}\lesssim 25$; Stellmach et al. 2014).

Table 1. Values of ${Ra}$, ${E}$, ${Ro \equiv E(}{{Ra/Pr)}}^{{1/2}}$, $\widetilde {{Ra}} \equiv {Ra}\,{{E}}^{{4/3}}$, and ${Ra}/{Ra_c}$ in the 16 numerical experiments. The ${Pr}$ value is fixed at 1. Note that Ra7 and E7 are identical.

Figure 2. The ${E}$–${Ra}$ parameter space investigated in this paper. The blue crosses denote experiments Ra1–Ra8. The red crosses denote experiments E1–E8. The black line denotes the critical Rayleigh number ${Ra_c}$ derived by Homsy & Hudson (1971).

The time interval of data output is $0.5$. The total simulation length is 120, which goes through the convective onset stage and finally reaches a statistically equilibrium stage. Before calculating any quantity, the data are interpolated from the vertically stretched mesh to a vertically uniform mesh with 300 layers using the cubic spline function. The motivation for the interpolation is to facilitate the calculation of numerical integral and finite difference. As readers will see, $S$ calculated with the interpolated data converges to zero in the ${Ro_g} \to 0$ limit, so the interpolation error should be negligible.

3. Simulation results

Unlike previous studies that investigate the vorticity skewness at different heights (Julien et al. 1996b; Kunnen et al. 2009), we investigate the contribution from different vertical eigenmodes. Vertical mode decomposition is a useful technique in studying linear and weakly nonlinear waves in a vertically confined domain (e.g. Vallis 2017). To our knowledge, it has not been applied to analyse the vorticity skewness in RRBC. For RRBC with stress-free boundary conditions, the vertical eigenfunction is the trigonometric function:

(3.1a,b)\begin{equation} \sin(n{\rm \pi} z)\quad\text{and}\quad\cos(n{\rm \pi} z),\quad n=0,1,2,\ldots. \end{equation}

The $n=0$ mode is called the barotropic mode, and $n\geq 1$ modes are called baroclinic modes. Variables with the Dirichlet and Neumann boundary conditions take the $\sin (n{\rm \pi} z)$ and $\cos (n{\rm \pi} z)$ shapes, respectively. Thus $u$, $v$, $\omega _z$, $p$ have a vertical structure of $\cos ( n {\rm \pi}z )$, and $w$, $T$, $\boldsymbol {\omega }_h$ (horizontal vorticity vector) have a vertical structure of $\sin ( n {\rm \pi}z )$.

In this section, we do not distinguish between baroclinic modes and simply decompose $\omega _z$ into the barotropic mode and a bulk baroclinic mode. The barotropic mode is simply the vertical average of $\omega _z$, defined as $\overline {\omega _z}$. The baroclinic mode is the vertical anomaly, defined as $\omega '_z$. Because only the baroclinic modes can be linearly unstable, $\overline {\omega _z}$ must be nonlinearly generated by baroclinic modes, and there should be $O(\overline {\omega _z}) \ll O(\omega '_z)$ in the weakly nonlinear regime. Substituting $\omega _z = \overline {\omega _z} + \omega '_z$ into (1.2a–c), we decompose $S$ into the barotropic contribution, the baroclinic contribution and the barotropic–baroclinic contribution:

(3.2) \begin{align} S &= \overbrace{ \frac{ \langle \overline{ \overline{\omega_z}^3 } \rangle }{\langle \overline{\omega^2_z} \rangle^{3/2} } }^{\textit{barotropic}} + \overbrace{ \frac{ \langle \overline{ {\omega'_z}^3 } \rangle }{\langle \overline{\omega^2_z} \rangle^{3/2} } }^{\textit{baroclinic}} + \overbrace{ 3\,\frac{ \langle \overline{ \overline{\omega_z}^2 {\omega'_z} } \rangle }{\langle \overline{\omega^2_z} \rangle^{3/2} } + 3\,\frac{ \langle \overline{ \overline{\omega_z} {\omega'_z}^2 } \rangle }{\langle \overline{\omega^2_z} \rangle^{3/2} } }^{\textit{barotropic--baroclinic}} \nonumber\\ &\approx \frac{ \langle \overline{ {\omega'_z}^3 } \rangle }{\langle \overline{\omega^2_z} \rangle^{3/2} } + 3\,\frac{ \langle \overline{ \overline{\omega_z} {\omega'_z}^2 } \rangle }{\langle \overline{ \omega^2_z } \rangle^{3/2} }. \end{align}

The purely barotropic contribution ${ \langle \overline { \overline {\omega _z}^3 } \rangle }/{\langle \overline {\omega ^2_z} \rangle ^{3/2}}$ is negligible because $O(\overline {\omega _z}) \,{\ll}\, O(\omega '_z)$. The $3({ \langle \overline { \overline {\omega _z}^2 {\omega '_z} } \rangle }/{\langle \overline {\omega ^2_z} \rangle ^{3/2} })$ term equals zero because $3 ({ \langle \overline { \overline {\omega _z}^2 {\omega '_z} } \rangle }/{\langle \overline {\omega ^2_z} \rangle ^{3/2} }) = 3( { \langle \overline {\omega _z}^2 \overline { {\omega '_z} } \rangle }/ \langle \overline {\omega ^2_z} \rangle ^{3/2}) =0$, where we have used $\overline {\omega '_z}=0$. This simplification is validated with DNS (figures 3 and 4). The only two terms remaining are the purely baroclinic term ${ \langle \overline { {\omega '_z}^3 } \rangle }/{\langle \overline {\omega ^2_z} \rangle ^{3/2} }$ and the barotropic–baroclinic term $3 ({ \langle \overline { \overline {\omega _z} {\omega '_z}^2 } \rangle }/{\langle \overline { \omega ^2_z } \rangle ^{3/2} })$. Figures 3(e–h) and 4(e–h) show that at the convective onset stage, the purely baroclinic term is negative, and the barotropic–baroclinic term is positive. The magnitude of the former is approximately half of the latter, except for the Ra1 and E1 experiments that have relatively large $\widetilde {{Ra}}$ and ${Ro}$.

Figure 3. (a–d) The time evolution of the vorticity skewness $S \equiv { \langle \overline { \omega ^3_z } \rangle }/{\langle \overline { \omega ^2_z } \rangle ^{3/2} }$ (black line) and its decomposed parts. The blue line shows ${ \langle \overline { \overline {\omega }^3_z } \rangle }/{\langle \overline {\omega ^2_z} \rangle ^{3/2} }$, the red line shows ${ \langle \overline { {\omega '_z}^3 } \rangle }/{\langle \overline {\omega ^2_z} \rangle ^{3/2} }$, the yellow line shows $3({ \langle \overline { \overline {\omega _z}^2 {\omega '_z} } \rangle }/{\langle \overline {\omega ^2_z} \rangle ^{3/2} })$, and the purple line shows $3 ({ \langle \overline { \overline {\omega _z} {\omega '_z}^2 } \rangle }/{\langle \overline {\omega ^2_z} \rangle ^{3/2} })$. The blue shading marks the convective onset stage, which ends at the time of maximum ${Ro_g}$. The red shading marks the equilibrium stage, which is between $t=90$ and $t=120$, and will be studied in § 5. The four plots are for experiments (a) Ra1, (b) Ra3, (c) Ra5, and (d) Ra8. (e–h) Comparison of the purely baroclinic term ${ \langle \overline { {\omega '_z}^3 } \rangle }/{\langle \overline {\omega ^2_z} \rangle ^{3/2} }$ and the barotropic– baroclinic term $3 ({ \langle \overline { \overline {\omega _z} {\omega '_z}^2 } \rangle }/{\langle \overline {\omega ^2_z} \rangle ^{3/2} })$, with the colour of the dots denoting time ($0< t<20$, from blue to yellow). The time interval between two dots is 0.5. For (e), only a small portion of data points lie in the scope of plotting. The solid black line is the $2:-1$ reference line. Here, (a,e) $\widetilde {Ra}=58.0$, (b,f) $\widetilde {Ra}=23.2$, (c,g) $\widetilde {Ra}=14.5$, and (d,h) $\widetilde {Ra}=10.5$.

Figure 4. The same as figure 3, but for experiments E1, E3, E5 and E8 that change ${E}$ and fix ${Ra}$. For (e,f), only a small portion of data points lie in the scope of plotting. Here, (a,e) $E=5.00\times 10^{-4}$, (b,f) $E=2.50\times 10^{-4}$, (c,g) $E=1.67\times 10^{-4}$, and (d,h) $E=9.09\times 10^{-5}$.

Equation (3.2) indicates that the skewness originates from the correlation between the baroclinic or barotropic vorticity with $\omega '^2_z$. The region with a high $\omega '^2_z$ is essentially the rich-vorticity region. We define the lines crossing the centre of the rich-vorticity region as the vortex axes. At the convective onset stage, they coincide with the axes of strongest vertical motions. Therefore, the vorticity structure along the vortex axis is crucial for modelling the skewness.

Figures 5 and 6 show $S \propto {Ro}_{{g}}$ for almost all simulations, with proportional factor approximately 2.5. The exceptions are the Ra1 and E1 experiments, where the ${Ro}\ll 1$ and $\widetilde {{Ra}} \lesssim O(10^1)$ conditions are violated. Note that ${Ro_g}$ and ${Ro}$ have different physical meanings: ${Ro_g}$ is a diagnostic quantity that represents the instantaneous vorticity amplitude, while ${Ro}$ is an a priori estimation of ${Ro_g}$ at the equilibrium state. Section 4 focuses on the convective onset stage. We use asymptotic analysis and vertical mode decomposition to derive a Galerkin model, which proves $S \propto {Ro_g}$ and explains why the barotropic–baroclinic and purely baroclinic terms make opposite contributions to $S$. They are the basis for understanding the more complicated equilibrium stage in § 5.

Figure 5. The dots show the dependence of the volumetric vorticity skewness $S$ on the global Rossby number ${Ro_{g}}$ during $0< t<20$. The data sampling time interval is 0.5. (a–h) Plots for experiments Ra1–Ra8 with a fixed ${E}=10^{-4}$. The colour of the dots denotes time ($0< t<20$, from blue to yellow). The solid black lines show the $S = 2.5 {Ro_g}$ fitting relation. Here, (a) $\widetilde {Ra}=58.0$, (b) $\widetilde {Ra}=38.7$, (c) $\widetilde {Ra}=23.2$, (d) $\widetilde {Ra}=16.6$, (e) $\widetilde {Ra}=14.5$, (f) $\widetilde {Ra}=12.2$, (g) $\widetilde {Ra}=11.6$, and (h) $\widetilde {Ra}=10.5$.

Figure 6. The same as figure 6, but for experiments E1–E8 that change ${E}$ and fix ${Ra}$: (a) $E=5.00\times 10^{-4}$, (b) $E=3.33\times 10^{-4}$, (c) $E=2.50\times 10^{-4}$, (d) $E=2.00\times 10^{-4}$, (e) $E=1.67\times 10^{-4}$, (f) $E=1.25\times 10^{-4}$, (g) $E=1.00\times 10^{-4}$, and (h) $E=9.09\times 10^{-5}$.

4. The vorticity skewness at the convective onset stage

4.1. Vertical mode decomposition

This subsection studies what controls $S$ at the convective onset stage. The main idea is to perform a vertical mode decomposition, and study the contribution from each mode. We decompose velocity $(u,v,w)$, horizontal vorticity $\boldsymbol {\omega }_h$, vertical vorticity $\boldsymbol {\omega }_z$, pressure $p$, and temperature $T$ into different vertical modes:

(4.1)\begin{equation} \left.\begin{gathered} u = u_0 + u_1 + u_2 +\cdots, \\ v = v_0 + v_1 + v_2 +\cdots, \\ w = w_0 + w_1 + w_2 +\cdots, \\ \boldsymbol{\omega}_h=\boldsymbol{\omega}_{h,0} +\boldsymbol{\omega}_{h,1}+ \boldsymbol{\omega}_{h,2} +\cdots, \\ \omega_z = \omega_{z,0} + \omega_{z,1} + \omega_{z,2} +\cdots, \\ p = p_0 + p_1 + p_2 +\cdots, \\ T = T_0 + T_1 + T_2 +\cdots, \end{gathered}\right\} \end{equation}

where the subscript $n$ denotes the $n\textrm {th}$ vertical mode, with vertical structure $\sin (n {\rm \pi}z)$ or $\cos (n {\rm \pi}z)$. For the weakly nonlinear regime, only the $n=1$ mode is susceptible to linear instability, so the $n=0$ and $n=2$ modes are driven nonlinearly by the $n=1$ quantities. Thus the $n=0$ and $n=2$ modes are smaller than the $n=1$ modes by an order of ${Ro}$ (strictly speaking, ${Ro_g}$). The $n \geq 3$ modes are higher-order nonlinear effects neglected in this analysis.

We express the third-order moment of vorticity, which is the numerator of $S$, with the $n=0, 1, 2$ modes:

(4.2)\begin{align} \left \langle \overline{\omega^3_z} \right \rangle & \approx \left \langle \overline{ \left( \omega_{z,0} +\omega_{z,1} +\omega_{z,2} \right)^3 } \right \rangle \nonumber\\ &= \langle \overline{\omega^3_{z,0}} + \overline{\omega^3_{z,1}} + \overline{\omega^3_{z,2}} \rangle + 6 \langle \overline{\omega_{z,0} \omega_{z,1} \omega_{z,2}} \rangle \nonumber\\ &\quad + 3 \langle \overline{\omega^2_{z,0} \omega_{z,1}} + \overline{\omega_{z,0} \omega^2_{z,1}} + \overline{\omega^2_{z,0} \omega_{z,2}} + \overline{\omega_{z,0} \omega^2_{z,2}} + \overline{\omega^2_{z,1} \omega_{z,2}} + \overline{\omega_{z,1} \omega^2_{z,2}}\rangle \nonumber\\ &\approx 3 \langle \overline{\omega^2_{z,1} \omega_{z,0}} \rangle + 3 \langle \overline{\omega^2_{z,1} \omega_{z,2}} \rangle. \end{align}

The first term of the fourth line is essentially the $3 \langle \overline { \overline {\omega _z} {\omega '_z}^2 } \rangle$ term in (3.2), and the second term is essentially the $\langle \overline { {\omega '_z}^3 } \rangle$ term. Many terms have been neglected, as explained below. Some terms are exactly zero:

(4.3)\begin{equation} \left.\begin{gathered} \langle \overline{ \omega^3_{z,1} } \rangle \propto \overline{\cos^3 \left( {\rm \pi}z \right)}=0, \quad 6 \langle \overline{\omega_{z,0} \omega_{z,1} \omega_{z,2}} \rangle \propto \overline{\cos \left( {\rm \pi}z \right) \cos \left( 2 {\rm \pi}z \right) }=0,\\ 3 \langle \overline{\omega^2_{z,0} \omega_{z,1}} \rangle \propto \overline{ \cos \left( {\rm \pi}z \right) }=0, \quad 3 \langle \overline{\omega^2_{z,0} \omega_{z,2}} \rangle \propto \overline{ \cos \left( 2 {\rm \pi}z \right) }=0,\\ 3 \langle \overline{\omega_{z,1} \omega^2_{z,2}} \rangle \propto \overline{ \cos \left( {\rm \pi}z \right) \cos^2 \left( 2 {\rm \pi}z \right) }=0. \end{gathered}\right\} \end{equation}

The rest of the neglected terms, i.e. $\langle \overline {\omega ^3_{z,0}} \rangle$, $\langle \overline {\omega ^3_{z,2}} \rangle$ and $3 \langle \overline {\omega _{z,0} \omega ^2_{z,2}} \rangle$, are order ${Ro}^2$ smaller than the retained terms. This section is devoted to understanding the relationship between $\omega _{z,0}$, $\omega _{z,1}$ and $\omega _{z,2}$, which involves the full coupling between dynamical and thermodynamic processes.

4.2. A Galerkin model

We derive a Galerkin model of $n=0, 1, 2$ modes from a simplified set of equations. As a simplification from Veronis (1959), we apply two assumptions that work for the rapidly rotating and weakly nonlinear regime.

1. (i) The divergence equation uses geostrophic balance, which requires ${Ro} \ll 1$. This assumption is crucial for deriving the quasi-geostrophic omega equation in § 4.3, which helps us to understand the origin of $S$.

2. (ii) The $n=1$ mode is the only unstable vertical mode, which requires $8.7 \lesssim \widetilde{{Ra}} \lesssim 21.9$. This is because the critical $\widetilde {{Ra}}$ of the $n$th mode in the rapidly rotating limit approximately obeys $8.7 n^{4/3}$ (Chandrasekhar 1953). We take it as $\widetilde {{Ra}} \lesssim 20$ for simplicity.

A natural consequence of the two assumptions is ${E}={Ro}^3\, \widetilde {{Ra}}^{-3/2}\,{Pr}^{3/2} \ll 1$. Because the most unstable horizontal total wavenumber is $K_m = ({\rm \pi} ^2/2)^{1/6} {E}^{-1/3}$, we must have $K_m \gg 1$, and the vortices must have a small width-to-depth ratio (Chandrasekhar 1953). This allows us to neglect viscosity and thermal diffusion in the vertical direction. The simplified governing equation set is shown below, with linear terms on the left-hand side and nonlinear terms on the right-hand side:

(4.4)$$\begin{gather} \left( \frac{\partial}{\partial t} - {Ra}^{-1/2}\,\nabla_h^2 \right) \omega_z - \frac{1}{{Ro}}\,\frac{\partial w}{\partial z} =- \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla}\omega_z + \boldsymbol{\omega} \boldsymbol{\cdot}\boldsymbol{\nabla} w, \end{gather}$$

(4.5)$$\begin{gather}\nabla^2_h p - \frac{1}{{Ro}}\,\omega_z = 0, \end{gather}$$

(4.6)$$\begin{gather}\left( \frac{\partial}{\partial t} - {Ra}^{-1/2}\,\nabla_h^2 \right) w + \frac{\partial p}{\partial z} - T =- \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}w, \end{gather}$$

(4.7)$$\begin{gather}\left( \frac{\partial}{\partial t} - {Ra}^{-1/2}\,\nabla_h^2 \right) T =- \boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}T, \end{gather}$$

(4.8)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u} = 0, \end{gather}$$

where $\boldsymbol {\nabla }_h \equiv \boldsymbol {i}\,\partial /\partial x +\boldsymbol {j}\,\partial /\partial y$ denotes the horizontal gradient. Equation (4.5) is the divergence equation using the geostrophic balance approximation. Also, note the $\nabla ^2 \approx \nabla ^2_h$ approximation used in the viscosity and diffusion terms.

The simplified equation set is similar to yet different from the non-hydrostatic quasi-geostrophic (NHQG) equation proposed by Sprague et al. (2006), which is derived in the regime of ${Ro} \ll 1$ and small width-to-depth ratio, using multiple time scale asymptotic expansion. The rigorous NHQG equation yields zero $S$. We retain a few terms neglected in the NHQG scaling, which are crucial for generating vorticity asymmetry.

1. (i) In the vertical vorticity equation (4.4), the advection of $\omega _z$ by the toroidal velocity, the stretching of $\omega _z$, and tilting terms are retained. The toroidal velocity includes the vertical velocity and the curl-free component of horizontal velocity.

2. (ii) In the $w$ equation (4.6), the advection of $w$ by the toroidal velocity is retained.

3. (iii) In the $T$ equation (4.7), the advection of perturbation temperature $T$ by the toroidal velocity is retained. This term is partially represented in the NHQG equation by allowing the horizontally averaged temperature perturbation $\langle T \rangle$ evolve.

We will show that factor (i) contributes positively to $S$, and factors (ii) and (iii) contribute negatively to $S$, at the convective onset stage. Substituting (4.1) into the simplified equation set (4.4)–(4.8), we obtain a Galerkin model of the $n=1$, $n=2$ and $n=0$ modes relevant for calculating skewness.

The $n=1$ mode equation is

(4.9)$$\begin{gather} \left( \frac{\partial}{\partial t} - {Ra}^{-1/2}\,\nabla_h^2 \right) \omega_{z,1} - \frac{1}{{Ro}}\,\frac{\partial w_1}{\partial z} = 0, \end{gather}$$

(4.10)$$\begin{gather}\nabla^2_h p_1 - \frac{1}{{Ro}}\,\omega_{z,1}= 0, \end{gather}$$

(4.11)$$\begin{gather}\left( \frac{\partial}{\partial t} - {Ra}^{-1/2}\,\nabla_h^2 \right) w_1 + \frac{\partial p_1}{\partial z} - T_1 = 0, \end{gather}$$

(4.12)$$\begin{gather}\left( \frac{\partial}{\partial t} - {Ra}^{-1/2}\,\nabla_h^2 \right) T_1 - w_1 = 0, \end{gather}$$

(4.13)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_1 = 0. \end{gather}$$

The eigenvalue of the $n=1$ equation is the growth rate $\sigma$ of the linear instability:

(4.14)\begin{equation} \sigma = \left( \frac{K^2 - {Ro}^{-2}\,{\rm \pi}^2}{K^2} \right)^{1/2} - {Ra}^{-1/2}\,K^2. \end{equation}

The $n=2$ mode governing equation is

(4.15)$$\begin{gather} \left( \frac{\partial}{\partial t} - {Ra}^{-1/2}\,\nabla_h^2 \right) \omega_{z,2} - \frac{1}{{Ro}}\,\frac{\partial w_2}{\partial z} = 0, \end{gather}$$

(4.16)$$\begin{gather}\nabla^2_h p_2 - \frac{1}{{Ro}}\,\omega_{z,2} = 0, \end{gather}$$

(4.17)$$\begin{gather}\left( \frac{\partial}{\partial t} - {Ra}^{-1/2}\,\nabla_h^2 \right) w_2 + \frac{\partial p_2}{\partial z} - T_2 =- \boldsymbol{u}_1 \boldsymbol{\cdot} \boldsymbol{\nabla} w_1, \end{gather}$$

(4.18)$$\begin{gather}\left( \frac{\partial}{\partial t} - {Ra}^{-1/2}\,\nabla_h^2 \right) T_2 - w_2 =- \boldsymbol{u}_1 \boldsymbol{\cdot}\boldsymbol{\nabla} T_1, \end{gather}$$

(4.19)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}_2 = 0. \end{gather}$$

For the $n=0$ mode, only the vertical vorticity equation is needed:

(4.20)\begin{align} \left( \frac{\partial}{\partial t} - {Ra}^{-1/2}\,\nabla_h^2 \right) \omega_{z,0} &=- \underbrace{ w_1\,\frac{\partial \omega_{z,1}}{\partial z} }_{\textit{vertical advection}} + \underbrace{\omega_{z,1}\,\frac{\partial w_1}{\partial z}}_{\textit{stretching}} \nonumber\\ &\quad - \underbrace{ \left( \boldsymbol{u}_1 \boldsymbol{\cdot} \boldsymbol{\nabla}_h \right) \omega_{z,1} }_{\textit{horizontal advection}} + \underbrace{\left( \boldsymbol{\omega}_{h,1} \boldsymbol{\cdot} \boldsymbol{\nabla}_h \right) w_1 }_{\textit{tilting}}. \end{align}

The stretching of solid-body vorticity cannot produce barotropic vorticity because the barotropic mode is purely two-dimensional. We call the right-hand side terms of (4.20) ageostrophic effects because they do not exist in the NHQG scaling (Sprague et al. 2006).

The Galerkin model has been greatly simplified by the assumption that the instability is dominated by a single mode with a unique horizontal and vertical wavenumber. The $- \boldsymbol {u}_1 \boldsymbol {\cdot } \boldsymbol {\nabla } w_1$ and $- \boldsymbol {u}_1 \boldsymbol {\cdot } \boldsymbol {\nabla } T_1$ terms only project onto the $n=2$ mode:

(4.21)\begin{equation} - \boldsymbol{u}_1 \boldsymbol{\cdot} \boldsymbol{\nabla} w_1,\ - \boldsymbol{u}_1 \boldsymbol{\cdot} \boldsymbol{\nabla} T_1 \propto \sin ( {\rm \pi}z ) \cos ( {\rm \pi}z ) \sim \sin ( 2 {\rm \pi}z ). \end{equation}

The sum of the vertical advection and stretching of $\omega _{z,1}$ only projects onto the $n=0$ mode vorticity equation:

(4.22)\begin{equation} - w_1\,\frac{\partial \omega_{z,1}}{\partial z} + \omega_{z,1}\,\frac{\partial w_1}{\partial z} \propto \sin^2 ( {\rm \pi}z ) + \cos^2 ( {\rm \pi}z ) \sim 1. \end{equation}

In Appendix B, we prove that with an additional axisymmetric assumption, which approximately works for vortices at the convective onset stage, the sum of the horizontal advection and tilting terms also only project onto the $n=0$ mode vorticity equation.

Next, we study the two parts of $S$, which is from the interaction between the $n=1$ and $n=2$ modes, and from the interaction between the $n=1$ and $n=0$ modes:

(4.23)\begin{equation} S\approx 3\,\frac{\langle \overline{\omega^2_{z,1} \omega_{z,2} } \rangle}{ \langle \overline{\omega^2_{z,1}} \rangle^{3/2} } + 3\,\frac{\langle \overline{ {\omega^2_{z,1}} \omega_{z,0} } \rangle}{ \langle \overline{\omega^2_{z,1}} \rangle^{3/2} }. \end{equation}

Note that the self-interaction of the $n=1$ mode does not generate vorticity skewness, because $\langle \overline { \omega ^3_{z,1} } \rangle =0$, as shown in (4.3).

4.3. The negative contribution to $S$ from the $n=2$ mode

The $n=2$ vertical vorticity equation (4.15) has an intriguing property: it is linear. Thus $\omega _{z,2}$ is only driven by the stretching of solid-body vorticity by $w_2$. The quasi-geostrophic approximation on the divergence equation (4.16) allows us to derive the diagnostic equation of $w_2$, essentially the quasi-geostrophic omega equation (e.g. Holton 2004) in the non-hydrostatic state with unstable stratification. Combining (4.15)–(4.19), we obtain

(4.24)\begin{align} &\left[ \underbrace{ \frac{1}{{Ro}^2}\,\frac{\partial^2}{\partial z^2} }_{\textit{from}\ p_2,\,<0} \underbrace{ - \nabla^2_h }_{\textit{from}\ T_2\ \textit{adv},\,>0} + \underbrace{ \left( \frac{\partial}{\partial t} - {Ra}^{-1/2}\,\nabla_h^2 \right)^2 \nabla_h^2}_{\textit{non-hydrostatic},\,<0} \right] w_2 \nonumber\\ &\quad =- \left( \frac{\partial}{\partial t} - {Ra}^{-1/2}\,\nabla_h^2 \right) \nabla_h^2 \left( \boldsymbol{u}_1 \boldsymbol{\cdot}\boldsymbol{\nabla} w_1 \right) - \nabla_h^2 \left( \boldsymbol{u}_1 \boldsymbol{\cdot}\boldsymbol{\nabla} T_1 \right). \end{align}

The operator on the left-hand side includes three important terms.

1. (i) The pressure gradient part $({1}/{{Ro}^2})({\partial ^2}/{\partial z^2})$, which suppresses $w_2$. This is because convergence produces a cyclone and therefore a low-pressure anomaly, and vice versa for divergence. Thus pressure gradient force always points from a divergent zone to a convergent zone, suppressing $w_2$.

2. (ii) The temperature advection part $\nabla _h^2$, which amplifies $w_2$. This term originates from the $w_2$ term in (4.18). It amplifies $w_2$ because a higher $w_2$ increases the buoyancy ($T_2$) and accelerates an updraft, and similarly for a downdraft. This term has an opposite sign in the traditional omega equation applied to a stably stratified fluid.

3. (iii) The non-hydrostatic term, which suppresses $w_2$.

The omega equation helps us to infer the trend of $w_2$ from the forcing terms. First, we analyse the left-hand-side operator of (4.24). The first and third terms yield a multiplier with a minus sign, and the second term yields a positive sign. We can determine the sign of the left-hand-side operator without solving the omega equation. If the $n=1$ variables were substituted into the left-hand side operator,

(4.25a–c)\begin{equation} n=1,\quad -\nabla^2_h=K_m^2, \quad \frac{\partial}{\partial t} = \sigma, \quad \frac{\partial^2}{\partial z^2}=-{\rm \pi}^2, \end{equation}

then the left-hand side would be zero, providing a reference. Now we substitute in the $n=2$ variables (e.g. $w_2$). Because the right-hand-side forcing terms of (4.24) are product terms of the $n=1$ mode, we must have $\partial /\partial t=2\sigma$, and the horizontal structure of the $n=2$ mode must be approximately twice as fine-grained as the $n=1$ mode. Thus we set $-\nabla ^2_h \approx \alpha ^2 K^2_m$, where $\alpha$ is an unknown parameter, diagnosed to be $1<\alpha <2$ in Appendix C (figure 12). The vertical structure of the $n=2$ mode is also more fine-grained, yielding $\partial ^2/\partial z^2=-4{\rm \pi} ^2$. This is summarized as

(4.26a–c)\begin{equation} n=2,\quad -\nabla^2_h \approx \alpha^2 K_m^2,\quad \frac{\partial}{\partial t} = 2 \sigma, \quad \frac{\partial^2}{\partial z^2}=-4{\rm \pi}^2. \end{equation}

Using (4.25a–c) and (4.26a–c), we calculate the change of magnitude of each left-hand-side term when switching from $n=1$ to $n=2$:

(4.27)\begin{equation} \left.\begin{gathered} \frac{1}{{Ro}^2}\,\frac{\partial^2}{\partial z^2}\times 4, \\ - \nabla^2_h \times \alpha^2,\\ \left( \frac{\partial}{\partial t} - {Ra}^{-1/2}\,\nabla_h^2 \right)^2 \nabla_h^2 \times \underbrace{\frac{(2 \sigma + \alpha^2 {Ra}^{-1/2}\,K^2_m)^2 \alpha^2}{(\sigma+{Ra}^{-1/2}\,K^2_m)^2}}_{> \alpha^2}. \end{gathered}\right\} \end{equation}

Because $1<\alpha <2$, the first and third terms of (4.27) are amplified more than the second term when switching from $n=1$ to $n=2$. Thus they dominate the sign of the left-hand-side operator, and the operator must correspond to a negative multiplier.

The above analysis shows that along the vortex axis, $w_2$ must have the same sign as $-\boldsymbol {u}_1 \boldsymbol {\cdot } \boldsymbol {\nabla } w_1$ and $-\boldsymbol {u}_1 \boldsymbol {\cdot } \boldsymbol {\nabla } T_1$, which causes $w$ to deform downwind vertically. The downwind deformation denotes the shift of the zero divergence height towards the outflow zone, essentially an intensification of outflow (figure 7a). As a result, the convergent zone is more diluted than the divergent zone, making the magnitude of the cyclone produced by stretching the solid-body vorticity smaller than the magnitude of the anticyclone produced by squashing the solid-body vorticity. This explains why the $n=2$ mode contributes negatively to $S$.

Figure 7. A schematic illustration of the vorticity skewness generation along the vortex axis at the convective onset stage. An updraft is used as an example. (a) The vertical velocity downwind deformation (outflow intensification) enhances the squashing of solid-body vorticity, generating a strong anticyclone (AC). It contributes negatively to $S$ via the interaction between the $n=2$ and $n=1$ modes. (b) The vertical advection and stretching of $\omega _z$ produce a barotropic cyclone (CC) along the vortex axis. It contributes positively to $S$ via the interaction between the $n=0$ and $n=1$ modes. The solid blue lines denote (a) the $n=1$ mode vertical velocity and (b) the vertical vorticity profile along the vortex axes. The solid red lines denote (a) the bulk vertical velocity and (b) the vertical vorticity profile. The dashed black lines are zero-value reference lines.

Here are physical explanations for the downwind deformation. We take the updraft as an example. Here, $-\boldsymbol {u}_1 \boldsymbol {\cdot } \boldsymbol {\nabla } w_1$ denotes the inertia of a parcel. Buoyancy remains positive along the air column, so the parcel keeps accelerating until steered to divergence near the upper plate, causing a downwind deformation of the vertical velocity profile. Meanwhile, the upward advection of temperature makes the temperature lapse rate smaller at the lower level and larger at the upper level. This reduces the low-level convective instability and enhances the upper-level convective instability. Thus the peak height of vertical velocity must shift to a higher level, also causing a downwind deformation. An extreme example of downwind deformation of vertical velocity is a buoyant plume, which keeps entraining (converging) fluids unless a ceiling forces it to diverge (e.g. Rooney & Linden 2012).

4.4. The positive contribution to $S$ from the $n=0$ mode

Equation (4.20) shows that along the vortex axis, the barotropic vorticity is produced by the vertical advection term $-w_1({\partial \omega _{z,1}}/{\partial z})$ and the stretching term $\omega _{z,1}({\partial w_1}/{\partial z})$. For an updraft, there is a low-level cyclone and an upper-level anticyclone. The down-gradient advection produces an overall cyclonic anomaly. The stretching on the low-level cyclone and the squashing on the upper-level anticyclone produce an overall cyclonic anomaly. Thus both terms generate cyclonic anomaly along the vortex axis, producing a barotropic cyclone, as illustrated in figure 7(b). The cyclonic bias produced by these mechanisms has been identified in many previous studies of rotating fluids (e.g. Schubert & Alworth 1987; Morize et al. 2005; Majda, Mohammadian & Xing 2008). The only discussion in the context of RRBC seems to be from Guervilly et al. (2014). They explained the cyclonic bias as the concentration of cyclonic vorticity by the convergent flow and the dilution of anticyclonic vorticity by the divergent flow. It is an explanation based on the flux-form vorticity equation (e.g. Haynes & McIntyre 1987), consistent with the explanation based on the stretching of $\omega _z$.

Readers might ask: what is the role of the horizontal advection and tilting terms? They are zero along the vortex axis. Tory, Montgomery & Davidson (2006) made a budget analysis of a simulated tropical cyclone precursor vortex. They found a core–shield structure of barotropic vorticity, with the horizontal advection and tilting terms producing an anticyclonic shield around a cyclonic core. The shield structure in the RRBC problem is theoretically confirmed in Appendix B and visible in the simulation result (figure 11(c) of Appendix C). Note that this shield is a nonlinearly generated depth-averaged structure that needs to be compared cautiously with the RRBC literature, where the shield (or sleeve) is more loosely defined as a ring of opposite-sign vorticity to the vortex core that is also seen in the linear instability (e.g. Sprague et al. 2006; Portegies et al. 2008; Grooms et al. 2010; Shi et al. 2020). Because the anticyclonic shield renders $\omega _{z,0}<0$, it produces negative vorticity skewness via $\langle \overline {\omega _{z,0} \omega ^2_{z,1}} \rangle$. However, because the shield is off the vortex axis, its negative contribution is minor compared to the positive contribution by vertical advection and stretching along the vortex axis. In fact, the bulk effect of the horizontal and vertical advection of $\omega _z$ does not produce volumetric vorticity skewness, even though it could generate local vorticity skewness along the vortex axis, as will be shown in the budget analysis of the $\langle \overline {\omega _z^3} \rangle$ equation in § 5. This is because the positive skewness produced by the vertical advection along the vortex axis is offset by the negative skewness produced by the horizontal advection in the anticyclonic shield. Given that advection does not produce volumetric vorticity skewness, we do not highlight the role of vertical advection as a generator of vorticity skewness.

In summary, $S$ at the convective onset stage is controlled by two factors. The downwind deformation of $w$ is inherently non-hydrostatic, contributing negatively to $S$. The production of barotropic vorticity by vortex tube stretching and tilting is inherently ageostrophic, contributing positively to $S$. Despite the rich insights, the above qualitative analysis cannot explain the $-1:2$ relative contribution of the $n=2$ and $n=0$ modes in DNS. It also cannot explain why $S \propto {Ro_g}$ and determine the proportional factor, which does not depend on ${Ra}$ and ${E}$, and is approximately 2.5 in DNS. In Appendix C, we perform an order-of-magnitude estimation of $S$ via a single-scale approximation on the Galerkin model, and by using the vorticity skewness along the vortex axis to estimate the volumetric skewness. The contributions of the $n=2$ and $n=0$ modes to $S$ are confirmed theoretically to be of the same order of magnitude. We show theoretically that $S/{Ro_g}\sim 1.5$, which is close to the $S/{Ro_g}\approx 2.5$ result in DNS.

5. The vorticity skewness at the equilibrium stage

The convective onset stage has a relatively simple flow pattern with erect vortices, and is a convenient starting point. This section studies the equilibrium stage and investigates how much of the mechanism discussed at the convective onset stage still holds. We define the equilibrium stage as the period where flow statistics do not vary significantly with time. In practice, we use the time-averaged data between $t=90$ and $t=120$ to diagnose its behaviour. The slot is marked as the red shading in figures 3 and 4. To test the convergence of the averaging, figure 8, which is associated with the equilibrium stage, is recalculated between $t=90$ and $t=105$, and between $t=105$ and $t=120$, and deposited in the supplementary material. No qualitative differences are found, so the averaging between $t=90$ and $t=120$ should be an adequate representation of the equilibrium-state statistics.

Figure 8. Statistics of the equilibrium stage (between $t=90$ and $t=120$) for experiments (a,c,e) Ra1–Ra8 and (b,d,f) E1–E8. (a,b) The budget of the $\langle \overline {\omega ^3_z} \rangle$ equation. The blue, red, yellow and purple lines show the $\omega _z$ stretching term, the tilting term, the solid-body vorticity stretching term, and the viscous term of (5.2). Each term has been normalized by the absolute value of the viscous term. (c,d) The volumetric skewness of $\partial w/\partial z$. (e,f) The vorticity skewness $S$ versus ${Ro_g}$. The blue lines show the diagnosed $S$. The solid black lines show the convective onset stage fitting: $S=2.5 {Ro_g}$. The dashed black lines show the heuristic equilibrium stage theory: $S=2.5\,{Ro_g}\exp (-0.25\,{Ro_g}/{Ro})$.

5.1. Budget analysis of the $\langle \overline {\omega ^3_z} \rangle$ equation

For the equilibrium stage, we have tried the analysis with the barotropic–baroclinic decomposition, but found a conversion between the barotropic–baroclinic term and the purely baroclinic term due to the unsteady flow, which complicates the attribution of physical mechanisms. Thus we take another route by making a budget analysis of the $\langle \overline {\omega ^3_z} \rangle$ equation. Multiplying by $3 \omega ^2_z$ on both sides of the vertical vorticity equation

(5.1)\begin{equation} \frac{\partial \omega_z}{\partial t} =- \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \omega_z + \boldsymbol{\omega} \boldsymbol{\cdot} \boldsymbol{\nabla} w + \frac{1}{{Ro}}\,\frac{\partial w}{\partial z} + \left( \frac{{Pr}}{{Ra}} \right)^{1/2} \nabla^2 \omega_z \end{equation}

and making a volumetric average, we get

(5.2) \begin{align} \frac{{\rm d} \langle \overline{\omega^3_z} \rangle}{{\rm d} t} &= \overbrace{\underbrace{ 3 \left\langle \overline{ \omega^3_z\,\frac{\partial w}{\partial z} } \right\rangle }_{\omega_z\ \textit{stretching}} + \underbrace{ 3 \langle \overline{\omega^2_z \left( \boldsymbol{\omega}_h \boldsymbol{\cdot}\boldsymbol{\nabla}_h w \right) } \rangle }_{\textit{tilting}} }^{\textit{ageostrophic}}\nonumber\\ &\quad + \overbrace{ \underbrace{ 3 \left\langle \overline{ \frac{1}{{Ro}}\,\omega^2_z\, \frac{\partial w}{\partial z} } \right\rangle }_{\textit{solid-body stretching}}}^{\textit{non-hydrostatic}} \underbrace{{}- 6 \left( \frac{{Pr}}{{Ra}} \right)^{1/2} \langle \overline{ \omega_z\, |\boldsymbol{\nabla}\omega_z|^2 } \rangle }_{\textit{viscosity}}. \end{align}

The right-hand side has four terms. The $\omega _z$ stretching and tilting terms carry the ageostrophic effects. The solid-body vorticity stretching term carries the non-hydrostatic effect, which causes the asymmetry between inflow and outflow. Figures 8(a,b) show that the sum of the four terms is slightly positive for ${Ro_g}\gtrsim 1$, a quantity that should equal zero at the equilibrium stage. We have carefully examined the budget calculation and infer it to be due to the numerical diffusion of the fifth-order advection scheme that could dissipate $\langle \overline{\omega^3_z} \rangle$.

The primary balance of the right-hand side is between the ageostrophic terms ($\omega _z$ stretching and tilting) that generate skewness and the viscosity that dissipates skewness. The negative contribution of tilting is likely due to the anticyclonic shield around a cyclonic core (Appendix B). The solid-body vorticity stretching term has an intriguing transition behaviour. It reduces skewness in the ${Ro_g}\lesssim 0.5$ range, and slightly increases it for a higher ${Ro_g}$. Because this term carries the asymmetry between inflow and outflow, the transition indicates the suppression of outflow strength as the flow gets more unsteady (figure 9a). Figures 8(c,d) confirm the drastic change in the inflow–outflow asymmetry by showing that the skewness of $\partial w/\partial z$ first drops and then increases as ${Ro_g}$ increases. This can be explained as a well-known process in plume dynamics: a stronger crossflow enhances entrainment dilution and makes an inflow turn to outflow at a lower depth (e.g. Lavelle 1997; Devenish et al. 2010). In the context of RRBC, this agrees qualitatively with the existing understanding that the turbulent dilution of the outflow generates positive vorticity skewness in the turbulent regime (e.g. Vorobieff & Ecke 2002), though we should be cautious that the equilibrium-state skewness of $\partial w/\partial z$ remains negative in our DNS. Our analysis reveals a competition in generating vorticity skewness between the stiffening role of the non-hydrostatic effect, which intensifies the outflow and drives the inflow–outflow asymmetry, and the disturbing role of the unsteady flow, which dilutes the outflow and weakens the inflow–outflow asymmetry.

Figure 9. A schematic diagram for the break of vertical coherence by the vertical shear induced by other convective vortices. The shear could (a) suppress the intensity of outflow and (b) cause misalignment between the plume and the vortex axis.

5.2. A heuristic skewness generation efficiency

The above analysis shows that the break of vortex coherency by the unsteady shear flow is the main difference between the onset and equilibrium stages. The unsteady flow suppresses outflow intensification (figure 9a). We infer that the shear could also tilt the plume, causing a misalignment between the plume and the vortex axis that makes the stretching of $\omega _z$ less efficient in producing a barotropic cyclone (figure 9b). Thus both the positive and negative skewness generation factors at the convective onset stage might be suppressed. The DNS result (figures 8e,f) also indicates that the equilibrium-state $S$ drops below the onset stage scaling $S \approx 2.5\,{Ro_g}$. Thus skewness might be generated at the highest efficiency at the convective onset stage, where the flow is stationary.

We introduce a heuristic ‘skewness generation efficiency’ to parametrize the deviation from $S \approx 2.5\,{Ro_g}$ at the equilibrium stage due to the unsteady flow. For simplicity, we assume that the vorticity stretching and outflow intensification are suppressed equally. We let the efficiency be a decreasing function of the vertical overturning time scale over the horizontal overturning time scale. A faster vertical overturning allows less time for eddies in the horizontal plane to tilt the plume. The vertical overturning time scale is estimated as the free-fall time scale (unity in this non-dimensional framework), and the horizontal overturning time scale is estimated as the inverse of vorticity scale $\langle \overline {\omega ^2_z} \rangle ^{-1/2}={Ro}/{Ro_g}$. We propose a heuristic expression for $S$ that works for the ${Ro_g}\lesssim 1$ regime at the equilibrium stage:

(5.3)\begin{equation} S = 2.5\,{Ro_g} \exp\left(-0.25\,\frac{{Ro_g}}{{Ro}} \right). \end{equation}

Here, the exponential function is a heuristic choice, and the $0.25$ factor is obtained by fitting. Equation (5.3) agrees well with DNS in the ${Ro_g}\lesssim 1$ regime (figures 8e,f). The next step is to solidify the theoretical basis of (5.3) by studying the resilience of a convective vortex in random vertical and horizontal shear, a topic also of interest to tropical cyclone researchers (e.g. Reasor, Montgomery & Grasso 2004; Tao & Zhang 2015).

6. Conclusion and discussion

Rotating Rayleigh–Bénard convection (RRBC) is a prototype model of geophysical and astrophysical convection. Like many other rotating fluid systems, cyclones tend to have higher vorticity magnitudes than anticyclones, rendering a cyclonic bias. The bias has been explained with three mechanisms. First, $\omega _z$ is stretched more efficiently in the cyclonic region than squashed in the anticyclonic region due to the influence of $\omega _z$ on the absolute vorticity (e.g. Morize et al. 2005; Guervilly et al. 2014). Second, turbulent mixing dilutes the outflow of a plume, making the vertical motion preferentially stretch the solid-body vorticity rather than squash it (e.g. Julien et al. 1996a; Vorobieff & Ecke 2002). Third, an anticyclone with negative absolute vorticity is unstable to centrifugal instability (e.g. Kunnen et al. 2010b). Despite the progress in identifying physical mechanisms, the current understanding of the cyclonic bias in rotating convection is not systematic. It remains unclear which mechanism dominates at which stage or regime.

A useful metric of cyclonic bias is the vorticity skewness. An important question is how the skewness depends on the Rossby number, which has not been solved in any regime or stage of RRBC. The convective onset stage of rapidly rotating RRBC, which is weakly nonlinear (finite amplitude) and has an organized flow structure, provides a straightforward starting point. To obtain a concise result, properly defining the skewness and Rossby number is crucial. Using DNS with a fixed ${Pr}=1$ and varying ${Ra}$ and ${E}$, we find that the volumetric vorticity skewness $S$ is proportional to the global Rossby number ${Ro_g}$ by a constant factor approximately 2.5. It works for the convective onset stage of the ${Ro} \ll 1$ and $\widetilde {{Ra}} \lesssim 20$ regime.

To understand this phenomenon, we perform a vertical eigenmode decomposition and find that $S$ is produced by the interaction between the $n=0$, $n=1$, and $n=2$ modes. We derive a Galerkin model following the asymptotic analysis procedure of Veronis (1959), but with several simplifications. A crucial simplification is applying the quasi-geostrophic assumption to the divergence equation, consistent with the NHQG model of Sprague et al. (2006). Our model also retains some terms neglected in the NHQG model that generate skewness. They include the advection of $\omega _z$, $T$ and $w$ by toroidal velocity, as well as the vertical stretching and tilting terms in the $\omega _z$ equation. An approximate solution of the Galerkin model confirms $S\propto {Ro_g}$ and shows a relation $S/{Ro_g}\sim 1.5$ that does not depend on ${Ra}$ and ${E}$, not far from the 2.5 value in the DNS.

Figure 10 summarizes the mechanism of vorticity skewness at the convective onset stage. The linearly unstable $n=1$ mode does not generate vorticity skewness itself. It nonlinearly drives the $n=0$ and $n=2$ modes, and interacts with them to generate skewness. The key mechanisms include:

1. (i) the interaction between the $n=1$ and $n=2$ modes, a non-hydrostatic effect;

2. (ii) the interaction between the $n=1$ and $n=0$ modes, an ageostrophic effect.

By deriving the quasi-geostrophic omega equation of the $n=2$ mode, we find the downwind deformation of vertical velocity (in other words, outflow intensification) due to the nonlinear advection of $w$ and $T$ makes a negative contribution to $S$. This results from non-hydrostatic effects, a factor that has not received much attention. The mechanisms associated with the $n=0$ mode are ageostrophic. The down-gradient vertical advection of $\omega _z$ and the more efficient stretching of absolute vorticity in the cyclonic region produces a barotropic cyclone along the vortex axis, contributing positively to $S$. Meanwhile, the horizontal advection of vertical vorticity and tilting of horizontal vorticity produce a barotropic anticyclonic shield, contributing negatively to $S$. Because the shield is off the vortex axis where the vorticity magnitude is the largest, its influence on $S$ is minor. In addition, we note that the bulk effect of the horizontal and vertical advection of $\omega _z$ cannot produce volumetric vorticity skewness. The DNS show that the positive contribution from ageostrophic effects ($n=0$ mode) is twice as strong as the negative contribution from non-hydrostatic effects ($n=2$ mode), making $S$ positive. By solving the quasi-geostrophic omega equation with a single-scale approximation, we theoretically confirm that these two terms are proportional, though not accurately enough to prove the $2:-1$ contribution ratio in the DNS.

Figure 10. A schematic diagram showing how the interaction of the $n=0$, $n=1$ and $n=2$ modes of the flow produces the vorticity skewness $S$ at the onset stage of rapidly rotating RRBC. It is the basis for understanding the equilibrium stage, where the flow is more unsteady.

Finally, we extend the theory from the convective onset to the equilibrium stage, where the flow is more unsteady. The budget analysis of the $\langle \overline { \omega ^3_z } \rangle$ equation shows that the stretching of $\omega _z$ is the main producer of vorticity skewness, mainly balanced by viscous dissipation. The non-hydrostatic effect in generating negative vorticity skewness is significantly suppressed for ${Ro_g}\gtrsim 0.5$, owing to the suppression of outflow intensification by the unsteady flow. Given that the unsteady vertical shear suppresses the outflow intensification and could also disturb the stretching of $\omega _z$ by tilting the plume, we heuristically parametrize the influence of shear on skewness as an efficiency factor, which is a decreasing function of the vertical overturning time scale over the horizontal overturning time scale. A unit efficiency represents the fully coherent state corresponding to the convective onset stage. By fitting one parameter, the agreement of the heuristic skewness theory with DNS is good in the ${Ro_g}\lesssim 1$ regime. An important future task is solidifying the physical explanation at the equilibrium stage by studying the vortex evolution under random vertical and horizontal shear.

Note that the theory only derives the $S\unicode{x2013} {Ro_g}$ relation and is not a fully predictive model because ${Ro_g}$ uses the diagnosed vorticity. To predict how $S$ depends on ${Ra}$ and ${E}$ at the equilibrium state, we can couple the $S\unicode{x2013} {Ro_g}$ relation with a model of how ${Ro_g}$ depends on ${Ra}$ and ${E}$. Sakai (1997) derived a model of the Rossby number defined with the horizontal velocity scale and length scale for the cellular regime, which might be a candidate model for ${Ro_g}$.

Examining the trend of the $S-{Ro_g}$ relation in a more turbulent regime is the next step. It remains unclear how the skewness transfers across scales in the turbulent RRBC and whether it could help to explain a critical problem: the formation of the large-scale vortex (LSV). The LSV is a barotropic turbulent cyclone produced by upscale energy transfer (Julien et al. 2012; Favier et al. 2014; Guervilly et al. 2014; Rubio et al. 2014; Guzmán et al. 2020). The LSV is not seen in our simulations. It typically takes much longer than the convective overturning time to form, and its formation criterion is still not well understood. Even within the finite-amplitude rapidly rotating regime, the investigation of $S$ is far from complete. Examining how much of the result still holds for ${Pr} \neq 1$ and a no-slip boundary is important. A no-slip boundary condition leads to Ekman pumping, a factor inferred to enhance the cyclonic bias via intensifying the inflow (e.g. Kunnen et al. 2006). Despite the limitations of the current theory, the analysis method in this paper could be generalized. First, studying the convective onset stage is a convenient intermediate step, which provides a short time slot to observe skewness at a relatively high ${Ro_g}$ without unsteady flow. At the onset stage, decomposing the vorticity skewness into the contributions from purely barotropic, purely baroclinic and barotropic–baroclinic modes is a useful diagnosis method. Second, focusing on the vorticity-rich region, i.e. the vortex axis, might be a useful simplification for studying vorticity skewness, given that the rotating flow tends to self-organize into columnar structures under the constraint of the Taylor–Proudman theorem. Third, using the quasi-geostrophic omega equation to quantify the asymmetry of convergent and divergent flow is a useful method. For example, it can be applied to evaluate the secondary circulation driven by Ekman pumping when the boundary is no-slip.