On the asymmetry of cyclones and anticyclones in the cellular regime of rotating Rayleigh-Bénard convection

Rotating Rayleigh-Bénard convection (RRBC) denotes the free convection between two parallel plates with a fixed temperature difference, placed in a rotating reference frame. It is a prototype model of geophysical and astrophysical convection. Rotation breaks the symmetry on its rotating axis, making the cyclones and anticyclones unequal in size and magnitude. Such an asymmetry has long been observed in experiments and simulations, but has not been explained with any theoretical model. A theory of such vorticity asymmetry is proposed specifically for the cellular regime, where background rotation is important and convection is weak. The property that columnar updraft and downdraft plumes are densely packed is shown to make the vertical vorticity profile at the vortex center approximately linear with height via thermal wind relation. This simplification of morphology enables a linkage between the vorticity strength of a plume which is quantified by vorticity Rossby number RoV, and the vorticity magnitude difference between the cyclonic and anticyclonic ends of plumes which is quantified with a nondimensional asymmetry factor δ. The lowest order relationship between δ and RoV is found to be constrained by vertical vorticity equation alone. An approximate analytical solution is found using asymptotic expansion, which shows that the asymmetry is generated mainly by the vertical advection and stretching of vertical vorticity in fluid interior, and is modified by the Ekman layer dynamics.

momentum [3]. Methods to quantify the asymmetry include investigating the probability distribution function (PDF) of vertical vorticity ω z and its skewness [4], the statistics of the angle between vorticity vector and the ambient rotation vector [5], and a direct vorticity magnitude comparison of the CC/AC ensemble if the vortex has a well-defined shape [6,7].
In the single layer shallow water equation which is the simplest paradigm of the large aspect ratio (width over height) stably-stratified rotating flow, ω z skews towards negative (stronger AC), because an AC area has thicker water depth, which leads to larger Rossby deformation radius that increases the vortex interaction range. This enhances mergers between ACs and therefore makes them less susceptible to the damping by viscosity [4,8]. When the relative vorticity is smaller compared to the ambient vorticity (small Rossby number), and the scale of motion becomes comparable to or smaller than Rossby deformation radius (weak layer thickness gradient), the motion is quasi-geostrophic and no CC/AC asymmetry exists [9]. The dominance of AC is also observed in free-evolving continuously stratified rotating flow due to the more efficient merger and vertical alignment that prevent substantial vorticity dissipation [10]. For 3D homogeneous rotating non-stratified turbulence, however, a CC is generally more long-lived and therefore stronger than an AC, mainly because an AC tends to be centrifugally unstable when the magnitude of ω z is larger than the background vorticity, as well as the effect of elliptical instability [5]. Centrifugal instability also exists in shallow water setup, but it is less effective in destabilizing anticyclones, probably because the vertical motion is weak in such a large aspect ratio setup [11].
In this paper, we focus on the asymmetry in 3D unstably-stratified rotating free convection, which has implications for the dynamics of open ocean convection, rotating cumulus cloud in tropical cyclone, Earth dynamo and star interior [12][13][14][15]. The strong vertical motion leads to vigorous vertical advection and stretching of vertical vorticity which is quite different from the stably-stratified flow. We study rotating Rayleigh-Bénard convection (RRBC) which is perhaps the simplest paradigm of rotating convection. It is Boussinesq convection between a pair of parallel warm lower plate and cold upper plate in a rotating frame [16]. RRBC is governed by three independent nondimensional parameters [17]: where β is the thermal expansion coefficient, g is the gravitational acceleration, ∆T is the temperature difference between the bottom and top, H is the depth of the fluid, f is Cori-olis parameter, ν and κ are the kinematic viscosity and diffusivity respectively. Rayleigh number Ra depicts the relative strength of the destabilizing buoyancy to the viscous and diffusive damping. Ekman number E represents the relative strength between viscous diffusion and rotational constraint, and (2E) 1/2 denotes the ratio of Ekman layer depth to the total depth. Prandtl number Pr is the ratio of kinematic viscosity to thermal diffusivity.
Combining the above dimensionless numbers, a convective Rossby number can be defined as Ro = E(Ra/Pr) 1/2 , which is regarded as a measure of the relative strength of thermal forcing to rotational effect [18].
Vorobieff and Ecke [1] showed that in the case of rapidly rotating convection, where Ro 1, the vertically-aligned columnar vortices dominate the flow. This phenomenon has been observed in many laboratory experiments [7,[19][20][21][22][23] and numerical simulations [18,[24][25][26][27][28]. The regime of RRBC can be further classified with a reduced Rayleigh number Ra= RaE 4/3 . It represents the extent to which the convective system is supercritical to neutral stability [16,17]. When Ra is small (less than 25 as is indicated by Stellmach et al. [17]), the flow is dominated by quasi-steady densely packed columnar vortices which barely interact with each other, called cellular or geostrophic regime. As Ra increases (but less than 70) it comes to convective Taylor column regime. The vortices are packed less densely, more significantly shielded, and vortex interaction remains weak. For higher Ra (over 70), organized vortices break into unsteady plumes and could organize into barotropic large-scale vortices due to upscale energy cascade [17].
RRBC is elegant because both its geometry and boundary condition are quite symmetric.
When the two plates are both no-slip or both free-slip, the only asymmetric factor is rotation.
Some laboratory experiments reported the skew towards anticyclone due to the centrifugal acceleration of the rotating apparatus [6,7], a factor that will not be considered in this work. Aside from this, the skew towards positive vorticity has been vaguely explained as being forced near the boundary where cyclones are ejected near the departure plate due to angular momentum convergence (or vorticity stretching) associated with the Ekman pumping and buoyant vertical motion, and are diluted in the fluid interior by turbulent mixing before reaching the destination plate [1,6,18,30,32,33]. Another factor is inertial instability which restricts the stable anticyclonic vertical vorticity magnitude to be smaller than f , and is especially clear-cut for the cellular regime where vortices are close to circular [32]. We are not aware of any theoretical model that quantitatively predicts the asymmetry in any regime. The relative contribution of Ekman pumping and the fluid interior factors to the asymmetry remains unknown.
We focus on the cellular regime which is weakly nonlinear, nonturbulent and has small Ro. We avoid inertial instability by considering only the weakly-convective (smaller Ro and Ra) cases where ω z > −f is always met. Portegies et al. [26] showed that a breaking of ω z symmetry does exist: within an updraft or downdraft in the cellular regime, the CC part has larger vorticity magnitude than the AC part, but this feature has not been explained. There are two candidate analytical vortex models for the cellular regime. One is by Portegies et al. [26] who solved the neutrally-stable mode of a linearized RRBC stability problem in cylindrical coordinate, and the other one is by Grooms et al. [28] who obtained a steady solution of non-hydrostatic quasi-geostrophic equation (NHQG) [25] and considered the nonlinear vertical advection by columnar vortices in shaping the basic state temperature.
Both used no-slip vertical boundary condition. They indicated that the main balance of vertical vorticity equation in fluid interior is the stretching of background vorticity and the horizontal diffusion of vertical relative vorticity. The lowest order balance of horizontal vorticity is close to but not the exact thermal wind balance due to the nonhydrostatic effect in such a small-aspect ratio vortex [21,25]. Neither of the models have considered CC/AC asymmetry. For Portegies et al. [26], the asymmetry naturally vanishes in the smallamplitude problem. For Grooms et al. [28], NHQG as a finite-amplitude model retains the horizontal advection of ω z but neglects the vertical advection, stretching and tilting of ω z .
This treatment, which is derived from rigorous multiscale expansion, can capture the basic flow pattern [17]. However, the opportunity for understanding the ω z asymmetry is missed.
To find the asymmetry, a natural path is using the first order solution (linear instability) to derive the complete next order one. However, the first order solution itself is only semianalytical, so this path is hard to provide any insight. We notice that the ω z profile of the vortex center has been found quasi-linear with height [26], but why? If it is robust, can we adopt this linear ω z profile as an approximation? With this in mind, we establish an analytical model of asymmetry which uses the vertically-averaged CC/AC vortex center ω z as the asymmetry norm. Three types of boundary conditions are considered to delineate the relative roles of the fluid interior dynamics (vertical advection and stretching of ω z ) and Ekman layer dynamics: both plates no-slip (NN), both free-slip (FF), and a mixed one (NF). A series of direct numerical simulations (DNS) of RRBC are performed to test the theoretical model.
The paper is organized in the following way. Section II introduces the governing equations and the DNS settings. Section III introduces the theoretical model of vorticity asymmetry and its comparison with DNS. Section IV discusses the physical understanding, and section V concludes the paper.

II. THE GOVERNING EQUATIONS AND DNS SETUP
The dimensional variables are denoted with "*". Let x, y, z be the unit vectors of Cartesian coordinate, x * be the position vector, u * = (u * , v * , w * ) be the velocity, p * be the pressure potential, θ * be the disturbance temperature that has subtracted a diffusiveequilibrium linear temperature profile, ω * = (ω * x , ω * y , ω * z ) be the vorticity with ω * = ∇ * ×u * , where ∇ * = x∂/∂x * + y∂/∂y * + z∂/∂z * is the gradient operator. Following [26] The convection is between a warm plate at z = −1/2 and a cold plate at z = 1/2, with a doubly-periodic side boundary. The flow obeys the Boussinesq equation: where ∇ = x∂/∂x + y∂/∂y + z∂/∂z is the nondimensional gradient operator. Note we have convective Rossby number Ro ≡ W/(f H) = E(Ra/Pr) 1/2 . The disturbance temperature boundary condition is Dirichlet: The impermeable velocity boundary condition is: We study three types of tangential velocity boundary conditions, no-slip-no-slip (NN), freeslip-free-slip (FF) and no-slip-free-slip (NF): The NF case may have implications for penetrative convection in the atmosphere where only the lower surface provides frictional drag, i.e., the tropospheric deep convection penetrating into the stably-stratified stratosphere, or the boundary layer convection penetrating into the capping inversion layer [34]. Approximating the stably-stratified upper layer as a free-slip lid is a usual theoretical practice [35,36]. However, Moeng and Rotunno [37] argued that a , κ e = (PrRa) −1/2 .  Here the black "+" denotes a detected and successfully paired updraft plume center, and the white "x" denotes a corresponding downdraft plume center.
part is generally stronger and compacter than the anticyclonic part. As the convective amplitude is relatively small, the vortices are approximately erect, axisymmetric and have approximately variable-separation structure in the vertical and radial direction [26]. Thus, understanding the asymmetry can be split into two tasks: understanding the asymmetry of the vortex center and the radial structure respectively. This paper focuses on the former but also takes the latter into account.
The asymmetry at the vortex center was first revealed by Portegies et al. [26] who made an ensemble of vortex center ω z and w vertical profile to validate their symmetric vortex model. We follow this approach and presented the ensemble average profiles diagnosed from our DNS in Fig. 3 for all the three boundary types. The vortex center detection algorithm used here is introduced in Appendix B. Fig. 3 shows that the ensemble standard deviation of ω z is relatively small and that of w is larger. Both are smaller for the smaller Ro test which has a weaker convective amplitude in the equilibrium state. Note that the ensemble- the horizontal Laplace of vertical vorticity ∇ 2 h ω z = ∂ 2 ω z /∂x 2 + ∂ 2 ω z /∂y 2 . The rank of sharpness is: ∇ 2 h ω z > ω z ∼ w > θ , because ∇ 2 h operator has a sharpening effect. The sharpness difference of θ and ω z is evident from the thermal wind vorticity equation averaged vortex is not a "real vortex" that exactly obeys N-S equation, because the nonlinear advection terms in the governing equation do not commute with the average operator. However, as the budget of ω z equation (shown in Fig. 4) is approximately balanced, we assume the ensemble-average vortex to obey N-S equation and use it as a reference for theoretical modeling. Fig. 1 and Fig. 2 show that the horizontal position of the vorticity extremums generally coincide with those of w and θ , so "vorticity center" and "updraft/downdraft center" will be regarded as the same object in this paper and is named "vortex center". row denotes ω z for updraft, the second row denotes ω z for downdraft, the third row denotes w for updraft, and the fourth row denotes w for downdraft. The solid lines denote the ensemble-averaged profiles from our DNS, the shadow denotes the ±1 standard deviation of the ensemble, and the dashed lines denote those reconstructed from the theory which uses w = ∆ω z w (0) + Ro V w (1) and We use the steady state ω z equation at the vortex center as the main tool: As the vortex center line is approximately vertical, there is no horizontal advection and tilting. Eq. (10) is a powerful constraint of ω z and w, though not enough to determine the full dynamics, especially convective strength. The vortices' small aspect ratio makes horizontal vorticity diffusion much larger than vertical diffusion in fluid interior, so ν e ∂ 2 ω z /∂z 2 can be neglected there [25]. For an updraft, the positive ω z is produced by the stretching of (ω z +f e ) at the lower layer. As the parcel moves up, ω z is dissipated by horizontal diffusion to be close to 0 at the middle plane, and it further decreases to negative value at the upper-layer due to the squeezing of (ω z + f e ). For the small Ro case where ω z f e , the vertical advection and stretching of ω z are small terms, so the primary fluid interior balance is: The remaining balance between ambient vorticity stretching and horizontal diffusion tends to render a vertically antisymmetric ω z and a symmetric w [39]. Should the CC/AC asymmetry be due to the finite Ro effect in Eq. (10)? With this in mind, we further simplify the ω z equation and its boundary condition.
For no-slip boundary, there is an Ekman layer that links the nearly geostrophic interior with large ω z to the zero ω z at the boundary. Vertical diffusion of ω z approximately balances the stretching of background vorticity, yielding a non-growing boundary layer depth For the Ro 1 case, both the vertical advection and stretching of ω z in the Ekman layer are small terms. As Ekman layer is very thin (E 1/2 1), Eq. (12) leads to a parameterization using Ekman pumping relation that links the w and ω z at the Ekman layer top [17,39,40], as if the boundary is penetrative. All of them use the linear Ekman layer theory for Ro 1 case, which renders symmetric CC/AC and does capture the lowest order dynamics. However, the small vertical advection and stretching of ω z also contribute to the asymmetry in Ekman layer, enhancing suction by AC and weakening pumping by CC. Ishida and Iwayama [41] used asymptotic expansion to show that the stretching of ω z strengthens the existing ω z and the vertical advection of ω z weakens it for both CC and AC. As the stretching dominates the advection, their net ω z production shares the burden of the f e stretching for CC but add some for AC. As a result, the vertical integral of |f e ∂w/∂z| within the Ekman layer and therefore the Ekman layer top |w| is smaller for CC than AC. The asymmetry could also be related to the centrifugal acceleration that decelerates (accelerates) the converging (diverging) fluid [42]. We adopt the first order truncation of the asymptotic Ekman pumping formula for a circular vortex derived by Hart [43]: Here z = ±H e /2 denotes the boundary layer top rather than the solid wall. The γ = 1/5 is a constant; b and t are used to denote no-slip and free-slip boundary condition at the lower and upper plate respectively: The derivation and validation of Eq. (13) is shown in Appendix C.  [26] for NN boundary type and has been reproduced by the neutral mode of their linear stability analysis. For FF type, the linear stability analysis [16] shows ω z ∼ ±sin(πz), agreeing with our DNS result in Fig. 3. It is also not far from a linear profile. The NF type profile looks like a mixture of the NN and FF ones. In Appendix D, the quasi-linear shape is explained with the advective-diffusive equilibrium of temperature field and thermal wind relation, as is also illustrated in the schematic diagram in Fig. 5.
This feature inspires us to fit the vortex center ω z (outside of the boundary layer) using a linear function of z with two unknown parameters ω z and ∆ω z : Here w is used to signify updraft or downdraft, ω z is the intercept which denotes the vertical average of ω z at fluid interior, ∆ω z represents the slope which denotes the vertical vorticity magnitude. An asymmetry factor δ which is ω z normalized by ∆ω z is defined to quantify the asymmetry of the vortex center ω z profile: The ∆ω z contains information of the convective strength that cannot be solely constrained by the ω z equation. Geometrically, δ is the shift of ω z = 0 height from the middle plane z = 0. It is also linked to the ratio of AC to CC vorticity magnitude at the Ekman layer top where their contrast is the largest: Our goal is to solve δ. In our work, ∆ω z is directly diagnosed from DNS. To theoretically estimate ∆ω z , one can try the finite-amplitude vortex model which tells the amplitude and wavenumber at the same time [28]. See section III C for further discussion.
Now we consider the radial ω z structure in a vortex which determines the vortex center ∇ 2 h ω z . This involves the CC/AC asymmetry in the radial structure. Fig. 2(c) shows the ∇ 2 h ω z at z = −0.289H e which is within the lower layer. It's not surprising that the ∇ 2 h ω z for updraft vortices is negative, indicating that horizontal viscosity reduces ω z magnitude.
However, the ∇ 2 h ω z of downdraft vortices is close to 0. The highest ∇ 2 h ω z is not at the updraft centers, but surrounding them in a petal shape. The budget of ω z equation at the vortex center in Fig. 4 shows that the updraft core's ∇ 2 h ω z of the reference test is also quasi-linear with height, but it is mostly negative. To explain this, we note that for doubly periodic domain the horizontally integrated ω z at any height should be 0: This makes cyclones which are larger in magnitude have smaller fractional area, and vice versa for anticyclones, as was first noticed by Julien et al. [18] and is illustrated here in the left panel of Fig. 5 as a schematic diagram. We speculate that the contraction (expansion) of vortex area at a CC (AC) might be partly driven by the pushing of the horizontally convergent (divergent) flow. The petal shape ∇ 2 h ω z patch could be the vortex shield of the strong cyclones which has been predicted in the vortex models of Portegies et al. [26] and Grooms et al. [28], and has been experimentally verified by Shi et al. [7]. It is not a closed donut in Fig. 2(c) due to the straining of the neighboring vortices that breaks the axisymmetry. As anticyclones are weak, the cyclones' shield ω z (<0) can outweigh and surround the anticyclones' core vorticity ω z (<0), making ∇ 2 h ω z remain negative. We present an approximation of the ∇ 2 h ω z profile by grasping its upper and lower ends. There, the vortex center vorticity still outweighs the neighbor's shield vorticity, with ∇ 2 h ω z and ω z having opposite sign. At both ends, we let the length scale for ∇ 2 h at a CC (AC) center be L + (L − ). Using Eq. (15), ∇ 2 h ω z is estimated as: where h ω z at the cyclonic (anticyclonic) end: Here α is a nondimensional shape factor. Portegies et al. [26] and Shi et al. [7] have shown that the horizontal structure of ω z in the cellular regime approximately obeys zeroth order Bessel function J 0 (kr) where k is the eigenvalue and r is the radial distance from vortex center. The corresponding Bessel equation yields ∇ 2 h ω z = −k 2 ω z at the vortex center. Let the first zero point (kr ≈ 2.4) of ω z be half of L d , we get kL d /2 ≈ 2.4, and therefore α ≈ 23 which is used in this model. L + and L − depends on CC/AC fractional area which is linked to CC/AC vorticity magnitude via ω z dxdy = 0: Suppose the number density of AC and CC are identical, the vortex distance scale L d can be interpreted as the length scale corresponding to the mean fractional area (proportional to the square of length scale) of AC and CC: Substituting Eqs. (18), (19), and (20) into Eq. (17), and using δ 1, we get an approximate expression of ∇ 2 h ω z which has O(δ 2 ) error: Thus, −∇ 2 h ω z ∼ (2δ − w z/H e ) skews more towards positive than ω z ∼ (δ − w z/H e ), due to the additional effect of the changing radial structure. This indicates that the ∇ 2 h ω z = 0 height is farther from the mid-plane than that of ω z = 0, in agreement with the DNS result which is illustrated in the budget plotting in Fig. 4 for updrafts. Note that our argument in deriving Eq. (21) requires the symmetry between an updraft and downdraft, which is not strictly met for the NF type. The error is not significant in Fig. 4(c) and (f), but it is indeed an unresolved flaw.
Substituting Eqs. (15) and (21) into Eq. (10), the vertical vorticity equation is simplified to a first order ODE about vortex center w with one unknown parameter ω z (or δ) and two given parameters ∆ω z and L d which are diagnosed from DNS. This, together with the two boundary conditions which is obtained by substituting Eq. (15) into Eq. (13), render a deterministic problem: One might speculate: how can a first order ODE satisfy two boundary conditions? In fact, such an overdetermined problem can be avoided by setting δ a proper value, which renders a nonlinear eigenvalue problem. In other words, the δ which stands for CC/AC asymmetry can be solved from the solvability condition, as is derived in section III B below. The asymmetry is produced by both the ODE which depicts the fluid interior dynamics and the Ekman layer boundary conditions.
B. Solving the asymmetry factor δ Because the vortex magnitude which is measured by ∆ω z is an input parameter for the asymmetry problem, w and ω z are rescaled with ∆ω z , yielding the new variables w + and ω + z : They will be used in the theoretical derivation in this section. Eq. (22) becomes: Here Ro V is vorticity Rossby number which denotes vorticity magnitude, and Γ is the vortex aspect ratio: Here the superscripts (0) , (1) and (2) denote the terms with different Ro V orders. Substituting Eq. (26) into Eq. (24), we get the zeroth order ODE and boundary condition: It describes the balance between background vorticity stretching and the horizontal diffusion of vorticity. δ (0) is solved by enforcing the solvability condition: Equation (28) The αE/Γ 2 essentially denotes the relative strength of horizontal diffusion to background vorticity stretching. The last term denotes the extra vertical velocity provided by Ekman pumping. For NN or FF boundary condition where δ (0) = 0, w (0) is a symmetric parabolic line which has the largest magnitude at the middle plane z = 0.
The first order equation and boundary conditions are It describes how the vertical advection (the left-hand-side) and stretching of ω z (the first term on the right-hand-side, RHS) are balanced by the asymmetry in ω z and w through the stretching of f e (the second term on the RHS), horizontal diffusion (the third term on the RHS) and the anomalous Ekman pumping/suction. The solvability condition leads to δ (1) which is generally positive: The last term in the numerator is the contribution from the nonlinear effects in Ekman layer which makes pumping less efficient in CC than suction in AC. The corresponding w (1) is: where the coefficients a, b and c are: Fig. 6 shows the w (0) and w (1) for Ro = 0.158 case (EXP 4) of the three boundary types.
w (1) is an odd function for NN and FF type, which is obvious from a symmetry analysis of Eq. (30). The w (1) tends to shift the extremum of w towards the anticyclonic side, agreeing with the DNS in Fig. 3. However, Fig. 3 also shows that the theory underestimates the magnitude of w for EXP 4 (large Ro test), as well as the asymmetric feature of w as a whole.
As will be shown in the validation against the numerical nonlinear solution in Fig. 8, the first order asymptotic solution renders sufficient accuracy, even for Ro V 1 cases. We conclude that the difference in boundary condition leads to zeroth order vorticity asymmetry, and a finite Ro V leads to first order asymmetry.
Here we introduce the method of diagnosing Ro V = ∆ω z /f e , δ = ω z /∆ω z and Γ=L d /H e from DNS. The ∆ω z and ω z are diagnosed from the ensemble-averaged vortex center ω z profile. Their values, as well as Ro V and δ are diagnosed separately for updrafts and downdrafts. The ∆ω z is set as the difference between the maximum and minimum ω z on the profile. The ω z is set as the numerical vertical average of the ω z profile. For convenience, it includes the Ekman layer which has tiny contribution to the integral due to its thin depth. and NF type lying between them. To diagnose L d and therefore Γ, we follow Sakai (1997) to count the vortex number, as is shown in Appendix B. Fig. 7 [16]. Physically, a stronger background rotation leads to a stronger tilting of f e to the tangential direction, which requires a larger baroclinic torque to maintain thermal wind balance. This leads to a smaller vortex size. The dependence on Ra is likely due the nonlinear upscale growth (e.g. vortex merger), where the horizontal scale is larger for a larger convective amplitude. The Γ for NN type is a bit larger than the FF and NF types, probably due to the difference in boundary condition and the convective amplitude.
Except for the Γ ∼ E 1/3 scaling, we will not invoke other theoretical determination of Γ and Ro V in this paper. Though Sakai [21] provided a scaling theory of vortex size, it is suitable for more vigorous convection and requires a well-developed thermal boundary layer.
The finite-amplitude model of Grooms et al. [28] is also a candidate for determining Γ and  7. (a) The dependence of Ro V on Ra for all the data. The "*" denotes NN boundary type, the "+" denotes FF type, and the "o" denotes NF type. The red marks denote ascending (warm) vortices, and the blue marks denotes descending (cold) vortices. (b) The mean vortex aspect ratio Γ of updraft and downdraft vortices for changing E (EXP 1-10) tests. The blue "*" denotes the Γ diagnosed from DNS for NN type; the blue "+" denotes FF type; the "o" denotes NF type.
The solid blue line denotes the E 1/3 law which is the most unstable wavelength predicted with the linear stability analysis of FF type: λ c = 2π/k c = 2 7/6 π 2/3 E 1/3 [16]. For a checkerboard vortex pattern, Γ should be compared to λ c /2. (c) The same as (b), but for Γ of the changing Ra tests . The λ c is not plotted.
Ro V , but it does not yield an explicit expression, therefore we leave the coupling of their model to our vorticity asymmetry theory for future work. is no statistical difference between the updraft and downdraft vortices, due to the symmetry in boundary condition. For both types, the DNS agrees with our prediction that δ ∝ Ro V , even when Ro V is as large as 1. The DNS confirms that the δ of NN type is a bit larger than FF type, as will be explained in section IV. This result is robust for different sampling time, though in general the difference between the FF and NN types is hard to identify at Ro V < 0.5. This is probably technically due to the uncertainty in vortex detection, or theoretically due to the slightly different shape of the ω z vertical profile that blurs the δ−Ro V relation. For NF type, DNS not only confirms the quasi-linear δ − Ro V trend, but also the prediction that δ is smaller for updraft and larger for downdraft by roughly a constant. The slope (δ (1) part) is captured well, but the theory seems to underestimate |δ (0) |. The deviation probably lies in Eq. (20) which requires the symmetry of updraft and downdraft but is not met for NF type. Note that as ω z is not a perfectly linear profile, different measures of ∆ω z and ω z can lead to different δ magnitude. For example, if δ is diagnosed as the height where the ensemble-averaged ω z curve crosses ω z = 0, the δ magnitude is systematically larger than the current method by around 50%, but the trend does not change.
The model is further validated by comparing the ω z equation budget of the DNS data and theory, as is shown in Fig. 4. The agreement is good within the fluid interior (outside of the boundary layer) where the theory works, and is even better as Ro goes smaller.

IV. UNDERSTANDING δ FOR THREE TYPES OF BOUNDARY CONDITIONS
In this section, we try to interpret the expression δ ≈ δ (0) +Ro V δ (1) for the three boundary types: δ FF , δ NN and δ NF .

A. Lower free-slip and upper free-slip (FF)
The asymmetry in FF type solely depends on the interior dynamics. Here δ is proportional to Ro V , and does not depend on E and Γ: The w (0) is proportional to αE/Γ 2 , so the vertical advection and stretching of ω z , which are asymmetric factors and are also associated with w (0) , are also proportional to αE/Γ 2 . In the first order equation (30), such asymmetry is balanced by the asymmetric horizontal diffusion which is also proportional to αE/Γ 2 . Thus, the αE/Γ 2 in δ (1) cancels out. Physically, the asymmetry is due to the larger absolute vorticity stretching in the cyclonic region that favors positive vorticity production, as well as the vertical advection of ω z that "pushes" the ω z profile toward the anticyclonic side to produce a shift of δ. That the vortex becomes thinner at the cyclonic side provides a stronger horizontal diffusion to balance the extra positive vorticity production there, helping reduce the asymmetry.

B. Lower no-slip and upper no-slip (NN)
The δ here clearly shows the dominant role of the fluid interior contribution (Ro V /12 ≈ 0.083Ro V ) which is already seen in the FF type, and a smaller correction (0.024Ro V for EXP 4) due to the Ekman layers which introduces the dependence on E and Γ: For given Ro V , δ NN is always larger than δ FF . On one hand, the lower and upper Ekman layer systematically raise the updraft strength, which enhance the vertical advection of ω z .
Such an asymmetric tendency is partly balanced by the horizontal diffusion whose strength is proportional to Γ −2 , so a thinner vortex with a smaller Γ has a lower δ NN .
On the other hand, the stronger vorticity and therefore the Ekman pumping anomaly at the cyclonic side suppresses the asymmetry (the 2(E/2) 1/2 part in the denominator of the correction term), because it makes the fluid interior w increment which is responsible for stretching smaller than the w decrement for squeezing. This effect also has its limitation. For a cyclone, the centrifugal effect (γ) reduces the efficiency for pumping mass with given vorticity, so the difference between the w increment near the lower boundary and the decrement near the upper boundary is mitigated.
The correction term can be simplified by substituting in the Γ ∼ E 1/3 law, and dividing both the numerator and denominator by E 1/2 . Now, only the first term in the denominator has a E −1/6 dependence on E. As "−1/6" is quite smooth a slope, the correction term is quite close to a constant.
In summary, the gross effect of the two Ekman layers is enhancing asymmetry. The direct effects of Γ and E in Eq. (35) are much less influential than Ro V .

C. Lower no-slip and upper free-slip (NF)
One may consider the NF type as an extreme case of a differential Ekman pumping/suction of the lower and upper plates, as is mentioned in analyzing the NN type. Here the approximate expression of δ under δ (0) 1 is: where δ (0) can be simplified using E ∼ 10 −4 and Γ ∼ 0.2, as well as the Γ ∼ E 1/3 law : Thus, δ (0) is a quite fixed quantity, with the EXP 4 case rendering δ (0) ≈ −0.03 w . The δ (0) is positive for downdraft and negative for updraft. To explain this, first consider an updraft.
The bottom Ekman pumping gives the parcel some w (0) , but when it approaches the upper plate, w (0) needs to decrease all the way to 0. As a result, the stretching of background vorticity is weaker than squeezing, which makes its ω z more negative. Oppositely, there is extra stretching for downdraft, which makes its ω z more positive.
As for δ (1) which is in the large bracket of Eq. Another possibly more accurate way to calculate δ (0) is solving the neutral mode of the linear stability problem [16,26] with NF type boundary condition, but an explicit expression is unlikely to be found. This is left for future work.

V. CONCLUSION
We study the asymmetry of cyclones and anticyclones for rotating convection between two parallel plates, specifically in its "cellular regime" where convection is weak and rotation is strong. The rotation shapes convection into densely packed vertically-aligned vortices, which are of comparable magnitude and in quasi-steady state. The updraft (downdraft) vortices are cyclonic (anticyclonic) at the lower level and anticyclonic (cyclonic) at the upper level.
Portegies et al. [26] derived a linear model of convective vortex between two no-slip plates.
They showed that the vertical vorticity ω z at the vortex center is quasi-linear with height.
A clear shift of the vortex center ω z profile toward the cyclonic side (i.e. the magnitude of positive ω z is larger than the negative ω z ) was observed but remains unexplained. In this work, a simple model is presented to explain such asymmetry in the cellular regime.
The observation that vortex center ω z profile is quasi-linear with height is explained with the thermal wind relation and the advection-diffusion of temperature. This constraint enables us to focus on ω z equation. The complex dynamic-thermodynamic coupling is only used for determining the convective strength which is characterized by the vorticity difference ∆ω z between the cyclonic and anticyclonic ends. The parameter ∆ω z is considered to be known a priori, and is represented as vorticity Rossby number Ro V = ∆ω z /f e , which can be either diagnosed from direct numerical simulation (DNS) or potentially from the existing theoretical models of convective amplitude [21,28]. The parameter Ro V , together with the vortex aspect ratio Γ and Ekman number E, renders an analytical model of the asymmetry factor δ = ω z /∆ω z (ω z is the vertically averaged ω z ) based on the balanced dynamics of ω z alone.
Physically, the asymmetry in the cellular regime is contributed from both the fluid interior and Ekman layer dynamics. As for the fluid interior within a vortex column, the downgradient vertical advection of ω z "pushes" its profile toward the anticyclonic side to enhance the cyclonic part. Meanwhile, the stretching of absolute vertical vorticity (ω z + f e ) is more efficient at the cyclonic side. Such positive ω z production makes the cyclonic side stronger and compacter. A compacter cyclone has a stronger horizontal diffusion of vorticity, so the feedback on vortex shape offsets some of the asymmetry. As a result, we get δ ≈ Ro V /12 for a vortex between two free-slip plates (FF type), without the influence of Ekman layer. Such quasi-linear dependence of δ on Ro V is valid even when Ro V is as large as 1. The boundary layer effect is studied for a vortex between two no-slip plates (NN type).
It turns out to be a minor modification to the asymmetry problem compared to the interior dynamics. The Ekman pumping (suction) at the cyclonic (anticyclonic) side elevates the vertical velocity strength, which enhances the interior asymmetry produced by the vertical advection of ω z . However, the pumping and suction strength difference also matters. Because a cyclone is generally stronger, the stronger pumping (weaker suction) leaves less w increment for vorticity stretching than squeezing in fluid interior, so the original cyclonic tendency is offset a bit. Another difference is due to the relatively weak nonlinear vertical advection and stretching of ω z in the Ekman layer, which makes cyclones pump less ver-tical velocity than anticyclones, and therefore slightly favors the asymmetry. An extreme case is the vortex between a no-slip and a free-slip plate (NF type). The stretching and squeezing outside of the Ekman layer are asymmetric even for the small-amplitude problem where Ro V → 0. This makes the updraft (downdraft) stretches (squeezes) less and squeezes (stretches) more, and therefore renders a baseline negative (positive) δ that depends on the aspect ratio Γ and E but not Ro V . We argue that the previous asymmetry explanations [1,6,18,30,32,33], which emphasize the role of Ekman pumping in generating vorticity at the cyclonic side and the turbulent mixing in diluting vorticity at the anticyclonic side, is not suitable for the non-turbulent cellular regime.
The DNS result is consistent with our theory of δ, though a little underestimation of the baseline δ magnitude for NF type (named δ (0) ) exists. The DNS agrees with the theoretical result that w peaks at the anticyclonic side for NN and FF types, but the theory underestimates the magnitude of this asymmetric feature. One flaw of the theory is the treatment of vorticity horizontal diffusion for NF type.
In future, whether the theory of δ is robust for different Prandtl number Pr needs to be checked, though Pr does not appear in the analytical model. The theory might be extended to consider non-Boussinesq effects such as the variable density (anelastic) effect [44,45] and the dependence of viscosity on temperature [46]. In a macroscopic view, the modeling framework is only a small step toward understanding cyclone/anticyclone asymmetry, because it only applies to the non-turbulent quasi-steady regime where the asymmetry factor δ and vorticity Rossby number Ro V are well-defined. By linking δ to other more general asymmetry norms such as skewness, a connection to the general asymmetry problem in more turbulent regimes would be possible.  Fig. 8 shows that the average distance between a cold and warm vortex is at least 0.15H e which corresponds to 12 grid points, so the horizontal resolution is enough. The vertical mesh generation function is: The CM1 code is modified from its configured "Rayleigh-Bénard convection case". Specifically, the following subroutines are modified: the vertical mesh is set in "param.F", the buoyancy expression is set in "solve.F" and the basic state profile is set in "base.F". The code has been benchmarked with critical Rayleigh number test for the FF type, using EXP 14 -EXP 16. Firstly, for each slice we use the Q criterion to circle out some "vortical area" and find all the vorticity extremum within it. The vortical region is defined as the place where the horizontal velocity gradient tensor ∇ h u h = ∂ i u j , i, j ∈ {1, 2} has complex eigenvalues [1,27]. This is equivalent to finding the Q 2D > 0 place where: Here "det" means the determinant, "Tr" means the trace.
Secondly, within Q 2D > 0 region, we pick out all the maximum ω z position as candidate CC centers and minimum ω z position as AC centers. If the distance of two vortices is smaller than 0.0375, it is judged as an overlap event and one of the vortex is discarded.
Thirdly, the candidate vortex centers at the two horizontal slices are paired: a CC (AC) at the lower slice and an AC (CC) at the upper slice constitute an updraft (downdraft). The criterion for pairing is that the two vortices' distance projected onto the horizontal plane is smaller than 0.06H e which is a substantial portion of the vortex width. The detected vortex centers that fail to meet the paring criterion are not used for generating profiles.
Fourthly, for each vortex pair, the vortex center line is chosen as the straight line that crosses their centers. Then the w, ω z and ∇ 2 h ω z are sampled on these lines. A byproduct of vortex detection is the vortex number density, which is used by Sakai [21] to calculate the characteristic distance between an updraft and downdraft L d . On each horizontal slice, L d is related to the total vortex number N from the second step (updraft plus downdraft, and no need to be paired) via: Here L 0 = 2.5 is the domain width. In data processing, L d is further calculated as the average between the L 0 /N 1/2 at the lower and upper slices.
The ensemble-averaged vortex profile of Portegies et al. [26] at Ra = 2.5 × 10 6 , E = 10 −4 and Pr = 1 shown in their Fig. 6 is compared with the detection result of EXP 4 in Fig. 3(a) and (g). The agreement is good, so both the detection method and the DNS solver with no-slip boundary condition are validated.
the first order, it is: Here h E = (2ν e /f e ) 1/2 = (2E) 1/2 H e is the characteristic depth of Ekman layer, V is the fluid interior tangential velocity, w b is the pumping vertical velocity at an infinite height and r is the distance to the vortex center. We are only interested in the relation at the vortex center.
The V and its derivatives are replaced by vorticity via assuming the vorticity radial profile to be zeroth order Bessel function J 0 (kr) which is suitable for the columnar vortices in the cellular regime of RRBC [26]. As the flow close to the vortex center can be approximated as rigid body rotation, we have: As ∂ω z /∂r ∼ dJ 0 (kr)/dr| r=0 = −kJ 1 (0) = 0, there is: Substituting Eqs. (C2) and (C3) into Eq. (C1), we get: Here γ = 1/5 is a constant that signifies the nonlinear effect. The derivation for the upper plate Ekman layer is similar, with a flip of sign at the right-hand-side of Eq. (C4).
The formula for the bottom Ekman layer is validated with Ro = 0.158 (EXP 4) and Ro = 0.132 (EXP 8), NN boundary type simulations, using vortex center ω z and w profile Here we validate Eq. (C4) with the DNS result. Fig. 9 shows that w and ω z keep changing with height above the Ekman layer (especially for downdraft vortices), so a strict comparison with formula is hard. We choose the w and ω z data at the height of maximum |ω z | (within 3.4h E , or 19 grid points from the bottom) for comparison, as is shown in Fig. 10. This height is chosen because the ∆ω z is calculated as the difference of ω z maximum and minimum. and hydrostatic balance which is not fully satisfied at updraft centers in the fluid interior [39]: Note that if the centrifugal acceleration is taken into account (gradient wind balance), the right part of Eq. (D1) is only valid right at the vortex center where tangential velocity is zero. As warm core updrafts and cold core downdrafts surround each other, the ∇ 2 h θ depends on the θ difference between a pair of neighboring updraft and downdraft, as well as the characteristic vortex distance L d : where θ ↑ (z) denotes the perturbed updraft center temperature profile and θ ↓ (z) denotes that of a downdraft. Meanwhile, as the updraft and downdraft are governed by the same dynamics for a symmetric boundary condition such as FF and NN (and a bit less for NF), it yields: Substitute Eqs. (D2) and (D3) into Eq. (D1): (D4) Thus, ∂ω z /∂z should be close to an even function of z. In an updraft, θ ↑ (z) is positive at all levels due to the upward transport of warm fluid, so ∂ω z /∂z is negative at all z, and vice versa for a downdraft.
Next, we explore θ structure at small and large w limits more closely. The vortex center temperature equation is a steady advection-diffusion equation: Here the vertical diffusion has been neglected due to the small aspect ratio of the vortex. The Here W e is the effective vertical velocity scale, which can be estimated with the balance of background vorticity stretching and the horizontal diffusion of vorticity in the fluid interior: Here α = 23 is the same as that used in Eq. (21), and ∆ω z is the vorticity difference between the CC/AC ends. Substituting Eq. (D7) into Eq. (D6), we get: Pe ∼αRo V Pr.
For cellular regime, Ro V O(1), so Pe O(10 1 ) for our Pr = 1 case, thus advection is dominant unless Ro V < 0.1 where convection is very weak.
The small Pe limit which corresponds to Ro V 1 and small Ra is a small-amplitude regime, where the normal mode solution of Portegies et al. [26] works. The main balance is between the vertical advection of basic state temperature and horizontal diffusion: Thus, the vertical profile of θ should be proportional to w whose largest magnitude is at the middle layer. Suppose w ∼ cos(πz) in this linear regime, Eq. (D7) shows that the vertical variation of ω z at small w limit is controlled by sin(πz). If there are Ekman layers, w will be more uniform in the fluid interior because Ekman pumping and suction have systematically elevated the magnitude of w. The more vertically uniform w is, the more vertically linear θ and therefore ω z should be. The solution to Eq. (24) that only satisfies the lower boundary condition is: where the integral constant C is: