Results from a set of three-dimensional numerical experiments of a hot Jupiter atmosphere

TL;DR

This study uses three-dimensional numerical simulations to investigate the dynamics of a hot Jupiter atmosphere, revealing the robustness of super-rotating equatorial jets and complex mechanisms behind jet acceleration.

Contribution

It provides new insights into the stability and driving mechanisms of atmospheric jets on hot Jupiters through extensive 3D modeling.

Findings

Super-rotating equatorial jet is robust across parameter changes.

Deep atmosphere forcing can diminish the jet.

Jet acceleration involves complex eddy and mean flow interactions.

Abstract

We present highlights from a large set of simulations of a hot Jupiter atmosphere, nominally based on HD 209458b, aimed at exploring both the evolution of the deep atmosphere, and the acceleration of the zonal flow or jet. We find the occurrence of a super-rotating equatorial jet is robust to changes in various parameters, and over long timescales, even in the absence of strong inner or bottom boundary drag. This jet is diminished in one simulation only, where we strongly force the deep atmosphere equator-to-pole temperature gradient over long timescales. Finally, although the eddy momentum fluxes in our atmosphere show similarities with the proposed mechanism for accelerating jets on tidally-locked planets, the picture appears more complex. We present tentative evidence for a jet driven by a combination of eddy momentum transport and mean flow.

Tables

Table 1. Table 1: Value of the standard parameters for the temperature forced (TF) and radiative transfer (RT) simulations.

Quantity	TF	RT
Horizontal resolution	${144}_{\lambda }$, ${90}_{\varphi }$
Vertical resolution, $N_z$	66
Dynamical Timestep (s)	1200	30
Radiative Timestep (s)	-	150
Initial inner boundary pressure, $p_{\max }$ (Pascals, Pa)	$220\times {10}^5$	$200\times {10}^5$
Rotation rate, $\Omega$ (s-1)	$2.06\times {10}^{-5}$
Radius, $R_p$ (m)	$9.44\times {10}^7$	$9.0\times {10}^7$
Radius to outer boundary, $R_{\mathrm{top}}$ (m)	$1.1\times {10}^7$	$9.0\times {10}^6$
Surface gravity, $g_p$ (ms-2)	9.42	10.79
Specific heat capacity (constant pressure), $c_p$ (Jkg-1K-1)	14 308.4	13 000.0
Ideal gas constant, $R$ (Jkg-1K-1)(1)	4593	3556.8
Diffusion setting, $K_{\lambda }$ (see Section 2.1 for applied coefficient, $K_{\lambda }^{\prime }$)	0.158
Horizontal vector laplacian damping coefficient, ${\nu }_{\lambda ,\varphi }$	0.1	0.0
Vertical vector laplacian damping coefficient, ${\nu }_r$	0.0
Vertical, ‘sponge’, damping coefficient $R_w$	0.15
Temperature characterising intrinsic luminosity of the planet (K), $T_{\mathrm{int}}$	-	100

Table 2. Table 2: Details, alongside short names, of the key simulations used in this work (a more extensive set of simulations were performed, but we only highlight those where results are directly used in this work). The heating/cooling scheme (temperature forced, TF, or radiative transfer, RT), total elapsed simulation time (rounded down to the nearest 100 days), variables or settings which are different from the standard models (see Table 1 for the basic parameters), and the sophistication of the dynamical equations is shown (see Mayne et al., 2014a , for definition) . The time sampling in all cases was ten days, with an additional run outputting at 1 day performed only over the first 100 days for each simulation. p 𝑝 p is the pressure, τ rad subscript 𝜏 rad \tau_{\rm rad} the radiative timescale and T eq subscript 𝑇 eq T_{\rm eq} the equilibrium temperature.

Short Name	Heating/Cooling	Length (Earth days)	Adjusted parameters		Equation set
Std Prim	TF	10 200	$R_w=0.20$		“Primitive”
Std Full	TF	13 300	-		“Full”
Std RT	RT	1 600	-		“Full”
Reduced $p_{\max }$	TF	15 400	$p_{\max }={10}^6$ Pa, $R_w=0.20$		“Full”
Deep $\Delta T_{\mathrm{eq}\to \mathrm{pole}}$	TF	10 800	$p>{10}^6\mathrm{Pa}$:	${\tau }_{\mathrm{rad}}\left(p\right)={\tau }_{\mathrm{rad}}\left({10}^6\mathrm{Pa}\right)$	“Full”
				$T_{\mathrm{eq}}=T_{\mathrm{eq}}+\mathrm{1\hspace{0.17em}000.0}{\sin }^2\left(\varphi \right)$

Equations

$$K_{\lambda }^{\prime }=\frac{\Delta t{K}_{\lambda }}{(r^2{\cos }^2\varphi (\Delta \lambda {)}^2)}.$$

$${\tau }_{\mathrm{th}}\sim \frac{C_{\mathrm{p}}T\Delta m}{L},$$

$$(\overline{\rho }\overline{u}{)}_{,t}$$

$$+2\Omega \overline{\rho v}\sin \varphi -2\Omega \overline{\rho w}\cos \varphi$$

$$-(\overline{{\rho }^{\prime }{u}^{\prime }}{)}_{,t}-\frac{{\left[\overline{(\rho v{)}^{\prime }{u}^{\prime }}{\cos }^2\varphi \right]}_{,\varphi }}{r{\cos }^2\varphi }-\frac{{\left[\overline{(\rho w{)}^{\prime }{u}^{\prime }}{r}^3\right]}_{,r}}{r^3}+$$

$$\overline{\rho G_{\lambda }},$$

$$\frac{Du}{Dt}$$

$$-\frac{uw}{r}+\frac{uv}{r}\tan \varphi +2\Omega v\sin \varphi -2\Omega w\cos \varphi$$

$$-\frac{1}{r\cos \varphi }\frac{1}{\rho }\frac{\partial p}{\partial \lambda }+G_{\lambda }$$

$$\frac{Du}{Dt}$$

$$\frac{\partial \rho }{\partial t}+\nabla .\left(\rho \text{u}\right)=0.$$

$$\rho \frac{Du}{Dt}$$

$$\frac{\partial \left(\rho u\right)}{\partial t}$$

$$\nabla .\left(\rho u\text{u}\right)$$

$$\rho \frac{Du}{Dt}$$

$$\nabla .\left(X\text{Y}\right)=$$

$$\frac{1}{r\cos \varphi }\frac{\partial \left(X{Y}_{\lambda }\right)}{\partial \lambda }+\frac{1}{r\cos \varphi }\frac{\partial \left(X{Y}_{\varphi }\cos \varphi \right)}{\partial \varphi }+\frac{1}{r^2}\frac{\partial \left(r^{2}X{Y}_r\right)}{\partial r},$$

$$\rho \frac{Du}{Dt}$$

$$\frac{\partial \left(\rho u\right)}{\partial t}+\frac{1}{r\cos \varphi }\frac{\partial \left(\rho u^2\right)}{\partial \lambda }$$

$$+\frac{1}{r\cos \varphi }\frac{\partial \left(\rho uv\cos \varphi \right)}{\partial \varphi }+\frac{1}{r^2}\frac{\partial \left(r^2\rho uw\right)}{\partial r}$$

$$\rho \frac{Du}{Dt}$$

$$-\frac{\rho uw}{r}+\frac{\rho uv}{r}\tan \varphi +2\Omega \rho v\sin \varphi -2\Omega \rho w\cos \varphi$$

$$-\frac{1}{r\cos \varphi }\frac{\partial p}{\partial \lambda }+\rho G_{\lambda },$$

$$\frac{\partial \left(\rho u\right)}{\partial t}+\frac{1}{r\cos \varphi }\frac{\partial \left(\rho u^2\right)}{\partial \lambda }+$$

$$\frac{1}{r\cos \varphi }\frac{\partial \left(\rho uv\cos \varphi \right)}{\partial \varphi }+\frac{1}{r^2}\frac{\partial \left(r^2\rho uw\right)}{\partial r}=$$

$$-\frac{\rho uw}{r}+\frac{\rho uv}{r}\tan \varphi +2\Omega \rho v\sin \varphi -2\Omega \rho w\cos \varphi$$

$$-\frac{1}{r\cos \varphi }\frac{\partial p}{\partial \lambda }+\rho G_{\lambda }.$$

$$\frac{\partial \left(\rho u\right)}{\partial t}+\frac{1}{r\cos \varphi }\frac{\partial \left(\rho u^2\right)}{\partial \lambda }+\frac{1}{r\cos \varphi }\frac{\partial p}{\partial \lambda }-\rho G_{\lambda }=$$

$$-\frac{1}{r\cos \varphi }\frac{\partial \left(\rho uv\cos \varphi \right)}{\partial \varphi }+\frac{\rho uv}{r}\tan \varphi +2\Omega \rho v\sin \varphi$$

$$-\frac{1}{r^2}\frac{\partial \left(r^2\rho uw\right)}{\partial r}-\frac{\rho uw}{r}-2\Omega \rho w\cos \varphi .$$

$$-\frac{1}{r\cos \varphi }\frac{\partial \left(\rho uv\cos \varphi \right)}{\partial \varphi }+\frac{\rho uv}{r}\tan \varphi =$$

$$-\frac{1}{r\cos \varphi }\frac{\partial \left(\rho uv\cos \varphi \right)}{\partial \varphi }+\frac{\rho uv\sin \varphi }{r\cos \varphi },$$

$$-\frac{1}{r\cos \varphi }\frac{\partial \left(\rho uv\cos \varphi \right)}{\partial \varphi }+\frac{\rho uv\sin \varphi }{r\cos \varphi }=$$

$$-\frac{1}{r\cos \varphi }\frac{\partial \left(\rho uv\cos \varphi \right)}{\partial \varphi }-\frac{\rho uv}{r\cos \varphi }\frac{\partial \left(\cos \varphi \right)}{\partial \varphi }=$$

$$-\frac{\cos \varphi }{r{\cos }^2\varphi }\frac{\partial \left(\rho uv\cos \varphi \right)}{\partial \varphi }-\frac{\rho uv\cos \varphi }{r{\cos }^2\varphi }\frac{\partial \left(\cos \varphi \right)}{\partial \varphi }.$$

$$-\frac{\cos \varphi }{r{\cos }^2\varphi }\frac{\partial \left(\rho uv\cos \varphi \right)}{\partial \varphi }-\frac{\rho uv\cos \varphi }{r{\cos }^2\varphi }\frac{\partial \left(\cos \varphi \right)}{\partial \varphi }=$$

$$-\frac{1}{r{\cos }^2\varphi }\frac{\partial \left(\rho uv{\cos }^2\varphi \right)}{\partial \varphi }.$$

$$-\frac{1}{r^2}\frac{\partial \left(r^2\rho uw\right)}{\partial r}-\frac{\rho uw}{r}=$$

$$-\frac{r}{r{r}^2}\frac{\partial \left(r^2\rho uw\right)}{\partial r}-\frac{r^2\rho uw}{r^{2}r}\frac{\partial r}{\partial r},$$

$$-\frac{r}{r{r}^2}\frac{\partial \left(r^2\rho uw\right)}{\partial r}-\frac{r^2\rho uw}{r^{2}r}\frac{\partial r}{\partial r}=$$

$$-\frac{1}{r^3}\frac{\partial \left(\rho uw{r}^3\right)}{\partial r}.$$

$$\frac{\partial \left(\rho u\right)}{\partial t}+\frac{1}{r\cos \varphi }\frac{\partial \left(\rho u^2\right)}{\partial \lambda }+\frac{1}{r\cos \varphi }\frac{\partial p}{\partial \lambda }-\rho G_{\lambda }=$$

Full text

1 Introduction

Three-dimensional (3D) general circulation models (GCMs) have become established tools used both to interpret current observations of exoplanets, and to predict those to be made by future instruments (see Kataria et al., 2016; Lee et al., 2016, for recent examples). Given the complex, non–linear and interacting physical mechanisms acting within a planetary atmosphere, combined with the incomplete nature of our observational access to exoplanets, GCMs are likely to become increasingly important in this field.

Currently, the most observationally constrained exoplanets are a subset termed hot Jupiters, being Jovian in size and orbiting close to their parent star (see discussion and review in Baraffe et al., 2010). Radial velocity and transit measurements provide estimates of the masses and radii of these planets, and numerous subsequent studies have inferred further atmospheric characteristics such as chemical compositions, temperature structures and even wind speeds (Kreidberg et al., 2015; Vidal-Madjar et al., 2011; Louden & Wheatley, 2015). The recent work of Sing et al. (2016) has also begun to classify these objects spectrally, as has been done for stars. The proximity of hot Jupiters to their parent star suggests that tidal forces are likely to rapidly evolve hot Jupiters into a tidally-locked state, with permanent day and night hemispheres (Baraffe et al., 2010). Therefore, although much progress has been made (and much more can yet be made) using 1D models, 3D models are required to truly unpick the observations, and extract robust physical meaning.

Several GCMs (or similar models) with varying levels of sophistication have been applied to hot Jupiters (see for example Cooper & Showman, 2005; Cho et al., 2008; Menou & Rauscher, 2009; Rauscher & Menou, 2010; Heng et al., 2011; Dobbs-Dixon & Agol, 2013; Parmentier et al., 2013; Showman et al., 2015; Helling et al., 2016; Kataria et al., 2016; Lee et al., 2016), including our own adaptation of the Met Office GCM termed the Unified Model (UM) (Mayne et al., 2014a, b; Amundsen et al., 2014; Helling et al., 2016; Amundsen et al., 2016, 2017; Boutle et al., 2017). However, much of the progress has been driven by application of a single GCM, the SPARC/MITgcm. In fact, a significant number of subsequent, detailed analyses of hot Jupiter atmospheres have been performed using output, or derivatives of the pioneering work of Showman et al. (2009). Studies based on the SPARC/MITgcm include exploring the presence of TiO/VO, advection driven non–equilibrium chemistry and changes in the dynamics over a range of hot Jupiters (Parmentier et al., 2013; Agúndez et al., 2014; Kataria et al., 2016). Despite this progress lessons learned from the Earth and solar system communities tell us that GCM results can be particularly model-dependent (e.g. Lebonnois et al., 2011). Simple model intercomparisons have been done (Heng et al., 2011), and efforts have also started to compare more complex models (Helling et al., 2016). In Amundsen et al. (2016) we presented simulations of HD 209458b using our own adapted GCM which is of commensurate sophistication to that of Showman et al. (2009). Amundsen et al. (2016) present qualitative and quantitative differences between their results and those of Showman et al. (2009), emphasising the need for further intercomparison of both the GCMs themselves and post-processing tools111We are working on a more direct comparison between the UM and SPARC/MITgcm, but such comparisons are difficult for models as complex as GCMs..

Previous results from GCM simulations of hot Jupiters hint at possible initial condition sensitivity, if extensive damping is not used at the bottom boundary (Liu & Showman, 2013; Cho et al., 2015), which could be caused by the deep, radiatively unforced atmosphere (hereafter the ‘deep’ atmosphere refers to the region where the pressure is in excess of $10^{6}$ Pa or 10 bar) dynamically evolving similarly to the oceans of Earth (Mayne et al., 2014a). Amundsen et al. (2016) also note an evolution of the deep atmosphere temperature–pressure profile in their simulations of HD 209458b. We have recently explored the evolution of the deep atmosphere using a 2D model assuming a steady state solution, the results of which will be discussed in an upcoming publication (Tremblin, et al., submitted).

Despite the possible uncertainties, several features are qualitatively, yet robustly reproduced, across most of the GCMs applied to hot Jupiters, in particular the occurrence of a super-rotating equatorial jet (coherent zonal flow). This jet is seemingly confirmed by observations of a shift in the brightest part of a hot Jupiter atmosphere away from the substellar point, observed via phase curves and suggested to be caused by wind-driven advection (Knutson et al., 2007; Zellem et al., 2014) by a super-rotating equatorial jet. Although, simulations suggest the cause is more subtle than this and that the initial temperature structure is setup by fast moving waves driven by the irradiation, subsequently acting to accelerate the jet itself (Showman & Polvani, 2011).

Although jets are commonplace amongst Solar system objects, and simple 2D arguments can be used to characterise their likely breadth in latitude (see reviews by Showman et al., 2008, 2011), analytical description of the mechanisms which accelerate the jet by the convergence of prograde momentum (or pumping), is much more challenging. The difficulty is that the diagnosis of such mechanisms requires an understanding of the interactions of the mean flow with eddies or perturbations. Despite this difficulty the mechanism for the acceleration and support of Earth’s mid–latitude jets is relatively well understood. The result comes directly from the conservation of wave activity (Vallis, 2006). Atmospheric Rossby waves, or vortices, excited by baroclinic instability at mid–latitudes, travel towards the pole (or equator). The characteristics of the Rossby wave dispersion relation lead to the transport of eastward (or positive) angular momentum into the excitation location and westward (negative) angular momentum into the dissipation site, which in the case of Earth is at low latitudes and at high latitudes. The total angular momentum in the system is conserved, and this is why Earth’s mid–latitude jets are coupled to the circumpolar and equatorial retrograde flow (Vallis, 2006). This effect can be diagnosed by exploring the Eliassen-Palm flux, a vector quantity representing the relative strengths of the eddy heat and momentum fluxes (Eliassen & Palm, 1960, 1961; Vallis, 2006).

Showman & Polvani (2010) explore the mechanism for generating super-rotating flows at the equator, extending the 2D analytical results of Matsuno (1966) (representing positive and negative heating on opposing hemispheres) and Gill (1980) (representing positive and zero forcing on opposing hemispheres). Showman & Polvani (2010) performed tests using a simple shallow–water model, simulating superrotation at Earth’s tropics and with the addition of momentum exchange between the upper, ‘active’ atmosphere and lower, ‘quiescent’ atmosphere super-rotating flows were accelerated. This indicates strongly that the mechanism which exchanges angular momentum vertically in an atmosphere is critical, in balance with the horizontal interactions, to the generation of a super-rotating jet at the equator, or in other words the process is truly a 3D one. As discussed by Showman & Polvani (2011), however, an Earth–like mechanism is unlikely to operate in the atmospheres of hot Jupiters (see discussion in Showman & Polvani, 2011). Showman & Polvani (2011) develop an alternative theory, still involving mean flow–eddy interaction, but reliant upon planetary scale, equatorially trapped, standing Rossby and Kelvin waves (which are responses to the large scale forcing as explained in Showman & Polvani, 2011). This mechanism, as for the previous more Earth-like case, relies on a balance between the vertical and meridional eddy momentum fluxes in the atmosphere, with angular momentum being extracted from a deep atmosphere reservoir (Showman & Polvani, 2011). Showman & Polvani (2011) suggest that the jet in their simulations of HD 209458b, is accelerated in the first few tens of days (all references within this work to days refer to Earth days i.e. 86 400 seconds) and reaches an equilibrium where the acceleration terms balance to zero at the equator. The equations and simulations of Showman & Polvani (2011) are, however, based on the primitive equations of motion (a simplified version of the equations of motion for an atmosphere, see Mayne et al., 2014a, for a full discussion in relation to hot Jupiters), incorporating the assumption of vertical hydrostatic equilibrium, a shallow atmosphere and constant gravity. Further testing of this mechanism has been performed using results from a more dynamically complete model by Tsai et al. (2014). Tsai et al. (2014) adopt the $\beta$-plane approximation and investigate the momentum fluxes in 3D, in an assumed steady state. In the $\beta$-plane approximation the Coriolis force is taken to vary linearly with latitude, $\phi$ (actual variation is $\propto\cos\phi$). Tsai et al. (2014) found a picture consistent with the ideas of Showman & Polvani (2011).

In this work we explore the dynamical form and acceleration of the super-rotating equatorial jet within various simulations based nominally on HD 209458b. A key feature of this study is that our model is a non-hydrostatic, deep atmosphere model, and we do not invoke a drag at the bottom boundary. The layout of this paper is as follows: in Section 2 we detail the model used and the simulations we have performed, referring to previous works for the details. Section 3 then details our results. Firstly, in Section 3.1 we present the morphology of the dominant, zonal advection over long simulation timescales, at various levels of approximation to the dynamical equations, treatments of the heating and treatment of the deep atmosphere thermal profile. Our results reveal a robust recurrence of the equatorial super-rotating jet. Then in Section 3.2 we explore the evolution of the deep, high pressure atmosphere, highlighting the gradual evolution of the kinetic energy and the fact that it has not reached a steady state in any of our simulations. In order to test the effect of this on the flows in the upper atmosphere (i.e. at low pressures), we perform a simple limiting numerical experiment, strongly forcing the deep atmosphere over long timescales, revealing an eventual weakening of the super rotating equatorial jet. Finally, in Section 3.3 we explore the acceleration and maintenance of the zonal flows. Section 3.3.1 presents the initial response (from rest) of the simulated atmosphere to the heating, featuring similar eddy patterns as observed in Matsuno (1966); Gill (1980); Showman & Polvani (2011); Tsai et al. (2014). Then, in Section 3.3.2 gradients of the momentum fluxes in the simulated atmospheres are used to reveal the mechanisms maintaining the zonal flows. We find similarities between momentum transport, and the previous work of Showman & Polvani (2011) and Tsai et al. (2014) but highlight additional complexities which merit further study beyond the scope of this work (Debras et al., in prep). In particular, the mean flow contributions to the momentum transport are non-negligible. Our conclusions are stated in Section 4. Finally, the Appendix contains extended detail of the flow structure (Appendix A), results from a ‘shallow-hot Jupiter’ test (Appendix B), the derivation of the eddy-mean flow interaction equation (Appendix C) and plots of additional terms within this equation for our simulations (Appendix D).

2 Model setup

Results from this work are from simulations using the basic setups of Mayne et al. (2014a) and Amundsen et al. (2016), which are differentiated by the treatment of heating/cooling and we refer to these papers for the basic numerical details. In Mayne et al. (2014a), the basic equations solved by the UM are introduced and results from simulations including a simple Newtonian relaxation, or radiative forcing scheme are presented, where the temperature is simply relaxed to a pre-calculated radiative equilibrium profile (hereafter termed TF simulations referring to ‘Temperature Forced’). This work also defines the various levels of simplification to the dynamical equations termed ‘primitive’, ‘shallow’, ‘deep’ and ‘full’ in order of the most simplified to the most complete. Essentially, the most approximated equations are the ‘’primitive’, assuming vertical hydrostatic equilibrium, gravity constant with height and a shallow–atmosphere. The ‘shallow’ equations are formed by relaxing the hydrostatic approximation, and the ‘deep’ by further relaxing the shallow–atmosphere approximation. Finally, the ‘full’ equations invoke none of these approximations. (see Mayne et al., 2014a, for details). Amundsen et al. (2016) features models using only the ‘full’ equations but incorporating a two-stream, dual band, correlated-k radiative transfer scheme, which has previously been tested in Amundsen et al. (2014, 2017), alongside minor updates in the treatment of diffusion and minor parameter changes (hereafter termed RT, referring to ‘Radiative Transfer’, simulation). The model uses SI units, and these are adopted throughout this work.

Although TF simulations are computationally much cheaper than their RT counterparts allowing us to reach extensive total integration times, they have several key disadvantages. As TF simulations linearly relax the temperature to a prescribed profile, they lack the inclusion of atmospheric interactions (i.e. thermal emission and absorption), and do not model the non–linear response of the atmosphere to significant levels of heating and cooling both of which are captured in the RT simulations. Additionally, the TF setups rely on initial calculation of equilibrium profiles, often from simple and physically incomplete 1D models, further reducing their accuracy and also flexibility. Our results exhibit clear differences between the TF and RT simulations (see Section 3), some elements of which will be caused by such issues.

2.1 Damping

Of particular note is the treatment of damping in the model. In a physical fluid eddies and turbulence will act to cascade energy to smaller scales, and ultimately convert kinetic energy into thermal energy (see discussion in Li & Goodman, 2010). Hot Jupiter GCM simulations do not correctly capture such small scale turbulence. Additionally, physical processes which are not explicitly modelled can also lead to further dissipation, for example magnetic braking in hot Jupiters (Menou, 2012). Explicit damping is often applied to account for these effects. However, we currently have no way of constraining the strength of the dissipation (Li & Goodman, 2010; Heng et al., 2011), except perhaps by matching simulated offsets in hot Jupiter hotspots from the substellar point with observed phase curve offsets, although processes not explicitly modeled complicate this (Zellem et al., 2014). Therefore, such damping and dissipation parametrisation are primarily used to achieve numerical stability in a physically plausible way, but are not robustly constrained. In particular the low pressure atmosphere, near the outer boundary is sensitive to perturbations which grow as they travel from higher pressure regions, and become unstable. For the inner, high pressure, boundary dissipation from the numerical scheme (i.e. discretisation) is more significant, and the fluid more stable. However, as the pressure (and therefore density) increases exponentially with depth into the hot Jupiter atmosphere, the angular momentum for a given flow rate will increase. Therefore, even very slow fluid flows at the inner boundary can significantly affect the angular momentum budget (Mayne et al., 2014a), and minor inaccuracies of the numerical solver can drastically affect the simulated dynamical structure (Cho et al., 2015).

In our model we include several forms of artificial dissipation. Firstly, as detailed in Mayne et al. (2014a) and updated in Amundsen et al. (2016) we include an explicit diffusion of the zonal flow (only). However, we incorrectly reported the coefficient for this scheme previously, in Amundsen et al. (2016), stating a single coefficient for the zonal diffusion operator ( $K_{\lambda}\sim 0.16$). Actually, due to simplification in the routine performing the diffusion several terms are implicitly included in the resulting applied diffusion coefficient, $K^{\rm\prime}_{\lambda}$, such that it is given by

[TABLE]

Where $\Delta t$ and $\Delta\lambda$ are the timestep and longitudinal grid spacing, respectively 222In practice the diffusion is set using a characteristic e–folding timescale, $K_{\lambda}=0.25\left(1-e^{\frac{-1}{t_{K}}}\right)$, where $t_{K}$ is the set number of simulation timesteps.. We also include a vertical ‘sponge layer’ as detailed in Mayne et al. (2014a) to damp vertical velocities close to the outer, low pressure boundary, and represent vertical waves propagating out of our modelled domain. For this work, to allow us more flexibility, we have also included a vector laplacian damping (to mimic an explicit viscosity) in the momentum equation. The vector Laplacian of the vector wind field is calculated, $\nabla^{2}(\mathbf{u})\equiv\nabla\left(\nabla\cdot\mathbf{u}\right)-\nabla\times\left(\nabla\times\mathbf{u}\right)$. The resulting components of this provide additional terms to equations (1), (2) and (3) of Mayne et al. (2014a), governing the $u$, $v$ and $w$ wind components, which are in the longitude ($\lambda$), latitude ($\phi$) and vertical directions ($r$), respectively. A separate multiplicative coefficient is then prescribed for the horizontal and vertical directions as $\frac{\nu_{\lambda,\phi}}{\rho}$ and $\frac{\nu_{\mathrm{r}}}{\rho}$, respectively, where $\rho$ is the density.

In effect the diffusion and vector Laplacian schemes apply the same, physical equation, but differ via their calculation of the coefficient. In practice neither of these schemes operate on the vertical component of velocity. Given our spatial resolution ($\sim 10^{5}$ & $\sim 10^{6}$ m in the vertical and horizontal directions, respectively), and maximum windspeeds ($\sim 10^{2}$, $\sim 10^{3}$ ms*-1*, in the vertical and horizontal directions, respectively), the vertical flow is much more accurately resolved. However, the vector laplacian scheme is applied to both the zonal and meridional directions, whereas the diffusion scheme is only applied in the zonal direction (see discussion in Amundsen et al., 2016). Given typical simulation parameters ($\Delta t\sim 120$ s, $r\sim 10^{8}$ m, $\Delta\lambda\sim 0.04$  radians, see Table 1) and densities throughout our domain which range from $\sim$1 to 2$\times 10^{-7}$ kg m*-3*, $K_{\lambda}^{\prime}\sim\frac{8\times 10^{-12}K_{\lambda}}{\cos^{2}\phi}$ and the applied vector laplacian coefficient, $\nu^{\prime}_{\lambda,\phi}\sim(10^{-7}\rightarrow 1)\,\nu_{\lambda,\phi}$. Therefore, the maximum applied diffusion coefficient (at the equator) is typically between 12 and five orders of magnitude smaller than the vector laplacian, for the same input value. The spatial variation of the coefficient also leads to the two schemes behaving differently. For the vector laplacian, the variation of the coefficient with density effectively means this process mimics kinematic visocity and will act most strongly on the horizontal wind in the upper, low pressure regions of the atmosphere. However, the diffusion coefficient varies strongly with latitude such that it reduces towards the pole, meaning this scheme effectively acts to suppress grid–scale noise in the dominant zonal flow. In Appendix B, we demonstrate the effect of varying the ‘strength’ of the diffusion on a shallow–hot Jupiter test case (introduced by Menou & Rauscher, 2009). It is important to reiterate that although generally physically motivated, these damping/dissipation mechanisms are poorly constrained for hot Jupiters, and are primarily used to achieve numerical stability.

Additional horizontal damping at the inner or high pressure boundary is also often employed (see discussion in Cho et al., 2015). This ‘bottom drag’ is employed similar to the Rayleigh friction schemes used to model the frictional damping of the Earth’s surface on horizontal winds, and is motivated by the potential for magnetic drag in the deep atmosphere, (Rogers & Komacek, 2014). However, the strength of this drag, and vertical profile are very poorly constrained, and it is not clear if a Rayleigh type drag is appropriate at all. Of course, when employing such a bottom drag GCM simulations are much more robust to conservation of total angular momentum, and less sensitive to initial conditions (Liu & Showman, 2013; Cho et al., 2015). In this work we do not include a bottom drag damping as we are, in part, interested in the evolution of the deep atmosphere, and its interaction with the lower pressure regions.

2.2 Model variations

In Table 1 we state the main parameters for both the TF and RT simulations. These differ slightly as the TF simulations were setup to match the work of Heng et al. (2011) and the RT simulations to be compared with Showman et al. (2009) but these differences are unlikely to affect our main conclusions. We have run a large set of simulations to explore various scenarios, but only directly report results here for a subset. We have also performed additional simulations of the ‘shallow-hot Jupiter’ setup of Menou & Rauscher (2009), the results of which are discussed in Appendix B. Table 2 details the elapsed simulation time, parameters adjusted from the default (shown in Table 1), dynamical equation set (following the nomenclature of Mayne et al., 2014a) and the ‘short name’ adopted throughout this manuscript.

Only a single RT simulation is included as the TF simulations allow greater, direct, control of the heating and run significantly faster. Although TF simulations are less physically accurate than their RT counterparts (see discussion in Showman et al., 2009; Amundsen et al., 2016), they have still proved useful in the study of hot Jupiter (and other) atmospheres. Recent examples of the use of TF simulations include exploration of the efficiency of dynamical redistribution of heat (Komacek & Showman, 2016) and the dependence on the atmospheric flow on bulk composition (Zhang & Showman, 2017). For the RT simulation we do not include TiO and VO formation as evidence for its presence has only been suggested for Wasp-121b (Evans et al., 2016) so far, and not HD 209458B.

Our TF simulations have been evolved for much longer than those of Mayne et al. (2014a), allowing us to explore the long term evolution for simulations adopting various levels of simplification to the dynamical equations (e.g., Std Prim & Std Full)333Note we have also performed simulations using the “shallow” and “deep” versions.. The profiles for the heating in the TF simulations are those of Mayne et al. (2014a), which are adjusted from Heng et al. (2011) which, in turn, are based on those of Iro et al. (2005). The radiative timescale is an (approximately) exponentially increasing function of pressure, ranging from $\sim 10^{3}-10^{8}$ over pressures of $\sim 10^{2}-10^{6}$ Pa, and is infinite for pressures of $10^{6}$ Pa or higher.

As discussed in Mayne et al. (2014a) evolution of the deep atmosphere is very gradual and is likely to require extremely long simulation integrations to reach a steady state. Mayne et al. (2014a) note a gradual increase in the difference between the equatorial temperature and that of the pole for the deep atmosphere. In several of their simulations the poles are warmed, and the equator cooled, which given the absence of forcing for the standard models can only be caused by compression/expansion or advection potentially by material lifting over the hotspot and falling towards the poles. A similar temperature evolution of the deep atmosphere can be seen in the RT simulations of Amundsen et al. (2016), Figure 7. However, for the models where the initial profile is hotter, presented in Amundsen et al. (2016), the latitudinal temperature gradient is reduced. This behaviour hints that the 3D model is evolving to a hotter equilibrium state than the one with which it is initialised, however this is beyond the scope of this work, where we are more interested in the effects of different scenarios for the unconstrained deep atmosphere, and will be discussed in an upcoming publication (Tremblin, et al., submitted). Therefore, we have simply performed a numerical experiment to explore the limits of this effect, and determine whether it will be significant for the flows in the upper atmosphere. In the simulation Deep $\Delta T_{\rm eq\rightarrow pole}$ we impose an additional, latitudinal temperature gradient in the deep atmosphere using a constant radiative timescale given as the value at $10^{6}$ Pa (10 bar). For this setup although the equilibrium temperature is increased by $(1\,000.0\,K)\sin^{2}(\phi)$ (where $\phi$ is latitude) the final temperature difference between the equator and pole will be much less as this must first compensate for the existing contrast enforced by the original equations of $\sim-500$ K (Heng et al., 2011; Mayne et al., 2014a). We have also run a simulation with the bottom boundary at $10^{6}$ Pa to explore the effect of omitting the deep region entirely (Reduced $p_{\rm max}$). It is important to note that these simulation setups are explorative and performed to investigate the interaction of the deep and low pressure atmosphere. We are unlikely to obtain observational constraints on the state of the deep atmosphere of hot Jupiters in the near future and so are forced to make choices regarding how to model this region.

All of our models are initialised with zero winds, and in solid body rotation with a hydrostatically balanced atmosphere. The initial temperature pressure profile used is either a midway profile between the hottest and coldest profiles for the TF simulations, or a globally averaged radiative equilibrium profile for the RT simulations (for further details refer to Mayne et al., 2014a; Amundsen et al., 2016), unless otherwise stated.

3 Results

In this section we present results for the dynamical state of our simulated atmospheres. In Section 3.1 we demonstrate that the occurrence of a super rotating equatorial jet is robust in our simulations, and present the flow regime for some examples. In Section 3.2 we demonstrate that the deep atmosphere is still evolving in our simulations, and show that this effect may disrupt the flow regime in the upper atmosphere, but only if strongly forced over long timescales. Finally, in Section 3.3 we present the mean flow and eddy interaction, demonstrating that the jet in our simulations is driven by a mix of eddy and mean flow momentum transport.

3.1 Zonal flow

Figure 1 shows the zonal and temporally averaged (mean) zonal wind, where red is prograde and blue is retrograde for the TF simulations where the completeness of the dynamical equations solved is varied (Std Prim and Std Full, left and right columns, respectively), at two different epochs (200–1 200 & 9 000–10 000 days as the top and bottom rows, respectively). It is clear that the jet evolves over the period of about 10 000 Earth days ($\sim$3 000 rotation periods). The breadth, in latitude, and depth range, as well as the peak prograde velocity all increase. Additionally, the morphology of the jet changes with the level of assumptions used in the dynamical equations. For the more simplified case the jet is generally slower and covers a narrower range of pressures, whilst being slightly broader in latitude444Results from the simulations of intermediate complexity (i.e. ‘shallow’ and ‘deep’ are consistent with this trend.. The difference found when moving from ‘full’ to ‘primitive’ is similar to what one might expect when reducing the rotation rate. However, despite these changes it is clear that, broadly and qualitatively speaking, the same jet structure is apparent; a prograde equatorial jet of a few kms*-1* flanked by retrograde jets, covering a broad range of pressures down to about $10^{6}$ Pa.

Figure 2 shows the results for the zonal flow for the Std RT simulation, in the same format as Figure 1, but for 600-1 600 days. However, the Std RT simulation bottom boundary is placed at lower pressures, and this combined with the subsequent evolution means the highest available pressure is slightly lower than that of the Std Prim or Std Full cases. Additionally, the much shorter total elapsed simulation time, compared to the TF simulations, means that the deeper atmosphere is likely to still be evolving. In fact, as we demonstrate in Section 3.2, and present in Figure 9, the kinetic energy of the atmosphere appears to still be evolving at pressure higher than $\sim 10^{6}$  Pa. As discussed in Showman et al. (2009) and Amundsen et al. (2016), moving from a temperature forcing to more complete radiative transfer scheme alters the simulated dynamics of the atmosphere. For the RT simulation the jet is slower and maintains a broader latitude profile at lower pressures when compared to the most comparable standard TF simulation (from simulations not presented here a similar result is found when reducing the gravity, or slowing the rotation speed of the planet). However, the qualitative result is still consistent with that of the Std Prim and Std Full simulations.

The set of results presented here, alongside further simulations (not presented), show broad qualitative agreement and are consistent with previously published results (e.g. Showman et al., 2009; Heng et al., 2011; Dobbs-Dixon & Agol, 2013), in that all produce a prograde equatorial jet. However, it is also clear that various parameter choices can affect the flow, and create differences which can be seen even in a relatively coarse measure such as the zonally and temporally averaged zonal wind plots

Figures 3 and 4 show the temperature (colour scale and contours) and horizontal wind (vector arrows) structure, after 1 200 and 10 000 days respectively. The flow is depicted on the two highest pressure isobaric surfaces as used in Heng et al. (2011) and subsequently Mayne et al. (2014a), namely 4.69$\times 10^{5}$Pa, and 21.9$\times 10^{5}$ Pa, as the top and bottom rows, respectively. The corresponding Figures for the lower pressure surfaces (213 and 21 600 Pa) are presented in Appendix A (and reveal little variation across time, or simulation setup). After 10 000 days the Std Prim simulation has a smooth, and relatively time invariant structure dominated by the homogenisation of temperature around the equator. However, the Std Full simulation, whilst still supporting a similar equatorial flow maintains vortices at 4.69$\times 10^{5}$Pa, and a more complex flow structure at 21.9$\times 10^{5}$ Pa. For the Std Full simulation the deepest isobaric slice shows time evolution, even after 10 000 days, as demonstrated by Figure 5, where this isobaric surface is shown at two further times, 8 000 and 13 000 days. The behaviour in this region is purely dynamical, driven by both circulations and adiabatic compressions/expansions from the upper atmosphere since the forcing timescale is infinite below $10^{6}$ Pa (Heng et al., 2011; Mayne et al., 2014a). The deepest layer of the Std Full simulation, 21.9$\times 10^{5}$ Pa after 10 000 days shows some asymmetry about the equator (see Figure 4). However, this is due to fluctuations of the still evolving fluid, as can be seen in Figure 5.

Figure 6, shows the same information as Figure 3, but for the Std RT simulation, after 1 600 days. Again the lower pressure surfaces are shown, and discussed in Appendix A. The Std RT ‘snapshots’ have been taken at the latest simulation time, 1 600 days, to capture the atmosphere in its most evolved state, as opposed to matching the 1 200 days of the first ‘snapshot’ of the Std Prim and Std Full simulations. The main differences between the thermodynamic and dynamical structure between RT and TF simulations have been discussed in the previous works of Showman et al. (2009) and Amundsen et al. (2016). Here we note that the deep atmosphere in this case is qualitatively different to the TF simulations. In particular, the slice at 21.9$\times 10^{5}$ Pa shows an equator to pole temperature gradient of $\sim+300$ K. As discussed in Amundsen et al. (2016) and Section 2.2 this may also be due to the 3D model adjusting from an incorrect, or inconsistent initial profile to a steady state (discussed further in, Tremblin et al., submitted).

3.2 Deep atmosphere evolution

As discussed in the previous section, as well as Mayne et al. (2014a) and Amundsen et al. (2016), the dynamical (and thermodynamic) state of the deep atmosphere is slowly evolving throughout the simulations (even out to several thousand rotation periods). Additionally, although the flows in upper, lower pressure atmosphere are broadly consistent across our simulations, as one would expect due to the strong forcing, the flows do vary for the deeper, higher pressure regions. Clearly, the evolutionary timescale for the deep atmosphere is likely to be very long. The thermal timescale, $\tau_{\rm th}$, can be estimated using (see for example Geroux et al., 2016),

[TABLE]

where $\Delta m$ is the mass in the layer, and $L$ the luminosity (estimated from $L=4\pi r^{2}\sigma T_{\rm int}^{4}$, where $\sigma$ is the Stefan-Boltzmann constant). Assuming a characteristic temperature of $\sim 100$ K, the timescale is $\sim 10^{4}$ years, at pressures of $10^{7}$ Pa. As the pressure, and therefore density, increase exponentially moving to higher pressures, this timescale will also increase exponentially (as the mass is linearly dependent on the density). Therefore, we have performed two numerical experiments to explore the possible impact of the deep atmosphere on the lower pressure regions accessible to observations, by artificially forcing, or omitting the deep regions entirely.

As demonstrated by Showman & Polvani (2011) equatorial superrotation can be achieved, in the context of a tidally-locked hot Jupiters, by the liberation of angular momentum from a deep atmosphere ‘reservoir’. Additionally, Mayne et al. (2014a) highlight an evolution in the equator-to-pole temperature difference of the high pressure atmosphere, in their simulations. Therefore, we have performed simulations omitting the radiatively inactive deep atmosphere, termed Reduced $p_{\rm max}$, or enforcing a strong equator-to-pole temperature gradient, termed Deep $\Delta T_{\rm eq\rightarrow pole}$. Figure 7 shows the zonal flow, in the same format as Figure 1, for the Deep $\Delta T_{\rm eq\rightarrow pole}$ and Reduced $p_{\rm max}$ (note the y-axis only extends to pressures of $\sim 10^{6}$ Pa or 10 bar) simulations, as the left and right columns, respectively. Interestingly, the absence of the deep atmosphere section does not disrupt the formation of a prograde equatorial jet. However, for the simulation forcing the equator-to-pole temperature gradient after 10 000 days the jet has significantly slowed and covers a much smaller region in both pressure and latitude. The angular momentum is balanced in this model by the prograde jets around the poles at high pressures, which are much faster than in the Std Full case (see Figure 1). Almost all of the simulations we have performed conserve angular momentum to better than $\sim$5% over periods of around 10,000 Earth days (most are better than 1%). The few exceptions include the simulation with the bottom boundary at $10^{6}$ Pa (10 bar) where the angular momentum almost doubles, and all are discussed later in this section.

Figure 8 shows the isobaric slices for the same pressures, and in the same format as Figure 4, but for the Deep $\Delta T_{\rm eq\rightarrow pole}$ simulation. The lower pressure slices are again presented in Appendix A, as Figure 19. The corresponding figures for the Reduced $p_{\rm max}$ are omitted as for the available pressures the results closely match those of the Std Full simulation in terms of morphology. For the Deep $\Delta T_{\rm eq\rightarrow pole}$ simulation the flow is strongly retrograde at mid to high latitudes, with a weak equatorial prograde flow, and subsequently weaker homogenisation of the temperature about the equator compared to the standard simulations (compare the 21 600 Pa slice from Figures 17 and 19, as well as the 4.69$\times 10^{5}$ Pa slice from Figures 4 and 8). Deeper in the atmosphere a temperature difference of $\sim+300$ K is achieved between the equator and the pole, similar to that seen in the Std RT model (compare Figures 8 and 6), and we see an extension of the retrograde flow.

To explore the evolution of the deep atmosphere in our simulations, Figure 9 shows the logarithm of the total kinetic energy ($\log_{10}$(KE [J]), including the atmospheric mass) as a function of time and pressure for the Std Prim, Std Full, Deep $\Delta T_{\rm eq\rightarrow pole}$ and Std RT simulations as the top left, top right, bottom left and bottom right panels, respectively, following the figure presented previously in Rauscher & Menou (2010). Note that the Std RT simulation only extends to 1 600 days, whereas the others are presented to 10 000 days. It is clear that the upper or lower pressure atmosphere accelerates very rapidly, with the KE growing rapidly and then reaches an almost steady value. The deeper atmosphere seems to evolve much more slowly in the Std Prim case than that of the Std Full simulation (although both are clearly far from completing their evolution after 10 000 days). The deep atmosphere of the Std RT simulation also appears to evolve very gradually, but of course the total simulation time here is much shorter. Mayne et al. (2014a) concluded that the deep atmosphere had quickly finished evolving in a simulation assuming a ‘shallow’ atmosphere (i.e. similar to our Std Prim case, but the assumption of vertical hydrostatic balance is relaxed), using the deep atmosphere equator to pole temperature contrast. However, Figure 9 actually shows that this might just be caused by an evolution which is much slower than the Std Full case. Additionally, it is clear that the KE peaks at around 1-10 bar, $1$–$10\times 10^{5}$ Pa except for the case of the Deep $\Delta T_{\rm eq\rightarrow pole}$, where the energy budget is dominated by the material at higher pressures, as one might expect given the flow morphology. Although this experiment of artificially forcing the deep atmosphere is somewhat arbitrary, it does illustrate the possible effects of deep atmosphere flows on the upper atmosphere. The final critical point is that the evolution of the deep atmosphere appears to require very lengthy elapsed simulation times, using Figure 9 as a guide. Therefore, it is quite likely that the final solutions presented by GCMs after current published simulation times may well be dependent on the initial state of the deep atmosphere, both dynamically, and thermodynamically (see discussion in Mayne et al., 2014a; Amundsen et al., 2016). However, it must be noted, our artificial forcing of the deep atmospheres thermodynamic state is quite strong (adopting the timescale assumed for $10^{6}$ Pa), and over 10 000 days has still not completely destroyed the prograde equatorial jet.

For the simulations presented in Figure 9 the total KE is still increasing towards the end of the simulation time, as one might expect given the increase in the deep atmosphere component. As the pressure, and therefore density, increases exponentially as one moves closer to the inner simulation boundary, a significant KE contribution can be achieved with very slow velocities. Generally, the maximum zonal velocity ($u_{\rm max}$) is used to determine a quasi-steady state in the literature for hot Jupiter simulations. Figure 10, left panel, shows $u_{\rm max}$ as a function of time for the Std Full, Std Prim, Deep $\Delta T_{\rm eq\rightarrow pole}$ and Reduced $p_{\rm max}$ simulations. As Figure 10 shows, the fact that the KE has not equilibrated is not detected in the maximum zonal velocity, with all of these simulations appearing to reach a quasi-steady solution. This is, of course driven by the fact that the deep atmosphere velocities are much slower than those in the upper, low pressure, regions. This suggests that assuming a steady state based on the maximum zonal velocity will not apply to the high pressure atmosphere.

The deep atmosphere also represents a huge reservoir of axial angular momentum (AAM), which can be used to accelerate zonal flows. Updrafts of material can effectively convert significant amounts of planetary angular momentum (simply due to the solid body rotation of the atmosphere) into atmospheric angular momentum (from the winds) (see discussion in Mayne et al., 2014a), as the total (i.e. the sum of the planetary and atmospheric components) is conserved. Additionally, as with the KE even very small velocities in the high pressure regions can make significant contributions to the total AAM budget. The right panel of Figure 10 shows the total normalised AAM, relative to the initial value, for the same simulations presented in the left panel of Figure 10, namely the Std Full, Std Prim, Deep $\Delta T_{\rm eq\rightarrow pole}$ and Reduced $p_{\rm max}$ simulations. The Reduced $p_{\rm max}$ shows a near doubling of AAM in the first $\sim$4 000 days. The Std Prim and Deep $\Delta T_{\rm eq\rightarrow pole}$ simulations both gain around 15% AAM after about 10 000 days, which is still a reasonable level of accuracy for such low resolution simulations run over such long periods. All of our remaining simulations conserve AAM to better than 5%, and in many cases 1%, and the Std Full model is presented as indicative of this. If angular momentum conservation issues were dominated solely by the fast dynamics in the upper atmosphere, one would expect the Reduced $p_{\rm max}$ simulation to show much poorer absolute AAM conservation. This is as this simulation does not include the vast AAM contribution from the deep atmosphere. The absolute change in AAM for the Reduced $p_{\rm max}$ simulation is $\sim 13.5\times 10^{32}$kg ms*-1* (from an initial total AAM of $\sim 15.4\times 10^{32}$kg ms*-1*), whereas the absolute change for the Std Full simulation is a factor $\sim$5 lower, $\sim 2.6\times 10^{32}$kg ms*-1* (from an initial total AAM of $\sim 33.9\times 10^{33}$kg ms*-1*) indicating that indeed the inaccuracies in the AAM conservation are dominated by the fast winds in the low pressure atmosphere. However, Cho et al. (2015) noted significant changes in AAM using models without bottom boundary drag, implying that velocities close to the bottom boundary lead to numerical losses. Here we do not employ a bottom boundary drag, and allow our deep atmosphere to slowly accelerate. For our Reduced $p_{\rm max}$ simulations fast winds are accelerated much closer to the bottom, or inner, boundary than in the other simulations. We are still investigating the reasons for this issue.

3.3 Jet pumping

Two of the main mechanisms for accelerating jets in planetary atmospheres are either the exchange of momentum via mean flows, or the convergence of momentum from the creation and destruction of eddies (see Showman et al., 2013, for a review in the terrestrial planet case). For hot Jupiter atmospheres the expected large scale of the atmospheric eddies (i.e. Rossby radius of deformation), and zonally asymmetric forcing (i.e. day and night side) led Showman & Polvani (2011) to propose an eddy momentum acceleration mechanism for the jets relying on the interaction of standing Kelvin and Rossby waves. Tsai et al. (2014) explored the validity of this mechanism in their simulations by presenting, in 3D adopting the $\beta$-plane (where the Coriolis force is simplified as varying linearly with latitude, instead of $\propto\cos\phi$), the eddy momentum fluxes as a function of pressure and a proxy for latitude. The simulations, and analysis of Showman & Polvani (2011) are performed in a pressure–based (hydrostatic) framework. Although in previous sections we have interpolated our prognostic variables onto pressure surfaces, the eddy momentum fluxes can not easily be translated in such a fashion. This is as the terms involve gradients, in the height–based framework, at constant $r$ in the longitude, $\lambda$, direction, which cancel to zero when zonally averaged. Simply interpolating the variables onto pressure surfaces invalidates this step, and the terms no longer vanish, introducing extra terms into the eddy–mean interaction equation. The derivation of this equation is explicitly stated in Appendix C, and this point highlighted. Therefore, to analyse the transport of momentum in the atmosphere we shift to using height, $z$, as a vertical coordinate. The simulations of Tsai et al. (2014), however, are performed using a height–based model, and the results interpolated onto a pressure grid, meaning they have omitted the extra contributions to the equation derived from the, resulting non–zero, zonally averaged vertical gradient terms.

3.3.1 Initial Response from Rest

Firstly, Figure 11 presents the evolution of the zonally summed linear zonal momentum ($\rho u$) for the Std Full and Std RT simulations, as the left and right columns, respectively. The top, middle and bottom rows show the zonal momentum after 1, 10 and 1 000 days, respectively.

Figure 11 shows a broadly similar evolution in the zonal momentum of the Std Full and Std RT simulations, when accounting for the vertical offset. Each model covers a slightly different pressure range, meaning the height axis will represent slightly different pressure regimes. The initial acceleration occurs very rapidly, as suggested by Showman & Polvani (2011), with a jet structure appearing in the first day. After around 10 days the acceleration of the flows in the low pressure, radiatively forced regions is largely complete for the Std TF simulation, with slow evolution out to 1 000 days. After 1 day the acceleration is broader and more localised to low pressures, for the Std RT simulation, however this may well simply be representative of transient evolution from slightly different initial conditions. The Std RT simulation also undergoes slower acceleration over the first 10 days. The later stages are broadly similar, but the Std RT presents a smoother jet profile. For the Std Full simulation, Figure 12 shows the same data as Figure 11 but after 10 000 days, near the end of our simulation. The overall structure of the zonal momentum has not dramatically changed, from 1 000 days, however the zonal momentum both in the jet (prograde), and at very high pressures (retrograde), has significantly increased.

3.3.2 Pseudo–Steady State Eddy Momentum Transport

For the rest of this section we focus on the eddy, and mean flow interaction in our Std Full and Std RT simulations. The work of Showman & Polvani (2011) built on the analytical solutions of Matsuno (1966) and Gill (1980). Matsuno (1966) and Gill (1980) studied the linear response (steady state) of an atmosphere at rest, in 2D, to non-axisymmetric heating as present in tidally-locked planets. Matsuno (1966), however, included a positive and negative heating, on opposing ‘hemispheres’, whereas Gill (1980) considered only the positive forcing, analogous to the “day side” heating. These differences yield a subtly different solution. In order to compare the structure of the perturbations in our eddies with the previous work of Showman & Polvani (2011), Matsuno (1966) and Gill (1980), we decompose the flow into a zonally–averaged mean flow and a perturbation from this mean. For example, $u\rightarrow\overline{u}+u^{\prime}$, where $\overline{u}$ is the zonally averaged zonal velocity. Note, these perturbations are departures from a zonal mean, not from a temporal mean flow.

Firstly, we explore the linear response of the atmosphere, from its initial rest state, to the forcing, after the first day. The work of Showman & Polvani (2011) incorporates simulations on a vertical pressure grid, whereas our model adopts height as the vertical coordinate (see discussion in Hardiman et al., 2010). Therefore, we interpolate our prognostic variables onto pressure surfaces to aid comparison of the eddy structures. However, for the divergence of the eddy momentum fluxes, discussed earlier in this section, such a transformation is much more complex so we retain our native height–based coordinate system.

Figure 13 shows the temperature and horizontal wind velocity (top row) alongside their eddy components (bottom row), as the contours and vectors, respectively, for the Std Full and Std RT simulation after 1 day (left and right columns, respectively), on the 30 mbar 3 000 Pa isobaric surface (approximately the infrared photosphere, Showman & Polvani, 2011). Figure 13 shows a similar picture for both the Std Full and Std RT simulations. Vortices, as presented in Showman & Polvani (2011) are seen either side of the equator, both to the west and east of the hot spot, with perturbations tilting westward, and thereby driving momentum to the equator (Vallis, 2006; Showman & Polvani, 2011). The pattern is apparent even after 1 day which is much faster than the advective timescale for a $\sim$kms*-1* wind traversing a planet with a radius of $\sim$108m. This indicates, as discussed in Showman & Polvani (2011) that fast waves rapidly setup the perturbation pattern which in turn accelerates the zonal flow. However, as we are starting from an unphysical, rest state, it is not clear that this activity is relevant to the study of hot Jupiters. The perturbations also show similar features around the equator, where peaks and troughs in the zonal velocity can be seen, corresponding to the equatorially trapped Kelvin or Rossby waves previously identified (Matsuno, 1966; Gill, 1980; Showman & Polvani, 2011)555Note Matsuno (1966); Gill (1980) plot pressure perturbations, whereas we plot temperature perturbations. However, we have checked density perturbations from our simulations, and the patterns match those in the temperature structure..

Figure 13 shows that the linear response of our simulations, under the radiative forcing, is similar to previous results (e.g. Showman & Polvani, 2011). However, given that our simulations are accelerating from a non-physical initial condition, it is not clear what relevance this has to real hot Jupiter atmospheres. Therefore, we now move to diagnosing the divergence of the eddy momentum fluxes, which tracks the transport of momentum by eddies in the simulated atmospheres. As discussed in Section 3.2 the deeper layers of our atmosphere are still evolving, therefore, we restrict this analysis to upper, low pressures regions. Diagnosing the momentum fluxes in the atmosphere is done using the eddy-mean flow interaction equation (see for example Hardiman et al., 2010). Showman & Polvani (2011) explore the meridional and vertical eddy terms of this equation in its primitive hydrostatic form. For our purposes, Hardiman et al. (2010) derive the same equation in height–based coordinates for our “Full” equation system. The derivation, and different forms of this equation are included in Appendix C, but for the simulations we analyse here we adopt (see Equation 23)

[TABLE]

where $r$ is radial position and $G_{\lambda}$ represents body forces acting in the zonal direction. The subscript preceded by a comma denotes a partial derivative, e.g,. $(X)_{,r}\rightarrow\frac{\partial X}{\partial r}$. As mentioned previously in this section we present results in our height–based system, as conversion to a pressure–based framework results in much more complex terms in the eddy–mean flow interaction equation.

In the mechanisms proposed by Showman & Polvani (2011) the main balance, over the equator, is between the meridional ($-\frac{\left[\overline{(\rho v)^{\prime}u^{\prime}}\cos^{2}\phi\right]_{,\phi}}{r\cos^{2}\phi}$) and vertical ($-\frac{\left[\overline{(\rho w)^{\prime}u^{\prime}}r^{3}\right]_{,r}}{r^{3}}$) eddy momentum transport terms i.e. those containing $v^{\prime}$ and $w^{\prime}$ (and gradients in the latitudinal, $\phi$ and the vertical, $r$, directions, respectively). The remaining important terms can be classified as mean flow terms involving the zonal velocity paired with either the meridional or vertical flow ($-\frac{(\overline{\rho v}\,\overline{u}\cos^{2}\phi)_{,\phi}}{r\cos^{2}\phi}$, $-\frac{(\overline{\rho w}\,\overline{u}r^{3})_{,r}}{r^{3}}$), or the mean Coriolis terms derived from the contributions of the Coriolis terms of the momentum equations ($+2\Omega\overline{\rho v}\sin\phi$, $-2\Omega\overline{\rho w}\cos\phi$). Tsai et al. (2014) test the presence in their simulations, of the proposed pumping mechanism of Showman & Polvani (2011) in the context of the $\beta$-plane where the Coriolis parameter is assumed to vary linearly with latitude (as opposed to $\propto\cos\phi$). Tsai et al. (2014) present figures of the meridional and vertical eddy momentum transport (strictly the divergence of the eddy momentum fluxes), supporting an eddy–driven jet, although they do not present contributions from mean flow terms.

Figure 14 shows the divergence of the eddy momentum fluxes for the Std Full and Std RT simulations, as the left and right columns, respectively, after 1 000 days. The latitudinal gradient in the meridional eddy momentum flux, the vertical gradient in the vertical eddy momentum flux, and their sum, are shown as the top, middle and bottom rows, respectively. This figure, and subsequently similar figures are truncated to the upper atmosphere i.e. $z>2\times 10^{6}$ m to aid interpretation666Contributions from transient perturbations, within these instantaneous plots, from the deeper regions (which have much higher densities, and therefore, momentum fluxes) can significantly alter the range of the plotted data. As discussed in Section 3.2, the deep atmosphere is still evolving there is not yet a mean flow, meaning the analysis we are performing is not meaningful in this region..

Figure 14 reveals some similarities in the momentum transport due to eddies, in our Std Full simulation, to those presented in Tsai et al. (2014). These structures are broadly consistent with the mechanism proposed by Showman & Polvani (2011), as the eddy momentum transport largely acts to bring prograde momentum into the jet region, however at around $4\times 10^{6}$ m, over the equator, a significant transport of retrograde momentum is present. The vertical term (middle left panel), shows predominantly prograde momentum being transport from lower altitudes (from below $5\times 10^{6}$ m) to the upper atmosphere and retrograde momentum at lower altitudes. The sum of these two terms, of the Std Full simulation, after 1 000 days, does not demonstrate clear balance between them (bottom left panel) and residual eddy momentum transport is present of the same order of the terms themselves. These patterns are steady in time, albeit gradually decreasing in magnitude (see Figure 21 in Appendix D). Interestingly, eddy momentum transport after 1 000 days is very different in the Std RT case (right column of Figure 14), compared to the Std Full simulation. The balance between the meridional and vertical component is closer, i.e. the magnitude of the summed momentum transports (bottom right) is much lower than the Std Full counterpart. The vertical term is the opposite to that shown for the Std Full simulation, as prograde momentum is transported into regions below an altitude of $4\times 10^{6}$ m. Additionally, the meridional component is transporting prograde eddy momentum in the flanks of the jet ($\phi=\pm 20^{\degree}$).

Figure 15 presents the mean flow contributions for the Std Full and Std RT simulations, as the left and right columns, respectively, after 1 000 days. The top, middle and bottom rows show the mean meridional ($-\frac{(\overline{\rho v}\,\overline{u}\cos^{2}\phi)_{,\phi}}{r\cos^{2}\phi}$) term, mean vertical ($-\frac{(\overline{\rho w}\,\overline{u}r^{3})_{,r}}{r^{3}}$) term, and their sum, respectively.

Figure 15 shows that the mean flow contributions are somewhat balanced throughout much of the atmosphere, but show a localised peak contribution in the sum over the equator at $4\times 10^{6}$ m, at $\sim 2\times 10^{6}$ m for the Std Full and Std RT simulations, respectively. These contributions are larger in magnitude than those provided by the divergence of the eddy momentum fluxes (compare with Figure 14), indicating a jet maintained by a mix of eddy momentum transport and the mean flow. The contributions from the remaining terms, deriving from the Coriolis terms of the momentum equations ($2\Omega\overline{\rho v}\sin\phi$ and $-2\Omega\overline{\rho w}\cos\phi$) contribute an order of magnitude lower momentum transport within the atmosphere (see Appendix D). Therefore, our results suggest that the jet driving mechanism, in our simulations, may actually be a combination of eddy and mean flow momentum transport.

4 Conclusions

In this study we have presented a set of simulations, nominally based on HD 209458b, aimed at investigating the robustness of the atmospheric dynamical morphology. These simulations have been designed to explore the dependence of the atmospheric flows on various assumptions made in standard models, primarily the state of the deeper, higher pressure atmosphere which is inaccessible to observations. We have presented the atmospheric dynamical structure found in this simulation set and for a subset, explored the eddy-mean interaction revealing momentum transport in the atmospheres.

Our results have shown that the super-rotating equatorial jet, found as a solution for hot Jupiter atmospheres from several models (see for example Cooper & Showman, 2005; Menou & Rauscher, 2009; Rauscher & Menou, 2010; Heng et al., 2011; Dobbs-Dixon & Agol, 2013; Parmentier et al., 2013; Showman et al., 2015; Helling et al., 2016; Kataria et al., 2016; Lee et al., 2016), is robust in our setup. Our model does not include a drag or friction at the bottom boundary potentially responsible for removing some of the dependence on initial conditions (see Liu & Showman, 2013; Cho et al., 2015). Only when artificially forcing the deep atmosphere over long timescales do we find the jet is significantly slowed, and reduced in latitudinal breadth. However, caution must still be exercised when interpreting results of simulations of hot Jupiters as, although the quantitative dynamical structures across our simulation set are consistent (apart from extreme cases where we strongly force the deep atmosphere over long timescales), the quantitative results are not. This variation in the dynamical structure, between our simulations (and indeed over time during the simulation) will in turn alter the thermodynamic profile of the atmosphere (this latter point warrants a more detailed follow-up study which we are engaged in, comparing results from the SPARC/MITgcm and the UM).

For simulations adopting a simplified radiative transfer scheme, our results show attributes of the proposed mechanism of Showman & Polvani (2011), however, the structure of the eddy momentum transport is complex, and there is a non-negligible contribution from the mean flow terms. Although the broad jet structure is similar between the Full and the RT simulations, we find very different patterns for the eddy momentum transport between these two simulations, suggesting indeed that the eddy momentum transport may not be the unique (or even the main) driver for the jet. Essentially, we find evidence for a jet driven by both contributions from eddies and mean flow interaction. A full investigation of the nature of the jet acceleration is beyond the scope of this work. However, we are performing a follow up study aimed at exploring the full, 3D solution to the standing linear problem studied by Showman & Polvani (2011) and Tsai et al. (2014) including an accurate representation of the heating rate (Debras et al., in prep).

Acknowledgements.

We thank the anonymous referee for comments which improved this manuscript. NM acknowledges funding from the Leverhulme Trust via a Research Project Grant. JM and CS acknowledge the support of a Met Office Academic Partnership secondment. DSA acknowledges support from the NASA Astrobiology Program through the Nexus for Exoplanet System Science. This work is partly supported by the European Research Council under the European Community’s Seventh Framework Programme (FP7/2007-2013 Grant Agreement No. 247060-PEPS and grant No. 320478-TOFU). The calculations for this paper were performed on the University of Exeter Supercomputer, a DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS, and the University of Exeter.

Appendix A Atmospheric flow structure

Figures 16 and 17 show the same information as Figures 3 and 4, shown in Section 3.1 but for the lower pressure surfaces at 213 and 21 600 Pa (as the top and bottom rows, respectively), after 1 200 and 10 000 days, respectively. These Figures reveal a broadly similar flow structure between the two times, and regardless of the completeness of the dynamical equations. This is to be expected given the strong radiative forcing at these pressures.

Figure 18 presents the results for the Std RT simulation, matching those presented in Figure 16, but after 1 600 days (the end of this simulation). Revealing differences similar to those previously discussed in Showman et al. (2009) and Amundsen et al. (2016).

Figure 19 presents the results for the $\Delta T_{\rm eq\rightarrow pole}$ simulation, matching those presented in Figure 17. These slices reveal several differences when compared with those presented in Figure 17. The lowest pressure slices, 213 Pa, show the diverging flow centred closer to $180^{\circ}$, and a sharper day–night temperature structure in the $\Delta T_{\rm eq\rightarrow pole}$, as opposed to the Std Full simulation. For the 21 600 Pa surface, the equatorial zonal flow is not uniformly prograde for the $\Delta T_{\rm eq\rightarrow pole}$ simulation, as also seen in the zonal and temporal mean of the zonal wind, Figure 7.

Appendix B Shallow-hot Jupiter

Mayne et al. (2014a) performed the Shallow-hot Jupiter (SHJ) test of Menou & Rauscher (2009), and compared to the results of Menou & Rauscher (2009) and Heng et al. (2011). Our simulations largely matched previous works except that the super rotating equatorial jet did not intersect the bottom, or inner, boundary. We assumed this to be caused by the difference in using a pressure-based (as is the case for Menou & Rauscher, 2009; Heng et al., 2011), or height-based model (the UM). However, the SHJ simulation has now been performed by Mendonça et al. (2016) using a height-based approach who were able to match the previous simulations, and obtain a jet intersecting the high pressure boundary. Therefore, we have revisited the SHJ test, and investigated this discrepancy further. Figure 20 shows the zonally and temporally averaged, from 200-1 200 days, zonal wind, in the same format as Figure 1 but for the SHJ case (see Mayne et al., 2014a, for details of the setup etc.). The left and right panels of Figure 20 show simulations applying ‘stronger’ ($K_{\lambda}\sim 0.158$, $K^{\prime}_{\lambda}\lesssim 1.2\times 10^{-12}$) and ‘weaker’ ($K_{\lambda}\sim 2.6\times 10^{-03}$, $K^{\prime}_{\lambda}\lesssim 2.0\times 10^{-14}$) diffusion, respectively. See Section 2.1 for explanation of how the diffusion coefficient is set ($K_{\lambda}$) and subsequently applied ($K^{\prime}_{\lambda}$). Figure 20 shows that the super rotating equatorial jet intersects the bottom boundary when using the larger diffusion coefficient, and therefore ‘stronger’ diffusion. The jet speeds are significantly slower when applying ‘weaker’ diffusion. This is counter–intuitive, and demonstrates clearly that not only the maximum wind/jet velocity, but also the shape of the atmospheric flow can be altered by the choice of diffusion settings (see also discussion in Heng et al., 2011; Li & Goodman, 2010). We are currently exploring this effect further, which is likely caused by the diffusion altering the eddy transport in the atmosphere.

Appendix C Derivation of the eddy mean interaction equation

Starting with the ‘full’ or ‘deep’ equations as defined in (see Mayne et al., 2014a). The derivation for the eddy–mean interaction equation, in height coordinates, is as follows: the zonal momentum equation is

[TABLE]

where $\frac{D}{Dt}$ is the material derivative. The mass continuity equation is

[TABLE]

Multiplying the momentum equation (Eqn 4) by $\rho$, and adding to the mass continuity equation (Eqn 6) multiplied by $u$, to the right hand side (which is equal to zero), gives

[TABLE]

Using the product rule,

[TABLE]

Therefore, Eqn 7 becomes,

[TABLE]

As,

[TABLE]

Eqn 9 can be expanded as,

[TABLE]

i.e. $X=\rho u$, and $\textbf{Y}=\textbf{u}=\left(u,v,w\right)$. Using the momentum equation (Eqn 4),

[TABLE]

i.e. multiply Eqn 4 by $\rho$. Therefore, combining Eqns 11 and 12, to eliminate $\rho\frac{Du}{Dt}$ yields,

[TABLE]

Rearranging by collecting all terms involving $v$ and $w$ on one side yields,

[TABLE]

Now we combine two of the $uw$ and $uv$ terms. Starting with the $uv$, expand the $\tan\phi$,

[TABLE]

rearranging and then multiplying the fractions by $\frac{\cos\phi}{\cos\phi}$, we can transform the $uw$ terms from Eqn 15 to

[TABLE]

This can then be compressed using the inverse product rule, i.e. $\partial\left(\left(\rho uv\cos\phi\right)\left(\cos\phi\right)\right)$ to yield,

[TABLE]

For the $uw$ terms from Eqn 14, multiplying the first by $\frac{r}{r}$ and the second by $\frac{r^{2}}{r^{2}}\frac{\partial r}{\partial r}$ yields

[TABLE]

which again can be compressed using the inverse product rule, $\partial\left(\left(r^{2}\rho wu\right)\left(r\right)\right)$,

[TABLE]

Now, Eqns 17 and 19 can be substituted into Eqn 14 to yield,

[TABLE]

We can now apply a zonal average to this equation, meaning that $\frac{\partial}{\partial\lambda}$ terms become zero. Critically, this only applies if a zonal average is performed along surfaces of equal $r$, as the derivatives are performed at constant $r$ (see discussion in Section 3.3),

[TABLE]

Then the wind terms involving the zonal wind, $u$ only, can be decomposed into a mean and perturbation using $u=\overline{u}+u^{\prime}$, $\overline{\rho u}=\overline{\rho}\,\overline{u}+\overline{\rho^{\prime}u^{\prime}}$, $\overline{\rho uv}=\overline{\rho v}\,\overline{u}+\overline{\left(\rho v\right)^{\prime}u^{\prime}}$ and $\overline{\rho uw}=\overline{\rho w}\,\overline{u}+\overline{\left(\rho w\right)^{\prime}u^{\prime}}$, yielding,

[TABLE]

Finally, rearranging this equation so that the mean terms are on the left and the perturbation terms on the right yields

[TABLE]

which matches that given as Eqn (6) in Hardiman et al. (2010).

Next in the ‘shallow’ or ‘primitive’ cases (see Mayne et al., 2014a).

The zonal momentum equation is

[TABLE]

where $r\rightarrow R_{\rm p}$, $\frac{uw}{r}\rightarrow 0$ and $w\cos\phi\rightarrow 0$. Similarly, Eqn 10 becomes,

[TABLE]

where $r\rightarrow z$ in the derivative.

Therefore, the key steps of combining the $uv$ and $uw$ terms, using Eqns 17 and 19 are simplified greatly. For the $uw$ there is only one term (as $\frac{uw}{r}$ has been removed), and for the $uv$ terms the process is exactly the same, just replacing $r$ with $R_{\rm p}$. Therefore, this yields,

[TABLE]

which can clearly be simplified further, but we leave in this form to provide comparison with Eqn 23.

Appendix D Eddy mean interaction terms

In Section 3.3.2, Figure 14 the latitudinal gradient of the meridional and the vertical gradient of the vertical eddy momentum fluxes, are presented after 1 000 days for the Std Full and Std RT simulations. Figure 21 presents the same data as the left column of Figure 14, i.e. the Std Full simulation, but after 10 000 days. Here, the lowest plotted height has been increased, to avoid the presence of fluctuations in the deep atmosphere dominating the plot (as discussed in Section 3.3.2). The structures are very similar, to the earlier time, albeit slight reduced in magnitude.

The remaining terms in the eddy-mean interaction equation, aside from those in Figures 14 and 15 in Section 3.3.2, are the mean flow Coriolis terms (see Section 3.3.2). The mean Coriolis terms, in the meridional ($+2\Omega\overline{\rho v}\sin\phi$) and vertical ($-2\Omega\overline{\rho w}\cos\phi$), as well as their sum, are shown as the top, middle and bottom rows of Figure 22 in the same format as Figure 14. The Std Full and Std RT simulations are shown as the left and right columns, respectively, after 1 000 days. Figure 22 shows that the remaining Coriolis terms have a significantly reduced contribution over the eddy and mean flow momentum flux transport, in the region of the jet (i.e. aside from fluctuations in the low altitude atmosphere).