# The effects of consistent chemical kinetics calculations on the pressure-temperature profiles and emission spectra of hot Jupiters

B. Drummond, Pascal Tremblin, I. Baraffe, D. Amundsen, J. Mayne, O.

Venot, J. Goyal

► To cite this version:

B. Drummond, Pascal Tremblin, I. Baraffe, D. Amundsen, J. Mayne, et al.. The effects of consistent chemical kinetics calculations on the pressure-temperature profiles and emission spectra of hot Jupiters. Astronomy & Astrophysics - A&A, 2016, 594, 10.1051/0004-6361/201628799. hal-01834127

HAL Id: hal-01834127

https://hal.science/hal-01834127v1

Submitted on 7 Nov 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# The effects of consistent chemical kinetics calculations on the pressure-temperature profiles and emission spectra of hot Jupiters

B. Drummond<fnref n="1" />, P. Tremblin<fnref n="1" />,<fnref n="2" />, I. Baraffe<fnref n="1" />,<fnref n="5" />, D. S. Amundsen<fnref n="1" />,<fnref n="3" />,<fnref n="4" />, N. J. Mayne<fnref n="1" />, O. Venot<fnref n="6" />, and J. Goyal<fnref n="1" />

Received 26 April 2016 / Accepted 13 July 2016

## ABSTRACT

In this work we investigate the impact of calculating non-equilibrium chemical abundances consistently with the temperature structure for the atmospheres of highly-irradiated, close-in gas giant exoplanets. Chemical kinetics models have been widely used in the literature to investigate the chemical compositions of hot Jupiter atmospheres which are expected to be driven away from chemical equilibrium via processes such as vertical mixing and photochemistry. All of these models have so far used pressure-temperature (P-T) profiles as fixed model input. This results in a decoupling of the chemistry from the radiative and thermal properties of the atmosphere, despite the fact that in nature they are intricately linked. We use a one-dimensional radiative-convective equilibrium model, ATMO, which includes a sophisticated chemistry scheme to calculate P-T profiles which are fully consistent with non-equilibrium chemical abundances, including vertical mixing and photochemistry. Our primary conclusion is that, in cases of strong chemical disequilibrium, consistent calculations can lead to differences in the P-T profile of up to 100 K compared to the P-T profile derived assuming chemical equilibrium. This temperature change can, in turn, have important consequences for the chemical abundances themselves as well as for the simulated emission spectra. In particular, we find that performing the chemical kinetics calculation consistently can reduce the overall impact of non-equilibrium chemistry on the observable emission spectrum of hot Jupiters. Simulated observations derived from non-consistent models could thus yield the wrong interpretation. We show that this behaviour is due to the non-consistent models violating the energy budget balance of the atmosphere.

**Key words.** planets and satellites: atmospheres – planets and satellites: composition

## 1. Introduction

Despite the discovery of ever smaller and more Earth-like rocky exoplanets (e.g. Berta-Thompson et al. 2015), hot Jupiters remain one of the most important classes of exoplanet due to the availability of follow-up characterisation observations. This allows for important comparisons between the various observables (transmission spectra (e.g. Sing et al. 2016), emission spectra (e.g. Knutson et al. 2008; Beaulieu et al. 2010; Diamond-Lowe et al. 2014; Evans et al. 2015) and phase curves e.g. Zellem et al. 2014; see Heng & Showman 2015, for a review) with the many atmosphere models in use throughout the community to understand the physical and chemical processes occuring in these atmospheres.

There exists a hierarchy of atmosphere models which have so far been applied to the study of exoplanets with each class of model having its own practical use. The three-dimensional (3D) general circulation models (GCMs) are required to gain insights into the atmospheric dynamics (e.g. Showman et al. 2009; Burrows et al. 2010; Heng et al. 2011; Rauscher & Menou 2012; Dobbs-Dixon & Agol 2013; Mayne et al. 2014; Polichtchouk et al. 2014) whilst the one-dimensional

(1D) radiative-convective models solve for single column pressure-temperature (P-T) profiles (e.g. Iro et al. 2005; Barman et al. 2005; Fortney et al. 2010; Spiegel & Burrows 2010). A further class of atmosphere model, the chemical kinetics models, determine precise chemical compositions (e.g. Moses et al. 2011; Venot et al. 2012; Zahnle & Marley 2014).

Early work using these chemical kinetics models focussed on the photochemical production and destruction processes of hydrocarbon hazes (Liang et al. 2004; Zahnle et al. 2009) and the production of atomic hydrogen through water photolysis (Liang et al. 2003). Later modelling efforts provided detailed chemical compositions of these atmospheres and considered the impact on the observable properties (Moses et al. 2011; Venot et al. 2012; Madhusudhan & Seager 2011).

The non-equilibrium processes of vertical mixing and photochemistry have been shown to have notable consequences on the simulated spectra of some hot exoplanet atmospheres (see Moses 2014, for a review). Vertical mixing can result in transport-induced quenching which increases or decreases the abundances of chemical species compared with their chemical equilibrium profiles. In addition, at low pressures high-energy photons can dissociate molecules into highly reactive daughter products.

<footnote n="1">¹ Astrophysics Group, University of Exeter, EX4 4QL, Exeter, UK
e-mail: bdrummond@astro.ex.ac.uk</footnote>

<footnote n="2">² Maison de la Simulation, CEA-CNRS-INRIA-UPS-UVSQ, USR 3441 Centre d'étude de Saclay, 91191 Gif-Sur-Yvette, France</footnote>

<footnote n="3">³ Department of Applied Physics and Applied Mathematics, Columbia University, New York, NY 10025, USA</footnote>

<footnote n="4">⁴ NASA Goddard Institute for Space Studies, New York, NY 10025, USA</footnote>

<footnote n="5">⁵ Univ. Lyon, ENS de Lyon, Univ. Lyon1, CNRS, CRAL, UMR 5574, 69007 Lyon, France</footnote>

<footnote n="6">⁶ Intituut voor Sterrenkunde, Katholieke Universiteit Leuven, Celestijnenlaan 200D, 3001 Leuven, Belgium</footnote>
These non-equilibrium processes have the potential to influence the simulated spectrum of the atmosphere by changing the amount of absorption or emission from these molecules. For example, Moses et al. (2011) found that increases in the abundances of methane and ammonia due to transport-induced quenching have important consequences for the emission spectrum of HD 189733b, decreasing the calculated eclipse depth over several wavelength regions.

One common feature amongst all of the chemical kinetics studies published so far is that none include a consistent approach to the model atmosphere. The background thermal structure of the atmosphere is treated as a fixed model input. Usually the supplied P-T profiles have been derived from 3D GCM results (e.g. Moses et al. 2011; Venot et al. 2012) or 1D radiative-convective models (e.g. Moses et al. 2013; Agúndez et al. 2014b). The equilibrium and non-equilibrium chemical abundances are then calculated using this fixed P-T profile leading to an inconsistency between the chemical abundances and the thermal structure. As the chemistry is driven away from chemical equilibrium, the opacity in the model atmosphere changes which should impact on the P-T profile.

Agúndez et al. (2014b) have attempted to take into account this effect by including one additional iteration in their radiative-convective model using the initial non-equilibrium chemical abundances they calculate for the hot Neptune GJ 436b. They found corrections to the temperature structure of <100 K. In addition, Hubeny & Burrows (2007) performed consistent non-equilibrium chemistry models of brown dwarf atmospheres using a timescale argument to quench the abundances of N₂/NH₃ and CO/CH₄, not using chemical kinetics, and found corrections to the P-T profile of 50-100 K.

In this study we use a 1D atmosphere code which solves for both hydrostatic and energy balance and includes a sophisticated chemistry scheme to solve for the non-equilibrium chemical abundances, including vertical mixing and photochemistry, consistently with the P-T profile. In Sect. 2 we describe our model setup. In Sect. 3 we present our results and, finally, in Sect. 4 we summarise and conclude.

## 2. Method and model description

We use the 1D atmosphere code ATMO (Amundsen et al. 2014; Tremblin et al. 2015, 2016) to generate models of two highly irradiated exoplanet atmospheres. The model has already been applied to the atmospheres of brown dwarfs (Tremblin et al. 2015, 2016) and we use essentially the same model setup here except that external irradiation is now included. In this section we present details of the model and our method for calculating non-equilibrium chemical abundances consistently with the P-T profile.

### 2.1. Radiative-convective equilibrium

ATMO solves for the pressure-temperature structure of an atmosphere by finding energy flux balance in each model level, i.e.,

$$ \int_{0}^{\infty} (F_{\text{rad}}(v) + F_{\text{star}}^{0}(v)e^{\tau_{v}/\mu_{\text{star}}}) dv + F_{\text{conv}} = \sigma T_{\text{int}}^{4} \quad (1) $$

where $F_{\text{rad}}(v)$ and $F_{\text{conv}}$ are the spectral radiative flux and the convective flux, respectively, and $\tau_{v}$ is the vertical monochromatic optical depth. $T_{\text{int}}$ is the internal temperature of the object corresponding to the surface flux at which the object cools in the

absence of irradiation; $T_{\text{int}}$ is equivalent to the effective temperature $T_{\text{eff}}$ in the absence of irradiation. $F_{\text{star}}^{0}(v)$ is the incoming irradiation flux from the star at the top of the atmosphere and $\mu_{\text{star}} = \cos \theta$ where $\theta$ is the angle of incoming radiation off the vertical; $\mu_{\text{star}}$ is negative. In addition, the pressure as a function of altitude is calculated using the equation of hydrostatic balance

$$ \frac{\partial P}{\partial z} = -\rho g \quad (2) $$

where $P$ is the pressure, $z$ is the altitude, $\rho$ is the density and $g$ is the gravity. An object's atmosphere is therefore completely determined by its internal temperature and gravity, often specified as $(T_{\text{int}}, \log(g))$, in the absence of irradiation. The convective flux is computed using the mixing length theory as presented by Henyey et al. (1965)

$$ F_{\text{conv}} = \frac{1}{2} \rho C_{\text{p}} T v_{\text{conv}} \frac{l}{H_{\text{p}}} \frac{\Gamma}{\Gamma + 1} (\nabla_{T} - \nabla_{\text{ad}}) \quad (3) $$

with $v_{\text{conv}}$ the convective velocity, $\Gamma$ the efficiency parameter, $\nabla_{T} = \partial \log(T)/\partial \log(P)$, $\nabla_{\text{ad}}$ the adiabatic gradient, $C_{\text{p}}$ the specific heat at constant pressure, $H_{\text{p}}$ the pressure scale height, and $l = \alpha H_{\text{p}}$ the mixing length. The details of the computation of $v_{\text{conv}}$ and $\Gamma$ can be found in Gustafsson et al. (2008) and we used the same standard parameters (temperature distribution in a convective element $y = 0.076$, energy dissipation by turbulent viscosity $\nu = 8$, mixing length<fnref n="1" /> $\alpha = 1.5$). The non-linear problem given by Eqs. (1) and (2) is solved by a Newton-Raphson algorithm to obtain the converged P-T profile of the atmosphere.

### 2.2. Radiative transfer

The radiative transfer equation is solved in 1D plane-parallel geometry including isotropic scattering. We include CH₄, H₂O, CO, CO₂, NH₃, TiO, VO, Na, K, Li, Cs and Rb and collision induced absorption due to H₂-H₂ and H₂-He as opacity sources in the atmosphere (see Tremblin et al. 2015; Amundsen et al. 2014, for details). ATMO can be used as a full line-by-line code at high spectral resolution (0.001 cm⁻<fnref n="1" />, evenly spaced in wavenumber), e.g. to post process high-resolution emission and transmission spectra, though this is not used in this study.

Presently ATMO is used with the correlated-k approximation using the random overlap method to compute the total mixture opacity (see Lacis & Oinas 1991; Amundsen et al., in prep.) for moderate resolution emission and transmission spectra with 500 or 5000 bands and, for a rapid computation of the radiative flux, with 32 bands in the Newton-Raphson iterations. The number of k-coefficients per band is not fixed but is instead computed based on a specified precision; the number of k-coefficients therefore varies depending on the band and the gas (this is explained in detail in Amundsen et al. 2014). The k-coefficients were derived for opacity tables as in Amundsen et al. (2014). The results obtained using the correlated-k approximation have been compared with the full line-by-line result, which agree very well.

The radiative transfer equation for the irradiation and thermal emission of the planet is split into a directly-irradiated part $I_{\text{star}}(v)$ and a thermal plus diffusely-irradiated part

<footnote n="1">¹ The choice of $\alpha$ is based on typical values, namely ~1.5–2, used in cool object atmosphere models (low mass stars and brown dwarfs) and based on calibration of 1D models with multi-D radiative hydrodynamical simulations (e.g. Baraffe et al. 2015). Such calibration has never been performed for hot Jupiter-like planets but we note that present hot Jupiter models are rather insensitive to this parameter for any value $\ge 1.5$.</footnote>
$I_{\text{rad}}(\nu)I_{\text{therm}}(\nu)+I_{\text{diff}}(\nu)$ (see Thomas & Stammes 1999, Chap. 6.2 for details). The direct part can be integrated straightforwardly to give

$$F_{\text{star}}(\nu) = 2\pi \int_{-1}^{1} I_{\text{star}}(\nu)\mu d\mu = F_{\text{star}}^{0}(\nu) e^{\tau_{\nu}/\mu_{\text{star}}} \quad (4)$$

where $F_{\text{star}}^{0}(\nu)$ is the top incoming flux defined as $F_{\text{star}}^{0}(\nu) = 4\pi H_{\text{star}}^{0}(\nu)\mu_{\text{star}}R_{\text{star}}^{2}/R_{\text{orbit}}^{2}$; here $H_{\text{star}}^{0}(\nu)$ is the flux at the surface of the star, and $R_{\text{star}}$ and $R_{\text{orbit}}$ are the stellar and orbital radius respectively. The diffuse and thermal part are computed using the radiative transfer equation with isotropic scattering (Thomas & Stammes 1999; Hubeny & Mihalas 2014)

$$\mu \frac{dI_{\text{rad}}(\nu)}{d\tau_{\nu}} = I_{\text{rad}}(\nu) - (1 - \epsilon_{\nu})J_{\text{rad}}(\nu) - (1 - \epsilon_{\nu})F_{\text{star}}(\nu)/4\pi - \epsilon_{\nu}B(\nu) \quad (5)$$

where $\epsilon_{\nu}$ is the photon destruction probability given as a function of absorption and scattering opacity $\kappa_{\text{abs}}/(\kappa_{\text{scat}} + \kappa_{\text{abs}})$, $B(\nu)$ is the Planck function, and

$$\begin{align}
J_{\text{rad}}(\nu) &= \frac{1}{2} \int_{-1}^{1} I_{\text{rad}}(\nu) d\mu \nonumber \\
F_{\text{rad}}(\nu) &= 2\pi \int_{-1}^{1} I_{\text{rad}}(\nu) \mu d\mu. \tag{6}
\end{align}$$

The ray directions specified by $\mu$ are sampled with a discrete-ordinate method using Gauss-Legendre quadrature (we usually use 16 rays in total). We solve the radiative transfer equation iteratively in its integral form following Bendicho & Bueno (1995) and using a Gauss-Seidel-type $\Lambda$-acceleration scheme for the scattering

$$\begin{align}
I_{\text{rad}}(\nu) &= I_{\text{rad}}^{0}(\nu) e^{-(\tau_{\nu}/\mu)} + \int_{0}^{\tau_{\nu}/\mu} S_{\nu}(t) e^{\tau' - (\tau_{\nu}/\mu)} d\tau' \nonumber \\
S_{\nu} &= (1 - \epsilon_{\nu}) J_{\text{rad}}(\nu) + (1 - \epsilon_{\nu}) F_{\text{star}}(\nu) / 4\pi + \epsilon_{\nu} B(\nu). \tag{7}
\end{align}$$

The radiative-transfer scheme has been benchmarked against the Met Office SOCRATES code (Amundsen et al. 2014). All results presented in this study include isotropic scattering.

## 2.3. Chemistry: gibbs energy minimisation

The Gibbs energy minimisation scheme is one of two chemistry schemes included in the model. Minimisation of the Gibbs energy is a popular method for obtaining chemical equilibrium abundances in atmosphere models (Burrows & Sharp 1999; Blecic et al. 2016). This technique does not rely on complicated chemical networks and is therefore much simpler to implement and compute than a chemical kinetics method.

We follow the method of Gordon & McBride (1994), using the same thermochemical data as Venot et al. (2012) in the form of NASA polynomial coefficients (McBride et al. 1993). The Gibbs mimimisation method allows for depletion of gas phase species due to condensation, however we do not consider condensation here and this is left for a future study.

At constant temperature and pressure, the Gibbs energy of the system, G, which is given by

$$G = \sum_{i=1}^{I} \mu_{i} n_{i}, \quad (8)$$

where $\mu_{i}$ and $n_{i}$ are the chemical potential and number of moles of chemical species i, must be minimised subject to the constraint of elemental conservation,

$$\sum_{i=1}^{I} a_{ij} n_{i} - b_{j}^{0} = 0, \quad (9)$$

where $a_{ij}$ is the number of atoms of element j per molecule i and $b_{j}^{0}$ is the total number of moles of element j in the mixture. This minimisation is performed using the method of Lagrange multipliers.

We include ~140 chemical species in total, with 105 of those species being the ones contained in the Venot et al. (2012) chemical network. In addition to those, we include Na, K, Li, Cs and Rb in monatomic form and their important molecular gas-phase species (e.g. NaCl, NaOH,...), which can deplete the abundance of the monatomic forms. Local chemical equilibrium abundances are acheived by minimising the Gibbs energy independently in each model level.

We benchmark the Gibbs minimisation scheme against the chemical equilibrium analytical formulae of Burrows & Sharp (1999) and Heng & Tsai (2016) in Appendix A.2.1.

## 2.4. Chemistry: chemical kinetics

We also include a chemical kinetics scheme which deals directly with chemical reactions. We solve for the chemical steady-state by integrating the continuity equation,

$$\frac{\partial n_i}{\partial t} = P_i - n_i L_i - \frac{\partial \Phi_i}{\partial z}, \quad (10)$$

so that $\frac{\partial n_i}{\partial t} \sim 0$. Here, $P_i$ and $L_i$ are the chemical production and chemical loss terms, derived from the system of chemical reactions in the chemical network, and the $\frac{\partial \Phi_i}{\partial z}$ term describes vertical mixing.

The vertical transport flux $\Phi_i$ is split into two components due to molecular diffusion and eddy diffusion and is described by (e.g. Gladstone et al. 1996),

$$\begin{align}
\Phi_i = -D_i \left( \frac{\partial n_i}{\partial z} + \frac{n_i}{H_i} + \frac{n_i(1+\alpha)}{T} \frac{dT}{dz} \right) \nonumber \\
- K_{zz} \left( \frac{\partial n_i}{\partial z} + \frac{n_i}{H_a} + \frac{n_i}{T} \frac{dT}{dz} \right). \tag{11}
\end{align}$$

where $D_i$ and $K_{zz}$ are the molecular and eddy diffusion coefficients, respectively, and $H_i$ and $H_a$ are the scale heights of each individual species and of the bulk atmosphere, respectively, and $\alpha$ is the thermal diffusion parameter. Equation (11) can be written more simply in terms of the mole fraction $f_i = n_i/n$

$$\Phi_i = -nD_i \left( \frac{\partial f_i}{\partial z} - \frac{f_i}{H_a} + \frac{f_i}{H_i} + \frac{f_i \alpha}{T} \frac{dT}{dz} \right) - nK_{zz} \left( \frac{\partial f_i}{\partial z} \right). \quad (12)$$

The molecular diffusion coefficient $D_i$ can be determined using the kinetic theory of gases (e.g. Wayne 1991) and is inversely proportional to the number density. Therefore, molecular diffusion becomes important at low pressures ($P < 10^{-5}$ bar) in the thermospheric regions of the atmosphere. At higher pressures, vertical transport is dominated by the eddy diffusion term for which the controlling coefficient $K_{zz}$ is far less well constrained, both observationally and theoretically. In previous studies, this term has been estimated using wind velocity fields derived from
3D general circulation models (e.g. Moses et al. 2011) or pa-
rameterised using the advection of passive tracers (Parmentier
et al. 2013); however if the turbulent diffusion is caused by small
scale processes linked to convection, overshooting of the con-
 vective region, gravity waves or fingering convection (Tremblin
et al. 2016), then care must be taken when deriving this parame-
ter from 3D simulations in which the distinction should be made
between the large scale vertical advection linked to the jet circu-
lation and the small scale turbulent diffusion processes that are
likely to be unresolved. Other studies treat the K<sub>zz</sub> parameter as
a free model parameter and test a range of plausible values (e.g.
Miguel & Kaltenegger 2014).

In this work, we choose two plausible values for the K<sub>zz</sub> pa-
rameter and show models for both cases. Values of the K<sub>zz</sub> pa-
rameter for hot Jupiter atmospheres used in other works vary
between ~10<sup>7</sup>-10<sup>12</sup> cm<sup>2</sup> s<sup>-1</sup> (e.g. Moses et al. 2011; Miguel &
Kaltenegger 2014). We take one case roughly in the middle of
this range (K<sub>zz</sub> = 10<sup>9</sup> cm<sup>2</sup> s<sup>-1</sup>) and one towards the upper limit
of this range (K<sub>zz</sub> = 10<sup>11</sup> cm<sup>2</sup> s<sup>-1</sup>). We impose zero-flux limits
on both the upper- and lower-boundaries, assuming that no mass
enters the atmosphere at the bottom or escapes the atmosphere
at the top. In addition, we chose the value of K<sub>zz</sub> and maximum
pressure level of our model such that the lower boundary will al-
ways remain in chemical equilibrium, meaning that our assump-
tion of zero-flux at the lower boundary is valid (Moses et al.
2011).

The chemical production and loss terms, *P*<sub>*i*</sub> and *L*<sub>*i*</sub>, are cal-
culated from the system of chemical reactions contained within
the chemical network employed. In this study, we adopt the
Venot et al. (2012) chemical network previously applied in sev-
eral studies of hot, hydrogen-dominated exoplanet atmospheres
(Agúndez et al. 2012, 2014a,b; Venot et al. 2014) and for brown
dwarf atmospheres (Tremblin et al. 2015, 2016).

We choose to use the original C₀–C₂ network of Venot et al. (2012) with hydrocarbon species of up to two carbon atoms, rather than the more recent update which includes higher-order hydrocarbons C₀–C₆ (Venot et al. 2015). The overall aim of this work is to study the effect of non-equilibrium chemistry on the background atmosphere and it is therefore most important to obtain accurate abundances for the species contributing to the opacity. The inclusion of higher-order hydrocarbons was found not to affect the abundances of the main species (e.g. CH₄, H₂O) (Venot et al. 2015). Therefore we choose to use the smaller and less computationally expensive C₀–C₂ network.

We use the LSODE solver<fnref n="2" /> (Hindmarsh 1983) using the DLSODES option for stiff system of ordinary differential equations to solve the system of continuity equations simultaneously (Eq. (10)) for each chemical species using the Backwards Differentiation Formulae method, until at least 10¹<fnref n="2" /> s when the chemistry has reached a steady-state. There are two options to choose the timestep of the chemistry iterations: a customised version which starts with an initial timestep of dt = 10⁻¹⁰ s which is progressively increased, by 10%, if the maximum relative variations of the molecular abundances do not exceed 10%, or we employ the timestep which is internally generated within the DLSODES routine. We have tested both methods, which give the same results for a given time. Since the timestep generated by DLSODES is typically larger than that using our custom method, we generally use the second option.

As well as vertical mixing, dissociations of molecules by en-
ergetic X-ray and UV (XUV) photons can drive the chemistry
away from local chemical equilibrium at low pressures where

the flux of these photons is high. The XUV photon flux is calcu-
lated as a function of pressure, using the same radiative transfer
scheme as described in Sect. 2.2, accounting for scattering due to
H₂ and He. We use the same molecule cross sections and quan-
tum yields as described in Venot et al. (2012).

We benchmark our chemical kinetics scheme with that of
Venot et al. (2012) in Sect. A.2.2.

2.5. Fully-consistent modelling

Since the chemical kinetics scheme is consistently coupled with
the radiative-convective scheme it allows the non-equilibrium
chemical abundances to be solved consistently with the P–T
profile. Previous chemical kinetics models of hot Jupiter atmo-
spheres supply the P–T profile as a fixed model input, on which
the chemical abundances are calculated (e.g. Line et al. 2010;
Moses et al. 2011; Venot et al. 2012; Kopparapu et al. 2012).
Our approach allows changes in opacity, due to non-equilibrium
chemistry, to feedback into the radiative-convective balance and
hence maintain a consistent temperature structure.

Initially, we solve for the *P*–*T* profile assuming local chem-
ical equilibrium. That is, we solve for radiative-convective bal-
ance using the Gibbs energy minimisation scheme to re-calculate
the abundances on every iteration of the Newton-Raphson solver.
Once radiative-convective balance (1) is found, a 1D *P*–*T* profile
and corresponding abundance profiles of all the chemical species
is obtained.

To find the *P*–*T* profile consistent with non-equilibrium chemistry the model is restarted, switching to the chemical kinetics scheme including vertical mixing and photochemistry. The chemical equilibrium abundances and corresponding *P*–*T* profile are used to initialise the model. Note that in the absence of vertical mixing and photochemistry the derived abundances from the Gibbs minimisation scheme and the chemical kinetics scheme agree very well; see Sect. A.2.2 and Venot et al. (2012). As the continuity equation (Eq. (10)) is solved, as outlined above, we periodically solve for the new radiative-convective balance using the up-to-date chemical abundances which may have changed due to vertical mixing and/or photochemical processes. We re-converge to find radiative-convective balance every 10 timesteps of the DLSODES solver; this value was chosen as the chemical abundances do not vary by a large amount within 10 iterations.

For convergence of the consistent model we use the criteria
that the error in the energy flux balance Eq. (1) is <10⁻⁴. For the
kinetics calculations we ensure that the maximum relative error
of all the chemical species (max(dnᵢ/nᵢ)) is <10⁻⁴ and the model
has reached an integration time of 10¹<fnref n="2" /> s. Overall model conver-
gence is deemed to be reached when both of these criteria are
satisfied simultaneously. We tested using a higher tolerance and
found that this is a good balance between accuracy and perfor-
mance. The resulting *P*–*T* profile and non-equilibrium chemical
abundance profiles are fully consistent with each other. In addi-
tion, we also perform a chemical kinetics calculation with the
*P*–*T* profile held fixed.

2.6. Test case planets

In this first application of ATMO to highly irradiated exoplanet
atmospheres we choose to model two hot Jupiters: HD 189733b
and HD 209458b. These atmospheres are well studied both ob-
servationally and theoretically. In this work we do not attempt to
match the various observations of these atmospheres. Instead we

<footnote n="2">² https://computation.lnl.gov/casc/odepack/</footnote>
<table><caption>**Table 1.** Parameters used in this study.</caption><thead><tr><th></th><th>HD 209458b</th><th>HD 189733b</th></tr></thead><tbody><tr><td>Mass <i>M</i><sub>J</sub></td><td>0.714</td><td>1.150</td></tr><tr><td>Radius <i>R</i><sub>J</sub></td><td>1.380</td><td>1.151</td></tr><tr><td>Orbital dist. AU</td><td>0.047</td><td>0.031</td></tr><tr><td>Stellar spectral type</td><td>G0V</td><td>K1-K2</td></tr><tr><td>Stellar radius <i>R</i><sub>Sun</sub></td><td>1.162</td><td>0.752</td></tr></tbody></table>

focus on the theoretical implications of non-equilibrium chemistry on both the temperature structure and the emission spectrum.

We model the pressure range $10^3$ to $10^{-5}$ bar with 100 vertical levels; we have tested with both higher and lower vertical resolutions and found that this a good balance between model accuracy (most importantly accurately capturing the location of the quenching point) and computational efficiency. We do not model lower pressures since the atmosphere transitions into the thermosphere around this region and a different modelling approach is required (e.g. Yelle 2004; Koskinen et al. 2013), requiring consideration of non-local thermodynamic equilibrium. Therefore, due to the lower pressure limit of our model, though we include photochemistry, we only capture the “edge” of the photodissociation zone; photodissociations become the dominant driver of the chemistry for $P < 10^{-5}$ bar. We assume an internal temperature of 100 K for both planets. For both cases we have assumed a solar composition with elemental abundances of Caffau et al. (2011).

The planetary, orbital and stellar parameters of Southworth (2010) are adopted in both cases and are summarised in Table 1. We use the Kurucz stellar spectra<fnref n="3" /> of both HD 189733 and HD 209458 to compute the short-wave irradiation. For HD 209458 we use an XUV spectrum based on the solar spectrum, as in Venot et al. (2012), where the solar XUV spectrum is scaled to the temperature and radius of HD 209458 with a scaling factor of $(R_{\text{star}}/R_{\text{o}})^2(T_{\text{star}}/T_{\text{o}})^4 \approx 1.8$. For the HD 189733 XUV flux we follow the same approach as Moses et al. (2011): the UV flux of a K2 V star, epsilon Eridani, taken from the CoolCAT database<fnref n="4" />, was used in range 115–230 nm, the solar UV flux for wavelengths below this range, and the solar flux divided by ten for longer wavelengths.

### 3. Results

In this section we present the *P*–*T* profiles, chemical abundances and simulated spectra derived from our 1D atmosphere code ATMO for the two test case hot Jupiter atmospheres outlined in the previous section.

For both planets we present a series of models with different approaches to the calculation of the chemical abundances. We show models which are consistent with local chemical equilibrium (EQ). We also include non-consistent non-equilibrium (NEQ) models, with chemical abundances derived from chemical kinetics calculations, including the processes of vertical mixing and photochemistry, where we hold the *P*–*T* profile fixed. Finally, we include consistent non-equilibrium (CNEQ) models where the abundances are, again, derived from chemical kinetics but the *P*–*T* profile is not fixed and is instead allowed to vary to maintain radiative-convective balance.

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Dayside average *P*–*T* profiles for HD 189733b assuming equilibrium chemistry (EQ model, dashed), consistent non-equilibrium chemistry with $K_{\text{zz}} = 10^{11}$ cm$^2$ s$^{-1}$ (CNEQ model, solid) and with $K_{\text{zz}} = 10^9$ cm$^2$ s$^{-1}$ (CNEQ model, dotted). Note that the NEQ models referred to in the text use the same EQ *P*–*T* profile as plotted in this figure.</figcaption></figure>

#### 3.1. HD 189733b

Figure 1 shows the *P*–*T* profiles of our EQ and CNEQ HD 189733b models using two different values for the $K_{\text{zz}}$ parameter, as described in Sect. 2. For the model with stronger vertical mixing there is a significant influence of non-equilibrium chemistry on the *P*–*T* profile. The CNEQ *P*–*T* profile is ~100 K warmer than the EQ profile for $P > 0.1$ bar. In the model with the smaller $K_{\text{zz}}$ parameter the CNEQ *P*–*T* profile is warmer than the EQ profile by about 15 K.

The equilibrium and consistent non-equilibrium chemical abundances for these HD 189733b models are shown in Figs. 2 and 3 for the models using $K_{\text{zz}} = 10^{11}$ and $10^9$ cm$^2$ s$^{-1}$, respectively. Qualitatively, we find similar behaviour to previous studies (Moses et al. 2011; Venot et al. 2012) where the chemistry remains in chemical equilibrium in the hot deep atmosphere, the mid-regions of the atmosphere are dominated by vertical mixing, and photochemistry begins to become important for $P \sim 10^{-5}$ bar.

The model using the stronger $K_{\text{zz}}$ shows significant increases in the abundance of both CH$_4$ and NH$_3$ compared to chemical equilibrium, which are quenched at around 10 bar and 100 bar, respectively. In chemical equilibrium, H$_2$O is more abundant than CO in the deep atmosphere, with CO becoming more abundant than H$_2$O at around 10 bar. This transition between the two molecules is removed with the inclusion of vertical mixing, as H$_2$O and CO are quenched below the transition, increasing the abundance of H$_2$O. The model with the smaller eddy diffusion coefficient shows a smaller increase in both CH$_4$ and NH$_3$ as their quench points are both shifted to lower pressures, reducing the quenched abundance. The effect on H$_2$O and CO is also much smaller in this model.

The temperature change due to non-equilibrium chemistry has a feedback impact on the temperature-dependent chemical abundances. Comparing the chemical abundances between the consistent and non-consistent (CNEQ and NEQ) models we see that for the strong vertical mixing case (Fig. 4) the abundances of CO and CH$_4$ are essentially reversed, due to an increase in CO and a decrease in CH$_4$ in the CNEQ model. Similarly, in the NEQ model we find that NH$_3$ should be the dominant nitrogen species throughout the atmosphere, whereas in the CNEQ model we find that N$_2$ is the dominant nitrogen species (below 100 bar).

<footnote n="3">³ http://kurucz.harvard.edu/stars.html</footnote>

<footnote n="4">⁴ http://casa.colorado.edu/~ayres/CoolCAT/</footnote>
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Chemical abundances of the major chemical species for the HD 189733b model with abundances from the EQ calculation (dashed) and abundances from the CNEQ calculation including vertical mixing ($K_{zz} = 10^{11}$ cm$^2$ s$^{-1}$) and photochemistry (solid).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Fig. 3.** Same as Fig. 2 but using $K_{zz} = 10^9$ cm$^2$ s$^{-1}$.</figcaption></figure>

The temperature change induces different abundances via two processes. Firstly, the chemistry in the deep atmosphere remains in chemical equilibrium, and as the temperature increases the chemistry reaches a new chemical equilibrium with abundances which are consistent with the new temperature.

The second effect is caused by a shifting of the quench point. The quench point occurs at the pressure level where the chemical timescale $\tau_{chem}$ is equal to the mixing timescale $\tau_{mix}$. Since $\tau_{chem}$ is dependent on temperature the quench point is shifted to a lower pressure level, in the case of a warmer atmosphere. This leads to a different quenched mole fraction effecting the quenched abundances for pressures below the quench point. This is complicated by the fact that the mole fractions of individual species have changed in the region of the quench point as they now exist in the new higher temperature chemical equilibrium, due to the first process explained.

This process also occurs in the model with weaker vertical mixing (Fig. 5) but to a smaller degree since the departure from chemical equilibrium and the induced temperature change is less.

Tests with only photochemistry included (i.e. without vertical mixing) show that the photochemistry has a negligible impact on the P–T profile. Photodissociations become important for pressures below $\sim 10^{-5}$ bar depending on the temperature and

<figure><img src="image_5.png" /><figcaption>**Fig. 4.** A comparison of the chemical abundances between the CNEQ model (solid) and the NEQ model (dashed) for the HD 189733b $K_{zz} = 10^{11}$ cm$^2$ s$^{-1}$ case.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Fig. 5.** Same as Fig. 4 but for the $K_{zz} = 10^9$ cm$^2$ s$^{-1}$ case.</figcaption></figure>

the UV flux. At these pressure levels the optical depth is small and changes in the chemical composition have negligible effect on the temperature structure. It is transport-induced quenching which is effective at higher pressures and higher optical depths which has the potential to alter temperature structure. However, it may be that photochemical production of chemical species not included in our model, or not included as opacity sources in our model, could contribute to heating at low pressures, for example ozone in the Earth atmosphere.

Overall, we find that for our model atmosphere of HD 189733b the process of transport-induced quenching causes temperature increases of up to 100 K between $10^{-1}$–$10^2$ bar. This temperature increase, in turn, effects the calculated mole fractions by 1) inducing a new chemical equilibrium consistent with the higher temperature in the deep atmosphere and 2) shifting the quench point ($\tau_{chem} = \tau_{mix}$) to lower pressures and altering the quenched abundances at low pressures; these two processes act simultaneously. For the strong vertical mixing case presented here, this results in an atmosphere where CO is the dominant carbon species for $P < 10$ bar. This contrasts with the non-consistent calculation where $CH_4$ is expected to be the dominant carbon species throughout the whole atmosphere. A similar process occurs for the $N_2$–$NH_3$ system leading to an $N_2$ dominated atmosphere in the CNEQ model but an $NH_3$ dominated atmosphere in the NEQ model.
<figure><img src="image_7.png" /><figcaption>**Fig. 6.** Emission spectra of the HD 189733b model with $K_{zz} = 10^{11}$ cm$^2$ s$^{-1}$, showing calculations based on the EQ calculation (blue), the NEQ calculation (red) and the CNEQ calculation (green).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Fig. 7.** Same as Fig. 6 but for the $K_{zz} = 10^9$ cm$^2$ s$^{-1}$ model.</figcaption></figure>

**Simulated emission spectra.** In this section we present the simulated emission spectra for this series of models of HD 189733b. Figure 6 shows the simulated emission spectrum for the stronger vertical mixing case for all three chemistry models (EQ, NEQ and CNEQ). Non-equilibrium chemistry has a strong impact on the simulated emission spectrum in the NEQ model. The NEQ spectrum has a significantly reduced flux ratio compared with the EQ model at almost all wavelengths. On the other hand, the CNEQ model shows a smaller discrepancy with the EQ model, except at around 4.5 µm where the CNEQ model shows a greater flux ratio than both the EQ and NEQ models.

The model with the lower $K_{zz}$ value (Fig. 7) shows similar trends, though the difference between the three chemistry cases is smaller as the departure from chemical equilibrium is not as strong. Our EQ and NEQ simulated spectra for this case agree well with the spectra of the “thermochemical model” and “pho-tochemical model” of Moses et al. (2011, their Fig. 11). We also find a reduction in flux at around 4 µm and at longer wavelengths for our NEQ model. Interestingly, however, performing the chemical kinetics calculations consistently (CNEQ model) completely removes this signature of non-equilibrium chemistry at 4 µm and also reduces the impact at longer wavelengths.

Most of the spectral features here are due to CH$_4$ (particularly around 3.6 µm) and, at longer wavelengths, to NH$_3$ whilst CO is the dominant absorber around 4.5 µm. Increases in the mole fractions of CH$_4$ and NH$_3$ due to transport-induced quenching increase the opacity in the wavelengths regions where they have absorption bands.

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** Dayside average $P-T$ profiles for HD 209458b assuming equilibrium chemistry (EQ models, dashed), consistent non-equilibrium chemistry with $K_{zz} = 10^{11}$ cm$^2$ s$^{-1}$ (CNEQ models, solid) and with $K_{zz} = 10^9$ cm$^2$ s$^{-1}$ (CNEQ models, dotted), with the model including TiO and VO (red) and without TiO/VO (black). Note that the NEQ models referred to in the text use the same EQ $P-T$ profile as plotted in this figure.</figcaption></figure>

## 3.2. HD 209458b

We now present a series of models for the atmosphere of HD 209458b, which is warmer than the atmosphere of HD 189733b due to its orbit around a warmer G-type star, see Table 1. There has been much debate about the presence of a thermal inversion in the atmosphere of HD 209458b, with early observations favouring the presence of an inversion (Knutson et al. 2008; Beaulieu et al. 2010). However, more recent re-analyses of these datasets suggest that this atmosphere does not contain a temperature inversion (Diamond-Lowe et al. 2014; Evans et al. 2015). In this study, we present both cases, as a way to explore a larger diversity of atmosphere types.

Figure 8 shows the $P-T$ profiles for the two different atmosphere types of HD 209458b. In each case we show $P-T$ profiles for both the EQ and CNEQ chemistry models and with two different strengths of vertical mixing: $K_{zz} = 10^{11}$ and $K_{zz} = 10^9$ cm$^2$ s$^{-1}$, as in the previous section. For the rest of this analysis, however, we discuss only the $K_{zz} = 10^{11}$ cm$^2$ s$^{-1}$ case, where the impact of non-equilibrium chemistry is larger, and instead focus on the difference between the temperature inversion and non-temperature inversion cases.

For the case without a temperature inversion there are very small differences between the EQ and CNEQ $P-T$ profiles. For pressures greater than ~0.1 bar the CNEQ model is marginally warmer than the EQ model; in the isothermal plateau region (1 bar < $P$ < 100 bar) the temperature is about <10 K warmer. On the other hand, the case with a temperature inversion shows a much greater discrepancy between the EQ and CNEQ models. The CNEQ $P-T$ profile is >100 K hotter than the EQ $P-T$ profile for pressures greater than 0.1 bar. At lower pressures the CNEQ $P-T$ profile is cooler than the EQ case by a similar amount. Interestingly, the position of the thermal inversion is also shifted to lower pressures.

Figure 9 shows the abundances of the major chemical species for the model without a temperature inversion. NH$_3$ and CH$_4$ are quenched at around 10 and 1 bar respectively. This has the effect of increasing their abundances with respect to chemical equilibrium by several orders of magnitude for pressures lower than the quench point. Despite this, their molar fractions do not exceed $1\times 10^{-6}$ and N$_2$ and CO remain the dominant nitrogen and carbon
<figure><img src="image_10.png" /><figcaption>**Fig. 9.** Chemical abundances of the major chemical species for the HD 209458b model without a temperature inversion with abundances from the EQ calculation (dashed) and abundances from the CNEQ calculation including vertical mixing and photochemistry (solid); $K_{zz} = 10^{11}$ cm$^2$ s$^{-1}$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Fig. 10.** Same as Fig. 9 but for the HD 209458b model with a temperature inversion.</figcaption></figure>

species, which are unaffected by vertical mixing processes, and retain constant mixing ratios below 10 bar. The effects of photochemistry can be seen in the very upper regions of the model, particularly by the dissociation of $NH_3$ at 0.1 mbar and photochemical production of atomic H from 10 mbar.

Figure 10 shows the abundances for the model with a temperature inversion which shows a much greater departure from chemical equilibrium. Though this model contains a hotter upper atmosphere, at depth the atmosphere is actually considerably cooler than the model without a temperature inversion. TiO and VO absorb visible photons at low pressures forming the temperature inversion. However, this reduces the flux of high energy photons which penetrate to depth and heat the lower atmosphere. This reduction of heating at high pressures leads to a cooler deep atmosphere. Due to the temperature dependence of the chemical timescale the quench point is shifted to higher pressures in the temperature inversion model, leading to larger quenched abundances of both $CH_4$ and $NH_3$. The quench points for $CH_4$ and $NH_3$ now lie at 10 and 100 bar, respectively. At around 1 mbar, where the temperature begins to increase again, the chemical timescale begins to speed up once more, and the species begin to move back towards their chemical equilibrium state, as seen in previous studies (Moses et al. 2011; Venot et al. 2012).

<figure><img src="image_12.png" /><figcaption>**Fig. 11.** A comparison of the chemical abundances between the CNEQ calculation (solid) and the NEQ calculation (dashed) for the HD 209458b model without a temperature inversion.</figcaption></figure>

The chemical abundances in the CNEQ and NEQ models (Fig. 12) show important differences. The NEQ model gives a $CH_4$ mole fraction around 5× larger than the CNEQ model. Similarly, the $NH_3$ mole fraction is ~3.5× larger in the NEQ model compared with the CNEQ model. These discrepancies between the two models occur in the pressure range where observations are available. The differences are, again, due to the increase in temperature at depth which 1) changes the chemical equilibrium abundances in the deep atmosphere and 2) shifts the location of the quench point.

There is negligible difference between the CNEQ and NEQ abundances for the model without a temperature inversion (Fig. 11), only minor decreases in $CH_4$ and $NH_3$ in the CNEQ case. This is to be expected since the temperature difference between the CNEQ and NEQ model is also very small.

We find that in the case of HD 209458b hosting a temperature inversion the deep atmosphere is cool enough to have significant non-equilibrium chemical abundances at depth. This leads to temperature increases of more than 100 K between $10^{-1}$–$10^2$ bar, with decreases in the temperature at lower pressures. The location of the temperature inversion is also shifted to lower pressures. These temperature changes have important consequences on the chemical abundances and we find significantly smaller mole fractions of $CH_4$ and $NH_3$ in the consistent (CNEQ) model compared with the non-consistent (NEQ) model.

For the model without a temperature inversion the deep atmosphere is much hotter and the quench point exists at much lower pressures. This limits the influence of non-equilibrium chemistry and only a small temperature change is seen even for the strong vertical mixing case. Consequently there is little difference between the mole fractions of the consistent and non-consistent models.

***Simulated emission spectra.*** In this section we show the simulated emission spectra for the series of HD 209458b models presented above.

For the model without a temperature inversion (Fig. 13) there is a very small difference in the flux between the EQ, NEQ and CNEQ cases. This is not surprising since the departure from chemical equilibrium is small; and consequently the induced change in the *P*–*T* profile is small.
<figure><img src="image_13.png" /><figcaption>**Fig. 12.** Same as Fig. 11 but for the HD 209458b model with a temperature inversion.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Fig. 13.** Emission spectra of the HD 209458b model without a temperature inversion, showing calculations based on the EQ calculation (blue), the NEQ calculation (red) and the CNEQ calculation (green); $K_{zz} = 10^{11} \text{ cm}^2 \text{ s}^{-1}$.</figcaption></figure>

However, the models including a temperature inversion do show important differences (Fig. 14). The NEQ spectrum shows a greater flux over most of the wavelength range compared with the EQ and CNEQ spectra. In particular there is a large increase in flux around 3.6 µm due to a large increase in the methane abundance due to transport-induced quenching. This is the opposite to what was found for the HD 189733b models where the flux was seen to decrease in the NEQ model. The primary difference between these models is that, in the HD 209458b model, at low pressures the temperature is increasing inversely with pressure, due to the presence of a temperature inversion.

On the other hand, the EQ and CNEQ spectra are remarkably similar, despite the fact the abundances of methane and ammonia are driven far from chemical equilibrium and the *P*–*T* profile is altered considerably. The increase in flux around 3.6 µm and at longer wavelengths present in the NEQ model are not apparent in the CNEQ model, removing the signatures of non-equilibrium chemistry. The cause of this will be explained in detail in the next section.

In this particular case, the overall effect of calculating the non-equilibrium chemistry consistently with a coupled temperature structure, rather than on a fixed *P*–*T* profile, is to reduce the influence of non-equilibrium chemistry on the emission spectrum, as found for the HD 189733b model, as the CNEQ spectrum tends back towards the EQ result.

<figure><img src="image_15.png" /><figcaption>**Fig. 14.** Same as Fig. 13 but for the HD 209458b model with a temperature inversion.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 15.** Atmosphere emission spectrum of the HD 189733b model with $K_{zz} = 10^{11} \text{ cm}^2 \text{ s}^{-1}$ showing calculations based on the EQ calculation (blue), the NEQ calculation (red) and the CNEQ calculation (green).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Fig. 16.** Same as Fig. 15 but for the $K_{zz} = 10^9 \text{ cm}^2 \text{ s}^{-1}$ model.</figcaption></figure>

### 3.3. Energy balance considerations

In this section we further investigate the differences resulting from a consistent treatment of calculating non-equilibrium chemical compositions.

Figures 15 and 16 show the thermal emission spectra of the atmosphere (not divided by $F_{star}$) for the HD 189733b models with $K_{zz} = 10^{11}$ and $10^9 \text{ cm}^2 \text{ s}^{-1}$, respectively. In both cases, the emission is lower at all wavelengths for the NEQ models, compared to both the EQ and CNEQ models. This discrepancy is larger for the model with stronger vertical mixing. By eye, one can already see that the total amount of energy emission of the atmosphere (i.e. the wavelength-integrated flux) is less in the NEQ model, compared with the EQ and CNEQ models.

Indeed, this is the case, as shown in Tables 2 and 3 which present the integrated top of atmosphere flux (not including the reflected component) and corresponding blackbody
<table><caption>**Table 2.** Integrated flux and corresponding blackbody temperatures for the HD 189733b model with $K_{zz} = 10^{11}$ cm$^2$ s$^{-1}$.</caption><thead><tr><th></th><th>Flux (kWm<sup>-2</sup>)</th><th>T<sub>BB</sub> (K)</th></tr></thead><tbody><tr><td>EQ</td><td>106.6</td><td>1171</td></tr><tr><td>NEQ</td><td>66.1</td><td>1039</td></tr><tr><td>CNEQ</td><td>106.4</td><td>1170</td></tr></tbody></table>

<table><caption>**Table 3.** Integrated flux and corresponding blackbody temperatures for the HD 189733b model with $K_{zz} = 10^{9}$ cm$^2$ s$^{-1}$.</caption><thead><tr><th></th><th>Flux (kWm<sup>-2</sup>)</th><th>T<sub>BB</sub> (K)</th></tr></thead><tbody><tr><td>EQ</td><td>106.6</td><td>1171</td></tr><tr><td>NEQ</td><td>94.6</td><td>1137</td></tr><tr><td>CNEQ</td><td>106.6</td><td>1171</td></tr></tbody></table>

temperatures for both $K_{zz}$ cases and for all chemistry models. In both $K_{zz}$ cases, the integrated flux for the EQ and CNEQ models agree well with each other, conserving the total amount of energy being emitted by the atmosphere. On the other hand, the NEQ models show strongly reduced integrated fluxes. The integrated flux in the NEQ models is ~38% and ~11% smaller than the EQ integrated flux for the $K_{zz} = 10^{11}$ cm$^2$ s$^{-1}$ and $K_{zz} = 10^{9}$ cm$^2$ s$^{-1}$ models, respectively.

These calculations show that the NEQ models do not conserve energy and the model atmosphere is not in a state of energy balance. The incoming energy (irradiation and internal heating) has not changed, only the chemical abundance profiles have changed, yet the atmosphere is emitting less energy. However, the CNEQ models do conserve the amount of energy being lost by the atmosphere, as the integrated flux is equivalent to that found for the EQ models.

To understand this further, we show the pressure level of peak emission (i.e. the photosphere) in Figs. 17 and 18 for the 3.6 µm and 8.0 µm Spitzer/IRAC channels, respectively. Here the pressure level of the photosphere is taken as the maximum of the contribution (or weighting) function (e.g. Knutson et al. 2009; Griffith et al. 1998). In both cases, the pressure level of the photosphere is shifted to lower pressures, and lower temperatures, in the NEQ model. This is a result of increased opacity due transport-induced quenching of CH$_4$ and NH$_3$.

Since the emission flux is strongly dependent on temperature, the shifting of the photosphere to lower temperatures results in a decreased emission in this wavelength band. Indeed, this occurs not just in this wavelength band but also at other points where CH$_4$/NH$_3$ absorb and is evident in the decreased integrated flux value previously shown for the NEQ models. For the CNEQ model, though the photosphere is shifted to lower pressures the temperature at this lower pressure is increased compared with the EQ/NEQ model. Here we see the P-T profile adapting to maintain radiative-convective equilibrium and energy balance in reaction to the changing chemical composition due to non-equilibrium chemistry.

The story is very similar for the models of HD 209458b. Figures 19 and 20 show the top of atmosphere emission spectrum for the HD 209458b models with and without a temperature inversion, respectively, in each case for the EQ, NEQ and CNEQ models. The model without a temperature inversion shows a negligible difference between all three cases. In the temperature inversion model, the NEQ case shows a greater flux at all wavelengths compared with the EQ case. In particular, there is a large increase in flux between 3 and 4 µm, which roughly

<figure><img src="image_18.png" /><figcaption>**Fig. 17.** Location of the peak emission in the 3.6 µm Spitzer band for the HD 189733b $K_{zz} = 10^{11}$ cm$^2$ s$^{-1}$ model with the EQ calculation (blue), NEQ calculation (red) and CNEQ calculation (green).</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Fig. 18.** As Fig. 17 for the Spitzer/IRAC 8.0 µm channel.</figcaption></figure>

<table><caption>**Table 4.** Integrated flux and corresponding blackbody temperatures for the HD 209458b temperature inversion model; $K_{zz} = 10^{11}$ cm$^2$ s$^{-1}$.</caption><thead><tr><th></th><th>Flux (kWm<sup>-2</sup>)</th><th>T<sub>BB</sub> (K)</th></tr></thead><tbody><tr><td>EQ</td><td>235.4</td><td>1427</td></tr><tr><td>NEQ</td><td>256.5</td><td>1458</td></tr><tr><td>CNEQ</td><td>235.4</td><td>1428</td></tr></tbody></table>

corresponds to the wavelength of peak emission. Methane dominates the absorption around 3.6 µm, which is increased in abundance by transport-induced quenching.

Again, by eye, it is possible to see that the NEQ models do not conserve the wavelength-integrated flux. Indeed, Table 4 shows the integrated flux of the atmosphere and the corresponding blackbody temperature for the temperature inversion models, where the integrated flux for the NEQ calculation is ~10% greater than the EQ model. On the other hand, the integrated flux for the CNEQ model is in excellent agreement with the EQ model. Likewise, Table 5 shows the same information for the HD 209458b model without a temperature inversion. In this case, the discrepancy between EQ and NEQ is smaller. The EQ and CNEQ models show consistent outgoing fluxes, however, the NEQ model shows a very small ~1% reduction in integrated flux.

For the model without a temperature inversion, the NEQ model shows a decreased integrated flux, similar to what was found for HD 189733b. However, the model with a temperature inversion shows an increased integrated flux for the NEQ model.
<figure><img src="image_20.png" /><figcaption>**Fig. 19.** Atmosphere emission spectrum of the HD 209458b model with a temperature inversion showing calculations based on the EQ calculation (blue), the NEQ calculation (red) and the CNEQ calculation (green); $K_{zz} = 10^{11} \text{ cm}^2 \text{ s}^{-1}$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>**Fig. 20.** Same as Fig. 19 but for the model without a temperature inversion; $K_{zz} = 10^{11} \text{ cm}^2 \text{ s}^{-1}$.</figcaption></figure>

<table><caption>**Table 5.** Integrated flux and corresponding blackbody temperatures for the HD 209458b model without temperature inversion, $K_{zz} = 10^{11} \text{ cm}^2 \text{ s}^{-1}$.</caption><thead><tr><th></th><th>Flux (kWm<sup>-2</sup>)</th><th><i>T</i><sub>BB</sub> (K)</th></tr></thead><tbody><tr><td>EQ</td><td>222.2</td><td>1407</td></tr><tr><td>NEQ</td><td>221.9</td><td>1406</td></tr><tr><td>CNEQ</td><td>222.7</td><td>1408</td></tr></tbody></table>

Figure 21 indicates the pressure level of the photosphere in the 3.6 µm *Spitzer*/IRAC band for the temperature inversion model. Similarly to the HD 189733b model, we see that transport-induced quenching pushes the photosphere to lower pressures by increasing the opacity. In this model however, where temperature is increasing with decreasing pressure, the photosphere moves to a higher temperature. This explains why we see an increased integrated flux in the NEQ model. However, for the CNEQ model, though the pressure level of the photosphere is still shifted to lower pressures the *P*–*T* profile has adapted so that the temperature at this level is now cooler, conserving the integrated flux to that of the EQ model, shown in Table 4.

For the model without a temperature inversion, a similar process occurs. However, because temperature decreases with altitude in this model, the new photosphere is both lower in pressure and lower in temperature resulting in a reduced emission flux. Therefore, to compensate, the *P*–*T* structure increases in temperature at this pressure level to increase the emission flux once again, and conserve energy balance.

This understanding allows us to explain why the simulated spectra from our CNEQ models are similar to the EQ model

<figure><img src="image_22.png" /><figcaption>**Fig. 21.** Location of the peak emission in the 3.6 µm Spitzer band for the HD 209458b model with a temperature inversion with the EQ calculation (blue), NEQ calculation (red) and CNEQ calculation (green); $K_{zz} = 10^{11} \text{ cm}^2 \text{ s}^{-1}$.</figcaption></figure>

spectra, reducing the impact of non-equilibrium chemistry. Our results show that it is not primarily the changing abundances due to non-equilibrium chemistry which affect the emission spectrum. It is instead the secondary effect of the non-equilibrium abundances shifting the location of the photosphere, by changing the opacity, which effects the calculated spectral flux, as the flux is now originating from a part of the atmosphere with a lower/higher temperature. This non-consistent method does not conserve energy balance in the atmosphere.

On the other hand, in the self-consistent non-equilibrium chemistry models, though the non-equilibrium abundances do also change the pressure level of the photosphere, the *P*–*T* profile adapts, by either heating up or cooling down, to maintain energy balance. Therefore, the resulting spectral flux for the CNEQ model shows a smaller discrepancy with the EQ model. Previous studies using non-consistent models may have overestimated the importance of non-equilibrium chemistry on the emission spectrum.

## 4. Conclusions

We have presented results based on a fully-consistent chemical kinetics model applied to the atmospheres of HD 189733b and HD 209458b.

Our simulations show that in cases of strong disequilibrium chemistry transport-induced quenching of absorbing species can induce changes in the *P*–*T* profile of up to 100 K. These temperature shifts can, in turn, have impacts on both the chemical abundances themselves and on the corresponding emission spectra. The chemical abundances are affected via two related processes: firstly, temperature changes at high pressures, where the chemistry remains in chemical equilibrium, induce new chemical equilibrium abundances, and secondly, the change in temperature shifts the quenching point which alters the quenched abundances at lower pressures.

For instance, in our model of HD 189733b we would conclude that CH$_4$ is more abundant than CO for the $K_{zz} = 10^{11} \text{ cm}^2 \text{ s}^{-1}$ case when not performing the calculation consistently. Instead, in the consistent approach, we find that CO is the most abundant carbon-bearing species. For the model of HD 209458b (with a temperature inversion) we find that the abundances of CH$_4$ and NH$_3$ are a factor of ~5 and ~3.5 lower in the consistent model compared with the non-consistent model, due to an increase in the temperature by more than 100 K in the deep atmosphere.
For both HD 189733b and HD 209458b models we find that consistent calculations of non-equilibrium chemistry reduce the overall impact of chemical disequilibrium on the emission spectrum. Our results show that in conventional chemical kinetics models, where the *P*-*T* profile is held fixed, the dominant mechanism for non-equilibrium chemistry to affect the emission spectrum is by changing the pressure level, and temperature, of the photosphere. The strong dependence of the emission flux on temperature will result in a very different simulated emission spectrum.

However, in our consistent model the *P*-*T* profile adapts to the new non-equilibrium chemical composition and to retain energy balance in the model atmosphere. The consequent temperature changes mitigate the effect of changing the location of the photosphere by either heating up or cooling down at the location of the new photosphere to preserve energy balance.

Based on these results, we urge caution when assessing the impact of non-equilibrium chemistry (transport-induced quenching and photochemistry) on the emission spectrum. Not including consistency between the chemical abundances and the temperature structure can lead to overestimates of the impact of non-equilibrium chemistry.

This work has only considered 1D (vertical) effects of non-equilibrium chemistry. Horizontal advection is expected to be very important in the atmospheres of tidally-locked exoplanets which possess very strong zonal wind velocities (Showman et al. 2009; Heng et al. 2011; Rauscher & Menou 2012; Mayne et al. 2014). In addition, these atmospheres can possess very large day-night temperature contrasts leading to large contrasts in horizontal chemical equilibrium abundances (Burrows et al. 2010; Kataria et al. 2016). Disequilibrium chemistry has already been suggested as a possible explanation to explain the discrepancies between the observed and model emission phase curves (Zellem et al. 2014). However, this work has shown that when performed consistently, transport-induced quenching has a smaller impact on the emission spectrum than previous studies suggest. It would therefore be very interesting and timely to study the process considered here including horizontal advection, by coupling a chemical kinetics scheme consistently to a 3D GCM.

**Acknowledgements.** This work is partly supported by the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement No. 247060-PEPS and grant No. 320478-TOFU). B.D. thanks the University of Exeter for support through a Ph.D. studentship. D.S.A. acknowledges support from the NASA Astrobiology Program through the Nexus for Exoplanet System Science, N.J.M. and J.G.'s contributions were in part funded by a Leverhulme Trust Research Project Grant, and in part by a University of Exeter College of Engineering, Mathematics and Physical Sciences studentship. This work used the DiRAC Complexity system, operated by the University of Leicester IT Services, which forms part of the STFC DiRAC HPC Facility. This equipment is funded by BIS National E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant ST/K0003259/1. DiRAC is part of the National E-Infrastructure. This work also used the University of Exeter Supercomputer, a DiRAC Facility jointly funded by STFC, the Large Facilities Capital Fund of BIS and the University of Exeter.

## References

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

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

* Agúndez, M., Venot, O., Selsis, F., & Iro, N. 2014b, ApJ, 781, 68

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

* Baraffe, I., Homeier, D., Allard, F., & Chabrier, G. 2015, A&A, 577, A42

* Barman, T. S., Hauschildt, P. H., & Allard, F. 2005, ApJ, 632, 1132

* Beaulieu, J. P., Kipping, D. M., Batista, V., et al. 2010, MNRAS, 409, 963

* Bendicho, P. F., & Bueno, J. T. 1995, ApJ

* Berta-Thompson, Z. K., Irwin, J., Charbonneau, D., et al. 2015, Nature, 527, 204

* Blecic, J., Harrington, J., & Bowman, M. O. 2016, ApJS, 225, 4

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

* Burrows, A., Rauscher, E., Spiegel, D. S., & Menou, K. 2010, ApJ, 719, 341

* Caffau, E., Ludwig, H.-G., Steffen, M., Freytag, B., & Bonifacio, P. 2011, Sol. Phys., 268, 255

* Diamond-Lowe, H., Stevenson, K. B., Bean, J. L., Line, M. R., & Fortney, J. J. 2014, ApJ, 796, 66

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

* Evans, T. M., Aigrain, S., Gibson, N., et al. 2015, MNRAS, 451, 680

* Fortney, J. J., Shabram, M., Showman, A. P., et al. 2010, ApJ, 709, 1396

* Gladstone, G. R., Allen, M., & Yung, Y. L. 1996, Icarus, 119, 1

* Gordon, S., & McBride, B. J. 1994, NASA Reference Publication, 1311

* Griffith, C. A., Yelle, R. V., & Marley, M. S. 1998, Science, 282, 2063

* Gustafsson, B., Edvardsson, B., Eriksson, K., et al. 2008, A&A, 486, 951

* Heng, K., Menou, K., & Phillipps, P. J. 2011, MNRAS, 413, 2380

* Heng, K. & Showman, A. P. 2015, Ann. Rev. Earth Planet. Sci., 43, 509

* Heng, K., & Tsai, S.-M. 2016, ApJ, 829, 104

* Henyey, L., Vardya, M. S., & Bodenheimer, P. 1965, ApJ, 142, 841

* Hindmarsh, A. 1983, Scientific Computing

* Hubeny, I., & Burrows, A. 2007, ApJ, 669, 1248

* Hubeny, I., & Mihalas, D. 2014, Theory of Stellar Atmospheres (Princeton University Press)

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

* Isaksen, I. S. A., Midtbø, K. H., Sunde, J., & Crutzen, P. J. 1977, Geophysica Norvegica, 31, 11

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

* Knutson, H. A., Charbonneau, D., Allen, L. E., Burrows, A., & Megeath, S. T. 2008, ApJ, 673, 526

* Knutson, H. A., Charbonneau, D., Cowan, N. B., et al. 2009, ApJ, 690, 822

* Kopparapu, R. K., Kasting, J. F., & Zahnle, K. J. 2012, ApJ, 745, 77

* Koskinen, T. T., Harris, M. J., Yelle, R. V., & Lavvas, P. 2013, Icarus, 226, 1678

* Lacis, A. A., & Oinas, V. 1991, J. Geophys. Res., 96, 9027

* Liang, M.-C., Parkinson, C. D., Lee, A. Y.-T., Yung, Y. L., & Seager, S. 2003, ApJ, 596, L247

* Liang, M.-C., Seager, S., Parkinson, C. D., Lee, A. Y.-T., & Yung, Y. L. 2004, ApJ, 605, L61

* Line, M. R., Liang, M. C., & Yung, Y. L. 2010, ApJ, 717, 496

* Madhusudhan, N., & Seager, S. 2011, ApJ, 729, 41

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

* McBride, B. J., Gordon, S., & Reno, M. A. 1993, NASA Technical Memorandum, 4513

* Michelangeli, D. V., Allen, M., Yung, Y. L., et al. 1992, J. Geophys. Res.: Atmospheres, 97, 865

* Miguel, Y., & Kaltenegger, L. 2014, ApJ, 780, 166

* Moses, J. I. 2014, Philos. Trans. Roy. Soc. London Series A, 372, 20130073

* Moses, J. I., Madhusudhan, N., Visscher, C., & Freedman, R. S. 2013, ApJ, 763, 25

* Moses, J. I., Visscher, C., Fortney, J. J., et al. 2011, ApJ, 737, 15

* Parmentier, V., Showman, A. P., & Lian, Y. 2013, A&A, 558, A91

* Polichtchouk, I., Cho, J. Y.-K., Watkins, C., et al. 2014, Icarus, 229, 355

* Rauscher, E., & Menou, K. 2012, ApJ, 750, 96

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

* Sing, D. K., Fortney, J. J., Nikolov, N., et al. 2016, Nature, 529, 59

* Southworth, J. 2010, MNRAS, 408, 1689

* Spiegel, D. S., & Burrows, A. 2010, ApJ, 722, 871

* Thomas, G. E., & Stammes, K. 1999, Radiative Transfer in the Atmosphere and Ocean (Cambridge University Press)

* Tremblin, P., Amundsen, D. S., Chabrier, G., et al. 2016, ApJ, 817, L19

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

* Venot, O., Agúndez, M., Selsis, F., Tessenyi, M., & Iro, N. 2014, A&A, 562, A51

* Venot, O., Hébrard, E., Agúndez, M., Decin, L., & Bounaceur, R. 2015, A&A, 577, A33

* Venot, O., Hébrard, E., Agúndez, M., et al. 2012, A&A, 546, A43

* Wayne, R. 1991, Chemistry of Atmospheres, 2nd edn. (Oxford University Press, New York: Oxford Science Publications)

* Yelle, R. V. 2004, Icarus, 170, 167

* Zahnle, K., Marley, M. S., Freedman, R. S., Lodders, K., & Fortney, J. J. 2009, ApJ, 701, L20

* Zahnle, K. J., & Marley, M. S. 2014, ApJ, 797, 41

* Zellem, R. T., Lewis, N. K., Knutson, H. A., et al. 2014, ApJ, 790, 53
# Appendix A: Model testing

In this appendix section we perform several comparisons with previously published works to benchmark both the radiative-convective and the chemistry modules of ATMO.

## A.1. Radiative-convective Scheme

In this section we examine the radiative-convective equilibrium component of ATMO by reproducing the P-T profiles of Iro et al. (2005) and Barman et al. (2005; hereafter I05 and B05) for HD 209458b (Fig. A.1). In both cases we match the planetary parameters (i.e. internal temperature, gravity etc.) as closely as possible to those used in the original studies. Remaining differences between our model and those of I05 and B05 stem from different opacity sources, the stellar irradiation spectrum, and the radiative-transfer and chemistry schemes.

I05 models an atmosphere with an internal temperature, $T_{\text{int}}$, of 100 K and 300 K and a planet-wide heat redistribution parameter, i.e. with an incident flux reduction factor of $\alpha = 0.25$ which accounts for emission from both the dayside and nightside of the planet.

Generally our P-T profile agrees very well with that of I05, though our model is ~100 K warmer throughout most of the modeled pressure range. This very likely stems from the different stellar irradiation used in the two models. I05 use a solar irradiation in their model, however we use the Kurucz simulated spectrum of HD 209458 which is several hundred Kelvin hotter than the Sun. Differences also likely stem from the treatment and source of opacities.

We also compare with the $\alpha = 0.25$ model of B05. Note that we adjusted the internal temperature parameter until the profiles reached the same adiabat at depth since we do not know the exact value of $T_{\text{int}}$ used in the B05 model. We found that $T_{\text{int}} = 380$ K matched the B05 profile well. We find that the ATMO profile is consistently lower in temperature with a maximum discrepancy at ~0.1 bar. Aside from the differences in opacity sources (we use more up-to-date sources with increased infrared opacity) between our model and the B05 model, the calculation of equilibrium chemistry is significantly different. B05 accounts for the “rainout” process whereby condensation of species removes material from the gas-phase, depleting the elemental abundance of those species at lower pressures. We do not include this process in ATMO, potentially leading to differences in chemical composition and hence opacity.

The ATMO models of HD 209458b are bracketed by those of I05 and B05. There are significant differences in the opacity sources, irradiation, chemical and radiative schemes between all three models. Considering these significant differences, the three models generally agree on the qualitative temperature structure, with quantitative differences which should be expected due to the differences in model inputs and schemes.

## A.2. Chemistry schemes

### A.2.1. Gibbs minimisation scheme

In this section we test the Gibbs minimisation scheme by comparing our chemical equilibrium abundances with those derived using the Burrows & Sharp (1999) and Heng & Tsai (2016) analytical formulae. Figures A.2 and A.3 show the abundances of CO, H₂O, CH₄, NH₃, N₂ and HCN derived using the Gibbs energy minimisation scheme and the two analytical methods for

<figure><img src="image_23.png" /><figcaption>Fig. A.1. A comparison of the P-T profiles for HD 209458b with ATMO results in solid lines and the models of Iro et al. (2005) and Barman et al. (2005) in dashed lines.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Fig. A.2. Chemical equilibrium abundances derived using the Gibbs energy minimisation scheme (coloured solid lines), the Burrows & Sharp (1999) analytical formula (coloured dashed lines) and the Heng & Tsai (2016) analytical formula (black lines) for a T = 1000 K isothermal profile. Note that HCN is not included in the Burrows & Sharp (1999) formula.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Fig. A.3. As Fig. A.2 for a T = 2000 K isothermal profile.</figcaption></figure>

$T = 1000$ K and 2000 K isothermal profiles, respectively; note that HCN is not included in the Burrows & Sharp (1999) scheme.

Our Gibbs minimisation scheme agrees very well with the analytical formula of Heng & Tsai (2016) for both temperature profiles. The agreement with the Burrows & Sharp (1999)
<figure><img src="image_26.png" /><figcaption>**Fig. A.4.** Chemical evolution of (from top left, clockwise) H, OH, NNH and C₂H₆ with different approaches to the reversal of certain chemical reactions (see Venot et al. 2012, for details) in the dashed and solid lines; to be compared with Venot et al. (2012, Fig. 2).</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Fig. A.5.** Chemical abundances calculated for HD 209458b with equilibrium abundances (dashed) and non-equilibrium abundances (solid) using the Moses et al. (2011) pressure-temperature profile (their Fig. 2, dayside average profile); to be compared with Venot et al. (2012, Fig. 5) and Moses et al. (2011, Figs. 5 and 6).</figcaption></figure>

formula is also very good. At low pressures in the T = 2000 K model the Gibbs minimisation abundances deviate at low pressures due to molecular hydrogen dissociation which is not included in the analytical formula.

### A.2.2. Chemical kinetics scheme

In this section we test the chemical kinetics scheme implemented in our model.

Firstly we test our chemical kinetics scheme in isolation (i.e. only calculation of the chemical rates) using a 0D (box) model by reproducing the chemical evolution plots of Venot et al. (2012, their Fig. 2) for H, OH, NNH and C₂H₆ (Fig. A.4). The evolution and final abundances of these species are, by eye, identical, validating the implementation of the Venot et al. (2012) chemical network in our model.

Secondly, we reproduce the 1D chemical abundances profiles for HD 209458b included in Venot et al. (2012) using the same model inputs; pressure-temperature profile, chemical network, UV spectrum and eddy diffusion coefficient profile. Venot et al. (2012) used the same model inputs as Moses et al. (2011), and compared the two different chemical networks. This allows us to compare with Venot et al. (2012), whose chemical network we use here, testing only differences in the model implementation, and also with Moses et al. (2011), who use both a different model and a different chemical network. The primary difference between the ATMO scheme and that of Venot et al. (2012) is the calculation of the UV flux. We calculate the UV flux as a function of pressure using the radiative transfer scheme as described in Sect. 2.2, whereas Venot et al. (2012) use the scheme of Isaksen et al. (1977), which has a more simple treatment of scattering.

Figure A.5 shows the chemical abundances calculated with ATMO, abundances calculated from the chemical kinetics scheme including vertical mixing and photochemistry. Overall, there is a very good agreement between our model and that of Venot et al. (2012). In the deep atmosphere we retain chemical equilibrium, due to the very fast chemical timescale due to high pressures and temperatures. The abundances also agree very well in the middle regions of the atmosphere, with accurate reproductions of the methane and ammonia quenching points.
There are, however, some disagreements for the abundances at low pressures. As noted previously, one of the most important differences between our model and that of Venot et al. (2012) is the calculation of the UV flux. In particular, we find a production of methane at low pressures, when photochemistry is included, which is not seen in the Venot et al. (2012) model. Interestingly, this feature is present in the Moses et al. (2011) model (see their Fig. 5). The Moses et al. (2011) model uses Feautrier radiative-transfer method (Michelangeli et al. 1992), taking into account multiple-scattering by H₂ and He. Venot et al. (2012) note that abundances at low pressures are sensitive to the treatment of Rayleigh scattering.

In summary, our Gibbs minimisation scheme and chemical kinetics scheme, using the Venot et al. (2012) chemical network, agree very well for chemical equilibrium. Our chemical equilibrium results agree well with Venot et al. (2012) and Moses et al. (2011), and with the analytical formula of Burrows & Sharp (1999). We also find a reasonably good agreement with Venot et al. (2012) when including vertical mixing and photochemistry, reproducing the quench points of the various species well. However, we find some discrepancies at low pressures where photochemistry is important, which we attribute to differences in the calculation of the UV flux.