Turbulence in large scale two-dimensional balanced hard sphere gas flow

In recent works we developed a model of balanced gas flow where the momentum equation possesses an additional mean field forcing term, which originates from the hard sphere interaction potential between the gas particles. We demonstrated that, in our model, a turbulent gas flow with a Kolmogorov kinetic energy spectrum develops from an otherwise laminar initial jet. In the current work, we investigate the possibility of a similar turbulent flow developing in a large scale two-dimensional setting, where a strong external acceleration compresses the gas into a relatively thin slab along the third dimension. The main motivation behind the current work is the following. According to observations, horizontal turbulent motions in the Earth atmosphere manifest in a wide range of spatial scales, from hundreds of meters to thousands of kilometers. Yet, the air density rapidly decays with altitude, roughly by an order of magnitude each 15-20 kilometers. This naturally raises the question as to whether or not there exists a dynamical mechanism which can produce large scale turbulence within a purely two-dimensional gas flow. To our surprise, we discover that our model indeed produces turbulent flows and the corresponding Kolmogorov energy spectra in such a two-dimensional setting.


Introduction
In his famous work, Reynolds [27] demonstrated that an initially laminar flow of a liquid consistently develops turbulent motions whenever the high Reynolds number condition is satisfied. Later, Kolmogorov [16,17,18] and Obukhov [23,24,25] observed that the time-averaged Fourier spectra of the kinetic energy of an atmospheric flow possess a universal decay structure, corresponding to the inverse five-thirds power of the Fourier wavenumber. With help of detailed measurements from the Global Atmospheric Sampling Program, Nastrom and Gage [22] estimated the energy spectra of meridional and zonal atmospheric winds, and found that the spectra exhibited the inverse cubic power scaling at large scales, and the inverse five-thirds power scaling at moderate and small scales. As of current, the physics of turbulence formation in a laminar flow, as well as the origin of power scaling of turbulent kinetic energy spectra, remain unknown.
In our recent work [3], we considered a system of many particles, each of mass m, interacting solely via a repelling short-range potential φ(r), with the convention that φ(r) → 0 as r → ∞. In the limit of infinitely many such particles, we obtained, via the standard Bogoliubov-Born-Green-Kirkwood-Yvon (BBGKY) formalism [5,6,15], the following Vlasov-type equation [30] for the mass-weighted distribution density of a single particle f (t, x, v): Above, the mean field potentialφ is given via The quantities ρ, p and θ are the mass density, pressure and kinetic temperature of f , respectively, given, together with the mean velocity u, via Further, computing the usual hierarchy of the transport equations for the velocity moments of f , and closing it under the infinite Reynolds number assumption, in [3] we arrived at the following equations for a balanced compressible gas flow: (1.4) ∂ρ ∂t + ∇ · (ρu) = 0, ∂(ρu) ∂t + ∇ · (ρu 2 ) + ∇p = −∇φ, ∂p ∂t + u · ∇p = 0.
The difference between the equations above and the conventional compressible Euler equations [4,10] is that the former preserve the pressure p along the stream lines (balanced flow), and have the mean field forcing −∇φ in the momentum equation.
In [3], we numerically simulated the equations in (1.4) for a gas of hard spheres with the mass and diameter corresponding to those of argon, flowing through a straight threedimensional pipe. For the initial condition of that simulation, we selected a straight laminar Bernoulli jet, which is a steady state in the conventional Euler equations. We observed, however, that in our model (1.4) such a jet quickly develops into a fully turbulent flow. We also examined the Fourier spectrum of the kinetic energy of the simulated flow, and found that its time average decayed with the rate of inverse five-third power of the wavenumber, which corresponded to the famous Kolmogorov spectrum.
The main goal of the current work is to develop and test a similar framework for a gas flow under a strong gravity acceleration, which, on a large scale, is effectively twodimensional, and, therefore, it may be impractical to introduce the vertical direction from the computational perspective. While it is obvious that such a two-dimensional flow, by its simplified nature, lacks many features of a full three-dimensional flow (even with the latter confined to a relatively thin slab), our goal here is to examine whether or not turbulent features could manifest in such a simplified two-dimensional model, which is severely restricted in other physical aspects by its own low dimensionality. This is important, in particular, for long-range climate prediction models, where the restricted dimension may be necessitated by the requirement of a long-range simulation.
The work is organized as follows. In Section 2 we derive the Vlasov equation which incorporates an external acceleration in a potential form. In Section 3 we assume that the strong constant external acceleration is applied in the vertical direction, and further arrive at the system of equations for a balanced flow, confined to a horizontal plane. In Section 4 we show the results of numerical simulations of a large scale flow in two scenarios -one is an inertial jet, and another is a cyclostrophic vortex. We observe that turbulent motions with Kolmogorov spectra manifest in both scenarios from laminar initial conditions. Section 5 summarizes the results of the work.

The Vlasov-type equation with external acceleration
Here we repeat the derivation of the Vlasov-type equation in (1.1) for the motion of particles under an external acceleration. In the presence of an external acceleration, the equations of motion for N particles interacting via a potential φ(r) are given via where g(x) is the external acceleration which acts on any particle at the location x, and depends only on the coordinates of that particle. Further, we will assume that g(x) is a potential acceleration, that is, there exists a potential function h(x) such that g(x) = −∂h(x)/∂x. As in our previous works [2,3], we concatenate Additionally, due to the presence of the external acceleration, here we concatenate In the new notations, the Liouville equation for the density of states F(t, X, V ) is given via Due to the presence of the external acceleration, the total momentum of the system is no longer a first integral. However, due to the fact that g(x) is a potential acceleration, the total energy E of the system remains an invariant of motion: Therefore, similarly to [2,3], any suitable function of the form is automatically a steady state of the Liouville equation in (2.4). Among those functions, the canonical Gibbs equilibrium state is given via where θ 0 is the equilibrium kinetic temperature of the system. For the purpose of the work, below we will assume that the gas is dilute (that is, the mean free path of a particle is much larger than the range of the potential interaction), and thus the contribution of Φ(X) to the integral for Z N is negligible: The latter leads to the following form of the canonical Gibbs state for a dilute gas:

Preservation of the Rényi divergences.
Despite the presence of the external acceleration, it can be easily shown that the Liouville equation in (2.4) still preserves the family of general Rényi divergences [26], as long as the integration by parts is not affected by the boundary conditions. Indeed, for two differentiable functions ψ 1 and ψ 2 , we have The Kullback-Leibler divergence [19] is a special case of the Rényi divergence with α = 1.
As we noted in [2,3], the preservation of the Rényi divergences plays a similar role to Boltzmann's H-theorem -that is, if the initial state of F is close to F 0 in the sense of any Rényi metric, then the solution will also remain a nearby state.

Joint two-particle marginal distribution of the Gibbs state.
In what follows, it is important to note the structure of the joint two-particle marginal distribution F G of the Gibbs state for a dilute gas in (2.9). Integrating F G in (2.9) over all particles but the first two, and taking advantage of the dilute gas assumption, one easily obtains (2.12) is the Gibbs state for a single particle. Here, note that, even at a thermodynamic equilibrium, the particles in a pair cannot be regarded as being independent, since there is a spatial correlation term through the interaction potential φ.
2.3. Transport equation for a single particle. We denote the marginal distribution for the first particle via f , that is, where, for convenience, we drop the subscripts from x 1 and v 1 . We then integrate the Liouville equation in (2.4) over dx 2 dv 2 . . . dx N dv N : where F 1,i is the joint distribution of the first and i-th particles, that is, . . dx N dv N , and y, w serve as dummy variables of integration for the coordinate and velocity, respectively, of the i-th particle. Observe that there cannot be any boundary effects in the velocity integration, because, realistically, F and its derivatives must vanish at infinite velocities. However, there could still be boundary effects in the spatial integration, since the spatial domain is generally bounded, and the presence of the external acceleration g introduces spatial anisotropy. Thus, we retain the terms with spatial derivatives of the joint two-particle distribution in the equation for now.

2.4.
A closure for the single particle transport equation. To obtain the transport equation for f alone, we need to express F (2) 1,i via f . For that, a standard assumption is that all pairs of particles are identically distributed [1,3], so that the transport equation becomes Following our recent work [3], we assume that F (2) has the same form as its corresponding equilibrium Gibbs distribution, that is .
Here, the kinetic temperature θ can no longer be regarded as being at equilibrium, and, instead, we take it at the midpoint between the two particles, for symmetry. Next, we renormalize f so that it becomes the mass density, that is, f → Nm f . This leads, as N → ∞, to where we integrated the product w f (w) over dw and obtained ρu, according to (1.3). Via the divergence theorem, the volume integral in the right-hand side can be expressed as the surface integral over the domain boundary S, where [ρu] in is the net inward momentum flux through the boundary. Above we note that x is located inside the domain, while y S is on the boundary, which allows us to set the exponential to 1 as long as the range of potential is sufficiently short. This leads to In a practical scenario, we can assume that the net momentum flux through the boundary is zero; indeed, if we have a horizontal slab domain with strong gravity and impenetrable bottom, the momentum fluxes through the top and bottom are zero (the former due to exponentially vanishing density, the latter due to impenetrability), and we can additionally assume that the net momentum flux through the sides is also zero if the total amount of the gas particles inside our domain remains fixed. This further leads to the following Vlasov-type equation for f with an external acceleration:

The two-dimensional moment equations in the presence of gravity
Here, we are going to assume that the vector of external acceleration g is constant, and that the reference frame is chosen so that g points in the opposite direction of the z-axis (that is, g = (0, 0, −g), and, as a consequence, h(x) = gz). In this situation, it is convenient to separate the variables into the horizontal and vertical ones; from now on, x = (x, y) and v = (v x , v y ) will denote the horizontal coordinate and velocity vectors, respectively, while z and w will refer to the vertical coordinate and velocity, respectively. Accordingly, ∇ = (∂/∂x, ∂/∂y) will refer to the horizontal spatial differentiation. In the new notations, the gravity-forced Vlasov-type equation in (2.21) becomes We now assume that the horizontal dimensions of the domain are much larger than the characteristic scale of the exponential decay of f in the vertical dimension. As a result, it can be assumed that, on the horizontal scale, the dynamics are largely two-dimensional, with the vertical dimension supplying certain parameterized effects. To describe this, we now integrate f over dz: which results in the advection term of the form as f vanishes exponentially with altitude. This leads to the following Vlasov-type equation forf :

The velocity moment equations.
For the zero-order moment we integrate the Vlasovtype equation forf in (3.4) over dvdw, with the notations This results in the following equation for the densityρ: where the latter identity is in effect since the bottom of the domain is impenetrable and thus u z | z=0 = 0. Observe that there is no advective coupling to u z , and only the horizontal momentum is present in the advection term of (3.6). Thus, it suffices to include only the horizontal momentum at the next level of the moment hierarchy. To obtain the transport equation for the horizontal momentum, we integrate (3.4) over vdvdw: where the horizontal pressure tensorP is given via The bottom boundary term in (3.7) can be expressed via where we use the fact that u z | z=0 = 0. Here, S xz | z=0 denotes the vector of the surface shear stress: To obtain the transport equation for the horizontal pressure tensorP, we integrate (3.4) over v 2 dvdw, which leads to the energy equation: In order to isolate the equation for the horizontal pressure tensorP alone, we observe that, from the horizontal momentum equation (3.7), Substituting the expression for the time derivative and noting that we arrive at the horizontal pressure tensor equation of the form The bottom boundary term in the above equation can be expressed as where Q x 2 z | z=0 is the surface skewness matrix, given via At this point, we define the pressurep in a standard way, that is, as a normalized trace ofP, and denote the corresponding horizontal shear stress deviator viaS: Next, we obtain the equation forp by taking the normalized trace of the equation forP in (3.14): Above,q and q z are the horizontal and surface heat fluxes, respectively: The equation for the horizontal shear stressS is obtained by subtracting the pressure equation from (3.14): Following [3], we split the horizontal skewness tensor into the appropriate combination of the horizontal heat fluxes and the traceless deviatorQ: In the same fashion as in our recent work [3], here we assume that both horizontal and surface shear stresses are negligible in comparison with the rest of the terms, which indicates a high Reynolds number flow. For the horizontal shear stressS to remain small, we, first, impose the balance of external forces and boundary effects in its equation (3.20), that is, where the latter equation is the trace of the former. Second, we also impose the balance between the pressure-velocity terms and the heat fluxes in the equation (3.20) for the horizontal shear stress, where, again, the latter equation is the trace of the former. Substituting the above relations into the pressure equation, we arrive at that is, in the same manner as in [3], the pressure is preserved along the stream lines (balanced flow). For the vanishing shear stress, the momentum equation in (3.7) becomes The density equation in (3.6), together with the pressure equation (3.24) and the momentum equation (3.25), constitute the two-dimensional balanced flow equations in the presence of the gravitational acceleration. For a numerical computation, it is more convenient to reformulate the pressure equation (3.24) in the form of a conservation law. For that, we can use the inverse kinetic temperature as a variable; indeed, denotingθ =p/ρ, one can, with the help of (3.6), rearrange (3.24) in the form Together, the equations (3.6), (3.25) and (3.26) form the system of nonlinear conservation laws for the densityρ, momentumρũ and inverse kinetic temperatureθ −1 , where the latter can be related to the standard temperature (in • K) via an appropriate gas constant.

Approximation for the vertically integrated mean field potential.
To estimate the effect of the potential forcing in the momentum equation (3.25), we assume that the gas particles interact via the hard sphere mean field potential -that is φ(r) is zero when r is greater than the diameter of the protection sphere [8] (that is, the sphere of minimal distance between the centers of the colliding hard spheres), and becomes infinite for r less than that. In such a case, the spatial integral in (1.2) equals the volume of this protection sphere regardless of the value of the kinetic temperature: Above, V prot and V HS denote the volumes of the protection and hard spheres, respectively, with the latter being eight times smaller than the former, because the size of the protection sphere is twice that of any of the two colliding hard spheres. Substituting V HS into (1.2) and denoting the density of the hard sphere via ρ HS = m/V HS , we arrive at a rather simple expression for the hard sphere mean field potentialφ HS : Note that the hard sphere mean field potential in (3.28) introduces the same correction into the momentum equation in (1.4) as does the Enskog correction in the Enskog-Euler equations (see Abramov [1] for details). Next, we assume that, in the vertical direction, the hydrostatic balance is achieved: We now express ρ from the hydrostatic balance relation and write For the integral of ρp in the vertical direction we thus have (3.31) Above, the variable η = η(z) is given via and we assume that θ remains bounded above and below with altitude (which is indeed the case in the troposphere and stratosphere). For the vertical integral of the hard sphere mean field potentialφ HS we thus have What remains to be done is to evaluate the surface pressure p| z=0 . For that, we recall that the surface pressure equals the total mass of the air column above the surface per unit area, multiplied by the gravity acceleration, which yields and, subsequently, As a result, the horizontal momentum equation for the hard sphere gas becomes

Numerical simulations
Here we present the results of numerical simulations of a balanced two-dimensional large scale gas flow for two scenarios -an inertial jet, and a cyclostrophic vortex. In both scenarios, we recreate (rather crudely) the main parameters of Earth's atmosphere. In particular, the acceleration constant in (3.36) is set to g = 9.81 m/s 2 , which corresponds to Earth's gravity, while the vertically integrated densityρ and temperatureθ of the gas are chosen to correspond to those typically observed in large scale atmospheric flows. 4.1. Density of the air-like hard sphere gas. The Earth's atmosphere consists of air, which is a mixture of various gases -primarily diatomic, such as the molecular nitrogen N 2 and molecular oxygen O 2 , but with a small amount of polyatomic gases, such as water H 2 O and carbon dioxide CO 2 . Our theory is, however, developed for a monatomic hard sphere gas. Thus, in order to numerically simulate the behavior of the Earth atmosphere using our model, we need to "create" a hypothetical hard sphere gas with a density which matches known physical properties of air. Here, we choose the two properties of our hard sphere gas to match with the air -the molar mass, and the viscosity.
As follows from our gas model above, the equations in (3.6), (3.24), (3.26) and (3.36) are derived from the assumption that the molecules interact in a fully time-reversible fashion. On the contrary, in the standard kinetic theory, the viscosity properties of a hard sphere gas are obtained directly from the time-irreversible Boltzmann collision integral [7,9,11,13]. This is not necessarily a contradiction -clearly, in a real gas, both types of interactions are present (for example, while potential interactions are timereversible, the quantum-mechanical Pauli repulsion is inherently stochastic), and can apparently manifest at different spatial scales. In particular, turbulence is not observed in flows with high Knudsen number (such as thin channels and capillaries) where the gas motion is dominated by viscous effects, whereas the situation reverses itself at low Knudsen numbers. Thus, the same gas can be treated in the context of the time-reversible model at large scales, and as a time-irreversible process at viscous scales.
Treating the collisions as time-irreversible at viscous scales, we recall the well known formula for the dynamic viscosity of the hard sphere gas [9,11,13]: where σ is the diameter of the hard sphere, N A = 6.02214 · 10 23 mol −1 is the Avogadro number, and M is the molar mass of the gas. In order to relate the expression above to the density of the hard sphere ρ HS , we observe that πN A σ 3 , where R = 8.31446 kg m 2 /K mol s 2 is the universal gas constant, and T is the usual temperature in • K. Expressing σ via ρ HS from the latter equation, and substituting the expression for θ, we obtain the formula for ρ HS via the viscosity µ and temperature T: . With (4.3), we can "create" a theoretical hard sphere analog of any gas with a prescribed molar mass M and known viscosity µ at a given temperature T. For air, M = 2.897 · 10 −2 kg/mol, and we take µ = 1.8194 · 10 −5 kg/m s at T = 293.15 K [14], for which (4.3) yields ρ HS = 1850 kg/m 3 . This value is used throughout all computations below.

Two simulated scenarios of a balanced flow.
Below we study the following two special cases of a balanced flow: • Inertial jet. In an inertial jet, the pressure is constant throughout the domain. Observe that, in our equations, the Coriolis acceleration of Earth rotation is not present, which corresponds roughly to the equatorial region. Thus, such inertial flow describes a special case of geostrophic flow near the equator. • Cyclostrophic vortex. In a cyclostrophic vortex, the centripetal force, acting on a rapidly rotating gas, is balanced by the pressure gradient, and the velocity is orthogonal to both. Again, given the absence of the Coriolis acceleration, this scenario corresponds to a large scale rapidly rotating flow, where the Coriolis acceleration is negligible in comparison with the centripetal acceleration -that is, a fully developed tropical cyclone.

Computation of initial and boundary conditions.
Below we describe how the initial and boundary values of the velocity, temperature and pressure are specified in each simulated scenario.
• First, we specify the velocityũ of the flow in the domain, and on those boundaries, where the Dirichlet boundary condition is specified. The way we specify the velocity is entirely scenario-dependent; in the inertial jet scenario, the initial velocity field forms a jet stream, while in the cyclostrophic vortex scenario the initial velocity field forms a rotating flow with an appropriate radial profile. • Onceũ is defined, we compute the kinetic temperatureθ of the flow using the Bernoulli law for a compressible gas: where the background temperature of the gas is set toT 0 = 250 K, that is, the approximate vertically averaged temperature of air throughout the troposphere.
Observe, however, that in present setting the gas is two-dimensional, and, therefore, the effective adiabatic exponent γ = 2. • Onceũ andθ are defined, we specify the pressurep (which automatically yields the density viaρ =p/θ). Below, we study two scenarios: an inertial jet, and a cyclostrophic vortex. In the inertial jet scenario, we set the pressure to a constant valuep 0 , whereρ 0 = 10 4 kg/m 2 , which constitutes the approximate mass of the air column per unit area of the Earth surface. In the cyclostrophic vortex scenario, the pressure gradient is balanced by the centripetal force: where r is the coordinate of the point relative to the center of rotation, Ω is the angular velocity of rotation, and the latter identity is valid for a planar rotation (that is, when r ⊥ Ω). Here, the velocityũ is orthogonal to both r and Ω, that is, (4.7)ũ = Ω × r, and, since r ⊥ Ω, ũ = Ω r .
Assuming, in turn, thatρ,ũ andp do not depend on the angle, and only depend on r , we denote r = r ,ũ = ũ , and arrive at the following explicit formula forp: (4.9) dp dr =ρũ Above, in the integration formula we assume thatũ vanishes sufficiently far away from the center of rotation. Also, in practice, one can simplify the integration by presuming thatθ varies weakly in comparison toũ 2 under the integral, and factor the kinetic temperature out of the integration (as we do further below).

Numerical methods and software.
In the current work, we use the same software as in our previous work [3] -namely, we use OpenFOAM [32] to perform all numerical simulations. Noting that the equations (3.6), (3.26) and (3.36), for the density, inverse temperature and momentum, respectively, comprise a system of nonlinear conservation laws, we simulate them with the help of an appropriately modified rhoCentralFoam solver [12], which uses the central scheme of Kurganov and Tadmor [20] for the numerical finite volume discretization, with the flux limiter due to van Leer [29]. The time-stepping of the method is adaptive, based on the 20% of the maximal Courant number.

Numerical simulation of an inertial jet.
In the inertial jet scenario, we simulate a large scale jet stream in a channel-like domain. The spatial domain is 2500 km in length, and 500 km in width. The spatial discretization is uniform in both directions, with a step of 2.5 km, which constitutes 1000 cells in the zonal direction, and 200 cells in the meridional direction. The domain has a 100 km-wide inlet in the middle of the western wall, while the outlet is the whole 500 km-wide eastern boundary. The initial velocity field is given via the shear jet where the vertical coordinate y is given relative to the central axis of the channel, d = 50 km, and u 0 = 35 m/s. As we can see, the speed of the flow is 35 m/s in the middle of the channel, and smoothly decays to zero at the distance of 50 km both to the north and south (which also corresponds to the boundaries of the inlet). The weak asymmetry parameter α = 0.05 is introduced in the same manner as by McCalpin and Haidvogel [21], to break up possible mirror symmetry effects of the flow in the two-dimensional domain.
Observe that the typical kinematic viscosity of air at normal conditions is ∼ 10 −5 m 2 /s, the characteristic size of our flow is no less than 10 5 m (the width of the jet), and the reference velocity is no less than 10 m/s, which yields the value of the Reynolds number Re ∼ 10 11 . Therefore, the assumptions made in the course of the derivation of our model are clearly valid for this set-up.
The boundary conditions are specified as follows. At the inlet, the velocity is set to (4.10) at all times, with the rest of the variables following the procedure set forth in Section 4.3. At the outlet, all variables are set to the free outflow conditions (that is, zero normal gradient). At the walls, which comprise the remainder of the domain boundary, the velocity is set to the no-slip condition, while the rest of the variables is set to the zero normal gradient condition (which corresponds to zero momentum flux for the pressure, and zero heat flux for the temperature).
We emphasize that these initial and boundary conditions correspond to a steady state for the conventional Euler equations, and even for our balanced flow equations taken without the mean field forcing in the momentum equation (3.25) (or its special case (3.36) for a hard sphere gas). Thus, any observed non-steady effects originate from the presence of the mean field forcing due to the molecular interaction potential.
Starting from the initial conditions described above, we integrate the system of equations in (3.6), (3.26)  from the CFL criterion, varied between 12 and 16 seconds. In Figure 1, we show the speed ũ of the flow, captured at 12, 18, 24 and 30 hours elapsed model time, on contour plots. Here we can see that small fluctuations in the jet manifest at 12 hours. At 18 hours, the jet visibly meanders roughly between 1000 and 2000 km. At 24 and 30 hours, the jet is completely broken up in the second half of the domain. In fact, visually, the snapshots of the flow speed at 24 and 30 hours resemble the famous Reynolds experiment [28]. After 30 hours, the general configuration of the flow does not seem to visibly change any further, and thus we do not show the snapshots past that time.
The vorticityω is given via the curl of velocityũ: (4.11)ω = ∇ ×ũ. Sinceũ is confined to the xy-plane,ω is orthogonal to this plane, and only has the z-component, given via In Figure 2, we show the contour plots ofω z for the inertial flow, also captured at 12, 18, 24 and 30 hours elapsed model time. The situation here is similar to what was observed in Figure 1 for the speed of the inertial flow; namely, small vorticity fluctuations are observed at 12 hours, and visible meandering manifests at 18 hours. At 24 hours, the structure of vorticity resembles a Kármán vortex trail [31]; at 30 hours the structure of the vorticity is completely lost in the second half of the channel, and the motion appears to be fully chaotic. In Figure 3, we show the time-averaged kinetic energy spectrum of the flow, computed in the rectangular sub-region of the domain, which extends from 1500 to 2500 km zonally, and from −100 to +100 km meridionally. This spectrum was computed as in [3]: first, the kinetic energy E(x) = u(x) 2 /2 was averaged across the channel (thus becoming the function of the x-coordinate only). Then, the linear trend was subtracted from the result in the same manner as was done by Nastrom and Gage [22], to ensure that there was no sharp discontinuity between the energy values at the western and eastern boundaries of the region. Finally, the one-dimensional discrete Fourier transformation was applied to the result. The subsequent time-averaging of the modulus of the Fourier transform was computed between t = 24 and t = 72 hours. As we can see in Fig Despite promising results in simulating turbulent motions of the velocity, our model also has its limitations. In particular, for a fully developed turbulent flow, we found that the densityρ may attain unrealistic values. In Figure 4, we show the contour plots of the density for the same times as the velocity and vorticity above, namely, at 12, 18, 24 and 30 hours. Observe that, in a fully developed chaotic flow, the density varies between 4 · 10 3 and 1.6 · 10 4 kg/m 2 (with the background value 10 4 kg/m 2 ), which, of course, does not happen in Earth's atmosphere.
The likely reason for this behavior is that our model is derived under the assumption that the flow remains balanced unconditionally, while in reality balanced flows manifest largely in the presence of relatively small density gradients. Should a large density gradient develop locally, the flow is likely to become isentropic in that region, and the usual compressible Euler equations would apply instead. Thus, in order to correct such a behavior of our model, one likely needs to introduce an appropriate mechanism of controlled switching to the compressible Euler equations locally, should the density gradient become "too large". In the present work, however, we report the results of our simulation without any corrections.

Numerical simulation of a cyclostrophic vortex.
In the scenario with a cyclostrophic vortex, we simulate a large scale rapidly rotating flow, resembling a tropical cyclone, in a square domain of size 1000×1000 km. The spatial discretization step is set to 2 km, which constitutes 500 cells in both horizontal directions. The boundary of the domain is impenetrable with appropriate boundary conditions, and the counter-clockwise vortex is confined entirely within it. Naturally, the direction of the velocity at each point is orthogonal to the direction towards the center of the vortex, and the speed of the flow is a function of the distance to the center, but not the angle. As a function of the distance to the center r, we set the speed of the flow to (4.13)ũ(r) =    u 0 4 ln(r/r 0 ) ln(R 0 /r) Here, the maximum speed u 0 = 40 m/s is achieved at the geometric mean distance r max u = √ r 0 R 0 , where r 0 and R 0 are the inner and outer radii of the vortex, respectively. For the simulation, we set R 0 = 10r 0 = 500 km (and thus, r max u ≈ 158 km). The reason whyũ(r) is chosen in this manner, is that it is, effectively, a function of ln r, which synergizes well with the computation of the integral in the cyclostrophic pressure profile formula (4.9). The Reynolds number of this scenario is similar to that of the inertial flow examined above. Upon integrating this initial set-up forward in time, we observed the following behavior. The time step, chosen from the CFL criterion, varied between 5 and 6 seconds. The flow remained largely laminar for the initial 150 minutes of the elapsed model time, after which turbulent motions rapidly developed around the region of the maximal flow speed. Shortly after 200 minutes of the elapsed model time, a numerical instability occurred in the solution, and a floating point exception was generated by the software. Thus, turbulent flow was observed between 150 and 200 minutes of the elapsed time.
A possible reason for the manifestation of the numerical instability is the apparent lack of any damping effects in the model, similarly to what we observed in our recent work [2]. In the inertial jet scenario above, the developing regions with large density gradients left the domain through the outlet before causing further problems (in a way, damping was created by the boundary conditions); here, however, the rotating flow is contained entirely within the domain, thus allowing instabilities to develop without restrictions.  In Figure 5, we show the speed ũ of the flow, captured at 160, 170, 180 and 190 minutes of the elapsed model time, on contour plots. Here we can see that small speed fluctuations manifest at 160 minutes, and become progressively larger at 170, 180 and 190 minutes. These fluctuations are, however, mainly confined to the region of maximal velocity of the vortex, with the flow at the "outskirts" remaining laminar.
In Figure 6, we show the vertical vorticity componentω z of the flow, captured at 160, 170, 180 and 190 minutes of the elapsed model time, on contour plots. Likewise, one can observe that the vorticity is largely confined to the region with maximal flow speed, and its fluctuations become progressively larger with elapsed time.
In Figure 7, we show the time-averaged kinetic energy spectrum of the flow. This spectrum was computed in the square sub-region of the domain, which extends from −400 to +400 km both meridionally and zonally, to reduce the effects of relatively "calm" corners of the domain. This computation was done in the same manner as above for the inertial flow, except, due to isotropy of the flow, here we computed the spectrum both in the zonal and meridional directions, and then took the average. The time averaging was done over the interval between 180 and 200 minutes of the elapsed model time, during which the most turbulent flow was observed. As we can see in Figure 7, the kinetic energy spectrum largely matches Kolmogorov's k −5/3 -decay throughout the whole range of the wavenumbers. In Figures 8 and 9, we show the contour plots of the pressurep and densityρ, respectively, of the cyclostrophic vortex, captured at 160, 170, 180 and 190 minutes of the elapsed model time. Note that while the pressure assumes the form of a "basin" with closed level curves (which is what naturally occurs in cyclones), the density exhibits the same unrealistic variations as above in Figure 4 for the inertial flow.

Summary
In the current work, we investigate the ability of our model of a balanced compressible hard sphere gas flow [3] to produce turbulent motions from a laminar initial condition in a purely two-dimensional setting, which corresponds to the large scale dynamics of the Earth atmosphere. The equations for such a flow are derived in the same manner as in our previous work, with the additional condition that a strong external acceleration in the downward direction, together with an impenetrable bottom boundary of the domain, compress the gas into a relatively thin horizontal slab.
We simulate two prototype scenarios of large scale atmospheric dynamics -an inertial jet, and a cyclostrophic vortex. In the inertial jet scenario, we observe that an initially laminar straight jet flow develops turbulent motions, bearing a close resemblance to Reynolds' famous experiment [28]. Moreover, the Fourier spectrum of the kinetic energy of the turbulent part of the inertial flow shows the Kolmogorov k −5/3 -decay at moderate scales, and the k −8/3 -decay at small scales. In the cyclostrophic vortex scenario, turbulent motions develop in the region of the maximal flow speed, whereas the Kolmogorov k −5/3 -decay of the kinetic energy spectrum is observed throughout the whole range of the Fourier wavenumbers.
The main result of the current work is the surprising discovery that, in our model of a balanced compressible hard sphere gas flow, turbulent motions with Kolmogorov kinetic energy spectra develop not only in a fully three-dimensional flow [3], but also in a restricted, two-dimensional setting. If our balanced gas flow equations indeed constitute a faithful model of the actual phenomenon of turbulence in gases, this result suggests that it should be possible to capture large scale turbulent atmospheric features even in simplified, lower-dimensional settings, such as those which are typically used in long range climate change predictions.