{'fraction_non_alphanumeric': 0.028493243433859002, 'fraction_numerical': 0.008350495105368578, 'mean_word_length': 4.515481562159639, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Mean flow generation in rotating anelastic two-dimensional convection

We investigate the processes that lead to the generation of mean flows in two-dimensional anelastic convection. The simple model consists of a plane layer that is rotating about an axis inclined to gravity. The results are two-fold: firstly we numerically investigate the onset of convection in three-dimensions, paying particular attention to the role of stratification and highlight a curious symmetry. Secondly, we investigate the mechanisms that drive both zonal and meridional flows in two dimensions. We find that, in general, non-trivial Reynolds stresses can lead to systematic flows and, using statistical measures, we quantify the role of stratification in modifying the coherence of these flows.


I. INTRODUCTION
Geophysical and astrophysical flows are often turbulent and characterised by the presence of a wide-range of temporal and spatial scales. It is often the case that systematic large-scale flows that vary on long timescales co-exist with shorter lived turbulent eddies. Famous examples of such large-scale flows are the zonal jets so visible at the surface of the gas giants 1,2 , the turbulent jet stream of the Earth 3 and the strong zonal and meridional flows in the interior of the Sun which are interpreted as the observed differential rotation and meridional circulations 4 .
A central question for fluid dynamicists is then to determine the role of the smaller-scale turbulent flows in modifying the systematic (or mean) flows. In some cases the turbulence will act purely as a dissipation or friction and act so as to damp any mean flows that occur. However in certain circumstances (usually when rotation is important) the turbulence may act as an anti-friction 5 and play a key role in driving and maintaining the mean flows. This mechanism is complicated by the tendency of the mean flows to act back on the turbulence and modify its form; a process that often determines the saturation amplitude of the mean flows. These are complicated interactions and although much progress has been made (as discussed briefly below) there is still much that is not understood.
Here we focus on a simple model where the turbulence is driven by a thermal gradient leading to convection. Owing to its importance in planetary and stellar interiors, this has been an extremely well studied problem, with many theoretical and numerical studies performed both in Cartesian (local) [6][7][8][9][10] and spherical (both shell and full sphere) geometries [11][12][13][14][15][16][17][18] . We note that, in order to capture some effects of vortex stretching in a spherical body, Busse 19 introduced an annulus model for its relative simplicity. This geometry has been used in attempts to model the zonal flow on Jupiter. For example, Jones et al. 20 used a rotating annulus model in a two-dimensional (2d) study and incorporated the possibility of boundary friction which allowed for the more realistic multiple jet solutions to be found more easily. Rotvig and Jones 21 examined this annulus model more extensively and identified a bursting mechanism that occurs in the convection in some cases.
We wish to focus on the role of stratification in altering the dynamics of the mean flow and so focus on a simple plane layer model in two dimensions allowing us to access some parameter regimes more easily than in other, more complicated geometries. The plane layer model, when the axis of rotation is allowed to vary from the direction of gravity, can be used to represent a local region at different latitudes of a spherical body. This paper builds on previous, largely Boussinesq, studies in a Cartesian domain which are summarised here -though this is by no means a complete review.
Much of the previous work using a local model makes the Boussinesq approximation so that density variations are neglected except in the buoyancy force 22,23 . For example, Julien and Knobloch 9 studied nonlinear convection cells in a rapidly rotating fluid layer with a tilted rotation vector. They performed an asymptotic analysis and concluded the orientation of the convection rolls affects the efficiency of mean flow generation. Hathaway and Somerville 6 performed three-dimensional (3d) simulations of convection in layers with tilted rotation vectors and no slip boundary conditions. They found that the horizontal component of the rotation vector gives a preference for cells aligned with the rotation axis which produces dynamical changes that drive a mean flow. There have also been a number of studies analysing the role of shear flow in a rotating plane layer. Hathaway and Somerville 8 considered Boussinesq convection in a rotating system with an imposed shear flow. The imposed flow was constant in depth, but varying in latitude. They found in the non-rotating case that the convection extracts energy from the mean flow and reduces the shear, but in the rotating case the convection can feed energy into the mean flow and increase the shear. Hathaway and Somerville 7 studied the interaction between convection, rotation and flows with vertical shear. They performed 3d simulations with no slip boundary conditions and found that in cases with vertical rotation the convection becomes more energetic by extracting energy from the mean flow. However, for cases with a tilted rotation vector the results depend on the direction of the shear. Saito and Ishioka 10 revisited the problem of the interaction of convection with rotation in an imposed shear flow. They were able to examine a larger region of parameter space than Hathaway and Somerville 8 and identified a feedback mechanism in which the convection interacts with the rotation in such a way that leads to an accelerated mean flow.
This mechanism operates when the sign of the shear flow is opposite to the vertical component of the rotation axis and relies upon the sinusoidal form of shear flow they imposed. Currie 24 studied the generation of mean flows by Reynolds stresses in Boussinesq convection both in the absence and in the presence of a thermal wind and showed whether convection acts to increase or decrease the thermal wind shear depends on the fluid Prandtl number and the angle of the rotation vector from the vertical.
In many astrophysical systems the fluid is strongly stratified (e.g., in stellar and some planetary interiors) and therefore density changes across the layer (that are neglected by the Boussinesq approximation) may play a significant role. Therefore, there exist studies where fully compressible convection has been simulated [25][26][27] but the full equations are computationally intensive to solve owing to the necessity of accurately tracking sound waves. However, for systems where there are a large number of scale heights involved but that remain close to being adiabatic, the anelastic equations are an improvement on the Boussinesq equations [28][29][30] . Furthermore, the anelastic equations allow for density stratification across the layer whilst still filtering out fast sound waves thus making studying a compressible layer more computationally accessible.
In this paper we use the anelastic approximation to investigate the effects of stratification on both linear and nonlinear convection, focussing on the role stratification plays in altering the dynamics of mean flows. In other words, we extend the work of Hathaway and Somerville 6 to include the important effects of stratification.
By considering density variations across the fluid layer one cannot take both the dynamic viscosity µ and the kinematic viscosity ν of the fluid to be constant (since µ = ρν where ρ is the fluid density) therefore, a choice has to made. Similarly, one cannot take both the thermal conductivity k and the thermal diffusivity κ of the fluid to be constant and again, a choice has to be made.
The results can depend on these choices (see, e.g., Glatzmaier and Gilman 31 ). In this paper we consider a formalism where ν and k are constant. Anelastic formalisms also differ depending on whether entropy or temperature is diffused in the energy equation 30 . If one diffuses entropy then temperature can be eliminated as a variable from the formulation; for simplicity, this is often the approach taken in nonlinear studies and is the choice we make here. Further differences between anelastic formalisms arise from whether one takes the vertical axis to be parallel or antiparallel to gravity. Typically, in Boussinesq formulations the vertical axis increases upwards, whereas in compressible studies it is taken downwards. In order to ease the comparison with Boussinesq models, we take the vertical axis to increase upwards in line with Mizerski and Tobias 32 . The existing studies of anelastic convection, linear and nonlinear, each consider a different formalism (discussed below) and therefore care has to be taken when comparing across the different models.
The onset of compressible convection in a local Cartesian geometry, using the anelastic approximation, has been studied in a number of papers. The earliest of these studies includes the work of Kato and Unno 33 who studied the onset of convection in an isothermal reference atmosphere.
As an alternative, it is common to assume a polytropic reference atmosphere, as we do here. For example, Jones et al. 34 considered a polytropic, constant conductivity model in a Cartesian geometry in which rotation, magnetic field and gravity were taken to be mutually perpendicular. This is different to the configuration in this paper as we consider rotation that is oblique to gravity.
More recently, Mizerski and Tobias 32 investigated the effect of compressibility and stratification on convection, using the anelastic approximation, in a rotating plane layer model where rotation and gravity were aligned; they compared a model with constant k with one with constant κ whilst keeping ν constant throughout and they chose to diffuse entropy. In this paper we essentially consider the constant conductivity case of Mizerski and Tobias 32 but we allow the rotation vector to tilt from the vertical in order to model behaviour at mid-latitudes of a spherical body. To date, the only other study examining the linear stability of compressible convection in this tilted f-plane geometry is that of Calkins et al. 35 who compared the onset of convection in compressible and anelastic ideal gases on a tilted f-plane. The majority of the results from that paper are concerned with comparing the anelastic equations with temperature diffusion with the fully compressible equations, and so those results are not directly comparable with ours. They do however, also consider an entropy diffusion model but in contrast to the work here, they take µ constant and have the vertical axis pointing downwards.
In contrast to the Boussinesq systems, there have been relatively few studies of mean flow generation by convection under the anelastic approximation in a local model. Rogers and Glatzmaier 36 did model penetrative convection in a system where a convective region is bounded below by a stable region and Rogers et al. 37 presented 2d simulations of turbulent convection using the anelastic approximation in a non-rotating system. However, most relevant to the work we undertake here is the study of Verhoeven and Stellmach 38 who used 2d anelastic simulations of rapidly rotating convection in the equatorial plane (gravity and rotation perpendicular) to support a compressional Rhines-type mechanism predicting the width of jets driven in such a rapidly rotating system. In their study, as we do, Verhoeven and Stellmach 38 take ν constant and consider an entropy diffusion model, however, in contrast to our work, Verhoeven and Stellmach 38 take κ constant.
Our study, using a local model, allows us to focus on the effect of stratification on mean flow generation without many of the complicating features of a spherical geometry, say. However our model still incorporates some important physical features that are expected to play a role in the determining the dynamics of many astrophysical objects, namely stratification on the f-plane.
This paper is organised as follows: in section II, the model and governing equations are presented. In section III, we aim to add to the existing literature discussed above by carrying out 3d linear stability analysis and present a new symmetry of which technical details are given in an appendix. In section IV, we analyse the role of stratification in mean flow generation at mid-latitudes using a 2d local model, focussing on the differences from Boussinesq models. This work complements previous Boussinesq and fully compressible studies in a Cartesian domain, as well as the global anelastic models considering mean flows.

II. MODEL SETUP AND EQUATIONS
We consider a local plane layer of convecting fluid rotating about an axis that is oblique to gravity, which acts downwards. The rotation vector lies in the y-z plane and is given by Ω = (0, Ω cos φ, Ω sin φ), where φ is the angle of the tilt of the rotation vector from the horizontal, so that the layer can be thought of as being tangent to a sphere at a latitude φ. In this case, the z-axis points upwards, the x-axis eastwards and the y-axis northwards.
We denote the fluid density, pressure, temperature, entropy and velocity by ρ, p, T , s and u = (u, v, w) respectively. The anelastic equations are then found by decomposing the thermodynamic variables into a reference state (denoted by the subscript ref ) and a perturbation. The reference state depends on z only and is assumed to be almost adiabatic. The departure from adiabaticity is measured by a small parameter, ǫ given by where d is the layer depth, the pressure scale height, c p is the specific heat at constant pressure, g is the acceleration due to gravity, the subscript ad indicates the value for an adiabatic atmosphere and a subscript r denotes a value taken at the bottom of the layer, z = 0.
ǫ will also be a measure of the relative magnitude of the perturbations and we assume so that the perturbations are small compared to the reference state.
It is useful to write the equations in a dimensionless form, using d as the length scale and the thermal diffusion time d 2 κr as the time scale, where κ is the thermal diffusivity. In this case, the leading order equations are given by where ς ij = ρ ref θ is the dimensionless superadiabatic temperature gradient and γ is the ratio of specific heats at constant pressure to constant volume. In this case the dimensionless parameters are the Rayleigh, Taylor and Prandtl numbers given by respectively. Again, note in our formalism that we assume the kinematic viscosity ν and the turbulent thermal conductivity k = ρ ref c p κ (where we now interpret κ as the turbulent thermal diffusivity) to be constant. Furthermore, we use a model that takes the turbulent thermal conductivity to be much larger than the molecular conductivity and so equation (5) contains an entropy diffusion term but not a thermal diffusion term; Braginsky and Roberts 30 discuss models including both terms. In addition, we have employed a technique introduced by Lantz 39

and Braginsky and
Roberts 30 to reduce the number of thermodynamic variables 29,40 .
In this paper, we consider a time-independent, polytropic reference state given by where m is the polytropic index and −1 < θ ≤ 0. We note in this model there is no adjustment of the reference state by any mean that may be generated.
The equations (3)-(7) are similar to those given in case (1) of Mizerski and Tobias 32 ; the key difference here however is the introduction of a tilted rotation vector, Ω = (0, cos φ, sin φ). In this formalism, the anelastic equations (3)-(7) reduce to the Boussinesq equations in the limit θ → 0 and so θ can be thought of as a measure of the degree of compressibility. However, a more intuitive measure is given by As χ is increased, the density contrast across the layer is also increased. Furthermore, it is common in the literature to use N ρ , the number of density scale heights, as a measure of the stratification.
Note, θ, χ and N ρ are related by the following relations: With the reference state specified, we solve the anelastic equations subject to impenetrable, stress free and fixed entropy boundary conditions.

III. LINEAR THEORY
Here we perform linear stability analysis of the 3d system in order to obtain information about the preferred modes of convection at onset. We linearise the anelastic equations (3)-(5) by perturbing about a trivial basic state and neglecting terms quadratic in the perturbation quantities. We then assume the perturbations are given by Here k and l are the x and y wavenumbers respectively such that a 2 = k 2 + l 2 and σ = σ R + iσ I = σ R + iω is the complex growth rate. S(z) is the amplitude function for s, W (z) is the amplitude function for w and Z(z) is the amplitude function for the vertical vorticity ζ = ∂v ∂x − ∂u ∂y . Taking the z-component of the curl and double-curl of the momentum equation to eliminate the pressure, along with the entropy equation, results in the following equations for the amplitude functions W , Z and S: σZ =T a 1 2 P r sin φ DW + mθ 1 + θz W + cos φilW to be solved subject to the boundary conditions: Note, to simplify notation, we have used D n to mean d n dz n . We solve this linear eigenvalue problem using the built-in bvp4c solver of Matlab.
Since we assume a reference state that is close to being adiabatic, we use m = 1.495 in all calculations. Berkoff et al. 41 demonstrated that the anelastic approximation gives a good approximation to fully compressible calculations even when the reference state is super-adiabatic, finding a 2% error even when ǫ ∼ 10; but as mentioned above, we prefer to remain close to the adiabatic state, where m = 1.5. . NS convection rolls are those whose axis is aligned in the y-direction (l = 0) and similarly, east-west (EW) convection rolls are those whose axis is aligned in the x-direction (k = 0). These are distinct when φ = π 2 , i.e., when rotation is oblique to gravity. For NS rolls, we observe that solutions with a positive frequency and those with a negative frequency are distinct unlike in the Boussinesq case 42 -we will examine this symmetry breaking in more detail in the next section. For small T a, the negative branch is preferred but this changes to the positive branch as T a is increased. For the cases shown, the positive solution always has the smaller critical wavenumber and critical frequency. In addition, the minimum Ra crit occurs for N ρ > 0, with both the minimum Ra crit and the N ρ at which it occurs, increasing with T a. Another feature of figure 1 (a)-(c) is the difference between the positive and negative frequency solutions; this can be seen in the larger variation of k crit with increasing N ρ for the positive case. In a similar study, Calkins et al. 35 found that the critical Rayleigh number is a monotonic function of stratification. This difference could arise because of a number of different modelling choices. First, as discussed in the introduction, Calkins et al. 35 take constant µ where we have constant ν. Second, the parameter regime considered is very different to that covered here; they choose to focus on the rapidly rotating case (T a ∼ 10 10 ) at P r = 1. Finally, and perhaps most importantly, the definition of Ra is different between our models (the difference resulting from the fact that we take z pointing upwards and ν constant in contrast to Calkins et al. 35   is no distinction to be made between the solutions with positive and negative frequency as they have the same critical values, hence we only plot the positive frequency solutions (we discuss this hidden symmetry in more detail in the next section). The behaviour is very similar to that in the NS case, but Ra crit is higher in the EW case, so that NS rolls are preferred, in agreement with Calkins et al. 35 and also with linear Boussinesq systems, e.g., Hathaway et al. 43 .

B. Symmetry considerations
As highlighted in the previous section, when χ = 1 (N ρ = 0) and l = 0 (NS rolls), there is a distinction to be made between solutions with a positive critical frequency and those with a negative critical frequency. However, when k = 0 (EW rolls), even when χ = 1, there is still a symmetry and the positive and negative branches have the same |ω crit |. We might expect that breaking the up-down symmetry of the system, via the introduction of a vertical density stratification, would cause a break in symmetry of the eigenvalue spectrum, and hence result in different frequencies for the positive and negative branches. Instead, when k = 0, the eigenvalues remain in complex conjugate pairs. We see that in figure 2 (a) and (c), the introduction of a vertical stratification across the layer has, as expected, broken the symmetry of the eigenvalue spectrum -they no longer appear in complex conjugate pairs. However, counter-intuitively, when k = 0 (subfigure(b)), the symmetry is not broken and the eigenvalues remain in complex conjugate pairs, in an analogous way to the Boussinesq case (χ = 1). Evonuk 44 and Glatzmaier et al. 45 describe a mechanism that is perhaps responsible for this difference between NS and EW rolls. The crux of their argument is that the vorticity equation (curl of equation (3)) contains a term proportional to Ω(∇ · u), which is in general, non-zero for anelastic convection. However, in our system, the x-component of this term is zero and so it does not have an effect on EW rolls, whereas, the y-component of this term

IV. NONLINEAR RESULTS
This section extends the work of section III to the nonlinear regime, allowing us to examine the mean flows driven by the system. For these nonlinear results we restrict ourselves to the 2d system which lies in the plane of the rotation vector, the y-z plane, i.e., ∂ ∂x ≡ 0. We solve the nonlinear equations (3)-(7) using a streamfunction, ψ(y, z), defined by and so ∇ · (ρ ref u) = 0 is automatically satisfied. We then write the equations in terms of ψ and and solve them for ω, u, s and ψ using a Fourier-Chebyshev pseudospectral method with a second order, semi-implicit, Crank-Nicolson/Adams-Bashforth time-stepping scheme; for details see Currie 24 and references within. v and w are then straightforward to obtain from ψ.
For diagnostic purposes, we decompose ω, u, s and ψ into means (horizontal averages) and fluctuations; where the mean (denoted by an overbar) is defined as, for example, where L is the length of the computational domain.

A. Bifurcation structure and large-scale solutions
In the limit θ → 0 (χ → 1), our anelastic system reduces to a Boussinesq system. As is typically seen in the Boussinesq case, if Ra is slowly increased from its critical value (with other parameters fixed) then the solutions in the anelastic system undergo a series of bifurcations. An example is shown in figure 4 which shows time series of the Nusselt number (Nu) at different Ra.
We define Nu to be the ratio of the total heat flux to the conductive heat flux in the basic state.
Clearly the system undergoes a number of bifurcations via steady, oscillatory, quasiperiodic and chaotic solutions, en route to chaos (Ruelle-Takens-Newhouse route to chaos 46,47 ). We note that some hysteresis of the solutions was observed, depending on the initial conditions used; however, a full investigation of the bistability is not examined here.
Whilst it is typical for the system to undergo this series of bifurcations, Currie 24 reported on a

B. Generation of mean flows
We investigate how the strength and direction of mean flows driven in our plane layer system are affected by stratification. In all simulations presented we fix the length of our computational box to be L = 5. To begin, we fix P r = 1, Ra = 2 × 10 5 , T a = 10 5 , φ = π 4 and considerū andv for three different stratifications (see figure 6). In (a), χ = 1.5 (close to Boussinesq) in (b), χ = 5 and in (c), χ = 10. Even though the critical Rayleigh number changes with χ and Ra is fixed in these cases, the degree of supercriticality is not vastly changed -ranging from approximately 5.8 to 6.2 times the onset value. A striking feature present in all of the flows shown in figure   6 is that of strong oscillations. These oscillations are likely to be inertial oscillations that arise here because the Rossby number is of order one and so rotation and convection are of roughly equal importance 50 . Such oscillations were also observed in the fully compressible calculations of Brummell et al. 26 . In figure 6 (a), where χ is close to one, i.e., almost Boussinesq, we see that the flows are almost symmetric but, when the strength of stratification is increased, the extent of the asymmetry of the mean flows is also increased. For example, the positive flow ofv in the upper half-plane only just penetrates down into the lower half-plane for χ close to one, but for stronger stratification it penetrates further into the layer.ū is more time-dependent and harder to interpret thanv, but the asymmetry is still evident. The asymmetry results from the fact that, with stress free boundaries, the horizontal mass flux must be zero. From figure 6, other effects of increasing the stratification appear to be that the maximum velocity achieved by the flow decreases as χ increases, but the flows become more systematic. By systematic we mean a flow with a definite mean in the sense that the velocity fluctuations (including oscillations) produce a significant mean when averaged in both horizontal space and time.
Our results share some common features with those found in previous 3d studies. For example, we see strong systematic shear flows driven when the rotation vector is oblique to gravity as were seen in the Boussinesq calculations of Hathaway and Somerville 6 and the fully compressible calculations of Brummell et al. 26 . Furthermore, the asymmetries introduced in the layer when stratification is added are also present in the flows of Brummell et al. 26 . We also find thatū andv are comparable in size; a similar feature is also found in Hathaway and Somerville 6 and Brummell et al. 26 . This is likely to be a result of the horizontal periodic boundary conditions artificially enhancing the meridional flow. We also note that there are some differences between the flows driven here and those in the previous work of Hathaway and Somerville 6 and Brummell et al. 26 . The most obvious difference being that the sense ofū is reversed (whilst the sense of v coincides). We also tend to find that the percentage of energy in the mean flows compared to the total energy is larger in our cases. We expect these differences result from the difference in parameter regimes considered and also the 2d nature of our system.
To quantify the flow properties described above, we consider the mean and variation of the flows in time, and see how they vary with χ and z. In figure 7, we plot the time-averaged mean forū andv along with error bars corresponding to the standard deviation (Σ) from that mean. The first comment we make is that the error bars are significant and the departure of the maximum flow speed (given by the colorbars in figure 6) from the average is also significant; this is because of the oscillations present in the flows (as discussed above). Further, in (a) we see that Σ(ū) is smallest near to mid-layer and grows as we move out towards the boundaries but in (b), Σ(ū) is smallest at a deeper layer. This behaviour is also seen in Σ(v), where for χ = 1.5, the standard deviation is fairly even across the layer but with its smallest value at approximately mid-layer; whereas for χ = 10, the smallest standard deviation is found at much smaller z. Note also, the mean ofū and v is close to zero at z = 0.5 in (a), but there is a significant flow at z = 0.5 in (b). These measures characterise the behaviour we saw in the time-dependent plots in figure 6. As a percentage of its mean, Σ(ū) is larger than Σ(v), indicative of the more time-dependent behaviour ofū we also observed in figure 6 . Figure 8 shows how the standard deviation varies with z for a number of different χ. In general, we see that for stronger stratification Σ is reduced. This behaviour is particularly evident in the lower depths of the layer (smaller z). It is also evident that for Σ(v), the minimum of the standard deviation occurs at a deeper level in the layer as χ is increased. For Σ(ū), the trend is not so clear, however, the flows corresponding to larger χ have a minimum at a lower z than the flows corresponding to smaller χ. Therefore, there are fewer fluctuations at lower levels with increasing χ, and it is this that results in the relatively large time-averaged mean at this level.
Reynolds stresses are known to drive mean flows 6,26 . To analyse their role in mean flow generation, we consider the mean equations obtained by horizontally averaging the x and y components where we have averaged in time and assumed a statistically steady state so that ∂ ∂t ū = ∂ ∂t v = 0. The quantities ρ ref uw, ρ ref vw are the Reynolds stresses terms, they measure the correlation between the horizontal and vertical velocity components. With a tilted rotation vector we might expect these correlations to be nonzero 6 . We note, from equations (19) and (20), that it is the z-derivative of ρ ref vw that drivesū and the z-derivative of ρ ref uw that drivesv. In what follows, for both equations (19) and (20), we refer to the term on the LHS as the Coriolis term, the first term on the RHS as the viscous term and the second term on the RHS as the Reynolds stress (RS) term.
The factor of ρ ref in the Coriolis terms of equations (19) and (20) means that, in theory, for two different χ, if the driving terms on the right-hand side are of the same size, then the case with the largest χ will yield the largestū andv, i.e., at any fixed z, if P rρ refū is the same for two different ρ ref (fixed P r) thenū will be larger for the smaller ρ ref . To see this, we plot each of the terms of equations (19) and (20) in figure 9; in addition, we plotū andv (without the P rρ ref factors). For χ = 1, the Coriolis term is equivalent to the mean flow since P r = 1, therefore no additional line is visible in this case. However, for χ = 1, there is a difference between the Coriolis term and the mean flow itself. In both (a) and (b), the strong dominance of the RS terms is clear. It is also evident that the viscous term is more important in determiningū than it isv, this is because the viscous term affecting each mean flow component depends on the gradient of the other mean flow component and we find that the gradient ofv tends to be larger than that ofū. It is clear that the RS terms are larger in the χ = 1 case and this results in the Coriolis terms being larger for χ = 1.
This effect is most prominent at the top of the layer, where the fluid mass is at its lowest.
In the Boussinesq case, the RS terms are symmetric about the mid-layer depth. As a result of this symmetry, and because the RS terms are the dominant terms in equations (19) and (20), the mean flows are also symmetric about the mid-layer depth. Moreover, if χ is increased, then the RS terms become asymmetric, leading to asymmetric mean flows. We comment that the symmetry of the RS terms can change both as a result of the small-scale turbulent interactions but also through the mean flow acting back on the turbulence.

V. CONCLUSION
The results in this paper can be categorised into two distinct sections. The first section was concerned with linear anelastic convection. There we demonstrated and proved the existence of a previously unknown, hidden symmetry in the equations present when EW rolls are considered.
We explained this by showing that the symmetry breaking term present in the vorticity evolution equation (introduced by the inclusion of a stratification) vanishes under certain circumstances.
The second part of the paper was concerned with the effect of stratification on nonlinear anelastic convection. We showed that efficient large-scale convection cells that have been shown to exist in 2d, Boussinesq Rayleigh-Bénard convection can also be found when stratification is introduced; although, unlike their Boussinesq counterparts, the anelastic solutions may not be time independent.
We went on to examine the effect of stratification on the generation of zonal and meridional mean flows. The use of an idealised 2d tilted plane layer model has allowed us to show the importance of correlations between velocity components resulting in Reynolds stresses that generate that is to be expected at P r = 1). However, we would expect the viscous forces to be much less important in a realistic setting as the diffusivities are much smaller. In fact, Currie 24 found small P r to play an important part in the dynamics of mean flow generation in Boussinesq convection.
In particular, for the same size Reynolds stresses a larger mean flow is driven at smaller P r. We would expect a similar effect to occur here as, from equations (19) and (20) Owing to the existence of magnetic fields in physical systems such as stars and planets and their interaction with convection and rotation, an obvious question to ask is how such magnetic fields modify the mean flows generated. This is a question we address in a subsequent paper. We conclude by acknowledging the limitations of our simple model. Clearly in two dimensions correlations may be over-exaggerated leading to the formation of strong mean flows. We are therefore currently investigating how extending the model to three dimensions weakens correlations and affects the turbulent driving of mean flows.

ACKNOWLEDGMENTS
L.K.C. is grateful to Science and Technology Facilities Council (STFC) for a PhD studentship.
We thank the anonymous referees for comments that helped to focus the content of this paper.

Appendix: Proof
The following is a proof of the symmetry of the spectrum of eigenvalues that exists when k = 0 (see section III B). The proof not only holds for the stress free boundary conditions considered above but is a more general result and holds for all natural boundary conditions. The proof is similar in nature to that of Proctor et al. 51 who prove a similar result. However, they consider a system with symmetric equations but break the symmetry through asymmetric boundary conditions. This is in contrast to this work, where we have asymmetric equations to begin with, and typically our boundary conditions are symmetric.
To begin the proof, we make a change of variables. Let Next, we define the inner product  So, if (σ,X 1 ) is an eigenvalue, eigenfunction pair for the system then so is (σ * ,Ỹ 1 ).
Hence, we have shown that as long as the boundary conditions onX andX * are the same, when k = 0, the eigenvalue spectrum is symmetric. This is in agreement with the numerical results we found in section III A.
If k = 0, then the final term on the right-hand-side of equation (A.5) is non-zero and must be added to the central entry of the matrices B and B † ; this results in a breakdown of the proof, as the last step (from equation (A.13) to equation (A.14)) can not be carried out. Therefore, when k = 0, the eigenvalue spectrum is not symmetric, again in agreement with the numerical results obtained in section III A.