<figure><img src="image_1.png" /></figure>

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

N. J. Mayne<fnref n="1" />, B. Drummond<fnref n="1" />, F. Debras<fnref n="2" />, E. Jaupart<fnref n="2" />, J. Manners<fnref n="3" />, I. A. Boutle<fnref n="3" />, I. Baraffe<fnref n="1" />,<fnref n="2" />, and K. Kohary<fnref n="1" />

Received 2018 September 10; revised 2018 November 9; accepted 2018 December 5; published 2019 January 22

## 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 daysides and nightsides. 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, and 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.

**Key words:** planets and satellites: atmospheres – planets and satellites: gaseous planets – planets and satellites: terrestrial planets

## 1. Introduction

The detection of exoplanets has revealed several classes of object without direct analogs 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_{\oplus} < R_{\text{p}} < 4R_{\oplus}$, where $R_{\oplus}$ and $R_{\text{p}}$ are the Earth and planetary radii, respectively. In addition, the Next Generation Transit Survey (Wheatley et al. 2013), Characterising ExoPlanet Satellite (Broeg et al. 2013), and *Transiting Exoplanet Survey Satellite* (TESS, Ricker 2014) will significantly increase the number of detections, and availability of characterization 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; Cáceres et al. 2014; Wilson 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.3R_{\oplus}$, 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 (Charnay et al. 2013; Leconte et al. 2013, respectively), alongside studies including surface effects (Lewis et al. 2018) and a dynamical ocean (e.g., Del Genio et al. 2017; Wolf 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

<footnote n="1">¹Astrophysics Group, University of Exeter, Exeter EX4 2QL, UK; n.j.mayne@exeter.ac.uk</footnote>

<footnote n="2">²Ecole normale supérieure de Lyon, CRAL, UMR CNRS 5574, 69364 Lyon Cedex 07, France</footnote>

<footnote n="3">³Met Office, Exeter EX1 3PB, UK</footnote>
included using simplified “relaxation methods” in 3D (Cooper & Showman 2006; Drummond et al. 2018c, 2018b) 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, Helling et al. 2016; Parmentier et al. 2016), and “prognostically” (where the clouds evolve and interact with the atmosphere, Lee et al. 2016; Lines et al. 2018a, 2018b).

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 MITgcm (with the main hot Jupiter adaptations presented in Showman et al. 2009) and the Laboratoire de Méthéorologie Dynamique zoom (LMDz) (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),<fnref n="4" /> 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 π Mensae c (HD 39091c) the latter being recently detected using the 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, summarize 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 UM, which we have adapted to enable modeling 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 minimization scheme (Drummond et al. 2018a) and chemical relaxation scheme (Drummond et al. 2018b, 2018c). We have also recently coupled a chemical kinetics scheme (originally developed in 1D by 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, 2018b), 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_{eq}$, over a parameterized timescale, $\tau_{rad}$. We adopt the $\tau_{rad}$ and $T_{eq}$ used by Zhang & Showman (2017, their Equations (8) and (9), respectively), where $T_{eq}$ is a function of longitude, latitude, and pressure, and $\tau_{rad}$ a function of pressure only. The dayside-nightside temperature contrast is controlled by the parameter $\Delta T_{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_{eq} = 600$ K where $p \le 10$ Pa and decrease linearly with $\ln p$ to $\Delta T_{eq} = 0$ K where $p \ge 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

<footnote n="4">⁴ This code explicitly solves for the centrifugal component, which is usually included in the gravitational term for GCMs (Showman et al. 2008).</footnote>
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, 2014b). The equations are detailed in full in several works, in particular, Mayne et al. (2014b). As described in Mayne et al. (2014b, 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 details 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_p\Omega^2/g$, where $R_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, 2017; Amundsen et al. 2016). 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_p$) and rotation rate ($\Omega$), among 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

<table><caption>**Table 1**
Values of the Standard Parameters for Simulations of GJ 1214b</caption><thead><tr><th>Quantity</th><th>Value</th></tr></thead><tbody><tr><td>Horizontal resolution, $N_\lambda$, $N_\phi$</td><td>144<sub>λ</sub>, 90<sub>φ</sub></td></tr><tr><td>Vertical resolution, $N_z$</td><td>66</td></tr><tr><td>Dynamical timestep, $\Delta t$ (s)</td><td>120</td></tr><tr><td>Initial inner boundary pressure, $p_{max}$ (Pa)</td><td>200 × 10<sup>5</sup></td></tr><tr><td>Rotation rate, $\Omega$ (s<sup>-1</sup>)</td><td>4.60 × 10<sup>-5</sup></td></tr><tr><td>Radius, $R_p$ (m)</td><td>1.45 × 10<sup>7</sup></td></tr><tr><td>Atmospheric height $z_{top}$ (m)</td><td>3.7 × 10<sup>6</sup></td></tr><tr><td>Surface gravity, $g_p$ (m s<sup>-2</sup>)</td><td>12.20</td></tr><tr><td>Specific heat capacity (constant pressure), $c_p$ (J kg<sup>-1</sup> K<sup>-1</sup>)</td><td>12300</td></tr><tr><td>Ideal gas constant, $R$ (J kg<sup>-1</sup> K<sup>-1</sup>)</td><td>3573.5</td></tr><tr><td>Temperature contrast, $\Delta T_{eq}$ (K)</td><td>600</td></tr><tr><td>Diffusion setting, $K_\lambda$ (see Mayne et al. 2017 for details)</td><td>0.158</td></tr><tr><td>Vertical, “sponge,” damping coefficient $R_w$</td><td>0.15</td></tr></tbody></table>

(dT+) we also present a simulation using the deep equation set, giving a total of 15 simulations in this work.

Each simulation is initialized 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 Figure 2. Each simulation was run for 1000 days (throughout this work days refer 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_{eq}$). Since our model is height-based, we adjust the top-of-atmosphere height ($z_{top}$) to achieve a similar pressure range between each simulation. Each simulation encompasses the approximate pressure range ~200 × 10<sup>5</sup> Pa to ~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 (~20%), and have no effect on our conclusions (see test and discussion in Section 3.2.3). The value of $z_{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).
<table><caption>**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 Details)</caption><thead><tr><th>Short Name</th><th>Adjusted Parameters</th><th>z<sub>top</sub> (m)</th><th>Equation Set</th></tr></thead><tbody><tr><td>Std Full</td><td>...</td><td>3.70 × 10<sup>6</sup></td><td>Full</td></tr><tr><td>Std Prim</td><td>...</td><td>3.00 × 10<sup>6</sup></td><td>Primitive</td></tr><tr><td>Std Hires Full</td><td>N<sub>λ</sub> = 288, N<sub>φ</sub> = 180, N<sub>z</sub> = 132, Δ t = 60 s</td><td>3.70 × 10<sup>6</sup></td><td>Full</td></tr><tr><td>Std Hires Prim</td><td>N<sub>λ</sub> = 288, N<sub>φ</sub> = 180, N<sub>z</sub> = 132, Δ t = 60 s</td><td>3.00 × 10<sup>6</sup></td><td>Primitive</td></tr><tr><td>R<sub>p</sub>—Full</td><td>R<sub>p</sub> = 6.0 × 10<sup>6</sup> m, Δ t = 80 s</td><td>5.50 × 10<sup>6</sup></td><td>Full</td></tr><tr><td>R<sub>p</sub>—Prim</td><td>R<sub>p</sub> = 6.0 × 10<sup>6</sup> m, Δ t = 80 s</td><td>3.00 × 10<sup>6</sup></td><td>Primitive</td></tr><tr><td>R<sub>p</sub> + Full</td><td>R<sub>p</sub> = 1.0 × 10<sup>8</sup> m</td><td>3.00 × 10<sup>6</sup></td><td>Full</td></tr><tr><td>R<sub>p</sub> + Prim</td><td>R<sub>p</sub> = 1.0 × 10<sup>8</sup> m</td><td>2.80 × 10<sup>6</sup></td><td>Primitive</td></tr><tr><td>CO<sub>2</sub> Full</td><td>c<sub>p</sub> = 10 123 J kg<sup>-1</sup> K<sup>-1</sup>, R = 188.9 J kg<sup>-1</sup> K<sup>-1</sup>, Δ t = 60</td><td>1.50 × 10<sup>5</sup></td><td>Full</td></tr><tr><td>CO<sub>2</sub> Prim</td><td>c<sub>p</sub> = 10 123 J kg<sup>-1</sup> K<sup>-1</sup>, R = 188.9 J kg<sup>-1</sup> K<sup>-1</sup>, Δ t = 60 s</td><td>1.50 × 10<sup>5</sup></td><td>Primitive</td></tr><tr><td>Ω+ Full</td><td>Ω = 9.2 × 10<sup>-5</sup> s<sup>-1</sup></td><td>3.70 × 10<sup>6</sup></td><td>Full</td></tr><tr><td>Ω+ Prim</td><td>Ω = 9.2 × 10<sup>-5</sup> s<sup>-1</sup></td><td>3.00 × 10<sup>6</sup></td><td>Primitive</td></tr><tr><td>dT+ Full</td><td>ΔT<sub>eq</sub> 800 K</td><td>3.30 × 10<sup>6</sup></td><td>Full</td></tr><tr><td>dT+ Deep</td><td>ΔT<sub>eq</sub> 800 K, g(r) = g<sub>p</sub></td><td>2.65 × 10<sup>6</sup></td><td>Deep</td></tr><tr><td>dT+ Prim</td><td>ΔT<sub>eq</sub> 800 K</td><td>2.65 × 10<sup>6</sup></td><td>Primitive</td></tr></tbody></table>

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 toward 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 nightside will lead to stronger homogenization 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, normalized 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 Equation (18) of Zhang & Showman 2017). The phase curve for a given pressure level is then normalized by (Zhang & Showman 2017, their Equation (18)),

$$ \bar{F}(\delta) = \frac{F(\delta) - \langle F(\delta) \rangle}{\langle F(\delta) \rangle}, \qquad (1) $$

where $\bar{F}(\delta)$ is the normalized phase curve, $F(\delta)$ is the emission flux as a function of phase ($\delta$), and $\langle F(\delta) \rangle$ is the mean of the flux over the entire phase. We note that the normalized phase curve is negative where $F(\delta) < \langle F(\delta) \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

Figure 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 color scales are varied between simulation “pairs” i.e., each matching pair of primitive and full simulations have the same color 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 ~1.5× 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 (bottom 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 (color scale) and horizontal wind (vector arrows) at isobaric surfaces (in $\lambda$ and $\phi$) of 100 and 3000 Pa at 1000 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,
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Figure showing the zonal and temporal mean of the zonal wind (m s⁻¹) as a function of latitude (φ) and pressure (log₁₀(*p* [Pa]), for the primitive and full versions of the Std and Std Hires simulations. Note the modeled pressure domain extends down to 200 × 10⁵ Pa, but only the relatively, dynamically active region of the atmosphere is shown here. The simulation short name (see Table 2 for explanation of simulation names) and temporal averaging period are given below each subfigure.</figcaption></figure>

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 substellar 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 3000 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 nightside.

The amplitude of the simplified phase curve is proportional to T⁴ cos φ, heavily weighted to the equatorial regions.<fnref n="5" /> 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 substellar 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

<footnote n="5">⁵ We note that we show normalized phase curves (see Equation (1)) which means that relative differences in the amplitude of the normalized phase curve between simulations can be larger than might be expected based on the relative temperature differences.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Figure showing the temperature (color scale, K) and horizontal wind (vector arrows) as a function of $\lambda$ and $\phi$, for isobaric surfaces at 100 and 3000 Pa at 1000 days for the Std Full simulation and the difference with the Std Prim (in the sense Std Full–Std Prim), as well as the simplified thermal phase curves at both pressures for both simulations (see Table 2 for explanation of simulation names). Note the change in the vertical axes for the bottom panels.</figcaption></figure>
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

$$ \frac{Du}{Dt} = \frac{uv \tan \phi}{r} - \frac{uw}{r} + fv - f'w - \frac{c_p \theta}{r \cos \phi} \frac{\partial \Pi}{\partial \lambda} + \mathbf{D}(u), \quad (2) $$

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

$$ f = 2\Omega \sin \phi, \quad (3) $$

and

$$ f' = 2\Omega \cos \phi. \quad (4) $$

The terms in Equation (2) can be “grouped” into those associated with rotation, i.e., $f$ and $f'$, 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_p \Omega^2/g$ is small (Bánard 2014). For our simulations, Tables 1 and 2 show $R_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_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 Equation (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,

$$ \frac{Du}{Dt} = \frac{uv \tan \phi}{R_p} + fv - \frac{c_p \theta}{R_p \cos \phi} \frac{\partial \Pi}{\partial \lambda} + \mathbf{D}(u), \quad (5) $$

with the full equation (Equation (2)). Specifically, for Equations (2) and (5), the shallow-fluid approximation results in $r$ being replaced with $R_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'$ (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 scale height, $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_p$, where $r$ can be expressed as $R_p + z_U$ we require $z_U \ll R_p$. For a super Earth $R_p$ is fixed as it has a solid surface, but for a small Neptune $R_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 modeled 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_p^+$ and $R_p^-$, respectively, see Table 2), to explore this effect. Additionally, we perform a simulation adopting the standard planetary radius, but with a CO$_2$-dominated atmosphere, following Zhang & Showman (2017), leading to an increased mean molecular weight, and subsequent vertical compression of the atmosphere itself.

Figure 3 shows the same information, in the same format as Figure 1 but for the $R_p^-$, $R_p^+$, and CO$_2$ simulations. Clearly, for the $R_p^-$ simulation differences remain between the primitive and full simulations, whereas these differences are almost completely removed when the planetary radius is
<figure><img src="image_4.png" /><figcaption>(a) $R_p$- Prim: 800-1 000 days</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>(b) $R_p$- Full: 800-1 000 days</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>(c) $R_p$+ Prim: 800-1 000 days</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>(d) $R_p$+ Full: 800-1 000 days</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(e) CO₂ Prim: 800-1 000 days</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(f) CO₂ Full: 800-1 000 days</figcaption></figure>

**Figure 3.** Figure similar to Figure 1 but for simulations adopting a decreased ($R_p$−) and increased ($R_p$+) planetary radius, and a CO₂ dominated atmosphere leading to a significantly reduced atmospheric scale height.
increased ($R_p$+). They are also significantly reduced when the atmospheric scale height is reduced ($CO_2$) 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_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 scale height 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_p$+ and $CO_2$ 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 difference is, at its maximum, ~20% (compared to a maximum of 10% in the $R_p$+ simulation).

For the $R_p$+ the temperature changes between the primitive and full simulations, at both 100 and 3000 Pa, are generally less than ~2 K. Therefore, figures of the temperature structure for these simulations are omitted. In the case of the $CO_2$ 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 toward lower pressures. Figures 4 and 5 shows the same information as Figure 2 but for the $R_p$- and $CO_2$ simulations.

The pattern of the change, at 100 Pa, when moving to the more complete equations of motion for the $R_p$- case is similar to our standard GJ 1214b setups, but enhanced. The regions east of the substellar point are heated and those to the west are cooled, with a net dayside heating and nightside cooling. For the higher pressure, 3000 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 dayside amplitude, and reduction in the nightside 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 $CO_2$ 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 superrotating jet at the equator, although some more complex changes occur at high latitudes for the lower pressure surface of the $CO_2$ 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_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 do not 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 behavior of the atmosphere and understand the simulation results. First, 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$ to $10^5$ Pa (or 1 mbar to 1 bar), it can reasonably be assumed that:

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

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

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

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

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

$$ \frac{D\rho}{Dt} = 0, \quad (7) $$

$$ \frac{DT}{Dt} = \frac{T_{eq} - T}{\tau_{\text{rad}}}, \quad (8) $$

and

$$ p = \rho RT, \quad (9) $$

respectively, where $\mathbf{v}$, $\nabla$, $\Omega$ and $\mathbf{g}$ are merely the vector forms of the variables and operator previously defined. Using the vector identity $\nabla(v^2/2) = (\mathbf{v} \cdot \nabla)\mathbf{v} + \mathbf{v} \wedge (\nabla \wedge \mathbf{v})$ and
<figure><img src="image_10.png" /><figcaption>(a) $R_p$ - Full: 100 Pa</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(b) $R_p$ - Full: 3000 Pa</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(c) $R_p$ - Full-Prim: 100 Pa</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(d) $R_p$ - Full-Prim: 3000 Pa</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(e) $R_p$ - Prim/Full: 100 Pa</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(f) $R_p$ - Prim/Full: 3000 Pa</figcaption></figure>

**Figure 4.** Figure similar to Figure 2, but for the reduced planetary radius simulations, $R_p$ - (see Table 2 for an explanation of the simulation names). Note the change in the vertical axes for the bottom panels.
<figure><img src="image_16.png" /><figcaption>(a) CO₂ Full: 100 Pa</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(b) CO₂ Full: 3 000 Pa</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(c) CO₂ Full-Prim: 100 Pa</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(d) CO₂ Full-Prim: 3 000 Pa</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(e) CO₂ Prim/Full: 100 Pa</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(f) CO₂ Prim/Full: 3 000 Pa</figcaption></figure>

**Figure 5.** Figure similar to Figure 2, but for the CO₂ simulations (see Table 2 for an explanation of the simulation names). Note the change in the vertical axes for the bottom panels.
taking the scalar product of Equation (6) with $v$ we obtain

$$v \cdot \nabla \left( \frac{v^2}{2} \right) = -\frac{1}{\rho} \frac{Dp}{Dt} - wg. \quad (10)$$

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 $\nabla \cdot v = 0$ and therefore $\frac{1}{r} \frac{\partial u}{\partial \lambda} + \frac{1}{r} \frac{\partial v}{\partial \phi} + \frac{\partial w}{\partial r} = 0$ suggesting the standard order of magnitude estimate of

$$\frac{U}{L_U} \sim \frac{W}{H_W}, \quad (11)$$

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$,

$$-\frac{1}{\rho} \frac{Dp}{Dt} - wg = -\frac{1}{\rho} (v \cdot \nabla_{\perp}) p \sim \frac{1}{\rho} \frac{Dp}{Dt}, \quad (12)$$

where $\nabla_{\perp}$ 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 Equation (9) using (8) we obtain

$$\frac{Dp}{Dt} = \rho R \frac{DT}{Dt} = \rho R \frac{(T_{eq} - T)}{\tau_{rad}}. \quad (13)$$

Combining Equations (10), (12), and (13) yields

$$v \cdot \nabla \left( \frac{v^2}{2} \right) \sim R \frac{(T_{eq} - T)}{\tau_{rad}}, \quad (14)$$

which can be simplified to

$$\frac{U^3}{2L_U} \sim R \frac{\Delta T}{\tau_{rad}}, \quad (15)$$

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_{eq,day} + T_{eq,night})/2$, where $T_{eq,day}$ and $T_{eq,night}$ are the average day and nightside equilibrium temperatures. Therefore, on average the difference between the steady-state temperature and equilibrium temperature is $\Delta T \sim (T_{eq,day} - T_{eq,night})/2$.

For our simulations the characteristic length scale for the zonal jet is about half the planetary circumference or $\sim \pi R_p$ and we expect $\Delta T \sim 0.5 \Delta T_{eq}$ (or half the total day-night contrast, see Section 2), which is 300 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 Equation (15) yields

$$U \sim \sqrt[3]{\frac{2\pi R_p R \Delta T}{\tau_{rad}}} \sim 1.4 \times 10^3 \text{ m s}^{-1}. \quad (16)$$

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. Equation (16) can then be substituted into Equation (11) to provide an estimate for $W$,

$$W \sim \frac{H}{L} \sqrt{2\pi R_p R \frac{\Delta T}{\tau_{rad}}}, \quad (17)$$

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 Equation (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, Equation (6), assuming that the Coriolis term is dominant can be written as

$$\frac{Dv}{Dt} \sim 2\Omega U \sin \phi, \quad (18)$$

which can be approximated by

$$\frac{UV}{L_V} \sim 2\Omega U \sin \phi, \quad (19)$$

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

$$\frac{UV}{L_V} \sim 2\Omega U \sin \phi V \sim 2L_V \Omega \sin \phi. \quad (20)$$

In our simulations which $L_V \sim \pi/4R_p$, yielding

$$V \sim \frac{\pi}{2} R_p \Omega \sin \phi \sim 400 \sin \phi \text{ ms}^{-1}. \quad (21)$$

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$,

$$V \tan \phi \sim \frac{\pi}{2} R_p \Omega \frac{\sin \phi^2}{\cos \phi}, \quad (22)$$

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

The equations we have derived (Equations (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. Equations (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. Equation (22) then requires the further assumption that the $V$ component is driven, chiefly, by the Coriolis forces, where Equation (17) requires the assumption that the flow is divergence-free. The derivation of these equations, however, does not rely on adopting either the
<figure><img src="image_22.png" /><figcaption>**Figure 6.** Figure similar to Figure 1, but showing the relative lack of difference between the primitive and full equations for a faster rotating planet (see Table 2 for an explanation of the simulation names).</figcaption></figure>

primitive or full equations explicitly. In Appendix B we explicitly verify these equations using our simulation outputs.

Although very simplified, without a proper prescription of *H*, *L*, or $\Delta T$, Equations (17) and (22) allow us to draw insight into the nature of the flow. First, 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, Equations (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 Equation (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_p$ remains constant, for the terrestrial regime. If $\Omega R_p$ is constant, Equation (22) is clearly constant for a given $\phi$. Equation (17) can be expressed as $W \propto \frac{R_p^{1/3}}{L}$, which using $L \sim R_p$ becomes $W \propto \frac{1}{R_p^{2/3}}$, which at constant $\Omega R_p$ can be expressed as $W \propto \Omega^{2/3}$. Therefore, at constant $\Omega R_p$ as $\Omega$ is increased $W$ increases, while $V$ is constant, meaning the condition $w \ll v \tan \phi$ becomes increasingly less valid.

Equations (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. (2018b).

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 CO₂ 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 CO₂ simulation the ratio of scale height of the atmosphere, to the planetary radius is similar to that of the $R_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 CO₂ case, which are not observed to the same extent in the $R_p$+ simulation, are likely caused by the change in the mean molecular weight itself, driving a change in the $W$ field, as shown in Equation (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), Equation (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 a 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
<figure><img src="image_23.png" /><figcaption>(a) Ω+ Full: 100 Pa</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(b) Ω+ Full: 3000 Pa</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>(c) Ω+ Full-Prim: 100 Pa</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(d) Ω+ Full-Prim: 3000 Pa</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(e) Ω+ Prim/Full: 100 Pa</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(f) Ω+ Prim/Full: 3000 Pa</figcaption></figure>

**Figure 7.** Figure similar to Figure 2, but for the Ω+ simulations (see Table 2 for an explanation of the simulation names). Note the change in the vertical axes for the bottom panels.
<figure><img src="image_29.png" /><figcaption>(a) dT+ Deep: 800-1 000 days</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(b) dT+ Prim: 800-1 000 days</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(c) dT+ Full: 800-1 000 days</figcaption></figure>

**Figure 8.** Figure similar to Figure 1, but for the case with an increased day-night temperature contrast, dT+ simulation, but including the version solving the deep equations (i.e., full but with constant gravity, see Table 2 for an explanation of the simulation names).

of the simplified phase curve is found, as temperature changes at
this level are around a few degrees.

*Forcing/advection and gravity*—Equation (17) shows that as the temperature contrast is increased *W* increases, leading to a weakening of the *w* ≪ *v* tan *φ* condition. Additionally, the acceleration due to gravity does not appear in our estimate of *V* tan *φ*, but does enter the estimate for *W* via the scale height (*H* ∝ *RT*/*g*) (Equations (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 an 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 ~30%, with the difference rapidly reducing
with altitude, ∝1/*r*². Additionally, the *z*top, and therefore, the
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
<figure><img src="image_32.png" /><figcaption>**Figure 9.** Figure similar to Figure 2, but for the dT+ simulations (see Table 2 for an explanation of the simulation names). The phase curves include the additional deep simulation. Note the change in the vertical axes for the bottom panels.</figcaption></figure>
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.

Figure 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_p$ becomes a 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_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 modeled 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_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_{\text{atm}}$, is

negligible compared to the mass of the bulk planet comprising the unmodeled material within the inner boundary, $M_{p}$. For our basic setup the adopted surface gravity, $g_{p} = 12.20 \text{ m s}^{-2}$ and planetary radius, $R_{p} = 1.45 \times 10^{7}$ m, imply $M_{p} \sim 3.8 \times 10^{25}$ kg and the standard simulations indicate $M_{\text{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_{\text{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_{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 modeled 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 intercomparison 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 Gyr, 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_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_s \sim \sqrt{\gamma RT}$, with $\gamma = 1.5$ and temperatures in the jet region of $\sim 500$ K this yields $c_s \sim 1600 \text{ m s}^{-1}$ for our simulations ($\sim 400 \text{ m s}^{-1}$ for the CO$_2$ 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 that 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, and Std Hires Prim), which are set up 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_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. 2018b, 2018c), and include the cloud scheme adopted in Lines et al. (2018b). 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.

N.J.M. is partly funded by a Leverhulme Trust Research Project Grant and gratefully acknowledges their support. This work is partly funded by the ERC grant Nos. 336792 and 787361-COBOM. J.M. and I.A.B. acknowledge the support of a Met Office Academic Partnership secondment. This work was partly supported by a Science and Technology Facilities Council Consolidated Grant (ST/R000395/1). Material was 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. The research data supporting this publication are openly available from the University of Exeter's institutional repository at: https://doi.org/10.24378/exe.1023, using http://hdl.handle.net/10871/35127.

## Appendix A
### Radiative Transfer

All of our simulations analyzed 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 analyzing 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 set up to match the Std Prim and Std Full simulations, with the complete setup (including radiative transfer elements) detailed in
<figure><img src="image_33.png" /><figcaption>**Figure 10.** Figure similar to Figure 1 for simulations matching the Std Prim and Std Full, but including full radiative transfer (as described in Drummond et al. 2018a) as opposed to temperature forcing (see Table 2 for an explanation of the simulation names). Note the shorter simulation time and averaging period.</figcaption></figure>

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 cannot 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 Equations (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 100$ Pa, 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_p$+ Full, dT+ Full, and $\Omega$+ Full simulations (see Table 2 for an explanation) after 1000 days.

Figure 11 shows a broad region extending to 40° 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_p$ simulation results in a significantly reduced region violating the $w \ll v \tan \phi$ consistent with our results, and interpretation.
<figure><img src="image_34.png" /><figcaption>(a) Std Full</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>(b) R<sub>p</sub> + Full</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(c) dT + Full</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(d) Ω + Full</figcaption></figure>

**Figure 11.** Figure showing the sign of $\frac{v \tan \phi}{10w} - 1$ as a function of longitude and latitude at a height corresponding to an equatorial pressure of ~100 Pa for four simulations, after 1000 days. Yellow and purple regions show positive ($w \ll v \tan \phi$) and negative ($w \gtrsim v \tan \phi$) values.

ORCID iDs

N. J. Mayne https://orcid.org/0000-0001-6707-4563

B. Drummond https://orcid.org/0000-0001-7589-5484

F. Debras https://orcid.org/0000-0002-3113-4840

References

Agúndez, M., Parmentier, V., Venot, O., Hersant, F., & Selsis, F. 2014, A&A, 564, A77

Agúndez, M., Venot, O., Iro, N., et al. 2012, A&A, 548, A73

Amundsen, D. S., Baraffe, I., Tremblin, P., et al. 2014, A&A, 564, A59

Amundsen, D. S., Mayne, N. J., Baraffe, I., et al. 2016, A&A, 595, A36

Amundsen, D. S., Tremblin, P., Manners, J., Baraffe, I., & Mayne, N. J. 2017, A&A, 598, A97

Bánard, P. 2014, QJRMS, 141, 195

Bean, J. L., Désert, J.-M., Kabath, P., et al. 2011, ApJ, 743, 92

Bean, J. L., Miller-Ricci Kempton, E., & Homeier, D. 2010, Natur, 468, 669

Berta, Z. K., Charbonneau, D., Désert, J.-M., et al. 2012, ApJ, 747, 35

Boutle, I. A., Mayne, N. J., Drummond, B., et al. 2017, A&A, 601, A120

Broeg, C., Fortier, A., Ehrenreich, D., et al. 2013, European Physical Journal Web of Conferences, 47, 03005

Burrows, A., & Sharp, C. M. 1999, ApJ, 512, 843

Cáceres, C., Kabath, P., Hoyer, S., et al. 2014, A&A, 565, A7

Carter, J. A., Winn, J. N., Holman, M. J., et al. 2011, ApJ, 730, 82

Charnay, B., Forget, F., Wordsworth, R., et al. 2013, JGRD, 118, 10

Charnay, B., Meadows, V., & Leconte, J. 2015a, ApJ, 813, 15

Charnay, B., Meadows, V., Misra, A., Leconte, J., & Arney, G. 2015b, ApJL, 813, L1

Cooper, C. S., & Showman, A. P. 2006, ApJ, 649, 1048

Crossfield, I. J. M., Barman, T., & Hansen, B. M. S. 2011, ApJ, 736, 132

Del Genio, A. D., Way, M. J., Amundsen, D. S., et al. 2017, arXiv:1709.02051

de Mooij, E. J. W., Brogi, M., de Kok, R. J., et al. 2012, A&A, 538, A46

Deming, D., Harrington, J., Laughlin, G., et al. 2007, ApJL, 667, L199

Désert, J.-M., Bean, J., Miller-Ricci Kempton, E., et al. 2011, ApJL, 731, L40

Dobbs-Dixon, I., & Agol, E. 2013, MNRAS, 435, 3159

Drummond, B., Mayne, N. J., Baraffe, I., et al. 2018a, A&A, 612, A105

Drummond, B., Mayne, N. J., Manners, J., et al. 2018b, arXiv:1810.09724

Drummond, B., Mayne, N. J., Manners, J., et al. 2018c, ApJL, 855, L31

Drummond, B., Tremblin, P., Baraffe, I., et al. 2016, A&A, 594, A69

Elkins-Tanton, L. T., & Seager, S. 2008, ApJ, 685, 1237

Fraine, J. D., Deming, D., Gillon, M., et al. 2013, ApJ, 765, 127

Fressin, F., Torres, G., Charbonneau, D., et al. 2013, ApJ, 766, 81

Fromang, S., Leconte, J., & Heng, K. 2016, A&A, 591, A144

Gerkema, T., Zimmerman, J. T. F., Maas, L. R. M., & van Haren, H. 2008, RvGeo, 46, RG2004

Helling, C., Lee, G., Dobbs-Dixon, I., et al. 2016, MNRAS, 460, 855

Huang, C. X., Burt, J., Vanderburg, A., et al. 2018, arXiv:1809.05967

Hubbard, W. B., Fortney, J. J., Lunine, J. I., et al. 2001, ApJ, 560, 413

Iro, N., Bézard, B., & Guillot, T. 2005, A&A, 436, 719

Kataria, T., Showman, A. P., Fortney, J. J., Marley, M. S., & Freedman, R. S. 2014, ApJ, 785, 92

Kataria, T., Sing, D. K., Lewis, N. K., et al. 2016, ApJ, 821, 9

Knutson, H. A., Benneke, B., Deming, D., & Homeier, D. 2014, Natur, 505, 66
Knutson, H. A., Lewis, N., Fortney, J. J., et al. 2012, ApJ, 754, 22

Kreidberg, L., Bean, J. L., Désert, J.-M., et al. 2014, Natur, 505, 69

Leconte, J., Forget, F., Charnay, B., Wordsworth, R., & Pottier, A. 2013, Natur, 504, 268

Lee, G., Dobbs-Dixon, I., Helling, C., Bognar, K., & Woitke, P. 2016, A&A, 594, A48

Lewis, N. T., Lambert, F. H., Boutle, I. A., et al. 2018, ApJ, 854, 171

Lines, S., Manners, J., Mayne, N. J., et al. 2018a, MNRAS, 481, 194

Lines, S., Mayne, N. J., Boutle, I. A., et al. 2018b, A&A, 615, A97

Lopez, E. D., & Fortney, J. J. 2014, ApJ, 792, 1

Madhusudhan, N., Agúndez, M., Moses, J. I., & Hu, Y. 2016, SSRv, 205, 285

Mayne, N. J., Baraffe, I., Acreman, D. M., et al. 2014a, GMD, 7, 3059

Mayne, N. J., Baraffe, I., Acreman, D. M., et al. 2014b, A&A, 561, A1

Mayne, N. J., Debras, F., Baraffe, I., et al. 2017, A&A, 604, A79

Mendonça, J. M., Grimm, S. L., Grosheintz, L., & Heng, K. 2016, ApJ, 829, 115

Menou, K. 2012, ApJL, 744, L16

Owen, J. E., & Jackson, A. P. 2012, MNRAS, 425, 2931

Parmentier, V., Fortney, J. J., Showman, A. P., Morley, C., & Marley, M. S. 2016, ApJ, 828, 22

Pollack, J. B., Hubickyj, O., Bodenheimer, P., et al. 1996, Icar, 124, 62

Ricker, G. R. 2014, JAVSO, 42, 234

Showman, A. P., Fortney, J. J., Lian, Y., et al. 2009, ApJ, 699, 564

Showman, A. P., Menou, K., & Cho, J. Y.-K. 2008, in ASP Conf. Ser. 398, Extreme Solar Systems, ed. D. Fischer et al. (San Francisco, CA: ASP), 419

Showman, A. P., & Polvani, L. M. 2011, ApJ, 738, 71

Sing, D. K., Désert, J.-M., Fortney, J. J., et al. 2011, A&A, 527, A73

Todorov, K. O., Deming, D., Burrows, A., & Grillmair, C. J. 2014, ApJ, 796, 100

Tokano, T. 2013, GeoRL, 40, 4538

Tort, M., Dubos, T., & Melvin, T. 2015, QJRMS, 141, 3056

Tremblin, P., Amundsen, D. S., Mourier, P., et al. 2015, ApJL, 804, L17

Tremblin, P., Chabrier, G., Mayne, N. J., et al. 2017, ApJ, 841, 30

Tsai, S.-M., Dobbs-Dixon, I., & Gu, P.-G. 2014, ApJ, 793, 141

Turbet, M., Leconte, J., Selsis, F., et al. 2016, A&A, 596, A112

Vallis, G. K. 2006, Atmospheric and Ocean Fluid Dynamics (Cambridge: Cambridge Univ. Press)

Venot, O., Drummond, B., Miguel, Y., et al. 2017, arXiv:1711.08433

Wheatley, P. J., Pollacco, D. L., Queloz, D., et al. 2013, European Physical Journal Web of Conferences, 47, 13002

White, A. A., & Bromley, R. A. 1995, QJRMS, 121, 399

Wilson, P. A., Colón, K. D., Sing, D. K., et al. 2014, MNRAS, 438, 2395

Wolf, E. T., Shields, A. L., Kopparapu, R. K., Haqq-Misra, J., & Toon, O. B. 2017, ApJ, 837, 107

Zhang, X., & Showman, A. P. 2017, ApJ, 836, 73