1 Introduction

The goal of this section is to compute the natural convection patterns in the Rama atmosphere, as well as the wind speed. Indeed, we assumed in the previous section that the atmosphere was at rest, but this is not compatible with the temperature variations we found in the end: these variations will inevitably lead to winds, due to the buoyancy force. The winds will modify these temperature variations through temperature advection, as well as the density profile, through matter advection. This, in turn, will modify how short wave and long wave radiations are scattered inside the atmosphere, thus modifying the temperature gradients, and then the buoyancy force.

The above discussion shows that, in theory, radiative transfer and fluid dynamics are completely coupled here. In practice, however, we found in the previous section that the temperature variations, in an atmosphere at rest, are very small. And we can expect the winds to reduce the temperature gradients even more, by transporting energy from the hot to the cold areas. Then, using the ideal gas law, the density variations should be equally small. We can thus safely assume that the density profile is the one computed in the axisymmetric case, with very small variations that will not notably influence the radiative transfer solution (i.e. the energy $\bar{J}$ received at each point in the atmoshere by radiative transfer).

In other words, we can simulate the winds by using the fluid dynamics equations alone, with a heat source term $\bar{J}$ given by our radiative transfer computations for an atmosphere at rest, computed in the previous sections. This decouples the two sets of equations, thereby greatly simplifying their resolution.

This is what we do in the following. We first briefly recall the Euler equations for ideal fluids. We continue by writing these equations for the Rama case, and then simplify them by using the Boussinesq approximation. We then discuss the numerical simulation method we use to solve these equations, and finally present our results.

2 Euler equations

The general fluid dynamics equations are the Navier-Stokes equations, but we assume here that the atmosphere is an ideal fluid without viscosity (and without thermal conduction either). In this case these equations reduce to the Euler equations. In an inertial reference frame, with an external force density $\mathbf{F}$ and an external energy flux $\mathbf{E}$ (in $W.m^{-2}$), these equations are the mass continuity, the momentum and the internal energy equations (here in conservation form): \begin{align} \frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\bu)&=0\label{eq:cont0}\\ \frac{\partial\rho\bu}{\partial t}+\nabla\cdot(\rho\bu\otimes\bu)&=-\nabla p+\rho\mathbf{F}\label{eq:momentum0}\\ \frac{\partial\rho e}{\partial t}+\nabla\cdot(\rho\bu e)&=-p\nabla\cdot\bu-\nabla\cdot\mathbf{E}\label{eq:internal0} \end{align} where $\bu$ is the velocity, and $e$ is the specific internal energy. These equations must be completed with boundary conditions, which are specific to each case. In a volume $V$ enclosed by walls, like the Rama atmosphere, the most important boundary condition is that the fluid can not flow through the walls, which means that $\bu\cdot\bn=0$ on the boundary $\partial V$ of $V$. And for a viscuous fluid, an even stronger no-slip constraint holds: $\bu=0$ on $\partial V$.

Energy conservation

Note that the Euler equations conserve total energy. Indeed, by multiplying the momentum equation \eqref{eq:momentum0} by $\bu$ and by adding the internal energy equation \eqref{eq:internal0} to the result, we get \begin{align*} \frac{\partial\rho(\frac{1}{2}\bu\cdot\bu+e)}{\partial t}+\nabla\cdot\left(\rho\bu\left(\frac{1}{2}\bu\cdot\bu+e\right)\right)&=-\nabla \cdot(p\bu)+\rho\mathbf{F}\cdot\bu-\nabla\cdot\mathbf{E} \end{align*} where we used the identity $f\frac{\partial fg}{\partial x}=\frac{1}{2}\frac{\partial f^2g}{\partial x}+\frac{f^2}{2}\frac{\partial g}{\partial x}$ and the continuity equation \eqref{eq:cont0} to rewrite the $\bu\cdot\nabla\cdot(\rho\bu\otimes\bu)$ term. Integrated over a volume $V$ enclosed by walls, and using the divergence theorem, this equation gives: \begin{align*} \frac{\partial}{\partial t}\int_V\rho\left(\frac{1}{2}\bu\cdot\bu+e\right)\diff v&=\int_V(\rho\mathbf{F}\cdot\bu-\nabla\cdot\mathbf{E})\diff v \end{align*} Thus, in the absence of external forces and external energy fluxes, the integral of the sum of the kinetic energy and of the internal energy remains constant over time. Also, if the forces are conservative, i.e. if $\mathbf{F}=-\nabla f$ with $f$ independent of $t$, then we have, using the continuity equation, $\rho\mathbf{F}\cdot\bu=-\nabla\cdot(\rho f\bu)-\frac{\partial \rho f}{\partial t}$. Thus energy is also conserved in this case, if we define it as the sum of the kinetic energy $\frac{1}{2}\rho\bu\cdot\bu$, internal energy $\rho e$ and potential energy $\rho f$.

Finally, in the absence of external forces, the Euler equations also conserve the total momentum $\int_V\rho\bu\diff v$ (as the name of the "momentum conservation" equation implies) and the total angular momentum $\int_V\bp\times\rho\bu\diff v$.

3 Euler equations for Rama

In order to apply the Euler equations to Rama, we must compute the external forces $\mathbf{F}$ and the external energy flux $\mathbf{E}$ acting on the atmosphere, as well as the internal energy $e$:

in an inertial frame, there are no external forces $\mathbf{F}$ acting on the atmosphere. However, we want to write the Euler equations in the Rama rotating frame. This frame is not an inertial frame, and writing the coordinate transform from a non rotating to the rotating frame yields additional, fictituous forces which are the centrifugal force $\rho r\omega^2\be_r$ and the Coriolis force $-2\rho\omega\be_z\times\bu$.

the energy flux $\mathbf{E}$ comes from the radiance in the short and long wave range. More precisely $\mathbf{E}$ is equal to the radiative vector flux $\mathbf{L}$ defined in the Atmosphere section. The divergence of this flux has already been computed in this section, and was found equal to $\nabla\cdot\mathbf{E}=\nabla\cdot\mathbf{L}=-4\pi k_e\bar{J}+4k_e\sigma T^4$. Note that without winds this term must be $0$ in a steady case, as can be seen from the internal energy equation, and as we assumed in the previous sections. In the presence of winds, however, it may not be $0$.
As in the previous sections we assume that the atmosphere is an ideal gas. We also assume here that the atmosphere is dry, i.e. that it does not contain water vapor. In this case the internal energy $e$ is equal to $c_v T$, where $c_v$ is the heat capacity at constant volume of air.

Reporting this in Eqs. \eqref{eq:cont0}-\eqref{eq:internal0} gives the Euler equations for Rama: \begin{align} \frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\bu)&=0\label{eq:cont1}\\ \frac{\partial\rho\bu}{\partial t}+\nabla\cdot(\rho\bu\otimes\bu)&=-\nabla p+\rho r\omega^2\be_r-2\omega\be_z\times\rho\bu\label{eq:momentum1}\\ \frac{\partial\rho c_v T}{\partial t}+\nabla\cdot(\rho\bu c_v T)&=-p\nabla\cdot\bu+4\pi k_e\bar{J}-4k_e\sigma T^4\label{eq:internal1} \end{align} where $\bar{J}$ was computed in the previous section (as discussed in the introduction, in theory, $\bar{J}$ depends on $\rho$ and $T$ which are the unknowns, with $\bu$ and $p$, of this system of equations. In practice, however, we showed above that we can use the value of $\bar{J}$ computed for the isothermal hydrostatic density $\rho_0$ and for an atmosphere at rest).

Although we assumed a fluid without viscosity to simplify the equations, we do use the no-splip boundary condition $\bu=0$ on $\partial V$. We do this in order to exclude the solution where the atmosphere would be at rest in an inertial frame, with Rama rotating around it without friction ($\bu=-r\omega\be_{\theta},\rho=cst,p=cst$). Finally, note that since Rama is invariant under rotations of 120°, the solution to these equations should also be invariant under the same rotations. This implies, in particular, that the velocity $\bu$ on the rotation axis must be $0$.

Energy conservation

Note that the centrifugal force is conservative, because $r\omega^2\be_r=\nabla \frac{1}{2}r^2\omega^2$. Note also that the Coriolis force does not change the kinetic energy, since $-(2\omega\be_z\times\rho\bu)\cdot\bu=0$. Thus, in the absence of radiative heating, i.e. if $4\pi k_e\bar{J}-4k_e\sigma T^4=0$, energy is conserved.

Since there are no external forces in the non rotating frame, momentum and angular momentum are also conserved in this frame. Then, since angular momentum in the inertial and rotating frames only differ by a constant, angular momentum is preserved in the Rama rotating frame as well. Finally, if the total momentum is 0 in the inertial frame, then it is also conserved in the rotating frame.

4 Anelastic Boussinesq approximation

The system of equations in Eqs. \eqref{eq:cont1}-\eqref{eq:internal1} is hard to simulate numerically because it allows sound waves, which propagate very fast. Then a very small time step must be used to solve the equations without numerical issues like oscillations or instabilities. To avoid this problem, and since sound waves play no role in natural convection, we use the Boussinesq approximation to filter them out. In other words we note $p=p_0+\delta p$, $\rho=\rho_0+\delta\rho$, and $T=T_0+\delta T$, where $p_0$, $\rho_0$, $T_0$ is the solution of isothermal hydrostatic equilibrium (i.e. $T_0$ constant independent of $\bp$ and $\nabla p_0=\rho_0r\omega^2\be_r$) and we assume that $|\delta p|\ll p_0$, $|\delta \rho|\ll \rho_0$, $|\delta T|\ll T_0$. We then approximate $\rho$ by $\rho_0$ eveywhere, except in the centrifugal force term — see The Anelastic and Boussinesq Approximations [Randall 2013] for a detailed derivation of this approximation. The equations become \begin{align} \nabla\cdot(\rho_0\bu)&=0\\ \frac{\partial\rho_0\bu}{\partial t}+\nabla\cdot(\rho_0\bu\otimes\bu)&=-\nabla p+\rho r\omega^2\be_r-2\omega\be_z\times\rho_0\bu\label{eq:momentum2}\\ \frac{\partial\rho_0 c_v T}{\partial t}+\nabla\cdot(\rho_0\bu c_v T)&=-p\nabla\cdot\bu+4\pi k_e\bar{J}-4k_e\sigma T^4 \end{align}

Using the above hypotheses, the first two right hand side terms in the momentum equation \eqref{eq:momentum2} can be rewritten as (following the derivation of [Randall 2013]): \begin{align*} -\nabla p+\rho r\omega^2\be_r&=-\rho\left(\frac{1}{\rho_0+\delta\rho}\nabla(p_0+\delta p)- r\omega^2\be_r\right)\\ &=-\rho\left(\frac{1}{\rho_0+\delta\rho}\nabla\delta p-\left(\frac{\rho_0}{\rho_0+\delta\rho}-1\right)r\omega^2\be_r\right)\\ &\approx-\rho_0\left(\frac{1}{\rho_0}\nabla\delta p-\frac{\delta\rho}{\rho_0}r\omega^2\be_r\right)=-\rho_0\left(\nabla\frac{\delta p}{\rho_0}+\frac{\delta p}{\rho_0^2}\nabla\rho_0-\frac{\delta\rho}{\rho_0}r\omega^2\be_r\right)\\ \end{align*} Now, taking the gradient of the logarithm of the perfect gas equation, we get $\frac{1}{p_0}\nabla p_0=\frac{1}{\rho_0}\nabla \rho_0+\frac{1}{T_0}\nabla T_0=\frac{1}{\rho_0}\nabla \rho_0$ since $T_0$ is constant. Also, equating the first order terms in $p_0+\delta p=(\rho_0+\delta\rho)R_s(T_0+\delta T)$ yields $\frac{\delta p}{p_0}=\frac{\delta\rho}{\rho_0}+\frac{\delta T}{T_0}$. Substituting these relations in the above equation gives the gradient of the dynamic pressure and the buoyancy force: \begin{align*} -\nabla p+\rho r\omega^2\be_r&=-\rho_0\nabla\frac{\delta p}{\rho_0}-\rho_0r\omega^2\frac{\delta T}{T_0}\be_r \end{align*}

Finally, using $T=T_0+\delta T$, the fact that $T_0$ is constant, and the continuity equation, the internal energy equation becomes: \begin{align*} \frac{\partial\rho_0 c_v \delta T}{\partial t}+\nabla\cdot(\rho_0\bu c_v \delta T)&=-p\nabla\cdot\bu+4\pi k_e\bar{J}-4k_e\sigma T^4\\ &\approx-p\nabla\cdot\bu+4\pi k_e\delta J-16k_e\sigma T_0^3\delta T \end{align*} where $\delta J\eqdef\bar{J}-\sigma T_0^4/\pi$.

Note that the $-p\nabla\cdot\bu$ term is problematic here. Indeed, in the full, non approximated equations, it normally combines, in the total energy equation, with the $-\nabla p\cdot\bu$ term coming from the kinetic energy equation to give $-\nabla\cdot(p\bu)$, which is conservative. Here, however, due to the above approximations, the pressure gradient $\nabla p$ has a new form which does not combine with $-p\nabla\cdot\bu$ into a conservative term. The result is that the above equations do not conserve energy (as shown in [Randall 2013]).

This non conservation is a problem for numerical simulations, especially for long quasi static simulations required to reach a steady state. Indeed, it means that energy can accumulate or dissipate over time, and thus exceed any other error due to other approximations. Therefore, in order to enforce energy conservation, as shown below, we simply ignore the $-p\nabla\cdot\bu$ term in the internal energy equation.

Putting everything together, we obtain the 3 anelastic Boussinesq equations that we want to simulate in order to find the atmospheric convection patterns in Rama (for the 3 unknowns $\bu$, $\delta T$ and $\delta p$): \begin{align} \nabla\cdot(\rho_0\bu)&=0\label{eq:cont3}\\ \frac{\partial\rho_0\bu}{\partial t}+\nabla\cdot(\rho_0\bu\otimes\bu)&=-\rho_0\nabla\frac{\delta p}{\rho_0}-\rho_0r\omega^2\frac{\delta T}{T_0}\be_r-2\omega\be_z\times\rho_0\bu\label{eq:momentum3}\\ \frac{\partial\rho_0 c_v \delta T}{\partial t}+\nabla\cdot(\rho_0\bu c_v \delta T)&=4\pi k_e\delta J-16k_e\sigma T_0^3\delta T\label{eq:internal3} \end{align}

Or, equivalently, in cylindrical coordinates: \begin{align} &\frac{1}{r}\frac{\partial r \rho_0 u_r}{\partial r} +\frac{1}{r}\frac{\partial\rho_0 u_{\theta}}{\partial \theta}=0\label{eq:cont4}\\ &\begin{split} \frac{\partial \rho_0u_r}{\partial t} +\frac{1}{r}\frac{\partial r\rho_0 u_r u_r}{\partial r} &+\frac{1}{r}\frac{\partial\rho_0 u_{\theta} u_r}{\partial\theta} -\frac{1}{r}\rho_0 u_{\theta} u_{\theta}\\ &=-\rho_0\frac{\partial}{\partial r}\frac{\delta p}{\rho_0} -\rho_0r\omega^2\frac{\delta T}{T_0} +2\omega\rho_0 u_{\theta} \end{split}\label{eq:momentumr}\\ &\begin{split} \frac{\partial \rho_0u_{\theta}}{\partial t} +\frac{1}{r}\frac{\partial r\rho_0 u_r u_{\theta}}{\partial r} &+\frac{1}{r}\frac{\partial\rho_0 u_{\theta} u_{\theta}}{\partial\theta} +\frac{1}{r}\rho_0 u_{\theta} u_r\\ &=-\rho_0\frac{1}{r}\frac{\partial}{\partial \theta}\frac{\delta p}{\rho_0} -2\omega\rho_0 u_r \end{split}\label{eq:momentumt}\\ &\begin{split} \frac{\partial \rho_0 c_v \delta T}{\partial t} +\frac{1}{r}\frac{\partial r \rho_0 u_r c_v \delta T}{\partial r} &+\frac{1}{r}\frac{\partial\rho_0 u_{\theta} c_v \delta T}{\partial \theta}\\ &=4\pi k_e\delta J-16k_e\sigma T_0^3\delta T \end{split}\label{eq:internal4} \end{align} with $u_r=u_{\theta}=0$ at $r=0$ and $r=r_s$, and Neumann boundary conditions $\partial\delta T/\partial r=0$ for the temperature.

Energy conservation

We can check that the above anelastic Boussinesq equations conserve energy, as follows. Multiplying the momentum equation \eqref{eq:momentum3} by $\bu$ and using the continuity equation \eqref{eq:cont3} yields: \begin{align*} \frac{\partial\rho_0\frac{1}{2}\bu\cdot\bu}{\partial t}+\nabla\cdot\left(\rho_0\bu\left(\frac{1}{2}\bu\cdot\bu\right)\right)&=-\nabla\cdot(\delta p\bu)-\rho_0r\omega^2\frac{\delta T}{T_0}u_r\\ \end{align*} Similarly, multiplying the internal energy equation \eqref{eq:internal3} by $\frac{1}{2}r^2\omega^2\frac{1}{T_0}$ yields (in the absence of external radiative heating): \begin{align*} \frac{\partial\frac{1}{2}r^2\omega^2\rho_0 c_v\frac{\delta T}{T_0}}{\partial t}+\nabla\cdot\left(\frac{1}{2}r^2\omega^2\rho_0\bu c_v\frac{\delta T}{T_0}\right)&=\rho_0r\omega^2\frac{\delta T}{T_0}u_r\\ \end{align*} Finally, summing these two equations and the internal energy equation gives (in the absence of external radiative heating): \begin{align*} \frac{\partial e_T}{\partial t}+\nabla\cdot(\bu e_T)&=-\nabla\cdot(\delta p\bu)\\ \end{align*} where the total energy $e_T$ is the sum of the kinetic energy, the internal energy and the potential energy $\frac{1}{2}r^2\omega^2\rho_0 c_v\frac{\delta T}{T_0}$: \begin{equation} e_T\eqdef \frac{1}{2}\rho_0\bu\cdot\bu + \rho_0c_v\delta T + \frac{1}{2}r^2\omega^2\rho_0 c_v\frac{\delta T}{T_0} \end{equation} This shows, as we wanted, that the integral of the total energy $e_T$ over the whole atmosphere remains constant over time, in the absence of external radiative heating.

It is also possible to show that the angular momentum conservation property ($\int_Vr\rho_0u_{\theta}dv=cst$) is also preserved by the anelastic Boussinesq equations, despite the approximations used to derive them. This can be done by multiplying Eq. \eqref{eq:momentumt} by $r$ and by rewriting the $\frac{\partial r\rho_0 u_r u_{\theta}}{\partial r}$ term as $\frac{1}{r}\frac{\partial r^2\rho_0 u_r u_{\theta}}{\partial r}-\rho_0u_ru_{\theta}$. The second term then cancels with the $\rho_0u_{\theta}u_r$ convection term, while the remaining terms are in conservation form and integrate to 0 on $V$.

5 Numerical simulation

Introduction

The non linear equations \eqref{eq:cont4}-\eqref{eq:internal4} can't be solved analytically, except in very specific cases (such as hydrostatic equilibrium conditions, which require $\delta J=0$). We must therefore use some numerical method to solve them. Many such methods exist. Some methods model the fluid with discrete particles moving along the flow, like the smoothed particle hydrodynamics method. But the majority of methods use a structured or unstructured grid, which does not move with the flow, to discretize the fluid volume. Examples of such methods include finite difference, finite element and finite volume methods.

In the case of Rama the fluid satisfies the anelastic incompressibility constraint $\nabla\cdot(\rho_0\bu)=0$, which is not easy to enforce with particle-based methods. We therefore use a grid-based method. And since the domain $V$ is a simple 2D disc, a rectangular structured grid in the $r,\theta$ polar coordinates space seems the most natural choice (because it represents the boundary $\partial V$ exactly; this is not possible with structured or unstructured grids in the $x,y$ cartesian coordinates space).

Two methods can be used to discretize the fluid properties in a structured or unstructured grid: the collocated and staggered schemes. In the collocated scheme all the variables are discretized at the grid cell centers, whereas in the staggered scheme the scalar variables are discretized at the cell centers, but the (normal component of the) velocity is discretized at cell face centers. The collocated scheme is simpler but suffers from numerical instability issues (the so called odd-even decoupling). We therefore use a staggered grid here, with a uniform spacing for simplicity.

For better stability and accuracy, we want a discretization scheme which conserves energy, like the continuous equations. Morinishi et al. propose such a scheme for structured staggered grids in cylindrical coordinates in their 2003 paper Fully conservative finite difference scheme in cylindrical coordinates for incompressible flow simulations. We therefore choose this numerical method to solve our problem. This method, however, is limited to fluids with a uniform density, which is not the case in Rama. Also, this paper does not present how to discretize the buoyancy force, the Coriolis force, and the internal energy equation. Extending the Morinishi et al. method to handle this is relatively straightforward, though, and we show how to do this below. On the other hand, Morinishi et al. explain how to handle singularities at the axis, but we don't have this issue here since, by symmetry, the velocity is 0 on the Rama axis.

Staggered grid

For the reasons presented above, we use a uniform grid in $r,\theta$ space, with the radial velocity $u_r$ discretized at "vertical" cell face centers, the tangential velocity $u_{\theta}$ discretized at "horizontal" cell face centers, and the pressure and temperature discretized at cell centers (see Fig. 1). In addition, constant values such as $r$ and $\rho_0$ are evaluated exactly at any point where we need them.

Figure 1: An example of a staggered grid used in our method, for $N_r=N_{\theta}=3$.

Thanks to the symmetry of the problem under rotations of 120° we discretize a $2\pi/3$ interval only, and use periodic boundary conditions for the $\theta=0$ and $\theta=2\pi/3$ boundaries. We divide the domain in $N_r \times N_{\theta}$ cells of the same size $\Delta r\eqdef r_s/N_r$ and $\Delta \theta\eqdef 2\pi/3N_{\theta}$. We note $r_i\eqdef (i-1)\Delta r$ and $\theta_j\eqdef (j-1)\Delta \theta$, where $i$ and $j$ can be half-integers, and note a discretized quantity $f$ with $i,j$ indices, e.g. $f_{i,\,j+\h}\eqdef f(r_i,\theta_{j+\h})$. Note that we have $N_rN_{\theta}$ values for $u_{\theta}$, $\delta p$ or $\delta T$, but that we have $(N_r+1)N_{\theta}$ values for $u_r$.

As in [Morinishi et al. 2003], we define the following central difference operators to approximate $\partial / \partial r$ and $\partial / \partial \theta$, as well as the following averaging operators, used to interpolate quantities at points where they are not discretized: \begin{align*} (\partial_r f)_{i,\,j}&\eqdef \frac{f_{i+\h,\,j}-f_{i-\h,\,j}}{\Delta r}&(\overline{f}^r)_{i,\,j}\eqdef\frac{f_{i+\h,\,j}+f_{i-\h,\,j}}{2}\\ (\partial_{\theta} f)_{i,\,j}&\eqdef \frac{f_{i,\,j+\h}-f_{i,\,j-\h}}{\Delta \theta}&(\overline{f}^{\theta})_{i,\,j}\eqdef\frac{f_{i,\,j+\h}+f_{i,\,j-\h}}{2} \end{align*}

Note that applying these operators on the boundaries of the grid requires values outside the grid. The values of the "top" and "bottom" rows are determined from the periodic boundary conditions in $\theta$. Similarly, the "left" and "right" columns are determined from the boundary conditions:

The Dirichlet boundary condition $u_{\theta}(r_s,\theta)=0$ imposes $u_{\theta,N_r+\frac{3}{2},j}=-u_{\theta,N_r+\h,j}$. Similarly, $u_{\theta,\h,j}=-u_{\theta,\frac{3}{2},j}$.
The Neumann boundary condition $\partial{\delta T}/{\partial r}(r_s,\theta)=0$ imposes $\delta T_{N_r+\h,j}=\delta T_{N_r-\h,j}$. Similarly, $\delta T_{\h,j}=\delta T_{\frac{3}{2},j}$.
Finally, the values $u_{r,0,j}$ and $u_{r,N_r+2,j}$ are determined by requiring that the continuity equation is satisfied in the cells outside the grid. This gives $u_{r,N_r+2,j}=u_{r,N_r,j}\frac{\rho_{0,N_r+\frac{3}{2}}}{\rho_{0,N_r+\h}}\frac{r_{N_r}\rho_{0,N_r}}{r_{N_r+2}\rho_{0,N_r+2}}$ and a similar relation for $u_{r,0,j}$.

Numerical method

Given the discretization scheme and the notations presented above, we discretize the spatial components in Eqs. \eqref{eq:cont4}-\eqref{eq:internal4} as follows, based on [Morinishi et al. 2003], extended here to handle a constant non uniform density $\rho_0$, the buoyancy and Coriolis forces, and the internal energy equation. The precise form of each term has been carefully chosen to preserve energy conservation, as shown in the next sub section (the continuity and internal energy equations are evaluated at cell centers, while the momentum equations are evaluated at cell face centers): \begin{align} &\frac{1}{r}\partial_r(r\rho_0 u_r)+\frac{1}{r}\partial_{\theta}(\rho_0 u_{\theta})=0\label{eq:cont5}\\ &\begin{split} \frac{\partial \rho_0u_r}{\partial t}+ \frac{1}{r}\partial_r \overline{r\rho_0 u_r}^r \overline{u_r}^r &+\frac{1}{r}\partial_{\theta}\overline{\rho_0 u_{\theta}}^r \overline{u_r}^{\theta} -\frac{1}{r}\overline{\overline{\rho_0 u_{\theta}}^{\theta} \overline{u_{\theta}}^{\theta}}^r\\ &=-\rho_0\partial_r\frac{\delta p}{\rho_0} -\rho_0r\omega^2\frac{\overline{\delta T}^r}{T_0} +2\omega\overline{\rho_0 \overline{u_{\theta}}^{\theta}}^r \end{split}\label{eq:momentumr5}\\ &\begin{split} \frac{\partial \rho_0u_{\theta}}{\partial t} +\frac{1}{r}\partial_r\overline{r\rho_0 u_r}^{\theta}\overline{u_{\theta}}^r &+\frac{1}{r}\partial_{\theta}\overline{\rho_0 u_{\theta}}^{\theta} \overline{u_{\theta}}^{\theta}+ \frac{1}{r}\overline{\overline{\rho_0 u_{\theta}}^{\theta}\overline{u_r}^r}^{\theta}\\ &=-\rho_0\frac{1}{r}\partial_{\theta}\frac{\delta p}{\rho_0} -2\omega\frac{1}{r}\overline{\rho_0 \overline{r u_r}^r}^{\theta} \end{split}\label{eq:momentumt5}\\ &\begin{split} \frac{\partial \rho_0 c_v \delta T}{\partial t} +\frac{1}{r}\partial_r(r \rho_0 u_r c_v &\overline{\delta T}^r) +\frac{1}{r}\partial_{\theta}(\rho_0 u_{\theta} c_v \overline{\delta T}^{\theta})\\ &=4\pi k_e\delta J-16k_e\sigma T_0^3\delta T \end{split}\label{eq:internal5} \end{align} Using the same notations as in [Morinishi et al. 2003], this can be summarized as: \begin{align*} (\mathrm{Cont.})&=0\\ \frac{\partial \rho_0u_r}{\partial t}+(\mathrm{Conv.})_r&=-(\mathrm{Pres.})_r-(\mathrm{Buo.})_r-(\mathrm{Corio.})_r\\ \frac{\partial \rho_0u_{\theta}}{\partial t}+(\mathrm{Conv.})_{\theta}&=-(\mathrm{Pres.})_{\theta}-(\mathrm{Corio.})_{\theta}\\ \frac{\partial \rho_0c_v\delta T}{\partial t} +(\delta T\mathrm{Conv.})&=(\mathrm{\delta T Source}) \end{align*}

We integrate these equations in time with a fourth order Runge-Kutta method. For each sub step inside a single Runge-Kutta time step $\Delta t$ we compute the velocity and temperature at $t'=t+\delta t$ using the semi-implicit scheme \begin{align*} \bu^{t'}&=\bu^t-\frac{\delta t}{\rho_0}\left((\mathrm{Conv.})^t+(\mathrm{Pres.})^{t'}+(\mathrm{Buo.})^t+(\mathrm{Corio.})^t\right)\\ \delta T^{t'}&=\delta T^t-\frac{\delta t}{\rho_0c_v}\left((\delta T\mathrm{Conv.})^t-(\mathrm{\delta T Source})^t\right) \end{align*} In order to compute the pressure term at $t'$ we require that $\bu^{t'}$ satisfies the continuity equation. This leads to the following Poisson equation: \begin{equation*} \begin{split} \partial_r\left(r\rho_0\partial_r \frac{\delta p}{\rho_0}\right) +\frac{1}{r}\partial_{\theta}\left(\rho_0\partial_{\theta} \frac{\delta p}{\rho_0}\right) =&\partial_r r\left(\frac{\rho_0u_r^t}{\delta t}-(\mathrm{Conv.})_r^t-(\mathrm{Buo.})_r^t-(\mathrm{Corio.})_r^t\right)\\ &+\partial_{\theta}\left(\frac{\rho_0u_{\theta}^t}{\delta t}-(\mathrm{Conv.})_{\theta}^t-(\mathrm{Corio.})_{\theta}^t\right)\eqdef g \end{split} \end{equation*} As in "A Simple Compact Fourth-Order Poisson Solver on Polar Geometry", Ming-Chih Lai 2001, we solve it by transforming it into a set of independent tri-diagonal linear systems by using a fast Fourier transform in $\theta$. Indeed, with the notations $q_{i,j}\eqdef(\delta p/\rho_0)_{i+\frac{3}{2},j+\frac{3}{2}}$ and $f_{i,j}\eqdef g_{i+\frac{3}{2},j+\frac{3}{2}}$, this equation can be written in expanded form as: \begin{equation} a_i q_{i-1,j} + b_i q_{i,j} + c_i q_{i+1,j} + d_i q_{i,j-1} + d_i q_{i,j+1} = f_{i,j}\label{eq:poisson1} \end{equation} where \begin{align*} a_i&\eqdef \frac{r_{i+1}\rho_{0,i+1}}{\Delta r^2}\\ b_i&\eqdef -(a_i+c_i+2d_i)\\ c_i&\eqdef \frac{r_{i+2}\rho_{0,i+2}}{\Delta r^2}\\ d_i&\eqdef \frac{\rho_{0,i+\frac{3}{2}}}{r_{i+\frac{3}{2}}\Delta\theta^2} \end{align*} Consider now the Fourier transforms $\tilde{q}_{i,j}\eqdef\sum_{k=0}^{N_{\theta}-1}q_{i,k}w^{-jk}$ and $\tilde{f}\eqdef\sum_{k=0}^{N_{\theta}-1}f_{i,k}w^{-jk}$ where $w=e^{2i\pi/N_{\theta}}$. The above equation can then be rewritten by using the inverse Fourier transform of $\tilde{q}$ and $\tilde{f}$, yielding: \begin{equation*} \sum_{k=0}^{N_{\theta}-1}\left(a_i \tilde{q}_{i-1,k} + b_i \tilde{q}_{i,k} + c_i \tilde{q}_{i+1,k} + d_i \tilde{q}_{i,k}w^{-k} + d_i \tilde{q}_{i,k}w^k\right)w^{jk} = \sum_{k=0}^{N_{\theta}-1}\tilde{f}_{i,k}w^{jk} \end{equation*} from which we deduce \begin{equation} a_i \tilde{q}_{i-1,k} + (b_i + 2d_i\cos\frac{2\pi k}{N_{\theta}}) \tilde{q}_{i,k} + c_i \tilde{q}_{i+1,k}=\tilde{f}_{i,k}\label{eq:poisson2} \end{equation} Each value of $k$ gives an independent tri-diagonal linear system. Solving these independent systems gives the values of $\tilde{q}$, and their inverse Fourier transform then give $q$, i.e. the pressure divided by the density.

Note: Eq. \eqref{eq:poisson1} involves pressure values outside the fluid, like $\delta p_{N_r+\frac{3}{2},j+\h}$. We can compute these values by using the radial momentum equation Eq. \eqref{eq:momentumr5} at $r=0$ and $r=r_s$, where $\partial \rho_0 u_r/\partial t=0$. Using the momentum equation at $r=r_s$, for instance, gives \begin{equation*} c_{N_r-1}q_{N_r,j}=c_{N_r-1}q_{N_r-1,j}-\frac{r_s}{\Delta r}\left((Conv.)_r+(Buo.)_r+(Corio.)_r\right)_{N_r+1,j+\frac{3}{2}} \end{equation*} Substituting this in Eq. \eqref{eq:poisson1}, we see that must add new terms in the definitions of $b_{N_r-1}$ and $f_{N_r-1,j}$, respectively $c_{N_r-1}$ and $r_s((Conv.)_r+\ldots)/\Delta r$. For the same reason, the definitions of $b_0$ and $f_{0,j}$ must be modified as well, using similar relations.

Energy conservation

In the absence of external radiative heating, the above discretized equations conserve the following discretized total energy expression (using the kinetic energy formulation of [Morinishi et al. 2003]): \begin{equation} e_T\eqdef \frac{1}{2r}\left(\overline{r\rho_0u_ru_r}^r+\overline{r\rho_0u_{\theta}u_{\theta}}^{\theta}\right)+\rho_0c_v\delta T+\frac{1}{2}\rho_0r^2\omega^2\frac{\delta T}{T_0} \end{equation}

This can be shown as follows. First of all the internal energy $\rho_0c_v\delta T$ is conserved because its time derivative, given by Eq. \eqref{eq:internal5}, is equal to a divergence term of the form $\frac{1}{r}\partial_r f+\frac{1}{r}\partial_{\theta} g$ where $f$ is $0$ on the $r$ boundaries, and $g$ is periodic in $\theta$. Thus the integral of the internal energy $e$ over the whole volume, $\sum_V e\, r \Delta r \Delta \theta=\Delta r \Delta \theta\sum_V \partial_r f+\partial_{\theta} g$ is equal to $0$. This is the discrete analog of the divergence theorem, which shows that the integral of any term of this form is $0$.

We can then compute the time derivative of the kinetic and potential energy terms, and show that they cancel each other. First, to get the time derivative of the kinetic energy, we compute $(\overline{ru_r\partial\rho_0u_r/\partial t}^r+\overline{ru_{\theta}\partial\rho_0u_{\theta}/\partial t}^{\theta})/r$ using Eqs. \eqref{eq:momentumr5}-\eqref{eq:momentumt5}:

The $(\mathrm{Conv.})$ terms become (following [Morinishi et al. 2003]): \begin{equation} \begin{split} &\frac{1}{r}\left(\overline{u_r\partial_r \overline{r\rho_0 u_r}^r \overline{u_r}^r}^r +\overline{u_r\partial_{\theta}\overline{\rho_0 u_{\theta}}^r \overline{u_r}^{\theta}}^r -\overline{u_r\overline{\overline{\rho_0 u_{\theta}}^{\theta} \overline{u_{\theta}}^{\theta}}^r}^r\right)\\ +&\frac{1}{r}\left(\overline{u_{\theta}\partial_r\overline{r\rho_0 u_r}^{\theta}\overline{u_{\theta}}^r}^{\theta} +\overline{u_{\theta}\partial_{\theta}\overline{\rho_0 u_{\theta}}^{\theta} \overline{u_{\theta}}^{\theta}}^{\theta} +\overline{u_{\theta}\overline{\overline{\rho_0 u_{\theta}}^{\theta}\overline{u_r}^r}^{\theta}}^{\theta}\right) \end{split}\label{eq:kinetic1derivative} \end{equation} Using the identity $\sum_V \overline{f\overline{g}^x}^x=\sum_V \overline{f}^xg$ (if $f$ and $g$ are periodic and/or $0$ on the boundary), we see that the integrals of the right most terms cancel each other. Using this same identity, the identity $f\partial_x(\overline{f}^x g)=\frac{1}{2}\partial_x(\widetilde{ff}^x g)+\frac{f^2}{2}\partial_x g$, where $\widetilde{fg}^r\eqdef(f_{i+\h,\,j}g_{i-\h,\,j}+f_{i-\h,\,j}g_{i+\h,\,j})/2$ — which is the discrete analog of $f\frac{\partial fg}{\partial x}=\frac{1}{2}\frac{\partial f^2g}{\partial x}+\frac{f^2}{2}\frac{\partial g}{\partial x}$, and the fact that $\partial_x \overline{f}^x=\overline{\partial_x f}^x$, the integral of the four remaining terms can be rewritten as the integral of a term in divergence form plus the integral of the product of $\frac{1}{2}\bu\cdot\bu$ by the continuity equation. In other words, the integral of Eq. \eqref{eq:kinetic1derivative} over $V$ is $0$.

The $(\mathrm{Pres.})$ terms become (following [Morinishi et al. 2003]): \begin{equation} \frac{1}{r}\left(\overline{r\rho_0u_r\partial_r\frac{\delta p}{\rho_0}}^r+\overline{\rho_0u_{\theta}\partial_{\theta}\frac{\delta p}{\rho_0}}^{\theta}\right) \end{equation} Using the identity $\partial_x(f\overline{g}^x)=\overline{f\partial_x g}^x+(\partial_x f)g$, which is the discrete analog of the product rule this can be rewritten as a term in divergence form, plus the product of the pressure by the continuity equation. In other words, the integral of this term is also $0$.

The $(\mathrm{Corio.})$ terms become: \begin{equation} \frac{1}{r}\left(\overline{2ru_r\omega\overline{\rho_0 \overline{u_{\theta}}^{\theta}}^r}^r -\overline{2u_{\theta}\omega\overline{\rho_0 \overline{r u_r}^r}^{\theta}}^{\theta}\right) \end{equation} Using again the identity $\sum_V \overline{f\overline{g}^x}^x=\sum_V \overline{f}^xg$, we see that the integral of the two terms cancel each other. In other words, as in the continuous case, the Coriolis force conserve the kinetic energy.

Finally, the $(\mathrm{Buo.})$ term becomes: \begin{equation} \frac{1}{r}\overline{\rho_0r^2u_r\omega^2\frac{\overline{\delta T}^r}{T_0}}^r \end{equation}

In summary, the time derivative of the kinetic energy is equal to $\Delta r \Delta \theta\frac{\omega^2}{T_0}\sum_V \overline{\rho_0r^2u_r\overline{\delta T}^r}^r$. We now want to compute the time derivative of the potential energy, and show that it cancels with this term. For this we multiply the internal energy equation by $\frac{r^2\omega^2}{2c_vT_0}$, yielding the derivative: \begin{equation} -\frac{r}{2}\partial_r\left(r\rho_0 u_r \omega^2\frac{\overline{\delta T}^r}{T_0}\right) -\frac{1}{2r}\partial_{\theta}\left(\rho_0 r^2 u_{\theta} \omega^2\frac{\overline{\delta T}^{\theta}}{T_0}\right) \end{equation} The integral of the second term over $V$ is $0$ since this is a term in divergence form. The integral of the first term gives $-\Delta r \Delta \theta\frac{\omega^2}{T_0}\frac{1}{2}\sum_Vr^2\partial_r\left(r\rho_0 u_r \overline{\delta T}^r\right)$ which, using the identity $\frac{1}{2}\sum_V x^2\partial_x f=\sum_V\overline{x f}^x$ — valid for a uniform grid with $f$ periodic and/or equal to $0$ on $\partial V$, is exactly the opposite of the time derivative of the kinetic energy. This shows that the time derivative of the total energy is $0$, as we wanted.

Note that in the above proof we computed continuous time derivatives, using rules such as $u_r\frac{\partial u_r}{\partial t}=\frac{1}{2}\frac{\partial u_r u_r}{\partial t}$. However, this equation is not verified when using discrete time derivatives: $u_r^t\frac{u_r^{t'}-u_r^t}{\delta t}\ne\frac{1}{2}\frac{u_r^{t'} u_r^{t'}-u_r^tu_r^t}{\delta t}$. Thus, to conserve energy in practice, i.e. to minimize the errors due to the time integration, it is necessary to use a high order time integration scheme. This is why we use the fourth order Runge-Kutta method.

Finally, it should be noted that, unfortunately, the discretized equations do not conserve the discretized angular momentum $r\rho_0\overline{u_{\theta}}^{\theta}$, unlike the continuous ones. Indeed, if we multiply Eq. \eqref{eq:momentumt5} by $r$ and rewrite the $\partial_r\overline{r\rho_0 u_r}^{\theta}\overline{u_{\theta}}^r$ term with the discretized product rule, we get a term in conservation form plus $-\frac{1}{r}\overline{\overline{r\rho_0u_r}^{\theta}\overline{u_{\theta}}^r}^r$ which does not cancel with the third convection term $\overline{\overline{\rho_0 u_{\theta}}^{\theta}\overline{u_r}^r}^{\theta}$ (unlike in the continuous case). We can ensure angular momentum conservation by using the former instead of the latter in Eq. \eqref{eq:momentumt5}, but then we loose the energy conservation property.

6 Results

Our implementation of the above numerical method, presented in the Appendix, gives the following results.

6.1 Atmosphere heating

Here we want to simulate what happens when the light sources in Rama are turned on, if the atmosphere is initially at hydrostatic equilibrium. Due to radiative transfer, the atmosphere is heated at different rates in different regions, which generates buoyancy forces and convection cells. The numerical results for this case are the following.

First, if we use the heat source term $\delta J$ obtained in the Coupled lighting and temperature profiles section, which has vey narrow spikes near the light sources, we quickly get a lot of turbulence. This is most probably due to an insufficient grid resolution for the very narrow spikes. To avoid this problem, we use a $\delta J$ term without the spikes, which would correspond to the case of linear light sources inside light bulbs (note that with or without spikes, the integral of this term over the whole atmosphere is 0; because of this we expect the internal energy to be conserved, despite the presence of this heat source term). The result in this case is shown below, for the first 6 hours (=120 rotations) after the lights are turned on:

t = 10m

where:

the left image shows the streamlines of the velocity field $\bu$. The red streamlines rotate in the same direction as Rama (i.e. counter clock-wise in these figures), while the blue ones rotate in the opposite direction. The background shows the magnitude of the velocity, from 0 (black) to the maximum velocity (white).
the middle image shows the isobars of the dynamic pressure field $\delta p$. Their color varies from blue for the lowest pressure to purple for the average pressure and to red for the highest pressure.
the right image shows the isotherms of the temperature field $\delta T$. Their color varies from blue for the lowest temperature to purple for the average temperature and to red for the highest temperature.

We can see from these images that:

the streamlines and the isobars are very similar, meaning that the velocity is nearly perpendicular to the dynamic pressure gradient. This is called a geostrophic flow.
the high pressure areas, or anticyclones, are to the right of the flow, and rotate in the opposite direction as Rama, like anticyclones in the Northern Earth hemisphere. Conversely, the low pressure areas, or cyclones, rotate in the same direction as Rama.
the high (resp. low) pressure areas correspond more or less to the high (resp. low) temperature areas.
at $t=3h$ the convection cells have become mostly horizontal, and stay like this afterwards. This is due to the very strong Coriolis force (this effect is greatly reduced in numerical tests using a smaller angular speed $\omega$ for Rama).
at $t\approx 4h$ the temperature gradient tends to $\infty$ at one point on the surface, and after that turbulence starts to appear near this point in all fields. This effect also happens, at the same time, with a doubled grid resolution ($128\times 256$ instead of $64\times 128$). It also happens, a bit later, with a smaller angular velocity $\omega$. This suggests that this effect is probably "real", as opposed to a simulation artefact, but this is not sure.

The plot of the kinetic, internal, potential and total energy, as well as of the angular momentum, during these 6 simulated hours, is shown in Fig. 2. This plot shows that internal energy is conserved, as expected, while the total energy is not (which is expected too since there is an external heat source). In fact the total energy increases almost proportionally to $t^4$. Finally, the angular momentum is not conserved, while it should be: this is a limitation of the numerical method, as discussed above.

Figure 2: The kinetic, internal, potential and total energy per meter (in $J.m^{-1}$), and the angular momentum per meter (in $N.s$), as a function of $t$ in hours.

We can also plot the maximum velocity, the maximum temperature variation and the maximum pressure variation as a function of time. The result is shown in Fig. 3. This shows that the temperature and pressure variations are very small compared to their hydrostatic values, as they are supposed to be for the anelastic Boussinesq approximation to be valid. The velocity is also quite small, about $1m.s^{-1}$. Finally, this graph also confirms that the simulation becomes somewhat chaotic after $t\approx 4h$, but does not tell whether this is a simulation artefact or a real effect.

Figure 3: The maximum velocity (in $m.s^{-1}$), temperature variation (in $K$) and pressure variation (in $Pa$), as a function of $t$ in hours.

If we let the simulation run for 18 simulated hours, these graphs become the ones shown in Fig. 4. They show a stabilization of the maximum velocity at about $3 m.s^{-1}$, of the maximum temperature variation at about $0.1 K$, and of the maximum pressure at about $45 Pa$. However, the turbulence which appeared at $t=4h$ has then spread to the whole atmosphere, as shown in Fig. 5, so these results can probably not be trusted (the velocity spike at $t\approx 8h$ is certainly an artefact).

Figure 4: The maximum velocity (in $m.s^{-1}$), temperature variation (in $K$) and pressure variation (in $Pa$), as a function of $t$ in hours.

Figure 5: The simulated atmosphere state a $t=18h$.

6.2 Atmosphere heating, 2x version

Doing the same simulation as above, using a grid resolution of $128\times 256$ instead of $64\times 128$ gives the results below. These results are very similar to the above ones, except for the turbulence at $t>4h$ (showing that, even if the emergence of this turbulence is a real phenomenon, it can't be precisely simulated at this resolution).

t = 10m

The quantitative results for the total energy are also very similar, as shown in Fig. 6:

Figure 6: The kinetic, internal, potential and total energy per meter (in $J.m^{-1}$), and the angular momentum per meter (in $N.s$), as a function of $t$ in hours.

6.3 Energy conservation

Here we suppose that the atmosphere has reached the radiative equilibrium computed in the Coupled lighting and temperature profiles section, and that we then turn the lights off. In other words we use the initial conditions $\bu=0$ and $\delta T=\pi\delta J/4\sigma T_0^3$, i.e. we use the temperature such that the right hand side of the internal energy equation is 0. We then force this right hand side to 0 at any time $t>0$. Note that this case is not very realistic, as it would require to prevent the atmosphere from moving while the light sources are on. However, this is an interesting case to check our numerical method (energy should be conserved since there is no external heat source). The result is the following (for 1 simulated hour):

t = 1m 40s

This result is qualitatively very similar to the "atmosphere heating" case (geostrophic flow, convection cells becoming nearly horizontal because of the Coriolis force, etc), including the apparition of an infinite temperature gradient at the same point, leading short after to turbulence in all fields. This time, however, total energy is strictly conserved, as expected and as shown in Fig. 7.

Figure 7: The kinetic, internal, potential and total energy per meter (in $J.m^{-1}$), and the angular momentum per meter (in $N.s$), as a function of $t$ in hours.

The maximum velocity, maximum temperature variation and maximum pressure variation are shown in Fig. 8, from which we can see that the maximum velocity "converges" towards about $5m.s^{-1}$ (for $t\rightarrow\infty$ the atmosphere should normally go back to its hydrostatic equilibrium in this case, meaning a 0 velocity).

Figure 8: The maximum velocity (in $m.s^{-1}$), temperature variation (in $K$) and pressure variation (in $Pa$), as a function of $t$ in hours.

6.4 Angular momentum

In order to check whether the previous results are affected by the fact that our numerical method does not conserve angular momentum (unlike the continuous equations), we also tested the modified $(Conv.)_{\theta}$ term providing angular momentum conservation (at the price of a non conservation of energy). We did this on the "atmosphere heating" case, and the result is the following:

t = 10m

In other words the result is almost the same whether we choose to conserve energy or angular momentum (whose conservation in this case can be checked in Fig. 9).

Figure 9: The kinetic, internal, potential and total energy per meter (in $J.m^{-1}$), and the angular momentum per meter (in $N.s$), as a function of $t$ in hours.

7 Conclusion

Assuming that the atmosphere is an ideal fluid (i.e. no viscosity, no thermal conduction), an ideal gas, and that it is dry, assuming as well that the Boussinesq approximations hold ($p$, $\rho$ and $T$ very close to their hydrostatic equilibrium values), that forcing energy conservation in these approximated equations is valid, and finally that natural convection does not significantly change the radiative transfer properties, we obtain the following results.

During the first 4 hours after the light sources are turned on, and starting from the hydrostatic equilibrium, 3 pairs of convection cells rotating in opposite directions appear and become progressively almost horizontal (i.e. parallel to the surface) because of the Coriolis force. The flow is always quasi geostrophic. The maximum velocity reaches about 0.5 $m.s^{-1}$ (1.8 $km.h^{-1}$) and the temperature and pressure variations remain very small (about 0.05° and 20 $Pa$). After 4 hours the results are not conclusive: turbulence appears at one point on the surface, and quickly expands to the whole atmosphere, but this could be an artefact of the simulation.

In all simulations, using different initial conditions, grid resolutions, simulated durations, and conservation properties, the wind speed never goes above 5 $m.s^{-1}$ (18 $km.h^{-1}$), so this is probably a good estimate of the maximum wind speed in Rama. Likewise, in all these cases, the convection pattern is very similar, with 3 pairs of convection cells that quickly become almost horizontal (they are then replaced with turbulence but, again, this might be a simulation artefact).