THOR 2.0: Major Improvements to the Open-source General Circulation Model

With new technology and data analysis techniques, we are entering an era in which 3D models of exoplanet atmospheres can be tested and validated. As observations improve, it will be important to test a variety of models, all of which make various assumptions in their representations of physical processes, to create the most accurate interpretations of the data.

Numerous theoretical and computational studies have shown that hot Jupiters have large day/night temperature contrasts and equatorial superrotation (Showman & Guillot 2002; Cooper & Showman 2005; Dobbs-Dixon & Lin 2008; Showman et al. 2009; Dobbs-Dixon et al. 2010; Rauscher & Menou 2010, 2012; Heng et al. 2011a, 2011b; Dobbs-Dixon & Agol 2013; Perez-Becker & Showman 2013; Kataria et al. 2015; Amundsen et al. 2016). These features are broadly consistent across a wide range of models and have been validated by observations (Snellen et al. 2010; Knutson et al. 2012; Louden & Wheatley 2015). The general consensus states that the superrotation is the product of interacting equatorial Rossby and Kelvin waves, which allow angular momentum to be transported toward the equator (Showman & Polvani 2011; Tsai et al. 2014; Hammond & Pierrehumbert 2018; Mendonça 2019).

Other works have explored the importance of clouds and hazes (Heng et al. 2012; Helling et al. 2016; Lee et al. 2016; Roman & Rauscher 2017; Mendonça et al. 2018a), atmospheric chemistry (Cooper & Showman 2006; Parmentier et al. 2013; Kataria et al. 2016; Drummond et al. 2018a, 2018b; Mendonça et al. 2018b), and shock physics (Goodman 2009; Li & Goodman 2010; Heng & Workman 2014; Fromang et al. 2016; Koll & Komacek 2018). Still others have focused on initial conditions and the details of numerical techniques (Thrastarson & Cho 2010, 2011; Watkins & Cho 2010; Liu & Showman 2013; Polichtchouk et al. 2014). These additional features are thought to be important, but the data are less conclusive in this regard (Heng & Showman 2015).

Further studies have focused on the atmospheric dynamics of other types of planets, including cooler Neptune-size planets (Charnay et al. 2015; May & Rauscher 2016; Mayne et al. 2019) and terrestrial planets (Merlis & Schneider 2010; Carone et al. 2014, 2016, 2018; Kaspi & Showman 2015; Guendelman & Kaspi 2018, 2019), particularly for the purposes of understanding habitability (Williams & Pollard 2002, 2003; Abe et al. 2011; Leconte et al. 2013a, 2013b; Yang et al. 2013, 2014; Kopparapu et al. 2016; Wolf et al. 2017; Way et al. 2018; Jansen et al. 2019). Atmospheric constraints for smaller and cooler planets, however, remain scarce because of the comparative difficulty of observation.

Even though constraining the atmospheric processes remains a challenge for extrasolar planets, the wide variety of orbits, masses, and sizes of planets indicates that the atmospheres will be quite diverse. Most studies of exoplanet atmospheres using general circulation models (GCMs) have adapted codes developed for Earth or solar system planets. The virtue of this method is that it relies on models that have been well tested; however, planet-specific tunings are often built into the code and can be difficult to find and generalize. In the development of THOR, we have chosen the opposite path: to develop a code from scratch that would be completely free of planet-specific tunings and would thus provide a flexible tool for the study of a diverse range of atmospheres. The additional benefit of this path is that we develop an intimate understanding of the physical processes at work and how these are represented in the code. The challenge that remains is the workload associated with development and testing new components.

Presented in Mendonça et al. (2016) and Mendonca et al. (2018) and utilized in Mendonça et al. (2018a, 2018b), THOR⁴ is a nonhydrostatic, fully global, 3D general circulation model developed specifically for the study of exoplanets. As such, it is free from the Earth and solar system tunings that often make use of 3D GCMs a challenge for exoplanets. However, because it is a young model developed from scratch, much development remains in order to make the model applicable to all types of planets. This work represents a step forward along this path.

The goals of this paper are to consolidate descriptions of improvements that have been made to the model since Mendonça et al. (2016), clarify the model framework, validate the new physics, and compare results from the model using different approximations, with implications for the general circulation of hot Jupiters.

The principal equations solved in THOR are the flux forms of the Euler equations:

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho }{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho {\boldsymbol{v}}}{\partial t}+{\rm{\nabla }}\cdot (\rho {\boldsymbol{v}}\otimes {\boldsymbol{v}})=-{\rm{\nabla }}P-\rho g\hat{{\boldsymbol{r}}}-2\rho {\boldsymbol{\Omega }}\times {\boldsymbol{v}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho \theta }{\partial t}+{\rm{\nabla }}\cdot (\rho \theta {\boldsymbol{v}})=0,\end{eqnarray} \tag{ 3 }$$

where ρ is the density, v is the velocity, P is the pressure, g is the acceleration due to gravity (assumed to be constant), ${\boldsymbol{\Omega }}$ is the planet's rotation vector, and θ is the potential temperature. Potential temperature is defined as

$$\begin{eqnarray}&&\theta \equiv T{\left(\displaystyle \frac{P}{{P}_{\mathrm{ref}}}\right)}^{-{R}_{d}/{C}_{P}}.\end{eqnarray} \tag{ 4 }$$

The system of equations (Equations (1)–(3)) is closed by the ideal gas law, P = ρR_dT. Additionally, a relationship between the pressure and potential temperature is necessary to update the pressure for use in Equation (2). From the ideal gas law and the definition of potential temperature, we have

$$\begin{eqnarray}&&P={P}_{\mathrm{ref}}{\left(\displaystyle \frac{{R}_{d}\rho \theta }{{P}_{\mathrm{ref}}}\right)}^{{C}_{P}/{C}_{V}}.\end{eqnarray} \tag{ 5 }$$

THOR solves the Euler equations using a finite-volume method on an icosahedral grid (Tomita et al. 2001; Tomita & Satoh 2004; Mendonça et al. 2016). The horizontal resolution is controlled by a single parameter, g_level, the number of times the icosahedral grid is refined (i.e., the number of times the sides of the icosahedron are subdivided into smaller triangles). The average angular size of the control volumes is given by

$$\begin{eqnarray}&&\bar{\theta }=\sqrt{\displaystyle \frac{2\pi }{5}}\displaystyle \frac{1}{{2}^{{g}_{\mathrm{level}}}}.\end{eqnarray} \tag{ 6 }$$

The lowest value of g_level used in this work is 4, which results in an angular resolution of $\bar{\theta }\approx 4^\circ$. Every increase of g_level by 1 decreases $\bar{\theta }$ by half. The average size (in m) of the control volumes is simply

$$\begin{eqnarray}&&\bar{d}={r}_{0}\bar{\theta },\end{eqnarray} \tag{ 7 }$$

where r₀ is the planet radius. The value of $\bar{d}$ is used to scale the numerical diffusion coefficients (Section 2.4).

A full description of the THOR algorithm was presented in Mendonça et al. (2016); however, there were a number of typographical errors in that paper and some of the details have changed, so we include a description here.

As described in Mendonça et al. (2016), we use a time-splitting algorithm based on Wicker & Skamarock (2002), Tomita & Satoh (2004), and Skamarock & Klemp (2008). In this scheme, the fluid equations are split into fast and slow modes, and the fast modes are integrated using a smaller time step than the slow modes. The time-stepping loop consists then of an outer loop (the large time step, which has variable length) and an inner loop (small time step, with length designated Δτ).

During the inner loop (at time τ or τ + Δτ), the deviation of any quantity from its large time step value (at time t) is

$$\begin{eqnarray}&&{{\rm{\Phi }}}^{\star [\tau ]}={{\rm{\Phi }}}^{[\tau ]}-{{\rm{\Phi }}}^{[t]},\end{eqnarray} \tag{ 8 }$$

or

$$\begin{eqnarray}&&{{\rm{\Phi }}}^{\star [\tau +{\rm{\Delta }}\tau ]}={{\rm{\Phi }}}^{[\tau +{\rm{\Delta }}\tau ]}-{{\rm{\Phi }}}^{[t]},\end{eqnarray} \tag{ 9 }$$

where the ${}^{\star }$ superscript indicates the deviation and the superscript in square brackets indicates how frequently the values are updated: slow modes (${}^{[t]}$) are updated every large time step, and fast modes (^[τ] or ${}^{[\tau +{\rm{\Delta }}\tau ]}$) are updated every small time step. Broadly, the fast modes are those terms that are associated with acoustic waves, and the slow modes are everything else. This time-splitting method allows for a less stringent constraint on the time step and thus a moderate boost in performance, as many of the terms do not have to be computed as frequently (the advection terms are particularly costly). Further, the 3D operators are split into horizontal and vertical components (Mendonça et al. 2016)

$$\begin{eqnarray}&&{\rm{\nabla }}{\rm{\Phi }}={{\rm{\nabla }}}_{h}{\rm{\Phi }}+\displaystyle \frac{\partial {\rm{\Phi }}}{\partial r}\hat{{\boldsymbol{r}}}\end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{\Phi }}={{\rm{\nabla }}}_{h}\cdot {\boldsymbol{\Phi }}+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}({r}^{2}{{\rm{\Phi }}}_{r}).\end{eqnarray} \tag{ 11 }$$

The vector $\hat{{\boldsymbol{r}}}$ represents the vertical direction, and r is the radial distance from the center of the planet. The altitude is given by z = r − r₀.

Equations (1)–(3) are then discretized as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\rho }^{\star [\tau +{\rm{\Delta }}\tau ]}-{\rho }^{\star [\tau ]}}{{\rm{\Delta }}\tau }+{{\rm{\nabla }}}_{h}\cdot {\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}{r}^{2}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}\\ \quad \ =\ -{{\rm{\nabla }}}_{h}\cdot {\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{[t]}-\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}{r}^{2}{\left(\rho {v}_{r}\right)}^{[t]}+{{ \mathcal F }}_{\rho }^{[t]},\end{array}\,\end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}-{\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{\star [\tau ]}}{{\rm{\Delta }}\tau }+{{\rm{\nabla }}}_{h}{P}^{\star [\tau ]}\\ \quad \ =-{{\rm{\nabla }}}_{h}{P}^{[t]}-\ {{\boldsymbol{ \mathcal A }}}_{h}^{[t]}-{{\boldsymbol{ \mathcal C }}}_{h}^{[t]}+{{\boldsymbol{ \mathcal F }}}_{{{\boldsymbol{v}}}_{h}}^{[t]}+{{\boldsymbol{ \mathcal G }}}_{{{\boldsymbol{v}}}_{h}}^{[\tau ]},\end{array}\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}-{\left(\rho {v}_{r}\right)}^{\star [\tau ]}}{{\rm{\Delta }}\tau }+\displaystyle \frac{\partial }{\partial r}{P}^{\star [\tau +{\rm{\Delta }}\tau ]}+{\rho }^{\star [\tau +{\rm{\Delta }}\tau ]}g\\ \quad \ =\ -\displaystyle \frac{\partial }{\partial r}{P}^{[t]}-{\rho }^{[t]}g-{{ \mathcal A }}_{r}^{[t]}-{{ \mathcal C }}_{r}^{[t]}+{{ \mathcal F }}_{{v}_{r}}^{[t]},\end{array}\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{\left(\rho \theta \right)}^{[\tau +{\rm{\Delta }}\tau ]}-{\left(\rho \theta \right)}^{[\tau ]}}{{\rm{\Delta }}\tau }+{{\rm{\nabla }}}_{h}\cdot {\theta }^{[t]}{\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{[\tau +{\rm{\Delta }}\tau ]}\\ \quad \ +\ \displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}{r}^{2}{\theta }^{[t]}{\left(\rho {v}_{r}\right)}^{[\tau +{\rm{\Delta }}\tau ]}=0.\end{array}\end{eqnarray} \tag{ 15 }$$

The terms ${{\boldsymbol{ \mathcal A }}}_{h}$ and ${{ \mathcal A }}_{r}$ represent the horizontal and vertical components of the advection term, ${\rm{\nabla }}\cdot (\rho {\boldsymbol{v}}\otimes {\boldsymbol{v}})$, and the terms ${{\boldsymbol{ \mathcal C }}}_{h}$ and ${{ \mathcal C }}_{r}$ represent the same for the Coriolis, $2\rho {\boldsymbol{\Omega }}\times {\boldsymbol{v}}$. The terms ${{ \mathcal F }}_{\rho }$, ${{\boldsymbol{ \mathcal F }}}_{{{\boldsymbol{v}}}_{h}}$, and ${{ \mathcal F }}_{{v}_{r}}$ represent fluxes from the "slow" drag or numerical dissipation mechanisms, in this case, hyperdiffusion and Rayleigh friction. The additional term ${{\boldsymbol{ \mathcal G }}}_{{{\boldsymbol{v}}}_{h}}$ represents the 3D divergence damping, which needs to be evaluated on the small time step.

An important note regarding the coordinate system used in Equation (13): the 2D spherical surface represented by this equation is transformed into a 3D Cartesian coordinate system centered on the planet's core and rotating with the planet. For example, the horizontal velocity is defined as

$$\begin{eqnarray}&&{{\boldsymbol{v}}}_{h}={v}_{1}{\hat{{\boldsymbol{e}}}}_{1}+{v}_{2}{\hat{{\boldsymbol{e}}}}_{2}+{v}_{3}{\hat{{\boldsymbol{e}}}}_{3}-{\boldsymbol{v}}\cdot \hat{{\boldsymbol{r}}},\end{eqnarray} \tag{ 16 }$$

where $\hat{{{\boldsymbol{e}}}_{i}}$ represent the axes of this coordinate system and v_i are the total (horizontal and vertical) velocities in the corresponding directions. The radial unit vector is related to the Cartesian coordinate system by

$$\begin{eqnarray}&&\hat{{\boldsymbol{r}}}=\cos \phi \cos \lambda \,\hat{{{\boldsymbol{e}}}_{1}}+\cos \phi \sin \lambda \,\hat{{{\boldsymbol{e}}}_{2}}+\sin \phi \,{\hat{{\boldsymbol{e}}}}_{3},\end{eqnarray} \tag{ 17 }$$

where ϕ is latitude and λ is longitude. In essence, Equation (16) defines the horizontal velocity as the total velocity minus the radial component. The advection and Coriolis terms, ${{\boldsymbol{ \mathcal A }}}_{h}$ and ${{\boldsymbol{ \mathcal C }}}_{h}$, are defined in the same fashion. The use of a 3D Cartesian coordinate system allows the horizontal and vertical components of the advection and Coriolis terms (in Equation (2)) to be cleanly separated and is also used in the NICAM GCM (Tomita & Satoh 2004). The separation, in turn, is key to allowing explicit integration in time in the horizontal dimensions and implicit integration in the vertical dimension (see Section 2.3).

Our prognostic variables are ρ, P, ρv_r, and the three Cartesian components of ρ v_h, as well as their deviations. Despite using a thermodynamic equation for potential temperature, pressure is used as a prognostic rather than θ or ρθ. As noted in Mendonça et al. (2016), in THOR we do not calculate the deviation of ρθ from the large time step, because such calculation can lead to numerical instability. Instead, we calculate the new value, ${\left(\rho \theta \right)}^{[\tau +{\rm{\Delta }}\tau ]}$, from Equation (15) and use this to update the pressure deviation, P^⋆. Similarly, the dissipation and heating terms are applied to the pressure deviation, so that

$$\begin{eqnarray}\begin{array}{rcl}{P}^{\star [\tau +{\rm{\Delta }}\tau ]} & = & {P}_{\mathrm{ref}}{\left(\displaystyle \frac{{R}_{d}{\left(\rho \theta \right)}^{[\tau +{\rm{\Delta }}\tau ]}}{{P}_{\mathrm{ref}}}\right)}^{{C}_{P}/{C}_{V}}\\ & & -{P}^{[t]}+\ {{ \mathcal F }}_{P}^{[t]}{\rm{\Delta }}\tau +\displaystyle \frac{{R}_{d}}{{C}_{V}}{q}_{\mathrm{heat}}{\rm{\Delta }}\tau ,\end{array}\end{eqnarray} \tag{ 18 }$$

where q_heat represents all of the additional sources of heating or cooling (currently only radiation or Newtonian cooling).

Table 1 gives an overview of the variables used during integration of the dynamical core, the roles of each, and their properties. All quantities are defined at the horizontal centers of the control volumes, though there is some staggering of the grid in the vertical.

Table 1.  Summary of Physical Quantities in Equations (12)–(18) and (35)–(36)

Variable	Role	Location	Update
ρ	Prog.	Center	Large step
P	Prog.	Center	Large step
ρ vh	Prog.	Center	Large step
ρvr	Prog.	Midpoint	Large step
ρ⋆	Prog.	Center	Small step
P⋆	Prog.	Center	Small step
${\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{\star }$	Prog.	Center	Small step
${\left(\rho {v}_{r}\right)}^{\star }$	Prog.	Midpoint	Small step
θ	Diag.	Center	Large step
ρθ	Diag.	Center	Small step
${{\boldsymbol{ \mathcal A }}}_{h}$	Diag.	Center	Large step
${{ \mathcal A }}_{r}$	Diag.	Center	Large step
${{\boldsymbol{ \mathcal C }}}_{h}$	Diag.	Center	Large step
${{ \mathcal C }}_{r}$	Diag.	Center	Large step
h	Diag.	Center	Large step
$\tilde{g}$	Diag.	Center	Large step
${{ \mathcal F }}_{{\rm{\Phi }}}$	Diag.	Center	Large step
${{\boldsymbol{ \mathcal G }}}_{{{\boldsymbol{v}}}_{h}}$	Diag.	Center	Large step
qheat	Diag.	Center	Master step

Note. The "Role" column indicates whether the variable is a primary variable of the integration (prognostic) or a secondary variable (diagnostic). The "Location" column indicates whether variables are defined at the center of the vertical layers or the midpoint between two layers (though many quantities are interpolated to the midpoints to solve Equation (35)). The "Update" column indicates when each quantity is updated (Master step refers to an update outside the dynamical core loop; see Section 2.5).

Download table as:  ASCIITypeset image

Table 4.  Model Parameters for Gray RT Simulations

Symbol	Description	Units	Earth	HD 189733 b
r0	Planet radius	m	6371,000	79,698,540
g	Gravity	m s−2	9.8	21.4
Ω	Rotation rate	rad s−1	7.292 × 10−5	3.279 × 10−5
Rd	Gas constant	J K−1 kg−1	287	3779
CP	Atmospheric heat capacity	J K−1 kg−1	1005	13226.8
Pref	Reference pressure (bottom boundary)	bars	1	220
Tinit	Initial temperature of atmosphere	K	300	1600
ΔtM	Time step	s	500/300	150/100
ztop	Altitude of model top	m	36000	3.6 × 106
glevel	Grid refinement level	⋯	5/6	4/5
vlevel	Number of vertical levels	⋯	32	40
Dhyp	Hyperdiffusion coefficient	⋯	4.8 × 10−3	9.97 × 10−3
Ddiv	Divergence damping coefficient	⋯	4.8 × 10−3	9.97 × 10−3
S0	Incident stellar flux	W m−2	1367	467072
A0	Albedo	⋯	0.3135	0.18
norb	Orbital mean motion	rad s−1	1.98 × 10−7	3.279 × 10−5
nlw	Power-law index for long-wave optical depth	⋯	4	2
fl	Strength of well-mixed absorber (long-wave)	⋯	0.1	0.5
τlw,eq	Long-wave optical depth at Pref at equator	⋯	6	4680
τlw,pole	Long-wave optical depth at Pref at poles	⋯	1.5	4680
nsw	Power-law index for short-wave optical depth	⋯	2	1
τsw	Short-wave optical depth at Pref	⋯	0.2	1170
${ \mathcal D }$	Diffusivity factor	⋯	1	2
Csurf	Heat capacity of surface (lower boundary)	J K−1 m−2	107	⋯
ksurf	Friction coefficient of lower boundary	s−1	1.1574 × 10−5	⋯
σb	Boundary layer top (fraction of Pref)	⋯	0.7	⋯
ksp	Sponge layer strength	s−1	⋯	10−3
ηsp	Bottom of sponge layer (fraction of ztop)	⋯	⋯	0.8
nlats	Number of latitude bins used in sponge layer	⋯	⋯	20

Download table as:  ASCIITypeset image

Rather than solving Equation (14) explicitly in time for the vertical momentum, ${\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}$, we follow Tomita & Satoh (2004) in combining the continuity equation, vertical momentum equation, and thermodynamic equation to form a 1D Helmholtz equation that can be solved implicitly. The implicit solution has the advantage of stabilizing the model without having to resolve the timescale associated with vertically propagating acoustic waves. The resulting Helmholtz equation was presented in Mendonça et al. (2016); however, there were a number of typographical errors, and some steps of the derivation were not particularly clear. We take the opportunity to reproduce the full derivation and to correct the prior mistakes here. We must stress that the typos present in Mendonça et al. (2016) did not propagate to the code itself—in other words, the model was coded correctly, despite typos in the manuscript.

2.3.1. Preparing the Thermodynamic Equation

The use of the entropy equation, once discretized (Equation (15)), does not result in the Helmholtz equation presented in Mendonça et al. (2016) (Equation (37) in that paper). Rather, it is the energy form of the thermodynamic equation that is used (Equation (3) of Tomita & Satoh 2004):

$$\begin{eqnarray}&&\displaystyle \frac{\partial \rho e}{\partial t}+{\rm{\nabla }}\cdot (h\rho {\boldsymbol{v}})={\boldsymbol{v}}\cdot {\rm{\nabla }}P+{q}_{\mathrm{heat}},\end{eqnarray} \tag{ 19 }$$

where ρ is the density, e is the specific internal energy, h is the specific enthalpy, v is the velocity of the fluid, P is the pressure, and q_heat is the diabatic heating rate. It can be shown that this equation is equivalent to Equation (3), aside from the heating term (see, e.g., Section 1.6 of Vallis 2006). However, once discretized, the discrete form for the pressure flux cannot be easily derived from the entropy version, so we begin here with the energy form. For reference, the specific internal energy and specific enthalpy are defined as

$$\begin{eqnarray}&&e={C}_{V}T,\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&h={C}_{P}T.\end{eqnarray} \tag{ 21 }$$

The internal energy can be written as E_int = ρe, which is also related to the pressure via the adiabatic gas index, ${\gamma }_{\mathrm{ad}}$:

$$\begin{eqnarray}&&{E}_{\mathrm{int}}=\displaystyle \frac{P}{{\gamma }_{\mathrm{ad}}-1}.\end{eqnarray} \tag{ 22 }$$

This allows us to write our thermodynamic equation as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial P}{\partial t}+({\gamma }_{\mathrm{ad}}-1){\rm{\nabla }}\cdot (h\rho {\boldsymbol{v}})=({\gamma }_{\mathrm{ad}}-1){\boldsymbol{v}}\cdot {\rm{\nabla }}P\\ \quad \ +\ ({\gamma }_{\mathrm{ad}}-1){q}_{\mathrm{heat}},\end{array}\end{eqnarray} \tag{ 23 }$$

or, since ${\gamma }_{\mathrm{ad}}-1={R}_{d}/{C}_{V}$,

$$\begin{eqnarray}&&\displaystyle \frac{\partial P}{\partial t}+\displaystyle \frac{{R}_{d}}{{C}_{V}}{\rm{\nabla }}\cdot (h\rho {\boldsymbol{v}})=\displaystyle \frac{{R}_{d}}{{C}_{V}}{\boldsymbol{v}}\cdot {\rm{\nabla }}P+\displaystyle \frac{{R}_{d}}{{C}_{V}}{q}_{\mathrm{heat}},\end{eqnarray} \tag{ 24 }$$

where R_d is the specific gas constant and C_V is the heat capacity at constant volume.

Now, we designate "slow" and "fast" quantities and discretize as we did for Equations (12)–(15). Denoting the small time step τ and the large time step t, Equation (24) can be written as

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{P}^{\star [\tau +{\rm{\Delta }}\tau ]}-{P}^{\star [\tau ]}}{{\rm{\Delta }}\tau }+\displaystyle \frac{{R}_{d}}{{C}_{V}}\left[{{\rm{\nabla }}}_{h}\cdot {h}^{[t]}{\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{[\tau +{\rm{\Delta }}\tau ]}\right.\\ \quad \ +\left.\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}{h}^{[t]}{\left(\rho {v}_{r}\right)}^{[\tau +{\rm{\Delta }}\tau ]}\right)\right]=\displaystyle \frac{{R}_{d}}{{C}_{V}}\left[\displaystyle \frac{{\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{[\tau +{\rm{\Delta }}\tau ]}}{{\rho }^{[t]}}\cdot {{\rm{\nabla }}}_{h}{P}^{[t]}\right.\\ \quad \ +\ \left.\displaystyle \frac{{\left(\rho {v}_{r}\right)}^{[\tau +{\rm{\Delta }}\tau ]}}{{\rho }^{[t]}}\displaystyle \frac{\partial }{\partial r}{P}^{[t]}+{q}_{\mathrm{heat}}+{{ \mathcal F }}_{P}^{[t]}\right].\end{array}\end{eqnarray} \tag{ 25 }$$

In this equation we also include the hyperdiffusive pressure flux, ${{ \mathcal F }}_{P}^{[t]}$.

Using the fact that ${\left(\rho {v}_{r}\right)}^{[\tau +{\rm{\Delta }}\tau ]}={\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}+{\left(\rho {v}_{r}\right)}^{[t]}$ (from Equation (9)), Equation (25) becomes

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{P}^{\star [\tau +{\rm{\Delta }}\tau ]}-{P}^{\star [\tau ]}}{{\rm{\Delta }}\tau }+\displaystyle \frac{{R}_{d}}{{C}_{V}}\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}{h}^{[t]}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}\right)\\ \quad \ -\ \displaystyle \frac{{R}_{d}}{{C}_{V}}\displaystyle \frac{{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}}{{\rho }^{[t]}}\displaystyle \frac{\partial {P}^{[t]}}{\partial r}=\displaystyle \frac{{R}_{d}}{{C}_{V}}\left[-{{\rm{\nabla }}}_{h}\cdot \left({h}^{[t]}{\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{[\tau +{\rm{\Delta }}\tau ]}\right)\right.\\ \quad \ -\ \displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}{h}^{[t]}{\left(\rho {v}_{r}\right)}^{[t]}\right)+\displaystyle \frac{{\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{[t]}}{{\rho }^{[t]}}\cdot {{\rm{\nabla }}}_{h}{P}^{[t]}\\ \quad +\ \left.\displaystyle \frac{{\left(\rho {v}_{r}\right)}^{[t]}}{{\rho }^{[t]}}\displaystyle \frac{\partial {P}^{[t]}}{\partial r}+{q}_{\mathrm{heat}}+{{ \mathcal F }}_{P}^{[t]}\right],\end{array}\end{eqnarray} \tag{ 26 }$$

where we have collected horizontal velocity terms and large time step (t superscripted) terms on the right-hand side. Following Tomita & Satoh (2004), we evaluate the pressure gradient force (third term on the right-hand side) and the buoyancy force (fourth term on the right-hand side) at the large time step. This is not the optimal choice for the conservation of energy (Satoh 2013), but the resulting error is small compared to the errors introduced by the diffusion schemes, and it allows the resulting Helmholtz equation (Equation (35)) to be solved implicitly. To achieve a more concise form (and stay consistent in notation with Mendonça et al. 2016), we will introduce the effective gravity,

$$\begin{eqnarray}&&{\tilde{g}}^{[t]}=-\displaystyle \frac{1}{{\rho }^{[t]}}\displaystyle \frac{\partial {P}^{[t]}}{\partial r},\end{eqnarray} \tag{ 27 }$$

and define

$$\begin{eqnarray}\begin{array}{rcl}{S}_{P} & \equiv & \displaystyle \frac{{R}_{d}}{{C}_{V}}\left[-{{\rm{\nabla }}}_{h}\cdot \left({h}^{[t]}{\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{[\tau +{\rm{\Delta }}\tau ]}\right)-\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}{h}^{[t]}{\left(\rho {v}_{r}\right)}^{[t]}\right)\right.\\ & & +\ \left.\displaystyle \frac{{\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{[t]}}{{\rho }^{[t]}}\cdot {{\rm{\nabla }}}_{h}{P}^{[t]}-{\tilde{g}}^{[t]}{\left(\rho {v}_{r}\right)}^{[t]}+{q}_{\mathrm{heat}}+{{ \mathcal F }}_{P}^{[t]}\right].\end{array}\end{eqnarray} \tag{ 28 }$$

Again, note that there are several typos in Equations (39) and (40) of Mendonça et al. (2016), which we have corrected in our Equations (27) and (28).

Then, Equation (26) becomes

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{{P}^{\star [\tau +{\rm{\Delta }}\tau ]}-{P}^{\star [\tau ]}}{{\rm{\Delta }}\tau }+\displaystyle \frac{{R}_{d}}{{C}_{V}}\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}\left({r}^{2}{h}^{[t]}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}\right)\\ \quad \ +\ \displaystyle \frac{{R}_{d}}{{C}_{V}}{\tilde{g}}^{[t]}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}={S}_{P}.\end{array}\end{eqnarray} \tag{ 29 }$$

Next, we must solve for ${P}^{\star [\tau +{\rm{\Delta }}\tau ]}$ and take its derivative in r to replace the pressure term in our vertical momentum equation. Noting that

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial r}\left[\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}{r}^{2}\right]=\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {r}^{2}}{r}^{2}-\displaystyle \frac{2}{{r}^{3}}\displaystyle \frac{\partial }{\partial r}{r}^{2},\end{eqnarray} \tag{ 30 }$$

and taking the derivative with respect to r, we have

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial r}{P}^{\star [\tau +{\rm{\Delta }}\tau ]}=\displaystyle \frac{\partial }{\partial r}{P}^{\star [\tau ]}-{\rm{\Delta }}\tau \displaystyle \frac{{R}_{d}}{{C}_{V}}\\ \quad \ \times \ \left[\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {r}^{2}}{r}^{2}{h}^{[t]}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}-\displaystyle \frac{2}{{r}^{3}}\displaystyle \frac{\partial }{\partial r}{r}^{2}{h}^{[t]}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}\right]\\ \ \ -\ {\rm{\Delta }}\tau \displaystyle \frac{{R}_{d}}{{C}_{V}}\displaystyle \frac{\partial }{\partial r}{\tilde{g}}^{[t]}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}+{\rm{\Delta }}\tau \displaystyle \frac{\partial }{\partial r}{S}_{P}.\end{array}\end{eqnarray} \tag{ 31 }$$

The vertical derivative is taken here in order to eliminate the pressure at time [τ + Δτ] from the vertical momentum equation, as shown below.

2.3.2. Preparing the Continuity Equation

We begin with the discretized form of the continuity Equation (1) and define

$$\begin{eqnarray}\begin{array}{rcl}{S}_{\rho } & \equiv & -{{\rm{\nabla }}}_{h}\cdot {\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}-{{\rm{\nabla }}}_{h}\cdot {\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{[t]}\\ & & -\ \displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}{r}^{2}{\left(\rho {v}_{r}\right)}^{[t]}+{{ \mathcal F }}_{\rho }^{[t]},\end{array}\end{eqnarray} \tag{ 32 }$$

which, like S_P, incorporates the terms evaluated at the large time step and the horizontal momentum term (which is already evaluated for the current small time step).

To eliminate the density at time [τ + Δτ] from the vertical momentum equation, we simply solve for ρ^⋆[τ+Δτ]:

$$\begin{eqnarray}&&{\rho }^{\star [\tau +{\rm{\Delta }}\tau ]}={\rho }^{\star [\tau ]}-{\rm{\Delta }}\tau \displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}{r}^{2}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}+{\rm{\Delta }}\tau {S}_{\rho }.\end{eqnarray} \tag{ 33 }$$

2.3.3. Solving the Vertical Momentum Equation

For the vertical momentum Equation (14), we again collect the terms evaluated at the large time step into a single term defined as

$$\begin{eqnarray}&&{S}_{{v}_{r}}\equiv -\displaystyle \frac{\partial }{\partial r}{P}^{[t]}-{\rho }^{[t]}g-{{ \mathcal A }}_{r}^{[t]}-{{ \mathcal C }}_{r}^{[t]}+{{ \mathcal F }}_{{v}_{r}}^{[t]},\end{eqnarray} \tag{ 34 }$$

where ${{ \mathcal C }}_{r}$ is the vertical component of the Coriolis acceleration and ${{ \mathcal A }}_{r}$ is the vertical component of the advection term ${\rm{\nabla }}\cdot (\rho v\otimes v)$. We now substitute Equations (31) and (33) into Equation (14) and arrange on the left-hand side all terms involving ${\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}$ (the deviation of the vertical momentum from the large time step [t] at time [τ + Δτ]). To achieve the final desired form, we divide through by ΔτR_d/C_V, resulting in

$$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {r}^{2}}{r}^{2}{h}^{[t]}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}+\displaystyle \frac{2}{{r}^{3}}\displaystyle \frac{\partial }{\partial r}{r}^{2}{h}^{[t]}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}\\ \ \ -\ \displaystyle \frac{\partial }{\partial r}{\tilde{g}}^{[t]}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}-\displaystyle \frac{{C}_{V}}{{R}_{d}}\displaystyle \frac{g}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}{r}^{2}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}\\ \ \ +\ \displaystyle \frac{{C}_{V}}{{R}_{d}{\rm{\Delta }}{\tau }^{2}}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}={C}_{0},\end{array}\end{eqnarray} \tag{ 35 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{C}_{0} & = & \displaystyle \frac{{C}_{V}}{{R}_{d}}\left[\displaystyle \frac{1}{{\rm{\Delta }}{\tau }^{2}}{\left(\rho {v}_{r}\right)}^{\star [\tau ]}+\displaystyle \frac{1}{{\rm{\Delta }}\tau }\left({S}_{{v}_{r}}-\displaystyle \frac{\partial }{\partial r}{P}^{\star [\tau ]}-{\rho }^{\star [\tau ]}g\right)\right.\\ & & -\ \left.\displaystyle \frac{\partial }{\partial r}{S}_{P}-{S}_{\rho }g\right].\end{array}\end{eqnarray} \tag{ 36 }$$

Equation (35) is the final Helmholtz equation for the vertical momentum. Note that the version printed in Mendonça et al. (2016) contains several typos that have been corrected here.

2.3.4. The Hydrostatic and Shallow Approximations

The hydrostatic approximation that is commonly used in atmospheric models is derived by assuming that the vertical pressure gradient and gravitational force are much larger than the other terms in the vertical momentum equation, in other words,

$$\begin{eqnarray}&&\left(\left|\displaystyle \frac{\partial }{\partial r}P\right|,\left|\rho g\right|\right)\gg \left(\left|\displaystyle \frac{\partial \rho {v}_{r}}{\partial t}\right|,\left|{{ \mathcal A }}_{r}\right|,\left|{{ \mathcal C }}_{r}\right|\right).\end{eqnarray} \tag{ 37 }$$

Making this assumption results in the well-known equation for hydrostatic equilibrium:

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial r}P=-\rho g.\end{eqnarray} \tag{ 38 }$$

In the work of Tomita & Satoh (2004), this assumption was made possible by maintaining a factor α next to all the terms in the vertical momentum equation except the two above (the pressure gradient and gravitational force), which could then be set to unity for nonhydrostatic or zero for hydrostatic, leaving the algorithm otherwise unchanged. This should be used with caution, however, as White et al. (2005) found that the application of hydrostatic balance solely to the vertical momentum equation produces an inconsistency in the representation of potential vorticity within the model unless the "shallow" approximation is also applied. Because of the mixture of coordinate systems in the Tomita & Satoh (2004) algorithm (Cartesian for the horizontal momentum and spherical for the vertical) and the form of the differential operators on the icosahedral grid, it is not simple to find exact correspondence between equations used in THOR and those studied in White et al. (2005). However, because the approximations discussed here refer implicitly to a spherical coordinate system, we believe that the problem described in that paper is relevant here. Hence, we follow their example in applying the hydrostatic and shallow approximations.

The first approximation, which we will call "quasi-hydrostatic" following the terminology used by White et al. (2005), involves neglecting the material derivative of the vertical velocity, in other words,

$$\begin{eqnarray}&&\displaystyle \frac{{{Dv}}_{r}}{{Dt}}=0.\end{eqnarray} \tag{ 39 }$$

In practice, this involves not only the first term (the partial derivative in time) of Equation (14) but also the advection term, ${{ \mathcal A }}_{r}$. The method used to calculate ${ \mathcal A }$ in Cartesian coordinates from the vector momentum retains the effects of curvature (which result in the well-known "metric" terms in a fully spherical coordinate system); thus, we cannot merely discard ${{ \mathcal A }}_{r}$ from Equation (14). However, the metric part of ${{ \mathcal A }}_{r}$ can be simply calculated from the horizontal momentum, as

$$\begin{eqnarray}&&{{ \mathcal A }}_{r}^{\mathrm{QH}}=\displaystyle \frac{\rho {{\boldsymbol{v}}}_{h}\cdot {{\boldsymbol{v}}}_{h}}{r},\end{eqnarray} \tag{ 40 }$$

where v_h is the horizontal momentum vector in Cartesian coordinates. Then, under the quasi-hydrostatic, deep (QHD) assumption, the vertical momentum equation becomes

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial }{\partial r}{P}^{\star [\tau +{\rm{\Delta }}\tau ]}+{\rho }^{\star [\tau +{\rm{\Delta }}\tau ]}g=-\displaystyle \frac{\partial }{\partial r}{P}^{[t]}\\ \quad \ -\ {\rho }^{[t]}g-{\left({{ \mathcal A }}_{r}^{\mathrm{QH}}\right)}^{[t]}-{{ \mathcal C }}_{r}^{[t]}.\end{array}\end{eqnarray} \tag{ 41 }$$

As before, we use the thermodynamic and continuity equations to substitute the terms on the left-hand side, which results in the modified Helmholtz equation:

$$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {r}^{2}}{r}^{2}{h}^{[t]}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}+\displaystyle \frac{2}{{r}^{3}}\displaystyle \frac{\partial }{\partial r}{r}^{2}{h}^{[t]}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}\\ \quad -\ \displaystyle \frac{\partial }{\partial r}{\tilde{g}}^{[t]}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}-\displaystyle \frac{{C}_{V}}{{R}_{d}}\displaystyle \frac{g}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}{r}^{2}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}={C}_{0}^{\mathrm{QH}},\end{array}\end{eqnarray} \tag{ 42 }$$

where

$$\begin{eqnarray}&&{C}_{0}^{\mathrm{QH}}=\displaystyle \frac{{C}_{V}}{{R}_{d}}\left[\displaystyle \frac{1}{{\rm{\Delta }}\tau }\left({S}_{{v}_{r}}^{\mathrm{QH}}-\displaystyle \frac{\partial }{\partial r}{P}^{\star [\tau ]}-{\rho }^{\star [\tau ]}g\right)-\displaystyle \frac{\partial }{\partial r}{S}_{P}-{S}_{\rho }g\right],\end{eqnarray} \tag{ 43 }$$

where ${S}_{{v}_{r}}^{\mathrm{QH}}$ is computed using the QH version of the advection term, ${{ \mathcal A }}_{r}^{\mathrm{QH}}$. Note that the hyperdiffusive flux term ${{ \mathcal F }}_{{v}_{r}}$ is also zero in this approximation. We then solve equations in the same fashion as in the nonhydrostatic model. Everywhere else, the equations remain identical to their nonhydrostatic versions.

The second approximation (the shallow approximation) is more involved. All horizontal differential operators are defined at the reference pressure (the bottom of the model), and these are proportional to 1/r₀, where r₀ is the radius of the planet (the distance from the center to the location of the reference pressure). To account for curvature, the horizontal operators are multiplied by r₀/r, where r = r₀ + z and z is the altitude. Simply, this scales the horizontal area of the control volumes with altitude. In the shallow approximation, r = r₀, and the scaling is removed from the horizontal operators. The radial derivative in the divergence also changes, as follows:

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial ({r}^{2}\psi )}{\partial r}\to \displaystyle \frac{1}{{r}_{0}^{2}}\displaystyle \frac{\partial ({r}_{0}^{2}\psi )}{\partial r}=\displaystyle \frac{\partial \psi }{\partial r},\end{eqnarray} \tag{ 44 }$$

where ψ is just a representative quantity. In this way, the second derivative (Equation (30)) becomes

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial r}\left[\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}{r}^{2}\right]\to \displaystyle \frac{\partial }{\partial r}\left[\displaystyle \frac{1}{{r}_{0}^{2}}\displaystyle \frac{\partial }{\partial r}{r}_{0}^{2}\right]=\displaystyle \frac{{\partial }^{2}}{\partial {r}^{2}}.\end{eqnarray} \tag{ 45 }$$

Using this in the thermodynamic equation and again constructing the Helmholtz equation, we have

$$\begin{eqnarray}&&\begin{array}{l}-\displaystyle \frac{{\partial }^{2}}{\partial {r}^{2}}{h}^{[t]}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}-\displaystyle \frac{\partial }{\partial r}{\tilde{g}}^{[t]}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}\\ \ \ -\ \displaystyle \frac{{{gC}}_{V}}{{R}_{d}}\displaystyle \frac{\partial }{\partial r}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}+\displaystyle \frac{{C}_{V}}{{R}_{d}{\rm{\Delta }}{\tau }^{2}}{\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}={C}_{0}^{{\rm{S}}},\end{array}\end{eqnarray} \tag{ 46 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{C}_{0}^{{\rm{S}}} & = & \displaystyle \frac{{C}_{V}}{{R}_{d}}\left[\displaystyle \frac{1}{{\rm{\Delta }}{\tau }^{2}}{\left(\rho {v}_{r}\right)}^{\star [\tau ]}+\displaystyle \frac{1}{{\rm{\Delta }}\tau }\left({S}_{{v}_{r}}^{{\rm{S}}}-\displaystyle \frac{\partial }{\partial r}{P}^{\star [\tau ]}-{\rho }^{\star [\tau ]}g\right)\right.\\ & & -\ \left.\displaystyle \frac{\partial }{\partial r}{S}_{P}^{{\rm{S}}}-{S}_{\rho }^{{\rm{S}}}g\right].\end{array}\end{eqnarray} \tag{ 47 }$$

In this case, ${S}_{{v}_{r}}^{{\rm{S}}}$, ${S}_{P}^{{\rm{S}}}$, ${S}_{\rho }^{{\rm{S}}}$, and all advection terms are calculated using the shallow operators described above. One additional change is made to ensure that the Coriolis acceleration is consistently represented by the horizontal and vertical equations. The Coriolis vector, ${\boldsymbol{ \mathcal C }}$, is now

$$\begin{eqnarray}\begin{array}{rcl}{{\boldsymbol{ \mathcal C }}}^{{\rm{S}}} & = & 2\rho {\boldsymbol{\Omega }}\times {\boldsymbol{v}}=2{\rm{\Omega }}\rho [\left({v}_{3}\sin \phi \cos \phi \sin \lambda -{v}_{2}{\sin }^{2}\phi \right){\hat{{\boldsymbol{e}}}}_{1}\\ & & -\ \left({v}_{3}\sin \phi \cos \phi \cos \lambda -{v}_{1}{\sin }^{2}\phi \right){\hat{{\boldsymbol{e}}}}_{2}\\ & & +\ \left({v}_{2}\sin \phi \cos \phi \cos \lambda -{v}_{2}\sin \phi \cos \phi \sin \lambda \right){\hat{{\boldsymbol{e}}}}_{3}],\end{array}\end{eqnarray} \tag{ 48 }$$

where the unit vectors ${\hat{{\boldsymbol{e}}}}_{{\boldsymbol{i}}}$ represent the rotating Cartesian coordinate system fixed at the center of the planet (the spin axis is aligned with $\hat{{{\boldsymbol{e}}}_{3}}$), and the velocities v_i are the total (horizontal plus vertical) in the corresponding direction. Equation (48) is equivalent to the form of Coriolis that appears in the primitive equations (i.e., the horizontal component of ${\boldsymbol{\Omega }}$ is neglected). The radial component, ${{ \mathcal C }}_{r}$, is now equal to zero.

2.3.5. Discretizing and Solving the 1D Helmholtz Equation

The time-discretized 1D Helmholtz equation (Equation (35)) is now spatially discretized in the following way. The vertical momentum is solved for at the midpoints between layers, in contrast to the horizontal momentum, pressure, and density, which are solved for at the center of the layer. Denoting quantities defined at the center of layers with subscript c and quantities defined at the midpoint between layers (the interfaces) with subscript m, we write the first term in Equation (35) containing the second derivative in r as

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{-1}{{r}^{2}}\displaystyle \frac{{\partial }^{2}}{\partial {r}^{2}}{r}^{2}{h}_{m}{W}_{m} & \approx & \displaystyle \frac{-1}{{\left({r}_{m}^{i}\right)}^{2}}\left[\displaystyle \frac{{\left({r}_{m}^{i+1}\right)}^{2}\,{h}_{m}^{i+1}{W}_{m}^{i+1}-{\left({r}_{m}^{i}\right)}^{2}\,{h}_{m}^{i}{W}_{m}^{i}}{{\rm{\Delta }}{r}_{m}^{i+1}{\rm{\Delta }}{r}_{c}^{i}}\right.\\ & & -\ \left.\displaystyle \frac{{\left({r}_{m}^{i}\right)}^{2}\,{h}_{m}^{i}{W}_{m}^{i}-{\left({r}_{m}^{i-1}\right)}^{2}\,{h}_{m}^{i-1}{W}_{m}^{i-1}}{{\rm{\Delta }}{r}_{m}^{i}{\rm{\Delta }}{r}_{c}^{i}}\right],\end{array}\end{eqnarray} \tag{ 49 }$$

for the ith midpoint (between layers), i.e., at ${r}_{m}^{i}$. Superscripts indicate the layer/midpoint at which each quantity is defined and are ordered such that the ith midpoint is at the bottom of the ith layer. Here, ${h}_{m}^{i}={h}^{[t]}$ and ${W}_{m}^{i}={\left(\rho {v}_{r}\right)}^{\star [\tau +{\rm{\Delta }}\tau ]}$ for the ith midpoint. For shorthand, the spatial separations are defined as

$$\begin{eqnarray}\begin{array}{rcl}{\rm{\Delta }}{r}_{c}^{i} & = & {r}_{c}^{i}-{r}_{c}^{i-1},\\ {\rm{\Delta }}{r}_{m}^{i} & = & {r}_{m}^{i}-{r}_{m}^{i-1},\\ {\rm{\Delta }}{r}_{m}^{i+1} & = & {r}_{m}^{i+1}-{r}_{m}^{i}.\end{array}\end{eqnarray} \tag{ 50 }$$

First derivatives are calculated by performing a first-order finite difference across the layers i and i − 1 and then interpolating to the midpoint. The resulting procedure, for terms 2–4 in Equation (35), is

$$\begin{eqnarray}\begin{array}{rcl}\displaystyle \frac{\partial }{\partial r}{F}_{m}^{i} & \approx & \displaystyle \frac{{r}_{m}^{i}-{r}_{c}^{i-1}}{{\rm{\Delta }}{r}_{m}^{i+1}{\rm{\Delta }}{r}_{c}^{i}}{F}_{m}^{i+1}+\left[\displaystyle \frac{{r}_{c}^{i}-{r}_{m}^{i}}{{\rm{\Delta }}{r}_{m}^{i}{\rm{\Delta }}{r}_{c}^{i}}-\displaystyle \frac{{r}_{m}^{i}-{r}_{c}^{i-1}}{{\rm{\Delta }}{r}_{m}^{i+1}{\rm{\Delta }}{r}_{c}^{i}}\right]{F}_{m}^{i}\\ & & -\ \displaystyle \frac{{r}_{c}^{i}-{r}_{m}^{i}}{{\rm{\Delta }}{r}_{m}^{i}{\rm{\Delta }}{r}_{c}^{i}}{F}_{m}^{i-1},\end{array}\end{eqnarray} \tag{ 51 }$$

where ${F}_{m}^{i}$ represents the various arguments inside the derivatives for the ith midpoint.

The resulting spatially discretized equations for each midpoint i form a system of equations that can be written as

$$\begin{eqnarray}&&{a}_{i}{W}_{m}^{i+1}+{b}_{i}{W}_{m}^{i}+{c}_{i}{W}_{m}^{i-1}={d}_{i}\ \mathrm{for}\ (1\lt i\lt n-1),\end{eqnarray} \tag{ 52 }$$

which we then solve using Thomas's algorithm for tridiagonal matrices and the boundary conditions $({W}_{m}^{0},{W}_{m}^{n})=0$ (n represents the index of the top layer). After applying the discretization process and rearranging terms, the resulting coefficients are

$$\begin{eqnarray}\begin{array}{rcl}{a}_{i} & = & \displaystyle \frac{1}{{\rm{\Delta }}{r}_{m}^{i+1}{\rm{\Delta }}{r}_{c}^{i}}\left[-\displaystyle \frac{{\left({r}_{m}^{i+1}\right)}^{2}\,{h}_{m}^{i+1}}{{\left({r}_{m}^{i}\right)}^{2}}+({r}_{m}^{i}-{r}_{c}^{i-1})\right.\\ & & \times \ \left.\left(\displaystyle \frac{2{\left({r}_{m}^{i+1}\right)}^{2}\,{h}_{m}^{i+1}}{{\left({r}_{m}^{i}\right)}^{3}}-{\tilde{g}}_{m}^{i+1}-\displaystyle \frac{{C}_{V}g{\left({r}_{m}^{i+1}\right)}^{2}}{{R}_{d}{\left({r}_{m}^{i}\right)}^{2}}\right)\right],\\ {b}_{i} & = & \displaystyle \frac{1}{{\rm{\Delta }}{r}_{m}^{i+1}{\rm{\Delta }}{r}_{c}^{i}}\left[{h}_{m}^{i}-({r}_{m}^{i}-{r}_{c}^{i-1})\left(\displaystyle \frac{2{h}_{m}^{i}}{{r}_{m}^{i}}-{\tilde{g}}_{m}^{i}-\displaystyle \frac{{C}_{V}g}{{R}_{d}}\right)\right]\\ & & +\ \displaystyle \frac{1}{{\rm{\Delta }}{r}_{m}^{i}{\rm{\Delta }}{r}_{c}^{i}}\left[{h}_{m}^{i}-({r}_{m}^{i}-{r}_{c}^{i})\left(\displaystyle \frac{2{h}_{m}^{i}}{{r}_{m}^{i}}-{\tilde{g}}_{m}^{i}-\displaystyle \frac{{C}_{V}g}{{R}_{d}}\right)\right]\\ & & +\ \displaystyle \frac{{C}_{V}}{{R}_{d}{\rm{\Delta }}{\tau }^{2}},\\ {c}_{i} & = & \displaystyle \frac{1}{{\rm{\Delta }}{r}_{m}^{i}{\rm{\Delta }}{r}_{c}^{i}}\left[-\displaystyle \frac{{\left({r}_{m}^{i-1}\right)}^{2}\,{h}_{m}^{i-1}}{{\left({r}_{m}^{i}\right)}^{2}}+({r}_{m}^{i}-{r}_{c}^{i})\right.\\ & & \times \ \left.\left(\displaystyle \frac{2{\left({r}_{m}^{i-1}\right)}^{2}\,{h}_{m}^{i-1}}{{\left({r}_{m}^{i}\right)}^{3}}-{\tilde{g}}_{m}^{i-1}-\displaystyle \frac{{C}_{V}g{\left({r}_{m}^{i-1}\right)}^{2}}{{R}_{d}{\left({r}_{m}^{i}\right)}^{2}}\right)\right],\\ {d}_{i} & = & {C}_{0}^{i}.\end{array}\end{eqnarray} \tag{ 53 }$$

As before, the enthalpy is defined on the large time step at the ith midpoint, i.e., ${h}_{m}^{i}={h}^{[t]}$, and likewise for ${\tilde{g}}_{m}^{i}$ and ${C}_{0}^{i}$.

Of course, in the case that the height grid is uniformly spaced, ${\rm{\Delta }}{r}_{c}^{i}={\rm{\Delta }}{r}_{m}^{i}={\rm{\Delta }}{r}_{m}^{i+1}$ and the above equations can be greatly simplified. The current version of THOR utilizes only a uniform grid; however, the model has been coded according to the above equations so that nonuniform grids can be utilized in the future. A further simplification can be made in the case of the shallow approximation, in which ${r}_{m}^{i+1}\approx {r}_{m}^{i}\approx {r}_{m}^{i-1}$. This simplification is carried out in the model when the shallow approximation is used.

The flux-form hyperdiffusion terms that are applied to Equations (12)–(14) and (18) are

$$\begin{eqnarray}&&{{ \mathcal F }}_{\rho }=-{{\rm{\nabla }}}_{h}^{2}{K}_{\mathrm{hyp}}{{\rm{\nabla }}}_{h}^{2}\rho ,\end{eqnarray} \tag{ 54 }$$

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal F }}}_{{{\boldsymbol{v}}}_{h}}=-{{\rm{\nabla }}}_{h}^{2}\rho {K}_{\mathrm{hyp}}{{\rm{\nabla }}}_{h}^{2}{{\boldsymbol{v}}}_{h},\end{eqnarray} \tag{ 55 }$$

$$\begin{eqnarray}&&{{ \mathcal F }}_{{v}_{r}}=-{{\rm{\nabla }}}_{h}^{2}\rho {K}_{\mathrm{hyp}}{{\rm{\nabla }}}_{h}^{2}{v}_{r},\end{eqnarray} \tag{ 56 }$$

$$\begin{eqnarray}&&{{ \mathcal F }}_{P}=-{R}_{d}{{\rm{\nabla }}}_{h}^{2}\rho {K}_{\mathrm{hyp}}{{\rm{\nabla }}}_{h}^{2}T.\end{eqnarray} \tag{ 57 }$$

The divergence damping term in the horizontal momentum equation is

$$\begin{eqnarray}&&{{\boldsymbol{ \mathcal G }}}_{{{\boldsymbol{v}}}_{h}}=-{K}_{\mathrm{div}}{{\rm{\nabla }}}_{h}{{\rm{\nabla }}}_{h}^{2}\left({{\rm{\nabla }}}_{h}\cdot (\rho {{\boldsymbol{v}}}_{h})+\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{\partial }{\partial r}(\rho {v}_{r}{r}^{2})\right).\end{eqnarray} \tag{ 58 }$$

Note that the order of the gradient and Laplacian operators was incorrectly reversed in Mendonça et al. (2016); the model is coded as written in our Equation (58). Divergence damping is necessary to eliminate noise produced by the time-splitting integration scheme (Skamarock & Klemp 1992; Mendonça et al. 2016).

The diffusion coefficients, K_hyp and K_div, have the same functional dependence on the grid resolution and time step size but can be individually adjusted. These are

$$\begin{eqnarray}&&{K}_{{\rm{hyp}}}={D}_{{\rm{hyp}}}\displaystyle \frac{{\bar{d}}^{4}}{{\rm{\Delta }}t}\end{eqnarray} \tag{ 59 }$$

$$\begin{eqnarray}&&{K}_{{\rm{div}}}={D}_{{\rm{div}}}\displaystyle \frac{{\bar{d}}^{4}}{{\rm{\Delta }}t},\end{eqnarray} \tag{ 60 }$$

where $\bar{d}$ is the average width of the control volumes given by Equation (7). The hyperdiffusion fluxes are updated on the large time step, while the divergence damping fluxes are updated every small time step.

The boundary conditions for the top and bottom of the model are that the vertical velocity must equal zero; this is the simplest assumption that allows for conservation of energy and axial angular momentum (Staniforth & Wood 2003). Unfortunately, this causes the boundary to act as a node for vertically propagating waves, causing them to reflect and potentially amplify. An additional dissipation mechanism is often needed to eliminate these reflecting waves. In a real atmosphere, vertical propagating waves will break in the upper layers; however, in the model the artificial reflection of these waves becomes an additional source of noise that can trigger numerical instabilities and cause the model to crash. Methods used to damp reflecting waves at the boundaries are often called "sponge layers." In THOR, we use the method based on Skamarock & Klemp (2008) and described in Mendonça et al. (2018b), which we briefly reiterate here.

The zonal, meridional, and vertical winds are all damped toward their zonal averages using Rayleigh friction. In principle, damping toward the zonal averages, rather than toward zero, will selectively damp waves while allowing the general flow to persist. In practice, the method is imperfect and the effect of the sponge layer on the flow at the highest layers is discernible as a decrease in the zonal wind speed. For this reason, we minimize the strength and size of the sponge layer in our simulations.

The damping takes the form

$$\begin{eqnarray}&&\displaystyle \frac{d{\boldsymbol{v}}}{{dt}}=-{k}_{\mathrm{sp}}(\eta )\left({\boldsymbol{v}}-\bar{{\boldsymbol{v}}}\right),\end{eqnarray} \tag{ 61 }$$

where v represents the vector velocity and $\bar{{\boldsymbol{v}}}$ represents the zonal mean of the components, η = z/z_top is the fractional altitude, and k_sp(η) is the damping strength as a function of the fractional altitude (with units of s⁻¹).

The damping strength is a function of altitude and takes the form

$$\begin{eqnarray}{k}_{\mathrm{sp}}(\eta )=\left\{\begin{array}{ll}0 & \eta \lt {\eta }_{\mathrm{sp}},\\ {k}_{\mathrm{sp}}({\eta }_{\mathrm{sp}}){\sin }^{2}\left(\tfrac{\pi }{2}\tfrac{\eta -{\eta }_{\mathrm{sp}}}{1-{\eta }_{\mathrm{sp}}}\right) & \eta \gt {\eta }_{\mathrm{sp}},\end{array}\right.\end{eqnarray} \tag{ 62 }$$

where η_sp, the fractional height at which the sponge layer begins, and k_sp(η_sp) are values set by the user.

To calculate the zonal mean on our icosahedral grid, we divide the sphere into a user-set number of latitude bins, n_lat, and compute the average within each bin. This is then considered as the average at the center latitude of each bin. At any given latitude, the zonal mean $\bar{{\boldsymbol{v}}}$ is determined by linearly interpolating from the center of the respective latitude bin. This allows the zonal-mean velocities to vary more smoothly with latitude.

Equation (61) is calculated for the zonal, meridional, and vertical wind speeds. Since version 2.3 of the code, there are several options for the time integration of the sponge layer friction. In "implicit mode," which is used in the simulations here, the calculation of the rate of change in wind speed is done in the ProfX step (see Section 2.5), and the velocities are updated directly, implicitly, during this step. In "explicit mode," the rate of change is converted into fluxes that are passed to the dynamical core. For the vertical winds, this flux takes the form ρ dv_r/dt (the vertical component of Equation (61)), which is added to Equation (56). For the horizontal winds, we first compute the zonal and meridional components of ρ dv_h/dt, convert those to Cartesian coordinates, and then add them to Equation (55).

For clarity, we outline and summarize the flow of integration. At the top level, each time step contains two components: the dynamical core (THOR) and physics modules (ProfX); see Figure 1. The "dynamical core" refers to the solution of the Euler equations, and "physics modules" refers to any additional processes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The integration of the dynamical core proceeds as follows. For any prognostic quantity, Φ, the value at time t + Δt is calculated using three large time steps of variable length (a third-order Runge–Kutta scheme; Wicker & Skamarock 2002; Tomita & Satoh 2004; Klemp et al. 2007; Mendonça et al. 2016). Each large time step is broken into smaller time steps of variable length: the first large time step consists of one small time step of length Δτ = Δt/3, the second large step consists of n_max/2 small steps of length Δτ = Δt/n_max, and the third large step consists of n_max steps of length Δτ = Δt/n_max.

A sketch of one time step for variable Φ follows, for n_max = 6 (see Figure 2). The beginning of the time step is time t. First, we calculate the slow terms in the derivative ∂Φ/∂t; these are the terms in Equations (12)–(15) and (18) designated with ^[t]. Next, we calculate the fast terms (designated with ${}^{[\tau ]}$ or ${}^{[\tau +{\rm{\Delta }}\tau ]}$) at times t and t + Δτ = t + Δt/3. The deviation, Φ^⋆, is defined with respect to time t; thus, Φ^⋆(t) = 0. The fast and slow terms are used to calculate Φ^⋆(t + Δt/3), the deviation of Φ at time t + Δt/3. Then, we compute Φ(t + Δt/3) = Φ(t) + Φ^⋆(t + Δt/3). This completes the first large step.

The second large step begins. First, we recompute the deviations, which are now defined with respect to t + Δt/3, so that Φ^⋆(t) = Φ(t) − Φ(t + Δt/3) $\ne$ 0. The slow terms (superscript ${}^{[t]}$) are recomputed at time t + Δt/3. The fast terms are recomputed at times t and t + Δτ = t + Δt/6 (again, n_max = 6 in this example). The fast terms at (t, t + Δt/6) and slow terms at t + Δt/3 are used to advance the deviations by a time step Δt/6, resulting in the value Φ^⋆(t + Δt/6). Fast terms are recomputed at t + Δt/6 and t + 2Δt/6 = t + Δt/3, and these are used with the slow terms to advance one more step of size Δt/6, resulting in Φ^⋆(t + Δt/3). We recompute the fast terms a third time and use them to compute Φ^⋆(t + Δt/2), and finally Φ(t + Δt/2) = Φ(t) + Φ^⋆(t + Δt/2). This completes the second large step.

The third large step begins. Again, the deviations are recalculated, this time with respect to t + Δt/2. Slow terms are recalculated at t + Δt/2, fast terms at t and t + Δt/6. These are used to calculate the deviation Φ^⋆(t + Δt/6). The fast terms are then recomputed and the deviation updated for a total of n_max = 6 times. We then have Φ^⋆(t + Δt), from which we calculate Φ(t + Δt) = Φ(t) + Φ^⋆(t + Δt). This completes the Runge–Kutta loop.

A more detailed outline of a single time step is as follows:

• 1.  
  ProfX step (additional physics)

  • (a)  
    Compute benchmark forcing if applicable. Typically, for benchmark tests, the prognostic variables are updated implicitly or explicitly during this step, rather than computing fluxes that are included in step 2.
  • (b)  
    Compute radiative transfer fluxes (Section 3.3). These are passed to the dynamical core as q_heat, rather than updating thermodynamic variables directly during this step.
  • (c)  
    Update sponge layer quantities (zonal-mean winds and resulting drag, Equation (61)). These can be used to implicitly update the wind speed during this step (implicit mode) or passed as fluxes to the dynamical core and added to the hyperdiffusive terms ${{\boldsymbol{ \mathcal F }}}_{{{\boldsymbol{v}}}_{h}}$ and ${{ \mathcal F }}_{{v}_{r}}$ (explicit mode).
• 2.  
  THOR step (dynamical core): solving fluid equations (Equations (12)–(15) and (18))

  • (a)  
    Begin large time step: three steps total, where the first step advances the prognostic variables to t + Δt/3, the second to t + Δt/2, and the third to t + Δt.

    • i.  
      Compute advection and Coriolis terms ${{\boldsymbol{ \mathcal A }}}_{h}$, ${{ \mathcal A }}_{r}$, ${{\boldsymbol{ \mathcal C }}}_{h}$, ${{ \mathcal C }}_{r}$ at time t, t + Δt/3, or t + Δt/2 (for the first, second, and third steps).
    • ii.  
      Compute enthalpy h (Equation (21)), effective gravity $\tilde{g}$ (Equation (27)), and potential temperature θ (Equation (4)) at time t, t + Δt/3, or t + Δt/2.
    • iii.  
      Compute hyperdiffusive and divergence damping fluxes ${{ \mathcal F }}_{\rho }$, ${{\boldsymbol{ \mathcal F }}}_{{{\boldsymbol{v}}}_{h}}$, ${{ \mathcal F }}_{{v}_{r}}$, ${{ \mathcal F }}_{P}$, ${{\boldsymbol{ \mathcal G }}}_{{{\boldsymbol{v}}}_{h}}$ (Equations (54)–(58)) at time t, t + Δt/3, or t + Δt/2. Add sponge layer drag (Equation (61)), if explicit mode is used for sponge layer.
    • iv.  
      Compute slow modes: sums of [t] terms (steps 2(a)i–2(a)iii) in Equations (12)–(18), including ${{\boldsymbol{ \mathcal G }}}_{{{\boldsymbol{v}}}_{h}}$ and RT fluxes (as q_heat) from ProfX step.
    • v.  
      Second and third large time steps only: update deviations ρ^⋆, ${\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{\star }$, ${\left(\rho {v}_{r}\right)}^{\star }$, P^⋆ (Equation (8)). Deviations are equal to zero on the first step.
    • vi.  
      Begin small time step: n_step steps, where n_step = (1, n_max/2, n_max) and Δτ = (Δt/3, Δt/n_max, Δt/n_max) for the first, second, and third large time step, respectively. For the nth iteration, the current time is τ = t + nΔτ.

      • A.  
        Update divergence damping ${{\boldsymbol{ \mathcal G }}}_{{{\boldsymbol{v}}}_{h}}$ (Equation (58)) at time τ.
      • B.  
        Compute horizontal momentum deviation ${\left(\rho {{\boldsymbol{v}}}_{h}\right)}^{\star }$ (Equation (13)) at time τ + Δτ.
      • C.  
        Compute S_P, S_ρ, ${S}_{{v}_{r}}$ (Equations (28), (32), and (34)). These terms encapsulate the slow modes plus the horizontal momentum deviations.
      • D.  
        Compute vertical momentum deviation ${\left(\rho {v}_{r}\right)}^{\star }$ (Equation (35)) using Thomas's algorithm described in Section 2.3.5, at time τ + Δτ.
      • E.  
        Compute density deviations ρ^⋆ (Equation (12)) at time τ + Δτ.
      • F.  
        Compute potential temperature density (ρθ) (Equation (15)) at time τ + Δτ.
      • G.  
        Compute pressure deviation P^⋆ (Equation (18)) at time τ + Δτ.
    • vii.  
      End small time step.
    • viii.  
      Update prognostic variables ρ, (ρ v_h), (ρv_r), and P using final deviations from small time step loop. These are now defined at times t + Δt/3, t + Δt/2, or t + Δt, for the first, second, or third large step, respectively.
  • (b)  
    End large time step.

By default, the model now outputs additional diagnostic quantities: total energy, mass, angular momentum, and entropy. The mass, angular momentum, total energy (kinetic + internal + potential), and entropy of the ith grid point and jth vertical level are, respectively,

$$\begin{eqnarray}\begin{array}{rcl}{m}_{{ij}} & = & {\rho }_{{ij}}{V}_{{ij}}\\ {{\boldsymbol{l}}}_{{ij}} & = & {\rho }_{{ij}}{{\boldsymbol{r}}}_{{ij}}\times ({{\boldsymbol{v}}}_{h}^{{ij}}+{\boldsymbol{\Omega }}\times {{\boldsymbol{r}}}_{{ij}}){V}_{{ij}}\\ {E}_{{ij}}^{\mathrm{tot}} & = & {\rho }_{{ij}}\left(\displaystyle \frac{1}{2}{{\boldsymbol{v}}}_{{ij}}^{2}+{C}_{V}{T}_{{ij}}+{{gz}}_{j}\right){V}_{{ij}}\\ {s}_{{ij}} & = & {\rho }_{{ij}}{C}_{P}\mathrm{log}({\theta }_{{ij}}){V}_{{ij}}.\end{array}\end{eqnarray} \tag{ 63 }$$

Above, ρ is the density, V is the volume of the control volume, r is the Cartesian vector position on the sphere, ${\boldsymbol{\Omega }}$ is the vector rotation rate, v_h is the horizontal wind vector, v is the total wind vector (horizontal + vertical), C_V is the specific heat at constant volume, T is the temperature, g is the gravity (assumed to be constant), z is the altitude, C_P is the specific heat at constant pressure, and θ is the potential temperature. Vectors are defined in a Cartesian coordinate system centered on the center of the planet and rotating about the ${\hat{{\boldsymbol{e}}}}_{3}$-axis with rotation rate Ω. The vertical momentum can be ignored in the calculation of l_ij, as it is parallel to r_ij by definition. When integrated over the sphere, the nonaxial (${\hat{{\boldsymbol{e}}}}_{1}$ and ${\hat{{\boldsymbol{e}}}}_{2}$) components of the angular momentum should vanish; in practice, they are not identically zero because of numerical noise, so these components can provide a useful test of numerical accuracy.

The global total of each quantity is calculated by summing over i grid points and j vertical levels. In the deep model, the volume of the ith grid point and jth vertical level is

$$\begin{eqnarray}&&{V}_{{ij}}=\displaystyle \frac{{A}_{i}}{3{r}_{0}^{2}}\left[{\left({r}_{m}^{j+1}\right)}^{3}-{\left({r}_{m}^{j}\right)}^{3}\right],\end{eqnarray} \tag{ 64 }$$

where A_i is the area of the ith control volume at the lowest boundary of the model, r₀ is the planet radius, and ${r}_{m}^{j+1}$ and ${r}_{m}^{j}$ are the radial coordinates of the top and bottom boundaries of the jth layer. Figure 3 illustrates the calculation of the size of the control volume, V_ij. The area, A_i, of each control volume is calculated during the grid construction (and adjusted by the spring dynamics process). This is calculated by decomposing the hexagon- or pentagon-shaped control volume into triangles formed by the center and two adjacent vertices (see Figure 2 of Mendonça et al. 2016) and summing the areas of the triangles. The areas of each triangle are calculated using the formula for spherical triangles ${A}_{\mathrm{tri}}=(\angle a+\angle b+\angle c-\pi ){r}_{0}^{2}$, where the angles $\angle a$, $\angle b$, and $\angle c$ are calculated from the vector locations of the center and corresponding vertices. When the shallow approximation is used, the control volumes are treated as flat, and so the volume is

$$\begin{eqnarray}&&{V}_{{ij}}={A}_{i}\left({r}_{m}^{j+1}-{r}_{m}^{j}\right).\end{eqnarray} \tag{ 65 }$$

Zoom In Zoom Out Reset image size

Here we provide a description of the dry convective adjustment scheme, for completeness. Nonhydrostatic models require a parameterization for convection at the coarse resolutions we typically use for exoplanet modeling; typically, scales less than tens of kilometers must be resolved to capture convection with no parameterization in Earth simulations (Arakawa 2004; Jung & Arakawa 2004; Rio et al. 2019). This scheme is based on Hourdin et al. (1993). After the dynamical code time step, but prior to the calculation of radiative transfer/Newtonian cooling, each vertical column is searched for unstable layers. Static stability is given by the condition

$$\begin{eqnarray}&&\displaystyle \frac{\partial \theta }{\partial r}\geqslant 0,\end{eqnarray} \tag{ 66 }$$

where θ is the potential temperature. When an unstable layer is detected (meaning Equation (66) is violated), we define a "mixed" potential temperature, θ_mixed, equal to the average potential temperature in the layer. First, we integrate upward through the unstable layer to find the enthalpy h:

$$\begin{eqnarray}&&h={C}_{P}{\int }_{{P}_{B}}^{{P}_{T}}\theta {\left(\displaystyle \frac{P}{{P}_{0}}\right)}^{{R}_{d}/{C}_{P}}{dP},\end{eqnarray} \tag{ 67 }$$

where P₀, P_B, and P_T are the pressure at bottom of the entire column, the pressure at the bottom of the unstable layer, and the pressure at the top of the unstable layer, respectively. Then, the potential temperature across the entire layer is set equal to θ_mixed, given by

$$\begin{eqnarray}&&{\theta }_{\mathrm{mixed}}=\displaystyle \frac{h}{{C}_{P}{\int }_{{P}_{B}}^{{P}_{T}}{\left(\tfrac{P}{{P}_{0}}\right)}^{{R}_{d}/{C}_{P}}{dP}}=\displaystyle \frac{{\int }_{{P}_{B}}^{{P}_{T}}\theta {\left(\tfrac{P}{{P}_{0}}\right)}^{{R}_{d}/{C}_{P}}{dP}}{{\int }_{{P}_{B}}^{{P}_{T}}{\left(\tfrac{P}{{P}_{0}}\right)}^{{R}_{d}/{C}_{P}}{dP}},\end{eqnarray} \tag{ 68 }$$

which is effectively like mixing the entropy across the unstable layer and enforcing an adiabatic profile in that region. After θ_mixed is calculated, the adjacent layers are tested for static stability again. If the adjacent layers are statically unstable with the new θ_mixed (e.g., the layer below, altitude-wise, has θ > θ_mixed), the entire process is repeated, including the additional unstable layers, until the entire column is statically stable.

The algorithm for radiative transfer is based on Lacis & Oinas (1991) and was described in Mendonça et al. (2018a). We do not reproduce the entire algorithm in this work, but we make several points of clarification.

We have reverted to using the diffusivity factor, ${ \mathcal D }$, instead of integrating intensities over angle using Gaussian quadrature as was done in Mendonça et al. (2018a). In the double-gray case, this approximation makes very little difference in the calculation of interlayer fluxes. When using multiwavelength radiative transfer, as in Lacis & Oinas (1991), more accurate integration over angle is probably warranted, but the double-gray approximation we use here is likely a much cruder assumption.

Our double-gray scheme is thus a true two-stream approximation, in which the diffusivity factor is used to calculate a characteristic angle for the path of the radiation, and radiation is assumed to be isotropic when the integral over angle is performed. In this approximation, the Planck functions in Equations (A1)–(A6) of Lacis & Oinas (1991) and Equations (3)–(8) of Mendonça et al. (2018a) are replaced by the angle-integrated flux calculated from the Stefan-Boltzmann law. The angle factor μ in those equations, as well as Equations (A9)–(A10) of Lacis & Oinas (1991) and Equations (9)–(10) of Mendonça et al. (2018a), is set to $\mu =1/{ \mathcal D }$, where ${ \mathcal D }$ is the diffusivity factor.

The optical depths are calculated using the form suggested by Frierson et al. (2006) and Heng et al. (2011a). For the short wave, we have a single power law:

$$\begin{eqnarray}&&{\tau }_{\mathrm{sw}}={\tau }_{\mathrm{sw},0}{\sigma }^{{n}_{\mathrm{sw}}},\end{eqnarray} \tag{ 69 }$$

where τ_sw,0 is the optical depth at P = P_ref, σ = P/P_ref, and n_sw is a tuneable factor meant to control the vertical distribution of absorbers. For example, n_sw = 1 would represent a uniformly mixed absorber. A value for n_sw > 1 represents absorbers that are denser in the lower atmosphere.

The long-wave optical depth is given by

$$\begin{eqnarray}&&{\tau }_{\mathrm{lw}}={\tau }_{\mathrm{lw},{\rm{w}}}\sigma +{\tau }_{\mathrm{lw},{\rm{s}}}{\sigma }^{{n}_{\mathrm{lw}}},\end{eqnarray} \tag{ 70 }$$

where τ_lw,w and τ_lw,s represent the surface optical depths of well-mixed absorbers and vertically segregated absorbers, respectively, and n_lw is again a tuneable factor controlling the vertical distribution of the segregated absorbers. With a factor f_l representing the percent of the surface optical depth attributable to the uniformly mixed absorbers, the optical depths above are given by

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}{\tau }_{{\rm{lw}},{\rm{w}}} & = & {\tau }_{{\rm{lw}},0}\,{f}_{l},\\ {\tau }_{{\rm{lw}},{\rm{s}}} & = & {\tau }_{{\rm{lw}},0}(1-{f}_{l}).\end{array}\end{array}\end{eqnarray} \tag{ 71 }$$

The surface optical depth τ_lw,0 can be assumed to be constant (hot Jupiter cases) or given a horizontal distribution. For our Earth-like double-gray case, we give τ_lw,0 a latitudinal dependence given by

$$\begin{eqnarray}&&{\tau }_{\mathrm{lw},0}={\tau }_{\mathrm{lw},\mathrm{eq}}+({\tau }_{\mathrm{lw},\mathrm{pole}}-{\tau }_{\mathrm{lw},\mathrm{eq}}){\sin }^{2}\phi ,\end{eqnarray} \tag{ 72 }$$

where τ_lw,eq and τ_lw,pole are the surface optical depths at the equator and the poles, respectively. This latitudinal dependence approximates the effect of decreased water vapor concentration in the polar regions (Frierson et al. 2006).

The total fluxes passing through each layer are calculated from Equations (1)–(10) of Mendonça et al. (2018a). The heating in the nth layer, used in Equations (18) and (28), is

$$\begin{eqnarray}&&{q}_{\mathrm{heat},n}=\displaystyle \frac{{F}_{n+1}^{\downarrow \mathrm{sw}}-{F}_{n}^{\downarrow \mathrm{sw}}-{F}_{n+1}^{\uparrow \mathrm{lw}}+{F}_{n}^{\uparrow \mathrm{lw}}+{F}_{n+1}^{\downarrow \mathrm{lw}}-{F}_{n}^{\downarrow \mathrm{lw}}}{{\rm{\Delta }}{z}_{n}},\end{eqnarray} \tag{ 73 }$$

where ${F}^{\downarrow \mathrm{sw}}$ is the downward-propagating short-wave (stellar) flux, ${F}^{\uparrow \mathrm{lw}}$ is the upward-propagating long-wave (thermal) flux, ${F}^{\downarrow \mathrm{lw}}$ is the downward-propagating long-wave flux, and Δz_n is the vertical thickness of the layer. Here q_heat has units of W m⁻³, equivalent to kg m⁻¹ s⁻³ or Pa s⁻¹, in line with Equations (23) and (24).

As in Heng et al. (2011a), the surface (when used) is treated as a slab with a constant heat capacity, C_surf. The temperature is modeled using the relation

$$\begin{eqnarray}&&{C}_{\mathrm{surf}}\displaystyle \frac{\partial {T}_{\mathrm{surf}}}{\partial t}={F}_{0}^{\downarrow \mathrm{sw}}-{F}_{0}^{\uparrow \mathrm{lw}}+{F}_{0}^{\downarrow \mathrm{lw}},\end{eqnarray} \tag{ 74 }$$

where ${F}_{0}^{\downarrow \mathrm{sw}}$, ${F}_{0}^{\downarrow \mathrm{lw}}$, and ${F}_{0}^{\uparrow \mathrm{lw}}$ are the short-wave and long-wave fluxes passing downward from and upward into the lowest atmospheric layer.

When no surface is included, as in our hot Jupiter simulations, the flux into and out of the lower boundary is zero, unless flux due to internal heating is included. We do not include flux from the interior in any of the simulations in this work; however, in the case where this is desired, it can be included by specifying an internal flux temperature, T_int.

The user can also set the planetary albedo, A₀. In this double-gray scheme, the albedo represents a top-of-atmosphere albedo and is thus only applied as a modification of the incoming stellar flux, Q.

A benchmark test case for a synchronously rotating Earth using prescribed thermal forcing (Newtonian cooling) was suggested by Heng et al. (2011b) for comparison with the model of Merlis & Schneider (2010), which used gray radiative transfer. We repeat this test simulation here, with input parameters given in Table 2.

The temperature is forced toward an equilibrium profile given by

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}=\max (200\ {\rm{K}},{T}_{\mathrm{HS}}),\end{eqnarray} \tag{ 75 }$$

where

$$\begin{eqnarray}\begin{array}{rcl}{T}_{\mathrm{HS}} & = & \left[315\ {\rm{K}}+{\rm{\Delta }}{T}_{\mathrm{EP}}\cos (\lambda -180^\circ )\cos \phi \right.\\ & & \,-\ \left.{\rm{\Delta }}{T}_{z}\mathrm{ln}\left(\displaystyle \frac{P}{{P}_{\mathrm{ref}}}\right){\cos }^{2}\phi \right]{\left(\displaystyle \frac{P}{{P}_{\mathrm{ref}}}\right)}^{{\kappa }_{\mathrm{ad}}},\end{array}\end{eqnarray} \tag{ 76 }$$

where ΔT_EP is the temperature difference between the substellar and antistellar points (rather than the equator-to-pole difference), ΔT_z is a characteristic scale for vertical temperature differences, and κ_ad = R_d/C_P is the adiabatic coefficient. This formulation is identical to the Held–Suarez test except for the second term of T_HS, which now depends on longitude, λ, to emulate the effect of having the substellar point permanently located at λ = 180° and ϕ = 0°. For the simulation here, ΔT_EP = 60 K and ΔT_z = 10 K, the same as in the Held–Suarez test. The damping timescale for temperature is also identical to the Held–Suarez test and is given by

$$\begin{eqnarray}&&{\tau }_{T}={\left[{k}_{a}+({k}_{s}-{k}_{a})\max \left(0,\displaystyle \frac{\sigma -{\sigma }_{b}}{1-{\sigma }_{b}}\right){\cos }^{4}\phi \right]}^{-1},\end{eqnarray} \tag{ 77 }$$

where k_a = 1/40 day⁻¹, k_s = 1/4 day⁻¹, σ = P/P_surf, and σ_b = 0.7. The surface pressure, P_surf, is calculated at the lower boundary of the lowest layer.

The horizontal winds are damped toward zero with the timescale

$$\begin{eqnarray}&&{\tau }_{v}={\left[{k}_{\mathrm{surf}}\max \left(0,\displaystyle \frac{\sigma -{\sigma }_{b}}{1-{\sigma }_{b}}\right)\right]}^{-1},\end{eqnarray} \tag{ 78 }$$

where k_surf = 1 day⁻¹.

Figure 4 shows results from the simulation with g_level = 5 (a horizontal resolution of ∼2°), plotted on two isobaric surfaces. Results are very similar for g_level = 4. At P = 0.9 bars, near the surface, we see convergence toward the substellar point (located at longitude λ = 180°), associated with the rising motion due to the intense heating at this location. At P = 0.25 bars, the flow diverges from the substellar point and then converges on the night side of the planet.

Zoom In Zoom Out Reset image size

Figure 5 shows several properties as a function of latitude. On the night side, the temperature is highest at P ∼ 0.8 bars, well above the surface. On the day side, the temperature is highest at the equatorial surface. The transport of heat upward and away from the substellar point is apparent from the temperature distribution centered on λ = 180° and the stream function, which shows a large equator-to-pole Hadley cell in each hemisphere. The potential temperature distribution shows that the atmosphere is marginally unstable near the substellar point. Convective adjustment was not included here.

Zoom In Zoom Out Reset image size

For this test, we have not included heating and cooling of the surface. The pressure in the lowest atmospheric layer varies from ∼0.93 bars near the antistellar point to ∼0.91 bars near the substellar point. Extrapolating the atmospheric temperature to the reference pressure (roughly the location of the implied surface) for plotting purposes produces poor results because of the steep gradients at the bottom of the model that are not well resolved by our uniform vertical mesh; hence, we do not attempt to extend our figures to this region.

In general, the results compare well with Heng et al. (2011b) and Merlis & Schneider (2010). The temperatures and wind speeds are similar to those of Heng et al. (2011b), with a temperature peak at ∼320 K at the equator, near the surface, and max wind speeds of ∼13 m s⁻¹ and ∼25 m s⁻¹ at P = 0.9 and P = 0.25 bars, respectively. The fact that the nightside temperature in Figure 4 is warmer than that which appears in Figure 3 of Heng et al. (2011b) is due to the difference in pressure level, owing to the poor vertical resolution near the surface in our simulation. The temperatures are ∼20 K lower across the globe in the simulation of Merlis & Schneider (2010), which utilized a radiative transfer scheme; the fact that our temperatures agree with those in Heng et al. (2011b) suggests that this is due to the tuning of the temperature forcing in Equations (75) and (76) rather than an error in the code. The temperature distribution is otherwise similar to Figure 9 in Merlis & Schneider (2010), with a maximum near 0.8 bars on the night side. The Hadley cells are similar in size to theirs, though it is about a factor of 2 weaker in our simulation.

Here we attempt to reproduce the deep hot Jupiter benchmark test from Heng et al. (2011b). Input parameters are given in Table 2. Like the Held–Suarez test, the benchmark is run with an idealized forcing to the temperature (the winds, in this case, are unforced). As noted in Mayne et al. (2014), this benchmark has several challenges. First, in a nonhydrostatic model, where the vertical coordinate is altitude rather than pressure, the night side of the planet tends to extend to several orders of magnitude lower pressure than the day side. The exact temperature–pressure profiles suggested by Heng et al. (2011b) tend to lead to runaway cooling on the night side at the start of the simulation, causing the model to crash. Mayne et al. (2014) successfully mitigated this issue in the UK Met Office model by increasing the temperatures at low pressures. We have attempted the same here, with less success. The second issue is the discontinuity in the temperature–pressure profiles at 10 bars, which can lead to numerical instabilities. To avoid the issue, we refit the temperature–pressure profiles used in Heng et al. (2011b) with a new set of polynomials, excluding the pressures near 10 bars. This produces a new profile that varies smoothly in this region. The new polynomial fits are shown in Figure 6 and presented in Appendix A.

Figure 7 shows the zonal-mean temperature and zonal wind for the resulting simulation. Unfortunately, our model domain is limited to pressures ≳10 mbar on the day side of the planet. Raising the model top any further causes the nightside instability (noted by Mayne et al. 2014 and described above) to occur. At present, it is unclear why the UK Met Office model is able to successfully extend the model domain to lower pressures while THOR is not. However, based on the issues discussed here and in Mayne et al. (2014), it appears that the deep hot Jupiter benchmark is a challenge to reproduce in altitude-grid models.

We have made several attempts to extend the model domain, including initializing the atmosphere from the average T_eq (blue curve in Figure 6) rather than isothermal conditions, including a sponge layer (see Section 2.4) at the top of the model, tuning the strength of divergence damping and hyperdiffusion, and initializing with a zonal wind profile to prompt dayside-nightside mixing. None of the above efforts were successful at preventing the runaway cooling on the night side that leads to a model instability.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Even with the model domain limited to P ≳ 10 mbar on the day side, we do reproduce many of the features of the experiment done in previous works. We see a zonal-mean temperature that is similar in magnitude and structure to that of Heng et al. (2011b) and Mayne et al. (2014). We see equatorial superrotation and return flow at the midlatitudes. The jet speed we find here is weaker (∼3600 m s⁻¹) than in the finite-volume simulation of Heng et al. (2011b; ∼5000 m s⁻¹) and the simulation in Mayne et al. (2014; ∼6000 m s⁻¹). This may be caused by the lower location of the model top (the velocities tend to increase with altitude), though it is also possibly explained by the low horizontal resolution we used for the simulation here (we revisit the problem of resolution in Section 5).

Here we demonstrate the THOR GCM's representation of acoustic waves. The purpose of this test is to characterize how well the model is able to represent the propagation of acoustic waves and provide a tool to isolate coding errors that may be hard to diagnose in more complicated scenarios. The setup is almost identical to Section 4.1 of Tomita & Satoh (2004), but we describe it here for completeness. The atmosphere is initialized with a background isothermal state, for an Earth-radius planet with no rotation and no forcing (radiation or Newtonian cooling). At the beginning of the first time step, a pressure perturbation is applied over a spatial distribution centered on longitude λ₀ = 0° and latitude ϕ₀ = 0°. The distribution in latitude ϕ and longitude λ is given by

$$\begin{eqnarray}\xi (\lambda ,\phi )=\left\{\begin{array}{ll}\tfrac{1}{2}(1+\cos (\pi x/L)) & x\lt L,\\ 0 & x\gt L,\end{array}\right.\end{eqnarray} \tag{ 79 }$$

where L is a representative horizontal length and x is the horizontal distance along a great circle from the point (λ₀, ϕ₀) and is given by

$$\begin{eqnarray}&&x={r}_{0}{\cos }^{-1}\left[\sin {\phi }_{0}\sin \phi +\cos {\phi }_{0}\cos \phi \cos (\lambda -{\lambda }_{0})\right].\end{eqnarray} \tag{ 80 }$$

The vertical distribution is

$$\begin{eqnarray}&&\zeta (z)=\sin \left(\displaystyle \frac{{n}_{v}\pi z}{{z}_{\mathrm{top}}}\right),\end{eqnarray} \tag{ 81 }$$

where z is the altitude, z_top is the height of the model top, and n_v is the vertical wave mode. The total initial perturbation field is

$$\begin{eqnarray}&&p^{\prime} (t=0)=\delta p\,\xi (\lambda ,\phi )\zeta (z),\end{eqnarray} \tag{ 82 }$$

where δp is the amplitude of the perturbation. Here, as in Tomita & Satoh (2004), we set δp = 100 Pa, r₀ = 6371 km, L = r₀/3, n_v = 1, and z_top = 10 km. Further input parameters are given in Table 3.

We compare results from two simulations. In the first, both hyperdiffusion and divergence damping are omitted (by "divergence damping" we are referring to the terms that depend on the divergence of momentum in Equations (44)–(46) of Mendonça et al. 2016, Equation (58) in this work). In the second, only divergence damping is enabled, with a coefficient D_div = 0.02 (hyperdiffusion is still disabled). In Figure 8, we show the resulting pressure field in the lowest altitude level at different times for the case with divergence damping. These compare well with Figure 3 of Tomita & Satoh (2004). We note, however, that the amplitude of the pressure perturbation in their figure appears much larger than in ours, peaking at p' ≈ 1000 Pa, which seems inconsistent with an initial perturbation of δp = 100 Pa as stated in the text. We also see that the amplitude decreases as the wave spreads away from the point of origin (evident in the ranges of the color scales in Figure 8), which is not seen in the Tomita & Satoh (2004) figure, unless their color scale is mislabeled.

Zoom In Zoom Out Reset image size

In Figure 9, we plot the globally integrated total energy, internal plus potential energy, and kinetic energy as a function of time. Compare with Figure 4 of Tomita & Satoh (2004)—though note that their figure has different units. There appears to be a slight mismatch between the exact hydrostatic initial conditions and the THOR algorithm's representation of hydrostatic balance. This leads to a jump in the total energy (on the order of 10¹⁴ J) on the first time step as all columns of the atmosphere adjust very slightly. At present the origin of the discrepancy is unclear, but in any case, the error is quite small compared to the total energy of the atmosphere and should not noticeably affect the results of simulations that include forcing, which produces a much larger change in the overall energy budget compared to the initial conditions. The pressure perturbation is applied at the end of the first time step, reaches the opposite side of the planet in ≈17 hr, and returns to the original location in ≈33 hr, indicating a sound speed of ∼337 m s⁻¹, as found in Tomita & Satoh (2004). The theoretical sound speed is ${c}_{s}=\sqrt{\gamma {R}_{d}T}\approx 347$ m s⁻¹, where γ = C_P/C_V.

Zoom In Zoom Out Reset image size

From Figure 9, we can see that the variation of kinetic energy associated with the sound wave is well compensated by the variation in internal plus potential energy. In the case without divergence damping, the total energy begins to increase erroneously after ∼35 hr. The simulation crashes not long after this time. This is the result of grid-scale noise that would be well eliminated by the divergence damping terms. A comparison between the damped and undamped cases of the global pressure field at 48.5 hr is shown in Figure 10. In the undamped case, we see spurious waves originating from the pentagonal grid points, but these waves are fully eliminated in the damped case.

Zoom In Zoom Out Reset image size

The case with divergence damping conserves the energy much better, though there is a slight energy loss over the course of the simulation. Comparing with Tomita & Satoh (2004), THOR does not conserve energy as well as their model NICAM. Our choice of entropy as the fundamental quantity in the thermodynamic equation rather than total energy means that the total energy is not conserved as precisely (Satoh 2002, 2003). Mendonça et al. (2016) tested an energy correction scheme in THOR based on Williamson et al. (2009) that, while improving total energy conservation, had little impact on the overall behavior of the model. Thus, until we can develop a more robust general approach that is physically well described at coarse resolutions, we choose to live with the gradual loss of energy to the dissipation scheme.

A further benchmark compares the representation of gravity waves in THOR with NICAM. The motivation behind this type of simulation is to characterize the model's representation of gravity waves. This test is originally presented in Tomita & Satoh (2004). We set up the atmosphere in the same manner as in that paper and in the acoustic wave experiment of the previous section, using the same planet radius, model top, and shape of the perturbation. See Table 3 for the input parameters. The key differences between these simulations and that in Section 4.3 are that the atmosphere is given a thermal profile corresponding to a constant Brunt–Väisälä frequency, instead of an isothermal profile, and the perturbation is applied to the potential temperature, given by

$$\begin{eqnarray}&&\theta ^{\prime} (t=0)=\delta \theta \,\xi (\lambda ,\phi )\zeta (z),\end{eqnarray} \tag{ 83 }$$

where θ' is the potential temperature perturbation function and δθ is its amplitude. The horizontal and vertical distributions, ξ and ζ, are again given by Equations (79) and (81). In all cases, δθ = 10 K.

The atmosphere for this case is initialized with a constant Brunt–Väisälä frequency, N. This can be defined in terms of the potential temperature as

$$\begin{eqnarray}&&{N}^{2}=\displaystyle \frac{g}{\theta }\displaystyle \frac{\partial \theta }{\partial z}.\end{eqnarray} \tag{ 84 }$$

In order to reduce numerical instabilities and to minimize motion introduced by the initial conditions, we initialize the atmosphere from hydrostatic equilibrium. For an atmosphere with a constant N, we must numerically solve the hydrostatic equation for the pressure of each layer as a function of the layer below. In our case, we begin at the reference pressure P_ref, which is the pressure at the bottom of the lowest layer at the beginning of the simulation. From there, we progress upward, numerically determining the pressure of each successive layer. More specifically, we use Newton–Raphson iteration to find the pressure in the ith layer that satisfies the equation

$$\begin{eqnarray}&&\displaystyle \frac{{P}_{i}-{P}_{i-1}}{{\rm{\Delta }}z}=\displaystyle \frac{1}{2}({\rho }_{i}+{\rho }_{i-1})g.\end{eqnarray} \tag{ 85 }$$

The relationship between the pressure and density is given by the ideal gas law. An additional constraint must then be made on the temperature, which can be derived from the definition of N:

$$\begin{eqnarray}&&{T}_{i}=\displaystyle \frac{1+\beta (N)}{1-\beta (N)}{T}_{i-1},\end{eqnarray} \tag{ 86 }$$

where

$$\begin{eqnarray}&&\beta (N)=\displaystyle \frac{1}{2}\displaystyle \frac{{N}^{2}}{g}{\rm{\Delta }}z+\displaystyle \frac{{R}_{d}}{{C}_{P}}\displaystyle \frac{{P}_{i}-{P}_{i-1}}{{P}_{i}+{P}_{i-1}}.\end{eqnarray} \tag{ 87 }$$

The quantities P_i−1, ρ_i−1, and T_i−1 are the pressure, density, and temperature of the layer below the ith, except in the case of the lowest layer, where P_i−1 = P_ref, ${\rho }_{i-1}={P}_{\mathrm{ref}}/({R}_{d}{T}_{\mathrm{init}})$, and T_i−1 = T_init. T_init is a user-set initial reference temperature. The value of Δz is the distance between the centers of the layers except in the case of the lowest layer, where it is simply the height of the lowest layer. After making the initial guess that ρ_i = ρ_i−1, we do one Newton–Raphson step to find a new value of P_i, update T_i and ρ_i (via Equation (86)), and then repeat the process until the change in P_i is below some threshold (we use 10⁻⁸).

In Figure 11, we show the results of three simulations. The first has (N, n_v) = (0.01 s⁻¹, 1), the second (N, n_v) = (0.02 s⁻¹, 1), and the third (N, n_v) = (0.01 s⁻¹, 2). The number of vertical levels is 20 in the n_v = 1 cases and 40 in the n_v = 2 case. The color contours show the temperature perturbation, ΔT, at the equator after 48 hr of integration. Tomita & Satoh (2004) give a theoretical estimate for the gravity wave speed, c_g = Nz_top/πn_v, which is 31.8, 63.7, and 15.9 m s⁻¹ in the three respective cases. The leading peaks in the three cases are located at λ ≈ (55°, 95°, 30°), indicating speeds c_g ≈ (35, 61, 19) m s⁻¹. Tomita & Satoh (2004) noted, in particular, the larger relative error in their n_v = 2 case (c_g ≈ 18 m s⁻¹) and theorized that this was caused by insufficient vertical resolution—they used 20 vertical levels in this case. However, we have run the same test with 20 levels and 40 levels (the bottom panel of Figure 11 shows the 40 level case), and the locations of the wave peaks are virtually identical in both cases. Most likely, the error in the speed is dominated simply by the difficulty in estimating the locations of the wave fronts.

Zoom In Zoom Out Reset image size

Figure 12 shows the evolution of the total, internal plus potential, and kinetic energy. As in the acoustic wave case, the kinetic energy mirrors the change in internal plus potential energy. The oscillation in the total energy is orders of magnitude smaller than the others, indicating that the total energy is conserved well. This compares well with Figure 6 in Tomita & Satoh (2004); however, as noted in the acoustic wave case, we plot the absolute energy (in J).

Zoom In Zoom Out Reset image size

Here we present a dry Earth-like case to compare to the Held–Suarez test and validate the double-gray radiative transfer scheme. Input parameters are given in Table 4. We utilize optical depth profiles for the short-wave and long-wave radiation in the same style as Frierson et al. (2006) and Heng et al. (2011a), with n_sw = 2, n_lw = 4, τ_sw = 0.2, τ_lw,eq = 6, and τ_lw,pole = 1.5. In this case, the insolation, or the distribution of incident solar/stellar radiation at the top of the atmosphere, is given by

$$\begin{eqnarray}&&Q={S}_{0}\cos \phi \cos (\lambda -\alpha ),\end{eqnarray} \tag{ 88 }$$

where S₀ = 1367 W m⁻², ϕ is the latitude, λ is the longitude, and the angle α at time t is

$$\begin{eqnarray}&&\alpha =({\rm{\Omega }}-{n}_{\mathrm{orb}})t.\end{eqnarray} \tag{ 89 }$$

Here Ω = 7.292 × 10⁻⁵ is the rotation rate and n_orb = 1.98 × 10⁻⁷ is the orbital mean motion. This insolation pattern assumes a zero eccentricity orbit and zero obliquity but resolves the diurnal cycle. With zero eccentricity, the mean motion n is simply the angular velocity of the planet about the Sun. In a future work, we will generalize the insolation to arbitrary orbital and rotation states.

The first case is nonhydrostatic, deep (NHD) without dry convective adjustment (Figure 13). The second case is NHD with dry convective adjustment (Figure 14). Comparing the two cases, we can see that without dry convection there is clear spurious heating at ∼20°–30° latitude. With convective adjustment enabled, we produce a temperature profile much more similar to the cases in Heng et al. (2011a). As in that work, we see that the potential temperature profile, in the case without convective adjustment, is unstable near the surface at latitudes ≲50°. The convective adjustment scheme rectifies the situation and cools the near-surface region near ∼20°–30° latitude.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In both cases, we see jet streams form at ϕ ∼ 30° and P ∼ 0.4 bars (top right panels of Figures 13 and 14). There is a hint of a second, weaker stream in each hemisphere at ∼60°. The Eulerian mass stream function (bottom right panels of Figures 13 and 14) is similar in the two cases. We see strong thermally direct cells (Hadley cells) at the equator, and weaker indirect cells adjacent to these. These overturning cells are narrower in latitude than the real Hadley and Ferrel cells seen on Earth. We speculate that this may be a consequence of the lack of a hydrologic cycle in the current version of our model, as narrow Hadley cells were also observed in the dry simulations of Heng et al. (2011b), compared with, for example, Merlis & Schneider (2010).

We calculate the Eulerian mass stream function using the following relation:

$$\begin{eqnarray}&&{\rm{\Psi }}=\displaystyle \frac{2\pi {r}_{0}\cos \phi }{g}{\int }_{0}^{P}\overline{[v]}{dP},\end{eqnarray} \tag{ 90 }$$

where Ψ is the stream function, r₀ is planet radius (at the model bottom), g is the gravity, and $\overline{[v]}$ is the time and zonally averaged meridional velocity. Implicit in this definition is hydrostatic balance; in the nonhydrostatic model, this condition is not exactly satisfied. However, the simulation is at all times near enough to hydrostatic equilibrium that there is little difference between in Ψ as calculated above compared to a calculation based on Equation (4.1) of Peixóto & Oort (1984), which does not assume hydrostatic balance.

We compare the model run using the QHD and hydrostatic, shallow (HSS) approximations to the NHD case, all with convective adjustment enabled. Figure 15 shows the difference in temperature and zonal wind for the QHD and HSS simulations from the NHD case. Over most of the model domain, the differences are relatively small, on average. The largest differences occur in wind speeds in the upper atmosphere between the HSS case and NHD cases. This suggests that nonhydrostatic effects and the geometric corrections for a deep atmosphere are unimportant in this regime, at least in terms of the global average state of the atmosphere, but that the geometric corrections are more important than nonhydrostatic effects.

Zoom In Zoom Out Reset image size

We have run a single NHD case at g_level = 6, corresponding to a horizontal resolution of ∼1°. The output (not shown) is qualitatively very similar to the g_level = 5 case, though, as we discuss below, the simulation conserves mass to slightly less precision than the lower-resolution cases.

Lastly, we compare several diagnostics for the Earth-like simulations. Figure 16 shows the evolution of the global atmospheric mass, energy, and axial angular momentum. We hope to conserve the total mass; the energy and angular momentum evolve owing to input from the stellar radiation and are thus more diagnostic of convergence. The simulations at resolution ∼2° all conserve mass to a similar precision, a few parts in 10¹² over 1200 days. The simulation at resolution of ∼1° does slightly worse in this respect. In all cases, the model reaches a steady state in ∼200–300 days, after which the total energy and axial angular momentum stay roughly constant.

Zoom In Zoom Out Reset image size

Here we show a comparison of QHD and NHD simulations of the hot Jupiter HD 189733 b. The input parameters are summarized in Table 4. Simulations including the shallow approximation (HSS) are very similar to the QHD simulations; thus, we stick here to a comparison of QHD and NHD. The largest differences in flow for this planet occur between the QHD and NHD simulations, indicating that the primary source of the differences described below is the term Dv_r/Dt, the Lagrangian derivative of the vertical velocity.

Figure 17 shows the zonally and temporally averaged temperature and zonal wind speed during the last 1000 days of the NHD and QHD simulations. The overall temperature structure is very similar between the two cases. The zonal wind speed plots indicate the presence of superrotation, as expected, and the overall structure is similar. However, the maximum velocities of the jets differ by ∼5%, broadly consistent with the comparison between the "Prim" and "Full" simulations of Mayne et al. (2017). While that study compared simulations using the primitive equations with the full nonhydrostatic equations, the sole difference between our simulations here is the neglect of the material derivative of the vertical velocity, Dv_r/Dt, and the hyperdiffusive term, ${{ \mathcal F }}_{{v}_{r}}$, in the QHD case. This results in a difference of the maximum jet speed, likely because the 3D representation of waves has been modified. Though the difference is small at this resolution, it becomes more pronounced in the higher-resolution simulations, as we describe shortly.

Zoom In Zoom Out Reset image size

Figure 18 shows the zonally and temporally averaged potential temperature and Eulerian mass stream function of the NHD and QHD cases. Though the dry convection scheme is enabled in these simulations, it likely has very little effect, as everywhere the atmosphere is quite stable, as indicated by the potential temperature structure. Though this quantity is averaged over longitude, local profiles of the potential temperature at different longitudes near the equator also show stability. In the plots for the mass stream function, we see thermally indirect overturning cells at the base of the equatorial jet, as also seen in Charnay et al. (2015) and Mendonça et al. (2018a). Thermally direct cells appear adjacent to the indirect cells at higher latitudes. Overturning in the upper atmosphere is weak in comparison and indiscernible on this scale.

Zoom In Zoom Out Reset image size

Figure 19 shows snapshots (at 10,000 days) of the temperature and horizontal wind speeds along isobars for the NHD and QHD simulations. At 0.1 bars and above (in altitude), in the region of the superrotational jet, we see the characteristic chevron shape and eastward hot spot offset associated with hot Jupiters. The wind vectors indicate eastward motion and that the gas is pushed away from the substellar point (at 0° longitude). At 1 mbar, we see return flow on the night side at high latitudes, which results in convergence on the equator at the eastern terminator and the recognizable chevron shape. The NHD and QHD simulations show only minor differences in temperature, but the velocities are slightly higher in the NHD case. Standing Rossby waves are easily discernible at these pressure levels. At the 10 bar level, flow at the equator is retrograde, though comparatively sluggish (∼10 m s⁻¹).

Zoom In Zoom Out Reset image size

Figure 20 shows the temperature and zonal wind speed, as in Figure 17, but for simulations with a resolution of ∼2°. Temperatures shows a similar pattern to the ∼4° resolution simulations, but the superrotational jet has increased in speed and has broadened in latitude, pushing the return flow toward higher latitudes. The increase in velocity is likely due to better resolution of waves that carry angular momentum toward the equator. The peak velocity of the NHD simulation is ∼20% greater than in the QHD simulation, similar to the effect seen at ∼4° resolution. The potential temperature and stream function for the $\sim {2}^{\circ }$ simulations are shown in Figure 21. The temperature and wind speeds along isobars are shown in Figure 22.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 23 shows the zonal and vertical velocities in a latitudinal band along the equator, averaged over the last 1000 days, for the NHD and QHD simulations at ∼4° resolution. Here we have used altitude as our vertical coordinate to avoid extrapolation at the top of the model. The values are averaged over a 20° latitudinal band weighted by $\cos \phi$, where ϕ is the latitude. Zonal winds are highest on the night side (longitudes 90°–270°) toward the western terminator. As the flow approaches the day side, it slows because of the increase in pressure. Vertical velocities are slow everywhere in comparison to the horizontal winds, though the figures show a lot of structure. Upwelling and downwelling occur, for example, side by side along the western terminator. Upwelling dominates on the day side of the planet, with downwelling on the night side, most prominently east of the eastern terminator. Figure 24 shows equatorial temperatures in the same style, comparing the NHD and QHD cases.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The same quantities are shown for the ∼2° simulations in Figures 25 and 26. While the temperature structure is largely unchanged from the ∼4° cases, the velocities have changed more significantly. In particular, there is now a strong difference between the peak zonal winds in the NHD and QHD cases. The zonal wind speeds have increased, the vertical wind speeds have slightly decreased, and the flow now extends deeper into the atmosphere, as also seen in Figure 20.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size