EDP Sciences
Free Access
Issue
A&A
Volume 595, November 2016
Article Number A36
Number of page(s) 11
Section Planets and planetary systems
DOI https://doi.org/10.1051/0004-6361/201629183
Published online 26 October 2016

© ESO, 2016

1. Introduction

Hot Jupiters, Jupiter-sized planets orbiting close to their parent stars, have the most observationally constrained atmospheres of all exoplanets. Transmission spectroscopy has been used to detect sodium, potassium and water (Charbonneau et al. 2002; Snellen et al. 2008; Redfield et al. 2008; Sing et al. 2011, 2012, 2015; Deming et al. 2013; Wakeford et al. 2013; McCullough et al. 2014; Evans et al. 2016) as well as identifying a continuum of hot Jupiter atmospheres ranging from “cloudy” to “clear” (Sing et al. 2016). Temperature contrasts and brightness temperature maps have been derived from phase curves (Knutson et al. 2007, 2009, 2012; Maxted et al. 2013; Zellem et al. 2014; Stevenson et al. 2014), and finally wind velocities have been estimated (Snellen et al. 2010; Louden & Wheatley 2015). These measurements have revealed that many hot Jupiters are subject to significant heat redistribution between their day and night sides, with hot spots shifted eastward of the substellar point, which can not be modelled consistently with one-dimensional (1D) models.

This has motivated the adaptation of global circulation models (GCMs), which to date have been applied to study the atmospheric circulation of several hot Jupiters (Showman & Guillot 2002; Cooper & Showman 2005; Showman et al. 2009; Lewis et al. 2010; Rauscher & Menou 2010; Thrastarson & Cho 2010; Polichtchouk & Cho 2012; Dobbs-Dixon & Agol 2013; Kataria et al. 2013, 2015, 2016). More recently, GCMs have also been used to study some aspects of cloud formation and evolution in hot Jupiter atmospheres (Parmentier et al. 2016; Helling et al. 2016; Lee et al. 2016). GCMs are three-dimensional (3D) models and include a dynamical core solving the equations of motion for the fluid, combined with a radiation scheme for treating stellar heating and thermal cooling of the atmosphere. Many GCMs applied to hot Jupiters solve the primitive equations (see e.g. Showman & Guillot 2002; Showman et al. 2009; Rauscher & Menou 2010; Kataria et al. 2013, 2015), which are an approximation to the full Navier-Stokes equations assuming that the atmosphere is shallow compared to the radius of the planet, in hydrostatic equilibrium and has a gravity constant with height. The exceptions are Dobbs-Dixon & Agol (2013) and Mayne et al. (2014a), who solved the full Navier-Stokes equations. As the vertical extent of hot Jupiter atmospheres can be about 10 % of the planet radius, the validity of the primitive equations is questionable (Mayne et al. 2014a).

Radiation schemes in initial hot Jupiter GCMs employed Newtonian forcing, where the temperature is relaxed linearly towards equilibrium PT profiles (Showman & Guillot 2002; Cooper & Showman 2005, 2006; Showman et al. 2008; Rauscher & Menou 2010). Such approaches have many disadvantages as radiative heating and cooling are not treated self-consistently: (i) appropriate equilibrium profiles are difficult to obtain from 1D models; (ii) the temperature relaxation is linear while in reality it may be non-linear for large deviations from the equilibrium profiles; (iii) atmospheric interactions due to exchange of radiative energy such as emission and absorption of thermal radiation are ignored and (iv) the model flexibility is poor since for each new planet modelled, the forcing must be changed.

More recent hot Jupiter GCMs have adopted radiation schemes using the two-stream approximation with grey (Rauscher & Menou 2012) or average opacity schemes (Dobbs-Dixon & Agol 2013), which has recently been shown to yield inaccurate heating rates when considering molecular absorption in these atmospheres (Amundsen et al. 2014). The most sophisticated radiation scheme employed to date is that presented in Showman et al. (2009) and later used in Lewis et al. (2010) and Kataria et al. (2013, 2015, 2016), which adopts the two-stream approximation (Thomas & Stamnes 2002) for the stellar component and the two-stream source function technique (Toon et al. 1989) for the thermal component, combined with the correlated-k method (Lacis & Oinas 1991) for treating opacities. This has been shown to yield significantly better agreement with observations compared to using Newtonian forcing (Showman et al. 2008, 2009). Amundsen et al. (2014) showed that the two-stream approximation and correlated-k method give accurate fluxes and heating rates for hot Jupiter atmospheres.

We have adapted the UK Met Office GCM, the Unified Model (UM), for the study of hot Jupiters. The UM dynamical core solves the full 3D Navier-Stokes equations with a height-varying gravity, and the radiation scheme is state-of-the-art using the two-stream approximation and correlated-k method to treat opacities. The adaptation of the dynamical core (Mayne et al. 2014a,b) and radiation scheme (Amundsen et al. 2014) have been presented in previous publications. Preliminary results from our model have been presented in Helling et al. (2016), however, here we use an updated opacity database and present the model and results in more detail.

The goal of the present work is to (i) provide the technical details of the coupling between the UM dynamical core and adapted radiation scheme for hot Jupiters for future reference; (ii) provide the first comparison between two hot Jupiter GCMs (ours and that of Showman et al. 2009) with similar state-of-the-art radiation schemes and investigate the robustness of these GCMs; (iii) evaluate differences in resulting synthetic observations calculated from model output as the model of Showman et al. (2009) has already been used extensively in the literature (see e.g. Agúndez et al. 2014; Fortney et al. 2010; Kataria et al. 2015; Moses et al. 2011; Wong et al. 2016; Zellem et al. 2014); and (iv) provide a guide for future more in-depth intercomparisons of these models.

To ease comparison with earlier GCMs and the models of Showman et al. (2009), we use parameters similar to those of HD 209458b. This is the first attempt to compare results obtained with two different hot Jupiter GCMs with sophisticated radiation schemes. Both dynamical cores1 and radiation schemes (Ellingson et al. 1991; Collins et al. 2006; Oreopoulos et al. 2012) are tested thoroughly through intercomparison projects for Earth-like conditions. Intercomparison of hot Jupiter GCMs will in the future become crucially important as the quality of observations improve.

In contrast to Showman et al. (2009), we do not include TiO and VO in our model. Unfortunately the model becomes unstable due to the intense heating of the upper atmosphere caused by these molecules, a consequence of their large opacity at visible wavelengths. This prevents a detailed comparison to the model of HD 209458b in Showman et al. (2009), however, at present there is no evidence for TiO and VO, or even a temperature inversion, in the atmosphere of this planet (Diamond-Lowe et al. 2014; Hoeijmakers et al. 2015; Evans et al. 2015; Schwarz et al. 2015; Line et al. 2016).

This paper is organised as follows: in Sect. 2 we briefly describe the model, including the dynamical core and radiation scheme. In Sect. 3 we discuss results from running the model of Mayne et al. (2014a) which uses Newtonian forcing in place of an accurate radiation scheme before discussing results from our coupled model with sophisticated radiative transfer in Sect. 4. Synthetic observations are calculated for all models presented and compared to available observations of HD 209458b in the literature. Our conclusions are presented in Sect. 5.

2. Model description

In this section we briefly describe the dynamical core (Sect. 2.1) and radiation scheme (Sect. 2.2), but refer to Mayne et al. (2014a), Amundsen et al. (2014), Amundsen (2015) and Amundsen et al. (2016) for more details.

2.1. Dynamics

We used the Met Office UM with the Even Newer Dynamics for General Atmospheric Modelling of the Environment (ENDGAME) dynamical core (Wood et al. 2014) to solve the non-hydrostatic, deep-atmosphere Navier-Stokes equations for planetary atmospheres with a height-varying gravity. The equations are solved on a latitude-longitude-height grid using a semi-Lagrangian semi-implicit scheme. We use free-slip impermeable upper and lower boundaries located at a fixed height. We have previously applied ENDGAME to HD 209458b successfully using a Newtonian forcing scheme (Mayne et al. 2014a), and we refer to this paper and references therein for more details on the dynamical core. Here we exclusively solve the “full” equation set solving the non-hydrostatic deep-atmosphere equations with a height-varying gravity.

The diffusion scheme is described in Mayne et al. (2014a,b), and includes separate components in the longitude (Kλ) and latitude (Kφ) directions. The UM numerical scheme does not explicitly enforce axial angular momentum conservation (AAM), however, in the results presented here AAM is conserved to better than 98%. AAM conservation was maximised by using Kλ ~ 0.16 and Kφ = 0.0, which is applied in all simulations presented in this work and those presented in Mayne et al. (2014a), where the latter incorrectly reported the values of these constants. As described in Mayne et al. (2014a,b), due to the particular horizontal grid staggering adopted at the pole, we do not need to use a polar filter, but note that the diffusion scheme has some aspects in common with a polar filter.

2.2. Radiation scheme

To calculate radiative heating rates we used the Suite of Community Radiative Transfer codes based on Edwards and Slingo (SOCRATES) scheme2 (Edwards & Slingo 1996; Edwards 1996), which uses the two-stream approximation combined with the correlated-k method for both the stellar and thermal components. We have presented the adaptation and testing of this radiation scheme for hot Jupiter-like atmospheres in Amundsen et al. (2014), where it was found to yield accurate fluxes and heating rates by comparing to discrete ordinate line-by-line calculations. Here we review modifications made to the radiation scheme since it was presented in Amundsen et al. (2014).

2.2.1. The radiative transfer equation

For the thermal component of the radiation we solve the two-stream equations as formulated by Zdunkowski & Korb (1985) and Edwards (1996) with no scattering and a diffusivity D = 1.66, which was found to yield the most accurate fluxes and heating rates in Amundsen et al. (2014). For the stellar component we solve the two-stream equations as formulated by Zdunkowski et al. (1980), which uses D = 2. Rayleigh scattering by H2 and He is included, with refractive indices for H2 and He from Leonard (1974) and Mansfield & Peck (1969), respectively3, and an anisotropy factor ρn = 0.02 (Penndorf 1957). Refractive indices for H2 and He are combined using the Lorentz-Lorentz relation (Heller 1965).

Vertically, the dynamical core defines potential temperatures in layers and exner pressures on levels (layer interfaces). The radiation scheme uses the same grid, with fluxes calculated at the levels defined by the exner pressure, and layer properties are set based on the potential temperature. As temperatures are also required at the exner pressure levels by the radiation scheme they are interpolated linearly in height from the layer values. Layer radiative heating rates can then be calculated and applied directly by differencing fluxes at neighbouring levels without interpolation.

2.2.2. Opacities

Our opacity database includes absorption by H2O, CO, CH4, NH3, and H2–H2 and H2–He collision induced absorption (CIA) as described in Amundsen et al. (2014), using the newest ExoMol line lists (Tennyson & Yurchenko 2012; Tennyson et al. 2016) where available, with the recent addition of the alkali metals Li, Na, K, Rb and Cs. Alkali metal opacities are included using transition probabilities and broadening coefficients from VALD3 (Heiter et al. 2008). Voigt profiles are used for all lines with a line cut-off at 4000 cm-1 from the line centres except for the Na and K D lines for which we use line profiles from the PHOENIX atmosphere code (Allard et al. 1999, 2003, 2007, Derek Homeier, priv. comm.). We have updated our CH4 opacities to use the YT10to10 CH4 line list (Yurchenko & Tennyson 2014).

Opacities are treated using the correlated-k method (Lacis & Oinas 1991) as described in Amundsen et al. (2014), but k-coefficients are here computed individually for each gas. Our 32 bands are defined in Amundsen et al. (2014), and in each band the main absorber is found by comparing transmissions using the maximum equilibrium abundance for each gas. k-coefficients for individual gases are combined on-the-fly in the UM using equivalent extinction (Edwards 1996), where all gases except the strongest absorber in each band is taken into account through a grey absorption. The direct stellar component, however, is computed directly by multiplying transmissions for each gas, which assumes absorption lines for different gases are randomly overlapping (Lacis & Oinas 1991), using all k-coefficients for all gases. We have found this approach to be more accurate than using a pre-computed table of k-coefficients for the gas mixture as used in Amundsen et al. (2014) and Showman et al. (2009) as the use of such tables involves interpolating both mixing ratios and gas opacities in temperature and pressure, not only the opacities of individual absorbers. For a more detailed discussion and a comparison of different treatments of overlapping gaseous absorption, see (Amundsen et al. 2016).

2.2.3. Abundances

Abundances were calculated as in Amundsen et al. (2014) using the analytical chemical equilibrium abundance formulas for H2O, CO, CH4 and NH3 from Burrows & Sharp (1999). The alkali metal abundances are approximated by assuming they are in atomic gaseous form above the chemical transformation temperature, for alkali metal i, and that for their atomic gaseous abundance is negligible. We take the chemical transformation curves for the alkali metal chlorides from Burrows & Sharp (1999) and apply an additional smoothing of the form (1)where φi(T) is the normalised abundance, is the chemical transformation temperature and is the characteristic scale over which the abundance changes for species i. We adopt 20 K for alkali metals. Physically this is a primitive way of taking into account the transition between for example Na and NaCl and avoids numerical problems associated with non-continuous abundance changes. We note, however, that the particular functional form of this smoothing has no physical basis, but was chosen as it is symmetric about , and it and all its derivatives are continuous functions of temperature.

2.2.4. Boundary conditions

An extra layer is included in the radiation scheme to account for absorption, emission and scattering above the dynamically modelled domain. The layer extends up to zero pressure, and extensive testing has shown that the absorption is accurately taken into account with temperatures extrapolated linearly in log pressure and the temperature at zero pressure set to the smallest temperature in our PT grid (70 K). At the lower boundary we impose a net intrinsic flux , thereby taking into account heat escaping from the planet interior. The thermal upward flux at the lower boundary surface, , is then given by (2)where is the downward flux at the lower boundary. To ease implementation we use the value of in Eq. (2) from the previous radiation time step. Consequently, the value of Tsurf, and therefore the upward surface flux used in the lower boundary condition, will lag one radiative time step behind the radiative transfer calculation. We have found the temporal variations in to be very small compared to the radiative time step, which ensures the validity of this approximation. More details on the boundary conditions are presented in Amundsen (2015).

2.3. Synthetic observations

We have calculated synthetic observations from UM output using our 1D discrete ordinate radiation code ATMO (Amundsen et al. 2014; Tremblin et al. 2015, 2016; Drummond et al. 2016). It uses the same opacity sources as described in Sect. 2.2.2, and to compute synthetic observations we use high resolution k-tables with 5000 bands with band limits evenly spaced at 10 cm-1 intervals. We note that the use of high resolution k-tables is necessary as a line-by-line approach is too computationally expensive and a reduced line-by-line resolution of ~1 cm-1 yields very large errors in band-integrated fluxes. Chemical equilibrium abundances are calculated using a Gibbs energy minimisation scheme (Drummond et al. 2016). Our calculations of emission spectra and phase curves from UM output are detailed below.

2.3.1. Emission spectra

The emission from a planet as measured on Earth is given by (Seager 2010) (3)where Rp,TOA is the planet radius at the top of the atmosphere, Do is the distance to the observer, and Is(θ,φ,ϑo,ϕo) is the intensity at the top of the atmosphere at the location defined by the polar angle θ and azimuth angle φ, which can be directly related to the latitude and longitude, in the direction of the observer (ϑo,ϕo), where ϑo is the polar angle and ϕo is the azimuth angle. The definitions of these angles are illustrated in Fig. 1. The coordinate system in which both the location (θ,φ) and direction (ϑ,ϕ) of the radiation are defined is placed so that the z-axis always points towards the observer, that is, ϑo = 0 and the angle between the planet surface normal at (θ,φ) and the direction of the observer is θ.

thumbnail Fig. 1

Illustration of the definition of the angles used to calculate the hemispherically integrated emission, see Eq. (3). (θ,φ) and (ϑ,ϕ) denote the position and the direction of the radiation, respectively, with the coordinate system placed such that the z-axis points towards the observer.

Open with DEXTER

The intensity Is is calculated at 16 discrete angles determined by the Gauss-Legendre points of the discrete ordinate method for all atmospheric columns from the UM and interpolated to obtain the intensity in the direction of the observer. These intensities, Is(θ,φ,ϑo,ϕo), are then integrated according to Eq. (3) to obtain the observed emitted flux.

2.3.2. Phase curves

Phase curves are the emission from the planet, as viewed from Earth, as a function of time or orbital phase angle. The integrated emission as a function of orbital phase is given by Eq. (3) for different observer directions, which are given by the phase angle α ∈ [0°,360°), where α = 0 is primary eclipse and α = 180° is secondary eclipse. Assuming the planet is tidally locked and in a steady state the intensity at the top of the atmosphere for a given latitude and longitude will be constant as a function of time. This simplification enables us to calculate Is(θ,φ,ϑ,ϕ) only once for the entire phase curve, greatly decreasing the computation time. As in Showman et al. (2009) we ignore the small inclination of the orbit (Fortney et al. 2006).

2.4. Model setup and parameters

Our model setup is similar to that used in Mayne et al. (2014a) with a few modifications. We provide our adopted model parameters in Table 1. Compared to Mayne et al. (2014a) the planet radius at the lower boundary of the model, Rp, has been decreased to account for the vertical extent of the atmosphere, the specific heat capacity, cP, has been changed to be in agreement with Showman et al. (2009), the specific gas constant, ℛ = R/ barm, has been changed to account for the mean molecular weight barm =2.3376 g mol− 1 used in Amundsen et al. (2014), the pressure at the lower boundary, Pbottom, has been chosen to be in agreement with Showman et al. (2009) and the height of the upper boundary has been slightly adjusted to account for the modified atmospheric scale height. We find that, for both numerical stability and accuracy, we need to use a dynamical time step of 30 s, much smaller than the time step used in Mayne et al. (2014a), but in agreement with time steps used by Showman et al. (2009).

It is common practice in GCMs to call the radiation scheme, that is, update fluxes and heating rates, less frequently than every dynamical time step. This is done mainly for computational efficiency, and is possible as changes in fluxes and heating rates may be small between dynamical time steps. We have tried several different radiative time steps, from calling the radiation scheme every dynamical time step to calling it every ten dynamical time steps, and have found that calling it every five dynamical time steps is a good compromise between numerical accuracy and computational cost. This leads to a radiative time step of 150 s.

We initialise our model with a PT profile from our radiative-convective equilibrium code ATMO (Amundsen et al. 2014; Tremblin et al. 2015, 2016). Profiles are calculated using μ0 = cosθ0 = 0.5, where θ0 is the star zenith angle, which corresponds to a day-side average and reduces the stellar flux at the top of the atmosphere by a factor of 1/2. To obtain a globally averaged PT profile the top of the atmosphere flux is further reduced by a factor 1/2 to account for redistribution to the night-side. All models are initialised with zero winds.

3. Model with Newtonian forcing from Mayne et al. (2014a)

Before discussing results from the coupled model we briefly summarise the results from running the model of Mayne et al. (2014a), which uses the Newtonian forcing scheme described in Cooper & Showman (2005, 2006), Rauscher & Menou (2010) and Heng et al. (2011). This enables us to compare the UM and SPARC/MITgcm without the additional complication of the radiation schemes used. The model setup used here is identical to Mayne et al. (2014a), with equilibrium PT profiles and timescales are from Iro et al. (2005). We initialise the model using an average between the day and night side PT profiles with zero winds. In Mayne et al. (2014a) results are averaged temporally from 200 d to1200 d (d denotes Earth days) as prescribed by the Heng et al. (2011) benchmarks. We have run the model for >1000 d, and, in contrast to Mayne et al. (2014a) who computed temporal averages from 200 d to1200 d, we generally show model results after 1600 d with no temporal averaging as is most common in studies applying GCMs to hot Jupiters (Showman et al. 2008, 2009; Kataria et al. 2013, 2015). We observe temperatures to reach an approximate steady-state after 1000 d for P<105 Pa, which is expected to be the observable part of the atmosphere (see Sect. 3.2). Longer simulation times will be needed to study the deeper P>105 Pa regions. This is in agreement with previously published hot Juptier GCMs with simplified forcing (Showman et al. 2008).

Table 1

Model parameters adopted for HD 209458b.

3.1. Results

We show in Fig. 2 (left column) the temperature and horizontal wind at various atmospheric depths as a function of longitude and latitude after 1600 d. At 100 Pa winds diverge from the hotspot located at substellar point (180° longitude, latitude). It is worth noting that, due to the very small radiative timescale at 100 Pa, the temperature is almost identical to the equilibrium temperature, which causes a large temperature contrast (>900 K) between the day and night side. For increasing pressures, the dynamical regime is dominated by a super-rotating equatorial jet spanning all longitudes. Dynamical processes redistributing the heat away from the substellar point become more dominant, which causes the temperature difference between the day and night side to decrease.

thumbnail Fig. 2

Horizontal wind velocity as arrows and temperature [K] as colours and contours from our models of HD 209458b after 1600 d. Left column: results from the model with Newtonian forcing discussed in Sect. 3 at 100 Pa, 3 × 103 Pa, 3 × 104 Pa, 1 × 105 Pa (from top to bottom, from Mayne et al. 2014a). Right column: results from the coupled model discussed in Sect. 4 at 3 Pa, 3 × 103 Pa, 3 × 104 Pa, 1 × 105 Pa (from top to bottom).

Open with DEXTER

We show in the left-hand panel of Fig. 3 the zonal mean of the zonal wind as a function of pressure and latitude. The zonal jet in the eastward direction mentioned above is clearly seen, and it reaches its maximum strength at about 103 Pa with a velocity of about 7 km s-1. At higher latitudes the mean flow is in the opposite (westward) direction, and much weaker in amplitude, with a maximum of about 1.2 km s-1.

thumbnail Fig. 3

Zonal mean of the zonal wind velocity [m s− 1] after 1600 d for the models of HD 209458b. Left: model with Newtonian forcing from Mayne et al. (2014a) discussed in Sect. 3. Right: coupled model discussed in Sect. 4. Red indicates a prograde wind, blue indicates a retrograde wind.

Open with DEXTER

In Fig. 4 we plot PT profiles for several different latitudes and longitudes. The temperature varies significantly across the globe, with night side temperatures down to ~600 K and day side temperatures up to ~1500 K at 103 Pa. A dynamically induced temperature inversion is even seen on the day side of the planet, which is caused by strong heating at the top of the atmosphere due to the short radiative timescale and the equatorial jet bringing cold material from the night side to the day side cooling the atmosphere down at larger pressures.

3.2. Discussion

These results are similar to those obtained with other hot Jupiter GCMs using Newtonian forcing schemes (see e.g. Cooper & Showman 2005; Showman et al. 2008; Heng et al. 2011; Rauscher & Menou 2010), as discussed in Mayne et al. (2014a). Using ATMO we have calculated synthetic day-side emission spectra and phase curves, which are shown in Fig. 5 and 6, respectively, together with observational data points from the literature.

The day-side emission spectrum agrees reasonably well with observations. This is rather surprising as the forcing profiles were estimated from the globally averaged PT profile of Iro et al. (2005; Cooper & Showman 2005; Rauscher & Menou 2010; Heng et al. 2011), and are therefore not expected to be particularly accurate. The amplitudes of the 4.5 μm phase curve show reasonably good agreement with the observed phase curve (Zellem et al. 2014), but the significant phase offset in the observed phase curve is lacking. The latter could be due to an underestimate of the radiative timescale, which would lead to a smaller offset of the hottest point in the atmosphere from the substellar point.

Fortney et al. (2006) presented synthetic day side emission spectra and phase curves using results from the GCM presented in Cooper & Showman (2005, 2006), which uses the same forcing scheme as the one adopted here (Mayne et al. 2014a). We have plotted the synthetic observations from Fortney et al. (2006) in Figs. 5 and 6 as dashed lines to ease comparison with our synthetic observations. Our day-side emission is somewhat larger than that obtained by Fortney et al. (2006), while differences in night side fluxes are smaller. There is also a noticeable difference in the phase offsets of the peak flux from 180° between the models, we obtain a significantly smaller phase offset than Fortney et al. (2006).

thumbnail Fig. 4

PT profiles around the globe after 1600 d for the model of HD 209458b with Newtonian forcing from Mayne et al. (2014a). Red solid lines and blue dashed-dotted lines are day and night side profiles, respectively, at latitude. Magenta dashed lines and cyan dotted lines are profiles between and 90° latitude for longitudes 180° and , respectively.

Open with DEXTER

It is difficult to pinpoint the exact causes of these differences, but there are several factors that might contribute: slight discrepancies in the temperature may be caused by numerical details of the GCMs. The model of Cooper & Showman (2005, 2006) solve the primitive equations on a pressure based grid using a gravity constant with height, while we solve the full 3D Navier-Stokes equations with a height-varying gravity on a height-based grid. In addition the numerical schemes, grids and resolutions are different. In our model, however, solving the primitive equations and assuming a gravity constant with height only has a minor effect on the emission compared to solving the full 3D Navier-Stokes equations with a height-varying gravity.

The main differences in emission may therefore be caused by differences in the tools used for post-processing such as different line list and line width sources, and slightly different elemental abundances, resulting in the calculation of somewhat different opacities. The fact that GCMs with the same simplified forcing scheme give such different results emphasise the need to include the post-processing tools in intercomparison studies.

Before discussing results from the coupled model we briefly reiterate the disadvantages of using Newtonian forcing schemes, as in e.g. Cooper & Showman (2005), Heng et al. (2011), Mayne et al. (2014a), Menou & Rauscher (2009), Rauscher & Menou (2010), Showman & Guillot (2002) and Showman et al. (2008), to treat the radiation:

  • 1.

    Equilibrium PT profiles needed by the forcing schemes, which vary as a function of latitude and longitude, are difficult to obtain from 1D models.

  • 2.

    The temperature relaxation is linear while in reality it is non-linear for large deviations from the equilibrium profiles.

  • 3.

    Atmospheric interactions due to exchange of radiative energy such as emission and absorption of thermal radiation are ignored.

  • 4.

    The model flexibility is poor since for each new planet modelled the forcing prescription must be changed.

Global circulation models with Newtonian forcing schemes can still be useful for exploring underlying dynamical processes, see e.g. Showman & Polvani (2011) and Komacek & Showman (2016). Including a proper treatment of radiative heating and cooling are essential, however, in order to improve model flexibility and, as we demonstrate in the next section (and discuss in Sect. 1), improve agreement with observations.

thumbnail Fig. 5

Observed (points) and synthetic (lines) emission spectra for HD 209458b models with Newtonian cooling. The solid line has been calculated from the model presented in Mayne et al. (2014a) and Sect. 3 using ATMO, the dashed line is the synthetic emission spectrum from Fortney et al. (2006), which is based on the models of Cooper & Showman (2005, 2006) using the same Newtonian forcing scheme. The black points are observations from Swain et al. (2008) (●), Crossfield et al. (2012) (▼), Deming et al. (2005) (◄), Zellem et al. (2014) (►), Diamond-Lowe et al. (2014) (■) and Evans et al. (2015) (♦).

Open with DEXTER

thumbnail Fig. 6

Synthetic Spitzer IRAC phase curves from the Newtonian forcing model. The solid lines are calculated from the model presented in Mayne et al. (2014a) and Sect. 3 using ATMO, the dashed lines are the synthetic phase curves from Fortney et al. (2006), which are based on the models of Cooper & Showman (2005, 2006) using the same Newtonian forcing scheme. The models have been integrated over the IRAC bands using the filter functions. The data points are from Zellem et al. (2014) (►), Diamond-Lowe et al. (2014) (■) and Evans et al. (2015) (♦). The best fit to the observed 4.5 μm phase curve from Zellem et al. (2014) is shown as a solid black line, the grey shaded area is the 1σ uncertainty for the offset of the observed planet to star flux ratio.

Open with DEXTER

4. Full coupled model results

Here we present results from simulations using the UM incorporating the full radiation scheme discussed in Amundsen et al. (2014), including the modifications described in Sect. 2.2. In the literature it is usually not explicitly stated how long simulations have been run for (Showman et al. 2009; Kataria et al. 2013, 2014, 2015). We have run simulations for >1000 d. For the observable part of the atmosphere, that is, pressures 105 Pa, the atmosphere has approximately reached a steady-state, while at larger pressures the atmosphere is still evolving and much longer integration timescales would be needed to study the evolution and its ramifications.

4.1. Results

Horizontal wind velocities and temperatures are plotted in Fig. 2 (right column) after 1600 d at several different pressures, and can be compared to the left column obtained with the model with Newtonian forcing. General features are similar to those found in the model with Newtonian forcing. At low pressures (100 Pa) the flow is again diverging from the substellar point, with a hotspot shifted eastward from the substellar point.

In the right panel of Fig. 3 we show the zonal mean zonal wind velocity after 1600 d as a function of pressure and latitude. As in the model with Newtonian forcing the eastward equatorial jet is prominent, with a mean westward flow at higher latitudes.

In Fig. 7 we show the variation in PT profiles across the globe. A large variation is evident, and some night side profiles have higher temperatures than some day side profiles. This is due to the strong eastward advection causing the terminator at 270° longitude to be much warmer than that at 90° longitude. For pressures 105 Pa profiles at latitude are dominated by the equatorial jet, causing very small temperature variations as a function of longitude. At other latitudes, however, temperature variations are larger. The deep atmosphere between 105 and 107 Pa is generally much hotter than the initial PT profile, with temperatures approaching 2000 K. This region is, due to the long dynamical timescale, evolving very slowly, and has not yet converged to a steady-state (Mayne et al. 2014a). It is clear, however, that temperatures are slowly increasing across the globe at these deep levels.

thumbnail Fig. 7

PT profiles around the globe after 1600 d for the coupled model of HD 209458b discussed in Sect. 4. Red solid lines and blue dashed-dotted lines are day and night side profiles, respectively, at latitude. Magenta dashed lines and cyan dotted lines are profiles between and 90° latitude for longitudes 180° and , respectively. The black line is the initial PT profile adopted.

Open with DEXTER

We also run simulations using a hotter initial temperature pressure profile, increasing from the standard global average 1D profile, uniformly, by 400 K and 800 K. The results are shown in Figs. 7 and 8, where it becomes clear that the initial PT has not converged to a steady-state for P105 Pa. In fact, the atmosphere may be converging towards temperatures significantly hotter than estimated from a 1D global average at P ~106 Pa. Unfortunately, due to computational limitations, we are unable to run our model for significantly longer timescales. However, these initial results suggest further work is required with regard to the sensitivity of hot Jupiter GCM results to such changes in the initial profile.

thumbnail Fig. 8

PT profiles around the globe after 1600 d for the coupled model of HD 209458b initialised with a PT profile that is 400 K (left) and 800 K (right) hotter than the global 1D mean. Lines are as in Fig. 7.

Open with DEXTER

thumbnail Fig. 9

Temperature difference ΔT between the PT profiles in the right-hand panel of Fig. 8, which are from the model with a 800 K hotter initial condition compared to the global 1D mean, and the PT profiles in Fig. 7, which are from the model initialised with the global 1D mean PT profile. Temperature differences are small for P105 bar, while the models have clearly not reached a steady-state for P105 bar. Results are similar for the case with a 400 K hotter initial condition. Lines are as in Fig. 7.

Open with DEXTER

4.2. Discussion

Our results are in qualitative agreement with those of Showman et al. (2009). The model exhibits a strong eastward equatorial jet, and hotspot shifted eastward of the substellar point. Unfortunately SPARC/MITgcm results for HD 209458b without TiO and VO have not been published in detail, which prohibits a more detailed comparison of temperature and wind fields. One characteristic of the hot Jupiter simulations presented in Showman et al. (2009) is what the authors term a “vertical coherency” of temperatures. This is particularly apparent in the solar metallicity HD 187733b setup, which excludes TiO and VO opacities, and is the most well matched to our HD 209458b simulations in terms of opacities. Vertical coherency describes the fact that the horizontal position of the hottest and coldest parts of the atmosphere vary only modestly between 102 Pa and 105 Pa. This was not seen, nor expected, in previous simulations adopting Newtonian forcing, as the radiative timescale varies by about two orders of magnitude over these depths. Therefore, one might expect the balance between the radiative forcing and advection to change with depth and lead to a significant change in the horizontal temperature distribution. Showman et al. (2009) proposed that the observed vertical coherency was caused by the vertical interaction of thermal radiation reducing vertical temperature gradients. This effect is self-consistently included in the models of Showman et al. (2009), but not included in those adopting Newtonian forcing.

We do not see vertical coherency in our simulations, despite self-consistently treating the thermal radiation: the position of the hottest and coldest points vary significantly with pressure. This is particularly noticeable at 105 Pa = 1 bar where both winds and temperatures are dominated by the eastward equatorial jet, and a weak retrograde flow at higher latitudes. Our longitudinal temperature variations are very small (<100 K), in contrast to the models in Showman et al. (2009) where temperatures vary by up to 500 K at high latitudes. The reason for this discrepancy is unclear, but we have run our model significantly longer, giving the system time to equilibrate at higher pressures, and we do not assume the atmosphere to be shallow. This may help explain the weaker vertical coherency in our model, but more in-depth comparisons are needed to understand these differences in more detail.

The fact that the deep layers are heating up compared to the initial PT profile is intriguing. Even though the temperatures have not converged to a steady-state, it suggests that the deep atmosphere, in equilibrium, should be hotter than predicted by simple 1D models. Further analysis is required to understand this feature.

We show in Fig. 10 and 11 the synthetic day side emission spectrum and phase curves calculated using ATMO. The day side emission spectrum provides a reasonably good fit to the observations, and so does the offset of the peak emission from 180° in the 4.5 μm phase curve, while the offset of the minimum emission is larger than observed.

thumbnail Fig. 10

Same as Fig. 5, but the synthetic emission spectrum has been calculated using results from our coupled model (Sect. 4) with ATMO.

Open with DEXTER

thumbnail Fig. 11

Same as Fig. 6, but the synthetic emission spectrum has been calculated using results from our coupled model (Sect. 4) with ATMO. The dashed line is the model 4.5 μm phase curve from Zellem et al. (2014) obtained using the GCM of Showman et al. (2009) without TiO and VO.

Open with DEXTER

As mentioned above, Showman et al. (2009) do not provide results for their model of HD 209458b without TiO and VO, which prevents direct comparison. The 4.5 μm phase curve from their model without TiO and VO is, however, presented in Zellem et al. (2014), and we have plotted it in Fig. 11 as a dashed line. The model 4.5 μm phase curves agree reasonably well, with minor differences in the night side emission. Interestingly, both models significantly overestimate the night side flux, indicating that this is a common feature of current GCMs. Zellem et al. (2014) suggest that this may be due to non-equilibrium carbon chemistry, specifically vertical quenching of CH4, leading to larger CH4 abundances and consequently more efficient night side cooling. Another potential explanation could be horizontal quenching of CO increasing the abundance of CO on the night side relative to that predicted by equilibrium. As CO has a strong absorption feature at 4.5 μm this could potentially decrease the emission at this wavelength.

We have previously found that the effect of non-equilibrium chemistry has likely been overestimated in previous studies due to the inconsistent treatment of non-equilibrium chemistry and its feedback on the PT profile Drummond et al. (2016). In addition, most previous studies are limited to considering vertical non-equilibrium effects due to the 1D models used. To investigate the effect of horizontal mixing in these atmospheres, which could dominate over vertical mixing due to the large horizontal wind speeds, we are currently coupling a non-equilibrium chemistry network to the GCM, which may provide a solution to our current inability to reproduce the observed night side emission.

The lack of a temperature inversion potentially caused by TiO and VO in the atmosphere of HD 209458b has been suggested to be due to a cold-trap (Hubeny et al. 2003; Spiegel et al. 2009; Parmentier et al. 2013, 2016). TiO and VO could potentially be advected from the day side to the night side, condense, rain out and be trapped on night side at high pressures. We note that, as these calculations were based on GCM results where the model, like ours, has not equilibrated for pressures P105 Pa. Temperatures at these pressures may be higher than estimated by the 1D global average used to initialise the models, see Figs. 7 and 8, and we therefore emphasise the need to study the long-term evolution of the deep atmosphere of hot Jupiters, both to study potential cold-traps and radius inflation.

We note that all temperatures in our model initialised with the 1D globally averaged PT profile are below the condensation temperature of TiO/VO. Consequently, we would not expect gaseous TiO and VO to form even if we included their opacity. They may form, however, shortly after the simulations start as the models are initialised with zero advection, causing the day side to initially heat up to temperatures potentially above the condensation temperature, and then cool down as the advection becomes more efficient. Studying the formation of TiO and VO with varying initial conditions in these models would therefore be beneficial in order to understand the robustness of the presence of these gases in hot Jupiter atmospheres.

5. Conclusions

We have presented results from the first application of the UM, including a sophisticated and accurate radiation scheme, to hot Jupiters. We have performed comparisons to both the model of Mayne et al. (2014a) which employs Newtonian forcing, to the SPARC/MITgcm of Showman et al. (2009) and to existing observations. Our main conclusions are

  • Our models with Newtonian forcing and sophisticated radiativetransfer have qualitatively similar wind and temperaturepatterns: a hotspot shifted eastward of the substellar point, a broadeastward equatorial jet and a westward mean flow at higherlatitudes. This is in good qualitative agreement with resultsobtained with the SPARC/MITgcm (Showmanet al. 2009).

  • The day side emission and phase offset of the coupled model fit the observed 4.5 μm phase curve reasonably well, while the night side emission is about a factor of two too large. This is similar to the results of Showman et al. (2009), and potential explanations are effects of non-equilibrium chemistry and super-solar metallicities. We are in the process of investigating this by coupling a non-equilibrium chemistry scheme to our GCM.

  • We do not see the vertical coherence of temperatures, that is, the small variation in the location of the hottest and coldest points in the atmosphere with pressure, seen in the SPARC/MITgcm. The cause of this difference is unclear, but the vertical coherence is not expected from simple time scale arguments.

  • The deep atmosphere has not converged to a steady-state even though we have run the model for 1600 d, and we see evidence for significant deviations from the globally averaged 1D radiative-convective equilibrium profile used to initialise the model. This highlights the need to use very long runs when evaluating the feasibility of cold-traps and studying mechanisms for radius inflation.

  • The observed differences between the UM and the SPARC/MITgcm highlight the importance of model intercomparisons, which are needed to separate physically robust results from model degeneracies. There is a need to compare the dynamical cores in more detail, but also both the radiation schemes and post-processing tools in isolation. In addition it will be important study the sensitivity of the model to the initialisation, particularly when including strong absorbers of stellar radiation with high condensation temperatures such as TiO and VO.

Despite the difficulties discussed here, GCMs have seen rapid improvement in their application to hot Jupiters. Combined with ever improving observations, as well as extensive intercomparison exercises, they will enable us to significantly improve our understanding of these exotic planets.


3

Data collected from http://refractiveindex.info/

Acknowledgments

We would like to thank Jonathan Tennyson and Travis Barman for insightful discussions. 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). DSA acknowledges support from the NASA Astrobiology Program through the Nexus for Exoplanet System Science. 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. DH acknowledges funding from the DFG through the Collaborative Research Centre SFB 881 “The Milky Way System”. 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.

References

All Tables

Table 1

Model parameters adopted for HD 209458b.

All Figures

thumbnail Fig. 1

Illustration of the definition of the angles used to calculate the hemispherically integrated emission, see Eq. (3). (θ,φ) and (ϑ,ϕ) denote the position and the direction of the radiation, respectively, with the coordinate system placed such that the z-axis points towards the observer.

Open with DEXTER
In the text
thumbnail Fig. 2

Horizontal wind velocity as arrows and temperature [K] as colours and contours from our models of HD 209458b after 1600 d. Left column: results from the model with Newtonian forcing discussed in Sect. 3 at 100 Pa, 3 × 103 Pa, 3 × 104 Pa, 1 × 105 Pa (from top to bottom, from Mayne et al. 2014a). Right column: results from the coupled model discussed in Sect. 4 at 3 Pa, 3 × 103 Pa, 3 × 104 Pa, 1 × 105 Pa (from top to bottom).

Open with DEXTER
In the text
thumbnail Fig. 3

Zonal mean of the zonal wind velocity [m s− 1] after 1600 d for the models of HD 209458b. Left: model with Newtonian forcing from Mayne et al. (2014a) discussed in Sect. 3. Right: coupled model discussed in Sect. 4. Red indicates a prograde wind, blue indicates a retrograde wind.

Open with DEXTER
In the text
thumbnail Fig. 4

PT profiles around the globe after 1600 d for the model of HD 209458b with Newtonian forcing from Mayne et al. (2014a). Red solid lines and blue dashed-dotted lines are day and night side profiles, respectively, at latitude. Magenta dashed lines and cyan dotted lines are profiles between and 90° latitude for longitudes 180° and , respectively.

Open with DEXTER
In the text
thumbnail Fig. 5

Observed (points) and synthetic (lines) emission spectra for HD 209458b models with Newtonian cooling. The solid line has been calculated from the model presented in Mayne et al. (2014a) and Sect. 3 using ATMO, the dashed line is the synthetic emission spectrum from Fortney et al. (2006), which is based on the models of Cooper & Showman (2005, 2006) using the same Newtonian forcing scheme. The black points are observations from Swain et al. (2008) (●), Crossfield et al. (2012) (▼), Deming et al. (2005) (◄), Zellem et al. (2014) (►), Diamond-Lowe et al. (2014) (■) and Evans et al. (2015) (♦).

Open with DEXTER
In the text
thumbnail Fig. 6

Synthetic Spitzer IRAC phase curves from the Newtonian forcing model. The solid lines are calculated from the model presented in Mayne et al. (2014a) and Sect. 3 using ATMO, the dashed lines are the synthetic phase curves from Fortney et al. (2006), which are based on the models of Cooper & Showman (2005, 2006) using the same Newtonian forcing scheme. The models have been integrated over the IRAC bands using the filter functions. The data points are from Zellem et al. (2014) (►), Diamond-Lowe et al. (2014) (■) and Evans et al. (2015) (♦). The best fit to the observed 4.5 μm phase curve from Zellem et al. (2014) is shown as a solid black line, the grey shaded area is the 1σ uncertainty for the offset of the observed planet to star flux ratio.

Open with DEXTER
In the text
thumbnail Fig. 7

PT profiles around the globe after 1600 d for the coupled model of HD 209458b discussed in Sect. 4. Red solid lines and blue dashed-dotted lines are day and night side profiles, respectively, at latitude. Magenta dashed lines and cyan dotted lines are profiles between and 90° latitude for longitudes 180° and , respectively. The black line is the initial PT profile adopted.

Open with DEXTER
In the text
thumbnail Fig. 8

PT profiles around the globe after 1600 d for the coupled model of HD 209458b initialised with a PT profile that is 400 K (left) and 800 K (right) hotter than the global 1D mean. Lines are as in Fig. 7.

Open with DEXTER
In the text
thumbnail Fig. 9

Temperature difference ΔT between the PT profiles in the right-hand panel of Fig. 8, which are from the model with a 800 K hotter initial condition compared to the global 1D mean, and the PT profiles in Fig. 7, which are from the model initialised with the global 1D mean PT profile. Temperature differences are small for P105 bar, while the models have clearly not reached a steady-state for P105 bar. Results are similar for the case with a 400 K hotter initial condition. Lines are as in Fig. 7.

Open with DEXTER
In the text
thumbnail Fig. 10

Same as Fig. 5, but the synthetic emission spectrum has been calculated using results from our coupled model (Sect. 4) with ATMO.

Open with DEXTER
In the text
thumbnail Fig. 11

Same as Fig. 6, but the synthetic emission spectrum has been calculated using results from our coupled model (Sect. 4) with ATMO. The dashed line is the model 4.5 μm phase curve from Zellem et al. (2014) obtained using the GCM of Showman et al. (2009) without TiO and VO.

Open with DEXTER
In the text

Current usage metrics show cumulative count of Article Views (full-text article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.

Data correspond to usage on the plateform after 2015. The current usage metrics is available 48-96 hours after online publication and is updated daily on week days.

Initial download of the metrics may take a while.