The Limits of the Primitive Equations of Dynamics for Warm, Slowly Rotating Small Neptunes and Super Earths

TL;DR

This study reveals that using simplified primitive equations instead of full Navier-Stokes equations significantly affects atmospheric flow simulations of warm, slowly rotating small Neptunes and super Earths, especially under certain conditions.

Contribution

It demonstrates the limitations of primitive equations for modeling atmospheres of small, slowly rotating exoplanets, highlighting the need for more comprehensive models.

Findings

Differences in atmospheric flow patterns between primitive and full Navier-Stokes simulations.

Primitive equations become less accurate as rotation slows and temperature contrasts increase.

Significant impact on heat redistribution and atmospheric temperature distribution.

Abstract

We present significant differences in the simulated atmospheric flow for warm, tidally-locked small Neptunes and super Earths (based on a nominal GJ 1214b) when solving the simplified, and commonly used, primitive dynamical equations or the full Navier-Stokes equations. The dominant prograde, superrotating zonal jet is markedly different between the simulations which are performed using practically identical numerical setups, within the same model. The differences arise due to the breakdown of the so-called `shallow-fluid' and traditional approximations, which worsens when rotation rates are slowed, and day-night temperature contrasts are increased. The changes in the zonal advection between simulations solving the full and simplified equations, give rise to significant differences in the atmospheric redistribution of heat, altering the position of the hottest part of the atmosphere and temperature contrast between the day and night sides. The implications for the atmospheric chemistry and, therefore, observations need to be studied with a model including a more detailed treatment of the radiative transfer and chemistry. Small Neptunes and super Earths are extremely abundant and important, potentially bridging the structural properties (mass, radius, composition) of terrestrial and gas giant planets. Our results indicate care is required when interpreting the output of models solving the primitive equations of motion for such planets.

Figures

Tables

Table 1. Table 1: Values of the standard parameters for simulations of GJ 1214b.

Quantity	Value
Horizontal resolution, $N_{\lambda }$, $N_{\varphi }$	${144}_{\lambda }$, ${90}_{\varphi }$
Vertical resolution, $N_z$	66
Dynamical Timestep, $\Delta t$ (s)	120
Initial inner boundary pressure, $p_{\max }$ (Pa)	$200\times {10}^5$
Rotation rate, $\Omega$ (s-1)	$4.60\times {10}^{-5}$
Radius, $R_p$ (m)	$1.45\times {10}^7$
Atmospheric height $z_{\mathrm{top}}$ (m)	$3.7\times {10}^6$
Surface gravity, $g_p$ (ms-2)	12.20
Specific heat capacity (constant pressure), $c_p$ (Jkg-1K-1)	12 300
Ideal gas constant, $R$ (Jkg-1K-1)	3573.5
Temperature contrast, $\Delta T_{\mathrm{eq}}$ (K)	600
Diffusion setting, $K_{\lambda }$ (see Mayne et al., 2017, for details)	0.158
Vertical, ‘sponge’, damping coefficient $R_w$	0.15

Table 2. Table 2: Short names of the simulations presented in this work of GJ 1214b, with variations in parameters from those shown in Table 1 , atmospheric height and the dynamical equations set used (see Mayne et al., 2014b , for explanation and detail) .

Short Name	Adjusted parameters	$z_{\mathrm{top}}$ (m)	Equation set
Std Full	-	$3.70\times {10}^6$	Full
Std Prim	-	$3.00\times {10}^6$	Primitive
Std Hires Full	$N_{\lambda }=288$, $N_{\varphi }=180$, $N_z=132$, $\Delta t=60$ s	$3.70\times {10}^6$	Full
Std Hires Prim	$N_{\lambda }=288$, $N_{\varphi }=180$, $N_z=132$, $\Delta t=60$ s	$3.00\times {10}^6$	Primitive
$R_p-$ Full	R${}_p{}^{}=6.0\times {10}^6$ m, $\Delta t=80$ s	$5.50\times {10}^6$	Full
$R_p-$ Prim	R${}_p{}^{}=6.0\times {10}^6$ m, $\Delta t=80$ s	$3.00\times {10}^6$	Primitive
$R_p+$ Full	R${}_p{}^{}=1.0\times {10}^8$ m	$3.00\times {10}^6$	Full
$R_p+$ Prim	R${}_p{}^{}=1.0\times {10}^8$ m	$2.80\times {10}^6$	Primitive
CO2 Full	$c_p=\mathrm{10\hspace{0.17em}123}$ Jkg-1K-1, $R=188.9$,Jkg-1K-1, $\Delta t=60$ s	$1.50\times {10}^5$	Full
CO2 Prim	$c_p=\mathrm{10\hspace{0.17em}123}$ Jkg-1K-1, $R=188.9$ Jkg-1K-1, $\Delta t=60$ s	$1.50\times {10}^5$	Primitive
$\Omega +$ Full	$\Omega =9.2\times {10}^{-5}$ s-1	$3.70\times {10}^6$	Full
$\Omega +$ Prim	$\Omega =9.2\times {10}^{-5}$ s-1	$3.00\times {10}^6$	Primitive
dT$+$ Full	$\Delta T_{\mathrm{eq}}$ 800 K	$3.30\times {10}^6$	Full
dT$+$ Deep	$\Delta T_{\mathrm{eq}}$ 800 K, $g\left(r\right)=g_p$	$2.65\times {10}^6$	Deep
dT$+$ Prim	$\Delta T_{\mathrm{eq}}$ 800 K	$2.65\times {10}^6$	Primitive

Equations

$$\bar{F}\left(\delta \right)=\frac{F\left(\delta \right)-⟨F\left(\delta \right)⟩}{⟨F\left(\delta \right)⟩},$$

$$\frac{Du}{Dt}=\frac{uv\tan \varphi }{r}-\frac{uw}{r}+fv-f^{{}^{\prime }}w-\frac{c_p\theta }{r\cos \varphi }\frac{\partial \Pi }{\partial \lambda }+\text{D}(u),$$

$$f=2\Omega \sin \varphi ,$$

$$f^{{}^{\prime }}=2\Omega \cos \varphi .$$

$$\frac{Du}{Dt}=\frac{uv\tan \varphi }{R_{\mathrm{p}}}+fv-\frac{c_p\theta }{R_{\mathrm{p}}\cos \varphi }\frac{\partial \Pi }{\partial \lambda }+\text{D}(u),$$

$$(\vec{v}\cdot \vec{\nabla })\vec{v}-2\vec{\Omega }\wedge \vec{v}=-\frac{1}{\rho }\vec{\nabla }p+\vec{g},$$

$$\frac{\mathrm{D}\rho }{\mathrm{D}t}=0,$$

$$\frac{\mathrm{D}T}{\mathrm{D}t}=\frac{T_{\mathrm{eq}}-T}{{\tau }_{\mathrm{rad}}},$$

$$p=\rho RT,$$

$$\vec{v}\cdot\vec{\nabla}\bigg{(}\dfrac{v^{2}}{2}\bigg{)}=-\dfrac{1}{\rho}\dfrac{\mathrm{D}p}{\mathrm{D}t}-wg.$$

$$\frac{U}{L_U}\sim \frac{W}{H_W},$$

$$-\frac{1}{\rho }\frac{Dp}{Dt}-wg=-\frac{1}{\rho }(\vec{v}\cdot {\vec{\nabla }}_{\perp })p\sim \frac{1}{\rho }\frac{Dp}{Dt},$$

$$\frac{Dp}{Dt}=\rho R\frac{DT}{Dt}=\rho R\frac{(T_{\mathrm{eq}}-T)}{{\tau }_{\mathrm{rad}}}.$$

$$\vec{v}\cdot\vec{\nabla}\bigg{(}\dfrac{v^{2}}{2}\bigg{)}\sim R\dfrac{(T_{\mathrm{eq}}-T)}{\tau_{\mathrm{rad}}}\\ ,$$

$$\frac{U^3}{2L_U}\sim R\frac{\Delta T}{{\tau }_{\mathrm{rad}}},$$

$$U\sim \sqrt[3]{2\pi R_{\mathrm{p}}R\frac{\Delta T}{{\tau }_{\mathrm{rad}}}}\sim 1.4\times 10^3{\mathrm{ms}}^{-1}.$$

$$W\sim \frac{H}{L}\sqrt[3]{2\pi R_{\mathrm{p}}R\frac{\Delta T}{{\tau }_{\mathrm{rad}}}},$$

$$\frac{Dv}{Dt}\sim 2\Omega U\sin \varphi ,$$

$$\frac{UV}{L_V}\sim 2\Omega U\sin \varphi ,$$

$$\begin{gathered} \frac{UV}{L_V}\sim 2\Omega U\sin \varphi \\ V\sim 2L_V\Omega \sin \varphi . \end{gathered}$$

$$V\sim \frac{\pi }{2}{R}_{\mathrm{p}}\Omega \sin \varphi \sim 400\sin \varphi {\mathrm{ms}}^{-1}.$$

$$V\tan \varphi \sim \frac{\pi }{2}{R}_{\mathrm{p}}\Omega \frac{\sin {\varphi }^2}{\cos \varphi },$$

Full text

editorials, notices — miscellaneous — catalogs — surveys

††journal: ApJ

1 Introduction

The detection of exoplanets has revealed several classes of object without direct analogues in our Solar system. One important example, termed hot Jupiters, Jovian type planets in short period orbits provide the strongest atmospheric signatures via observations such as emission spectra (e.g., Todorov et al., 2014), transmission spectra (e.g., Sing et al., 2011) and emission as a function of orbital phase, or phase curves (e.g., Knutson et al., 2012). However, these same techniques are beginning to be applied to a second class of planet, those between the size of Neptune and Earth, termed Super Earths and small Neptunes. Importantly, this class of planet has been shown to be the most abundant, when adopting the distinctions of Fressin et al. (2013), namely 1.25$R_{\earth}<R_{\rm p}<$4$R_{\earth}$, where $R_{\earth}$ & $R_{\rm p}$ are the Earth and planetary radii, respectively. In addition, the Next Generation Transit Survey (NGTS, Wheatley et al., 2013), CHaracterising ExoPlanet Satellite (CHEOPS, Broeg et al., 2013) and Transiting Exoplanet Survey Satellite (TESS, Ricker, 2014) will significantly increase the number of detections, and availability of characterisation observations, in this size range in the near future. For example, the recent detection of $\pi$ Mensae c (HD 39091c) using data from TESS (Huang et al., 2018).

Aside from the abundance of potential targets, Super Earths and small Neptunes also inhabit a potentially critical region of planetary parameter space. These planets likely bracket the point at which runaway accretion of a primary gaseous atmosphere occurs in the core accretion model (Pollack et al., 1996). Therefore, they bridge the structures of giant planets with thick hydrogen/helium dominated atmospheres, to terrestrial planets with much thinner ‘secondary’ atmospheres (Lopez & Fortney, 2014), as well as being in the size range where irradiative evaporation becomes significant (Owen & Jackson, 2012). Finally, as these planets are generally cooler than their hot Jupiter counterparts, the characteristic timescale to reach chemical equilibrium for many species increases, leading to stronger potential effects from chemical kinetics driven by transport and photochemistry (see Madhusudhan et al., 2016, for review).

Observations of several targets in this regime have demonstrated the complex nature of their atmospheres. The Super Earth GJ 1214b (equilibrium temperature $\sim$500 K), has been intensively observed returning a flat, featureless spectrum via both ground-based (Bean et al., 2010, 2011; Crossfield et al., 2011; de Mooij et al., 2012; Wilson et al., 2014; Cáceres et al., 2014) and space-based instruments (Désert et al., 2011; Berta et al., 2012; Fraine et al., 2013; Kreidberg et al., 2014). Kreidberg et al. (2014) effectively ruled out a cloud–free atmosphere, but did so without constraint on the bulk composition (i.e. metallicity) of the gas–phase. Additionally, for the ‘warm Neptune’ GJ 436b (slightly larger than the small Neptune limit at $\sim$4.3$R_{\earth}$, Deming et al., 2007) observations indicate a condensate–rich atmosphere (Knutson et al., 2014).

The main focus of the theoretical modeling of exoplanets has been either the well observed hot Jupiter objects, or potentially habitable terrestrial planets. For terrestrial planets 3D General Circulation Models (GCMs) adapted from those used to study Earth, have been applied to explore the potential climates of, for example, Proxima Centauri b (Turbet et al., 2016; Boutle et al., 2017), the response of an Earth–like climate to weakened and intensified stellar irradiation (Leconte et al., 2013; Charnay et al., 2013, respectively), alongside studies including surface effects (Lewis et al., 2018) and a dynamical ocean (e.g., Wolf et al., 2017; Del Genio et al., 2017). For hot Jupiters, gas–phase chemical equilibrium simulations have been performed, in 3D using GCMs, across a range of targets (e.g., Kataria et al., 2016). Treatments of the chemical kinetics have also been included using simplified “relaxation methods” in 3D (Cooper & Showman, 2006; Drummond et al., 2018c, b) or more complete chemical networks in pseudo-2D (Agúndez et al., 2012, 2014). Finally, the treatment of clouds has been added both ‘diagnostically’ (where feedbacks of the cloud presence on the atmosphere are neglected, Parmentier et al., 2016; Helling et al., 2016), and ‘prognostically’ (where the clouds evolve and interact with the atmosphere, Lee et al., 2016; Lines et al., 2018a, b).

The atmosphere of the super Earth GJ 1214b is one of the most extensively simulated atmospheres, using a range of different 3D models (Menou, 2012; Kataria et al., 2014; Charnay et al., 2015a; Drummond et al., 2018a) often combining radiative transfer calculations with equilibrium chemistry. Zhang & Showman (2017) took a more simplified modeling approach but explored the effect of a wide–range of bulk compositions. Charnay et al. (2015b) included a simplified treatment of condensation of KCl and ZnS clouds.

Considering only the solution to the dynamical equations, much progress has been made by GCMs solving the simplified or “primitive” equations of motion (invoking hydrostatic balance, the assumption of a thin or shallow atmosphere and a gravity constant with height, see discussion in Mayne et al., 2014a), such as derivatives of the MIT–GCM (with the main hot Jupiter adaptations presented in Showman et al., 2009) and the LMD$z$ (used, for example, in Charnay et al., 2015b). Concerns over the validity of the “primitive” equations for “thick” atmospheres, where the atmospheric scale height becomes significant compared to the assumed total planetary radius have been raised for Titan and Venus (see Tokano, 2013). Additionally, work has been done assessing the limits of the primitive equations for Earth (see Tort et al., 2015, and references therein), and the traditional approximation specifically (Gerkema et al., 2008). Several 3D models have been applied to exoplanets which solve the full Navier–Stokes equations, such as the 3D radiation–hydrodynamics model of Dobbs-Dixon & Agol (2013)111This code explicitly solves for the Centrifugal component, which is usually included in the gravitational term for GCMs (Showman et al., 2008), the recently developed THOR dynamical core (Mendonça et al., 2016), and of course our own work (e.g., Mayne et al., 2014b). However, significant differences between simulations of exoplanets using “primitive” and more complete dynamical cores has yet to be found.

In our own work we have applied the Met Office Unified Model (or UM) to exoplanets (first in Mayne et al., 2014a), which is able to solve the dynamical equations adopting various levels of simplification within the same numerical framework. In Mayne et al. (2014b) we applied this capability to hot Jupiters, finding differences in the evolution of the deep atmosphere between the primitive and more complete equation set, but not in the qualitative dynamics of the upper, observable atmosphere.

In this paper we continue our work using the UM and investigate the effect of the simplifications made to the dynamical equations in the Super Earth or Small Neptune regime, for example, GJ 1214b, or $\pi$ Mensae c (HD 39091c) the latter being recently detected using the Transiting Exoplanet Survey Satellite (TESS) (Huang et al., 2018). Focusing our work on GJ 1214b we find significant differences in the resulting dynamical structure of the simulated atmosphere between the “primitive” and more complete equations. This is caused by the breakdown of the ‘shallow–fluid’ and traditional approximations, similar to that found for a terrestrial planet by Tort et al. (2015). The rest of this paper is structured as follows: in Section 2 we introduce the model we are using, summarise its evolution and describe the setups and parameters adopted. In Section 3 we present results from the simulations, starting with a “standard” setup for GJ 1214b (Section 3.1), followed by an exploration of the limits of the primitive equations (Section 3.2, where Section 3.2.3 details the derivation of order of magnitude estimates for the wind speeds, which are verified for our simulations in Appendix B). We then discuss the limitations of our work in Section 4, and conclude in Section 5.

2 Model Description

For this study we use the Met Office Unified Model (UM) which we have adapted to enable modelling of a variety of planets and exoplanets. Initial adaptations and the implementation of a Newtonian relaxation, or temperature forcing scheme where the temperature is relaxed to a prescribed radiative equilibrium temperature over a specified timescale, are detailed in Mayne et al. (2014a), followed by adaptations for gas giants or planets with extended atmospheres in Mayne et al. (2014b). The radiative transfer component was adapted in Amundsen et al. (2014), and further in Amundsen et al. (2017), with simulations including full radiative transfer used in Helling et al. (2016) and Amundsen et al. (2016). The model chemistry options include an analytical chemical equilibrium scheme (Burrows & Sharp, 1999; Amundsen et al., 2016) as well as both a Gibbs energy minimisation scheme (Drummond et al., 2018a) and chemical relaxation scheme (Drummond et al., 2018c, b). We have also recently coupled a chemical kinetics scheme (originally developed in 1D Tremblin et al., 2015; Drummond et al., 2016) to the UM which will be the focus of a future work. The same model has been coupled to a ‘prognostic’ cloud scheme (Lines et al., 2018a, b), and used to explore the evolution of the deep atmosphere of hot Jupiters, combined with a 2D code (Tremblin et al., 2017), alongside exploration of the dynamical acceleration mechanisms (Mayne et al., 2017). Finally, we have also adapted the surface schemes to model terrestrial exoplanets (Boutle et al., 2017; Lewis et al., 2018).

In this study we adopt a temperature forcing scheme (implementation described in Mayne et al., 2014a) where the temperature is forced to a prescribed equilibrium temperature, $T_{\rm eq}$, over a parameterised timescale, $\tau_{\rm rad}$. We adopt the $\tau_{\rm rad}$ and $T_{\rm eq}$ used by Zhang & Showman (2017, their Eqn. 8 and 9, respectively), where $T_{\rm eq}$ is a function of longitude, latitude, and pressure, and $\tau_{\rm rad}$ a function of pressure only. The dayside–nightside temperature contrast is controlled by the parameter $\Delta T_{\rm eq}$, also a function of pressure ($p$), which Zhang & Showman (2017) set at 600 K at the top of the atmosphere, decreasing linearly with $ln\,p$ toward zero at the bottom of the atmosphere. In our height–based model we set $\Delta T_{\rm eq}=600$ K where $p\leq 10$ Pa and decrease linearly with $ln\,p$ to $\Delta T_{\rm eq}=0$ K where $p\geq 200\times 10^{5}$ Pa. In Appendix A we show that our conclusions are independent of the choice of this temperature forcing scheme and hold when using a more complete radiative transfer scheme. The planet is effectively modeled as a gas giant, as the inner boundary does not include a surface treatment of the land or ocean and radiative absorption/emission. However, the inclusion of a terrestrial surface, in the case of vertically extended atmospheres, will not likely change our main results. As in the study of Mayne et al. (2014b) no ‘drag’ is applied near the bottom boundary, as it is not required for stability and if present (due to, for example, magnetic drag in the case of a gas giant or surface friction for a terrestrial planet) the form and magnitude are poorly constrained.

The main equations solved by the UM include the conservation of momentum in each of the three directions, longitude, latitude and height, coordinates $\lambda$, $\phi$ and $r$, respectively, with the corresponding wind coordinates $u$ (zonal), $v$ (meridional) and $w$ (vertical), respectively. Additionally, the thermodynamic equation is solved adopting potential temperature, $\theta$ and the Exner Function, $\Pi$, where $\Pi=\left(\frac{p}{p_{0}}\right)^{R/c_{p}}=\frac{T}{\theta}$ and $c_{p}$ is the specific heat capacity, $R$ is the specific gas constant, $p_{0}$ is a chosen reference pressure and $T$ the normal temperature in Kelvin. This equation set is then closed via the ideal gas equation, and heating injected into the system using the Newtonian relaxation or temperature forcing scheme via the energy equation (see, for example Mayne et al., 2014a, b). The equations are detailed in full in several works, in particular Mayne et al. (2014b). As described in Mayne et al. (2014b) and Mayne et al. (2017) a vertical ‘sponge’ layer is used to damp vertical velocities near the upper, low pressure boundary, and a diffusion scheme is applied to remove grid–scale noise and aid numerical stability. SI units are used throughout the code, and adopted throughout this manuscript aside from when explicitly stated.

The UM, as applied to exoplanets in Mayne et al. (2014b), is capable of solving varying levels of the dynamical equations within the same model framework, from the most simplified ‘primitive’ equations which assume the atmosphere is in vertical hydrostatic equilibrium, is a ‘shallow–fluid’ and gravity is constant with height, to the ‘full’ equations, which do not make these same assumptions. For this work, we adopt the nomenclature of Mayne et al. (2014b), i.e. “primitive” (most simplified) and “full” (least simplified), and refer the reader to Mayne et al. (2014b) (and Section 3.2) for the detail of the exact terms within the equations used. However, as discussed in Mayne et al. (2014b), the “full” equations still include approximations. Notably the assumption of a spherical geoid, which is valid for small values of $R_{\rm p}\Omega^{2}/g$, where $R_{\rm p}$ and $\Omega$ are the planetary radius and rotation rate, respectively (Bánard, 2014). For our simulations of large radii, slow–rotating planets, this parameter is small (see Section 3.2). However, on Earth it can be larger and errors associated with this approximation can be comparable to those associated with the ‘shallow–fluid’ and traditional approximations (Bánard, 2014).

Model parameters for our standard simulation are shown in Table 1 where the planet specific parameters correspond to those for GJ 1214b Carter et al. (2011). The parameters relating to numerical settings are the same as those used in previous studies using the same model (see Mayne et al., 2014b; Amundsen et al., 2016; Mayne et al., 2017). We note that while our model setup specifically corresponds to GJ 1214b we expect our results to be representative of other highly–irradiated, Super Earth/Warm Neptune atmospheres with similar properties.

In addition to the standard simulation of GJ 1214b, we design eight additional setups that investigate the effect of varying parameters such as the planet radius ($R_{\rm p}$) and rotation rate ($\Omega$), amongst others; the motivation for each of these setups will be described in a later section. Each setup is described in Table 2, where we show the parameters that are adjusted from the standard ones, shown in Table 1. We also give each simulation setup a short name, which we use to refer to the setups in the text. For each setup we present simulations using both the full and primitive equation sets, while for one setup (dT$+$) we also present a simulation using the deep equation set, giving a total of 15 simulations in this work.

Each simulation is initialised at rest and in hydrostatic balance, with a zonally and meridionally homogeneous temperature–pressure profile. The substellar point is set at 180*∘* longitude. For the initial temperature profile we use the mean temperature ($T_{0}$) profile from Zhang & Showman (2017, see their Fig. 2). Each simulation was run for 1000 days (throughout this work days refers to Earth days). For all simulations, the maximum zonal wind velocity and the mean zonal wind structure has ceased to evolve after a few hundred days, reaching a pseudo–steady state. The deep, high pressure atmosphere continues to slowly evolve after 1000 days, however as shown for hot Jupiters by Mayne et al. (2017), this does not appear to have a significant effect on the flow in the lower pressure regions.

The pressure–altitude structure will vary between the full and primitive (or deep) equation sets due to the assumption of a gravity that is constant with height in the latter. In addition, the pressure–altitude structure also depends on various model parameters (e.g. $R$, $\Delta T_{\rm eq}$). Since our model is height–based, we adjust the top–of–atmosphere height ($z_{\rm top}$) to achieve a similar pressure range between each simulation. Each simulation encompasses the approximate pressure range $\sim$200$\times 10^{5}$ Pa to $\sim 10$ Pa. This means that the vertical resolution in height will be slightly different between different simulations though the vertical resolution in pressure will be approximately the same. However, these differences are typically small ($\sim$20%), and have no effect on our conclusions (see test and discussion in Section 3.2.3). The value of $z_{\rm top}$ used for each simulation is shown in Table 2.

3 Results

In this section, we first present results from our baseline or standard simulations (Std, see Table 2), with parameters matching those derived from observations of GJ 1214b (Section 3.1). These simulations demonstrate a clear difference in the resolved atmospheric structure between the more simplified primitive equations and their more complete counterparts. Therefore, we follow this with an exploration of the fundamental limits of the primitive equations, and a demonstration of the effects of exceeding these limits on the simulated flows using our remaining simulation set (Section 3.2).

The majority of simulations of tidally–locked atmospheres return a prograde, superrotating equatorial jet i.e. a coherent zonal flow in the direction of the planetary rotation peaking in speed towards low latitudes (see discussion in Showman & Polvani, 2011; Tsai et al., 2014; Mayne et al., 2017). Therefore, our analysis of the bulk dynamical structure of the atmosphere is performed via comparison of the zonal–mean, temporal–mean, zonal wind as a function of latitude and pressure (using linear interpolation to convert quantities from a height to pressure surfaces, as in Mayne et al., 2014b, 2017). For tidally-locked planets in close orbits, with extended atmospheres, the thermal structure in the upper, low pressure region of the atmosphere is dominated by the radiative forcing. In the temperature forcing setup, this is the regime of a very short radiative timescale, where the temperature is rapidly relaxed to the equilibrium temperature (Iro et al., 2005). As the radiative timescale increases with pressure the advection can begin to alter the temperature structure, and drive it from radiative equilibrium (see for example, Showman et al., 2009; Mayne et al., 2014b; Zhang & Showman, 2017). Therefore, changes in the advection between simulations will lead to effects on the thermal structure effectively weighted by the depth at which the flow is occurring. Basically, faster flows from the day to night side will lead to stronger homogenisation of the zonal temperature structure which, in turn, will have an increasing effect with increasing pressure.

To explore differences in the temperature structure between simulations we follow the approach of Zhang & Showman (2017) and present simple, normalised thermal phase curves. This is simply the blackbody thermal emission from an isobaric surface integrated over the observable hemisphere as a function of phase (see Eq. (18) of Zhang & Showman, 2017). The phase curve for a given pressure level is then normalised by (Zhang & Showman, 2017, their Eq. (18)),

[TABLE]

where $\bar{F}\left(\delta\right)$ is the normalised phase curve, $F\left(\delta\right)$ is the emission flux as a function of phase ($\delta$) and $\langle F\left(\delta\right)\rangle$ is the mean of the flux over the entire phase. We note that the normalised phase curve is negative where $F\left(\delta\right)<\langle F\left(\delta\right)\rangle$.

We use these simple phase curves as a tool to explore trends in the dayside–nightside temperature contrast as well as the longitudinal offset of the hot spot. We stress that these phase curves are not intended to represent the real emission from the planet, which would require a full radiative transfer approach to capture the pressure and wavelength dependence of the emission flux.

3.1 Standard Case

Figures 1, shows the zonal–mean, temporal–mean, zonal wind as a function of latitude and pressure, for the Std Prim, Std Full, and their high spatial resolution counterparts (Std Hires Prim and Std Hires Full, respectively, where the horizontal and vertical resolution have been doubled). Note, for these and subsequent similar figures the contour lines for all figures are the same, but the colour scales are varied between simulation ‘pairs’ i.e. each matching pair of primitive and full simulations have the same colour scales. Figure 1 shows a clear difference between the simulation using the primitive equations and the more complete full versions (top panels). The zonal prograde flow is decelerated, and spread across a larger latitude range and shallower pressure/height range for the simulation using the full equations. The difference is marked with a $\sim 1.5\times$ increase in the maximum zonal wind speed, and significant prograde zonal velocities penetrating an order of magnitude deeper in pressure for the primitive case, but ‘sharpened’ to a significantly more peaked latitudinal profile. This difference is also recovered in the higher spatial resolution simulations (middle panels), demonstrating that this is not an issue of poor resolution or resolution differences between simulation pairs (to the author’s knowledge these are the highest resolution GCM simulations published for an exoplanet atmosphere to date). For this work, all of our simulations have been run using a simplified temperature forcing scheme to model the atmospheric heating. However, as shown in Appendix A our conclusions are unchanged when moving to the use of a more sophisticated radiative transfer scheme.

Figure 2 shows the temperature (colour scale) and horizontal wind (vector arrows) at isobaric surfaces (in $\lambda$ and $\phi$) of 100 and 3 000 Pa at 1 000 days for the Std Full simulation (top panels), and the differences in these fields for the Std Prim simulation (middle panels, where all differences are in the sense full simulation minus primitive simulation). Additionally, the simplified thermal phase curves for both simulations at the two depths are shown (bottom panels).

Figure 2 shows significant changes in the thermal structure, and subsequent phase curve between the primitive and full equations for our standard setup. The changes in the zonal flow (see Figure 1) clearly translate to alterations in the temperature structure. For the standard simulations at 100 Pa (left column, Figure 2), the full equations result in an increase in temperature near and to the east of the sub–stellar point (180*°*), with cooler temperatures to the west, effectively driven by a weaker superrotating jet, as can be seen from the difference in the horizontal velocity (vector arrows, and see Figure 1). This results, overall, in a warmer dayside and cooler nightside and, therefore, an increased amplitude in the simplified phase curve. For 3 000 Pa (right column), the regions flanking the equatorial jet, in latitude, and at the poles are warmed, with the most significant warming on the dayside. The jet, and mid-latitude regions are cooled, with the strongest cooling on the night side. The amplitude of the simplified phase curve is proportional to $T^{4}\cos\phi$, heavily weighted to the equatorial regions222We note that we show normalised phase curves (see Eqn. (1)) which means that relative differences in the amplitude of the normalised phase curve between simulations can be larger than might be expected based on the relative temperature differences.. Therefore, the overall amplitude is again increased, however here the peak of the curve is also shifted. As discussed the radiative timescale is longer at deeper pressures so changes in the flow can more easily affect the longitudinal temperature structure. The peak of the warming is close to the sub–stellar point, shifting the peak amplitude of the simplified phase curve back toward this point.

Clearly, for the case of a slowly rotating warm Small Neptune (or potentially Super Earth) such as GJ 1214b, care must be taken when interpreting results derived from simulations solving the primitive equations. As discussed there are several additional assumptions made when simplifying from the full to the primitive equations so it is important to isolate how each of them contributes to the differences in the resolved flow.

3.2 Limits of the Primitive Equations

In Mayne et al. (2014b) we detail the terms included for each equation set, but here we restate only the dominant zonal momentum equation to explore the impacts of the approximations made in the primitive equations.

The equation of zonal momentum conservation in the full equations is:

[TABLE]

where $\frac{D}{Dt}$ is the material derivative and D the diffusion operator (see Mayne et al., 2017). $f$ and $f^{{}^{\prime}}$ are the Coriolis parameters defined as,

[TABLE]

and

[TABLE]

The terms in Eqn. (2) can be ‘grouped’ into those associated with rotation, i.e. $f$ and $f^{\prime}$, so–called ‘metric’ or advective terms associated with flow along a curved surface, i.e. $\frac{uv\tan\phi}{r}$ and $\frac{uw}{r}$, and the pressure gradient term, here $\frac{c_{p}\theta}{r\cos\phi}\frac{\partial\Pi}{\partial\lambda}$. As discussed in Section 2 the full equations include the spherical–geoid approximation which is only justified when $R_{\rm p}\Omega^{2}/g$ is small (Bánard, 2014). For our simulations, Tables 1 and 2 show $R_{\rm p}\Omega^{2}/g$ is $\lesssim 0.01$ for all cases, and $\sim 0.0025$ for the standard case, meaning this approximation does not lead to significant errors. For faster rotating, and smaller radii planets this may, however, begin to introduce more significant errors (Bánard, 2014).

The key underlying assumption made in constructing the primitive equations is that the aspect ratio of the motion is small i.e. the flow in the vertical direction is significantly smaller in scale than that in the horizontal direction. This validates the three major assumptions (see, for example, Vallis, 2006; Mayne et al., 2014b). (1) The hydrostatic approximation: the atmosphere is in vertical hydrostatic balance, regulated by fast acoustic waves, (2) The shallow–fluid approximation: the atmosphere is thin in relation to the size of the planet meaning that $r$ can be replaced with $R_{\rm p}$ and $\partial/\partial r$ can be replaced by $\partial/\partial z$, and (3) The traditional approximation: where the smaller metric and Coriolis terms are neglected. In the case of the zonal momentum equation, for example, this requires that $w\ll v\tan\phi$ such that the term $\frac{uw}{r}$ is ensured to be negligible when compared to the term $\frac{uv\tan\phi}{r}$ (see Eqn. (2) and Mayne et al., 2014b). This condition also ensures that the term $2\Omega w\cos\phi$ is small compared to $2\Omega v\sin\phi$ and can similarly be dropped. The assumption of a shallow–fluid and the traditional approximation must be taken together to ensure that the equation set conserves angular momentum (White & Bromley, 1995), and also allow the further assumption; (4) gravity is constant with height. The approximations of hydrostasy and a gravity constant with height directly impact the vertical momentum equation, and indirectly impact the zonal momentum (see Mayne et al., 2014b, for explicit equations). However, the shallow–fluid and traditional approximations directly impact the zonal (and meridional) momentum equation, evident when comparing the zonal momentum equation for the primitive case,

[TABLE]

with the full equation (Eq. (2)). Specifically, for Eqns (2) and (5), the shallow–fluid approximation results in $r$ being replaced with $R_{\rm p}$ for the first (metric) and final (pressure gradient) terms. The traditional approximation then leads to the omission of the second (metric) $\frac{uw}{r}$ term and $f^{\prime}$ (or $2\Omega w\cos\phi$, a rotational or Coriolis term) as discussed.

In this section we review each of the assumptions made in constructing the primitive equations, for the regime of warm small Neptunes (or Super Earths), and explore manifestations of their limitations using the simulation set presented in Table 2.

3.2.1 The Hydrostatic Approximation

Clearly, hydrostasy is enforced in the primitive simulations, however for all our simulations solving the full equations vertical hydrostatic equilibrium remains a good, global, approximation. Small departures from hydrostatic balance are found in small regions of the highest, low pressure, atmospheric layers, but do not significantly alter the bulk flow. Typically, hydrostasy will hold for flows with a horizontal scale larger than the vertical scaleheight, $H=RT/g\sim 1.5\times 10^{5}$ m, in our case, which is less than 0.2% of the planetary circumference.

3.2.2 The Shallow–fluid Approximation

The shallow–fluid approximation is evidently contingent on the condition that the height above the inner boundary of the dominant zonal flow, or the vertical extent of the dynamically active atmosphere ($z_{U}$) is much smaller than the planetary radius i.e. in order to replace $r$ with $R_{\rm p}$, where $r$ can be expressed as $R_{\rm p}+z_{U}$ we require $z_{U}\ll R_{\rm p}$. For a Super Earth $R_{\rm p}$ is fixed as it has a solid surface, but for a small Neptune $R_{\rm p}$ could simply be shifted to lower pressures if regions of the deep atmosphere were quiescent and not affecting the overall dynamics. However, from our simulations we see that the flow extends throughout the majority of our modelled domain, and the inner boundary could not be raised significantly without affecting the flow. In the Std Full and Std Prim cases, shown in Figure 1, the vertical extent of the dynamically active atmosphere is $\sim 20\%$ of the planetary radius. We have performed simulations matching the standard setup yet with an increased, or decreased planetary radius ($R_{\rm p}+$ and $R_{\rm p}-$, respectively, see Table 2), to explore this effect. Additionally, we perform a simulation adopting the standard planetary radius, but with a CO2 dominated atmosphere, following Zhang & Showman (2017), leading to an increased mean molecular weight, and subsequent vertical compression of the atmosphere itself.

Figures 3 shows the same information, in the same format as Figure 1 but for the $R_{\rm p}-$, $R_{\rm p}+$ and CO2 simulations. Clearly, for the $R_{\rm p}-$ simulation differences remain between the primitive and full simulations, whereas these differences are almost completely removed when the planetary radius is increased ($R_{\rm p}+$). They are also significantly reduced when the atmospheric scaleheight is reduced (CO2) compared with the standard setup. However, some differences remain, which are discussed in Section 3.2.3. It is important to note that due to the assumption of a constant gravity in the primitive equations, as opposed to reducing proportional to $(R_{\rm p}/r)^{2}$ in the full case, differences in $g$ between the primitive and full cases will be reduced as the planetary radius is increased, or atmospheric scaleheight reduced. This could be affecting our results and contributing to the increasing similarity in the resolved flow between the full and primitive cases for the $R_{\rm p}+$ and CO2 simulations. However, we show in Section 3.2.3 that for the standard case the difference in gravity between the primitive and full case does not significantly affect the flow (see Figure 8), where this differences is, at its maximum, $\sim 20$% (compared to a maximum of $10$% in the $R_{\rm p}+$ simulation).

For the $R_{\rm p}+$ the temperature changes between the primitive and full simulations, at both 100 and 3 000 Pa, are generally less than $\sim 2$ K. Therefore, Figures of the temperature structure for these simulations are omitted. In the case of the CO2 simulation, the differences in the zonal wind structure resolved using the primitive or full equations are reduced, over the standard setup, but are not negligible. However, the differences between the two flows does reduce towards lower pressures. Figures 4 and 5 shows the same information as Figure 2 but for the $R_{\rm p}-$ and CO2 simulations.

The pattern of the change, at 100 Pa, when moving to the more complete equations of motion for the $R_{\rm p}-$ case is similar to our standard GJ 1214b setups, but enhanced. The regions east of the sub–stellar point are heated and those to the west are cooled, with a net day side heating and night side cooling. For the higher pressure, 3 000 Pa, the pattern is slightly different, with heating close to the poles to the west of the substellar point, and east of the substellar point at the equator. Again, as the simplified phase curve amplitude is weighted by $\cos\phi$ we have an increase in the peak day side amplitude, and reduction in the night side minimum. The offset between the peak at the substellar point also reduces, as for the standard simulations, caused by overall less efficient redistribution of heat in the more complete equation case. This is shown, at 100 Pa, by a net reduction in the prograde or superrotating flow. For the CO2 simulation, in Figure 5, more significant temperature changes are found at 100 Pa, but as these occur very close to the poles they contribute negligibly to the resulting simplified phase curve, and the full and primitive versions closely match. For the deeper pressure level, the absolute changes are smaller, but distributed closer to the equator resulting in a slight shift in the simplified phase curve. The change in the horizontal wind, shown by the vector arrows, is again a reduction in the strength of the superrorating jet at the equator, although some more complex changes occur at high latitudes for the lower pressure surface of the CO2 simulation.

The simulations in this section demonstrate that the flow recovered from simulations adopting either the primitive or full equations of motion differ significantly for ‘thick atmospheres’, where $z_{U}/R_{\rm p}\sim 20\%$ or more. The limitation of the shallow–fluid approximation, within the primitive equations, is well known, but here we demonstrate that an important class of exoplanet, warm small Neptunes or Super Earths, potentially inhabit this problematic regime.

3.2.3 Traditional

Approximation

As discussed in Mayne et al. (2014b) the physical justification for the traditional approximation is relatively weak, and concerns have been raised over its validity for thick atmospheres (Tokano, 2013; Mayne et al., 2014b), and terrestrial planets (Gerkema et al., 2008; Tort et al., 2015). Also as discussed the condition $w\ll v\tan\phi$ indicates whether the terms omitted by this approximation are negligible or not. It is hard to diagnose this condition a priori as we don’t know the magnitude of the winds. However, this condition can clearly never be met at the equator where $\tan\phi=0$. Therefore, for this approximation to hold, and the primitive equations to correctly capture the flow, the equatorial region where $w\gtrsim v\tan\phi$ must be sufficiently restricted, and connected smoothly to the flow at mid–latitudes where the condition can more easily be met.

In order to develop an order of magnitude estimate for the typical $v\tan\phi$ and $w$ values we invoke several assumptions, before exploring the results from the simulations themselves. The resulting equations are highly simplified, but allow us to explore the approximate behaviour of the atmosphere and understand the simulation results. Firstly, we assume that our standard case is a reasonable approximation of the global dynamic and thermodynamic structure of GJ 1214b and similar warm small Neptunes or Super Earths. Given this assumption, within the dynamically ‘active’ region from $10^{2}-10^{5}$ Pa (or 1 mbar to 1 bar), it can reasonably be assumed that:

• •

The atmosphere is globally superrotating: this always holds at the equator, and in most cases is also true elsewhere in the atmosphere.

• •

The longitudinal variation in temperature is smaller than the equilibrium day–night temperature contrast (i.e. the jet is acting to redistribute heat and homogenise the temperature).

• •

The meridional winds are predominantly driven by the Coriolis effect on the zonal wind, which follows from our previous two assumptions.

• •

The gas is incompressible, and in vertical hydrostatic equilibrium where the former is reasonable as the flow speed is negligible compared to the sound speed for the majority of the atmosphere (see Section 4), and as we have discussed the former holds for all but the lowest pressure regions (Section 3.2.1).

For an inviscid incompressible atmospheric flow in a steady state the Euler, mass continuity, energy and ideal gas equations can be written as:

[TABLE]

and

[TABLE]

respectively. Where $\vec{v}$, $\vec{\nabla}$, $\vec{\Omega}$ and $\vec{g}$ are merely the vector forms of the variables and operator previously defined. Using the vector identity $\vec{\nabla}(v^{2}/2)=(\vec{v}\cdot\vec{\nabla})\vec{v}+\vec{v}\wedge(\vec{\nabla}\wedge\vec{v})$ and taking the scalar product of Eqn.(6) with $\vec{v}$ we obtain

[TABLE]

This equation simply shows that the advection of kinetic energy is balanced by the advection of pressure and bouyancy in the atmosphere.

For an incompressible gas $\vec{\nabla}\cdot\vec{v}=0$ and therefore $\dfrac{1}{r\cos\phi}\dfrac{\partial u}{\partial\lambda}+\dfrac{1}{r}\dfrac{\partial v}{\partial\phi}+\dfrac{\partial w}{\partial r}=0$ suggesting the standard order of magnitude estimate of

[TABLE]

where $U$ is an order of magnitude for $u$, $W$ for $w$, $L_{U}$ is a characteristic horizontal length for $u$ and $H_{W}$ a characteristic scale height for $w$. Hydrostatic balance then implies that the vertical component of the advection of $p$ is balanced by $wg$,

[TABLE]

where $\vec{\nabla_{\bot}}$ is the horizontal gradient. The second order of magnitude equality arises from the fact that $U$ is dominating the other winds, and therefore we do not expect the total advection to differ strongly from the horizontal advection. We have verified this point in our simulations. Further by substituting for $p$ using Eqn. (9) using (8) we obtain

[TABLE]

Combining Eqns.(10), (12) and (13) yields

[TABLE]

which can be simplified to

[TABLE]

where $\Delta T$ is the order of magnitude of the difference between the actual temperature in the atmosphere and the corresponding temperature at radiative equilibrium. As we assume that the temperature variations with longitude are negligible compared to the equilibrium day–night contrast, we can further state that the steady–state temperature around the planet is approximately constant at $T\sim(T_{\mathrm{eq,day}}+T_{\mathrm{eq,night}})/2$, where $T_{\mathrm{eq,day}}$ and $T_{\mathrm{eq,night}}$ are the average day and night side equilibrium temperatures. Therefore, on average the difference between the steady state temperature and equilibrium temperature is $\Delta T\sim(T_{\mathrm{eq,day}}-T_{\mathrm{eq,night}})/2$.

For our simulations the characteristic length scale for the zonal jet is about half the planetary circumference or $\sim\pi R_{\rm p}$ and we expect $\Delta T\sim 0.5\Delta T_{\rm eq}$ (or half the total day–night contrast, see Section 2) which is $300$ K and $400$ K for the Std and dT$+$ simulations, respectively. At 100 Pa, where the maximum zonal velocity is found the radiative timescale is $\sim 10^{4}s$. Therefore, using these values in Eqn. (15) yields

[TABLE]

Our maximum zonal velocities are close to, or within approximately a factor 2 of this value which is a remarkable level of agreement given the approximations we have invoked to derive the estimated wind speeds, giving us confidence in our expressions. Eqn. (16) can then be substituted into Eqn. (11) to provide an estimate for $W$,

[TABLE]

where $H$ is now the typical height scale for $W$, $L$ the typical horizontal scale for the typical zonal velocity, $U$ (i.e. we have dropped the subscripts used in Eqn. (11)).

To estimate $V$, as stated we assume that the meridional motions are primarily driven by the Coriolis force. The meridional component of the momentum equation, Eqn. (6), assuming that the Coriolis term is dominant can be written as

[TABLE]

which can be approximated by

[TABLE]

where $L_{V}$ is a characteristic scale for $V$, and this expression can further be simplified to

[TABLE]

In our simulations which $L_{V}\sim\pi/4R_{\rm p}$, yielding

[TABLE]

Again this estimate matches the simulation results remarkably well to within a factor $\sim$2. A final, simple manipulation provides an estimate for $V\tan\phi$,

[TABLE]

which can be combined with Eqn. (17), to explore the global behaviour of the $w\ll v\tan\phi$ condition.

The equations we have derived (Eqns. (17) and (22)) to estimate the magnitude of the flow are valid in the upper atmosphere where the radiative timescale is short and temperature gradient high (and friction is low). This is also the region, as we have shown, where the zonal flow is dominant and the largest differences between the primitive and full equations of motion are found. Eqns. (22) and (17) are linked by a common assumption namely that the flow is dominated by the $U$ component which is balancing the pressure (temperature) gradient. Eqn. (22) then requires the further assumption that the $V$ component is driven, chiefly, by the Coriolis forces, where Eqn. (17) requires the assumption that the flow is divergence free. The derivation of these equations, however, does not rely on adopting either the primitive or full equations explicitly. In Appendinx B we explicitly verify these equations using our simulation outputs.

Although very simplified, without a proper prescription of $H$, $L$ or $\Delta T$, these Eqns. (17) and (22) allow us to draw insight into the nature of the flow. Firstly, as the rotation rate ($\Omega$) increases $V\tan\phi$ increases, strengthening the validity of the traditional approximation. This is the same effect observed in the simulations presented in Section 3.2.2. Essentially, the faster the rotation the more valid the primitive equations will become, all else being equal. The second insight provided by the order of magnitude estimates, Eqns (22) and (17), is that as the radiative forcing is increased i.e. $\Delta T$, which acts to increase the horizontal advection, the $W$ velocities will increase and the primitive equations become less valid. Additionally, the presence of the ideal gas constant in the estimate for $W$ shows that the composition of the atmosphere, and mean molecular weight will also play a role in the validity of the traditional approximation. Finally, the aspect ratio of the atmosphere, $H/L$, also appears in Eqn. (17), confirming the well known limit of the primitive equations to flows with small aspect ratios.

In addition Tort et al. (2015) have shown that the traditional approximation becomes increasingly less valid as, in their case, $\Omega$ is increased, yet $\Omega R_{\rm p}$ remains constant, for the terrestrial regime. If $\Omega R_{\rm p}$ is constant, Eqn (22) is clearly constant for a given $\phi$. Eqn (17) can be expressed as $W\propto\frac{R_{\rm p}^{1/3}}{L}$ which, using $L\sim R_{\rm p}$ becomes $W\propto\frac{1}{R_{\rm p}^{2/3}}$ which at constant $\Omega R_{\rm p}$ can be expressed as $W\propto\Omega^{2/3}$. Therefore, at constant $\Omega R_{\rm p}$ as $\Omega$ is increased $W$ increases, while $V$ is constant, meaning the condition $w\ll v\tan\phi$ becomes increasingly less valid.

Eqns. (22) and (17) also make interesting predictions for the atmospheric interaction with condensates. Larger vertical velocities are predicted for increasing levels of irradiation which could potentially lead to larger particles sizes being ‘lofted’ to high altitudes compared to less irradiated planets. Additionally, the meridional velocities increase with rotation, potentially leading to more meridional mixing of cloud material as found by Lines et al. (2018a).

In the remainder of this section, we explore the results from the simulations where the rotation rate, forcing or temperature contrast, mean molecular weight and gravity are varied demonstrating the regimes where the primitive equations become inaccurate.

Composition

results from the CO2 simulation, where the mean molecular weight of the atmosphere has been increased, have already been introduced in Figures 3 and 5, as part of the shallow–fluid discussion (Section 3.2.2). For the CO2 simulation the ratio of scaleheight of the atmosphere, to the planetary radius is similar to that of the $R_{\rm p}+$ simulation, being around 1% and 3%, respectively (see Tables 1 and 2). Therefore, the residual differences between the resolved primitive and full flows for the CO2 case, which are not observed to the same extent in the $R_{\rm p}+$ simulation, are likely caused by the change in the mean molecular weight itself, driving a change in the $W$ field, as show in Eqn (17). Figure 3 shows the peak velocities of prograde equatorial jet extending down to shallower depths, and therefore higher pressures in the full case, compared to that of the primitive case. Although the value of $R$ has changed by an order of magnitude moving from the standard simulation (compare values from Table 1 with those from Table 2), Eqn. (17) shows that $W\propto\sqrt{R}$, suggesting a less significant change in the flow is expected.

Rotation rate

Figure 6 shows the zonal–mean, temporal–mean zonal wind for the $\Omega+$ Prim and $\Omega+$ Full simulations. For these simulations the setup is the same as the standard case but with an increased rotation rate (see Table 2). As predicted both by our order of magnitude analysis the difference in the flow, between the primitive and full equations, is reduced when increasing the rotation rate, as the flow becomes dominated by Coriolis forces.

The differences in zonal flow between the primitive and full $\Omega+$ simulations (see Figure 6), result in changes to the temperature structure and simplified thermal phase curve which are shown in Figure 7. As for the standard case (Figure 2) a net warming is shown on the dayside and a net cooling on the nightside for the $\Omega+$ simulations, at both pressure levels, albeit at a reduced level when moving to the more complete equations. For the higher pressure, a very small shift in the peak amplitude of the simplified phase curve is found, as temperature changes at this level are around a few degrees.

Forcing/advection and gravity

Eqn (17) shows that as the temperature contrast is increased $W$ increases, leading to a weakening of the $w\ll v\tan\phi$ condition. Additionally, the acceleration due to gravity does not appear in our estimate of $V\tan\phi$, but does enter the estimate for $W$ via the scaleheight ($H\propto RT/g$) (Eqns (22) and (17), respectively). To demonstrate, and test the impact of this, we have run a set of simulations (dT$+$, see Table 2), matching the standard setup, but with an increased temperature contrast solving the primitive, deep and full equations of motion (see Mayne et al., 2014b, for explicit definition of equation sets). The deep and full equations differ only by the assumption of a gravity constant with height in the former.

Figure 8 shows the wind structure with an increased planetary temperature contrast (thereby increasing the zonal wind speed). The difference between the primitive and full simulations is enhanced when the temperature contrast is increased and the resulting zonal flow velocity increased, thereby moving the flow to a more advectively dominated regime (bottom panel). Figure 8(a) then shows the same setup as the dT$+$ Full and dT$+$ Prim cases but enforcing a constant gravity with height (termed the ‘deep’ equations, see Mayne et al., 2014b). Importantly, the resulting flow matches, more closely, the full equation case indicating that the changes are independent of the treatment of gravity. This is perhaps unsurprising given that the maximum variation in $g$ between the constant and varying case, found at the top of the atmosphere, is only $\sim 30$%, with the difference rapidly reducing with altitude, $\propto 1/r^{2}$. Additionally, the $z_{\rm top}$, and therefore vertical resolution in height is identical between the dT$+$ Prim and dT$+$ Deep, meaning that as the flow from the dT$+$ Deep simulation matches that of the dT$+$ Full simulation the vertical height resolution is not the cause of the differences between the primitive and full setups. Figure 9 then shows the temperature structure and phase curves, as in Figure 2 but for the dT$+$ case. A similar picture emerges, and it is important to note that the dT$+$ Deep simulation matches the dT$+$ Full simulation versions very closely as one would expect from Figures 1 and 8.

Figures 9 shows the resulting thermal structure, and simplified thermal phase curve for the dT$+$ Prim and Full simulations (the simulation solving the deep equations are omitted due to its similarity with the full case) in the same format as Figure 2. As in the standard case Figure 9 show significant changes in the thermal structure, and subsequent phase curve between the primitive and full equations for our setups with enhanced temperature contrast. The changes in the zonal flow (see Figure 8) clearly translate to alterations in the temperature structure, and the changes shown in Figure 9 are similar to those found in the standard setup (Figure 2).

In summary, the changes in the atmospheric dynamics are found when moving from solving the primitive to full dynamical equations for cases where the flow is strong (and rotation weak), and the atmosphere is thicker in relation to the planetary radius. These changes in the advection, result in changes to the temperature structure. Specifically, we have shown, using simplified phase curves, that moving from the primitive to full equations can affect the day-night temperature contrast and the location of the hot spot. This indicates a net change in the overall heat redistribution efficiency of these atmospheres when modeled using dynamical equations of differing simplicity. Essentially, the validity and applicability of the primitive equations of dynamics becomes questionable as the atmosphere becomes thicker (i.e. $z_{U}/R_{\rm p}$ becomes larger, shallow–fluid approximation) and $w\ll v\tan\phi$ is violated (the traditional approximation). In Appendix B, Figure 11 we show the direct violation of this limit from our simulations (namely the Std, $R_{\rm p}+$, dT$+$ and $\Omega+$ simulations). These conditions are well known, but we have demonstrated that an important class of exoplanet could well exist within the regime where the primitive equations are no longer valid namely warm/hot small Neptunes or Super Earths.

4 Assumptions and Limitations

The inner boundary is modelled as a friction free impermeable boundary with no heat or mass exchange (aside from the prescribed internal heat flux in the radiative transfer simulations presented in Appendix A, see Amundsen et al., 2016). Therefore, we have effectively modeled GJ 1214b as a gas giant planet, without a treatment of a potential rocky (or ocean) surface nor inclusion of a ‘drag’ (see discussion in Mayne et al., 2014b). However, as our inner boundary pressure is $200\times 10^{5}$ Pa, and the atmosphere vertically extended, the inclusion of a surface treatment is unlikely to alter our results as radiation does not significantly penetrate to the surface itself. Additionally, the near surface layers of the atmosphere are not flowing rapidly, and will therefore not be strongly affected by frictional effects. Furthermore, Mayne et al. (2014b, 2017) show that simulations without a bottom ‘drag’ return the same qualitative flow as those without. Additionally, Tort et al. (2015), have already shown that the traditional approximation weakens, for terrestrial planets, as $\Omega$ is increased but $\Omega R_{\rm p}$ is held constant, (see discussion in Section 3.2.3), where we have focused on varying the atmospheric extent.

Our model assumes that the self–gravity of the atmosphere is negligible, i.e. the mass of the modeled atmosphere, $M_{\rm atm}$ is negligible compared to the mass of the bulk planet comprising the unmodeled material within the inner boundary, $M_{\rm p}$. For our basic setup the adopted surface gravity, $g_{\rm p}=$12.20 ms*-2* and planetary radius, $R_{\rm p}=1.45\times 10^{7}$ m imply $M_{\rm p}\sim 3.8\times 10^{25}$ kg and the standard simulations indicate $M_{\rm atm}\sim 5.1\times 10^{21}$ kg, suggesting that this is a reasonable assumption. For the purposes of this study we have simply altered the planetary radius without a commensurate change in the gravity, while the atmospheric extent has also been altered to preserve a similar pressure range between simulations. For the reduced and increased planetary radii simulations $M_{\rm atm}\sim 9.7\times 10^{20}$ and $2.3\times 10^{23}$ kg, respectively, and the assumed planetary masses (as the surface gravity is not altered) are, $M_{\rm p}\sim 6.58\times 10^{24}$ and $1.83\times 10^{27}$ kg, respectively, indicating that neglecting self–gravity remains a reasonable approximation for all simulations. Additionally, although vertically extended the total atmospheric mass, in our modelled domain ($\sim 200\times 10^{5}$–$10$ Pa), remains a plausible fraction of the total planetary mass, i.e. less than 1% (Elkins-Tanton & Seager, 2008; Lopez & Fortney, 2014).

As discussed the vertical extent of our atmosphere is set to retain a similar pressure range across all simulations, corresponding to $\sim 200\times 10^{5}$ Pa to $\sim 10$ Pa, allowing both inter–comparison of our simulations, and comparison with the results of Zhang & Showman (2017). The model calculates hydrostatic balance for an input temperature–pressure profile, which is derived using a 1D radiative transfer code (selected as on the mean $T_{0}$ profile from Zhang & Showman, 2017), an initial inner boundary pressure (see Table 1) and the selected height. The selected height is set so that a similar pressure range for all simulations is achieved. The pressure structure will then evolve, including that at the inner and outer boundary, however the evolution does not result in a significantly different maximum or minimum pressures. This means that the height of our atmosphere is effectively set by the gravity and inner boundary pressure we select, where the gravity is approximately consistent with the measured planetary mass.

One of our main conclusions rests on the modeled atmospheric height becoming comparable to the defined planetary radius. Therefore, could the differences be erased by selecting a different planetary radius and atmospheric height? For our simulations the inner boundary is placed at pressures high enough so as to capture the entire ‘dynamically active’ region of the atmosphere. The atmospheric flow we are simulating often extends across the entire pressure or height range we simulated, with only small quiescent regions close to the inner boundary, most importantly for the Std Full and Std Prim cases. This means that restricting our height range will likely alter the dynamics, and not capture the atmospheric flow correctly. Our inner boundary pressure, chosen to be consistent with Amundsen et al. (2016) is also consistent with the typical radiative–convective boundary for small Neptune sized planets at a few Gyrs which is $10^{7}-10^{8}$ Pa (Lopez & Fortney, 2014). The planetary radius is usually measured using the optical transit, and can be estimated to be the radius at 20 mbar or $2\times 10^{3}$ Pa (Hubbard et al., 2001). Here we have set this radius as our inner boundary at a pressure of $\sim 200\times 10^{5}$ Pa, with the standard simulation reaching $2\times 10^{3}$ Pa at $\sim 2\times 10^{6}$ m. Therefore, our actual $R_{\rm p}$ could plausibly be reduced from $1.4\times 10^{7}$ to $\sim 1.2\times 10^{7}$ m for the standard setup i.e. fixing the radial distance out to a pressure of $2\times 10^{3}$ Pa as the measured value of $1.4\times 10^{7}$ m, and therefore, deriving a new value for our inner boundary location such that the height above the inner boundary of $2\times 10^{3}$ Pa matches our simulation value ($\sim 2\times 10^{6}$ m.). This reduction in assumed planetary radius, or radial distance to the inner boundary would act to enhance the differences we have found between the primitive and full equations.

An estimate of the sound speed is $c_{\rm s}\sim\sqrt{\gamma RT}$, with $\gamma=1.5$ and temperatures in the jet region of $\sim 500$ K this yields $c_{\rm s}\sim 1600$ ms*-1* for our simulations ($\sim 400$ ms*-1* for the CO2 simulations). Therefore the approximate Mach number in our simulations is less than unity throughout the atmosphere except, in some cases, the jet core, where it can reach $\sim 1.5$ in isolated regions. It is possible shocks might occur in these regions, dissipating energy from the flow, a mechanism which is not captured in our model. However, this study is concerned with the relative changes in flow between the primitive and full simulations, and the maximum wind speed for our simulations can effectively be set by the level of dissipation, which is a free parameter. Additionally, studies of the effect of shocks in the atmospheres of hot Jupiters, sharing several characteristics with the targets of this paper, have shown their effects to be negligible on the large–scale dynamics of the atmospheres (Fromang et al., 2016).

5 Conclusions

Modeling the 3D dynamics of planetary atmospheres can be performed at varying levels of simplification. Much progress has been made using models adopting the primitive equations, but more recent results have been derived using models adopting the more complete “full” equations (following the nomenclature of Mayne et al., 2014b).

In this work we have demonstrated with a set of self–consistent simulations that the dominant zonal flow recovered adopting the primitive or full dynamical equations differs markedly for GJ 1214b. Additionally, we have demonstrated that the change in the zonal–mean, zonal flow results in changes to the thermal structure, via a change in the efficiency of the heat redistribution. The day-night temperature contrast and location of the hot spot differ markedly between pairs of simulations, using the simplified and more complete dynamical equations for several cases. In other words the simplifications involved in the derivation of the primitive equations begin to break down, specifically the so–called ‘shallow–fluid’ approximation where the atmosphere is assumed to be small in vertical extent compared to the radius of the bulk planet, and the traditional approximation. We have shown that these differences are most apparent for slowly rotating planets, where advective terms dominate over the Coriolis terms, i.e. planets in short–period orbits which are strongly heated driving fast winds. Here we stress that our standard simulations (i.e. Std Full, Std Prim, Std Hires Full & Std Hires Prim) which are setup to be physically based on GJ 1214b, clearly exhibit this effect. Our remaining simulations which are less physically–based are then used to explore the parameter space.

Although the validity of the ‘shallow–fluid’ and traditional approximations is not a new consideration (see Gerkema et al., 2008; Tort et al., 2015, and references therein in the context of Earth), we have shown that for an exceedingly important class of exoplanet i.e., Super Earths or small Neptunes, with potentially large, vertically extended atmospheres in tidally–locked orbital configurations, the dynamical solutions obtained from treatments adopting the primitive equations are likely to become inaccurate. Indeed our conclusions complement those of Tort et al. (2015) who find the traditional approximation weakens for faster rotating terrestrial planets on the condition that $\Omega R_{\rm p}$ is held constant (see Section 3.2.3). Importantly, our numerical setups are practically identical between pairs of simulations solving the primitive and full equations. The fact that the hydrostatic approximation appears to hold for the atmospheres we have simulated also suggests that the quasi–hydrostatic deep atmosphere–equations (e.g., White & Bromley, 1995; Tort et al., 2015) should also yield comparable results to the ‘full’ equations.

Given the potential for these planets to show departures from chemical equilibrium driven by the advection (Madhusudhan et al., 2016; Venot et al., 2017), inaccuracies in the derived dynamical structure of the atmosphere could have significant impacts on the inferred underlying chemical and thermodynamic structures. This class of objects is also important as it represents the regime bridging giant planets with terrestrial planets (Lopez & Fortney, 2014).

The next steps for this study would be to move to a full radiative transfer solution (as used in Amundsen et al., 2016), relax the assumption of chemical equilibrium (Drummond et al., 2018c, b) and include the cloud scheme adopted in Lines et al. (2018a). These next steps would allow us to investigate the effect of the differences in the dynamics on both the chemistry and cloud structures, but also to derive more accurate synthetic observations for comparison with real data.

NJM is part funded by a Leverhulme Trust Research Project Grant and gratefully acknowledges their support. BD acknowledges funding from the European Research Council (ERC) under the European Unions Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement no. 336792. JM and IAB acknowledge the support of a Met Office Academic Partnership secondment. Material produced using Met Office Software. This work used the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility (www.dirac.ac.uk ). This equipment is funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National E-Infrastructure. This research made use of the ISCA High Performance Computing Service at the University of Exeter. We acknowledge the comments from the anonymous reviewer which improved this manuscript.

Appendix A Radiative Transfer

All of our simulations analysed in this work, and listed in Table 2 are performed using temperature forcing (see Mayne et al., 2014a, for implementation). Such temperature forcing schemes have proved exceedingly useful in analysing the flow regimes of planetary atmospheres, and are extremely computationally efficient. However, this approach is certainly less physically accurate than the use of more complete radiative transfer scheme (see discussion in Showman et al., 2009; Amundsen et al., 2014, 2016). To ensure that our conclusions are not affected by the use of a simplified treatment of the heating, we have run two additional simulations for a shorter period of 500 days employing the full radiative transfer scheme detailed in Amundsen et al. (2014, 2016). These simulations are setup to match the Std Prim and Std Full simulations, with the complete setup (including radiative transfer elements) detailed in Drummond et al. (2018a). The resulting temporal and zonal mean of the zonal wind for both these simulations are shown in Figure 10.

As shown in Figure 10, the difference between the primitive and full equations is also apparent when using a radiative transfer scheme, and qualitatively similar to that found in our temperature forced setups (see Figure 1). A full analysis of the radiative transfer simulations is beyond the scope of this paper, and we reserve this for an upcoming work which will focus on the implications for the atmospheric chemistry and observations (requiring radiative transfer) of our findings.

Appendix B Verification of Approximations

As discussed in Section 3.2.3 the traditional approximation, invoked to derive the primitive equations, requires $w\ll v\tan\phi$ to be valid. Also as discussed in Section 3.2.3 this condition can not be met at the equator, and a critical or limiting latitude for the condition can be defined $\phi_{c}$. Where $|\phi|<\phi_{c}$, $w\gtrsim v\tan\phi$, and $w\ll v\tan\phi$ elsewhere. Clearly, the smaller $\phi_{c}$ the more reliable the flow resolved using the primitive equations becomes. Our estimates from Eq.(22) and (17) allow us to understand the parameters controlling the value of $\phi_{c}$: the larger $W$ compared to $V\tan\phi$, the larger $\phi_{c}$ is expected to be. In Figure 11 we verify our estimations by showing the sign of $\frac{v\tan\phi}{10w}-1$ at $p\approx 100Pa$, where yellow and purple denote regions of $v\tan\phi\gg w$ and $v\tan\phi\lesssim w$, respectively. Figure 11 presents values for four of our simulations solving the full equations of motion, namely the Std Full, $R_{\rm p}+$ Full, dT$+$ Full and $\Omega+$ Full simulations (see Table 2 for explanation) after 1 000 days.

Figure 11 shows a broad region extending to $40$ degrees in latitude where $v\tan\phi\gg w$ for the Std Full and dT$+$ Full simulations, where the flow will be resolved differently in the full and primitive cases. However, as expected, when the rotation rate is increased ($\Omega+$ Full) this region is reduced somewhat as the rotational terms become stronger. Finally, the $R_{\rm p}$ simulation results in a significantly reduced region violating the $w\ll v\tan\phi$ consistent with our results, and interpretation.