{'fraction_non_alphanumeric': 0.02128700758502569, 'fraction_numerical': 0.009093874887855803, 'mean_word_length': 4.811137440758293, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The role of density stratification in generating zonal flow structures in a rotating fluid

Local generation of vorticity occurs in rotating density-stratified fluids as fluid parcels move radially, expanding or contracting with respect to the background density stratification. Thermal convection in rotating 2D equatorial simulations demonstrates this mechanism. The convergence of the vorticity into zonal flow structures as a function of radius depends on the shape of the density profile, with the prograde jet forming in the region of the disk where the greatest number of density scale heights occurs. The number of stable jets that form in the fluid increases with decreasing Ekman number and decreases with increasing thermal driving. This local form of vorticity generation via the density stratification is likely to be of great importance in bodies that are quickly rotating, highly turbulent, and have large density changes, such as Jovian planets. However, it is likely to be of lesser importance in the interiors of planets such as the Earth, which have smaller density stratifications and are less turbulent.


Introduction
Rotation plays a key role in fluid motions in both the liquid interiors of planets and planetary atmospheres. Much work has been done in the context of thin atmospheres and the Beta-effect (Williams 1979(Williams , 1985(Williams , 2003Cho & Polvani 1996;Allison 2000). The effect of rotation and its role in producing zonal flow patterns as a function of latitude, such as seen on the giant gas planets, has been studied by several groups in the Boussinesq, constant background density, approximation (Christensen 2001(Christensen , 2002Aurnou & Olson 2001;Aurnou & Heimpel 2004;Heimpel et al. 2005); as well as in the anelastic, stratified background density, approximation (Glatzmaier 2005). While studies in 3D and in the 2D meridional plane (Jones et al 2003;Rotvig and Jones 2006) reproduce jet structures similar to those seen in giant planets, these simulations are run in parameter regimes far from those present in planetary bodies. In particular, due to numerical constraints, these simulations are executed with viscosities and diffusivities that are too high, resulting in more diffusive and more laminar flows than should be present in the planets. While some of these simulations (Glatzmaier 2005;Heimpel et al. 2005) display some turbulent behavior they are still not as fully turbulent as we expect planets, such as Jupiter and Saturn, to be. This means that simulations in this more laminar regime interact with their boundaries more than they would in a more turbulent regime. Interactions with the boundaries lead to vorticity generation via the global mechanism of Busse rolls (Proudman 1916;Busse 1983Busse , 1994Busse , 2002, an effect that should become weaker as the fluid becomes more turbulent. Extrapolations have been made to project trends seen in simulations to higher Rayleigh numbers and lower Ekman numbers and thus into the parameter regime proper for planetary bodies (Christensen & Aubert 2006), however these extrapolations are made across many orders of magnitude and therefore may not hold. Further, the bulk of these studies have been conducted with constant background density while the planets exist in the density-stratified regime. As direct 3D simulations of fully turbulent bodies are not possible this paper seeks to expand the explored parameter regime in 2D and to look in particular at the effect of density stratification on the fluid behavior. These 2D simulations are in the equatorial plane without additional terms to add the effect of curvature of the boundaries in the northern and southern hemispheres. This essentially eliminates the effect of these boundaries and therefore any contribution to the vorticity from Busse roll structures. These simulations also span the full equatorial plane and do not include a core, eliminating the effect of a solid core on the fluid behavior. Therefore the vorticity and any jets formed will be purely the consequence of the local interaction of the rotation with the density stratification.
To review how jets form with latitude via interactions with the boundaries in 3D, I will first present vorticity generation via the geostrophic approximation. Next the effects of density will be discussed to show how the density stratification can provide a local source of vorticity. This local mechanism, in the context of the 2D equatorial plane, will be contrasted to the Boussinesq, constant density case. A suite of density stratified numerical simulations will be presented exploring a range of rotation rates (Ekman numbers) and driving rates (Rayleigh numbers) to explore the effect of rotation on the number of jets that are created, as a function of radius in the fluid. While the parameter regimes expected for planetary interiors cannot be reached even in these 2D simulations we will see that the results point towards the growing importance of this local mechanism of vorticity generation, especially in the regime appropriate for giant planets.

The Geostrophic Approximation
Of foremost importance in this discussion is the momentum equation, written here in terms of the anelastic approximation with t being time, ρ the density, u the velocity vector, P the pressure, δ ij the Kronecker delta, ν the viscous diffusion, e ij the viscous stress tensor, g the gravity,r the unit vector in the radial direction, and Ω the rotation rate. Quantities with an over bar indicate a background radial profile. If we assume that the pressure gradient and the Coriolis terms are the dominant terms in the momentum equation, and that the background density (ρ o ) is constant, then we get the following relation often referred to as the Geostrophic Approximation: ∇P = 2ρ o u × Ω. Taking the curl of this equation gives the Taylor-Proudman theorem ∂u/∂z = 0, which dictates to first order that the change in the velocity parallel to the axis of rotation is zero. Therefore vertical columns spanning the convecting fluid are expected to form. These types of columns have been observed in laboratory experiments (Carrigan & Busse 1983;Busse 1994) and have been implied to exist in the Earths core from magnetic observations (Gubbins & Bloxham 1987). How vorticity is generated in this situation can be viewed in the following way: as the column moves toward or away from the axis of rotation its height is forced to change due to the curvature of the bounds of the convective region. The columns change their width to conserve mass, as a column moves towards the axis of rotation it becomes taller and thinner, while as it moves away from the axis it becomes shorter and wider. This, coupled with the Coriolis force, generates positive vorticity in the case of the lengthening column and negative vorticity in the case of the shortening column. If we now add the sense of curvature of the boundary it can be shown (Busse 2002) that in the case of concave curvature, which applies to the outer boundaries of the planets, these columns will become tilted perpendicular to the axis of rotation, allowing the convergence of positive vorticity, prograde flow, on the equator-side of the columns, and negative vorticity, retrograde flow, on the side of the columns towards the axis. Pairs of these Busse rolls can then set up alternating prograde and retrograde flows at the surface as a function of latitude with the primary prograde flow at the equator similar to the flows observed in the surface clouds of the giant planets (Busse 1994). This flow pattern is seen in Boussinesq simulations (Heimpel et al. 2005) and arguably in anelastic simulations as well (Glatzmaier 2005). However, as the level of turbulence in the fluid rises it becomes increasingly unlikely that columns spanning the convective zone can be maintained. The Geostrophic Approximation will no longer be valid as the fluid becomes more turbulent and the inertial terms play a more important role.
However, there is another mechanism that acts locally to generate vorticity that may gain increasing importance as this global mechanism breaks down in the high turbulence regime.

The effect of density stratification
All convecting planetary interiors are density stratified to some extent. The density of the Earths outer metallic core decreases on the order of 23% from bottom to top (Boehler 1996). This translates to 0.2 density scale heights where the number of density scale heights is defined as N ρ = ln(ρ b /ρ t ), with ρ b and ρ t being the densities at the bottom and top of the convective zone respectively. Other planets experience much greater changes in density through their fluid interiors. The density in Jupiters hydrogen and helium envelope changes by over 200 times if the upper level of the convective zone is taken to be at a pressure of 1000 bars, or 10 8 Pa (Guillot 1999). This translates to at least 5.3 density scale heights; naturally, as more of the upper atmosphere is included a greater number of density scale heights will occur. This change in the density does not go unnoticed by the fluid. Simple 2D simulations in a box show that a significant asymmetry in the fluid flow occurs in the density-stratified case, while in the constant density case symmetric flow is observed (Evonuk & Glatzmaier 2004). As fluid plumes in the density-stratified case rise they move into less dense regions and expand, forming mushroom shaped plumes, while the downwelling plumes are compressed. Meanwhile, for constant density cases upwelling and downwelling plumes are in a sense indistinguishable. An additional asymmetry in the density-stratified case appears in the velocities. In denser regions the fluid moves more slowly than in less dense regions. When these asymmetries are coupled to the rotation some interesting effects become apparent (the following is a review of material to be covered in more detail in Glatzmaier et al. 2008, in preparation). In particular expansion and contraction velocities of rising and sinking fluid parcels act via the Coriolis force to produce vorticity. Rising expanding parcels generate negative vorticity while sinking contracting parcels generate positive vorticity. This mechanism is an entirely local process so there is no longer a need for coherent columns spanning the fluid. Instead, many small disconnected fluid columns can form aligned along the axis of rotation and their individual radial movements will generate vorticity. We can examine how this mechanism works in more detail in the limit of two dimensions in the equatorial plane. To look at the time rate of change of the vorticity due to rotation we take the curl of the Coriolis force, ∇ × (2u × Ω) = −2Ω(∇ · u). In the Boussinesq case where the density is constant there is no contribution from the Coriolis term to the vorticity as ∇ · u = 0. However, in the anelastic case, which includes a density-stratification, we see that a non-zero term is produced by this cross product: This term depends on the rate of rotation, the radial velocity, and the local inverse density scale height, h ρ = d lnρ/dr. For a stably stratified fluid where the density decreases with radius, the inverse density scale height is negative everywhere but changes in amplitude depending on the shape of the density curve. In a fluid where the density drops off more quickly towards the surface of the planet (Figure 1a) peak values of the inverse density scale height are seen near the surface of the planet. Conversely, if the density drops off more slowly near the surface than at depth, then peak values of the inverse density scale height occur at depth. This plays an important role in the convergence of the angular momentum. A simple schematic (Figure 1b, c) shows that fluid parcels leaving both the inner and outer regions in the equatorial plane both ultimately transfer negative vorticity away from the boundaries as the sinking contracting parcel and the rising expanding parcel both veer in the direction retrograde to that of rotation. If we consider a density profile where the peak negative values of the inverse density scale height occur at the upper boundary, or near the surface, then the transport of negative vorticity from the upper boundary will be dominant, resulting in a prograde flow near the surface of the disk and a retrograde flow at depth. Conversely, if the peak inverse density scale height occurs at depth then a retrograde flow will form near the surface.
While Busse rolls are unlikely to exist in the turbulent interiors of giant planets, in the more laminar anelastic simulations in 3D we expect to see both the global process of Busse rolls and the local process via the density stratification to play a role in the formation of zonal flow structures. Since these processes are difficult to separate in 3D it is of interest to look at the effect of the density change in 2D to eliminate Busse rolls and to approach more turbulent regimes.

Numerics
Thermal-convection simulations are performed using the Evonuk-Glatzmaier finite-volume code on a Cartesian grid in 2D (Evonuk & Glatzmaier 2006). The basic equations solved are the momentum equation (Eq. 1), conservation of mass (Eq. 2), and the energy equation in terms of the entropy S as shown below: (5) In this equation κ T is the thermal diffusivity, C P is the specific heat capacity at constant pressure, κ R is the radiative diffusivity, T is the temperature, and Q s is the heating function. In these simulations the rotation rate, Ω, and the rate of heating, Q s , were varied to explore a range of parameters. The use of the finite volume method allows the inner core to be easily removed from the simulations to focus on the effect of rotation and the density stratification as there are no numerical issues associated with the origin point. Over four density scale heights were included along the profile shown in Figure 1a with numerical values appropriate for the interior of a Jovian planet with no core (Guillot 1999) but with a shape that is likely to also be appropriate for an Earth-like planet without a core. In the context of the above discussion this will imply a convergence of positive vorticity or prograde flow near the surface. To prevent the convecting fluid from splashing against the outer boundary of the simulation or feeling the jaggedness of the outer boundary, a stable buffer zone was added in the outer 10% of the disk. Simulation results will be viewed looking down from the north so that eastward, or prograde motion, is in the counterclockwise direction. Simulations were performed with grid resolutions from 500 × 500 to 1600 × 1600.
There are three primary dimensionless numbers which can be used to characterize the equations. The Rayleigh number, provides a measure of the buoyancy terms versus the diffusive terms. Here it is written in terms of the change in entropy, ∆S, of the fluid in the convecting region. As there is no inner core, D is the radius of the disk and the total change in entropy is calculated from the output of the simulation. Higher Rayleigh numbers are generally indicative of higher driving and more turbulent flows. The Ekman number, provides a comparison of the viscous terms to the Coriolis term. Smaller Ekman numbers indicate fluids where the rotation is more dominant in determining fluid behavior. The Prandtl number, P r = ν/κ T , is the a ratio of the viscous diffusivity to the thermal diffusivity. In all of the simulations the Prandtl number was taken to be 1.0 and the diffusivities were prescribed to be constant. A fourth dimensionless number, the Reynolds number, Re = U D/ν, compares the inertial terms to the viscous terms, providing a non-dimensional way of looking at the time-averaged mean Euclidean norm of the velocity, U, of the fluid.

Results
The fluid flow patterns for three representative cases are shown in Figure 2. The middle panel shows a case with two jets, as hypothesized by the above discussion with an exterior prograde jet and an interior retrograde jet. A case with a smaller rate of rotation (also characteristic of a larger rate of driving) shows the flow moving directly through the interior forming a dipolar flow structure (Figure 2c). This flow pattern occurs when the effect of rotation is too small, with respect to the driving, to form jets. Conversely, at higher rotation rates or lower driving, a third jet can form (Figure 2a). The amplitudes of the jets in these representative cases are shown in Figure 3 which shows the time and longitudinally-averaged angular velocity (i.e. the tangential velocity divided by the radius) scaled to the rotation rate for the cases with three jets and two jets in Figure 2. These families, grouped by number of jets can be seen in Figure 4a showing the number of jets as a function of Ekman and Reynolds numbers. In agreement with Figure 2, as the rotation rate increases (Ekman number decreases) more jets are formed, whereas when the fluid is driven harder and the velocities increase (Reynolds number increases) fewer jets are formed.
It is also of interest to look the directionality of the velocities, that is, the ratio of the time-averaged maximum radial (u r ) to maximum tangential (u φ ) velocities. This ratio is shown in Figure 4b plotted against the Ekman number. As the rotation rate increases, the vertical convection is suppressed and a greater fraction of the fluid flow is aligned in the tangential direction, as would be expected when more jets are formed. As large-scale vertical convection is suppressed, smaller scale vortices develop to transport energy radially, an effect shown in Evonuk & Glatzmaier (2006). This graph shows a very clear division between the number of jets with decreasing Ekman numbers and the ratio of maximum radial velocity to tangential velocity.
The convective heat flux, F c = C Pρ T u r is a function of the temperature perturbation and the radial velocity, while the diffusive heat flux, F d = −C Pρ κ T ∂(T + T )/∂r depends primarily on the total temperature,T +T . The relationship between the Ekman number and the ratio of the mean convective heat flux to the mean diffusive heat flux is therefore very similar in distribution (though not in value) to that of the relationship between the Ekman number and the Reynolds number (Figure 4a), with larger ratios of convective heat flux to diffusive heat flux corresponding to higher Reynolds numbers. The Rayleigh number versus the ratio of the mean heat fluxes shows a very strong correlation between the two quantities as higher Rayleigh numbers mean larger driving and therefore larger radial velocities ( Figure  4c). In this plot there appears to be no meaningful distinction based purely on number of jets present in the simulations.
Lastly, Figure 4d shows the Ekman number versus the Rayleigh number. Again, as the Ekman number decreases, more jets develop, while as the Rayleigh number increases (buoyancy increases due to increased driving), fewer jets form. Since planets are in a parameter regime of very low Ekman number and very high Rayleigh number it is of interest to see where they would lie on this plot. Figure 5 expands this graph to include representative values for the Earths liquid outer core and Jupiter for comparison. It is obvious that even 2D simulations are far from the parameter regimes expected for these planets. Any trends chosen need to be extrapolated over many orders of magnitude. Two possible extrapolations have been included in Figure 6 to demonstrate this sensitivity. A conservative extrapolation places both the Earth and Jupiter in a regime of three jets, or a region where local vorticity generation is very important. However a trend could also be chosen to put them both in the transition area between the three-jet and two-jet regions, where this local mechanism becomes less important.
Also of importance when interpreting these results for individual planets is the extent of the density stratification. These simulations have been done with four density scale heights. Therefore we expect a graph appropriate for the Earth, with 0.2 density scale heights, to have the Earths location on the graph shifted to the right with respect to the observed pattern. This would shift the Earth into the region of parameter space where the densitystratification plays a smaller role in generating vorticity. Conversely, Jupiter has a greater number of density scale heights; this will shift Jupiter to the left with respect to this pattern placing Jupiter more firmly in the parameter space where local vorticity generation, and the inclusion of the planetary density stratification, is very important.

Conclusions
This paper presents a local form of vorticity generation via interactions of the fluid moving through the background density profile with the Coriolis term. Simulations in 2D show this mechanism organizing the fluid into zonal flows with respect to radius in the equatorial plane with more jets forming as the rotation rate is increased. As the rotation rate is increased the radial convection of the fluid is suppressed and the energy moves radially via small convective cells imbedded in the zonal flow, as highlighted in Evonuk & Glatzmaier (2006). High driving suppresses the effect of rotation and results in a lower number of jets. This form of local vorticity generation is likely to be less important in planets that have small density stratifications and less turbulent flow, such as in the Earths core. However in giant planets where there are a large number of density scale heights and the flow is very turbulent this form of vorticity generation is likely to play an important role. Therefore, in addition to asymmetries in the fluid flow structures and velocities, local vorticity generation provides another argument for including the density-stratification when modeling the fluid behaviors and magnetic field generation particularly in giant planets.
Lastly, it is important to keep in mind that in 3D, as well as in 2D, we are still far from the parameter regimes of planetary interiors and should therefore be cautious when extrapolating current modeling results to planetary interiors.
I would like to thank Jon Aurnou for several helpful discussions on the topic as well as Chris Finlay and Gary Glatzmaier. Thanks also to Philippe Cardin for inviting me to give a talk at IUGG and allowing me to attend the Les Houches Dynamo workshop, both of which helped me organize my thoughts on the subject. Computations were run on the UCSC Beowulf cluster, upsand, and on the CSCS cluster, palu, through ETH. Fig. 1.-Plot of a background density profile (solid line) and its corresponding inverse density scale heights, h ρ (dashed line), as a function of radius for a planet whose density drops off more quickly near the surface of its convective zone (a). Rising, expanding material generates negative vorticity due to Coriolis forces while sinking, contracting material generates positive vorticity (b). In (c) positive vorticity is advected towards both the inner and outer regions of the disk, however the dominant process occurs at the outer boundary due to that region's peak negative value of the inverse density scale height. Grey arrows indicate positive vorticity.   Lines have been added to show approximate divisions between the different groups. Plot of Rayleigh number versus the ratio of the mean convective heat flux to the mean diffusive heat flux (c). There is no apparent dependance on the number of jets present. Plot of Ekman number versus Rayleigh number (d). Increasing rotation rate leads to more jets while increasing Rayleigh number decreases the number of jets.  Figure 4d, but expanded to show the parameter space relevant to planetary interiors. Proposed parameter values for the Earths liquid outer core and Jupiter are plotted for comparison. Simulations are still far, even in 2D, from the regimes in which these planets reside. Possible extrapolations for the transition between the two and three jet regions are shown in light grey. Simulations were done with four density scale heights, so the patterns seen should be shifted for the Earth (0.2 density scale heights) and Jupiter (over five density scale heights). Dark grey arrows show that the Earth would plot further to the right in a region where the local vorticity generation is less important and Jupiter would plot further to the left in a region where local vorticity generation is more important.