Turbulence in Large-Scale Two-Dimensional Balanced Hard Sphere Gas Flow

Abstract

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. However, the air density rapidly decays with altitude, roughly by an order of magnitude each 15–20 km. 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.

1. Introduction

In his famous work, Reynolds [1] demonstrated that an initially laminar flow of a liquid consistently develops turbulent motions whenever the high Reynolds number condition is satisfied. Later, Kolmogorov [2,3,4] and Obukhov [5,6,7] 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 the help of detailed measurements from the Global Atmospheric Sampling Program, Nastrom and Gage [8] 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 yet, 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 [9], we considered a system of many particles, each of mass m, interacting solely via a repelling short-range potential $\varphi \left(r\right)$, with the convention that $\varphi \left(r\right)\to 0$ as $r\to \infty$. In the limit of infinitely many such particles, we obtained, via the standard Bogoliubov–Born–Green–Kirkwood–Yvon (BBGKY) formalism [10,11,12], the following Vlasov-type equation [13] for the mass-weighted distribution density of a single particle $f(t,\mathit{x},\mathit{v})$:

$$\frac{\partial f}{\partial t}+\mathit{v}·\frac{\partial f}{\partial \mathit{x}}=\frac{1}{\rho }\frac{\partial \overline{\varphi }}{\partial \mathit{x}}·\frac{\partial f}{\partial \mathit{v}}.$$

(1)

Above, the mean field potential $\overline{\varphi }$ is given via

$$\overline{\varphi }(t,\mathit{x})=\frac{1}{2m}\rho (t,\mathit{x})p(t,\mathit{x})\int \left[1-exp\left(-\frac{\varphi (\parallel \mathit{y}\parallel )}{\theta (t,\mathit{x})}\right)\right]\mathrm{d}\mathit{y}.$$

(2)

The quantities $\rho$, p and $\theta$ are the mass density, pressure and kinetic temperature of f, respectively, given, together with the mean velocity $\mathit{u}$, via

$$\rho =\int f\mathrm{d}\mathit{v},\rho \mathit{u}=\int \mathit{v}f\mathrm{d}\mathit{v},p=\rho \theta =\frac{1}{3}\int {\parallel \mathit{v}-\mathit{u}\parallel }^{2}f\mathrm{d}\mathit{v}.$$

(3)

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 [9], we arrived at the following equations for a balanced compressible gas flow:

$$\frac{\partial \rho }{\partial t}+\nabla ·\left(\rho \mathit{u}\right)=0,\frac{\partial \left(\rho \mathit{u}\right)}{\partial t}+\nabla ·\left(\rho {\mathit{u}}^2\right)+\nabla p=-\nabla \overline{\varphi },\frac{\partial p}{\partial t}+\mathit{u}·\nabla p=0.$$

(4)

The difference between the equations above and the conventional compressible Euler equations [14,15] is that the former preserve the pressure p along the stream lines (balanced flow), and have the mean field forcing $-\nabla \overline{\varphi }$ in the momentum equation.

In [9], we numerically simulated the equations in (4) for a gas of hard spheres with the mass and diameter corresponding to those of argon, flowing through a straight three-dimensional pipe. For the initial condition of this 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 (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 two-dimensional, 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.

2. The Vlasov-Type Equation with External Acceleration

Here, we repeat the derivation of the Vlasov-type equation in (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 $\varphi \left(r\right)$ are given via

$$\frac{\mathrm{d}{\mathit{x}}_i}{\mathrm{d}t}={\mathit{v}}_i,\frac{\mathrm{d}{\mathit{v}}_i}{\mathrm{d}t}=\mathit{g}\left({\mathit{x}}_i\right)-\frac{\partial }{\partial {\mathit{x}}_i}\sum _{\genfrac{}{}{0.0pt}{}{j=1}{j\ne i}}^N\varphi (\parallel {\mathit{x}}_i-{\mathit{x}}_j\parallel ),$$

(5)

where $\mathit{g}\left(\mathit{x}\right)$ is the external acceleration that acts on any particle at the location $\mathit{x}$, and depends only on the coordinates of that particle. Further, we will assume that $\mathit{g}\left(\mathit{x}\right)$ is a potential acceleration; that is, there exists a potential function $h\left(\mathit{x}\right)$ such that $\mathit{g}\left(\mathit{x}\right)=-\partial h\left(\mathit{x}\right)/\partial \mathit{x}$. As in our previous works [9,16], we concatenate

$$\mathit{X}=({\mathit{x}}_1,\dots ,{\mathit{x}}_N),\mathit{V}=({\mathit{v}}_1,\dots ,{\mathit{v}}_N),\Phi \left(\mathit{X}\right)=\sum _{i=1}^{N-1}\sum _{j=i+1}^N\varphi (\parallel {\mathit{x}}_i-{\mathit{x}}_j\parallel ).$$

(6)

Additionally, due to the presence of the external acceleration, here, we concatenate

$$H\left(\mathit{X}\right)=\sum _{i=1}^{N}h\left({\mathit{x}}_i\right).$$

(7)

In the new notations, the Liouville equation for the density of states $F(t,\mathit{X},\mathit{V})$ is given via

$$\frac{\partial F}{\partial t}+\mathit{V}·\frac{\partial F}{\partial \mathit{X}}=\frac{\partial (H+\Phi )}{\partial \mathit{X}}·\frac{\partial F}{\partial \mathit{V}}.$$

(8)

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 $\mathit{g}\left(\mathit{x}\right)$ is a potential acceleration, the total energy E of the system remains an invariant of motion:

$$E=\frac{{\parallel \mathit{V}\parallel }^2}{2}+H\left(\mathit{X}\right)+\Phi \left(\mathit{X}\right)=\mathrm{const}.$$

(9)

Therefore, similarly to [9,16], any suitable function of the form

$$F_0(\mathit{X},\mathit{V})=F_0\left(E\right)=F_0\left({\parallel \mathit{V}\parallel }^2/2+H\left(\mathit{X}\right)+\Phi \left(\mathit{X}\right)\right)$$

(10)

is automatically a steady state of the Liouville equation in (8). Among those functions, the canonical Gibbs equilibrium state is given via

$$F_G=\frac{1}{{\left(2\pi {\theta }_0\right)}^{3N/2}{Z}_N}exp\left(-\frac{{\parallel \mathit{V}\parallel }^2/2+H+\Phi }{{\theta }_0}\right),Z_N=\int e^{-(H+\Phi )/{\theta }_0}\mathrm{d}\mathit{X},$$

(11)

where ${\theta }_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 $\Phi \left(\mathit{X}\right)$ to the integral for $Z_N$ is negligible:

$$Z_N=\int e^{-(H+\Phi )/{\theta }_0}\mathrm{d}\mathit{X}\approx \int e^{-H/{\theta }_0}\mathrm{d}\mathit{X}={\left(\int e^{-h/{\theta }_0}\mathrm{d}\mathit{x}\right)}^N=Z_1^N.$$

(12)

The latter leads to the following form of the canonical Gibbs state for a dilute gas:

$$F_G(\mathit{X},\mathit{V})=e^{-\Phi \left(\mathit{X}\right)/{\theta }_0}\prod _{i=1}^N{f}_G({\mathit{x}}_i,{\mathit{v}}_i),f_G(\mathit{x},\mathit{v})=\frac{e^{-{(\parallel \mathit{v}\parallel }^2/2+h\left(\mathit{x}\right))/{\theta }_0}}{{\left(2\pi {\theta }_0\right)}^{3/2}{Z}_1}.$$

(13)

2.1. Preservation of the Rényi Divergences

Despite the presence of the external acceleration, it can be easily shown that the Liouville equation in (8) still preserves the family of general Rényi divergences [17], as long as the integration by parts is not affected by the boundary conditions. Indeed, for two differentiable functions ${\psi }_1$ and ${\psi }_2$, we have

$$\begin{array}{c}\frac{\partial }{\partial t}\int {\psi }_1\left(F\right){\psi }_2\left(F_0\right)\mathrm{d}\mathit{X}\mathrm{d}\mathit{V}=\int \frac{\partial {\psi }_1\left(F\right)}{\partial t}{\psi }_2\left(F_0\right)\mathrm{d}\mathit{X}\mathrm{d}\mathit{V}=\int {\psi }_2\left(F_0\right)(\frac{\partial (H+\Phi )}{\partial \mathit{X}}·\frac{\partial }{\partial \mathit{V}}\hfill \\ \hfill -\mathit{V}·\frac{\partial }{\partial \mathit{X}}){\psi }_1\left(F\right)\mathrm{d}\mathit{X}\mathrm{d}\mathit{V}=\int {\psi }_1\left(F\right)\left(\mathit{V}·\frac{\partial }{\partial \mathit{X}}-\frac{\partial (H+\Phi )}{\partial \mathit{X}}·\frac{\partial }{\partial \mathit{V}}\right){\psi }_2\left(F_0\right)\mathrm{d}\mathit{X}\mathrm{d}\mathit{V}=0.\end{array}$$

(14)

Choosing ${\psi }_1\left(x\right)=x^{\alpha }$, ${\psi }_2\left(x\right)=x^{1-\alpha }$ for some parameter $\alpha >0$ leads to the preservation of the Rényi divergence $D_{\alpha }(F,F_0)$:

$$D_{\alpha }(F,F_0)=\frac{1}{\alpha -1}ln\int F^{\alpha }{F}_0^{1-\alpha }\mathrm{d}\mathit{X}\mathrm{d}\mathit{V}.$$

(15)

The Kullback–Leibler divergence [18] is a special case of the Rényi divergence with $\alpha =1$. As we noted in [9,16], 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.

2.2. 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^{\left(2\right)}$ of the Gibbs state for a dilute gas in (13). Integrating $F_G$ in (13) over all particles but the first two, and taking advantage of the dilute gas assumption, one easily obtains

$$F_G^{\left(2\right)}({\mathit{x}}_1,{\mathit{v}}_1,{\mathit{x}}_2,{\mathit{v}}_2)=\int F_G\mathrm{d}{\mathit{x}}_3\mathrm{d}{\mathit{v}}_3\dots \mathrm{d}{\mathit{x}}_N\mathrm{d}{\mathit{v}}_N=f_G({\mathit{x}}_1,{\mathit{v}}_1)f_G({\mathit{x}}_2,{\mathit{v}}_2)e^{-\varphi (\parallel {\mathit{x}}_2-{\mathit{x}}_1\parallel )/{\theta }_0},$$

(16)

where $f_G(\mathit{x},\mathit{v})$ 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 $\varphi$.

2.3. Transport Equation for a Single Particle

We denote the marginal distribution for the first particle via f, that is,

$$f(t,\mathit{x},\mathit{v})=\int F\mathrm{d}{\mathit{x}}_2\mathrm{d}{\mathit{v}}_2\dots \mathrm{d}{\mathit{x}}_N\mathrm{d}{\mathit{v}}_N,$$

(17)

where, for convenience, we drop the subscripts from ${\mathit{x}}_1$ and ${\mathit{v}}_1$. We then integrate the Liouville equation in (8) over $\mathrm{d}{\mathit{x}}_2\mathrm{d}{\mathit{v}}_2\dots \mathrm{d}{\mathit{x}}_N\mathrm{d}{\mathit{v}}_N$:

$$\frac{\partial f}{\partial t}+\mathit{v}·\frac{\partial f}{\partial \mathit{x}}+\mathit{g}·\frac{\partial f}{\partial \mathit{v}}=\sum _{i=2}^N\int \left(\frac{\partial }{\partial \mathit{x}}\varphi (\parallel \mathit{x}-\mathit{y}\parallel )·\frac{\partial }{\partial \mathit{v}}-\mathit{w}·\frac{\partial }{\partial \mathit{y}}\right)F_{1,i}^{\left(2\right)}(\mathit{x},\mathit{v},\mathit{y},\mathit{w})\mathrm{d}\mathit{y}\mathrm{d}\mathit{w},$$

(18)

where $F_{1,i}^{\left(2\right)}$ is the joint distribution of the first and i-th particles, that is,

$$F_{1,i}^{\left(2\right)}(\mathit{x},\mathit{v},{\mathit{x}}_i,{\mathit{v}}_i)=\int F\mathrm{d}{\mathit{x}}_2\mathrm{d}{\mathit{v}}_2\dots \mathrm{d}{\mathit{x}}_{i-1}\mathrm{d}{\mathit{v}}_{i-1}\mathrm{d}{\mathit{x}}_{i+1}\mathrm{d}{\mathit{v}}_{i+1}\dots \mathrm{d}{\mathit{x}}_N\mathrm{d}{\mathit{v}}_N,$$

(19)

and $\mathit{y}$, $\mathit{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 $\mathit{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_{1,i}^{\left(2\right)}$ via f. For this, a standard assumption is that all pairs of particles are identically distributed [9,19], so that the transport equation becomes   

$$\frac{\partial f}{\partial t}+\mathit{v}·\frac{\partial f}{\partial \mathit{x}}+\mathit{g}·\frac{\partial f}{\partial \mathit{v}}=(N-1)\int \left(\frac{\partial }{\partial \mathit{x}}\varphi (\parallel \mathit{x}-\mathit{y}\parallel )·\frac{\partial }{\partial \mathit{v}}-\mathit{w}·\frac{\partial }{\partial \mathit{y}}\right)F^{\left(2\right)}(\mathit{x},\mathit{v},\mathit{y},\mathit{w})\mathrm{d}\mathit{y}\mathrm{d}\mathit{w}.$$

(20)

Following our recent work [9], we assume that $F^{\left(2\right)}$ has the same form as its corresponding equilibrium Gibbs distribution, that is

$$F^{\left(2\right)}(\mathit{x},\mathit{v},\mathit{y},\mathit{w})=f(\mathit{x},\mathit{v})f(\mathit{y},\mathit{w})exp\left(-\frac{\varphi (\parallel \mathit{x}-\mathit{y}\parallel )}{\theta \left(\frac{\mathit{x}+\mathit{y}}{2}\right)}\right).$$

(21)

Here, the kinetic temperature $\theta$ 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\to Nmf$. This leads, as $N\to \infty$, to

$$\frac{\partial f}{\partial t}+\mathit{v}·\frac{\partial f}{\partial \mathit{x}}=\left(\frac{1}{\rho }\frac{\partial \overline{\varphi }}{\partial \mathit{x}}-\mathit{g}\right)·\frac{\partial f}{\partial \mathit{v}}-\frac{f}{m}\int \frac{\partial }{\partial \mathit{y}}·\left[\rho \left(\mathit{y}\right)\mathit{u}\left(\mathit{y}\right)exp\left(-\frac{\varphi (\parallel \mathit{x}-\mathit{y}\parallel )}{\theta \left(\frac{\mathit{x}+\mathit{y}}{2}\right)}\right)\right]\mathrm{d}\mathit{y},$$

(22)

where we integrated the product $\mathit{w}f\left(\mathit{w}\right)$ over $\mathrm{d}\mathit{w}$ and obtained $\rho \mathit{u}$, according to (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,

$$\begin{array}{c}\int \frac{\partial }{\partial \mathit{y}}·\left[\rho \left(\mathit{y}\right)\mathit{u}\left(\mathit{y}\right)exp\left(-\frac{\varphi (\parallel \mathit{x}-\mathit{y}\parallel )}{\theta \left(\frac{\mathit{x}+\mathit{y}}{2}\right)}\right)\right]\mathrm{d}\mathit{y}\hfill \\ \hfill {=\int \left(\rho \mathit{u}\right)|}_{S}exp\left(-\frac{\varphi (\parallel \mathit{x}-{\mathit{y}}_S\parallel )}{\theta \left(\frac{\mathit{x}+{\mathit{y}}_S}{2}\right)}\right){·\mathit{n}\mathrm{d}S=\int \left(\rho \mathit{u}\right)|}_S·\mathit{n}\mathrm{d}S=-{\left[\rho \mathit{u}\right]}_{\mathrm{in}},\end{array}$$

(23)

where ${\left[\rho \mathit{u}\right]}_{\mathrm{in}}$ is the net inward momentum flux through the boundary. Above, we note that $\mathit{x}$ is located inside the domain, while ${\mathit{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

$$\frac{\partial f}{\partial t}+\mathit{v}·\frac{\partial f}{\partial \mathit{x}}=\left(\frac{1}{\rho }\frac{\partial \overline{\varphi }}{\partial \mathit{x}}-\mathit{g}\right)·\frac{\partial f}{\partial \mathit{v}}+\frac{f}{m}{\left[\rho \mathit{u}\right]}_{\mathrm{in}}.$$

(24)

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:

$$\frac{\partial f}{\partial t}+\mathit{v}·\frac{\partial f}{\partial \mathit{x}}=\left(\frac{1}{\rho }\frac{\partial \overline{\varphi }}{\partial \mathit{x}}-\mathit{g}\right)·\frac{\partial f}{\partial \mathit{v}}.$$

(25)

3. The Two-Dimensional Moment Equations in the Presence of Gravity

Here, we are going to assume that the vector of external acceleration $\mathit{g}$ is constant, and that the reference frame is chosen so that $\mathit{g}$ points in the opposite direction of the z-axis (that is, $\mathit{g}=(0,0,-g)$, and, as a consequence, $h\left(\mathit{x}\right)=gz$). In this situation, it is convenient to separate the variables into the horizontal and vertical ones; from now on, $\mathit{x}=(x,y)$ and $\mathit{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, $\nabla =(\partial /\partial x,\partial /\partial y)$ will refer to the horizontal spatial differentiation. In the new notations, the gravity-forced Vlasov-type equation in (25) becomes

$$\frac{\partial f}{\partial t}+\mathit{v}·\nabla f+w\frac{\partial f}{\partial z}=\frac{\nabla \overline{\varphi }}{\rho }·\frac{\partial f}{\partial \mathit{v}}+\left(g+\frac{1}{\rho }\frac{\partial \overline{\varphi }}{\partial z}\right)\frac{\partial f}{\partial w}.$$

(26)

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 $\mathrm{d}z$:

$$\tilde{f}(t,\mathit{x},\mathit{v},w)={\int }_0^{\infty }f(t,\mathit{x},z,\mathit{v},w)\mathrm{d}z.$$

(27)

which results in the advection term of the form

$${\int }_0^{\infty }\left(\frac{\partial f}{\partial t}+\mathit{v}·\nabla f+w\frac{\partial f}{\partial z}\right)\mathrm{d}z=\frac{\partial \tilde{f}}{\partial t}+\mathit{v}·\nabla \tilde{f}{-wf|}_{z=0},$$

(28)

as f vanishes exponentially with altitude. This leads to the following Vlasov-type equation for $\tilde{f}$:

$$\frac{\partial \tilde{f}}{\partial t}+\mathit{v}·\nabla \tilde{f}{=wf|}_{z=0}+\frac{\partial }{\partial \mathit{v}}·{\int }_0^{\infty }\frac{\nabla \overline{\varphi }}{\rho }f\mathrm{d}z+\frac{\partial }{\partial w}{\int }_0^{\infty }\left(g+\frac{1}{\rho }\frac{\partial \overline{\varphi }}{\partial z}\right)f\mathrm{d}z.$$

(29)

3.1. The Velocity Moment Equations

For the zero-order moment, we integrate the Vlasov-type equation for $\tilde{f}$ in (29) over $\mathrm{d}\mathit{v}\mathrm{d}w$, with the notations

$$\tilde{\rho }=\int \tilde{f}\mathrm{d}\mathit{v}\mathrm{d}w,\tilde{\rho }\tilde{\mathit{u}}=\int \mathit{v}\tilde{f}\mathrm{d}\mathit{v}\mathrm{d}w.$$

(30)

This results in the following equation for the density $\tilde{\rho }$:

$$\frac{\partial \tilde{\rho }}{\partial t}+\nabla ·\left(\tilde{\rho }\tilde{\mathit{u}}\right)=\left(\rho u_z\right){|}_{z=0}=0,$$

(31)

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 (31). 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 (29) over $\mathit{v}\mathrm{d}\mathit{v}\mathrm{d}w$:

$$\frac{\partial \left(\tilde{\rho }\tilde{\mathit{u}}\right)}{\partial t}+\nabla ·(\tilde{\rho }{\tilde{\mathit{u}}}^2+\tilde{\mathit{P}}){=\int \mathit{v}wf|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w-\nabla {\int }_0^{\infty }\overline{\varphi }\mathrm{d}z,$$

(32)

where the horizontal pressure tensor $\tilde{\mathit{P}}$ is given via

$$\tilde{\mathit{P}}=\int {(\mathit{v}-\tilde{\mathit{u}})}^2\tilde{f}\mathrm{d}\mathit{v}\mathrm{d}w.$$

(33)

The bottom boundary term in (32) can be expressed via

$$\begin{array}{c}{\int \mathit{v}wf|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w=\int (\mathit{v}-\mathit{u}{|}_{z=0}){wf|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w+\mathit{u}{{|}_{z=0}\int wf|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w\hfill \\ \hfill =\int (\mathit{v}-\mathit{u}{|}_{z=0})(w-u_z{|}_{z=0})f{{|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w={\mathit{S}}_{\mathit{x}z}|}_{z=0},\end{array}$$

(34)

where we use the fact that $u_z{|}_{z=0}=0$. Here, ${\mathit{S}}_{\mathit{x}z}{|}_{z=0}$ denotes the vector of the surface shear stress:   

$${\mathit{S}}_{\mathit{x}z}{|}_{z=0}={\left.\int \left(\begin{array}{c}(v_x-u_x)(w-u_z)\\ (v_y-u_y)(w-u_z)\end{array}\right)f\mathrm{d}\mathit{v}\mathrm{d}w\right|}_{z=0}.$$

(35)

To obtain the transport equation for the horizontal pressure tensor $\tilde{\mathit{P}}$, we integrate (29) over ${\mathit{v}}^2\mathrm{d}\mathit{v}\mathrm{d}w$, which leads to the energy equation:

$$\begin{array}{c}\frac{\partial }{\partial t}(\tilde{\rho }{\tilde{\mathit{u}}}^2+\tilde{\mathit{P}})+\nabla ·\left(\tilde{\rho }{\tilde{\mathit{u}}}^3+\tilde{\mathit{u}}\tilde{\mathit{P}}+{\left(\tilde{\mathit{u}}\tilde{\mathit{P}}\right)}^T+{\left(\tilde{\mathit{u}}\tilde{\mathit{P}}\right)}^{TT}\right)+\nabla ·\int {(\mathit{v}-\tilde{\mathit{u}})}^3\tilde{f}\mathrm{d}\mathit{v}\mathrm{d}w\hfill \\ \hfill =\int {\mathit{v}}^2{wf|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w-{\int }_0^{\infty }\left(\nabla \overline{\varphi }{\mathit{u}}^T+\mathit{u}\nabla {\overline{\varphi }}^T\right)\mathrm{d}z.\end{array}$$

(36)

In order to isolate the equation for the horizontal pressure tensor $\tilde{\mathit{P}}$ alone, we observe that, from the horizontal momentum equation (32),

$$\begin{array}{c}\frac{\partial \left(\tilde{\rho }{\tilde{\mathit{u}}}^2\right)}{\partial t}=\frac{\partial \left(\tilde{\rho }\tilde{\mathit{u}}\right)}{\partial t}{\tilde{\mathit{u}}}^T+\tilde{\mathit{u}}\frac{\partial \left(\tilde{\rho }{\tilde{\mathit{u}}}^T\right)}{\partial t}-{\tilde{\mathit{u}}}^2\frac{\partial \tilde{\rho }}{\partial t}\hfill \\ =\left({\int \mathit{v}wf|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w-{\int }_0^{\infty }\nabla \overline{\varphi }\mathrm{d}z\right){\tilde{\mathit{u}}}^T+\tilde{\mathit{u}}{\left({\int \mathit{v}wf|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w-{\int }_0^{\infty }\nabla \overline{\varphi }\mathrm{d}z\right)}^T\\ \hfill +\nabla ·\left(\tilde{\rho }\tilde{\mathit{u}}\right){\tilde{\mathit{u}}}^2-\nabla ·(\tilde{\rho }{\tilde{\mathit{u}}}^2+\tilde{\mathit{P}}){\tilde{\mathit{u}}}^T-\tilde{\mathit{u}}\nabla ·{(\tilde{\rho }{\tilde{\mathit{u}}}^2+\tilde{\mathit{P}})}^T.\end{array}$$

(37)

Substituting the expression for the time derivative and noting that

$$\nabla ·\left(\tilde{\rho }\tilde{\mathit{u}}\right){\tilde{\mathit{u}}}^2-\nabla ·\left(\tilde{\rho }{\tilde{\mathit{u}}}^2\right){\tilde{\mathit{u}}}^T-\tilde{\mathit{u}}\nabla ·{\left(\tilde{\rho }{\tilde{\mathit{u}}}^2\right)}^T=-\nabla ·\left(\tilde{\rho }{\tilde{\mathit{u}}}^3\right),$$

(38)

we arrive at the horizontal pressure tensor equation of the form

$$\begin{array}{c}\frac{\partial \tilde{\mathit{P}}}{\partial t}+\nabla ·\left(\tilde{\mathit{u}}\tilde{\mathit{P}}\right)+\tilde{\mathit{P}}\nabla \tilde{\mathit{u}}+\nabla {\tilde{\mathit{u}}}^T\tilde{\mathit{P}}+\nabla ·\int {(\mathit{v}-\tilde{\mathit{u}})}^3\tilde{f}\mathrm{d}\mathit{v}\mathrm{d}w\hfill \\ \hfill =\int {(\mathit{v}-\tilde{\mathit{u}})}^2{wf|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w-{\int }_0^{\infty }\left(\nabla \overline{\varphi }{(\mathit{u}-\tilde{\mathit{u}})}^T+(\mathit{u}-\tilde{\mathit{u}})\nabla {\overline{\varphi }}^T\right)\mathrm{d}z.\end{array}$$

(39)

The bottom boundary term in the above equation can be expressed as

$$\begin{array}{c}\int {(\mathit{v}-\tilde{\mathit{u}})}^2{wf|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w={\int (\mathit{v}-\mathit{u}|}_{z=0}{+\mathit{u}|}_{z=0}-\tilde{\mathit{u}}{{)}^{2}wf|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w\hfill \\ =\int {(\mathit{v}-\mathit{u}{|}_{z=0})}^{2}wf{{|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w+\int (\mathit{v}-\mathit{u}{|}_{z=0}){(\mathit{u}{|}_{z=0}-\tilde{\mathit{u}})}^{T}wf|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w\\ +\int (\mathit{u}{|}_{z=0}-\tilde{\mathit{u}}){(\mathit{v}-\mathit{u}{|}_{z=0})}^{T}wf{{|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w+\int {(\mathit{u}{|}_{z=0}-\tilde{\mathit{u}})}^{2}wf|}_{z=0}\mathrm{d}\mathit{v}\mathrm{d}w\\ \hfill ={\mathit{Q}}_{{\mathit{x}}^{2}z}{|}_{z=0}+{\mathit{S}}_{\mathit{x}z}{{|}_{z=0}{(\mathit{u}{|}_{z=0}-\tilde{\mathit{u}})}^T+(\mathit{u}{|}_{z=0}-\tilde{\mathit{u}}){\mathit{S}}_{\mathit{x}z}|}_{z=0}^T,\end{array}$$

(40)

where ${\mathit{Q}}_{{\mathit{x}}^{2}z}{|}_{z=0}$ is the surface skewness matrix, given via

$$\begin{array}{c}{\mathit{Q}}_{{\mathit{x}}^{2}z}{|}_{z=0}={\left.\int {(\mathit{v}-\mathit{u})}^{2}wf\mathrm{d}\mathit{v}\mathrm{d}w\right|}_{z=0}\hfill \\ \hfill ={\left.\int \left(\begin{array}{cc}{(v_x-u_x)}^2(w-u_z)& (v_x-u_x)(v_y-u_y)(w-u_z)\\ (v_x-u_x)(v_y-u_y)(w-u_z)& {(v_y-u_y)}^2(w-u_z)\end{array}\right)f\mathrm{d}\mathit{v}\mathrm{d}w\right|}_{z=0}.\end{array}$$

(41)

At this point, we define the pressure $\tilde{p}$ in a standard way—that is, as a normalized trace of $\tilde{\mathit{P}}$—and denote the corresponding horizontal shear stress deviator via $\tilde{\mathit{S}}$:

$$\tilde{p}=\frac{1}{2}\mathrm{tr}\left(\tilde{\mathit{P}}\right),\tilde{\mathit{S}}=\tilde{\mathit{P}}-\tilde{p}\mathit{I}.$$

(42)

Next, we obtain the equation for $\tilde{p}$ by taking the normalized trace of the equation for $\tilde{\mathit{P}}$ in (39):

$$\begin{array}{c}\frac{\partial \tilde{p}}{\partial t}+\nabla ·\left(\tilde{p}\tilde{\mathit{u}}\right)+\tilde{p}\nabla ·\tilde{\mathit{u}}+\tilde{\mathit{S}}:\nabla \tilde{\mathit{u}}+\nabla ·\tilde{\mathit{q}}\hfill \\ \hfill =q_z{{|}_{z=0}+{\mathit{S}}_{\mathit{x}z}|}_{z=0}·(\mathit{u}{|}_{z=0}-\tilde{\mathit{u}})-{\int }_0^{\infty }(\mathit{u}-\tilde{\mathit{u}})·\nabla \overline{\varphi }\mathrm{d}z.\end{array}$$

(43)

Above, $\tilde{\mathit{q}}$ and $q_z$ are the horizontal and surface heat fluxes, respectively:

$$\tilde{\mathit{q}}=\frac{1}{2}\int \parallel \mathit{v}-\tilde{\mathit{u}}{{\parallel }^2(\mathit{v}-\tilde{\mathit{u}})\tilde{f}\mathrm{d}\mathit{v}\mathrm{d}w,q_z|}_{z=0}=\frac{1}{2}{\left.{\int \parallel \mathit{v}-\mathit{u}\parallel }^2(w-u_z)f\mathrm{d}\mathit{v}\mathrm{d}w\right|}_{z=0}.$$

(44)

The equation for the horizontal shear stress $\tilde{\mathit{S}}$ is obtained by subtracting the pressure equation from (39):

$$\begin{array}{c}\frac{\partial \tilde{\mathit{S}}}{\partial t}+\nabla ·\left(\tilde{\mathit{u}}\tilde{\mathit{S}}\right)+\tilde{\mathit{S}}\nabla \tilde{\mathit{u}}+\nabla {\tilde{\mathit{u}}}^T\tilde{\mathit{S}}-(\tilde{\mathit{S}}:\nabla \tilde{\mathit{u}})\mathit{I}+\tilde{p}(\nabla \tilde{\mathit{u}}+\nabla {\tilde{\mathit{u}}}^T-(\nabla ·\tilde{\mathit{u}})\mathit{I})\hfill \\ +\nabla ·\int {(\mathit{v}-\tilde{\mathit{u}})}^3\tilde{f}\mathrm{d}\mathit{v}\mathrm{d}w-(\nabla ·\tilde{\mathit{q}})I={\mathit{Q}}_{{\mathit{x}}^{2}z}{{|}_{z=0}-q_z|}_{z=0}\mathit{I}\\ +{\mathit{S}}_{\mathit{x}z}{{|}_{z=0}{(\mathit{u}{|}_{z=0}-\tilde{\mathit{u}})}^T+(\mathit{u}{|}_{z=0}-\tilde{\mathit{u}}){\mathit{S}}_{\mathit{x}z}|}_{z=0}^T-({\mathit{S}}_{\mathit{x}z}{|}_{z=0}·(\mathit{u}{|}_{z=0}-\tilde{\mathit{u}}))\mathit{I}\\ \hfill -{\int }_0^{\infty }\left(\nabla \overline{\varphi }{(\mathit{u}-\tilde{\mathit{u}})}^T+(\mathit{u}-\tilde{\mathit{u}})\nabla {\overline{\varphi }}^T-((\mathit{u}-\tilde{\mathit{u}})·\nabla \overline{\varphi })\mathit{I}\right)\mathrm{d}z.\end{array}$$

(45)

Following [9], we split the horizontal skewness tensor into the appropriate combination of the horizontal heat fluxes and the traceless deviator $\tilde{\mathit{Q}}$:

$$\int {(\mathit{v}-\tilde{\mathit{u}})}^3\tilde{f}\mathrm{d}\mathit{v}\mathrm{d}w=\frac{1}{2}\left(\tilde{\mathit{q}}I+{\left(\tilde{\mathit{q}}I\right)}^T+{\left(\tilde{\mathit{q}}I\right)}^{TT}\right)+\tilde{\mathit{Q}},$$

(46a)

$$\nabla ·\int {(\mathit{v}-\tilde{\mathit{u}})}^3\tilde{f}\mathrm{d}\mathit{v}\mathrm{d}w-(\nabla ·\tilde{\mathit{q}})I=\frac{1}{2}\left(\nabla \tilde{\mathit{q}}+\nabla {\tilde{\mathit{q}}}^T-(\nabla ·\tilde{\mathit{q}})I\right)+\nabla ·\tilde{\mathit{Q}}.$$

(46b)

In the same fashion as in our recent work [9], 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 stress $\tilde{S}$ to remain small, we, first, impose the balance of external forces and boundary effects in its equation (45), that is,

$${\mathit{Q}}_{{\mathit{x}}^{2}z}{{|}_{z=0}={\int }_0^{\infty }\left(\nabla \overline{\varphi }{(\mathit{u}-\tilde{\mathit{u}})}^T+(\mathit{u}-\tilde{\mathit{u}})\nabla {\overline{\varphi }}^T\right)\mathrm{d}z,q_z|}_{z=0}={\int }_0^{\infty }(\mathit{u}-\tilde{\mathit{u}})·\nabla \overline{\varphi }\mathrm{d}z,$$

(47)

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 (45) for the horizontal shear stress,

$$\tilde{p}(\nabla \tilde{\mathit{u}}+\nabla {\tilde{\mathit{u}}}^T)+\frac{1}{2}(\nabla \tilde{\mathit{q}}+\nabla {\tilde{\mathit{q}}}^T)+\nabla ·\tilde{\mathit{Q}}=\mathit{0},\tilde{p}\nabla ·\tilde{\mathit{u}}+\frac{1}{2}\nabla ·\tilde{\mathit{q}}=0,$$

(48)

where, again, the latter equation is the trace of the former. Substituting the above relations into the pressure equation, we arrive at

$$\frac{\partial \tilde{p}}{\partial t}+\tilde{\mathit{u}}·\nabla \tilde{p}=0,$$

(49)

that is, in the same manner as in [9], the pressure is preserved along the stream lines (balanced flow). For the vanishing shear stress, the momentum equation in (32) becomes

$$\frac{\partial \left(\tilde{\rho }\tilde{\mathit{u}}\right)}{\partial t}+\nabla ·\left(\tilde{\rho }{\tilde{\mathit{u}}}^2\right)+\nabla \tilde{p}=-\nabla {\int }_0^{\infty }\overline{\varphi }\mathrm{d}z.$$

(50)

The density equation in (31), together with the pressure equation (49) and the momentum equation (50), 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 (49) in the form of a conservation law. For this, we can use the inverse kinetic temperature as a variable; indeed, denoting $\tilde{\theta }=\tilde{p}/\tilde{\rho }$, one can, with the help of (31), rearrange (49) in the form  

$$\frac{\partial \left({\tilde{\theta }}^{-1}\right)}{\partial t}+\nabla ·\left({\tilde{\theta }}^{-1}\tilde{\mathit{u}}\right)=0.$$

(51)

Together, the equations (31), (50) and (51) form the system of nonlinear conservation laws for the density $\tilde{\rho }$, momentum $\tilde{\rho }\tilde{\mathit{u}}$ and inverse kinetic temperature ${\tilde{\theta }}^{-1}$, where the latter can be related to the standard temperature (in ${}^{\circ }$K) via an appropriate gas constant.

3.2. Approximation for the Vertically Integrated Mean Field Potential

To estimate the effect of the potential forcing in the momentum equation (50), we assume that the gas particles interact via the hard sphere mean field potential—that is, $\varphi \left(r\right)$ is zero when r is greater than the diameter of the protection sphere [20] (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 (2) equals the volume of this protection sphere regardless of the value of the kinetic temperature:

$$\int \left[1-exp\left(-\frac{\varphi (\parallel \mathit{y}\parallel )}{\theta }\right)\right]\mathrm{d}\mathit{y}=V_{prot}=8V_{HS}.$$

(52)

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 (2) and denoting the density of the hard sphere via ${\rho }_{HS}=m/V_{HS}$, we arrive at a rather simple expression for the hard sphere mean field potential ${\overline{\varphi }}_{HS}$:

$${\overline{\varphi }}_{HS}=\frac{4\rho p}{{\rho }_{HS}}.$$

(53)

Note that the hard sphere mean field potential in (53) introduces the same correction into the momentum equation in (4) as does the Enskog correction in the Enskog–Euler equations (see Abramov [19] for details).

Next, we assume that, in the vertical direction, the hydrostatic balance is achieved:

$$\frac{\partial p}{\partial z}=-g\rho =-\frac{gp}{\theta },\frac{\partial lnp}{\partial z}=-\frac{g}{\theta }{,p=p|}_{z=0}exp\left(-g{\int }_0^z\frac{\mathrm{d}\xi }{\theta \left(\xi \right)}\right).$$

(54)

We now express $\rho$ from the hydrostatic balance relation and write

$$\rho p=-\frac{p}{g}\frac{\partial p}{\partial z}=-\frac{1}{2g}\frac{\partial \left(p^2\right)}{\partial z}=\frac{{\left(p{|}_{z=0}\right)}^2}{\theta }exp\left(-2g{\int }_0^z\frac{\mathrm{d}\xi }{\theta \left(\xi \right)}\right).$$

(55)

For the integral of $\rho p$ in the vertical direction, we thus have   

$${\int }_0^{\infty }\rho p\mathrm{d}z={\left(p{|}_{z=0}\right)}^2{\int }_0^{\infty }exp\left(-2g{\int }_0^z\frac{\mathrm{d}\xi }{\theta \left(\xi \right)}\right)\frac{\mathrm{d}z}{\theta \left(z\right)}=\frac{{\left(p{|}_{z=0}\right)}^2}{2g}{\int }_0^{\infty }{e}^{-\eta }\mathrm{d}\eta =\frac{{\left(p{|}_{z=0}\right)}^2}{2g}.$$

(56)

   Above, the variable $\eta =\eta \left(z\right)$ is given via

$$\eta \left(z\right)=2g{\int }_0^z\frac{\mathrm{d}\xi }{\theta \left(\xi \right)},$$

(57)

and we assume that $\theta$ 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 ${\overline{\varphi }}_{HS}$, we thus have

$${\int }_0^{\infty }{\overline{\varphi }}_{HS}\mathrm{d}z=\frac{4}{{\rho }_{HS}}{\int }_0^{\infty }\rho p\mathrm{d}z=\frac{2{\left(p{|}_{z=0}\right)}^2}{g{\rho }_{HS}}.$$

(58)

What remains to be done is to evaluate the surface pressure ${p|}_{z=0}$. For this, 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

$${p|}_{z=0}=g{\int }_0^{\infty }\rho \mathrm{d}z=g\int \left({\int }_0^{\infty }f\mathrm{d}z\right)\mathrm{d}\mathit{v}\mathrm{d}w=g\int \tilde{f}\mathrm{d}\mathit{v}\mathrm{d}w=g\tilde{\rho },$$

(59)

and, subsequently,

$${\int }_0^{\infty }{\overline{\varphi }}_{HS}\mathrm{d}z=\frac{2g{\tilde{\rho }}^2}{{\rho }_{HS}}.$$

(60)

As a result, the horizontal momentum equation for the hard sphere gas becomes

$$\frac{\partial \left(\tilde{\rho }\tilde{\mathit{u}}\right)}{\partial t}+\nabla ·\left(\tilde{\rho }{\tilde{\mathit{u}}}^2\right)+\nabla \left(\tilde{p}+\frac{2g{\tilde{\rho }}^2}{{\rho }_{HS}}\right)=\mathit{0}.$$

(61)

4. 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 the Earth’s atmosphere. In particular, the acceleration constant in (61) is set to $g=9.81$ m/s${}^2$, which corresponds to Earth’s gravity, while the vertically integrated density $\tilde{\rho }$ and temperature $\tilde{\theta }$ 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’s atmosphere using our model, we need to “create” a hypothetical hard sphere gas with a density that matches the 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 (31), (49), (51) and (61) 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 [21,22,23,24]. This is not necessarily a contradiction—clearly, in a real gas, both types of interactions are present (for example, while potential interactions are time-reversible, 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 a 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 [21,22,24]:

$$\mu =\frac{5}{16\sqrt{\pi }}\frac{m}{{\sigma }^2}\sqrt{\theta }=\frac{5}{16\sqrt{\pi }}\frac{M}{N_A{\sigma }^2}\sqrt{\theta },$$

(62)

where $\sigma$ is the diameter of the hard sphere, $N_A=6.02214\times {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 ${\rho }_{HS}$, we observe that

$$\theta =\frac{RT}{M},{\rho }_{HS}=\frac{m}{V_{HS}}=\frac{6m}{\pi {\sigma }^3}=\frac{6M}{\pi N_A{\sigma }^3},$$

(63)

where $R=8.31446$ kg m${}^2$/K mol s${}^2$ is the universal gas constant, and T is the usual temperature in ${}^{\circ }$K. Expressing $\sigma$ via ${\rho }_{HS}$ from the latter equation, and substituting the expression for $\theta$, we obtain the formula for ${\rho }_{HS}$ via the viscosity $\mu$ and temperature T:

$${\rho }_{HS}=\frac{384}{5\sqrt{5}}{\left(\frac{N_A^{2}M{\mu }^6}{\pi R^3{T}^3}\right)}^{1/4}$$

(64)

With (64), we can “create” a theoretical hard sphere analog of any gas with a prescribed molar mass M and known viscosity $\mu$ at a given temperature T. For air, $M=2.897·{10}^{-2}$ kg/mol, and we take $\mu =1.8194·{10}^{-5}$ kg/m s at $T=293.15$ K [25], for which (64) yields ${\rho }_{HS}=1850$ kg/m${}^3$. This value is used throughout all computations below.

4.2. 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’s 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.

4.3. 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 $\tilde{\mathit{u}}$ of the flow in the domain, and on those boundaries where the Dirichlet boundary condition is specified. The way in which 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 $\tilde{\mathit{u}}$ is defined, we compute the kinetic temperature $\tilde{\theta }$ of the flow using the Bernoulli law for a compressible gas:

  $$\tilde{\theta }=\frac{R{\tilde{T}}_0}{M}-\frac{\gamma -1}{2\gamma }{\parallel \tilde{\mathit{u}}\parallel }^2,$$

  (65)

  where the background temperature of the gas is set to ${\tilde{T}}_0=250$ K—that is, the approximate vertically averaged temperature of air throughout the troposphere. Observe, however, that in the present setting, the gas is two-dimensional, and, therefore, the effective adiabatic exponent $\gamma =2$.
• Once $\tilde{\mathit{u}}$ and $\tilde{\theta }$ are defined, we specify the pressure $\tilde{p}$ (which automatically yields the density via $\tilde{\rho }=\tilde{p}/\tilde{\theta }$). Below, we study two scenarios: an inertial jet, and a cyclostrophic vortex. In the inertial jet scenario, we set the pressure to a constant value ${\tilde{p}}_0$,

  $${\tilde{p}}_0=\frac{R{\tilde{\rho }}_0{\tilde{T}}_0}{M},$$

  (66)

  where ${\tilde{\rho }}_0={10}^4$ kg/m${}^2$, which constitutes the approximate mass of the air column per unit area of the Earth’s surface. In the cyclostrophic vortex scenario, the pressure gradient is balanced by the centripetal force:

  $$\nabla \tilde{p}=-\tilde{\rho }\Omega \times (\Omega \times \mathit{r})=\tilde{\rho }{\parallel \Omega \parallel }^2\mathit{r},$$

  (67)

  where $\mathit{r}$ is the coordinate of the point relative to the center of rotation, $\Omega$ is the angular velocity of rotation, and the latter identity is valid for a planar rotation (that is, when $\mathit{r}\perp \Omega$). Here, the velocity $\tilde{\mathit{u}}$ is orthogonal to both $\mathit{r}$ and $\Omega$, that is,

  $$\tilde{\mathit{u}}=\Omega \times \mathit{r},\mathrm{and},\sin \mathrm{ce}\mathit{r}\perp \Omega ,\parallel \tilde{\mathit{u}}\parallel =\parallel \Omega \parallel \parallel \mathit{r}\parallel .$$

  (68)
  Upon substitution, this leads to

  $$\nabla \tilde{p}=\frac{\tilde{\rho }{\parallel \tilde{\mathit{u}}\parallel }^2}{{\parallel \mathit{r}\parallel }^2}\mathit{r}.$$

  (69)
  Assuming, in turn, that $\tilde{\rho }$, $\tilde{\mathit{u}}$ and $\tilde{p}$ do not depend on the angle, and only depend on $\parallel \mathit{r}\parallel$, we denote $r=\parallel \mathit{r}\parallel$, $\tilde{u}=\parallel \tilde{\mathit{u}}\parallel$, and arrive at the following explicit formula for $\tilde{p}$:

  $$\frac{\mathrm{d}\tilde{p}}{\mathrm{d}r}=\frac{\tilde{\rho }{\tilde{u}}^2}{r}=\frac{\tilde{p}{\tilde{u}}^2}{\tilde{\theta }r},\mathrm{or}\tilde{p}\left(r\right)={\tilde{p}}_{0}exp\left(-{\int }_r^{\infty }\frac{{\tilde{u}}^2\left(\xi \right)}{\tilde{\theta }\left(\xi \right)}\frac{\mathrm{d}\xi }{\xi }\right).$$

  (70)
  Above, in the integration formula, we assume that $\tilde{u}$ vanishes sufficiently far away from the center of rotation. Moreover, in practice, one can simplify the integration by presuming that $\tilde{\theta }$ varies weakly in comparison to ${\tilde{u}}^2$ under the integral, and factor the kinetic temperature out of the integration (as we do further below).

4.4. Numerical Methods and Software

In the current work, we use the same software as in our previous work [9]—namely, we use OpenFOAM [26] to perform all numerical simulations. Noting that equations (31), (51) and (61), 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 [27], which uses the central scheme of Kurganov and Tadmor [28] 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 20% of the maximal Courant number.

4.5. 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

$$\tilde{\mathit{u}}{|}_{t=0}=(u(y/d),0),u\left(z\right)=u_0\frac{1+cos\left(\pi z\right)}{2}(1-\alpha z),$$

(71)

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 $\alpha =0.05$ is introduced in the same manner as by McCalpin and Haidvogel [30], 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 $\mathit{Re}\sim {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 (71) 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 are 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 (50) (or its special case (61) 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 (31), (51) and (61) forward for 72 h in model time. The time step, chosen from the CFL criterion, varied between 12 and 16 s. In Figure 1, we show the speed $\parallel \tilde{\mathit{u}}\parallel$ of the flow, captured at 12, 18, 24 and 30 h elapsed model time, on contour plots. Here, we can see that small fluctuations in the jet manifest at 12 h. At 18 h, the jet visibly meanders roughly between 1000 and 2000 km. At 24 and 30 h, 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 h resemble the famous Reynolds experiment [31]. After 30 h, the general configuration of the flow does not seem to visibly change any further, and thus we do not show the snapshots past this time.

The vorticity $\tilde{\omega }$ is given via the curl of velocity $\tilde{\mathit{u}}$:

$$\tilde{\omega }=\nabla \times \tilde{\mathit{u}}.$$

(72)

Since $\tilde{\mathit{u}}$ is confined to the $xy$-plane, $\tilde{\omega }$ is orthogonal to this plane, and only has the z-component, given via

$${\tilde{\omega }}_z=\frac{\partial {\tilde{u}}_y}{\partial x}-\frac{\partial {\tilde{u}}_z}{\partial y}.$$

(73)

In Figure 2, we show the contour plots of ${\tilde{\omega }}_z$ for the inertial flow, also captured at 12, 18, 24 and 30 h 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 h, and visible meandering manifests at 18 h. At 24 h, the structure of vorticity resembles a Kármán vortex trail [32]; at 30 h, 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 [9]: first, the kinetic energy $E\left(\mathit{x}\right)={\parallel \mathit{u}\left(\mathit{x}\right)\parallel }^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 [8], 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$ h. As we can see in Figure 3, Kolmogorov’s famous $k^{-5/3}$-decay of the kinetic energy spectrum occurs roughly between the wavenumbers 8 and 50; for smaller scales, the $k^{-8/3}$-decay is observed. Similar spectra, with two different powers of decay at different ranges of wavenumbers, have been registered in the Jovian atmosphere [33,34].

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 $\tilde{\rho }$ 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 h. Observe that, in a fully developed chaotic flow, the density varies between $4\times {10}^3$ and $1.6\times {10}^4$ kg/m${}^2$ (with the background value ${10}^4$ kg/m${}^2$), which, of course, does not happen in the 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 this region, and the usual compressible Euler equations would apply instead. Thus, in order to correct such 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.

4.6. 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

$$\tilde{u}\left(r\right)=\left\{\begin{array}{cc}{\displaystyle u_0\frac{4ln(r/r_0)ln(R_0/r)}{{ln}^2(R_0/r_0)},}\hfill & r_0\le r\le R_0,\hfill \\ \mathrm{zero},\hfill & \mathrm{otherwise}.\hfill \end{array}\right.$$

(74)

Here, the maximum speed $u_0=40$ m/s is achieved at the geometric mean distance $r_{\max u}=\sqrt{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}\approx 158$ km). The reason that $\tilde{u}\left(r\right)$ is chosen in this manner, is that it is, effectively, a function of $lnr$, which synergizes well with the computation of the integral in the cyclostrophic pressure profile formula (70). 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 s. The flow remained largely laminar for the initial 150 min of the elapsed model time, after which turbulent motions rapidly developed around the region of the maximal flow speed. Shortly after 200 min 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 min 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 [16]. 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 $\parallel \tilde{\mathit{u}}\parallel$ of the flow, captured at 160, 170, 180 and 190 min of the elapsed model time, on contour plots. Here, we can see that small speed fluctuations manifest at 160 min, and become progressively larger at 170, 180 and 190 min. 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 ${\tilde{\omega }}_z$ of the flow, captured at 160, 170, 180 and 190 min 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 that, 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 min 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. A similar decay of the kinetic energy spectrum at moderate and small scales has been observed in the Earth’s atmosphere [8].

In Figure 8 and Figure 9, we show the contour plots of the pressure $\tilde{p}$ and density $\tilde{\rho }$, respectively, of the cyclostrophic vortex, captured at 160, 170, 180 and 190 min 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.

5. Summary

In the current work, we investigate the ability of our model of a balanced compressible hard sphere gas flow [9] 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’s 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, compresses 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 [31]. 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 [9], 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.