Simulating Non-hydrostatic atmospheres on Planets (SNAP): formulation, validation and application to the Jovian atmosphere

TL;DR

This paper introduces a new non-hydrostatic atmospheric model built on Athena++, capable of simulating moist convection and cloud formation in planetary atmospheres, validated on Earth and Jovian conditions, and relevant to recent Juno observations.

Contribution

The paper presents a novel, unified atmospheric modeling framework that extends hydrodynamics to include vapors, clouds, and microphysics for planetary atmospheres, with advanced numerical methods.

Findings

Model accurately simulates Earth's atmospheric convection.

Jovian simulations reveal superadiabatic temperature gradients near water condensation levels.

Ammonia gas profiles show depletion consistent with observations.

Abstract

A new non-hydrostatic and cloud-resolving atmospheric model is developed for studying moist convection and cloud formation in planetary atmospheres. It is built on top of the Athena++ framework, utilizing its static/adaptive mesh-refinement, parallelization, curvilinear geometry, and dynamic task scheduling. We extend the original hydrodynamic solver to vapors, clouds, and precipitation. Microphysics is formulated generically so that it can be applied to both Earth and Jovian planets. We implemented the Low Mach number Approximate Riemann Solver (LMARS) for simulating low speed atmospheric flows in addition to the usual Roe and HLLC Riemann solvers. Coupled with a fifth-order Weighted Essentially Nonoscillatory (WENO) subgrid-reconstruction method, the sharpness of critical fields such as clouds is well-preserved, and no extra hyperviscosity or spatial filter is needed to stabilize the model. Unlike many atmospheric models, total energy is used as the prognostic variable of the thermodynamic equation. One significant advantage of using total energy as a prognostic variable is that the entropy production due to irreversible mixing process can be properly captured. The model is designed to provide a unified framework for exploring planetary atmospheres across various conditions, both terrestrial and Jovian. First, a series of standard numerical tests for Earth's atmosphere is carried out to demonstrate the performance and robustness of the new model. Second, simulation of an idealized Jovian atmosphere in radiative-convective equilibrium shows that 1) the temperature gradient is superadiabatic near the water condensation level because of the changing of the mean molecular weight, and 2) the mean profile of ammonia gas shows a depletion in the subcloud layer down to nearly 10 bars. Relevance to the recent Juno observations is discussed.

Tables

Table 1. Table 1: List of symbols

Notation	Meaning	Units
${\rho }_d$	density of dry air	$\mathrm{kg}/m^3$
${\rho }_i$	density of vapor $i$	$\mathrm{kg}/m^3$
${\rho }_{ij}$	density of the condensate $j$ associated with the vapor $i$	$\mathrm{kg}/m^3$
$\rho$	total density, $\rho ={\rho }_d+{\sum }_i{\rho }_i+{\sum }_{i,j}{\rho }_{ij}$	$\mathrm{kg}/m^3$
$q_d$	mass mixing ratio of dry air, $q_d=\frac{{\rho }_d}{\rho }$	1
$q_i$	mass mixing ratio (specific humidity) of vapor $i$, $q_i=\frac{{\rho }_i}{\rho }$	1
$q_{ij}$	mass mixing ratio of condensate $(i,j)$, $q_{ij}=\frac{{\rho }_{ij}}{\rho }$	1
$q_{ij}^{\ast }$	saturation mass mixing ratio of vapor $i$ over condensate $(i,j)$	1
$g$	gravitational acceleration	$m^2/s$
$x,y$	horizontal coordinates	m
$z$	vertical coordinate	m
$u,v$	horizontal velocities	$m/s$
$w$	vertical velocity	$m/s$
$w_{ij}^t$	sedimentation velocity of condensate $(i,j)$	$m/s$
$p_d$	partial pressure of dry air	$\mathrm{pa}$
$p_i$	partial pressure of vapor $i$	$\mathrm{pa}$
$p_{ij}^{\ast }$	saturation vapor pressure of vapor $i$ over condensate $(i,j)$	$\mathrm{pa}$
$p$	total pressure, $p=p_d+{\sum }_i{p}_i$	$\mathrm{pa}$
$e$	total specific energy (not including gravitational potential)	$J/\mathrm{kg}$
$𝐗$	primitive variables $𝐗={(\rho ,q_i,q_{ij},u,v,w,p)}^T$
$𝐘$	conserved variables $𝐘={({\rho }_d,{\rho }_i,{\rho }_{ij},\rho u,\rho v,\rho w,\rho e)}^T$
$c_{v,d}$	specific heat capacity of dry air at constant volume	$J/\left(\mathrm{kg}K\right)$
$c_{p,d}$	specific heat capacity of dry air at constant pressure	$J/\left(\mathrm{kg}K\right)$
$c_{v,i}$	specific heat capacity of vapor $i$ at constant volume	$J/\left(\mathrm{kg}K\right)$
$c_{p,i}$	specific heat capacity of vapor $i$ at constant pressure	$J/\left(\mathrm{kg}K\right)$
$c_{ij}$	specific heat capacity of condensate $(i,j)$	$J/\left(\mathrm{kg}K\right)$
$\Delta c_{ij}$	difference of specific heat capacity between condensate and vapor, $\Delta c_{ij}=c_{ij}-c_{p,i}$	$J/\left(\mathrm{kg}K\right)$
$L_{ij}$	latent heat of forming condensate $(i,j)$	$J/\mathrm{kg}$
$L_{ij}^r$	latent heat of forming condensate $(i,j)$ at reference temperature $T^r$	$J/\mathrm{kg}$
$T^r$	reference temperature	$K$
$p_{ij}^r$	saturation vapor pressure of vapor $i$ over condensate $(i,j)$ at reference temperature $T^r$	$\mathrm{pa}$
${\mu }_i$	chemical potential of vapor $i$	$J/\mathrm{kg}$
${\mu }_{ij}$	chemical potential of condensate $(i,j)$	$J/\mathrm{kg}$
$m_d$	molecular weight of dry air	$\mathrm{kg}/\mathrm{mol}$
$m_i$	molecular weight of vapor $i$	$\mathrm{kg}/\mathrm{mol}$
$\hat{R}$	ideal gas constant, $\hat{R}=8.3144598$ J/(mol K)	$J/\left(\mathrm{mol}K\right)$
$R_d$	ideal gas constant of dry air, $R_d=\hat{R}/m_d$	$J/\left(\mathrm{kg}K\right)$
$R_i$	ideal gas constant of vapor $i$, $R_i=\hat{R}/m_i$	$J/\left(\mathrm{kg}K\right)$
${ϵ}_i$	ratio of molecular weights between vapor and dry air, ${ϵ}_i=m_i/m_d$	1
${\sigma }_{v,i}$	ratio of specific heats between vapor and dry air, ${\sigma }_{v,i}=c_{v,i}/c_{v,d}$	1
${\sigma }_{p,i}$	ratio of specific heats between vapor and dry air, ${\sigma }_{p,i}=c_{p,i}/c_{p,d}$	1
${\sigma }_{v,ij}$	ratio of specific heats between condensate and dry air, ${\sigma }_{v,ij}=c_{ij}/c_{v,d}$	1
${\sigma }_{p,ij}$	ratio of specific heats between condensate and dry air, ${\sigma }_{p,ij}=c_{ij}/c_{p,d}$	1
$\gamma$	adiabatic index	1
${\gamma }_d$	adiabatic index of dry air, ${\gamma }_d=c_{p,d}/c_{v,d}$	1
$T$	temperature	$K$
$T_v$	virtual temperature, $T_v=T\left(1+{\sum }_i{q}_i\left(1/{ϵ}_i-1\right)-{\sum }_{i,j}{q}_{ij}\right)$	$K$
$\theta$	potential temperature	$K$
${\theta }_v$	virtual potential temperature, ${\theta }_v=\theta \left(1+{\sum }_i{q}_i\left(1/{ϵ}_i-1\right)-{\sum }_{i,j}{q}_{ij}\right)$	$K$

Equations

$$\frac{\partial {\rho }_d}{\partial t}+\frac{\partial ({\rho }_{d}u)}{\partial x}+\frac{\partial ({\rho }_{d}v)}{\partial y}+\frac{\partial ({\rho }_{d}w)}{\partial z}$$

$$\frac{\partial {\rho }_i}{\partial t}+\frac{\partial ({\rho }_{i}u)}{\partial x}+\frac{\partial ({\rho }_{i}v)}{\partial y}+\frac{\partial ({\rho }_{i}w)}{\partial z}$$

$$\frac{\partial {\rho }_{ij}}{\partial t}+\frac{\partial ({\rho }_{ij}u)}{\partial x}+\frac{\partial ({\rho }_{ij}v)}{\partial y}+\frac{\partial ({\rho }_{ij}w)}{\partial z}$$

$$\frac{\partial (\rho u)}{\partial t}+\frac{\partial (u\rho u+p)}{\partial x}+\frac{\partial (v\rho u)}{\partial y}+\frac{\partial (w\rho u)}{\partial z}$$

$$\frac{\partial (\rho v)}{\partial t}+\frac{\partial (u\rho v)}{\partial x}+\frac{\partial (v\rho v+p)}{\partial y}+\frac{\partial (w\rho v)}{\partial z}$$

$$\frac{\partial (\rho w)}{\partial t}+\frac{\partial (u\rho w)}{\partial x}+\frac{\partial (v\rho w)}{\partial y}+\frac{\partial (w\rho w+p)}{\partial z}$$

$$\rho ={\rho }_d+\sum _i{\rho }_i+\sum _{i,j}{\rho }_{ij}$$

$$\rho e={\rho }_d{c}_{v,d}T+\sum _i{\rho }_i{c}_{v,i}T+\sum _{i,j}{\rho }_{ij}{c}_{ij}T+\frac{1}{2}\rho (u^2+v^2+w^2)+\sum _i{\mu }_i{\rho }_i+\sum _{i,j}{\mu }_{ij}{\rho }_{ij},$$

$$L_{ij}(T)=L_{ij}^r-\Delta c_{ij}(T-T^r),$$

$$L_{ij}(T)=c_{p,i}T+{\mu }_i-(c_{ij}T+{\mu }_{ij})$$

$${\mu }_{ij}={\mu }_i-(L_{ij}(T)+\Delta c_{ij}T)={\mu }_i-(L_{ij}^r+\Delta c_{ij}{T}^r)$$

$$({\rho }_d{c}_{v,d}+\sum _i{\rho }_i{c}_{v,i}+\sum _{i,j}{\rho }_{ij}{c}_{ij})T=({\rho }_d{c}_{v,d}+\sum _i{\rho }_i{c}_{v,i}+\sum _{i,j}{\rho }_{ij}{c}_{ij}-\Delta \rho c_{v,k}+\Delta \rho c_{kl})(T+\Delta T)+({\mu }_{kl}-{\mu }_k)\Delta \rho ,$$

$$\Delta T=\frac{\Delta \rho (L_{kl}-R_{k}T)}{{\rho }_d{c}_{v,d}+{\sum }_i{\rho }_i{c}_{v,i}+{\sum }_{i,j}{\rho }_{ij}{c}_{ij}-\Delta \rho c_{v,k}+\Delta \rho c_{kl}}$$

$$\frac{\partial (\rho e)}{\partial t}+\frac{\partial [u(\rho e+p)]}{\partial x}+\frac{\partial [v(\rho e+p)]}{\partial y}+\frac{\partial [w(\rho e+p)]}{\partial z}=-\rho gw-\sum _{i,j}{\rho }_{ij}g{w}_{ij}^t-\sum _{i,j}\frac{\partial (w_{ij}^t{\rho }_{ij}{e}_{ij})}{\partial z},$$

$$p$$

$$=\rho R_{d}T(q_d+\sum _i{q}_i/{ϵ}_i)$$

$$\displaystyle=\rho R_{d}T\big{[}1+\sum_{i}q_{i}(1/\epsilon_{i}-1)-\sum_{i,j}q_{ij}\big{]}$$

$$=\rho R_d{T}_v,$$

$$T_{v}=T\big{[}1+\sum_{i}q_{i}(1/\epsilon_{i}-1)-\sum_{ij}q_{ij}\big{]}$$

$$\frac{\partial q_i}{\partial t}$$

$$\frac{\partial q_{ij}}{\partial t}$$

$$\frac{{\mathbf{X}}^{n+1,m}-{\mathbf{X}}^n}{\Delta t}=\mathbf{\Phi }({\mathbf{X}}^{n+1,m})\approx \mathbf{\Phi }({\mathbf{X}}^{n+1,m-1})+\frac{\partial \mathbf{\Phi }}{\partial \mathbf{X}}({\mathbf{X}}^{n+1,m}-{\mathbf{X}}^{n+1,m-1}),$$

$${\mathbf{X}}^{n+1,0}={\mathbf{X}}^n,$$

$$\mathbf{X}^{n+1,m}=\mathbf{X}^{n+1,m-1}+(\frac{\mathbf{I}}{\Delta t}-\frac{\partial\mathbf{\Phi}}{\partial\mathbf{X}})^{-1}\big{[}\mathbf{\Phi}(\mathbf{X}^{n+1,m-1})-\frac{\mathbf{X}^{n+1,m-1}-\mathbf{X}^{n}}{\Delta t}\big{]}$$

$${\mathbf{X}}^{n+1,m}$$

$$\frac{{\mathbf{X}}^{n+1}-{\mathbf{X}}^n}{\Delta t}$$

$$\frac{\partial q_i}{\partial t}$$

$$\frac{\partial q_{i1}}{\partial t}$$

$$\frac{\partial q_{i2}}{\partial t}$$

$$w_{i1}^t$$

$$w_{i2}^t$$

$$\frac{\partial \mathbf{\Phi }}{\partial \mathbf{X}}=\left(\begin{array}{ccc}-k_1& 0& 0\\ k1& -k_2-k_3{q}_{i2}& -k_3{q}_{i1}\\ 0& k_2+k_3{q}_{i2}& k_3{q}_{i1}\end{array}\right)$$

$$\frac{\partial \mathbf{\Phi }}{\partial \mathbf{X}}=\left(\begin{array}{ccc}-k_4{q}_{i2}& 0& 0\\ 0& -k_2-k_3{q}_{i2}& -k_3{q}_{i1}\\ k_4{q}_{i2}& k_2+k_3{q}_{i2}& k_3{q}_{i1}\end{array}\right)$$

$$Q_{i,j,k}^{n+1}=Q_{i,j,k}^n$$

$$-\frac{\Delta t}{\Delta y}(G_{i,j+1/2,k}^n-G_{i,j-1/2,k}^n)$$

$$-\frac{\Delta t}{\Delta z}(H_{i,j,k+1/2}^n-H_{i,j,k-1/2}^n)+\Delta t{\mathcal{F}}_{i,j,k}^n,$$

$$F_{i-1/2}=\mathcal{R}(Q_{i-1/2}^{-},Q_{i-1/2}^{+})$$

$$Q_{i-1/2}^{-}$$

Full text

cloud, atmosphere, Jupiter, cloud-resolving, Riemann solver, ammonia

††journal: ApJS

1 Introduction

Atmospheric science for Earth and planets has come to a point where a unified model is needed to study diverse climates on numerous exotic worlds. Earth’s atmosphere has been observed and studied extensively. Yet, from a broader perspective, the possible range of Earth’s climate only occupies a niche of the parameter space of all planetary atmospheres that have been observed in the universe. Mars, Pluto, and Triton have atmospheres made up mostly of condensible gases. Condensation of CO2 or N2 on to the planetary surface removes a significant portion of the total mass of the atmosphere, resulting in a global circulation driven by sublimation and condensation (Pierrehumbert & Ding, 2016; Ding & Pierrehumbert, 2018). The atmosphere of giant planets harbors multiple condensible species whose molecular weights are much larger than the ambient hydrogen atmosphere. The loading of heavy molecules stratifies the atmosphere and inhibits moist convection (Sugiyama et al., 2014; Li & Ingersoll, 2015; Leconte et al., 2017; Friedson & Gonzales, 2017). Beyond the solar system, nearly 3000 exoplanets have been discovered at the time of writing111Statistics are gathered from www.exoplanets.org. Many of them have no counterparts in our solar system (hot Jupiters, super Earth, lava planets, etc.). It is, therefore, challenging to picture their atmospheric circulation pattern using the knowledge from studying the solar system.

Global circulation models (GCM) are powerful tools for simulating atmospheric flows on planets. Limited by computational resources, the horizontal resolution of a GCM is usually around 100 kilometers for Earth (e.g. AM4 model: Zhao et al. 2018; Lin04) and 500 kilometers for Jovian planets (e.g. Jupiter GCM: Schneider & Liu 2009). Models with such a low resolution cannot resolve small-scale processes like moist convection or cloud formation. Therefore, they must employ parameterization schemes that approximate the effect of small-scale processes. The parameterization schemes for Earth are usually motivated by observations or by regional numerical models that resolve these processes (see Arakawa & Jung 2011 for a review). For planets other than Earth, the absence of granular observations precludes the possibility of using observations to constrain the parameterization scheme. Direct modeling of convection embedded in a large-scale environment becomes necessary to understand how the parameterization scheme should be for a planet GCM. Right now, many planet GCMs either borrow parameterization schemes designed for Earth (e.g. Schneider & Liu 2009), or neglect any parameterization of convection (e.g. Lian & Showman 2010). Such approaches may be valid for investigating problems like jet formation or global temperature distribution, but it is doubtful when one tries to use a low-resolution GCM to understand how tracers or clouds are distributed without an appropriate treatment of convection and diffusion, especially for atmospheres that do not resemble Earth’s.

The purpose of this paper is to introduce a newly developed three-dimensional non-hydrostatic model that simulates moist convection and cloud formation across various planetary conditions. Using a regional non-hydrostatic model to study moist convection is standard practice for Earth’s atmosphere, and many numerical models – named cloud-resolving models (CRM) or large eddy simulations (LES) – have been developed accordingly (Pielke et al., 1992; Khairoutdinov & Randall, 2003; Bryan & Fritsch, 2002; Pressel et al., 2015). Similar non-hydrostatic models have been developed for Mars (Rafkin et al., 2001), Jupiter (Sugiyama et al., 2014; Hueso & Sánchez-Lavega, 2001) and Saturn (Hueso & Sánchez-Lavega, 2004). But they lack a well-documented comparison against an accepted nonlinear solution, which seems to be the first step of developing a robust numerical model. Since there are no standard tests of non-hydrostatic models for planets other than Earth, we will primarily perform simulations against standard tests for Earth-like conditions. Then, we perform an idealized simulation of the Jovian atmosphere in radiative-convective equilibrium (RCE) to study moist convection in hydrogen atmospheres.

The structure of this work is organized as follows. Section 2 is devoted to the equation of motion, which includes dry air, vapors, clouds, and precipitation. Section 3 discusses the basic microphysics scheme employed in the model. Section 4 describes the numerical schemes that are designed for simulating atmospheric flows. Section 5 provides a series of benchmark tests against known solutions in the literature. Section 6 performs an idealized RCE simulation of the Jovian atmosphere. Section 7 concludes and outlines the future applications. To make the notations consistent and clear, major symbols used in this paper are summarized in Appendix A for reference.

2 Equation of motion

We solve the most generic form of the equation of motion, the fully compressible Euler equations, to avoid unnecessary assumptions that are only valid for a particular type of planet. Though Euler equations can be written in many mathematically equivalent forms, they are different when discretized and implemented in a numerical model. Bryan & Fritsch (2002) has compared five different ways of implementing the Euler equations, using potential temperature, equivalent potential temperature, and liquid water potential temperature. They concluded that the form of governing equations used in a numerical model has a profound effect on the simulation of a warm rising bubble. Satoh (2002) pointed out that in a model with diabatic forcing, the change in the domain integral of total energy is generally different from energy gain/loss from the boundaries if the potential temperature is used as the prognostic variable, which causes an error of the energy budget that accumulates with time. Furthermore, the expressions of potential temperature and equivalent potential temperature become complicated when the heat capacity of the atmosphere varies with temperature or when multiple condensing species exist in the atmosphere. Given the limitations of using potential temperature, we use total energy as the prognostic variable in implementing the Euler equations. We build our hydrodynamic solver on top of the Athena++ framework, an efficient, scalable and well-tested astrophysical code (White et al., 2016; Stone et al., 2008; Stone & Gardiner, 2009). Using the Athena++ framework as the software infrastructure not only saves us a tremendous amount of the development work but also enables the atmospheric part to utilize all the power of the Athena++ infrastructure such as the static/adaptive mesh-refinement, curvilinear geometry, and dynamic task scheduling.

To simulate a heterogeneous atmosphere under the gravitational field, we extend the original Euler equations to vapors, clouds, and precipitation. The continuity equation for each homogeneous component should be solved individually, leading to a family of continuity equations. A heterogeneous air parcel comprises dry air ($\rho_{d}$), vapors ($\rho_{i}$), and condensates ($\rho_{ij}$), where $i$ is the index for the vapor and $j$ is the index for the phase of the condensate. A condensate is a general description of any condensed material, which may represent a liquid cloud, ice cloud, graupel, rain, snow, etc. Two indices $(i,j)$ are needed to identify a condensate because of the presence of multiple condensing species. Several simplifications to the true physical processes are made to limit the complexity of this model. We assume that the gaseous components share the same mean motion velocities ($u,v,w$) and the condensates have additional terminal velocities ($w^{t}_{ij}$) with respect to the mean motion of the air parcel. We treat each component in the atmosphere equally and solve the evolution of its density tendency separately as shown in equations (1-3).

[TABLE]

Comparing to equations (1-2), the equation for condensates, equation (3), has an additional flux term $w^{t}_{ij}\rho_{ij}$ in the vertical direction due to the flux of sedimentation. Sedimentation of cloud droplets is found to be important for simulating stratocumulus cloud on Earth (Ackerman et al., 2004; Bretherton et al., 2007) and probably for simulating the folded filament cloud on Jupiter as well. The value of terminal velocity $w^{t}_{ij}$ should be provided externally by a microphysics package specific to the application. Changes to the continuity equations due to thermodynamics or microphysics will be discussed separately in Section 3.

Precipitation and sedimentation not only transfer mass but transfer momentum and energy as well. The momentum equations now read:

[TABLE]

where the total density $\rho$ is defined as:

[TABLE]

Note that the total density includes the contribution from all gases and condensates regardless of whether sedimentation occurs or not. If the condensate is aloft in the air, its mass is part of the total mass of the air parcel. Otherwise, the frictional force during sedimentation equals the gravity of the condensate. Then, the term $\sum_{i,j}\rho_{ij}g$ can be interpreted as the drag force acted upon the air parcel by sedimentation.

The total energy of a heterogeneous air parcel is defined as:

[TABLE]

where the first three terms represent the internal energies, the fourth term is the kinetic energy, and the last two terms represent the chemical potentials. Because the chemical potential is defined only to within an additive constant, it is possible to set the chemical potentials of all gases ($\mu_{i}$) to zero so long as the gases are chemically independent. Otherwise, different chemically potentials should be assigned to reflect the differences in their internal energies. For example, in a hydrogen atmosphere, hydrogen in the ortho state (odd rotational quantum numbers, $J$) has higher internal energy than hydrogen in the para state (even $J$). The difference increases with decreasing temperature. Brought from the deeps, the ortho hydrogen slowly converts to the para hydrogen at a shallower depth on a time scale of about $10^{8}$ sec (Conrath & Gierasch, 1984). The latent heat release associated with the conversion is thought to be important for the energy transport in hydrogen atmospheres (Conrath et al., 1998). If such effect were to be calculated, ortho hydrogen and para hydrogen should be treated as different species and with different chemical potentials.

The chemical potential of the condensate ($\mu_{ij}$) is related to the chemical potential of the corresponding gas ($\mu_{i}$) and the latent heat ($L_{ij}$). Using the Kirchhoff’s equation222There are a lot of Kirchhoff’s equations. What we mean is: $dL_{ij}(T)/dT=c_{p,i}-c_{ij}$, the latent heat as a function of temperature is approximately (Emanuel, 1994):

[TABLE]

where $L^{r}_{ij}$ is latent heat of condensate $(i,j)$ at reference temperature $T^{r}$ and $\Delta c_{ij}=c_{ij}-c_{p,i}$ is the difference between the specific heat capacities. Latent heat is defined as enthalpy difference between the vapor and the condensate:

[TABLE]

Therefore, the chemical potential of condensate $(i,j)$ is:

[TABLE]

Assuming that gas $k$ undergoes phase transition, in which the mass transfer is $\Delta\rho$ and the condensate is $(k,l)$, the equation of the conservation of energy is:

[TABLE]

where $\Delta T$ is the temperature change. Using equations (9), (11), and (12) to solve for $\Delta T$ yields:

[TABLE]

The numerator is the latent heat release, and the denominator is the heat capacity of the air parcel after condensation. Note that the latent heat release at constant volume is different from the latent heat release at constant pressure, which is typically defined and used in a large-scale GCM.

The equation of total energy for atmospheric flows has forcing terms from gravity, friction, and sedimentation, which is:

[TABLE]

where $e_{ij}=c_{ij}T+\frac{1}{2}(u^{2}+v^{2}+w^{2})+\mu_{ij}$ is the specific total energy of condensate $(i,j)$. The first term on the right-hand-side is the forcing from gravity, the second term is the frictional heating from sedimentation, and the third term is the divergence of the energy flux due to sedimentation.

Finally, the Euler equations are closed by the equation of state, which is simply $p=\rho R_{d}T$ for a homogeneous ideal gas. Yet, in the presence of gases and clouds, the equation of state is obtained by summing the partial pressure of all gaseous components. Let $\epsilon_{i}$ be the molecular weight ratio of the vapor to the dry air. The total pressure is:

[TABLE]

where the virtual temperature $T_{v}$ is define as:

[TABLE]

In an atmosphere laden with vapors and clouds, it is the virtual temperature rather than the actual temperature that determines the density of an air parcel. Thus, a moist air parcel loaded with water vapor in Jupiter’s atmosphere $(\epsilon_{i}>1)$ is denser than a dry air parcel at the same temperature and pressure.

3 Microphysics and hydrological cycle

Simulating the formation and dissipation of clouds is indispensable in a non-hydrostatic and cloud-resolving model. In the presence of abundant condensation nuclei, cloud forms when the partial pressure of the vapor exceeds its saturation vapor pressure. Then cloud particles undergo coagulation and coalescence to grow in size. When they become large enough to have appreciable falling velocities, precipitation occurs. The falling rain droplets or snowflakes re-evaporate into unsaturated warm air below the cloud, adding mass to the air parcel and reducing its temperature. The whole subject of microphysical processes regarding clouds is very complicated and is still not very well understood especially for how cloud forms on other planets. Moreover, most microphysical packages are developed specifically for water clouds on Earth (e.g. Kessler 1969; Lin et al. 1983; Seifert & Beheng 2001). How to adapt the existing parameterizations developed for Earth to planets is a grand topic open for future study (some advances can be found at Rossow 1978; Carlson et al. 1988). It should be made clear that the subject of this section is not to compare and contrast different microphysics schemes, nor to design a new one, but rather, to describe a unified framework of implementing microphysics in a dynamic model, that works for both terrestrial planets and giant planets. It may seem unnecessary at this point for the dynamic part (core) of a numerical model is usually separated from other physical packages in a traditional GCM. However, since we are using total energy (including chemical potentials) as the prognostic variable of the thermodynamic equation, the microphysics is so intimately related to dynamics that should be implemented in a way that is congruent with the formulation of the dynamics.

The fact that dynamic model stays with density of the dry air ($\rho_{d}$) and total energy ($\rho e$) as prognostic variables eases the implementation of the microphysics because those two fields are also conserved during phase change. What remains is to calculate the partition of the vapors ($\rho_{i}$) and the corresponding condensates ($\rho_{ij}$) subject to the specific formulation of the microphysics. Assuming that the reaction rates are described by abstract functions:

[TABLE]

where $q_{i}=\rho_{i}/\rho$ and $q_{ij}=\rho_{ij}/\rho$ are mass mixing ratios. The reaction rates are written as functions of primitive variables $\mathbf{X}=(\rho,q_{i},q_{ij},u,v,w,p)^{T}$ since microphysics or chemical reactions usually works with mixing ratios rather than densities. The conversion between the conserved variables $\mathbf{Y}=(\rho_{d},\rho_{i},\rho_{ij},\rho u,\rho v,\rho w,\rho e)^{T}$ and the primitive variables $\mathbf{X}$ is detailed in Appendix B.

The evolution equations (17) and (18) can be integrated using a standard ODE integrator, such as the first-order backward Euler integration method:

[TABLE]

where $\frac{\partial\mathbf{\Phi}}{\partial\mathbf{X}}$ is the Jacobian matrix and $m$ is the iteration number. Iteration is needed because a linear expansion is used to approximate $\mathbf{\Phi}(\mathbf{X}^{n+1})$. The iteration starts with:

[TABLE]

and progresses forward with:

[TABLE]

The iteration stops when $||\mathbf{X}^{n+1,m}-\mathbf{X}^{n+1,m-1}||$ is smaller than an arbitrary small number. Then, at the last iteration:

[TABLE]

which completes the backward Euler integration method.

Equation (21) requires both the reaction rate $\mathbf{\Phi}(\mathbf{X})$ and the Jacobian matrix of the reaction rate $\frac{\partial\mathbf{\Phi}}{\partial\mathbf{X}}$ to be evaluated at each time step. Although equation (21) is only first-order accurate in time, it is unconditionally monotone and is especially good at integrating stiff equations in which the chemical timescale (such as condensation) is much shorter than the dynamic time step. A second-order accurate and L-stable integration method like TR-BDF2 (Hosea & Shampine, 1996) can also be used to integrate equations (17) and (18). However, given the uncertainties in the reaction rate coefficients, using the backward Euler scheme is probably enough for simulating planetary atmospheres. The above procedure costs more in computational time than the simple Euler forward method, but it is more robust regarding numerical stability and can be generically applied to any chemical systems.

Here we give an example of how the Kessler scheme (Kessler, 1969) is implemented according to the above framework. The Kessler scheme is the simplest cloud microphysics scheme that represents the cycling of the three phases of a condensible gas (vapor, cloud, and precipitation). Although it is considered outdated for modeling Earth’s atmosphere, it is still of great value when applied to planetary atmospheres due to its simplicity. In fact, many cloud-resolving simulations for the Jovian atmosphere, including Yair et al. (1995); Nakajima et al. (2000); Hueso & Sánchez-Lavega (2004); Sugiyama et al. (2014), started with microphysics schemes similar to Kessler (1969). The original paper of Kessler (1969) was specifically written for Earth’s atmosphere, with hardwired physical constants in the equations. We rewrite the formula and summarize the major equations in the following:

[TABLE]

where $k_{1},k_{2},k_{3}$, and $k_{4}$ are reaction rate coefficients concerning the rate of condensation, the rate of autoconvection (cloud becomes precipitation), the rate of accretion (precipitation grows bigger by collecting cloud particles), and the rate of evaporation, respectively. $k_{5}$ is a physical constant describing the terminal velocity of the precipitation. Sedimentation of cloud particles was not considered in Kessler (1969). $q_{i1}$ is the mass mixing ratio of the cloud, $q_{i2}$ is the mass mixing ratio of the precipitation, and $q^{*}_{i1}$ is the saturation vapor mixing ratio with respect to cloud $q_{i1}$. The symbol $x^{+}=max(x,0)$.

Coefficient $k_{1}$ controls how fast supersaturated vapor condenses to form clouds. If the condensation nuclei are abundant, condensation can occur almost instantaneously when the vapor is saturated. Otherwise, spontaneous nucleation requires a high degree of supersaturation. A larger value of $k_{1}$ produces less supersaturation. Coefficient $k_{2}$ controls how fast clouds start to precipitate, and thus the amount of cloud present in the atmosphere. Coefficient $k_{3}$ functions similarly as $k_{2}$, but the accretion rate $k_{3}q_{i1}q_{i2}$ is a quadratic term depending both on the amount of cloud and the amount of precipitation. Coefficient $k_{4}$ controls how fast precipitation evaporates. The evaporation rate is proportional to the amount of precipitation $q_{i2}$ and the amount of saturation deficit $(q^{*}_{i1}-q_{i})^{+}$. The values of these coefficients depend largely on particle number density, which can be either empirically prescribed or physically calculated via a more sophisticated two-moment microphysics scheme such as Seifert & Beheng (2006). Often, the exact values of these coefficients have little impact on the dynamics as long as they are chosen within a reasonable range, which can be informed by Earth’s conditions. Sugiyama et al. (2014) has varied the autoconversion rate ($k_{2}$) for two orders of magnitude, and they concluded that intermittent outburst of convective activity found in their simulation did not change greatly with $k_{2}$. In addition to the microphysical reactions outlined in equations (24) – (26), precipitation ($q_{i2}$) may boil in the atmosphere when the saturation vapor pressure equals the ambient atmospheric pressure. Boiling of rain droplets is not likely to occur in the terrestrial environment but must happen in Jovian atmospheres if the evaporation rate is low (large droplet size), and if the precipitation survives evaporation before reaching the boiling level – 30 bars for water and 4 bars for ammonia in the current Jupiter’s atmosphere (shown in Figure 1). The dynamic effect of boiling rain droplets is also an interesting topic worth future investigation.

The Jacobian matrix of the reaction rate has two expressions depending on whether the atmosphere is saturated or not. For a saturated atmosphere ($q_{i}>q^{*}_{i1}$):

[TABLE]

For an unsaturated atmosphere ($q_{i}<q^{*}_{i1}$):

[TABLE]

Since the Jacobian matrix is only $3\times 3$ in this case, using the backward Euler integration method as described by equation (21) is almost as fast as the simple Euler forward method with significant improvement in the numerical stability.

4 Numerical schemes

The hydrodynamic equations are discretized using the finite-volume method, in which all conserved quantities are placed in a collocated grid, representing the spatial average of the finite volume. Let $Q_{i,j,k}^{n}$ be any discretized conserved quantity at the cell $(i,j,k)$ and the $n$-th time. The simplest first-order accurate evolution equation in time is:

[TABLE]

where $F_{i-1/2,j,k}^{n}$ is the numerical flux between cells $(i-1,j,k)$ and $(i,j,k)$ at the $n$-th time step, and $\mathscr{F}_{i,j,k}$ is the body force in the right-hand-side of the hydrodynamic equations. The meanings of $G$ and $H$ are similar, but denote the fluxes in the $y-$ and $z-$ directions. Since the fluxes in all directions are calculated by the Riemann solver, we will focus on the $x$-direction and omit the $j,k$ indices in the subscript and the $n$ in the superscript for clarity.

The numerical flux $F_{i-1/2}$ is the solution of an initial value problem, the Riemann problem $\mathcal{R}$, consisting two states that are at the left ($Q_{i-1/2}^{-}$) and the right ($Q_{i-1/2}^{+}$) side of a jump discontinuity.

[TABLE]

The left and right states are interpolated from cell-averaged quantities, a process called subgrid-reconstruction. A $(2m-1)$-th order reconstruction scheme takes the form of:

[TABLE]

where $\mathcal{I}$ is the interpolation function and is biased toward one side of the discontinuity. Equations (31) - (34) summarized the standard Godunov method of solving hyperbolic equations (LeVeque, 2002), which were implemented in the hydrodynamic solver of the Athena++ code (Stone et al. (2008); Stone & Gardiner (2009)). However, because the atmospheric flows usually are at low Mach number and are stratified in the vertical direction, the solver designed for supersonic astrophysical flows should be modified according to the characteristics of the atmospheric flow. We elaborate our extension of the numerical methods of Riemann solver, subgrid-reconstruction, and time-stepping to atmospheric flows in the following subsections.

4.1 Riemann solver

The standard Riemann solvers in the Athena++ model are Roe solver (Roe, 1981) and HLLC solver (Toro, 2013). Applying the default solvers to atmospheric problems will induce spurious errors because the gravitational acceleration nearly balances the vertical pressure gradient. The small imbalance is what drives the vertical motion, but would be obscured by the interpolation error if the stratification in pressure is not taken into account in the reconstruction step or the Riemann solver (demonstrated in Figure 2 and Figure 3). Many remedies have been proposed to alleviate the numerical difficulty in the quasi-steady problems including subtracting a stationary state (Dedner et al., 2001), recognizing the pressure difference in the Riemann solver (Bale et al., 2003), incorporating the hydrostatic equilibrium in the pressure reconstruction (Käppeli & Mishra, 2016), etc. Here we use the recently developed Low Mach number Approximate Riemann Solver (LMARS) (Chen et al., 2013) to solve equation (32).

The LAMRS solver has several advantages over other approaches for atmospheric problems. First, as the Mach number approaches zero, the compressible Euler equations converge to the incompressible limit. Yet, the standard Riemann solver that captures shock, such as the Roe solver, failed to produce the convergence (Guillard & Viozat, 1999). The LAMRS solver is designed for low Mach number problems but can also handle the case in which the flow speed approaches the sound speed, bring an ideal solution to atmospheric flows. Second, the formulation of LMARS allows a seamless transition from hydrostatic models to non-hydrostatic models, a feature that will be useful for an adaptive mesh refined model. Third, the LAMRS scheme is easier to implement and faster to calculate than the Roe or HLLC solver. Its key results are:

[TABLE]

where $p^{\pm}_{i\pm 1/2}$ denotes the reconstructed pressure in the horizontal direction or the reconstructed perturbation pressure in the vertical direction. Knowing the pressure $p_{i-1/2}$ and the velocity $u_{i-1/2}$ at the cell boundary, it is easy to calculate the flux $F_{i-1/2}$. More details of the scheme are presented in Chen et al. (2013) and not repeated here.

4.2 subgrid-reconstruction

The Athena++ model has implemented two reconstruction methods, the piecewise-linear method (PLM) and the piecewise-parabolic method (PPM, Colella & Woodward (1984)). The PLM model is only second-order accurate in space, which would cause large phase error in the wave propagation when collocated variables are used (Chen et al., 2018). The PPM model has a higher order of accuracy but does not have an asymmetry between the left and right side of the cell boundary, which is required in the LMARS scheme. Because of the aforementioned shortcomings, we implement the fifth-order Weighted Essentially Non-Oscillatory (WENO) scheme (Shu, 1998) in lieu of the original ones. The fifth-order WENO scheme is computed as convex combinations of three third-order stencils:

[TABLE]

The weights $\alpha_{0},\alpha_{1}$ and $\alpha_{2}$ are nonnegative and add up to one; they control which stencil to use. For a smooth field, $\alpha_{0}=0.1,\alpha_{1}=0.6,\alpha_{2}=0.3$, and equation (37) reduces to the fifth-order polynomial interpolation scheme:

[TABLE]

If a discontinuity exists, the weights are adjusted such that the stencil is placed toward the direction away from the discontinuity. The formulas for the weights are provided in Jiang & Shu (1996) and Shu (1998). Following the convention of the original Athena++ model, the subgrid-reconstruction are performed using primitive variables ($\mathbf{X}$) instead of the conserved variables ($\mathbf{Y}$).

4.3 time-stepping

Equation (31) is only first-order accurate in time. To achieve higher order accuracy, a multi-stage time stepping method shall be used. A widely used third-order total variation diminishing (TVD) Runge-Kutta method is (Shu & Osher, 1988):

[TABLE]

where $F_{i-1/2}^{\dagger}$, $F_{i-1/2}^{\ddagger}$ and $\mathscr{F}_{i}^{\dagger}$, $\mathscr{F}_{i}^{\ddagger}$ are fluxes and body forces evaluated using intermediate states $Q_{i}^{\dagger}$ and $Q_{i}^{\ddagger}$. By default, the above third-order scheme is used. Using the original second-order schemes such as the van Leer integrator is also possible.

5 Benchmark tests

To illustrate that our model formulation and numerical scheme are suitable for atmospheric simulation, we perform three benchmark tests. The first one, proposed by Straka et al. (1993), simulates a dense sinking bubble and its subsequent propagation on the lower boundary. This case is targeted to validate how well the numerical model handles a solid lower boundary. The second case, proposed by Robert (1993), simulates a buoyant rising bubble and the fully developed Kelvin-Helmholtz instability. The forcing is weak in this case so that the turbulence field is sensitive to the numerical scheme. The third case, proposed by Bryan & Fritsch (2002), is similar to the second one but includes moisture and phase change. This case is intended for testing the interaction between the dynamics and the thermodynamics.

5.1 Straka density current

A standard test of a non-hydrostatic model is to simulate a sinking bubble and the propagation of the resulting density current. It has been widely used as a test case for many numerical models (e.g. Skamarock & Klemp (1993), Ooyama (2001), Pressel et al. (2015)) since Straka et al. (1993) formulated this classic problem. The background atmosphere has a constant potential temperature of 300 K referenced at 1 bar. A cold bubble is put aloft in the air with maximum temperature difference -15 K. Then the bubble sinks and propagates after it hits the ground.

The original setup specified by Straka et al. (1993) applied viscous dissipations to both the momentum equations and the energy equation to ensure the convergence of the solution. Here, we test the viscous solution for the comparison purpose as well as the nearly inviscid solution, in which no explicit diffusion is applied, to demonstrate the stability and the robustness of the numerical scheme. Figure 2(a) shows the viscous solution using LMARS solver, which is almost identical to the reference solution provided by Straka et al. (1993)’s Figure 1. However, if the HLLC solver was used, the potential temperature at the surface has unrealistically increased by over 1 K ahead of the density current (Figure 2b). The maximum potential temperature anomaly at 900s is 1.84 K using HLLC solver, while it is only 0.12 K using LMARS solver. In an adiabatic atmosphere, the potential temperature is very close to a conserved quantity. Therefore, in the inviscid limit, the maximum potential temperature anomaly should be zero, and the minimum potential temperature anomaly should be -15 K, the same as the initial condition.

As pointed out in Ooyama (2001), it is possible to run the same simulation without any explicit diffusion terms, which results in the inviscid solution. Of course, a truly inviscid solution cannot be reached because numerical viscosity always exists as long as the mesh size is finite. The adjective “inviscid” refers to the formulation without explicit diffusion terms. The inviscid model is useful to test the stability of the numerical scheme. Figure 3 shows the inviscid solution, which generates more rotary eddies and sharper potential temperature gradient than the viscous one. The boundary condition becomes more problematic when explicit diffusion is disabled, and HLLC solver is used. HLLC solver produces unrealistic potential temperature perturbation of more than 12 K at 900s near the surface (Figure 3b), while potential temperature perturbation remains as low as 0.3 K using LMARS solver (Figure 3a). The solution we obtained using the LMARS solver at 600s is similar to Ooyama (2001)’s result in their Figure C3 and Guerra & Ullrich (2016)’s result in their Figure 12, confirming the veracity of the model.

5.2 Robert rising bubble

The second case tests the turbulence field generated by a warm rising bubble (Robert, 1993). The background atmospheric temperature is still an adiabat at 303.15 K and the surface pressure is 1 bar. A warm bubble, of a peak temperature excess 0.5 K, is placed near the bottom of the domain. Since the forcing is weak, the solution is more sensitive to the numerical scheme than the first case. The computational domain is 1 km wide and 1.5 km tall, with a spatial resolution of 5 m. Its evolution with time is displayed in Figure 4 (potential temperature anomaly) and Figure 5 (vertical velocity). The bubble rises and develops two rotors at both sides. Vertical velocity is positive at the arch of the rotor and negative both inside and outside. At about 18 minutes, the Kelvin-Helmholtz instability due to the vertical velocity shear destroys the smooth arch. Two rotors become turbulent and stop rising. The warmest part of the bubble concentrates at its head, protrudes out, and develops a secondary circulation. The turbulent potential temperature field at 18 minutes is identical to the reference solution calculated by Chen et al. (2013) and Guerra & Ullrich (2016). Though simple as it looks, many advanced numerical models failed to produce the correct turbulence field; some models generated more KH instabilities along the arch of the thermal bubble, and others, being too diffusive, did not produce any stability at all (e.g. Konor 2014; Flyer et al. 2016; Abdi et al. 2017). Passing this test indicates that our model can capture the details of the turbulence field induced by weak forcing.

A more difficult problem arises when the rising bubble has a uniform potential temperature. Weak forcing and having a discontinuity in the temperature field make the solution very susceptible to the details of the numerical treatment. Similar to the previous one, the temperature excess is 0.5 K, and the warm bubble is evolved to 10 mins. As far as we know, there is no consensus on what the true solution is. Various authors get similar but somewhat different results. Robert (1993)’s solution (their Figure 2) developed one rotor at 7 mins at the shoulder of the warm bubble while we have two (Figure 6). Chen et al. (2013) did not provide a solution at 5m resolution, but the turbulence field in our result at 10 mins is very similar to theirs.

5.3 Bryan moist bubble

Extending the dry simulation to a moist one with vapor and cloud, we perform the final test, which simulates a saturated bubble rising in a saturated environment. The initial temperature profile is a moist adiabat, and the total water mixing ratio, including vapor and cloud, is constant at all levels. The surface temperature and pressure are 289.85 K and 1 bar, respectively. A similar Gaussian-shaped bubble is placed near the bottom of the domain, but filled with water vapor. The bubble is warmer than the environment by less than 2 K. Thus, it contains slightly more water vapor and slightly less cloud so that it is both thermally and compositionally buoyant (background atmospheric composition is O2+N2). The temperature and the pressure of the bubble drop when it rises, causing water vapor to condense and latent heat released. The domain size is 10 km in the vertical and 20 km in the horizontal, with a spatial resolution of 100 m. Other details of the initial condition are described in Bryan & Fritsch (2002). For readers who are not familiar with moist thermodynamics, we summarize the key equations and definitions in Appendix C.

Because we are using total energy as the thermodynamic variable, the temperature of the background atmosphere is obtained by integrating the moist adiabatic temperature gradient (see equation C3). An adaptive mesh is used to assure that the background atmosphere has a uniform equivalent potential temperature profile as demanded by the test case. We tried two ways of implementing instantaneous condensation when supersaturation occurs. One is to choose a very large value of condensation rate, for example, $10^{9}$ s*-1*. The other is to calculate the equilibrium state iteratively, a process called saturation adjustment (see Appendix D for the numerical scheme). We found that both ways yield the same result.

Figure 7 shows the equivalent potential temperature anomaly and the vertical velocity. Our result is almost identical to the reference solution given by Bryan & Fritsch (2002) at 100m resolution, except that the maximum equivalent potential temperature anomaly (see its definition at equation C9) is 4.4 K at 1000s in our simulation whereas Bryan & Fritsch (2002) got 4.09 K. We suspect that the difference is due to the difference in formulating the thermodynamic equation. Because the equivalent potential temperature is a variable representing entropy, the increase of the maximum value in an adiabatic environment may result from irreversible mixing, which cannot be captured if entropy or potential temperature is used as the prognostic variable unless all the irreversible sources of entropy production are explicitly accounted for. Pressel et al. (2015) provided a list of irreversible entropy production due to mixing and precipitation process. However, fundamental difficulty arises from the estimation of the entropy generation by numerical diffusion due to the discretization of the governing equations. Although we cannot prove rigorously that the increase in the maximum equivalent potential temperature is caused by irreversible mixing, this exercise shows that using total energy as a prognostic variable obeys thermodynamic laws with less effort, and can potentially capture entropy production due to irreversible mixing, which may be an advantage over the traditional models using potential temperature as the prognostic variable.

6 Simulation of Jovian atmosphere in radiative-convective equilibrium

In terms of a real application to planets, we design an idealized model for the Jovian atmosphere, whose primary focus is to investigate how the internal heat is transported through the weather layer and finally radiates into the space at the optically thin level; for that reason, we are interested in the mean profiles of water, ammonia, and temperature in radiative-convective equilibrium (RCE) rather than the initial ones that trigger moist convection (Stoker, 1986; Yair et al., 1995; Hueso & Sánchez-Lavega, 2001). The computational domain is two-dimensional, covering 1667 km in the horizontal and 400 km in the vertical, with a spatial resolution of around 5.5 km in both directions. A doubly periodic horizontal boundary condition is assumed. The model stays one step ahead of the one-dimensional equilibrium condensation models (e.g. Weidenschilling & Lewis 1973) or one-dimensional models with parameterized eddy diffusion (e.g. Ackerman & Marley 2001) by directly simulating the dynamic effects.

We use a spatially homogeneous and temporally constant body cooling rate as a substitute for the radiative cooling process since our main focus is convection. The cooling rate is set to 1 K/day between 2 bars and 0.2 bar, which is about 100 times larger than what Jupiter emits realistically. Choosing a larger cooling rate shortens the period needed to equilibrate. Between 0.2 bar and the top boundary, there is a sponge layer absorbing kinetic energy. The total energy is not modified in the sponge layer, and therefore the damped kinetic energy turns into heat, mimicking the effect of wave breaking in the stratosphere. The sponge layer is formulated by Newtonian relaxation:

[TABLE]

where $\tau_{0}=10^{4}$ s and $p_{0}=0.2$ bars. The choice of the lower boundary condition matters for the dynamics since the lower boundary is the source of energy and is the place where momentum dissipates. Here, we apply a no-slip boundary condition at 200 bars to represent the damping of momentum by magnetohydrodynamic drag that acts deep in the atmosphere (Grote & Busse, 2001; Liu08a), and a fixed temperature of 800 K to represent the internal heat source. Although the actual momentum dissipation level is much deeper than what employed here, this drag provides the principle momentum dissipation mechanism on Jupiter. A thorough discussion of the role of the lower boundary condition shall be canvassed in a dedicated study later.

Two condensible gases, water and ammonia, as well as their clouds and precipitation are included in the model. The mass mixing ratios of water and ammonia at depth are 40 g/kg and 2.7 g/kg respectively, corresponding to a nominal water cloud base at 7 bars and ammonia cloud base at 0.7 bars. The initial temperature profile is a pseudo moist adiabat up to 0.2 bars, then isothermal above. Microphysical rate constants are largely uncertain in the Jovian atmosphere due to the unknown size distribution of clouds and precipitation. We choose nominal values of $k_{1}=10^{9}$ s*-1*, $k_{2}=10^{-4}$ s*-1*, $k_{3}=0$ s*-1*, $k_{4}=10^{-2}$ s*-1*, and $k_{5}=-20$ m/s (see Section 3 for the definition of the coefficients). The extremely large value of $k_{1}$ manifests that condensation occurs instantaneously when vapor is saturated. The autoconversion timescale is $10^{4}$ s, in the middle of the range studied by Sugiyama et al. (2014). The evaporation rate and terminal velocity are roughly estimated since water and ammonia precipitation will eventually boil at 30 bars and 4 bars respectively (Figure 1).

The most prominent and striking feature of the Jovian atmosphere in RCE is that the mean potential temperature profile decreases with altitude between 10 bars and 5 bars (solid red line in Figure 8), suggesting a superadiabatic temperature gradient. However, the seemingly unstable temperature gradient is balanced by a negative gradient of the concentration of water vapor (solid blue line in Figure 8) such that the density profile is still stably stratified, as indicated by a positive gradient of the virtual potential temperature profile (Figure 9). The amount of water and ammonia gas is also much less than their initial values, partly due to the cold atmosphere between 7 bars and 1 bar. This result contrasts the vertical profiles obtained by a one-dimensional equilibrium condensation model (Weidenschilling & Lewis, 1973) that the temperature profile is a moist adiabat and that ammonia and water are well-mixed in the subcloud layer. Since the temperature at 1 bar pressure level is measured as 165 K by Lindal et al. (1981), a superadiabatic temperature gradient indicates that Jupiter’s interior may be warmer than what has been assumed before based on adiabatic extrapolation. Though water vapor contributes to less than one percent of the total mass of the atmosphere in the weather layer, the condensation of water vapor significantly changes the mean molecular weight of the atmosphere, and thus the thermal stratification.

It is illuminating to see how the superadiabatic temperature gradient develops and maintains over time. We have integrated the model for over 1700 days and plotted the evolution of mean temperature anomaly (with respect to the initial temperature) and mean zonal wind in Figure 10. The superadiabatic temperature gradient develops soon after the model started. Along with it, negative wind appears above the cloud, and positive wind appears below the cloud. The wind shear across the water cloud breaks after 200 days, likely triggered by convective instability because a large reduction in temperature above 7 bars occurs simultaneously. An unsteady state featured by alternating zonal wind lasts for about 500 days. The period of the oscillation is about 50 days, but sometimes one oscillatory phase may be skipped, for example, between day 500 and day 600. After 800 days, the atmosphere reaches a steady state. The zonal wind is strongest in the water cloud, reaching 50 m/s, and decreases with depth toward zero. The final equilibrated temperature above 7 bars is slightly warmer than what has been in the unsteady state. No matter whether the steady state is reached or not, a colder atmosphere above the water condensation level is a robust feature of the simulated Jovian atmosphere in RCE.

A snapshot of the spatial distribution of water and ammonia is shown in Figure 11, where the vapors are plotted in colored contours, clouds in shaded contours, and precipitation in scattered dots. Dark areas in each panel indicate regions devoid of the vapor. Unlike the prediction by equilibrium condensation models, the main cloud water appears at about 2 bars, with patchy and cumulus-like cloud scattered around 5 – 7 bars. Large cloud concentration usually coincides with clustered precipitation. Some precipitation deposits to levels deeper than 10 bars but none exists below 30 bars because of boiling. Ammonia cloud, on the other hand, forms at 0.7 bar, similar to the prediction by equilibrium condensation models, but the concentration of ammonia gas shows large variability in the subcloud layer. Some regions bear more ammonia than the deep atmosphere in mass mixing ratios, due to the re-evaporation of the ammonia precipitation, others are significantly less because compensating downdrafts carry dry air downward. On average, the concentration of the ammonia gas increases with depth and reaches its deep value at about 10 bars, where the atmosphere turns neutrally stable. Similar depletion of ammonia gas in the subcloud level was also found in the simulation result of Sugiyama et al. (2014), although the level where ammonia recovers is deeper in our simulation than in Sugiyama et al. (2014)’s.

The recent observations from the Juno spacecraft revealed the thermal emission of Jupiter’s atmosphere from depth well below the water cloud (Bolton et al., 2017). Yet, the interpretation of its microwave spectra is confounded by the nonuniform distribution of ammonia gas (Li et al., 2017; Ingersoll et al., 2017). It remains a mystery that why ammonia, an essentially chemically inert species, is depleted in the atmosphere far below its cloud condensation level. The previous study attributed the desiccation of the ammonia gas after a single convective event to the compensating downdrafts that brought dry air down from higher up (Li & Ingersoll, 2015). However, it is unclear how this mechanism works in a steady state in which the atmosphere is continuously heated from below and cooled from the top. This simulation demonstrates that moist convection alone can cause the ammonia gas to deplete at levels well below its nominal cloud bottom. Therefore, the effects of dynamics must be considered in modeling the vertical structure of a planet’s atmosphere.

7 Conclusion and discussion

In this work, we presented a new non-hydrostatic and cloud-resolving model for planets. It is developed with the goal of simulating cloud formation on both terrestrial planets and giant planets. We build this model on top of the Athena++ framework, allowing it to exploit the features provided by the Athena++ code, such as static/adaptive mesh-refine, curvilinear geometry, and dynamic task scheduling. In order to correctly and efficiently simulate atmospheric flows, we have implemented several extensions to the original Athena++ solver. They are summarized in the following.

The continuity equations are extended to vapors, clouds and precipitation. 2. 2.

Latent heat terms are added to the energy equations. 3. 3.

The equation of state is modified to represent a heterogeneous mixture of gases and clouds. 4. 4.

Microphysical processes are formulated generically such that both terrestrial and Jovian planets can use the same framework. Specifically, we discussed the formulation of condensation, precipitation and evaporation. 5. 5.

The LMARS Riemann solver is used instead of the Roe or HLLC solver for it better preserves hydrostatic balance in the vertical. 6. 6.

A fifth-order WENO reconstruction method is implemented to achieve a high order of accuracy as well as to minimize the phase error in the wave propagation.

We pioneered the use of total energy as the prognostic variable for the thermodynamic equation while a more conventional choice is potential temperature. Potential temperature, a concept primary for interpretive purposes, becomes complicated when the heat capacity of the atmosphere varies with temperature or when multiple condensing species exist in the atmosphere as so in the Jovian atmosphere. Some authors prefer using entropy rather than potential temperature for Earth’s atmosphere (e.g. Hauf & Höller 1987; Pressel et al. 2015). However, fundamental difficulties arise when one wishes to account for all entropy production due to irreversible processes, including the entropy generation by numerical diffusion due to the discretization of the governing equations. Therefore, we argue that using total energy is a better choice since total energy is conserved regardless of whether the process is reversible or not. We carried out a series of numerical tests to demonstrate the veracity of the model. The simulation results are almost identical to the reference solutions in the literature.

Then, we designed an idealized simulation of the Jovian atmosphere in radiative-convective equilibrium, with water and ammonia included as condensible species. Our simulation result showed that the concentration of the ammonia gas is highly variable in the subcloud layer, with some region enriched in ammonia caused by re-evaporation and some region developed in ammonia caused by downdrafts. On average, the subcloud layer is depleted in ammonia gas and the depletion can reach 10 bars in pressure. Although, at this stage, we cannot claim that this is exclusively the cause of the observed depletion of ammonia in Jupiter’s atmosphere (Li et al., 2017) because we have not yet considered the effect of large-scale circulation, the Coriolis force, and the three-dimensionality, we can assure that moist convection alone can cause the ammonia gas to deplete at levels well below its nominal cloud bottom. As a result, the concentration profiles resulted from the equilibrium condensation calculation should be interpreted qualitatively but not quantitatively.

Moreover, in contrast to the assumed adiabaticity of Jupiter’s atmosphere, we found that the temperature gradient is superadiabatic near the water condensation level, due to the changing of the mean molecular weight. We are certain about the mechanism but less sure on the exact magnitude because we have not yet tested the effect of the lower boundary condition and the effect of the cooling rate, which we will leave for future exploration.

We expect many future applications of the cloud-resolving model to planetary atmospheres. First, the speed of sound on Mars is about 2/3 of what it is on Earth, but the pressure fluctuation is much larger. Particularly, Mars has been shown in the MACDA dataset to have transonic jet streaks (Lewis et al., 2017). Thus, our model, using the LAMRS Riemann solver, offers advantages to investigate the dynamics of the strong polar jet. Other planets with tenuous atmospheres, such as Io and Pluto, will also benefit on this point.

Second, many GCMs for giant planets have been developed with an active hydrological cycle and microphysics (Palotai & Dowling, 2008; Lian & Showman, 2010). But a GCM cannot directly resolve moist convection due to the coarse resolution and hydrostatic assumption. Our cloud-resolving model can complement the large-scale models by devising convective parameterization schemes. Conversely, large-scale models help our simulation by providing the context and forcing that drive the cloud-resolving model.

Third, advances in the observations of exoplanets have further broadened the potential application of a cloud-resolving model. Many exoplanets’ atmospheres are thought to harbor multiple layers of exotic clouds made of KCl, ZnS, MgSiO3, etc. (Sing et al., 2016; Morley et al., 2012; Crossfield et al., 2014). Understanding their vertical structure and extend is important for interpreting the transit spectra. However, most cloud models for exoplanets’ atmospheres are restricted to 1D, in which physical processes are mingled in an ad-hoc parameters, $K_{zz}$ (Ackerman & Marley, 2001; Zhang & Showman, 2018a, b). There is no physical basis behind which planet should have what values of $K_{zz}$ other than it fits the observation (Li et al., 2014). Able to resolve the turbulence field, a cloud-resolving model will lay a solid foundation for the physical processes governing the distribution and structure of the clouds.

We thank all the people that are working or have worked on the development of Athena++ code, and hosted it on an open source platform. It is supposed to be a software that solves astrophysical problems but we found no difficulty to extend the solver to atmospheric flows. We also thank Andrew Ingersoll, Xi Zhang, Zhaohuan Zhu, Zhihong Tan, Xianyu Tan and Kyle Pressel for stimulating discussions of the development of the model. C.L. is supported by the Juno mission.

Appendix A List of symbols

All important symbols used in the manuscript and their CGI units are collected in the following table for reference. \startlongtable

Appendix B Conversion between primitive variables and conserved variables

The primitive variables are :

[TABLE]

The conserved variables are:

[TABLE]

Converting the primitive variables to the conserved variables takes two steps. The first step converts specific humidities to densities and velocities to momentums.

[TABLE]

The second step converts total pressure to total energy. Total pressure is the sum of the partial pressure of each gaseous component:

[TABLE]

Total energy includes internal energy, kinetic energy, and chemical potentials:

[TABLE]

where $K=\frac{1}{2}\rho(u^{2}+v^{2}+w^{2})$ is the kinetic energy and $A=\sum_{i}\mu_{i}\rho_{i}+\sum_{i,j}\mu_{ij}\rho_{ij}$ is the chemical potential. Substituting equation (B9) into equation (B10) to replace the temperature gives:

[TABLE]

Converting the conserved variables to the primitive variables is the reverse of the previous steps. First, the specific humidities are obtained by letting:

[TABLE]

Define internal energy as:

[TABLE]

Then the formula for pressure is:

[TABLE]

Now we can define an effective adiabatic index $\gamma$ for an heterogeneous air parcel:

[TABLE]

It can be shown that this definition is consistent with the definition of potential temperature in equation (C6). In the special case in which there are no clouds ($q_{c_{i}}=0$) and the adiabatic indices of all vapors equal that of dry air:

[TABLE]

The adiabatic index of the air mixture equals that of the dry air: $\gamma=\gamma_{d}$, and equation (B19) reduces to the equation of state of a single-component ideal gas regardless of the differences in molecular weight:

[TABLE]

Appendix C Thermodynamic formulas

In this section, we derive the expression of several important thermodynamic variables. Some of the expressions are available in a standard textbook (e.g. Emanuel 1994). But often, these expressions only consider one condensing species, water vapor. We extend them to multiple condensing species and emphasize on the quantities that will be used to diagnose or initialize the model.

C.1 Moist static energy

A conserved quantity during adiabatic displacement in which the pressure change is hydrostatic is moist static energy $h_{m}$. The generalized expression is:

[TABLE]

If the atmosphere is at rest, using equations (8), (9) and (C1), it is easy to verify that:

[TABLE]

C.2 Moist adiabatic temperature gradient

The moist adiabatic temperature gradient in pressure coordinate for an atmosphere with multiple condensing species is derived in Li et al. (2018), which reads:

[TABLE]

where $r_{i}=q_{i}/q_{d}$ and $r_{ij}=q_{ij}/q_{d}$ are the mass mixing ratios of vapor $i$ and cloud $(i,j)$ with respect to the dry air respectively. Using the hydrostatic balance and the equation of state (equation 15), the temperature gradient in height coordinate is:

[TABLE]

Either equation (C1) or equation (C4) can be used to initialize a background atmosphere that is neutral to moist convection. In a model that uses height coordinate, directly integrating equation (C4) is often easier than an iterative approach to find a solution that satisfies the conservation of moist static energy. Refinement in the grid size is sometimes necessary in order to achieve high accuracy, especially in the grid where condensation occurs. The pros and cons of different approaches are discussed in Li et al. (2018).

C.3 Potential temperature

The potential temperature is the temperature that an air parcel would have when it is adiabatically displaced to the reference pressure without having phase change. Let the latent heat terms ($L_{ij}$) to be zero, equation (C3) reduces to the dry adiabatic temperature gradient:

[TABLE]

Therefore, the potential temperature is obtained by integrating equation (C5):

[TABLE]

Similar to the definition of virtual temperature, the virtual potential temperature is defined as:

[TABLE]

The gradient of the virtual potential temperature is approximately related to the gravity wave frequency, the Brunt-Väisälä frequency:

[TABLE]

C.4 Equivalent potential temperature

The equivalent potential temperature is a measurement of entropy taken the latent heat into account. It is conserved during adiabatic and reversible phase change. In a saturated atmosphere with a single condensible species, the formula of equivalent potential temperature is given by equation 4.5.11 in Emanuel (1994):

[TABLE]

where $q_{t}=q_{i}+q_{ij}$ is the total mass mixing ratio of the vapor and the cloud. Note that this expression is only valid for a single condensible species and a single condensate. So, the expression is rarely used in our model but is described for comparison purpose because one of the test case, the Bryan moist bubble, uses equivalent potential temperature as a diagnostic variable.

Appendix D Saturation adjustment scheme

It is easier to work with molar mixing ratios rather than mass mixing ratios when dealing with thermodynamics. Therefore, in this section, we will work with molar quantities. A variable with a hat symbol ($\hat{}$) on its top represents the molar equivalence of the usual meaning defined in Table 1. Thus, the thermodynamic state variables are:

[TABLE]

The following two equations relate the mass and molar representations:

[TABLE]

Vapor condenses to cloud when its partial pressure exceeds the saturation vapor pressure. In a numerical model that uses control volume, condensation process happens under the condition of constant volume and constant total energy instead of the usual assumption of isobaric process. The saturation adjustment scheme efficiently calculates the transfer of mass from vapor to cloud or vice versa that would maintain thermodynamic equilibrium between vapor and cloud. There are a number of ways to calculate the saturation vapor pressure including using a look-up table or using an empirical formula such as the Antoine equation. We directly integrate the Clausius-Clapeyron relation because we will use the gradient of the saturation curve, which is exactly the Clausius-Clapeyron relation, to derive the adjustment scheme. The Clausius-Clapeyron relation is:

[TABLE]

We have used equation (9) to correct the temperature dependence of the latent heat. Integrating the above Clausius-Clapeyron relation gives the saturation vapor pressure:

[TABLE]

where $p^{r}_{ij}$ is the saturation vapor pressure at reference temperature $T^{r}$, and $t,\beta_{ij},\delta_{ij}$ are non-dimensional parameters:

[TABLE]

Assuming vapor $i$ and its cloud $(i,j)$ undergo a phase change. The transfer of substance in molar amount is $\Delta\hat{q}$. The partial pressure of the vapor before is:

[TABLE]

where $g=\hat{q}_{d}+\sum_{j\neq i}\hat{q}_{j}$ is the sum of the molar mixing ratios of other gases. Ignoring the change of pressure for now, after condensation, the temperature changes to $T+\Delta T$ and the partial pressure becomes the saturation vapor pressure:

[TABLE]

Let $s(T)=p_{ij}^{*}(T)/p$, then:

[TABLE]

Expand the above equation near $T$ and keep the first order term:

[TABLE]

Using the Clausius-Clapeyron relation:

[TABLE]

Let $\alpha=\Delta T/(T\Delta\hat{q})$ during the phase change. Substituting $s^{\prime}(T)$ in equation (D12) with equation (D13) and $\alpha$ gives:

[TABLE]

Under isobaric process, enthalpy $h$ is conserved:

[TABLE]

Therefore,

[TABLE]

Under isochoric process, energy is conserved:

[TABLE]

Therefore,

[TABLE]

Iteration is needed in the presence of multiple condensing species. Starting from each iteration, equation (D14) is used to calculate the projected mass transfer. The value of $\Delta\hat{q}$ must be limited to ensure the positiveness of all substances. Then either equation (D15) or equation (D17) is used to update temperature depending on whether it is an isobaric or isochoric process, i.e. whether the vertical coordinate is pressure or height. If the phase change happens at constant volume, pressure needs to be updated accordingly:

[TABLE]

where

[TABLE]

is another expression of virtual temperature using molar mixing ratios (c.f. equation 16). Next iteration begins with the updated temperature, pressure and abundances until the solution converges. Because the first derivative of the saturation vapor pressure is taken care of, a quadratic convergence is expected. In real applications, we find that the iteration converges to machine precision in less than 5 iterations.