© ESO 2020

1 Introduction

Observations of exoplanets, particularly of transiting, tidally–locked hot Jupiters, have yielded the detection of several atmospheric species (e.g. Snellen et al. 2008; Wilson et al. 2015; Sing et al. 2015; Kreidberg et al. 2015), alongside evidence for aerosols (e.g. photochemical hazes or condensate clouds, see for example, Sing et al. 2016). Observations of some hot hydrogen-dominated atmospheres show clear absorption or emission features of absorbing gas-phase species (e.g. Nikolov et al. 2018); while for other planets, these features appear to be obscured, possibly by the presence of absorbing or scattering clouds and hazes (e.g., Deming et al. 2013).

Theoretical modelling of hot Jupiter atmospheres has recently shifted its focus to chemical composition, in terms of both the gas-phase chemistry and cloud formation. Dynamical mixing can alter the local abundances of chemical species in an atmosphere. If the transport of material occurs on a timescale that is faster than the timescale for the atmosphere to reach a state of chemical equilibrium, then advection of chemical species has a strong influence on the local chemical composition. This process of “quenching” was inferred several decades ago in the atmosphere of Jupiter since observed abundances of carbon monoxide (CO) did not agree with predictions based on local chemical equilibrium (e.g. Prinn & Barshay 1977). Vertical quenching, due to vertical mixing, has also been inferred from observations of young Jupiter analogues (e.g. Skemer et al. 2014) and brown dwarfs (see Saumon et al. 2003; Zahnle & Marley 2014, and references therein).

Until recently, studies of the effect of dynamical mixing on the chemical composition have mainly been restricted to the use of vertical diffusion in one-dimensional (1D) atmosphere models (e.g. Zahnle et al. 2009; Moses et al. 2011; Venot et al. 2012; Rimmer & Helling 2016; Drummond et al. 2016; Tsai et al. 2017), primarily because of the high computational expense of the chemistry calculations, making higher-dimensional simulations challenging. However, for short–period, tidally–locked exoplanets material is also expected to be mixed horizontally between a permanently irradiated dayside and an unirradiated nightside. For hot Jupiters, fast horizontal winds (e.g. Showman & Guillot 2002; Perez-Becker & Showman 2013) acting to redistribute heat from the dayside to nightside have been inferred from day–night temperaturecontrasts and shifting of the thermal emission, or ‘hot spot’, from the closest point to the star (e.g. Knutson et al. 2007, 2012; Zellem et al. 2014; Wong et al. 2016). These winds have also been observed more directly using high–resolution spectroscopy (e.g. Snellen et al. 2010), indicating “diverging” flows from the day to night side (Louden & Wheatley 2015). Observations of the emission as a function of orbital phase, or phase curves, have also been shown to depart from model predictions for the nightside, which has been proposed to be caused by horizontal mixing of chemical species (e.g. Stevenson et al. 2010; Knutson et al. 2012; Zellem et al. 2014).

Future facilities, such as the James Webb Space Telescope (JWST; e.g. Bean et al. 2018) and the Atmospheric Remote-sensing Infrared Exoplanet Large-survey (ARIEL; e.g. Tinetti et al. 2018), will provide more accurate phase curves, for an increasing number of targets, alongside the potential for more accurate measurements of chemical abundances (e.g. methane, CH₄) from emission and transmission spectra. Combining these observations with a more complete understanding of the mixing processes controlling the observed abundances, and phase curves, will provide a unique “window” into the atmospheric dynamics of tidally–locked exoplanets.

Cooper & Showman (2006) were the first to investigate the effect of three-dimensional (3D) atmospheric mixing on the composition of hot Jupiter atmospheres using a simplified chemistry scheme (chemical relaxation) and temperature forcing within a 3D hydrodynamics code. Focusing on HD 209458b, the abundance of CO was “relaxed” back to an equilibrium abundance over a timescale chosen to be indicative of the relevant chemical pathways. The results of Cooper & Showman (2006) indicated that vertical mixing was actually dominant, in terms of controlling the composition in the observable region of the atmosphere, compared to the horizontal mixing. Agúndez et al. (2014a) also investigated the role of vertical and horizontal mixing using a more complete chemistry scheme but a highly idealised approach to the dynamics. Their pseudo two-dimensional (2D) model approximated horizontal mixing with a time-varying pressure–temperature profile within a 1D atmosphere model, intended to represent a column of atmosphere rotating around the equator within the equatorial zonal jet. Agúndez et al. (2014a) concluded that horizontal mixing was also important in controlling the local chemical composition, finding that the nightside can be “contaminated” by material transported from the hotter dayside.

Recently, Drummond et al. (2018a,b) revisited the chemical relaxation scheme developed by Cooper & Showman (2006) and implemented it into a different 3D atmosphere model. The main improvement over the earlier work of Cooper & Showman (2006) was the coupling of the chemistry scheme to a sophisticated radiative transfer scheme (Amundsen et al. 2014, 2016). This improvement allows changes in the local composition, due to advection, for instance, to be consistently fed back to the heating rates and the thermal evolution of the atmosphere, as well as the calculation of synthetic observations (e.g. Lines et al. 2018a). Drummond et al. (2018a,b) found that both horizontal and vertical quenching play important roles in determining the abundance of CH₄ in the observable region of the atmosphere for both HD 209458b and HD 189733b. Mendonça et al. (2018) recently included an updated chemical relaxation scheme, developed in Tsai et al. (2018), in a 3D atmosphere model for the case of WASP–43b. In contrast to the results of Drummond et al. (2018a,b), who found that chemical gradients are very efficiently removed due to advection, they found that significant horizontal structure in the chemical composition remains. However, a fair comparison cannot be made since there are many differences in the implementation of the model and the physics schemes, as well as the planetary parameters assumed, between our earlier work and that of Mendonça et al. (2018).

Bordwell et al. (2018) considered the quenching of chemical species by modelling the (2D and 3D) fluid dynamics of the convective region of substellar atmospheres including a reactive chemical tracer. Their results suggest that the length scale associated with the quenching is dependent on the chemical properties of the species of interest, in agreement with the earlier 1D work of Smith (1998). Zhang & Showman (2018a,b) investigated the eddy diffusion coefficient (K_zz) using results from 3D simulations including a transported reactive chemical species and found that K_zz should in principle be different for different chemical species. They derived also an analytical theory for the quench point of chemical species in substellar atmospheres.

As well as the gas-phase chemistry, the study of the formation of clouds and their effects on the atmosphere have recently been a significant focus. In the framework of 3D modelling there are two main approaches to treating cloud formation. One approach is to include cloud formation as a post-processing step, using output from a cloud-free radiation-hydrodynamics code, effectively decoupling the dynamic and thermodynamic evolution of the atmosphere from the cloud treatment (e.g. Parmentier et al. 2016; Helling et al. 2016; Powell et al. 2019). Another option is to take a more complete approach by coupling treatments of the cloud formation and their radiative properties directly to an atmosphere model, allowing for cloud-radiative feedback on the atmosphere (e.g. Lee et al. 2016; Lines et al. 2018b,a, 2019). In this work we choose to ignore the presence of clouds and focus on the gas-phase chemical composition. Future works will need to treat the kinetics of both gas-phase and cloud chemistry to fully understand the composition of hot Jupiter atmospheres. However, coupling these two processes consistently in a numerical model is a very significant scientific and technical challenge.

We present a significant step towards a better understanding of the chemical composition of hot exoplanet atmospheres by focusing in detail on the gas-phase chemistry. We couple our 3D model, including sophisticated treatments of the radiative transfer and dynamics (as used in, e.g. Amundsen et al. 2016; Drummond et al. 2018c; Mayne et al. 2019), to a chemical kinetics scheme using a chemical network specifically developed for application in 3D atmosphere models (Venot et al. 2019). This model allows for the interactions and feedback between chemistry, radiation, and dynamics.

In Sect. 2 we introduce the model development to this point, focusing on the chemistry (Sect. 2.1), and detail the parameters used in the simulations that we present (Sect. 2.2). We then move to our results in Sect. 3 demonstrating that including 3D transport is crucial in understanding the atmospheric chemical structure of hot Jupiters, and that the effects of chemical transport have significant consequences for their predicted spectra. We discuss our results in the context of previous studies in Sect. 4 and finally, we state our conclusions in Sect. 5.

2 Methods and model description

We use an idealised configuration of the Met Office Unified Model (UM) to simulate the atmospheres of highly-irradiated gas-giant planets, namely HD 209458b and HD 189733b. The UM, originally developed for Earth atmosphere applications, has previously been used to simulate the atmospheres of hot Jupiters and warm Neptunes (Mayne et al. 2014a, 2017, 2019; Amundsen et al. 2016; Drummond et al. 2018a,b,c; Lines et al. 2018a,b, 2019), as well as terrestrial exoplanets (Mayne et al. 2014b; Boutle et al. 2017; Lewis et al. 2018; Yates et al. 2020). Our modifications to the model for this work were made in a development branch using the UM at version 11.4 (released June 2019) as its base.

The dynamical core of the UM, called ENDGame (Wood et al. 2014), has been validated for exoplanet applications by reproducing long-term, large-scale flows of idealised standard tests (Mayne et al. 2014b) and tested for the typical flows expected in hot Jupiter like atmospheres (Mayne et al. 2014a, 2017). ENDGame has the unique capability to seamlessly transition between the simplified “primitive” equations of motion and the “full” equations of motion, the latter not adopting a series of assumptions(e.g. hydrostatic balance, shallow-atmosphere) that are used to simplify the equation set (see Mayne et al. 2014a, for discussion). Using this, Mayne et al. (2019) demonstrated that the primitive equations produce a significantly different circulation, and hence also thermal structure, compared with the full equation set for warm small Neptune and super Earth atmospheres. For typical hot Jupiter atmospheres, however, the differences in the resulting flow are negligible when using the primitive or full equations (Mayne et al. 2014a, 2019). In this work we use the full equations of motion.

The UM uses a flexible and sophisticated radiative transfer scheme, the open-source Suite Of Community RAdiative Transfer codes based on Edwards and Slingo (SOCRATES) scheme¹ (Edwards & Slingo 1996; Edwards 1996), to calculate the radiative heating rates that, in part, drive the thermal evolution of the atmosphere. SOCRATES has been updated and tested for the high-temperatures and compositions relevant for hot Jupiter atmospheres (Amundsen et al. 2014, 2016, 2017) and is now being used in other models in the community (Way et al. 2017, 2018; Thomson & Vallis 2019).

SOCRATES adopts the two-stream approximation and treats opacities using the correlated-k approximation, the latter significantly improving computational efficiency while maintaining accuracy, compared with alternative methods (Lacis & Oinas 1991; Amundsen et al. 2014). The k-coefficients of individual gases are combined on-the-fly using the method of equivalent extinction (Edwards 1996; Amundsen et al. 2017), meaning that the composition of the gas does not need to be known a priori, unlike the commonly used alternative method of “pre-mixed” k-coefficient tables (e.g. Showman et al. 2009). In the case of pre-mixed k-coefficients a chemical equilibrium composition is typically assumed, with abundances tabulated as a function of pressure and temperature. When including 3D chemical mixing combining opacities of individual gases on-the-fly is extremely important as the abundances in a grid cell are not necessarily dependent on the local pressure and temperature, if the chemistry is not in equilibrium.

The model includes absorption opacity due to H₂O, CO, CO₂, CH₄, NH₃, HCN, Na, K, Li, Cs, and Rb, as well as collision-induced absorption (CIA) due to H₂ –H₂ and H₂–He. We also include Rayleigh scattering due to H₂ and He. We note that CO₂ and HCN have been added to the model for this work and were not included in our previous modelling studies of hot Jupiters. The implementation of opacities in the model is described in detail in Amundsen et al. (2014, 2016, 2017) and we refer the reader to Goyal et al. (2018) for the most up-to-date information on the sources of line-lists used in the model. Where available we use line lists from the ExoMol project².

The treatment of chemistry in the exoplanet configuration of the UM has been incrementally improved and presented in a series of recent papers. Initially, Amundsen et al. (2016) used a simple and extremely efficient analytical solution to equilibrium chemistry (Burrows & Sharp 1999) for the species H₂O, CO, CH₄, and NH₃, combined with a simple parameterisation for the abundances of alkali species (Na, K, Li, Cs, and Rb) based on their condensation pressure-temperature curves. Drummond et al. (2018c) coupled a Gibbs energy minimisation scheme to the UM, also a local chemical equilibrium scheme, increasing the flexibility of the chemistry by allowing the numerical calculation of local chemical equilibrium abundances of a very large number of chemical species over a wide range of pressure, temperature, and element abundances.

Most recently, Drummond et al. (2018a,b) implemented a chemical relaxation scheme, based on the earlier work of Cooper & Showman (2006), which considers the time-dependent evolution of CO, CH₄, and H₂O and allows therefore the assumption of local chemical equilibrium to be relaxed for these species. This enables effects such as 3D chemical transport to be included, which is not possible when using a local chemical equilibrium scheme. Using a chemical relaxation scheme involves the estimation of a chemical timescale that represents either the time for interconversion to take place between two chemical species (e.g. for CO and CH₄, Cooper & Showman 2006) or the time taken for a species to evolve back to its local chemical equilibrium, following a perturbation (e.g. Tsai et al. 2018). In Cooper & Showman (2006) the timescale for CO–CH₄ interconversion was derived using the rate-limiting step in the series of elementary reactions that link the two species. This approach was questioned in recent work by Tsai et al. (2018).

In this work, we further improve the treatment of chemistry in the UM by coupling a flexible chemical kinetics scheme that was previously developed within the framework of the 1D atmosphere code ATMO (Tremblin et al. 2015; Drummond et al. 2016). This improves on our earlier work (Drummond et al. 2018a,b) by allowing for the inclusion of more chemical species, most importantly those that are included as opacities in the model. In addition, the chemical network that we use (Venot et al. 2019) is derived from a larger network that has been experimentally validated across a wide-range of pressure and temperature (Venot et al. 2012), likely leading to a more accurate chemical evolution compared with the Cooper & Showman (2006) chemical relaxation formulation.

In Sect. 2.1 we describe the coupling of the chemical kinetics scheme to the UM. Then, in Sect. 2.2, we describe the setup of the model including the various planet-star specific parameters for each simulation.

2.1 Chemical kinetics in a 3D atmosphere model

Implementing a chemical kinetics scheme in a 3D model such as the UM requires its coupling with the dynamics and radiative transfer schemes, and the selection of a suitable chemical network.

2.1.1 Coupling the chemical kinetics scheme to the model

The chemical kinetics scheme that we couple to the UM was initially developed within the framework of the 1D atmosphere code ATMO. The scheme has previously been described in Drummond et al. (2016). The methods are similar to those taken in other chemical kinetics codes applied to hot exoplanet atmospheres presented in the literature (e.g. Moses et al. 2011; Kopparapu et al. 2012; Venot et al. 2012; Agúndez et al. 2014b; Tsai et al. 2017). We give a brief overview here.

A chemical kinetics scheme requires the solution of a system of coupled ordinary differential equations (ODEs) that describe the time evolution of a set of chemical species. A chemical equilibrium scheme, in contrast, requires only the solution of algebraic equations. The terms in the ODEs correspond to chemical production and loss terms that are calculated based on a chosen chemical network, which describes the chemical reactions that link each of the chemical species. One of the fundamental choices to make when developing a chemical kinetics scheme is how to solve this system of ODEs.

We use the publicly available Fortran library ODEPACK (Hindmarsh 1983) to solve the system of ODEs. Specifically, we use the DLSODES routine which is the variant of the package suitable for stiff systems, which uses an implicit Backward Differentiation Formulae (BDF) method (Radhakrishnan & Hindmarsh 1993). Several user-defined options are required inputs for DLSODES. We use a relative tolerance 1 × 10⁻³ and we use an order of 2 in the BDF method (i.e. the results from two previous iterations are used in the calculation for the current iteration). The step size of each iteration is internally generated by DLSODES, though we impose minimum and maximum timesteps of 1 × 10⁻¹⁰ s and 1 × 10¹⁰ s, respectively.The user must also set a method flag (MF) to control the operation of DLSODES. Here we use MF = 222 which is the option for stiff problems using an internally generated (numerically derived) Jacobian, as opposed to user-supplied analytical Jacobian.

Tracer advection is handled by the well-tested semi-implicit, semi-Lagrangian advection scheme of the ENDGame dynamical core (Wood et al. 2014). The tracers (one for each chemical species, totalling 30 in this case) are initialised to chemical equilibrium abundances at the start of the model integration, calculated using the coupled Gibbs energy minimisation scheme (Drummondet al. 2018c). The quantity that is advected is the mass mixing ratio.

We call the chemical kinetics scheme, that handles the chemical evolution, individually for each cell of the 3D model grid on a frequency that corresponds to a chemical timestep. The chemical timestep is longer than the dynamical timestep meaning that the chemistry iscalled less frequently than every dynamical timestep, motivated by improving the computational efficiency. Asimilar approach is taken for the radiative transfer calculations (Amundsen et al. 2016). Details of the dynamical, radiative, and chemical timesteps are given in Sect. 2.2, and these were selected to provide accuracy and efficiency. Since the chemistry calculations are done separately from the advection step (i.e. operator splitting), the process could be thought of as solving a series of 0D box models, one for each cell of the 3D grid.

The initial values for chemical integration are the current values of the advected tracers converted from mass mixing ratio into number densities using the total gas number of that grid cell. The time-stepping in the chemical kinetics iterations is done using an adaptive timestep, starting with an initial small timestep of 1 × 10⁻¹⁰ s. Subsequent timesteps are controlled by the output of the DLSODES routine, as described above. Chemical iterations continue until either the integration time is equal to the set chemical timestep or the chemical system has reached a steady-state. Steady-stateis assessed via the maximum relative change in the number densities of the individual species in an approach similar to Tsai et al. (2017). We find that continuing iterations significantly beyond the point at which the system has reached a steady-state can lead to numerical errors. This means that for grid cells in the deep high-pressure, high-temperature region where the chemistry reaches a steady-state relatively quickly (as the chemical timescale is short) we may integrate the chemistry for a time less than the chemical timestep. For cells at lower pressures and temperatures (where the chemical timescale is long) the chemistry iterations are more likely to be stopped by reaching an integration time equivalent to the chemical timestep.

2.1.2 Chemical network

There are now a large number of chemical networks that have been used to model the atmospheres of hot exoplanets (e.g. Liang et al. 2003; Zahnle et al. 2009; Moses et al. 2011; Venot et al. 2012; Rimmer & Helling 2016; Tsai et al. 2017). Since the early work of Liang et al. (2003) the trend has been for such networks to grow in size and complexity. However, more recently the focus has shifted to the development of smaller chemical networks (Tsai et al. 2017; Venot et al. 2019) that are less computationally demanding, making them feasible for use in a 3D atmosphere model.

We use the “reduced” chemical network recently developed by Venot et al. (2019). The network comprises 30 chemical species and 181 reversible reactions. This network was constructed by reducing the larger (105 species, ~ 1000 reversible reactions) network of Venot et al. (2012) using commercial software. The original network of Venot et al. (2012) was experimentally validated for a wide-range of pressure and temperature relevant to hot Jupiter atmospheres. The Venot et al. (2012) network, excluding photodissociations, was reduced by eliminating chemical species and reactions while maintaining a set level of accuracy for six target species: H₂O, CH₄, CO, CO₂, NH₃, and HCN. The reduced network was compared against the original network within a 1D chemical-diffusion model and results for the six target species agreed extremely well across a remarkably large parameter space (pressure, temperature, and element abundances). The network is publicly available via the KIDA database³ (Wakelam et al. 2012).

Recently, Venot et al. (2020) released updates to their original chemical network (Venot et al. 2012) with a focus on improvements to the chemistry of methanol. They also re-derived a reduced chemical network from the updated network. Negligible differences were found in the quench points of several chemical species, using a 1D photochemical model, between theoriginal and updated networks for hot Jupiters (specifically HD 209458b and HD 189733b were tested). On the other hand, significant differences were found for cooler atmospheres. As we only consider hot Jupiter atmospheres in this study, our use of the Venot et al. (2019) chemical network, based on the original Venot et al. (2012) network without updates to the methanol chemistry, does not affect our results.

We note that we only consider thermochemical kinetics and transport-induced quenching in this work and do not include photochemistry (photodissociations or photoionisations). Nonetheless the model developments presented in this work represent a significant step forwards in exoplanet climate modelling. Inclusion of photochemistry represents a significant further model development.

2.2 Model setup and simulation parameters

The basic model setup is similar to that used for previous hot, hydrogen-dominated atmosphere simulations with the UM (e.g. Amundsen et al. 2016). We use 144 longitude points and 90 latitude points giving a horizontal spatial resolution of 2.5° and 2° in longitude and latitude, respectively. The UM uses a geometric height-based grid, for which we use 66 levels, as opposed to a pressure-based grid (see Sect. 2.5 of Mayne et al. 2014a, for details). For each planet setup we adjust the height at the top-of-atmosphere (z_toa) to reach approximately the same minimum pressure range (~1 Pa) for the initial profile. The bottom boundary of the model has a pressure of ~ 2 ×10⁷ Pa. The model includes a vertical damping “sponge” layer to reduce vertically propagating waves resulting from the rigid boundaries of the model, and we use the default setup for this as described in Mayne et al. (2014a) and Amundsen et al. (2016).

For the planetary and stellar parameters we use the values from the TEPCat database⁴, which are given in Table 1. For the stellar spectra we use the Kurucz spectra for HD 209458 and HD 189733⁵.

For the element abundances we use a solar composition with values taken from Caffau et al. (2011). The model requires global values of the specific heat capacity (c_P) and the specific gas constant ($R=\overline{R}/\mu$, where $\overline{R}=8.314 \text{J} {\text{K}}^{-1} {\text{mol}}^{-1}$ is the molar gas constant and μ is the mean molecular weight). We calculatevertical profiles of c_P and R for the initial 1D thermal profiles and use the straight average of these profiles for the globally uniform values required as model inputs, as previously described in Drummond et al. (2018c).

We use a main model (dynamical) timestep of 30 s which we find is a good compromise between efficiency and stability (Amundsen et al. 2016). The radiative heating rates are updated every 5 dynamical timesteps, giving a radiative transfer timestep of 150 s. The chemistry scheme is called every 125 dynamical timesteps giving a chemical timestep of 3750 s or approximately 1 h. Using a more efficient model setup (replacing the radiative transfer with a Newtonian cooling approach) we found little difference in the model result after 1000 (Earth) days of simulation between a model using a chemical timestep of 10 dynamical timesteps and 125 dynamical timesteps. We therefore choose to use a longer chemical timestep here simply for improved computational efficiency.

The model is initialised at rest (i.e. zero wind velocities) and with a horizontally-uniform thermal structure. The initial vertical thermal profile is calculated using the 1D radiative-convective equilibrium model ATMO, using the same planetary and stellar parameters Table 1. We use a zenith angle of cos θ= 0.5 for the ATMO calculation of the initial thermal profile, representative of a dayside average (e.g. Fortney & Marley 2007). For the chemical kinetics simulations the chemical tracers are initialised to the chemical equilibrium values in each grid cell. Since the initial thermal field is horizontally uniform, this entails that the initial chemical tracer values are also horizontally uniform.

We integrate the model for 1000 Earth days (i.e. 8.64 × 10⁷ s) by which point the maximum wind velocities and chemical abundances have approximately ceased to evolve (see Appendix A). A common problem with simulations of gas-giant planet atmospheres is the slow evolution of the deep atmosphere. Here, where the atmosphere is initialised with a “cold” initial thermal profile, the radiative and dynamical timescales are too long to feasibly reach a steady-state due to computational limitations (Rauscher & Menou 2012; Mayne et al. 2014a; Amundsen et al. 2016; Sainsbury-Martinez et al. 2019). For HD 209458b we find no significant differences when initialising with a hotter thermal profile (Appendix C). Conservation of several quantities is presented in Appendix B.

Chemical kinetics calculations are typically computationally demanding because they require the solution of a system of stiff ODEs. In our simulations, including the chemical kinetics calculations increases the computational cost by at least a factor of two, compared with an otherwise similar setup that uses the Gibbs energy minimisation (chemical equilibrium) scheme. However, the simulations remain computationally feasible and required approximately one-week of walltime using ~ 200 cores on the DiRAC DIaL supercomputing facility⁶, as an example. Using a larger chemical network would increase the computational cost because (1) chemistry calculations would become more expensive and, (2) more advected tracers would be required.

In Sect. 3 we present results from four simulations in total. We model the atmospheres of HD 209458b and HD 189733b both under the assumption of local chemical equilibrium and including the effect of 3D transport. For the former we use the coupled Gibbs energy minimisation scheme (Drummond et al. 2018c) to calculate the chemical abundances and we refer to these simulations throughout simply as the “equilibrium” simulations. For the latter we use the newly coupled chemical kinetics scheme and include tracer advection and we refer to these simulations as the “kinetics” simulations.

Fig. 1 Zonal-mean zonal wind velocity for the equilibrium simulations of HD 209458b (top) and HD 189733b (bottom).

3 Results

We first focus on the overall circulation and thermal structure resulting from our simulations, before moving on to describe the chemical structure, the chemical-radiative feedback, and the impacts on the synthetic observations.

3.1 Circulation and thermal structure

The temporal-mean (800–1000 days) zonal-mean zonal wind for the equilibrium simulations of HD 209458b and HD 189733b are shown in Fig. 1. For both planets, the circulation is dominated by a superrotating equatorial zonal jet, which is typical of close-in, tidally-locked atmospheres and is shown to be present across a wide parameter-space (e.g. Kataria et al. 2016). The maximum zonal-mean zonal wind velocity is ~ 6 km s⁻¹ in both cases. At higher latitudes, and also at higher pressures, the flow is westward (retrograde) with maximum wind velocities of ~ 1 km s⁻¹.

The maximum zonal-mean zonal wind for HD 209458b in previous results in the literature varies between 5 km s⁻¹ and 6 km s⁻¹ (Showman et al. 2009; Kataria et al. 2014) while for HD 189733b it varies between 3.5 km s⁻¹ and 6 km s⁻¹ (Showman et al. 2009; Kataria et al. 2014; Lee et al. 2016; Flowers et al. 2019). The jet speed depends on the magnitude of dissipation implemented in the model (required for numerical stability) and is currently unconstrained by observation (Heng et al. 2011; Mayne et al. 2014a).

The atmospheric temperature on the 10⁴ Pa isobar for the equilibrium simulations of both planets is shown in Fig. 2. The thermal structure is qualitatively quite similar between the two planets, with a hot dayside and relatively much cooler nightside and a hotspot that is shifted eastward of the substellar point. HD 209458b is generally a few hundred Kelvin warmer than HD 189733b.

Figure 3 shows pressure-temperature profiles extracted from the 3D grid of equilibrium simulations of HD 209458b and HD 189733b, as well as dayside and nightside average profiles. The 1D radiative-convective equilibrium pressure-temperature profile used to initialise each 3D simulation is also shown. Both atmospheres have significant zonal temperature gradients for P ≲ 10⁵ Pa with maximum day-night temperaturecontrasts of ~500 K. There is also a significant latitudinal thermal gradient, at a given longitude and pressure, of several hundred Kelvin between the equator and mid-latitudes.

A high-pressure thermal inversion at P ~ 10⁶ Pa is present inboth equilibrium simulations. This is a result of the deep atmosphere slowly converging towards a temperature profile that is hotter than the 1D radiative-convective equilibrium profile used to initialise the model (shown in black in Fig. 3). Using a 2D steady-state circulation model, Tremblin et al. (2017) showed that the steady-state of the deep atmosphere is actually significantly hotter than predicted by 1D radiative-convective equilibrium models, likely due to the advection of potential temperature. It has also been noted that 3D simulations appear to slowly evolve towards hotter profiles in the deep atmosphere (Amundsen et al. 2016), however very long timescale simulations are required to actually reach the end state (as performed in Sainsbury-Martinez et al. 2019). Therefore, the results for pressures greater than ~ 10⁶ Pa are likely to be strongly effected by the initial condition.

The atmospheric circulation and thermal structure from the equilibrium simulations of HD 209458b and HD 189733b are almost identical to results presented in our previous work (Drummond et al. 2018a,b). Minor differences result from a number of changes in the model setup and chosen parameters, including, in no particular order: use of a different chemical equilibrium scheme (Gibbs energy minimisation versus analytical chemical equilibrium formulae), slight adjustments to some planetary and stellar parameters (compare the parameters for HD 209458b and HD 189733b in Table 1 with those presented in Drummond et al. 2018a,b, respectively), and the addition CO₂ and HCN as opacity sources which were not previously included. Differences in the wind velocities and temperatures between the equilibrium and kinetics simulations due to chemical-radiative feedback are discussed later in Sect. 3.3.

Fig. 2 Atmospheric temperature (colour scale and black contours) on the 1 × 104 Pa (0.1 bar) isobar with horizontal wind velocity vectors for the equilibrium simulations of HD 209458b (top) and HD 189733b (bottom).

Fig. 3 Pressure-temperature profiles extracted from the 3D grid of equilibrium simulations of HD 209458b (top) and HD 189733b (bottom). Solid grey lines show profiles at various longitude points around the equator (0° latitude) while dashed grey lines show profiles at various longitude points at a latitude of 60°. Area-weighted dayside average (red) and nightside average (blue) profiles are shown, along with the 1D radiative-convective equilibrium model profile used to initialise the 3D model (black).

3.2 Chemical structure

To compare the differing chemical structure resulting from either the equilibrium or kinetics approach, we focus on each planet separately.

3.2.1 HD 209458b

Figure 4 shows the mole fractions of a subset of chemical species (CH₄, CO, H₂O, CO₂, NH₃, and HCN) on the P = 10⁴ Pa isobar for the equilibrium and kinetics simulations of HD 209458b. Figure 5 shows both their equilibrium and kinetics abundances as vertical (pressure) mole fraction profiles at a series of equally spaced longitude points around the equator (i.e. at a latitude of θ = 0). We also show the vertical abundance profiles for a latitude of 45° in Appendix D. We show this particular subset of six molecules (out of a total of 30 species included in the chemical network) since each of these species are included as opacity sources in the model and they are also the six “species” for which the reduced chemical network was constructed (Venot et al. 2019).

The equilibrium abundances (left column of Fig. 4) clearly trace the temperature structure (Fig. 2), as they necessarily should since the chemical equilibrium composition only depends on local pressure, temperature, and element abundances. The abundances of CH₄, CO₂, and NH₃ are lower on the warm dayside compared with the cooler nightside. The equilibrium abundance of CH₄ is large enough in the coolest regions at this pressure level to capture a significant number of carbon atoms, thereby reducing the equilibrium abundance of CO, though CO remains the dominant carbon species. The decrease in CO leads to a small increase in the equilibrium abundance of H₂O as more oxygen atoms are available. HCN shows a very small horizontal variation, since its abundance does not strongly depend on temperature for T >1000 K (e.g. Heng & Tsai 2016). However, this molecule does show a significant vertical abundance gradient (Fig. 5). These trends are simply due to the temperature and pressure dependence of the Gibbs free energy of each chemical species. The equilibrium abundances of CH₄, CO, and H₂O are almost identical to those presented in Drummond et al. (2018a), since the temperature structure is approximately the same.

We now focus on the differences in the chemical abundances between the equilibrium and kinetics simulations to investigate the effect of advection on the atmospheric composition. On the P = 10⁴ Pa isobar the abundance fields from the kinetics simulation (right column of Fig. 4) differ from the equilibrium simulation for all of the presented species. The kinetics abundances of CO and H₂O become horizontally more uniform compared with their equilibrium abundances, losing the variations across the cold, nightside mid-latitudes. The quantitative differences between the two simulations, however, are small and CO and H₂O are the dominant carbon and oxygen species throughout the atmosphere, as in the equilibrium case. CO₂ is close to its chemical equilibrium abundance for most of the domain, on this pressure level, except in the coldest region (nightside mid-latitudes) where there is a small difference between the equilibrium and kinetics simulations. Here CO₂ is maintaining a “pseudo-equilibrium” with CO and H₂O, since the cycling between oxidised forms of carbon remains efficient (e.g. Moses et al. 2011; Tsai et al. 2018). Therefore, CO₂ only shows an obvious difference with its chemical equilibrium value in the region where CO and H₂O change, at the nightside mid-latitudes. In other regions at this pressure level, the abundances of CO and H₂O do not change significantly, and therefore neither does the abundance of CO₂.

The effect of advection on the horizontal abundance distribution of CH₄ at P = 10⁴ Pa is more complex. The zonal (longitudinal) abundance gradient is decreased compared with the equilibrium simulation but there is a large latitudinal abundance gradient. The abundance of CH₄ from the kinetics simulation is around an order of magnitude smaller in the equatorial region compared with at higher latitudes. A similar qualitative structure is found for NH₃ and HCN, with lower abundances in the equatorial region compared with higher latitudes.

From the vertical (pressure) profiles of the equatorial abundances of the six chemical species that we focus on here (Fig. 5), it is clear that H₂O and CO are the most abundant trace species (i.e. after H₂ and He), showing negligible abundance variations with both pressure and longitude (for P > 10⁶ Pa). We note that the small variations in the abundances between the equilibrium and kinetics simulations that are apparent in Fig. 4 are mainly present at mid-latitudes and relatively small differences are also less apparent on the log-scale of this figure (Fig. 5). Both the kinetics abundances of NH₃ and HCN diverge from their respective equilibrium abundances between 10⁵ and 10⁴ Pa, with their abundances becoming zonally and vertically well mixed towards lower pressures.

The vertical profiles of the equatorial abundance of CH₄ are the most interesting. In the pressure range 10⁴−10⁵ Pa the kinetics abundances of CH₄ tend towards its equilibrium abundance in the hottest part of the atmosphere. This has the effect of significantly reducing the zonal abundance gradient, in this pressure range, by decreasing the nightside abundance. However, there remains a significant vertical abundance gradient in the same pressure range. This structure can be explained by CH₄ undergoing efficient zonal mixing in the pressure range 10⁴−10⁵ Pa, while not being affected by vertical mixing. As a result, zonal abundance gradients are homogenised but vertical abundance gradients persist. For pressures less than 10⁴ Pa the equatorial abundance of CH₄ becomes approximately vertically uniform.

CO₂ is an interesting molecule in that it stays close to its chemical equilibrium abundance for much of the equatorial region, except for the coolest nightside profiles at low pressures. As previously discussed, CO₂ maintains a pseudo-equilibrium with CO and H₂O, even once those molecules themselves have quenched (Moses et al. 2011; Tsai et al. 2018). We attribute the departure of the abundances of CO₂ from chemical equilibrium as being due to its pseudo-equilibrium with CO and H₂O, rather than to quenching of CO₂ itself.

To understand the differences between the equilibrium and kinetics simulations, we compare the advection and chemical timescales. Where the advection timescale is much smaller than the chemical timescale (τ_adv < < τ_chem) the composition should be well-mixed. Conversely, where the chemical timescale is much smaller than the advection timescale (τ_adv > > τ_chem), local chemical equilibrium should hold. Where the two timescales are comparable (τ_adv ≈ τ_chem) both processes are important to consider and the situation is more complex.

We estimate the zonal, meridional and vertical advection timescales as ${\tau }_{\text{adv}}^u\approx \frac{2\pi R_{\text{P}}}{\left|u\right|}$, ${\tau }_{\text{adv}}^v\approx \frac{\pi R_{\text{P}}}{2\left|v\right|}$, and ${\tau }_{\text{adv}}^w\approx \frac{H}{\left|w\right|}$, respectively, where R_P is the planet radius, H is the vertical pressure scale height and u, v, and w are the zonal, meridional, and vertical components of the wind velocity. We sample the wind velocities directly from the 3D grid to capture the spatial variation of the advection timescale. We note that these can be considered only as zeroth-order estimates.

We estimate the chemical timescale using the methods described in Tsai et al. (2018). Since CH₄ shows a particularly interesting effect from advection, and because it is an important absorbing species, we focus on that molecule for this timescale comparison. However, we note that the chemical timescale, for a given pressure and temperature, varies significantly for different molecules (Tsai et al. 2018). The chemical timescale for CH₄ that we use in this analysis has been specifically calculated for the Venot et al. (2019) chemical network. In practice, we calculated the CH₄ chemical timescale for a grid of pressures and temperatures. We then interpolated these timescales on to the model 3D grid using the pressureand temperature of each grid cell.

Figure 6 shows the ratio of the chemical to advection timescales for CH₄ (α = τ_chem∕τ_adv) for each of the three components of the 3D wind (zonal, meridional, and vertical) for the kinetics simulation of HD 209458b. For the zonal component we show the timescales ratio as a function of pressure and longitude for a meridional mean between ± 30°, whereas for the meridional and vertical components we show the timescales ratio as a function of pressure and latitude for a pole-to-pole slice centred on the anti-stellar point (i.e. a longitude of ϕ =0°). Where α <1 (red colours in Fig. 6) the chemical timescale is smaller (faster) than the advection timescale and the chemical equilibrium is expected to hold. However, where α > 1 (blue colours in Fig. 6) the chemical timescale is larger (slower) than the advection timescale and the atmosphere is expected to be well-mixed (in the direction of the considered wind velocity). We reiterate that both the chemical and advection timescales are only estimates, which conservatively should be seen as zeroth-order estimates.

For thezonal component of the timescale ratio (top panel, Fig. 6), α decreases with increasing pressure, which is primarily due to the pressure and temperature dependence of the chemical timescale. Following the contour of α = 1, where the timescales are comparable, it is clear that there is a significant dependence of α on longitude. For a line of constant pressure, in between ~10³ and ~10⁴ Pa, α >1 on the nightside whereas α < 1 for a significant portion of the dayside, largely due to the dependence of the chemical timescale on the local atmospheric temperature. This has the consequence that the atmosphere is expected to hold local chemical equilibrium on the dayside, but not on the nightside in the presence of zonal advection. This relates to the idea of “contamination” of the nightside by the chemistry of the dayside as discussed by previous authors (e.g. Agúndez et al. 2014a). The dayside atmosphere is hot enough to remain in chemical equilibrium and the fast zonal winds can transport this material to the cooler nightside that is not hot enough to maintain chemical equilibrium.

This structure of the timescale ratios helps to explain vertical profiles of the equatorial abundance of CH₄ in Fig. 5 (top left panel) where we see its nightside abundances departing from chemical equilibrium and tending towards its equilibrium abundance in the hottest part of the dayside. Essentially this is horizontal quenching of CH₄ (Agúndez et al. 2014a). We note that the pressure range of this feature does not exactly match what we would expect from the timescale comparison shown in Fig. 6. We attribute this simply to the fact that these timescale estimates are only zeroth-order estimates that we do not expect to precisely match the results of the full numerical simulations. We use the estimated timescales only as an aid to understand the complex 3D chemical structure that results from 3D advection.

For the meridional component of the timescale ratio (middle panel, Fig. 6), the contour of α = 1 shows a significant latitudinal dependence and lies at lower pressures nearer to the equator. Again, this is primarily due to the temperature dependence of the chemical timescale. Near the equator the atmosphere is warmer, which results in a faster chemical timescale. The spatial variation of the meridional wind velocity also plays a role. From this figure we can see that for some pressure range the atmosphere is expected to be in local chemical equilibrium in the equatorial region, but not for higher-latitudes. This has the important consequence that, for some pressure levels, the equatorial region can become chemically isolated from the mid-latitudes and polar regions. This can explain why our results show different compositions in the equatorial region compared with higher latitudes on the P = 10⁴ Pa isobar for CH₄, and also for NH₃, and HCN, as shown in Fig. 4.

The vertical component of the timescale ratio (bottom panel Fig. 6) is shown as a latitudinal slice at a longitude of 0°. The structure of the vertical component of the timescale ratio is very similar to that of the meridional component. The contour of α = 1 lies at lower pressures in the equatorial region, where temperatures are higher, compared with higher latitudes. The transition lies at much lower pressures on the warmer dayside, in fact almost at the upper model boundary (see Fig. E.1). Because of this zonal variation in the ratio of the chemical to vertical advection timescale, we only expect vertical mixing to be important for the cooler nightside region. However, because the zonal advection timescale becomes faster than the chemical timescale at higher pressures, once a species becomes vertically homogenised on the nightside, that vertically quenched abundance can be subsequently transported zonally, leading to both horizontal and vertical homogenisation all around the planet. This explains the vertically uniform equatorial abundance of CH₄ shown in Fig. 5.

Fig. 4 Mole fractions (colour scale and white contours) of CH4, H2O, NH3, CO, CO2 and HCN and wind vectors (arrows) on the 1 × 104 Pa (0.1 bar) isobar for the equilibrium simulation (left column) and the kinetics simulation (right column) of HD 209458b.

Fig. 4 continued.

Fig. 5 Vertical mole fraction profiles for different longitude points around the equator, for the equilibrium simulation (dashed lines) and kinetics simulation (solid lines) of HD 209458b. Warmer colours indicate a profile that is closer to the substellar point, cooler colours are closer to the antistellar point.

Fig. 6 Ratios of the chemical to advection timescales for CH4 (α = τchem∕τadv), for the zonal (top), meridional (middle), and vertical (bottom) components of the wind for the kinetics simulation of HD 209458b. We show a meridional-mean between ± 30° latitude for the zonal component, and slices at a longitude of 0° (i.e. the antistellar point) for the meridional and vertical components. The colour scale shows α, with blue colours corresponding to α > 1 and red colours corresponding to α < 1. The black contours show the α = 0.1, 1, 10 values. The fine structure shown in the plots is a reflection of regions where the wind velocity changes signs (e.g. eastward to westward) and the wind velocity approaches a value of zero over a very small region at the transition.

3.2.2 HD 189733b

We now consider the simulations of the atmosphere of HD 189733b which, though overall cooler, shows a similar qualitative thermal structure and circulation to HD 209458b. Figure 7 (left column) shows the mole fractions of CH₄, CO, H₂O, CO₂, NH₃, and HCN on the P = 10⁴ Pa isobar for the equilibrium simulation of HD 189733b. The distribution of abundances are generally very similar to those for the warmer HD 209458b (Fig. 7) but with quantitative differences due to differences in the temperature (Fig. 2). Since the atmosphere is cooler than HD 209458b there is an overall increased relative abundance of CH₄ by one to two orders of magnitude, which is favoured thermodynamically for cooler temperatures. In particular, in the nightside mid-latitude regions the temperature is cold enough at P = 10⁴ Pa for CH₄ to become the most abundant carbon species, replacing CO.

Figure 7 (right column) shows the mole fractions of CH₄, CO, H₂O, CO₂, NH₃, and HCN on the P = 10⁴ Pa isobar for the kinetics simulation of HD 189733b. Here we see that for most molecules the effect of 3D advection is to largely remove the horizontal abundance gradients that are present in the chemical equilibrium simulation. For CH₄, CO, H₂O, HCN, and NH₃ the remaining horizontal abundance gradients are very small. However, it is clear that CO₂ does not change significantly from its chemical equilibrium abundance on the dayside. This means that temperature-dependent kinetics is still important for this molecule on the warm dayside and advection is not the dominant process in determining its abundance. Larger differences in the CO₂ abundance between the equilibrium and kinetics simulations are found on the cooler nightside. For CH₄, CO, and H₂O our results qualitatively agree with our earlier work using a more simplified chemical scheme (Drummond et al. 2018a), as we also find that these molecules are efficiently horizontally mixed.

Figure 8 shows both the equilibrium and kinetics abundances of CH₄, CO, H₂O, CO₂, NH₃, and HCN, as vertical (pressure) mole fraction profiles around the equator. This also shows that most of the molecules are zonally and vertically well-mixed over a large pressure range at the equator. Compared with HD 209458b (see Fig. 5), the molecules typically quench at higher pressures, which is not surprising since the temperature is lower and the chemical timescale is longer.

The most obvious difference between our results for HD 189733b and HD 209458b are in the vertical profiles of CH₄. For HD 209458b, over a significant pressure range the abundance of CH₄ appears to be zonally but not vertically quenched, with vertical mixing becoming important only for pressures less than 10⁴ Pa. For HD 189733b on the other hand we find that CH₄ becomes vertically and horizontally well-mixed (i.e. quenched) much deeper, approximately at the base of the equatorial jet (Fig. 1). In the pressure range 10⁵ Pa to 10⁴ Pa the equatorial abundance of CH₄ becomes larger than any of its equatorial equilibrium abundances, significantly increasing the relative abundance of CH₄ for this region of the atmosphere. This structure is very similar to that found previously using a more simplified chemical scheme (Drummond et al. 2018a,b). In those earlier works we argued that this feature is due to meridional transport of material from the cooler mid-latitudes, where CH₄ is more abundant at local chemical equilibrium, to the equatorial region.

In Drummond et al. (2018b) we conducted a simple tracer experiment that demonstrated that mass can be transported from the mid-latitudes into the equatorial region on a timescale of several hundred days. The tracer was initialised to an abundance of zero everywhere, except in a “source” region at high-pressure and high-latitudes. After several hundred days of model integration, the tracer had increased in abundance significantly at the equator with a clear equator-ward transport occurring between 10⁵ Pa and 10³ Pa. Once the material is transported into the equatorial region, the material is rapidly distributed zonally and vertically by the fast equatorial jet and large dayside upwelling. Since the circulation for both HD 209458b and HD 189733b is broadly similar between this work and that tracer experiment, we expect the same transport process to be occuring. The difference between the two cases is in the chemical timescale, with a faster chemical timescale leading to a lower pressure quenching in the case of HD 209458b, above the region where the meridional transport is important.

The abundance profiles of CO₂ are quite similar to those found for HD 209458b. CO₂ remains close to its chemical equilibrium value to much lower pressures than other molecules shown. For P > 10³ Pa, we attribute the departure of the abundances of CO₂ from chemical equilibrium as being due to its pseudo-equilibrium with CO and H₂O. But for P ≲ 10³ Pa, where CO₂ becomes approximately zonally uniform we expect this to be due to zonal quenching of CO₂ itself.

Figure 9 shows the ratio of the chemical to advection timescales for CH₄ for the kinetics simulation of HD 189733b. Generally, the point at which the chemical and advection timescales become comparable (black line) lies at higher pressures than for HD 209458b. In addition, the pressure level at which the advection and chemical timescales become comparable is much more similar between the three components (zonal, meridional and vertical) compared with for HD 209458b, with the transition occurring at P ~ 10⁵ Pa. Together this means that we expect the molecules to quench deeper, and for zonal, meridional and vertical quenching to act in approximately the same regions. This explains why we see such different behaviour for CH₄ between HD 209458b and HD 189733b.

Fig. 7 Mole fractions (colour scale and white contours) of CH4, H2O, NH3, CO, CO2, and HCN and wind vectors (arrows) on the 1 × 104 Pa (0.1 bar) isobar for the equilibrium simulation (left column) and the kinetics simulation (right column) of HD 189733b.

Fig. 7 continued.

3.3 Chemical-radiative feedback

The advected chemical abundances of CH₄, CO, CO₂, H₂O, NH₃, and HCN are used to calculate the total opacity in each grid cell which is then used to calculate the radiative heating rates. This allows for a chemical-radiative feedback process as changes in the chemical composition, due to the transport-induced quenching, can affect the thermal structure and circulation. Of course, the chemistry itself is strongly dependent on the temperature via the temperature-dependent rate coefficients.

Feedback between transport-induced quenching of chemical species and the temperature profile has previously been investigated in 1D considering only vertical diffusion (Drummond et al. 2016), showing temperature changes of up to 100 K. In 3D, Drummond et al. (2018b) found that advection of CH₄ can lead to changes in the temperature by up to 5–10% for HD 189733b, compared with the chemical equilibrium case. Steinrueck et al. (2019) artificially varied the CO/CH₄ ratio in their 3D model and found that significant temperature responses can result, compared with the assumption of local chemical equilibrium.

For HD 209458b we find negligible differences (<1%) in the atmospheric temperature and wind velocities between the equilibrium and kinetics simulations. This is in agreement with our earlier work using a simplified chemical relaxation scheme (Drummond et al. 2018a). While there are significant changes in the relative abundances of several molecules, as shown in the previous section, the abundances of these molecules remain relatively small and therefore make only a minor contribution to the total opacity.

On the other hand, for HD 189733b we find more significant differences (5− 10%) in the temperatures and wind velocities between the equilibrium and kinetics simulations. The results are very similar to those presented in Drummond et al. (2018b) for the same planet but using a simplified chemical relaxation scheme (only considering CH₄, CO, and H₂O) and we therefore do not show the results here. The fact that the results are so similar is not overly surprising since in Drummond et al. (2018b) we found that the changes in the temperature and winds were mainly due to the changes in the CH₄ abundance, and the CH₄ abundance profiles for HD 189733b are very similar between this work and Drummond et al. (2018b). This suggests that quenching of other chemical species considered in this work is not important for the radiative heating rates.

Fig. 8 As Fig. 5, but for the equilibrium simulation (dashed lines) and kinetics simulation (solid lines) of HD 189733b.

3.4 Spectra and phase curves

In this section we present synthetic spectra calculated directly from our 3D atmosphere simulations. We calculate both transmission spectra and secondary eclipse emission spectra, and search for signatures of the chemical changes between the equilibrium and kinetics simulations. We also present simulated JWST and ARIEL observations to assess the observational significance of any differences in the model spectra from the perspective of these up-coming space-based instruments.

PandExo (Batalha et al. 2017) was used to simulate JWST observations of the emission and transmission spectra. The simulated instruments used were the Near InfraRed Spectrograph (NIRSpec) and the Mid-Infrared Instrument (MIRI), using their G395H and LRS modes respectively. The stellar parameters were sourced from the TEPCat database (Southworth 2011). The results were computed assuming a single orbit/transit, and an equal amount of observing time in-and-out of transit, corresponding to a total observation time of 6.096 hours for HD 209458 and 3.614 hours for HD 189733.The ARIEL Noise Simulator will be described in a future work (Mugnai et al., in prep.). For the ARIEL simulated observations five eclipses where used for the emission spectra and ten transits for transmission spectra.

Fig. 9 As Fig. 6, but for the kinetics simulation of HD 189733b.

3.4.1 HD 209458b

Figure 10 shows the emission and transmission spectra calculated for both the equilibrium and kineticssimulations of HD 209458b. In emission it is clear that there are negligible differences between thetwo cases. This is not surprising since the most abundant absorbing species (CO and H₂O) show relatively little difference in abundance between the equilibrium and kinetics simulations. Whilst the abundances of CH₄, NH₃, and HCN change significantly, their relative abundances are typically below ~ 10⁻⁷ and therefore make only a small contributionto the total opacity. CO₂ remains close to its equilibrium abundance on the dayside, which the secondary eclipse is probing. The emission spectrum is also very sensitive to the temperature structure, though there are only very minor (< 1%) changes in the temperature.

In transmission, however, there are two clear spectral regions where the equilibrium and kinetics simulations show differences. Firstly, at ~4.5 μm the kinetics simulation shows a reduced transit depth compared with the equilibrium simulation. This is due to a reduction in the abundance of CO₂ in the terminator regions (see Fig. 5) which decreases the opacity. In addition, there is an increase in the transit depth between 10− 12μm, which we find is due to enhanced NH₃ absorption. The NH₃ abundance is increased significantly in the kinetics simulation compared with the equilibrium simulation, which increases the transit depth in this region where NH₃ absorbs.

The middle and bottom rows of Fig. 10 overlay the simulated JWST and ARIEL observations, respectively. The JWST simulated observations assume a single transit/eclipse while the ARIEL simulated observations include five transits for transmission and ten eclipses for emission. For a single transit the NIRSpec G395H instrument is able to resolve the feature at 4.5 μm, with additional eclipses expected to improve this. Ten transits with ARIEL are also enough to resolve this feature. The feature due to NH₃ at longer wavelengths is not well resolved by the MIRI LRS simulated observations, and is not included in the spectral coverage of ARIEL.

Figure 11 shows the simulated emission phase curves for HD 209458b. Considering the secondary eclipse emission spectra previously discussed, it is unsurprising that there is relatively small difference between the equilibrium and kinetics simulations. The most notable difference between the two cases occurs when the nightside of the planet is in view, with the kinetics simulation showing an increased flux ratio compared with the equilibrium simulation. This occurs in the 3–4, 7–8 and 8–9 μm bands, where CH₄ is a prominent absorber and is almost certainly due to the decreased nightside abundance of CH₄, due to zonal quenching (Fig. 5). The decreased CH₄ abundance shifts the photosphere to higher pressures and temperatures, resulting in a larger emission flux.

3.4.2 HD 189733b

Figure 12 shows the emission and transmission spectra calculated for both the equilibrium and kinetics simulations of HD 189733b. In emission there are much more obvious differences between the equilibrium and kinetics simulations for HD 189733b, compared with for HD 209458b. This is primarily due to an enhanced absorption by CH₄ which is increased in abundance by several orders of magnitude in the kinetics simulation. Additionally there are features due to NH₃ at ~6.5 μm and 10−12 μm, where increased abundances lead to increased opacity and therefore a reduced eclipse depth. In the middle and lower panels of Fig. 12 it is clear that a single eclipse observation with JWST and five eclipses with ARIEL are sufficient to clearly differentiate between the two models.

In transmission there are very significant differences between the equilibrium and kinetics simulations across a wide range of wavelengths. The most prominent features are due to CH₄ between 3−4 μm and 7−10 μm, where an increased CH₄ abundanceleads to an enhanced transit depth. There is also a significant contribution from NH₃ which increases the transit depth at ~6.5 and 10−15 μm. The pale orange line in Fig. 12 (top right) shows a calculation that excludes absorption due to NH₃ in the spectrum calculation only (it is still included in the main model integration) for the kinetics simulation. Comparing this with the spectra for the standard kinetics simulation clearly shows the importance of NH₃ in shaping the transmission spectrum for HD 189733b, which is greatly enhanced by 3D mixing. HCN also absorbs in this wavelength region, however a similar test that removes the contribution of HCN to the transit depth shows negligible difference with the kinetics spectrum that includes the full set of absorbers. The transit depth is decreased at 4.5 μm due to a decrease in the CO₂ abundance, similarly to HD 209458b.

The simulated JWST (single transit) clearly shows the enhanced transit depths in the region where CH₄ absorbs. However, the large features due to NH₃ at longer wavelengths are lost in the noise. Ten transits with ARIEL enables the clear differentiation between the two models in the regions of the enhanced CH₄ absorption and decreased CO₂ absorption. The enhanced transit depth at ~6.5 μm is also apparent in the simulated ARIEL observations.

Figure 13 shows the simulated emission phase curves for HD 189733b. There is a significant difference between the equilibrium and kinetics simulations for many of the wavelength bands presented. In most of the wavelength bands the simulation including chemical mixing shows a decreased dayside emission (as shown previously in the secondary eclipse emission spectra) and an increase in the nightside emission. This is primarily due to increases in the abundance of CH₄, which is particularly important on the dayside. However, the large increase in the NH₃ abundance also contributes in the 6–7 μm band. Increases in the abundances of these species increases the atmospheric opacity and leads to a lower temperature emission (since the temperature decreases with altitude). The 4–5 μm channel shows very little change in the dayside emission but a small increase in the nightside emission. This is likely due to the decrease in the nightside CO₂ abundance, leading to a deeper and hotter photosphere in this wavelength region.

The main effect of chemical mixing on the emission phase curve is to decrease the amplitude (the difference between the maximum dayside emission and minimum nightside emission). However, it is also apparent that there is a slight increase in the offset of the phase curve, with the peak emission being shifted further away from secondary eclipse. Thisis consistent with our previous simulations of HD 189733b using much simpler chemical scheme and is due to the temperature response of the atmosphere (Drummond et al. 2018b). The temperature response of the atmosphere is very similar to that already presented in Drummond et al. (2018b) and in Fig. 5 of that paper we show that the largest temperature increase occurs east of the substellar point, which shifts the hotspot further east.

Fig. 10 Transmission (left) and emission (right) spectra for the equilibrium (blue) and kinetics (red) simulations of HD 209458b. Top row: model spectra, middle row: simulated JWST observations and bottom row: simulated ARIEL observations.

Fig. 11 Emission phase curves for the equilibrium (solid) and kinetics (dashed) simulations of HD 209458b in a series of spectral bands. The y-axis is the ratio of the planet to stellar flux.

Fig. 12 As Fig. 10 but for HD 189733b.

Fig. 13 As Fig. 11 but for HD 189733b.

4 Discussion

Here we discuss the results presented in this work in the context of both our own previous studies, and those from the literature, with particular focus on two key works (Cooper & Showman 2006; Agúndez et al. 2014a).

4.1 Comparison with our earlier work

In our earlier work (Drummond et al. 2018a,b) we used a much simpler and more computationally efficient chemical scheme, by revisiting the chemical relaxation scheme developed by Cooper & Showman (2006) where only the advection and chemical relaxation of CO is considered. The chemical timescale was based on a rate-limiting reaction that was identified in the interconversion steps from CO to CH₄. The mole fractions of CH₄ and H₂O were subsequently found by assuming mass balance assuming that all carbon is in the form of CO and CH₄ and all oxygen is in the form of CO and H₂O.

For HD 189733b there is a very good agreement for the abundance distributions of CH₄, H₂O, and CO found using the Cooper & Showman (2006) chemical relaxation scheme (Drummond et al. 2018b) or using the coupled chemical kinetics scheme (this work). This implies that in the particular case of HD 189733b the parameterised chemical timescale developed by Cooper & Showman (2006) gives a similar chemical evolution to the reduced chemical network of Venot et al. (2019) for CO. The results and conclusions presented in Drummond et al. (2018b) are therefore validated here with the use of a more accurate chemical scheme. This important result demonstrates that, at least in some cases, simplified chemistry schemes can be used in 3D atmosphere models, greatly increasing computational efficiency while maintaining accurate results, compared with more sophisticated chemical schemes.

For HD 209458b the picture is quite different and there are significant differences in the abundance distributions of CH₄ between simulations using the chemical relaxation scheme and simulations using the coupled chemical kinetics scheme. Specifically, CH₄ remains in chemical equilibrium at lower pressures when using the coupled chemical kinetics scheme rather than the chemical relaxation scheme (Drummond et al. 2018a). This has the important consequence that when CH₄ does undergo mixing, the quench point lies in a different region of the atmosphere where the circulation is different, resulting in a rather different CH₄ distribution.

In Cooper & Showman (2006), the chemical timescale represents the time for the chemical interconversion of two molecules, namely CO and CH₄, to take place. This, combined with the fact that CH₄ is found by mass balance with CO, means that both of these molecules have the same chemical timescale by definition. In practice, this entails that both molecules quench at the same location in the atmosphere. In contrast, in the chemical relaxation scheme more recently developed by Tsai et al. (2018) the chemical timescale represents the time taken for a particular species to evolve back to its local chemical equilibrium, following a perturbation.

In Fig. 14 we show the chemical timescales for CH₄ and CO, as a function of pressure and temperature, calculated specifically for the Venot et al. (2019) network, using the same method as described in Tsai et al. (2018). Comparing the timescales for CH₄ and CO it is apparent that for some pressure-temperature regimes the timescales for these two molecules can be significantly different (see also Tsai et al. 2018). For example, at P = 10⁵ Pa and T = 1200 K the timescale for CH₄ is log₁₀(τ[s]) ~ 10 while for CO it is log₁₀(τ[s]) ~ 12.5, more than two orders of magnitude difference. This suggests that the assumption of a single timescale describing the quenching of both CO and CH₄ as used in previous work is not accurate (Cooper & Showman 2006; Drummond et al. 2018a,b). Cooper & Showman (2006) present their CO–CH₄ interconversion timescale in their Table 2. Comparing our CO and CH₄ chemical timescales to this at P = 10⁵ Pa and T = 1500 K we find a good agreement for CO (~10⁷ s) but a larger discrepancy for CH₄ (~ 10⁵ s for our CH₄ chemical timescale compared with ~10⁷ s for the Cooper & Showman (2006) CO-CH₄ interconversion timescale). This is expected since the Cooper & Showman (2006) interconversion timescale is based on the reduction of CO.

To be clear, we are not suggesting that the CO-CH₄ interconversion timescale of Cooper & Showman (2006) is incorrect. However, the above discussion does imply that it is not appropriate to use the CO–CH₄ interconversion timescale to find the quench point of CH₄. Finally, this explains why we find the CH₄ quench point to lie at lower pressures in this work compared with our previous results (Drummond et al. 2018a,b) using the Cooper & Showman (2006) interconversion timescale.

Fig. 14 Chemical timescales of CH4 (top) and CO (bottom) estimated for the Venot et al. (2019) chemical network as a function of pressure (in bar, y-axis) and temperature (in K, x-axis). The quantity shown is the common logarithm of the timescale in seconds. The colourscales are different between the two panels.

4.2 Comparison with Cooper & Showman (2006)

Cooper & Showman (2006) were the first to investigate 3D transport of chemical species in hot Jupiter atmospheres, using the simple chemical relaxation scheme described above when simulating the atmosphere of HD 209458b. The main conclusion of this work was that chemical equilibrium does not hold for pressures less than ~ 1 bar, leading to a significant increase in CO on the cold nightside, at the expense of CH₄. Based on their simulations, the authors argued that it is predominantly vertical mixing that drives the disequilibrium of CO and CH₄. This is due to vertical quenching from a deep atmospheric level that is horizontally uniform and very abundant in CO (and therefore not very abundant in CH₄). This vertical quenching occurs all across the planet, leading to horizontally uniform abundances in the photosphere.

In Drummond et al. (2018a) we compared our results, from a different 3D atmospheric dynamics model but using the same simple chemical relaxation scheme, with those of Cooper & Showman (2006) and found important quantitative and qualitative differences. Firstly, we found a significantly different chemical equilibrium composition, with CO being the dominant carbon-species everywhere (except for very high pressures), as opposed to the CO-rich dayside and CH₄-rich nightside atmosphere found by Cooper & Showman (2006). We also found that horizontal (specifically, meridional) quenching, in combination with vertical quenching, is important in determining the photospheric abundances of CH₄. We attributed these differences to being almost entirely due to differences in the way the atmospheric heating is taken into account in the model, leading to a very different thermal and chemical structure. Cooper & Showman (2006) used a simple temperatureforcing (or Newtonian cooling) scheme whereas Drummond et al. (2018a) used a coupled radiative transfer scheme to calculate the radiative heating rates consistently with the chemical composition. The temperature forcing scheme developed for HD 209458b significantly overpredicts the day-night temperature contrast with a dayside that is too hot and a nightside that is too cold, compared with the coupled radiative transfer scheme (e.g. Showman et al. 2009; Amundsen et al. 2016).

4.3 Comparison with Agúndez et al. (2014a)

Agúndez et al. (2014a) developed a model that includes a sophisticated chemistry scheme but makes a very crude approximation to horizontal chemical transport. The Agúndez et al. (2014a) model is in practice a 1D vertical column model with a time-varying pressure-temperature profile sourced from previous 3D atmosphere model output for different longitudes. This so-called “pseudo-2D” model represents therefore a column of atmosphere rotating around the equator (i.e. within the equatorial zonal jet) with a solid-body rotation. Vertical diffusion was also included and parameterised with an eddy diffusion coefficient. In terms of the chemistry, Agúndez et al. (2014a) included a chemical kinetics scheme and used the 105 species chemical network of Venot et al. (2012). The authors also included photochemical dissociation. Overall, compared with other works investigating horizontal chemical mixing in hot Jupiter atmospheres (Cooper & Showman 2006; Drummond et al. 2018a,b; Mendonça et al. 2018, and this work) the model developed by Agúndez et al. (2014a) includes a much more sophisticated approach to chemistry, and is so far the only one to consider photochemistry, but takes a much cruder approach to horizontal transport, and only considers the zonal direction.

The main conclusions from Agúndez et al. (2014a) are that horizontal quenching can occur for certain molecules leading to their abundances being homogenised with longitude, at values typical of their chemical equilibrium values on the dayside. This was found to have a particularly important effect on the abundance of CH₄, which is significantly decreased in their pseudo-2D model on the nightside relatively its local chemical equilibrium abundance. This is due to quenching from the dayside where CH₄ is much lower in abundance in chemical equilibrium. This agrees with our results for HD 209458b where we find zonal quenching of CH₄ over the pressure range 0.1 bar to 0.01 bar, decreasing the nightside abundance. However, for HD 189733b we find that deeper meridional transport has a more important impact on the CH₄ abundance, with equatorward transport of CH₄ from mid-latitudes increasing the equatorial abundance. Such a process would not be captured by the purely zonal-vertical pseudo-2D model of Agúndez et al. (2014a). We note that the HD 209458b model presented by Agúndez et al. (2014a) includes a thermal inversion which is likely to lead to differences in the chemical structure and transport, in addition to differences that result from the very different modelling approach.

Agúndez et al. (2014a) include photochemical dissociation in their model, which is found to have some impact on molecular abundances in the upper atmosphere. One of the principle effects is the increase the abundance of HCN. However, these photochemical effects are restricted to pressures less than ~ 1 Pa, above the upper domain of our model. Therefore, even if we were to include photochemical dissociation in our 3D model the effects on the chemical abundances in the pressure region covered by our domain are likely to be small.

Our results qualitatively agree with the conclusions of Agúndez et al. (2014a) for the case of HD 209458b. Zonal quenching appears to “contaminate” the cooler nightside with the composition from the warmer dayside. In this case, it may be possible that a 2D model only accounting for zonal and vertical transport can be sufficient to accurately model the equatorial regions of hot Jupiters. However, for HD 189733b we find that meridional transport, occurring deeper in the atmosphere, has a more important role in setting the quenched abundances, especially for CH₄. In this case, a full 3D model is required to capture all of the important mixing processes.

4.4 Comparison with other works

Several other recent works have also tackled the issue of chemical mixing in hot Jupiter atmospheres. Mendonça et al. (2018) coupled a chemical relaxation scheme to the THOR general circulation model to investigate the effect of 3D mixing in the atmosphere of WASP-43b. The chemical relaxation scheme was recently developed in Tsai et al. (2018) and includes more chemical species and a more accurate chemical timescale calculation than the earlier method of Cooper & Showman (2006). Mendonça et al. (2018) conclude that zonal quenching is the dominant process in the equatorial region of WASP-43b, with vertical and meridional mixing having only a secondary effect.

The results of Mendonça et al. (2018) are qualitatively consistent with our findings for HD 209458b but different to our findings for HD 189733b. For HD 209458b we find that not all horizontal abundance gradients are removed, in agreement with the conclusions of Mendonça et al. (2018) for WASP-43b, while for HD 189733b we find that the chemistry isvery efficiently horizontally homogenised. However, given that Mendonça et al. (2018) simulate the atmosphere of a different planet, with a significantly different circulation and thermal structure, it is not surprising that we should find differences. WASP-43b has a similar planetary equilibrium temperature to HD 209458b (~ 1450 K), but is considerably warmer than HD 189733b (~1200 K). The circulation and temperature structure is also likely to be quite different due to the much higher surface gravity (~ 55 m s⁻²) and slower rotation rate (~ 1.5 × 10⁻⁵ s⁻¹), compared with both HD 209458b and HD 189733b (see Table 1). In addition to these physical differences in the planet-stellar properties, there are many differences in the model implementation, level of assumption, and included model physics between our model and that of Mendonça et al. (2018). For example, Mendonça et al. (2018) use a relatively simple 2-band radiative transfer scheme compared with our non-grey radiative transfer scheme. Also, Mendonça et al. (2018) include the effect of additional cloud opacity on the nightside, while we choose to assume a cloud-free atmosphere. The very large number of differences, both in the details of the numerical model and also in the properties of the planetary system under investigation, means that a meaningful comparison is not possible.

One of the more interesting features discussed in Mendonça et al. (2018) is the presence of chemically distinct polar vortexes. We do not find similar features in our simulations, however this may again be simply due to the use of different planetary and stellar parameters associated with different planets. A more controlled comparison of our model with that of Mendonça et al. (2018), using model parameters and a setup that are as close as feasibly possible would be an interesting future exercise.

Steinrueck et al. (2019) recently investigated the feedback effect of non-equilibrium chemical abundances on the 3D circulation and thermal structure of HD 189733b. Assuming globally uniform chemical abundances (i.e. no interactive chemistry scheme) the authors varied the ratio of CO and CH₄. The assumption of globally uniform abundances for HD 189733b is likely a good one based on our simulation results presented here and in Drummond et al. (2018a), with both of these molecules having roughly uniform abundances for pressures less than ~ 10⁵ Pa.

Steinrueck et al. (2019) found that in cases where CO remains the dominant carbon-species, varying the CO/CH₄ ratio can alter the temperature by 50–100 K. In cases where CH₄ is more abundant than CO (a less likely scenario) a larger temperature response (200–400 K) is found for changes in the CO/CH₄ ratio. In the more likely CO-dominated regime, our results using a more sophisticated approach to the chemistry support the conclusions of Steinrueck et al. (2019) that disequilibrium chemical abundances can alter the thermal structure by up to 10% compared to abundances predicted assuming local chemical equilibrium. For the equatorial region of HD 189733b, our simulations predict that HD 189733b remains in the CO-dominated regime, but that 3D vertical mixing can decrease the CO/CH₄ ratio by several orders of magnitude, due to an increase the equatorial abundance of CH₄. This leads to temperature changes of 10% compared with the local chemical equilibrium simulation (Drummond et al. 2018b).

5 Conclusions

In this paper we presented results from a model that consistently couples a chemical kinetics scheme to a 3D radiation-hydrodynamics atmosphere model. We use the recently released 30 species carbon-oxygen-nitrogen-hydrogen chemical network of Venot et al. (2019). We investigated the effect of 3D advection on the chemical structure of hot Jupiter atmospheres, using the particular cases of HD 209458b and HD 189733b. Focusing on six important absorbing species (CH₄, CO, H₂O, CO₂, NH₃, and HCN) we found that 3D advection can lead to significant changes in the chemical composition compared with a simulation that assumes local chemical equilibrium. These chemical differences lead to important signatures in simulated emission and transmission spectra, primarily due to CH₄, CO₂, and NH₃. The model presented in this work represents an improvement on our earlier studies (Drummond et al. 2018a,b) which used a much simpler chemical relaxation scheme (that of Cooper & Showman 2006), while using essentially the same radiative transfer and hydrodynamics setup.

For HD 189733b, interestingly, we find qualitatively very similar results to our earlier work (Drummond et al. 2018b). The abundance profiles of CH₄, CO, and H₂O (the only three molecules considered in that earlier work) are very similar. The most important consequence of 3D mixing is that the abundance of CH₄ is significantly increased, particularly in the low pressure regions of the atmosphere, in comparison to the chemical equilibrium case, leading to a decrease in the eclipse depth and an increase in the transit depth, where CH₄ absorbs. Our new results that include more chemical species show that CO₂ and NH₃ also make important contributions to both the emission and transmission spectra, and that the abundance profiles of these species are strongly effected by 3D mixing, particularly NH₃. We therefore conclude that CH₄ and NH₃ are important targets for future observations of HD 189733b as they may help to constrain the chemical-dynamical interactions occurring in HD 189733b-like atmospheres.

For HD 209458b we find qualitatively quite different results when using the newly coupled chemical kinetics scheme, compared with our earlier results using a much simpler chemical relaxation scheme (Drummond et al. 2018a). In Drummond et al. (2018a) we found that CH₄ undergoes a quenching behaviour similar to that which we found for HD 189733b and which significantly increases the abundance of CH₄, at low pressures, in comparison to the chemical equilibrium case. Here, however, we find that CH₄ remains in equilibrium to lower pressures before quenching occurs, leading to a much smaller quenched abundance of CH₄. This has the result that we no longer find the significant effects of increased CH₄ absorption in the emission and transmission spectra as previously found in Drummond et al. (2018a).

An important conclusion of this work is that understanding the chemical structure of hot Jupiter atmospheres is not trivial, and differences in the winds or the temperature structure can have large impacts on the predicted chemical structure. There is likely to be a huge diversity of exoplanetary atmospheres, showing much more variety in the atmospheric circulation patterns and temperature structures than captured by the study of just two specific hot Jupiters, as featured in this work. The consequences of this on the predicted chemical composition would require a large parameter space study, which would be difficult to achieve with such computationally expensive models, requiring a step–change in the efficiency of the model framework.

Another important conclusion from this work is that the advection processes that determine the chemical structure are strongly three-dimensional, and these processes cannot be captured in a 1D, or even 2D, model. For HD 209458b our kinetics simulation shows that in some cases the chemistry can become efficiently zonally mixed but is not strongly effected by vertical mixing. This occurs because the horizontal advection timescale is typically a few orders of magnitude faster than the vertical advection timescale. At lower pressures, when vertical mixing eventually becomes important, the vertically quenched abundance is effectively set by the zonally quenched abundance that is determined at higher pressures. The note of caution, however, is that precise absolute wind velocities can not be easily predicted by such models (see discussion in Heng et al. 2011), and require observational constraint (e.g. Snellen et al. 2010; Louden & Wheatley 2015).

The quantitative wind velocities, and therefore the advection timescales, are partly determined by the dissipation included in the model (required for the numerical stability) the magnitude of which is currently unconstrained by observations Heng et al. (2011). Adjusting this model free parameter could lead to changes in the quenching behaviour of chemical species, as the balance of advection and chemistry will change. However, we do not expect that changes of a factor of a few in the wind velocities (and therefore also in the advection timescales) will cause a very large change in the quenching behaviour since the chemical timescale that depends on the temperature varies much more steeply spatially in the atmosphere. Differences could also arise between models based on the chosen chemical network. Different choices of reactions or their associated rate constants can lead to a variety of quenching points and quenched abundances for a given circulation (e.g. Venot et al. 2012; Moses 2014). The precise chemical structure may vary from one model to another, with different quench points and different quenched abundances, due to differences in the predicted wind and temperature fields and also the use of different chemical networks. An intercomparison of models will be required in future to validate our results, once a family of models with similar capability are available in the community.

Recently, Sainsbury-Martinez et al. (2019) confirmed, using a 3D dynamical atmosphere model, that the deep atmosphere is much hotterthan what is predicted by 1D radiative-convective equilibrium models, due to vertical advection of potential temperature(Tremblin et al. 2017). Integration times of much longer than 1000 (Earth) days, as used here, are required to achieve a steady-state in the deep atmosphere, when initialised from a “cold” initial profile. However, Sainsbury-Martinez et al. (2019) demonstrated that initialising the model with a hotter initial profile can reduce the integration time required for steady-state to be reached. Lines et al. (2018b), who coupled a cloud kinetics formation scheme to a 3D radiation-hydrodynamics atmosphere model, has shown that initialising the atmosphere with a “hot” initial profile can lead to significantly different cloud structures in the upper atmosphere, due to differences in the temperature profile at P > 10⁵ Pa. We performed a similar test for our HD 209458b case (see Appendix C) but found no significant differences in the atmospheric structure or composition between the hot and cold initial profile models. This is likely because most species quench at pressures P ≲ 10⁵ Pa while the temperature profile is only different for pressures greater than this. Since the chemistry quenches deeper in the case of HD 189733b it is possible that considering the hot deep atmosphere may lead to more significant differences in the upper atmosphere chemistry, if the quench points are deep enough. If the deep temperature profile has not yet reached a steady-state at the pressure level of quenching, then a wrong quenched abundance will be found, affecting the abundance at lower pressures. This potential issue should be considered when interpreting our results, and should be investigated in more detail in future by running the model for a much longer integration time and/or initialising with a hotter initial thermal profile Sainsbury-Martinez et al. (2019).

The model described in this paper represents a major step forward in the ability to study the interactions of chemistry, radiation and fluid dynamics in a 3D exoplanet atmosphere. However, there are a number of important physical elements that are clearly missing that are likely to have a significant impact. Firstly, though the model takes into account 3D transport of chemical species by the resolved winds we do not yet account for photochemical processes, as included routinely in 1D chemistry models (e.g. Moses et al. 2011; Venot et al. 2012). Photodissociation and photoionisation have the potential to significantly change the atmospheric composition, particularly at high altitudes (low pressures). This is not only due to the direct effect of additional loss terms for molecules and atoms via dissociation and ionisation. The products of photodissociation and photoionisation are often very reactive species themselves, which in turn impact the chemical balance even further. Photochemical products formed on the irradiated dayside may well be efficiently transported to the cooler nightside, exciting additional chemical pathways. However, 1D models typically show that photochemical effects only become important for pressures less than around  1 Pa (e.g. Moses et al. 2011; Drummond et al. 2016) which is typically the pressure of the upper boundary of 3D simulations. In order to study photochemical effects in a 3D atmosphere, the model may first have to be extended to lower pressures which is would require significant development and testing.

A further major limitation of this work is the absence of condensates and clouds in the simulations. Here we refer to condensate clouds rather than photochemical hazes. In this work we consider a gas-phase only composition. Condensates and clouds of various compositions are expected to form in hot Jupiter atmospheres and there is a wide-variety of approaches with varying complexity to model them in both 1D (e.g. Barstow et al. 2017; Morley et al. 2017; Goyal et al. 2018; Powell et al. 2019) and 3D models (e.g. Charnay et al. 2015; Lee et al. 2016; Lines et al. 2018b,a, 2019). The presence, distribution and composition of cloud can strongly effect the radiation field, which determines the thermal structure and observable properties (Lines et al. 2018a, 2019). The formation of cloud can also have an important effect on the gas-phase composition, as the gas becomes depleted in the elements that form the cloud condensates themselves (see discussion in Lee et al. 2015; Woitke et al. 2018; Lines et al. 2018b). In this work we focus on understanding the gas-phase chemistry which, fundamentally, is the necessary precursor to any cloud formation. However, fully understanding the composition of hot Jupiter atmospheres likely requires coupling gas-phase chemical kinetics and cloud kinetics formation schemes within the framework of a 3D model. This represents another very significant scientific and technical challenge.

Acknowledgements

The authors thank the anonymous referee for a very useful and constructive report that led to improvements to the manuscript. B.D., E.H. and N.J.M. acknowledge support from a Science and Technology Facilities Council Consolidated Grant (ST/R000395/1). N.J.M. was part funded by a Leverhulme Trust Research Project Grant during this project. This work also benefited from the 2018 Exoplanet Summer Program in the Other Worlds Laboratory (OWL) at the University of California, Santa Cruz, a program funded by the Heising-Simons Foundation. O.V. thanks the CNRS/INSU Programme National de Planétologie (PNP) and the Centre National d’Études Spatiales (CNES) for funding support. Q.C. acknowledges funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 758892, ExoAI), under the European Union’s Seventh Framework Programme (FP7/2007-2013)/ ERC grant agreement numbers 617119 (ExoLights) and the European Union’s Horizon 2020 COMPET programme (grant agreement No 776403, ExoplANETS A) and also by the Science and Technology Funding Council (STFC) grants: ST/K502406/1, ST/P000282/1, ST/P002153/1 and ST/S002634/1. J.M. acknowledges the support of a Met Office Academic Partnership secondment. P.T. acknowledges supports by the European Research Council under Grant Agreement ATMO 757 858. This work was performed using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which formspart of the STFC DiRAC HPC Facility (www.dirac.ac.uk). The equipment was funded by BEIS capital funding via STFC capital grants ST/K000373/1 and ST/R002363/1 and STFC DiRAC Operations grant ST/R001014/1. DiRAC is part of the National e-Infrastructure. This research also made use of the ISCA High Performance Computing Service at the University of Exeter. Additionally, we acknowledge use of the Monsoon system, a collaborative facility supplied under the Joint Weather and Climate Research Programme, a strategic partnership between the Met Office and the Natural Environment Research Council. Material produced using Met Office Software. The research data supporting this publication are available upon request to the authors.

Appendix A Model evolution and steady-state

Here we assess the model evolution in order to judge whether a steady-state has been reached. Figure A.1 shows the evolution of the globalmaximum in each of the three components (zonal, meridional and vertical) of the wind velocities. After starting from rest (0 km s⁻¹) the winds quickly spin up to several km s⁻¹ over a timescale of a few hundred days. By 1000 days, the maximum integration time of the simulations presented here, the acceleration of the global maximum winds is very small. The zonal wind component for HD 189733b shows a small gradient at 1000 days.

Figure A.2 shows the evolution of the mole fractions of CH₄ and HCN in a series of arbitrarily chosen grid cells for the kinetics simulation of HD 209458b. The grid cells correspond to three model vertical (altitude) levels (~10² Pa, ~ 10³ Pa, and ~ 10⁴ Pa) at two arbitrary coordinates (λ = ϕ = 0° and λ = 180°, ϕ =45°). The legends refer to the pressures of the corresponding grid cells at the end of the simulation (1000 days). Clearly, most of the changes in the mole fractions, from their initial values, occurs within the first few hundred days, reflecting somewhat the evolution of the global maximum wind velocities (Fig. A.1). At the end of the simulation the gradients are approximately flat. Figure A.3 shows the same information but for the kinetics simulation of HD 189733b. The evolution profile shapes are quite different. However, as for HD 209458b, the gradients of the chemical evolution are approximately flat by 1000 days.

Based on the evolution of the global maximum wind velocities and the mole fractions of CH₄ and HCN, all of which show very small gradients by 1000 days, we conclude that our simulations have reached a pseudo-steady state (for P > 10⁵ Pa) by 1000 days. However, we cannot rule out the possibility of long-term chemical evolution that might become apparent for much longer integration times. The challenge is that it is very difficult to run these complex models for very long integration times.

Fig. A.1 Evolution of the global maximum for the zonal (blue), meridional (green) and vertical (red) components of the wind velocities for the kinetics simulation of HD 209458b (top) and HD 189733b (bottom).

Fig. A.2 Evolution of the mole fractions of CH4 (top) and HCN (bottom) on three particular model vertical (altitude) levels (~ 102 Pa, ~ 103 Pa, and ~ 104 Pa) at two arbitrary coordinates (λ = ϕ = 0° and λ = 180°, ϕ = 45°) for the kinetics simulation of HD 209458b.

Fig. A.3 As Fig. A.2 but for the kinetics simulation of HD 189733b.

Appendix B Conservation

The UM does not explicitly enforce global axial angular momentum in the numerical scheme. However, all simulations presented conserve initial global axial angular momentum to within 99%.

An important quantity to conserve is the global (i.e. summed over all grid cells) mass of each element: hydrogen, helium, carbon, oxygen, and nitrogen. The global mass of each chemical species is not a conserved quantity due to chemical transformation. We calculate the global mass of each element from our simulations using Eq. (E10) presented in Drummond et al. (2018b) and show it as a function of model integration time in Fig. B.1 for both planets. The mass of each element shows a very small increase over the course of the simulation but importantly is conserved to within 99.9% at 1000 days, the end of the simulation.

Fig. B.1 Conservation of global element mass shown as fraction of the initial mass for the kinetics simulation of HD 209458b (top) and HD 189733b (bottom). The global element mass is conserved to within better than 99.9% for all elements.

Appendix C Hot initial profile

Figure C.1 shows the pressure-temperature profiles from the chemical equilibrium simulation of HD 209458b initialised with a hotter thermal profile. Specifically the initial profile simply has 800 K added to the temperature at every pressure level, compared with the nominal case. The temperature profiles in Fig. C.1 can be directly compared with the temperature profiles for the nominal HD 209459b simulation in Fig. 3. The only significant difference is in the deep atmosphere (P ≳ 10⁶ Pa) where the nominal simulation shows a much lower temperature as the model has not been evolved for long enough to reach a steady-state. At lower pressures differences in the temperature profiles between the simulations are minimal. There are also negligible differences in the zonal-mean zonal wind (not shown).

Figure C.2 shows abundance profiles of the six key molecules considered in the main body of the paper, for both the equilibrium and kinetics simulations with a hot initial profile. Comparing this figure with abundance profiles for the nominal HD 209459b simulations (Fig. 5) clearly shows that there are negligible differences between the simulations with a “hot” and “cold” initial profile for most of the model domain. Differences are apparent for P > 10⁶ Pa, where all six species are in equilibrium, which is simply a result of the aforementioned temperature difference. Importantly, the deep atmosphere does not effect the lower atmosphere in this particular case, as the quench point for all six molecules lies above (towards lower pressures) the region where the atmosphere has not reached a steady-state.

Fig. C.1 As Fig. 3 but for the simulation of HD 209458b initialised with a hot thermal profile.

Fig. C.2 Vertical mole fraction profiles for different longitude points around the equator, for the equilibrium simulation (dashed lines) and kinetics simulation (solid lines) of HD 209458b initialised with a hot temperature profile.

A further outcome of this experiment is that it shows the insensitivity of the model result to the initial conditions with respect to the chemical abundances. The initial chemical abundances correspond to the chemical equilibrium at the start of the simulation, which is defined by the initial thermal profile. Since the initial profile is uniformly 800 K hotter in this test, compared with the nominal case, this corresponds to a significantly different chemical equilibrium composition. Despite this, the simulations show a very similar result (except in the deep atmosphere)after 1000 days of integration. This suggests that the results are not strongly sensitive to the initial conditions.

Appendix D Abundance profiles at higher latitudes

In this section we present abundance profiles of the six species of interest at a latitude of 45°. These complement the abundance profiles at the equator shown in the main body of the text. The abundance profiles at 45° latitude are shown in Figs. D.1 and D.2 for HD 209458b and HD 189733b, respectively.

Fig. D.1 Vertical mole fraction profiles for different longitude points around the latitude of 45°, for the equilibrium simulation (dashed lines) and kinetics simulation (solid lines) of HD 209458b.

Fig. D.2 Vertical mole fraction profiles for different longitude points around the latitude of 45°, for the equilibrium simulation (dashed lines) and kinetics simulation (solid lines) of HD 189733b.

Appendix E Timescale profiles

In this section, we present profiles of the advection and chemical timescales (for CH₄ and CO) at specific points around the planet. The advection timescales are estimated using the pressure profile of the wind velocity (no spatial averaging) for each component, using the equations presented in Sect. 3.2.1. The chemical timescales are estimated for the Venot et al. (2019) network using the method described in Tsai et al. (2018).

Fig. E.1 Pressure profiles of the advection timescale (showing seperately the zonal (u), meridional (v), and vertical (w) components) and the chemical timescale for CH4 and CO for the HD 209458b simulation for different columns around the planet.

Fig. E.2 As Fig. E.1 but for the HD 189733b simulation.

References

All Tables

All Figures

	Fig. 1 Zonal-mean zonal wind velocity for the equilibrium simulations of HD 209458b (top) and HD 189733b (bottom).
In the text

	Fig. 2 Atmospheric temperature (colour scale and black contours) on the 1 × 104 Pa (0.1 bar) isobar with horizontal wind velocity vectors for the equilibrium simulations of HD 209458b (top) and HD 189733b (bottom).
In the text

Fig. 3 Pressure-temperature profiles extracted from the 3D grid of equilibrium simulations of HD 209458b (top) and HD 189733b (bottom). Solid grey lines show profiles at various longitude points around the equator (0° latitude) while dashed grey lines show profiles at various longitude points at a latitude of 60°. Area-weighted dayside average (red) and nightside average (blue) profiles are shown, along with the 1D radiative-convective equilibrium model profile used to initialise the 3D model (black).

In the text

	Fig. 4 Mole fractions (colour scale and white contours) of CH4, H2O, NH3, CO, CO2 and HCN and wind vectors (arrows) on the 1 × 104 Pa (0.1 bar) isobar for the equilibrium simulation (left column) and the kinetics simulation (right column) of HD 209458b.
In the text

	Fig. 4 continued.
In the text

	Fig. 5 Vertical mole fraction profiles for different longitude points around the equator, for the equilibrium simulation (dashed lines) and kinetics simulation (solid lines) of HD 209458b. Warmer colours indicate a profile that is closer to the substellar point, cooler colours are closer to the antistellar point.
In the text

Fig. 6 Ratios of the chemical to advection timescales for CH4 (α = τchem∕τadv), for the zonal (top), meridional (middle), and vertical (bottom) components of the wind for the kinetics simulation of HD 209458b. We show a meridional-mean between ± 30° latitude for the zonal component, and slices at a longitude of 0° (i.e. the antistellar point) for the meridional and vertical components. The colour scale shows α, with blue colours corresponding to α > 1 and red colours corresponding to α < 1. The black contours show the α = 0.1, 1, 10 values. The fine structure shown in the plots is a reflection of regions where the wind velocity changes signs (e.g. eastward to westward) and the wind velocity approaches a value of zero over a very small region at the transition.

In the text

	Fig. 7 Mole fractions (colour scale and white contours) of CH4, H2O, NH3, CO, CO2, and HCN and wind vectors (arrows) on the 1 × 104 Pa (0.1 bar) isobar for the equilibrium simulation (left column) and the kinetics simulation (right column) of HD 189733b.
In the text

	Fig. 7 continued.
In the text

	Fig. 8 As Fig. 5, but for the equilibrium simulation (dashed lines) and kinetics simulation (solid lines) of HD 189733b.
In the text

	Fig. 9 As Fig. 6, but for the kinetics simulation of HD 189733b.
In the text

	Fig. 10 Transmission (left) and emission (right) spectra for the equilibrium (blue) and kinetics (red) simulations of HD 209458b. Top row: model spectra, middle row: simulated JWST observations and bottom row: simulated ARIEL observations.
In the text

	Fig. 11 Emission phase curves for the equilibrium (solid) and kinetics (dashed) simulations of HD 209458b in a series of spectral bands. The y-axis is the ratio of the planet to stellar flux.
In the text

	Fig. 12 As Fig. 10 but for HD 189733b.
In the text

	Fig. 13 As Fig. 11 but for HD 189733b.
In the text

	Fig. 14 Chemical timescales of CH4 (top) and CO (bottom) estimated for the Venot et al. (2019) chemical network as a function of pressure (in bar, y-axis) and temperature (in K, x-axis). The quantity shown is the common logarithm of the timescale in seconds. The colourscales are different between the two panels.
In the text

	Fig. A.1 Evolution of the global maximum for the zonal (blue), meridional (green) and vertical (red) components of the wind velocities for the kinetics simulation of HD 209458b (top) and HD 189733b (bottom).
In the text

	Fig. A.2 Evolution of the mole fractions of CH4 (top) and HCN (bottom) on three particular model vertical (altitude) levels (~ 102 Pa, ~ 103 Pa, and ~ 104 Pa) at two arbitrary coordinates (λ = ϕ = 0° and λ = 180°, ϕ = 45°) for the kinetics simulation of HD 209458b.
In the text

	Fig. A.3 As Fig. A.2 but for the kinetics simulation of HD 189733b.
In the text

	Fig. B.1 Conservation of global element mass shown as fraction of the initial mass for the kinetics simulation of HD 209458b (top) and HD 189733b (bottom). The global element mass is conserved to within better than 99.9% for all elements.
In the text

	Fig. C.1 As Fig. 3 but for the simulation of HD 209458b initialised with a hot thermal profile.
In the text

	Fig. C.2 Vertical mole fraction profiles for different longitude points around the equator, for the equilibrium simulation (dashed lines) and kinetics simulation (solid lines) of HD 209458b initialised with a hot temperature profile.
In the text

	Fig. D.1 Vertical mole fraction profiles for different longitude points around the latitude of 45°, for the equilibrium simulation (dashed lines) and kinetics simulation (solid lines) of HD 209458b.
In the text

	Fig. D.2 Vertical mole fraction profiles for different longitude points around the latitude of 45°, for the equilibrium simulation (dashed lines) and kinetics simulation (solid lines) of HD 189733b.
In the text

	Fig. E.1 Pressure profiles of the advection timescale (showing seperately the zonal (u), meridional (v), and vertical (w) components) and the chemical timescale for CH4 and CO for the HD 209458b simulation for different columns around the planet.
In the text

	Fig. E.2 As Fig. E.1 but for the HD 189733b simulation.
In the text