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

# The 3D Thermal, Dynamical, and Chemical Structure of the Atmosphere of HD 189733b: Implications of Wind-driven Chemistry for the Emission Phase Curve

Benjamin Drummond¹, Nathan J. Mayne¹, James Manners¹,², Isabelle Baraffe¹,³, Jayesh Goyal¹, Pascal Tremblin⁴, David K. Sing⁵, and Krisztian Kohary¹

¹ Astrophysics Group, University of Exeter, EX4 2QL, Exeter, UK; b.drummond@exeter.ac.uk

² Met Office, Exeter, EX1 3PB, UK

³ Univ. Lyon, Ens de Lyon, Univ Lyon1, CNRS, CRAL, UMR5574, F-69007, Lyon, France

⁴ Maison de la simulation, CEA, CNRS, Univ. Paris-Sud, UVSQ, Universit Paris-Saclay, F-91191 Gif-Sur-Yvette, France

⁵ Department of Earth and Planetary Sciences, Johns Hopkins University, Baltimore, MD, USA

Received 2018 July 23; revised 2018 October 5; accepted 2018 October 22; published 2018 December 7

## Abstract

In this paper we present 3D atmospheric simulations of the hot Jupiter HD 189733b under two different scenarios: local chemical equilibrium and including advection of the chemistry by the resolved wind. Our model consistently couples the treatment of dynamics, radiative transfer, and chemistry, completing the feedback cycle between these three important processes. The effect of wind-driven advection on the chemical composition is qualitatively similar to our previous results for the warmer atmosphere of HD 209458b, found using the same model. However, we find more significant alterations to both the thermal and dynamical structure for the cooler atmosphere of HD 189733b, with changes in both the temperature and wind velocities reaching ~10%. We also present the contribution function, diagnosed from our simulations, and show that wind-driven chemistry has a significant impact on its 3D structure, particularly for regions where methane is an important absorber. Finally, we present emission phase curves from our simulations and show the significant effect of wind-driven chemistry on the thermal emission, particularly within the 3.6 µm Spitzer/IRAC channel.

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

## 1. Introduction

Interpreting observations of exoplanet atmospheres requires comparison with theoretical models in order to infer the properties of the atmosphere, such as the thermal structure and the chemical composition. Whether this is done by a retrieval method (e.g., Irwin et al. 2008; Madhusudhan & Seager 2009; Waldmann et al. 2015) or by comparison with forward models (e.g., Barman et al. 2005; Fortney et al. 2005; Madhusudhan et al. 2011; Moses et al. 2011; Goyal et al. 2018), these theoretical tools are most often in the form of a 1D model. However, the results of 3D radiative-hydrodynamic codes (Showman & Guillot 2002; Dobbs-Dixon & Lin 2008; Menou & Rauscher 2009; Showman et al. 2009; Heng et al. 2011; Mayne et al. 2014a; Amundsen et al. 2016) clearly show that the atmospheres of highly irradiated tidally locked planets are far from being horizontally symmetric.

For tidally locked hot Jupiters, intense stellar irradiation heats the dayside to thousands of kelvin, while the nightside is significantly cooler, driving fast zonal winds (~$10^3$ m s⁻¹) that redistribute heat around the planet (e.g., Showman et al. 2009; Dobbs-Dixon & Agol 2013; Amundsen et al. 2016). High-resolution transmission spectra revealing Doppler shifting of atomic absorption lines (Snellen et al. 2010; Louden & Wheatley 2015), as well as measurements of the emission as a function of orbital phase (e.g., Knutson et al. 2012; Zellem et al. 2014), have since confirmed these theoretical predictions. However, Dang et al. (2018) recently reported a westward shift in the location of the hot spot, which contradicts the eastward shifts unanimously predicted by current 3D models.

Recently, Dobbs-Dixon & Cowan (2017) demonstrated the importance of considering the 3D structure of the atmosphere when interpreting the emission phase curve. They showed that

the contribution function, the peak of which indicates the pressure level of the photosphere, varies with longitude. In general, they found that the photosphere typically lies at lower pressures on the dayside compared with the nightside, with the conclusion that the observed emission originates from different pressure levels of the atmosphere throughout the orbit of the planet.

Most 3D models of hot exoplanet atmospheres to date share a common limitation: the assumption of a cloud-free atmosphere that is in local chemical equilibrium. However, many investigations with 1D models have shown the importance of vertical transport and photochemistry that drive the chemistry away from equilibrium (e.g., Line et al. 2010; Moses et al. 2011; Venot et al. 2012; Zahnle & Marley 2014; Drummond et al. 2016; Tsai et al. 2017). Large dayside–nightside temperature contrasts naturally lead to large dayside–nightside contrasts in the chemical equilibrium composition. In addition, the properties of clouds have recently been investigated using 3D models of various complexity (Lee et al. 2016; Parmentier et al. 2016; Roman & Rauscher 2017; Lines et al. 2018b).

Cooper & Showman (2006) used simplified chemistry and radiative transfer schemes with a 3D hydrodynamics code and found that vertical transport from high-pressure regions, which are horizontally uniform, leads to homogenization of the chemistry at lower pressures. Agúndez et al. (2014) used a “psuedo-2D” code (in practice a 1D column with a time-varying temperature profile) and found that horizontal (zonal) transport has a more important effect than vertical transport and contaminates the nightside chemistry with that of the hot dayside.

Recently, we coupled the simplified chemical scheme of Cooper & Showman (2006) to our own 3D model, including a
state-of-the-art radiative transfer scheme (Amundsen et al. 2016), completing the feedback cycle between the dynamics, radiation, and gas-phase chemistry (Drummond et al. 2018b). We found that a combination of horizontal (zonal and meridional) and vertical transport ultimately determines the abundance of methane in the atmosphere of HD 209458b. The overall effect of 3D transport is to increase the mole fraction of methane, leading to more prominent methane absorption features in both the simulated transmission and emission spectra, compared with the chemical equilibrium case.

In this paper we present results from the same 3D radiative-hydrodynamics code, coupled with the same gas-phase chemical relaxation scheme (Cooper & Showman 2006; Drummond et al. 2018b), applied to the specific case of HD 189733b. We compare our results with our previous simulations of the warmer atmosphere of HD 209458b and present the response of the atmospheric temperature and circulation to wind-driven chemistry, in Sections 3 and 4. We then present contribution functions diagnosed from our 3D simulations and compare with the results of Dobbs-Dixon & Cowan (2017), where we find significant qualitative differences, in Section 5. Finally, we consider the effect of departures of the chemistry from chemical equilibrium on the emission phase curve in Section 6.

## 2. Model Description

### 2.1. The Unified Model (UM)

We use the Met Office UM to simulate the atmosphere of HD 189733b. The UM has been used in previous works to simulate the exoplanet atmospheres of HD 209458b (Mayne et al. 2014a, 2017; Amundsen et al. 2016; Drummond et al. 2018b; Lines et al. 2018b), GJ 1214b (Drummond et al. 2018a), and Proxima Centauri b (Boutle et al. 2017; Lewis et al. 2018).

The dynamical core of the UM (ENDGame) solves the deep-atmosphere, non-hydrostatic Navier-Stokes equations (Mayne et al. 2014b, 2017; Wood et al. 2014). The heating rates are computed using the open-source SOCRATES<fnref n="6" /> radiative transfer scheme (Edwards 1996; Edwards & Slingo 1996), which has been updated and tested for the high-temperature, hydrogen-dominated conditions of hot Jupiter atmospheres (Amundsen et al. 2014, 2016, 2017). The chemical composition is derived using the same methods as described in Drummond et al. (2018b): an analytical formula to derive the chemical equilibrium abundances and a chemical relaxation scheme based on Cooper & Showman (2006).

The chemical relaxation method describes the time dependence of a chemical species by relaxing the mole fraction toward a prescribed equilibrium profile on some chemical timescale (Cooper & Showman 2006; Tsai et al. 2018; Drummond et al. 2018b). The equilibrium profile is taken to be the mole fraction corresponding to chemical equilibrium, while the timescale is estimated or parameterized based on the elementary reactions involved in the interconversion of the relevant chemical species. The accuracy of the chemical relaxation method is primarily determined by the accuracy of the chosen chemical timescale (Tsai et al. 2018).

In Appendix C we compare the results of the chemical relaxation method against a full chemical kinetics method

<table><caption>**Table 1**
Key Model Parameters</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Mass, M<sub>P</sub></td><td>2.18 × 10<sup>27</sup> kg</td></tr><tr><td>Radius, R<sub>P</sub></td><td>8.05 × 10<sup>7</sup> m</td></tr><tr><td>Semimajor axis, a</td><td>0.031 au</td></tr><tr><td>Surface gravity, g<sub>surf</sub></td><td>22.49 ms<sup>-2</sup></td></tr><tr><td>Intrinsic temperature, T<sub>int</sub></td><td>100 K</td></tr><tr><td>Lower boundary pressure, P<sub>lower</sub></td><td>2 × 10<sup>7</sup> Pa</td></tr><tr><td>Rotation rate, Ω</td><td>3.28 × 10<sup>-5</sup> s<sup>-1</sup></td></tr><tr><td>Specific heat capacity, c<sub>P</sub></td><td>13 kJ kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Specific gas constant, R</td><td>3556.8 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Vertical damping coefficient, R<sub>w</sub></td><td>0.15</td></tr><tr><td>Vertical damping extent, η<sub>s</sub></td><td>0.75</td></tr><tr><td>Horizontal resolution</td><td>144 × 90</td></tr><tr><td>Vertical resolution</td><td>66</td></tr><tr><td>Dynamical time step</td><td>30 s</td></tr><tr><td>Radiative time step</td><td>150 s</td></tr></tbody></table>

**Note.** Stellar and planetary parameters for HD 189733b adapted from Southworth (2010).

within the framework of a 1D model and find a good agreement. We also demonstrate the sensitivity of our results to the chosen chemical timescale in Appendix D by artificially increasing/decreasing the value of the chemical timescale by a factor of 10. We find that varying the timescale does not significantly affect the final mole fraction of methane over a large pressure range, and therefore uncertainty in the precise value of the timescale is unlikely to affect our conclusions.

Tracer advection is handled using the extensively tested semi-implicit semi-Lagrangian scheme of the UM (Wood et al. 2014). In Appendix E we test the conservation of the global mass of elemental carbon and oxygen for our simulations. The model conserves the global mass of these elements to within better than 99.9%.

### 2.2. Model Parameters and Setup

We use the planetary and stellar parameters of HD 189733b from Southworth (2010), which we summarize in Table 1. The intrinsic temperature T<sub>int</sub> is the blackbody temperature of the net intrinsic flux at the lower boundary, accounting for heat escaping from the interior (Amundsen et al. 2016). A full description of the vertical damping terms (R<sub>w</sub> and η<sub>s</sub>) can be found in Mayne et al. (2014a).

For the stellar irradiation spectrum we use the Kurucz spectrum for HD 189733.<fnref n="7" /> We assume solar elemental abundances from Asplund et al. (2009). For a complete description of the included opacities see Amundsen et al. (2016). For the main model integration the radiation spectrum is divided into 32 bands for the heating rate calculations (Amundsen et al. 2016). However, for the simulated observation calculations we restart the model, starting from the state at 1000 days, and run for a small number of time steps using 500 spectral bands, to obtain a higher spectral resolution (as done in Boutle et al. 2017; Drummond et al. 2018a, 2018b; Lines et al. 2018b).

The model setup is broadly the same as in our previous coupled radiative transfer simulations of hot Jupiters (Amundsen et al. 2016; Drummond et al. 2018b). As in Drummond et al. (2018b),

<footnote n="6">⁶ https://code.metoffice.gov.uk/trac/socrates</footnote>

<footnote n="7">⁷ http://kurucz.harvard.edu/stars.html</footnote>
we perform two simulations that are identical except that one assumes local chemical equilibrium while the other includes the effect of advection, due to the large-scale resolved wind, and chemical evolution via the chemical relaxation method (Cooper & Showman 2006). We refer to these two simulations as the “equilibrium” and “relaxation” simulations, respectively.

The model is initialized with zero winds and with a horizontally uniform thermal profile. For the latter we use a radiative–convective equilibrium temperature profile from the 1D ATMO model using the same stellar and planetary parameters (Drummond et al. 2016). The simulations are integrated for 1000 Earth days (hereafter days refers to Earth days), by which point the maximum wind velocities have ceased to evolve; however, we note that the deep atmosphere ($P \gtrsim 10^6$ Pa) has not reached a steady state (Amundsen et al. 2016; Mayne et al. 2017), though this does not affect lower pressures. Unless otherwise stated, all figures show the simulations at 1000 days. The total axial angular momentum is conserved to within 99% for the simulations presented here.

We note that throughout the results sections of this paper we regularly calculate the difference of several quantities (e.g., temperature) between the relaxation and equilibrium simulations, to aid interpretation. Unless otherwise stated, the difference we refer to is the absolute difference, $A^{\text{relaxation}} - A^{\text{equilibrium}}$. In some cases, where more useful, we consider the relative difference, $(A^{\text{relaxation}} - A^{\text{equilibrium}})/A^{\text{relaxation}}$, instead.

### 3. Wind-driven Chemistry in the Atmosphere of HD 189733b

In this section we present the thermal, dynamical, and chemical structure of the atmosphere. We show the thermal and dynamical structure of the equilibrium simulation and compare the chemical structure of the equilibrium and relaxation simulations.

Figure 1 shows the zonal-mean zonal wind velocity, meridional wind velocity at $P = 5 \times 10^4$ Pa, and temperature structure for the chemical equilibrium simulation of HD 189733b. The wind velocities and temperatures are qualitatively similar to previous 3D simulations of HD 189733b (Showman et al. 2009; Dobbs-Dixon & Agol 2013; Rauscher & Menou 2013; Kataria et al. 2016). The circulation is characterized by an equatorial jet with a maximum wind velocity of ~6 km s⁻¹ with slower retrograde circulation at higher latitudes. The dayside–nightside temperature contrast (typically hundreds of kelvin) increases with decreasing pressure while the hot spot moves closer to the substellar point, due to the decreasing radiative timescale with decreasing pressure (e.g., Iro et al. 2005).

Figure 2 shows the mole fractions of carbon monoxide, water, and methane for the equilibrium simulation. The distribution of these species exactly traces the temperature structure (compare with Figure 1) since in chemical equilibrium the composition is entirely dependent on the local pressure and temperature, for a given mix of elements.

Carbon monoxide is more abundant than methane almost everywhere in the modeled domain, though the methane mole fraction varies significantly between ~$10^{-4}$ on the nightside and ~$10^{-10}$ on the warmer dayside. In the midlatitude region of the nightside (where the atmosphere is the coolest) methane becomes more abundant than carbon monoxide; there is a corresponding increase in the mole fraction of water.

For the relaxation simulation, shown in Figure 3, which includes advection due to the resolved wind, the chemistry is

homogenized both vertically and horizontally over a large pressure range ($P \lesssim 10^5$ Pa), similar to our results for HD 209458b (Drummond et al. 2018b). Overall this leads to a much larger dayside, equatorial methane abundance compared with the equilibrium simulation.

In Drummond et al. (2018b) we identified a mechanism whereby meridional transport leads to an increase in the equatorial methane abundance, compared with chemical equilibrium, for the case of HD 209458b. Using a simple tracer experiment, we demonstrated that mass is transported from higher latitudes to the equatorial region. Since the atmosphere is typically cooler at higher latitudes, compared with the equator, the equilibrium abundance of methane is larger. The net result is the transport of mass that has a relatively high methane fraction toward the equator. Figure 1 (top right panel) clearly shows that significantly large regions of equatorward flow exist. This 3D effect cannot be captured by 1D (e.g., Moses et al. 2011; Venot et al. 2012; Drummond et al. 2016) or 2D (Agúndez et al. 2014) models.

Figure 4 shows vertical profiles of the mole fractions of carbon monoxide, water, and methane for a number of longitude points around the equator. Methane becomes vertically quenched at $P \sim 10^5$ Pa with a mole fraction of ~4 × $10^{-5}$. Importantly, we find the same meridional transport effect (the sharp increase in methane abundance between $1 \times 10^5$ and $2 \times 10^5$ Pa) as identified in our simulations of HD 209458b (Drummond et al. 2018b), which is due to meridional transport from the midlatitudes toward the equator.

Overall, we find qualitatively the same trends as we previously found for HD 209458b (Drummond et al. 2018b). The main quantitative differences are the larger mole fractions of methane in the cooler atmosphere of HD 189733b.

### 4. Thermal and Dynamical Response of the Atmosphere

In this section we consider the thermal and dynamical response of the atmosphere to the changes in the local chemical composition due to wind-driven chemistry, by comparing the temperature structures of the equilibrium and relaxation simulations. The chemical composition and the thermal structure are linked via the radiative heating rates.

#### 4.1. Temperature Response

Figure 5 shows the absolute temperature difference between the equilibrium and relaxation simulations in various different perspectives. A positive difference indicates a higher temperature in the relaxation simulation.

For pressures less than ~$10^4$ Pa the atmosphere is generally cooler on the dayside and warmer on the nightside, for the relaxation simulation. In particular, there is a significant temperature increase of >80 K in the midlatitude regions of the nightside. For larger pressures, particularly between $10^4$ and $10^5$ Pa, the temperature increases significantly in the equatorial region, within the equatorial jet, for all longitudes.

We compared estimates of both the radiative ($\tau_{\text{rad}}$) and dynamical ($\tau_{\text{dyn}}$) timescales for the equilibrium simulation. For cases where $\tau_{\text{rad}} < \tau_{\text{dyn}}$ the atmosphere is expected to be in radiative equilibrium (the atmosphere is radiatively driven), while for $\tau_{\text{rad}} > \tau_{\text{dyn}}$ advection of heat is expected to be important (the atmosphere is dynamically driven). We used Equation (10) of Showman & Guillot (2002) to estimate $\tau_{\text{rad}}$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Top left panel: zonal-mean zonal wind velocity (color scale and black contours) temporally averaged between 800 and 1000 days. Top right panel: meridional wind velocity (color scale and black contours) on the 5 × 10⁴ Pa isobaric surface. Middle row: atmospheric temperature (color scale and black contours) on the 1 × 10² Pa isobaric surface (left) and on the 5 × 10⁴ Pa isobaric surface (right) with horizontal wind velocity vectors (black arrows). Bottom left: area-weighted meridional mean (between ±20° latitude) of the atmospheric temperature (color scale and black contours) with vertical-zonal wind velocity vectors (black arrows). Bottom right panel: atmospheric temperature (color scale and black contours) at a longitude of 0° with meridional-vertical wind vectors (black arrows).</figcaption></figure>

and for $\tau_{dyn}$ assumed (for the upper atmosphere)

$$ \tau_{dyn} \sim \frac{R_P}{u} \sim \frac{H}{w} \sim 10^4 \text{ s}, \qquad (1) $$

where $R_P \sim 10^7$ m is the planet radius, $H \sim 10^5$ m is the vertical scale height, and $u \sim 10^3$ m s⁻¹ and $w \sim 10^1$ m s⁻¹ are the horizontal and vertical wind velocities, respectively.

Comparing the timescales, we find that the atmosphere is expected to be radiatively driven ($\tau_{rad} < \tau_{dyn}$) for $P < 10^4$ Pa and dynamically driven ($\tau_{rad} > \tau_{dyn}$) for $P > 10^4$ Pa. This transition gives a hint as to why we find different trends in thermal responses above and below $P \sim 10^4$ Pa.

To unpick the combined effects of the three molecules (methane, carbon monoxide, and water), we performed an additional test simulation (not shown) that is identical to the relaxation simulation, except that the mole fractions of water and carbon monoxide used in the radiative transfer calculations (i.e., to calculate the heating rate) correspond to chemical

equilibrium, isolating the effect of methane. The resulting temperature and wind structure is almost identical to the nominal relaxation simulation, indicating that methane is the most important driver of these temperature changes, with water and carbon monoxide making a less important contribution. This is not surprising given the relatively small changes in water and carbon monoxide abundances between the equilibrium and relaxation simulations, despite both species being more abundant than methane.

## 4.2. Top-of-atmosphere (TOA) Radiative Flux

The energy balance of a close-in tidally locked atmosphere is dominated by stellar (shortwave) heating of the dayside atmosphere due to irradiation by the host star and thermal (longwave) cooling that occurs throughout the atmosphere. Throughout the rest of this paper we will refer to the stellar and thermal components as shortwave and longwave, respectively. The dominant source of heating for the nightside atmosphere is
<figure><img src="image_3.png" /><figcaption>Figure 2. Chemical equilibrium mole fractions (color scale and black contours) and wind velocity vectors (black arrows) of carbon monoxide (left column), water (middle column), and methane (right column) at a longitude of 0° (top row), on the 5 × 10⁴ Pa isobaric surface (middle row) and an area-weighted meridional mean (±20° latitude) around the equator (bottom row).</figcaption></figure>

advection of heat from the irradiated dayside. To obtain a global view of the energy balance, we consider the outgoing longwave radiative flux at the TOA, spectrally integrated over all wavelengths, shown in Figure 6.

The TOA longwave flux qualitatively traces the temperature structure at $P \sim 5 \times 10^4$ Pa (see Figure 1), with the largest emission occurring eastward of the substellar point, where the atmosphere is the warmest. Comparing the equilibrium and relaxation simulations, it is apparent that the dayside TOA flux is decreased in the relaxation simulation, while on the nightside it is increased. The greatest change occurs in the midlatitude regions of the nightside, where the TOA flux is increased by around 25%.

As the primary source of energy for the nightside atmosphere is advection of heat from the irradiated dayside, this indicates an overall increased efficiency of dayside-to-nightside heat transport. This trend agrees with the changes in the temperature explored in the previous section, where the atmosphere was generally found to be warmer on the nightside but cooler on the dayside, for the relaxation simulation, for pressure less than $10^4$ Pa. In the following sections we will show that this is due to changes in the radiative heating rates and an increase in

the speed of the equatorial jet, which themselves are not independent of each other.

### 4.3. Heating Rates

To further understand the temperature response of the atmosphere, we consider the radiative heating rates, which are shown in Figure 7 for the equatorial region. We show separately the shortwave heating (positive heating rate) and the longwave cooling (negative heating rate), as well as the net heating rate, which is the sum of the shortwave and longwave components.

We note that in Figure 7 we show the heating rate $\mathcal{H}$ in units of $\mathrm{W m}^{-3}$, the rate of change of energy per unit volume (Amundsen et al. 2014, Equation (17)). This is related to $\mathcal{H}$ in units of $\mathrm{K s}^{-1}$, the rate of change of temperature, by

$$ \mathcal{H}(\mathrm{W m}^{-3}) = \rho c_p \mathcal{H}(\mathrm{K s}^{-1}), \qquad (2) $$

where $\rho$ is the mass density and $c_p$ is the specific heat capacity. In these simulations $c_p$ is a global constant, shown in Table 1.

Figure 7 (left column) shows the heating rates for the equilibrium simulation, which are qualitatively very similar to those from the relaxation simulation (not shown). Naturally, the
<figure><img src="image_4.png" /><figcaption>Figure 3. Same as Figure 2, but for the relaxation simulation.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 4. Vertical profiles of the carbon monoxide (red), water (blue), and methane (black/gray) mole fractions for a number of columns equally spaced in longitude around the equator (at 0° latitude) for the chemical equilibrium simulation (dashed) and relaxation simulation (solid).</figcaption></figure>

shortwave heating (top left panel) is restricted to the dayside atmosphere and peaks at $P \sim 4 \times 10^4$ Pa at the substellar point. Longwave cooling (middle left panel), which occurs

across both the dayside and nightside, peaks at similar pressures, although shifted in longitude eastward of the substellar point. The net heating rate (bottom left panel) shows an overall positive heating rate (i.e., net heating) for the dayside and an overall negative heating rate (i.e., net cooling) for the nightside.

Figure 7 also shows the absolute difference in the heating rates between the relaxation and equilibrium simulations. On the dayside, there is a clear shift to lower pressures for the shortwave heating (top right panel), as stellar flux is absorbed higher in the atmosphere owing to the enhanced methane abundance. This results in a relative heating for $P < 2 \times 10^4$ Pa and a relative cooling for larger pressures, in the shortwave.

In the longwave (middle right panel), the peak of the cooling is also shifted to lower pressures, particularly for the region eastward of the substellar point, where the increase in the methane abundance is most significant. In Figure 7 this is indicated by a negative difference in the longwave heating rate (i.e., more cooling) for $P < 10^4$ Pa and a positive difference (i.e., less cooling) for $P > 10^4$ Pa. The change in the net heating rate (bottom right panel) shows that the shortwave effects are dominant on the dayside and the longwave effects are dominant on the nightside.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Temperature difference between the relaxation and equilibrium simulations at a longitude of 0° (top), on the $P = 5 \times 10^4$ Pa isobaric surface (middle), and an area-weighted meridional mean between ±20° latitude (bottom). A positive difference indicates a higher temperature in the relaxation simulation.</figcaption></figure>

The overall effect of the changes to the shortwave and longwave heating rates is to increase the shortwave heating on the dayside and increase the longwave cooling on the nightside, for the $10^4$ Pa < *P* < $10^5$ Pa region. This increased differential heating acts to drive a faster equatorial jet. In Figure 8 we show the absolute difference in the zonal-mean zonal wind between the relaxation and equilibrium simulations. There is clearly an

overall increase in the zonal wind velocity within the equatorial jet of 250–500 m s⁻¹. Comparing this with Figure 1, this corresponds to a 5%–10% increase compared with the equilibrium simulation.

The estimated dynamical and radiative timescales (Section 4.1) indicate that the atmosphere is radiatively driven for $P < 10^4$ Pa and dynamically driven for higher pressures. Advection of heat by the equatorial jet, aided by the increased wind velocities, leads to a heating between $10^4$ Pa < *P* < $10^5$ Pa all around the planet. For lower pressures, where the atmosphere is radiatively driven, the local temperature change is more spatially dependent.

We note that the longwave cooling is itself dependent on the atmospheric temperature and that the temperature changes between the equilibrium and relaxation simulations will result in different steady-state cooling rates. We further note that the net radiative heating rates shown in Figure 7 are approximately balanced by advection of heat in the latter stages of the simulation, where the atmosphere (for less than $P \sim 10^6$ Pa) is in an approximate steady state.

### 4.4. Summary

We find that for pressures less than $10^4$ Pa the atmosphere is generally cooler on the dayside and warmer on the nightside owing to wind-driven chemistry, compared with chemical equilibrium. For the pressure region $10^4$ Pa < *P* < $10^5$ Pa we find a significant warming all around the planet, in the equatorial region. An increased differential heating between the dayside and nightside drives a faster equatorial jet. These temperature changes lead to a decreased TOA radiative flux for the dayside and an increased TOA radiative flux for the nightside. Our testing shows that these changes are primarily due to changes in the abundance of methane.

## 5. Contribution Functions

To further understand the effect of wind-driven chemistry on the radiative properties of the atmosphere, we consider the contribution function, which quantifies the contribution of a layer to the upward intensity at the TOA. A peak in the contribution function effectively indicates the pressure level of the photosphere. The calculation of the contribution function is described and validated in Appendices A and B, respectively. We first discuss the structure of the contribution functions from the equilibrium simulation and compare with previous studies in Section 5.1. We then show the effect of wind-driven chemistry in Section 5.2.

### 5.1. Chemical Equilibrium

Figure 9 shows an area-weighted meridional mean (±20° latitude) of the contribution function as a function of longitude and pressure for the four *Spitzer*/IRAC channels centered on 3.6, 4.5, 5.8, and 8.0 µm. Methane is a prominent absorber in the 3.6 and 8.0 µm channels, while water and carbon monoxide are the primary absorbers in the 4.5 and 5.8 µm channels, respectively.

In the case of chemical equilibrium (left panels, Figure 9), it is immediately apparent that the contribution function varies significantly with longitude for the 3.6 and 8.0 µm channels. The peak of the contribution function in these channels generally occurs at higher pressures on the dayside compared
<figure><img src="image_7.png" /><figcaption>**Figure 6.** TOA thermal radiative flux for the equilibrium simulation (left) and the relative difference TOA flux between the relaxation and equilibrium simulations (right). A positive difference indicates a larger TOA flux in the relaxation simulation.</figcaption></figure>

with the nightside. On the other hand, the contribution function is relatively flat (isobaric) in the 4.5 and 5.8 µm channels.

The shapes of the 3.6 and 8.0 µm contribution functions are explained by the distribution of methane in chemical equilibrium. The dayside equilibrium abundance of methane is many orders of magnitude less than the nightside equilibrium abundance (Section 3), and therefore the opacity, within these spectral channels, is also smaller for the dayside. The smaller dayside opacity results in a smaller optical depth and hence a “deeper” photosphere. On the nightside, where opacity due to methane is significantly larger, the photosphere is shifted to lower pressures.

The contribution function in the 3.6 µm channel actually contains a “double peak” in the nightside region, with a secondary, slightly smaller peak at $P \sim 10^5$ Pa. This means that these two pressure regions both make a significant contribution to the intensity emerging from the TOA. We find that this double peak is due to the spectral dependence of the opacity within the 3.6 µm Spitzer/IRAC channel, which spans the range 3.08–4.01 µm and corresponds to ~75 spectral bands in the high spectral resolution setup of our model (see Section 2.2). For bands where methane has significant absorption the contribution function peaks at lower pressures, but for bands where methane absorption is not so prominent the contribution peaks at higher pressures. When combined in the 3.6 µm channel, this leads to two distinct peaks.

The 4.5 and 5.8 µm contribution functions are also consistent with the equilibrium composition. Both water and carbon monoxide show small variations with longitude, compared with methane, which varies by orders of magnitude. Therefore, the contribution functions in the spectral regions where these species are dominant absorbers are approximately isobaric around the equatorial region.

We now compare our results from the equilibrium simulation (Figure 9) with those of (Dobbs-Dixon & Cowan 2017, see their Figure 4), where we find several significant differences. We note that in our figures (Figure 9) the substellar point is located at 180° longitude, while in Figure 4 of Dobbs-Dixon & Cowan (2017) the substellar point is located at 0° longitude.

First, for the 3.6 and 8.0 µm channels, we find that the contribution function peaks at lower pressures on the nightside, compared with the dayside. Intuitively, we expect this to be the case given the larger abundance of methane on the nightside. However, the opposite appears to be true for the results of Dobbs-Dixon & Cowan (2017), with the contribution function

peaking at lower pressures on the dayside, compared with the nightside. Second, for the 4.5 and 5.8 µm channels Dobbs-Dixon & Cowan (2017) find that the contribution function (generally) shifts to a higher pressures on the nightside, in contrast to the approximately isobaric contribution functions that we find in the same channels. The temperature structure appears to be very similar between the two models (compare our Figure 1 with their Figure 4), and therefore the equilibrium chemical composition should also be similar. The cause of the discrepancy is therefore puzzling.

Whatever the cause of the difference, the opposing results have important consequences for the interpretation. Taking the 3.6 µm channel as an example, we find that the peak of the contribution function crosses many isotherms between the dayside and nightside of the planet, since the shape of the contribution function is approximately opposite to that of the isotherms. In contrast, the shapes of the contribution functions presented by Dobbs-Dixon & Cowan (2017) cause them to approximately follow the isotherms, meaning that the temperature contrast at the photosphere around the planet should be much smaller.

## 5.2. The Effect of Wind-driven Chemistry

The contribution functions for the relaxation simulations (right panels of Figure 9) show significant differences with the equilibrium simulations, particularly for the 3.6 and 8.0 µm channels. Most importantly, we find that the peaks of the contribution functions occur at lower pressures on the dayside than on the nightside, opposite to the chemical equilibrium case.

Since the chemical composition is approximately homogeneous in the relaxation simulation, the non-isobaric shape of the contribution functions is likely related to the temperature structure. To check, we performed a series of simple tests. First, we fixed the temperature at which the Planck function $B(T)$ (see Equation (8)) is evaluated, and the shapes of the contribution functions become flatter, though not isobaric. Second, we also fixed the temperature at which the absorption coefficient $\kappa$ is calculated. Fixing the temperature for both $B(T)$ and $\kappa$ yields an isobaric contribution function, indicating that it is the temperature dependence of these two terms that drives the shape of the contribution functions for the relaxation simulations.

Since the temperature at the peak of the contribution function (i.e., the photosphere) effectively determines the TOA flux,
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Left column: shortwave (top), longwave (middle), and net (bottom) heating rates (color scale and black contours) in $W m^{-3}$ for the equilibrium simulation. Right column: absolute differences in the heating rates (color scale and black contours) between the relaxation and equilibrium simulations for the shortwave (top), longwave (middle), and net (bottom). All panels show area-weighted meridional means between ±20° latitude. We note that, for the difference plots, a positive difference in the shortwave heating rate (top right panel) means a more positive heating rate (more heating), while a positive change in the longwave (middle right panel) means a less negative heating rate (less cooling). A positive difference in the net heating rate (bottom right panel) means a larger net heating rate (i.e., more net heating or less net cooling).</figcaption></figure>

changes in the shape of the contribution function will effect the observed emission from the atmosphere. As previously discussed, the shape of the contribution functions in the equilibrium case means that the photosphere cuts across many isotherms between the dayside and nightside. In contrast, for the relaxation case the contribution functions approximately follow the shape of the isotherms, leading to a much smaller temperature difference at the pressure level of the photosphere between the dayside and nightside.

With the above argument, this allows us to predict that we should find a larger phase amplitude in the simulated emission

phase curve, in the 3.6 µm channel, for the equilibrium simulation compared with the relaxation simulation.

### 5.3. Summary

We have presented contribution functions calculated for both the equilibrium and relaxation simulations of HD 189733b. In the chemical equilibrium case, we find relatively flat (isobaric) contribution functions for the 4.5 and 5.8 µm Spitzer/IRAC channels but non-isobaric contribution functions in the 3.6 and 8.0 µm channels, where methane is an important absorber. In
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Absolute difference in the zonal-mean temporal-mean (800–1000 days) zonal wind velocity between the relaxation and equilibrium simulations. A positive difference indicates an increased wind velocity in the relaxation simulation.</figcaption></figure>

the latter channels, the contribution function peaks at lower pressures on the nightside compared with the dayside, opposite to the trends found by Dobbs-Dixon & Cowan (2017). However, when including wind-driven chemistry, this trend reverses, with the contribution functions peaking at lower pressures on the dayside in the 3.6 µm channel. This has a significant effect on the temperature at the pressure level of the photosphere.

# 6. Emission Phase Curves

In this section we present the emission phase curves calculated from both the equilibrium and relaxation simulations of HD 189733b. We compare the results from the equilibrium simulation with results from previous studies and discuss the effect of wind-driven chemistry.

## 6.1. Chemical Equilibrium

Figure 10 shows the 3.6 and 4.5 µm emission phase curves from our simulations, as well as those from previous simulations of Showman et al. (2009) and Dobbs-Dixon & Agol (2013). We first compare our chemical equilibrium phase curves with those of previous models and focus on the normalized phase amplitude and the phase offset. The normalized phase amplitude quantifies the difference between the maximum and minimum flux ratio normalized by the flux ratio at secondary eclipse (e.g., Dobbs-Dixon & Cowan 2017), while the phase offset quantifies the angular shift between the maximum flux ratio in the phase curve and the secondary eclipse. The normalized phase amplitudes and phase offsets are shown separately in Figure 11, which is based on Figure 1 from Dobbs-Dixon & Cowan (2017).

First, for the 3.6 µm channel, we find that all three models show similar normalized phase amplitudes (see Figure 11). However, Figure 10 shows that the Dobbs-Dixon & Cowan (2017) phase curve has a significantly larger actual (i.e., non-normalized) phase amplitude, compared with both our result and that of Showman et al. (2009). There is a significant phase offset difference between all three models, with an offset of ~50° for Showman et al. (2009), ~35° for this work, and ~15° for Dobbs-Dixon & Cowan (2017).

Second, for the 4.5 µm channel, we obtain a similar normalized phase amplitude (~0.5) to that of Showman et al. (2009). The normalized phase amplitude of Dobbs-Dixon & Cowan (2017) is significantly larger (~0.8) than both our result and that of Showman et al. (2009), which is particularly apparent in Figure 10. The phase offsets for each model are similar to their values in the 3.6 µm channel.

The trends in the emission phase curves from our simulations, shown in Figures 10 and 11, are consistent with the trends in the contribution functions previously discussed in Section 4. The equilibrium 3.6 µm contribution function crosses many isotherms as it peaks at lower pressures on the nightside compared with the dayside, leading to a large normalized phase amplitude. In contrast, the 4.5 µm contribution function is approximately isobaric and crosses fewer isotherms, resulting in a smaller emission phase amplitude. However, the contribution functions in both channels peak at similar pressures on the dayside, near the location of the hot spot, giving a similar phase offset.

Figure 11 also shows the phase offset and normalized phase amplitude derived from observations of HD 189733b by Knutson et al. (2012). We note that none of the models reproduce the significant phase offset difference between the 3.6 and 4.5 µm channels.

## 6.2. The Effect of Wind-driven Chemistry

We now consider the effect of wind-driven chemistry on the emission phase curve, also shown in Figure 10. As expected from previous discussion of the contribution function, we find only a small difference in the 4.5 µm phase curve between the equilibrium and relaxation simulations. The phase offset is slightly increased and the normalized phase amplitude is slightly decreased compared with equilibrium. On the other hand, we find a significant difference in the 3.6 µm channel with a strongly decreased normalized phase amplitude compared with the equilibrium simulation.

This effect on the 3.6 µm channel phase curve is consistent with the previous discussion of the contribution functions (Section 5). Homogenization of the methane chemistry effectively inverts the shape of the contribution function (Figure 9) with the ultimate effect of significantly reducing the number of isotherms that the photosphere crosses, leading to a smaller emission phase amplitude. We note that the 3.6 µm normalized phase amplitude for the relaxation simulation provides a significantly poorer match to the observed normalized phase amplitude, compared with the equilibrium simulation.

## 6.3. Summary

The emission phase curves calculated from our equilibrium simulation indicate that we find a phase offset somewhere in between the previous results of Showman et al. (2009) and Dobbs-Dixon & Agol (2013). The normalized phase amplitude in the 3.6 µm channel is somewhat similar between all three models; however, Dobbs-Dixon & Agol (2013) obtain a significantly larger normalized phase amplitude in the 4.5 µm channel compared with both our results and those of Showman et al. (2009). The effect of wind-driven chemistry is to significantly decrease the normalized phase amplitude in the 3.6 µm channel, as well as to slightly increase the phase offset in both the 3.6 and 4.5 µm channels.
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Area-weighted meridional means (±20° latitude) of the normalized contribution function (color scale) for different Spitzer/IRAC channels: 3.6 µm (top row), 4.5 µm (second row), 5.8 µm (third row), and 8.0 µm (bottom row). The atmospheric temperature is also shown in black contours. Results from the equilibrium simulation are shown in the left column, and those from the relaxation simulation are shown in the right column.</figcaption></figure>
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Calculated emission phase curves from our simulations with both chemical equilibrium and chemical relaxation, as well as from Dobbs-Dixon & Agol (2013) and Showman et al. (2009) (both assuming chemical equilibrium), for the 3.6 µm (top) and 4.5 µm (bottom) *Spitzer*/IRAC channels.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Phase offset vs. the normalized phase amplitude for the 3.6 µm (blue) and 4.5 µm (green) *Spitzer*/IRAC channels from our UM simulations, as well as those of Dobbs-Dixon & Agol (2013) and Showman et al. (2009). Also shown are the phase offset and normalized phase amplitude derived from the observed phase curves of HD 189733b (Knutson et al. 2012).</figcaption></figure>

# 7. Discussion and Conclusions

## 7.1. Summary of Results

We have presented results from two simulations of the atmosphere of the hot Jupiter HD 189733b: one with the assumption of local chemical equilibrium and the other with the inclusion of wind-driven advection and chemical relaxation.

The trends in the chemical composition between the equilibrium and relaxation simulations are qualitatively similar to previous results for the warmer atmosphere of HD 209458b (Drummond et al. 2018b) and further demonstrate the importance of 3D modeling of exoplanet atmospheres. The main difference with HD 209458b is a larger equilibrium abundance of methane due to the lower temperature of the atmosphere, leading to larger methane mole fractions in the relaxation simulation. The net result of wind-driven advection is to increase the methane mole fraction by several orders of magnitude, compared with chemical equilibrium, throughout most of the modeled domain, with the largest effect on the dayside.

In Drummond et al. (2018b) we found relatively unimportant changes (~1%) to the thermal and dynamical structure due to wind-driven chemistry for the atmosphere of HD 209458b. However, for the simulations of HD 189733b we find more significant differences of up to ~10% in both temperature and wind velocities. For pressures less than $10^4$ Pa we estimate the radiative timescale to be faster than the dynamical timescale, and in this region we observe a local temperature decrease (increase) where the local methane abundance is increased (decreased), compared with chemical equilibrium. This is dominated by changes in the thermal (longwave) cooling.

For pressures greater than $10^4$ Pa, where advection of heat is expected to be more important, this trend no longer holds. We find a significant temperature increase between $10^4$ and $10^5$ Pa in the equatorial region due to an increase in the differential heating between the dayside and nightside. This increased differential heating drives a faster equatorial jet.

Our tests show that these changes to the temperature and circulation are due to the large increase in the methane abundance, compared with chemical equilibrium, with a minor contribution from the relatively small changes in the abundances of carbon monoxide and water.

The temperature changes found in this paper (up to 10%) using a 3D model are similar in magnitude to the temperature changes that we previously found using a 1D model (Drummond et al. 2016), including only vertical transport. One of the main effects of vertical transport in the 1D model is to increase the abundance of methane above chemical equilibrium.

When considering the TOA thermal radiative flux, it is clear that there is a total reduction in the amount of energy emitted from the dayside and an increase from the nightside when wind-driven chemistry is taken into account. Overall, this indicates an increased efficiency of heat transport from the dayside to nightside atmosphere.

We find significant quantitative and qualitative differences in the 3D contribution functions between our chemical equilibrium simulation and the results of Dobbs-Dixon & Cowan (2017). The contribution functions of Dobbs-Dixon & Cowan (2017) generally peak at lower pressures on the dayside compared with the nightside for all spectral channels considered. In contrast, our simulations show that the contribution functions in the 4.5 and
5.8 µm channels are approximately isobaric, while we find an inverse trend to Dobbs-Dixon & Cowan (2017) in the 3.6 and 8.0 µm channels with contribution functions that peak at lower pressures on the nightside. Since the temperature structure between our simulation and theirs appears to be similar, the cause of the discrepancy is undetermined. When including the effect of wind-driven chemistry, we find that the longitude dependence of the contribution functions in the 3.6 and 8.0 µm channels is essentially inverted, peaking at lower pressures on the dayside than on the nightside.

Differences in the shapes of the contribution functions explain the differences that we find in the simulated emission phase curves. The main effect of wind-driven chemistry is to significantly reduce the amplitude of the 3.6 µm phase curve, since the photosphere crosses fewer isotherms in the relaxation simulation compared with the equilibrium simulation. The 4.5 µm phase curve does not change significantly between the equilibrium and relaxation simulations.

Comparing the phase curves calculated from our equilibrium simulation with previous simulations of HD 189733b, we find fairly similar phase amplitudes to those of Showman et al. (2009) in both the 3.6 and 4.5 µm spectral channels, but larger differences in the phase amplitude from those of Dobbs-Dixon & Cowan (2017). The phase offset in both channels varies significantly between all three models. We find that none of the three models (this work; Showman et al. 2009; Dobbs-Dixon & Cowan 2017) match the observed phase curve characteristics well.

## 7.2. The Spectrum of HD 189733b and the Potential Effect of Clouds and Haze

The observed transmission spectrum of HD 189733b shows potential signatures of cloud or haze (Pont et al. 2008, 2013; Sing et al. 2016; Barstow et al. 2017). The presence of cloud or haze particles can affect the observed spectrum both directly and indirectly. The direct effect is to change the optical properties of the atmosphere and thereby change the observed transit or eclipse depth. The indirect effect occurs if the presence of cloud or haze particles alters the radiative heating rates and thereby also the thermal structure and atmospheric circulation.

The structure, composition, and effect of clouds on the background atmosphere have been investigated using a number of different 3D models with varying degrees of complexity and levels of approximation (Lee et al. 2016; Parmentier et al. 2016; Roman & Rauscher 2017; Lines et al. 2018b). Recently, Lines et al. (2018a) reported a flat transmission spectrum from their 3D simulations of HD 189733b that include a state-of-the-art treatment of cloud formation (Helling et al. 2008; Lee et al. 2016), contrary to the molecular absorption features that have been measured (Sing et al. 2016). This model–observation discrepancy suggests that current theoretical models are missing important physical mechanisms at play in the atmospheres of hot Jupiters (Lines et al. 2018a). Lee et al. (2017) investigated the effect of inhomogeneous cloud structures on the reflected light and thermal emission spectra of HD 189733b. Haze particles may also be produced in an HD 189733b-like atmosphere via photochemical processes, with their presence impacting the observed transmission spectrum and, potentially, the thermal profile of the atmosphere (Lavvas & Koskinen 2017).

The simulations presented in this paper assume that clouds and haze are not present in the atmosphere. The inclusion of clouds would likely alter the thermal and dynamical structure of the model atmosphere, and changes to the optical properties would also directly impact the emission phase curves that we present.

In this paper we choose to focus on the gas-phase chemistry and, in particular, how this interacts with the circulation and radiative heating. This is important because absorption from gas-phase chemical species (e.g., water) gives rise to the dominant features in many observed spectra of transiting exoplanets (e.g., Sing et al. 2016). In addition, models that include a sophisticated treatment of cloud formation currently assume a simple treatment of the gas-phase chemistry (i.e., local chemical equilibrium; Lee et al. 2016; Lines et al. 2018b). The simulations presented in this paper therefore represent an important step in the hierarchy of model complexity.

## 7.3. Future Prospects

In this study we found a larger effect of wind-driven chemistry on the thermal and dynamical structure compared with our previous results for the warmer case of HD 209458b (Drummond et al. 2018b), with differences in the temperature and wind velocities reaching ~10%. This is due to the cooler atmosphere of HD 189733b (compared with HD 209458b) having a larger equilibrium methane abundance and hence larger quenched methane abundance. A natural question to ask is whether this trend will continue with lower-temperature planets.

It is possible that a “sweet spot” exists where the dayside atmosphere is warm enough that carbon monoxide is dominant over methane, in chemical equilibrium, but the nightside is cool enough that methane is dominant over carbon monoxide. In this case, horizontal and/or vertical transport may lead to significant changes in the relative abundances of methane and carbon monoxide, with consequences for the opacity and heating rates. For cooler atmospheres still, it is likely that methane will be the dominant carbon species everywhere, with carbon monoxide present as a trace species. If this is true, we might expect a similar process to that presented in this paper, but with the roles of carbon monoxide and methane reversed.

Our simulations used a chemical relaxation scheme (Cooper & Showman 2006) to solve for the chemical evolution. The accuracy of this type of scheme is reliant on an accurate estimation of the chemical timescale (Tsai et al. 2018). While we have validated the Cooper & Showman (2006) chemical relaxation scheme against a full chemical kinetics calculation in a 1D model (Appendix C), a full kinetics calculation within a 3D framework has not yet been achieved in the literature. As we have previously discussed in Drummond et al. (2018b), there remain significant differences in the results between the 3D chemical relaxation approach (complex dynamics with simplified chemistry) and the pseudo-2D chemical kinetics approach (complex chemistry with simplified dynamics) of Agúndez et al. (2014). It is also unclear whether it is possible to include a treatment of photochemistry within the chemical relaxation method. For these reasons, a model consistently coupling chemical kinetics calculations within a 3D framework is required.

Using a 2D steady-state circulation model, Tremblin et al. (2017) found that the deep atmosphere ($P \gtrsim 1 \times 10^5$ Pa) of hot Jupiters may be significantly hotter than predicted by 1D radiative–convective equilibrium models, possibly due to the advection of potential temperature. General circulation models
also show a trend of converging toward a hotter profile in the deep atmosphere (Amundsen et al. 2016; Mayne et al. 2017), although, due to computational limitations, the models cannot be integrated for long enough to reach a steady state at these pressures.

A hotter deep atmosphere may directly impact atmospheric emission within the 4.5 µm Spitzer/IRAC channel, as the contribution function in this spectral channel peaks near this region in the simulations presented here. In addition, our results have shown that methane first departs from the equilibrium profile at pressures greater than $P = 1 \times 10^5$ Pa. Therefore, a hotter deep atmosphere may change the equilibrium abundance at the quench point and the pressure level of the quench point itself, with consequences for the chemical abundances at lower pressures.

Our results add to the evidence that solar composition gas-phase models fail to match the observed properties of the emission phase curve (i.e., the phase offset and normalized phase amplitude) within the Spitzer/IRAC 3.6 and 4.5 µm channels, both assuming local chemical equilibrium and with the effect of wind-driven chemistry. An important missing ingredient in the simulations presented in this paper is cloud formation (Section 7.2), which can have a significant impact on the thermal structure and observed spectra. Consistently coupling a cloud formation module with gas-phase wind-driven chemistry in a 3D model is a significant challenge, but this is likely required to understand the atmospheric composition of hot exoplanets. On the other hand, Nikolov et al. (2018) recently reported the clearest hot gas-giant atmosphere to date, which may prove to be a highly useful observational target for testing the theory of gas-phase chemistry.

The authors thank the anonymous referee for their report that helped to improve the quality of this paper. This work is partly supported by the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 grant agreement No. 336792-CREATES and No. 320478-TOFU). N.J.M. and J.G. are partially funded by a Leverhulme Trust Research Project Grant. J.M. acknowledges the support of a Met Office Academic Partnership secondment. This work was performed using the DiRAC Data Intensive service at Leicester, operated by the University of Leicester IT Services, which forms part 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 work also used the University of Exeter Supercomputer ISCA. Material produced using Met Office Software. The research data supporting this publication are openly available from the University of Exeter's institutional repository at: https://doi.org/10.24378/exe.884.

## Appendix A
### Deriving the Contribution Function

We start with the half-range (upward) intensity at the TOA (see Thomas & Stammes 1999, Section 5.4.3)

$$I_{\nu}^{+}(0, \mu, \phi) = I_{\nu}^{+}(\tau^{*}, \mu, \phi)e^{-\tau^{*}/\mu} + \int_{0}^{\tau^{*}} \frac{d\tau'}{\mu} e^{-\tau'/\mu} B_{\nu}(\tau'), \quad (3)$$

where $I_{\nu}^{+}$ is the hemispherical intensity at wavenumber $\nu$, $\mu = \cos\theta$, where $\theta$ is the zenith angle, $\phi$ is the azimuthal angle,

$B_{\nu}$ is the Planck function, and $\tau$ is the optical depth,

$$\tau = \int_{0}^{s} ds' \kappa(s'), \quad (4)$$

where $\kappa$ is the absorption coefficient and $s$ is the path.

The TOA is denoted $\tau = 0$, while the bottom of the atmosphere is denoted $\tau = \tau^{*}$. The first term on the right in Equation (3) is the upward intensity from the bottom of the atmosphere that has survived extinction, while the second term on the right is the thermal contribution from the atmosphere. For a gas-giant atmosphere with no solid surface we can assume $\tau^{*} = \infty$, which reduces Equation (3) to

$$I_{\nu}^{+}(0, \mu, \phi) = \int_{0}^{\infty} \frac{d\tau'}{\mu} e^{-\tau'/\mu} B_{\nu}(\tau'). \quad (5)$$

Changing the limits of the integral, we can write the intensity contribution of a particular layer (dropping the $\nu$ notation) as

$$I^{+}(0, \mu, \phi) \equiv \mathcal{CF}_{I}(\mu, \phi) = \int_{\tau_{1}}^{\tau_{2}} \frac{d\tau'}{\mu} e^{-\tau'/\mu} B(\tau'), \quad (6)$$

where we have defined the intensity contribution function $\mathcal{CF}_I$. Evaluating the integral, we obtain

$$\begin{aligned} \mathcal{CF}_I(\mu, \phi) &= B(\tau_{1,2})[e^{-\tau_1/\mu} - e^{-\tau_2/\mu}] \\ \mathcal{CF}_I(\mu, \phi) &= B(\tau_{1,2})d[e^{-\tau/\mu}], \end{aligned} \quad (7)$$

where $B(\tau_{1,2})$ is the Planck function of the layer $\tau_1 \rightarrow \tau_2$, which is assumed to be constant across the layer.

Discretizing Equation (7) onto a model grid would give the contribution of a model layer to the emergent intensity. To generalize this, we can introduce a factor $1/d \log P$ (e.g., Chamberlain & Hunten 1987) to give the contribution per decade in pressure. Finally, then, we have

$$\mathcal{CF}_I(\mu, \phi) = B(\tau_{1,2}) \frac{d[e^{-\tau/\mu}]}{d[\log P]} \quad (8)$$

It is clear that Equation (8) is equivalent to Equation (2) of Knutson et al. (2009) if $\mu = 1$. For the calculations presented in this paper $\mu = 1$. However, our equation is different from that presented in both Griffith et al. (1998) and Dobbs-Dixon & Cowan (2017), where the negative sign is missing in the exponential factor; we assume that this negative sign is taken care of elsewhere in their implementation.

Importantly, we note that $\mathcal{CF}_I$ is a physical quantity with units W m$^{-2}$ sr$^{-1}$. However, it is more common to present the normalized contribution function ($\mathcal{CF}_I$) such that the layer that contributes the most to the TOA intensity has a value $\mathcal{CF}_I = 1$. In the main body of this paper we present the quantity $\mathcal{CF}_I$.

## Appendix B
### Validating the Contribution Function Calculation

To validate the implementation of the contribution function calculation, we first derive the flux contribution function ($\mathcal{CF}_F$). Integrating Equation (7) over $\mu$ and $\phi$
<figure><img src="image_14.png" /><figcaption>Figure 12. Flux contribution function $CF_F$ and the net flux divergence $\nabla F$ across the layer for a 1000 K isothermal atmosphere.</figcaption></figure>

(see Thomas & Stammes 1999, Equation (5.48)),

$$ 
\begin{aligned}
CF_F &= \int_0^{2\pi} d\phi \int_0^1 d\mu \mu CF_I(\mu, \phi) \\
CF_F &= 2\pi B(\tau_{1,2}) \int_0^1 d\mu [\mu e^{-\tau/\mu} - \mu e^{-\tau_2/\mu}].
\end{aligned} 
\quad (9) $$

To evaluate the integral in Equation (9), we use the diffusivity approximation (see Thomas & Stammes 1999, Section 11.2.5)

$$ \int_0^1 d\mu \mu e^{-\tau/\mu} = \frac{e^{-D\tau}}{D}, $$

where D is the diffusivity factor. Finally, we now have

$$ 
\begin{aligned}
CF_F &= \frac{2\pi B(\tau_{1,2})}{D} [e^{-D\tau_1} - e^{-D\tau_2}] \\
CF_F &= \frac{2\pi B(\tau_{1,2})}{D} d[e^{-D\tau}].
\end{aligned} 
\quad (10) $$

$CF_F$ is the amount of flux that escapes to space from a particular layer of the atmosphere. For the special case of an isothermal atmosphere this will be equal to the divergence of the radiative flux ($\nabla F$) in the layer, as exchange with other layers of the atmosphere will be zero (see Thomas & Stammes 1999, Section 11.2.7). We note that $\nabla F$ is the main output of the radiative scheme and is extensively tested. In the isothermal case $D \sim 2$ (Edwards 1996).

We compare $CF_F$ and $\nabla F$ for a 1000 K isothermal atmosphere (where we have used $D = 2.0$) in Figure 12, and we find excellent agreement between the two quantities, validating our implementation of the contribution function calculation.

## Appendix C

### Validating the Chemical Relaxation Scheme

Chemical relaxation methods rely on the accurate estimation and parameterization of a chemical timescale that is based on results from a full chemical kinetics calculation (Tsai et al. 2018). To assess the accuracy of the chemical relaxation scheme that we use in this study (Cooper & Showman 2006),

<figure><img src="image_15.png" /><figcaption>Figure 13. Mole fraction of methane along a 1D profile with different values of the $K_{zz}$. We compare the steady-state abundances obtained from the Cooper & Showman (2006) chemical relaxation scheme (CS06 Relax) with the Tsai et al. (2018) chemical relaxation scheme (T17 Relax), as well as chemical kinetics calculations using both the Venot et al. (2012) and Tsai et al. (2017) chemical networks (V12 Kinetics and T17 Kinetics, respectively). The chemical equilibrium mole fraction is shown in dotted blue. The calculations are performed with the pressure-temperature profile and planetary parameters (e.g., surface gravity) of HD 189733b.</figcaption></figure>

we compare it against a full chemical kinetics calculation using two different chemical networks: Venot et al. (2012) and Tsai et al. (2017). We used the chemical equilibrium pressure-temperature profile and model parameters for HD 189733b from Drummond et al. (2016) and tested a range of eddy diffusion coefficients ($10^8 < K_{zz} < 10^{11} \text{ cm}^2 \text{ s}^{-1}$) that are constant with pressure.

Figure 13 shows the comparison of the Cooper & Showman (2006) chemical relaxation with the Venot et al. (2012) and Tsai et al. (2017) chemical kinetics schemes, as well as the Tsai et al. (2018) relaxation scheme. For the smallest $K_{zz}$ value, the Cooper & Showman (2006) relaxation result matches very well with the Tsai et al. (2017) kinetics result, while for larger values the Cooper & Showman (2006) relaxation result is somewhere in between the Tsai et al. (2017) and Venot et al. (2012) kinetics results.

We conclude that the Cooper & Showman (2006) chemical relaxation method has an acceptable level of accuracy, compared with a full chemical kinetics calculation, for the current application. Importantly, we note that the result obtained with the Cooper & Showman (2006) relaxation scheme typically lies somewhere in between the results obtained with the Tsai et al. (2017) and Venot et al. (2012) kinetics schemes, depending on the value of $K_{zz}$. Therefore, potential inaccuracies of the present chemical relaxation method are within the expected uncertainties of currently available chemical networks, for HD 189733b conditions.

## Appendix D

### Sensitivity to the Chemical Timescale

To test the sensitivity of the model to the value of the chemical timescale, we performed additional simulations that are identical to the relaxation simulation but with an artificially increased/decreased chemical timescale by a factor of 10. The
<figure><img src="image_16.png" /><figcaption>Figure 14. Vertical profiles of the chemical equilibrium (dashed) and chemical relaxation (solid) mole fractions of methane, carbon monoxide, and water for several longitude points around the equator, for simulations with the Cooper & Showman (2006) chemical timescale multiplied by a factor of 10 (left panel) and the Cooper & Showman (2006) chemical timescale divided by a factor of 10 (right panel). This figure should be compared with Figure 4, which shows the simulation with the nominal Cooper & Showman (2006) chemical timescale.</figcaption></figure>

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

vertical mole fraction profiles around the equator for these tests are shown in Figure 14 and should be compared with the results using the nominal chemical timescale in Figure 4.

We find that an increase (decrease) in the chemical timescale leads to an increase (decrease) in the pressure level of the “turnoff” point, where the methane abundance begins the increase with decreasing pressure. The location of the quench point depends on the ratio of the transport to chemical timescales. The chemical timescale typically increases with decreasing pressure owing to the pressure and temperature dependence of the reaction rate, while the transport timescale depends on the local wind velocities and the relevant length scale (see Drummond et al. 2018b, Figure 3). An increase in the chemical timescale will therefore shift the quench point to higher pressures, while decreasing the chemical timescale will shift the quench point to lower pressures.

While changing the chemical timescale clearly alters the pressure level at which meridional transport becomes important, the final vertically quenched methane abundance shows only a small variation between the three cases. The quenched mole fraction of methane using the nominal Cooper & Showman (2006) timescale is ~4.5 × 10⁻⁵, while in the ×10 and ÷10 cases this changes to ~6 × 10⁻⁵ and ~3.5 × 10⁻⁵, respectively. This test suggests that uncertainties or errors in the value of the chosen timescale do not lead to large differences in our results.

## Appendix E
Conservation of Global Mass of the Elements

As a basic test of the advection scheme, we check the conservation of the global mass of carbon and oxygen. We note that this test does not necessarily validate the accuracy of the advection scheme but does provide a basic check that the advection scheme is not artificially gaining or losing mass.

The total mass of the atmosphere M can written as

$$M = \sum_k \rho_k V_k, \qquad (11)$$

where ρ_k and V_k are the mass density and volume of the cell k, respectively, and the sum is over the total number of cells in the model grid. Similarly, the total mass of an element i can be

written as

$$M_i = \sum_k \sum_j \alpha_{i,j} w_{j,k} \rho_k V_k = \frac{1}{\mu} \sum_k \sum_j \alpha_{i,j} \mu_j f_{j,k} \rho_k V_k, \quad (12)$$

where w_{j,k} and f_{j,k} are the mass fraction and mole fraction of the species j, respectively, α_{i,j} is the fractional mass of element i in species j, μ_j is the molar mass of the species j, and μ is the mean molar mass. Specifically for the simulations presented in this paper, which include methane, water, and carbon monoxide, we can write the global mass of carbon as

$$\begin{aligned} M_C &= \frac{1}{\mu} \sum_k \rho_k V_k [\alpha_{C,CO} \mu_{CO} f_{CO,k} + \alpha_{C,CH_4} \mu_{CH_4} f_{CH_4,k}], \\ M_C &= \frac{1}{\mu} \sum_k \rho_k V_k [\alpha_{C,CO} \mu_{CO} f_{CO,k} \\ & \qquad + \alpha_{C,CH_4} \mu_{CH_4} (A_C - f_{CO,k})], \end{aligned} \qquad (13)$$

where we have replaced f_{CH_4,k} with f_{CO,k} using the mass balance approach and A_C is a constant (Cooper & Showman 2006; Drummond et al. 2018b). Similarly for oxygen, we have

$$\begin{aligned} M_O &= \frac{1}{\mu} \sum_k \rho_k V_k [\alpha_{O,CO} \mu_{CO} f_{CO,k} + \alpha_{O,H_2O} \mu_{H_2O} f_{H_2O,k}], \\ M_O &= \frac{1}{\mu} \sum_k \rho_k V_k [\alpha_{O,CO} \mu_{CO} f_{CO,k} \\ & \qquad + \alpha_{O,H_2O} \mu_{H_2O} (A_O - f_{CO,k})]. \end{aligned} \qquad (14)$$

As before, we have replaced f_{H_2O,k} with f_{CO,k} by assuming mass balance and that A_O is a constant. We note that it is the mole fraction of carbon monoxide (f_{CO}) that is advected as a tracer in these simulations.

Figure 15 shows the conservation of the global mass of carbon and oxygen as a function of simulation time. The global mass of carbon and oxygen is conserved to within significantly better than 99.9% over 1000 days. As a point of comparison,
<figure><img src="image_18.png" /><figcaption>**Figure 15.** Left panel: percentage of the initial global mass of carbon and oxygen throughout the simulation. Right panel: global mass of carbon monoxide throughout the simulation.</figcaption></figure>

Figure 15 also shows the change in the global mass of carbon monoxide, which increases by approximately 25%.

The fractional change in the mass of carbon monoxide is much larger than the errors in the conservation of the global mass of elemental carbon and oxygen. This suggests that the change in the mass is due to the physics schemes included in the model and not due to numerical gain or loss of mass from the advection scheme.

## ORCID iDs

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

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

Pascal Tremblin https://orcid.org/0000-0001-6172-3403

David K. Sing https://orcid.org/0000-0001-6050-7645

## References

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

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

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

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

Asplund, M., Grevesse, N., Sauval, A. J., & Scott, P. 2009, ARA&A, 47, 481

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

Barstow, J. K., Aigrain, S., Irwin, P. G. J., & Sing, D. K. 2017, ApJ, 834, 50

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

Chamberlain, J. W., & Hunten, D. M. 1987, Theory of Planetary Atmospheres. An Introduction to their Physics and Chemistry (Orlando, FL: Academic Press)

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

Dang, L., Cowan, N. B., Schwartz, J. C., et al. 2018, NatAs, 2, 220

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

Dobbs-Dixon, I., & Cowan, N. B. 2017, ApJL, 851, L26

Dobbs-Dixon, I., & Lin, D. N. C. 2008, ApJ, 673, 513

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

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

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

Edwards, J. M. 1996, JATS, 53, 1921

Edwards, J. M., & Slingo, A. 1996, QJRMS, 122, 689

Fortney, J. J., Marley, M. S., Lodders, K., Saumon, D., & Freedman, R. 2005, ApJL, 627, L69

Goyal, J. M., Mayne, N., Sing, D. K., et al. 2018, MNRAS, 474, 5158

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

Helling, C., Woitke, P., & Thi, W.-F. 2008, A&A, 485, 547

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

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

Irwin, P. G. J., Teanby, N. A., de Kok, R., et al. 2008, JQSRT, 109, 1136

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

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

Knutson, H. A., Lewis, N., Fortney, J. J., et al. 2012, ApJ, 754, 22

Lavvas, P., & Koskinen, T. 2017, ApJ, 847, 32

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

Lee, G. K. H., Wood, K., Dobbs-Dixon, I., Rice, A., & Helling, C. 2017, A&A, 601, A22

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

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

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

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

Louden, T., & Wheatley, P. J. 2015, ApJL, 814, L24

Madhusudhan, N., Burrows, A., & Currie, T. 2011, ApJ, 737, 34

Madhusudhan, N., & Seager, S. 2009, ApJ, 707, 24

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

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

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

Menou, K., & Rauscher, E. 2009, ApJ, 700, 887

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

Nikolov, N., Sing, D. K., Fortney, J. J., et al. 2018, Natur, 557, 526

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

Pont, F., Knutson, H., Gilliland, R. L., Moutou, C., & Charbonneau, D. 2008, MNRAS, 385, 109

Pont, F., Sing, D. K., Gibson, N. P., et al. 2013, MNRAS, 432, 2917

Rauscher, E., & Menou, K. 2013, ApJ, 764, 103

Roman, M., & Rauscher, E. 2017, ApJ, 850, 17

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

Showman, A. P., & Guillot, T. 2002, A&A, 385, 166

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

Snellen, I. A. G., de Kok, R. J., de Mooij, E. J. W., & Albrecht, S. 2010, Natur, 465, 1049

Southworth, J. 2010, MNRAS, 408, 1689

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

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

Tsai, S.-M., Kitzmann, D., Lyons, J. R., et al. 2018, ApJ, 862, 31

Tsai, S.-M., Lyons, J. R., Grosheintz, L., et al. 2017, ApJS, 228, 20

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

Waldmann, I. P., Tinetti, G., Rocchetto, M., et al. 2015, ApJ, 802, 107

Wood, N., Staniforth, A., White, A., et al. 2014, QJRMS, 140, 1505

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