A reduced chemical scheme for modelling warm to hot hydrogen-dominated atmospheres

TL;DR

This paper introduces a simplified chemical scheme with 30 species for 3D atmospheric modeling of warm to hot hydrogen-rich planets and brown dwarfs, enabling faster computations without significant accuracy loss.

Contribution

A new reduced chemical scheme for 3D atmospheric modeling that balances computational efficiency and accuracy for a wide range of warm planetary atmospheres.

Findings

Accurately reproduces vertical profiles of key observable species.

Valid for a broad temperature range, excluding hot C-rich atmospheres.

Applicable in both 1D and 3D models for exoplanets and brown dwarfs.

Abstract

Three dimensional models that account for chemistry are useful tools to predict the chemical composition of (exo)planet and brown dwarf atmospheres and interpret observations of future telescopes. Recent Juno observations of the NH3 tropospheric distribution in Jupiter also indicate that 3D chemical modelling may be necessary to constrain the deep composition of the giant planets of the Solar System. However, due to the high computational cost of chemistry calculations, 3D chemical modelling has so far been limited. A solution to include chemical kinetics in a 3D model is to use a reduced chemical scheme. In this view, we have developed a reduced scheme containing 30 species and 181 reversible reactions, from the full chemical scheme of Venot et al. 2012. This scheme reproduces accurately the vertical profiles of the observable species (H2O, CH4, CO, CO2, NH3, and HCN) and has a large range of validity. It can be used to study all kind of warm atmospheres, except hot C-rich ones, which contains a high amount of C2H2. It can be used in 1D models, for fast computations, but also in 3D models for hot giant (exo)planet and brown dwarf atmospheres.

Tables

Table 1. Table 1: For the nominal model, maximum ( Δ ​ m ​ a ​ x Δ 𝑚 𝑎 𝑥 \Delta max ) differences of abundances (in %) for each species for which the reduced scheme has been built. The pressure level (@level in mbar) where the maximum variation is reached is indicated within parentheses. These values are calculated within the [0.1–1000] mbar region only, which is probed by infrared observations.

Species	$\Delta max$
H2O	2$\times$10-4 (@6$\times$102)
CH4	7$\times$10-3 (@8$\times$102)
CO	1$\times$10-2 (@1$\times$10-1)
CO2	2$\times$10-2 (@3$\times$102)
NH3	2$\times$10-3 (@1$\times$10-1)
HCN	2$\times$10-3 (@8$\times$102)

Table 2. Table 2: For GJ 436b models, with various elemental abundances, maximum differences of abundances (in %) for each species for which the reduced scheme is designed. The pressure level (@level in mbar) where the maximum variation is reached is indicated in parenthesis. These values are calculated within the [0.1–1000] mbar region only, which is probed by infrared observations.

Species	Metallicity = 10	Metallicity = 100	C/O =0.5$\odot$	C/O = 3$\odot$
H2O	2$\times$10-2 (@1$\times$10-1)	2$\times$10-1 (@7)	4$\times$10-4 (@1$\times$103)	3$\times$10-3 (@5$\times$10-1)
CH4	6$\times$10-2 (@9$\times$102)	6$\times$10-1 (@6$\times$102)	2$\times$10-2 (@1$\times$103)	6$\times$10-3 (@7)
CO	1$\times$10-1 (@5$\times$102)	2$\times$10-1 (@5$\times$102)	8$\times$10-3 (@1$\times$103)	7$\times$10-2 (@6$\times$102)
CO2	2$\times$10-1 (@1$\times$10-1)	1$\times$10-1 (@7)	2$\times$10-2 (@2$\times$101)	9$\times$10-2 (@1$\times$102)
NH3	1$\times$10-3 (@1$\times$10-1)	2$\times$10-2 (@6$\times$102)	2$\times$10-3 (@1$\times$10-1)	1$\times$10-3 (@8$\times$102)
HCN	1$\times$10-1 (@1$\times$10-1)	3$\times$10-1 (@4$\times$102)	1$\times$10-2 (@1$\times$103)	3$\times$10-2 (@1$\times$10-1)

Table 3. Table 3: For HD 209458b, HD 189733b, GJ 1214b models, maximum variations of abundances (in %) for each species for which the reduced scheme is designed. The pressure level (@level in mbar) where the maximum variation is reached is indicated in parenthesis. These values are calculated within the regions probed by infrared observations ([0.1–1000] mbar).

Species	HD 209458b	HD 189733b	GJ 1214b
H2O	1$\times$10-1 (@1$\times$10-1)	6$\times$10-1 (@3$\times$10-1)	3$\times$10-1 (@3$\times$101)
CH4	6$\times$102 (@1$\times$10-1)	1$\times$10-1 (@3$\times$102)	6$\times$10-1 (@3$\times$102)
CO	7$\times$10-2 (@1$\times$10-1)	5$\times$10-1 (@6$\times$101)	1 (@2$\times$102)
CO2	1$\times$10-2 (@3$\times$10-1)	8$\times$10-1 (@1$\times$10-1)	9$\times$10-1 (@1$\times$10-1)
NH3	2$\times$103 (@1$\times$10-1)	2$\times$10-2 (@1$\times$10-1)	3$\times$10-2 (@1$\times$10-1)
HCN	6$\times$101 (@1$\times$10-1)	1 (@1$\times$102)	1 (@2$\times$102)

Table 4. Table 4: For Uranus, Neptune, and SD 1110 models, maximum variations of abundances (in %) for each species for which the reduced scheme is designed. The pressure level (@level in mbar) where the maximum variation is reached is indicated in parenthesis. These values are calculated within the region where quenching occurs for the Giant Planets ([10 6 –10 7 ] mbar) and within the region probed by infrared observations ([10 2 –10 4 ] mbar) for the Brown Dwarf.

Species	Uranus	Neptune	SD 1110
H2O	6$\times$10-1 (@8$\times$106)	2 (@9$\times$106)	2 (@9$\times$103)
CH4	1$\times$10-1 (@1$\times$106)	1 (@9$\times$106)	3 (@9$\times$103)
CO	3 (@1$\times$106)	6 (@1$\times$106)	2$\times$10-1 (@9$\times$103)
CO2	3 (@1$\times$106)	1$\times$101 (@1$\times$106)	3 (@1$\times$102)
NH3	7$\times$10-1 (@8$\times$106)	6 (@1$\times$106)	9$\times$10-2 (@1$\times$103)
HCN	4$\times$10-1 (@5$\times$106)	5 (@1$\times$106)	3 (@2$\times$103)

Table 5. Table 5: For HD 209458b with a C/O = 3 × \times and 6 × \times solar value, maximum variations of abundances (in %) for each species for which the reduced scheme is designed. The pressure level (@level in mbar) where the maximum variation is reached is indicated in parenthesis. These values are calculated within the regions probed by infrared observations ([0.1–1000] mbar).

Species	C/O =3$\odot$	C/O=6$\odot$
H2O	1$\times$102 (@1$\times$10-1)	9$\times$101 (@2$\times$10-1)
CH4	4$\times$105 (@1$\times$10-1)	8$\times$105 (@1$\times$10-1)
CO	1$\times$10-1 (@1$\times$10-1)	2$\times$10-1 (@1$\times$10-1)
CO2	1$\times$102 (@1$\times$10-1)	9$\times$101 (@2$\times$10-1)
NH3	1$\times$104 (@2$\times$10-1)	1$\times$104 (@2$\times$10-1)
HCN	3$\times$102 (@1$\times$10-1)	1$\times$102 (@1$\times$10-1)

Full text

*Aims. *Our goal is to develop a reduced chemical scheme from the full chemical scheme of Venot et al. (2012) able to reproduce accurately the vertical profiles of the observable species (H2O, CH4, CO, CO2, NH3, and HCN). This reduced scheme should have a size compatible with three dimensional models and be usable across a large parameter space (e.g. temperature, pressure, elemental abundance). The absence of C2H2 from our reduced chemical scheme prevents its use to study hot C-rich atmospheres.

*Methods. *We use a mechanism-processing utility designed for use with Chemkin-Pro to reduce a full detailed mechanism. ANSYS Chemkin-Pro Reaction Workbench allows the reduction of a reaction mechanism for a given list of target species and a specified level of accuracy. We take a warm giant exoplanet with solar abundances, GJ 436b, as a template to perform the scheme reduction. To assess the validity of our reduced scheme, we take the uncertainties on the reaction rates into account in Monte-Carlo runs with the full scheme, and compare the resulting vertical profiles with the reduced scheme. We explore the range of validity of the reduced scheme even further by applying our new reduced scheme to GJ 436b’s atmosphere with different elemental abundances, to three other exoplanet atmospheres (GJ 1214b, HD 209458b, HD 189733b), a brown dwarf atmosphere (SD 1110), and to the troposphere of two giant planets of the Solar System (Uranus and Neptune).

*Results. *For all cases except one, the abundances predicted by the reduced scheme remain within the error bars of the model with the full scheme. Expectedly, we found important differences that cannot be neglected only for the C-rich hot atmosphere. The reduced chemical scheme allows more rapid runs than the full scheme it derived from ($\sim$30 times faster).

*Conclusions. *We have developed a reduced scheme containing 30 species and 181 reversible reactions. This scheme has a large range of validity and can be used to study all kind of warm atmospheres, except hot C-rich ones, which contains a high amount of C2H2. It can be used in 1D models, for fast computations, but also in 3D models for hot giant (exo)planet and brown dwarf atmospheres.

1 Introduction

The next generation of space telescopes, such as the James Webb Space Telescope (JWST) and the Atmospheric Remote-Sensing Infrared Exoplanet Large-survey (ARIEL) will certainly revolutionise our knowledge and understanding of exoplanetary worlds. The high sensitivity of these future telescopes will permit the determination of horizontal and vertical variations of temperature and chemical composition in the atmospheres of hot and warm giant planets (e.g. Bean et al. 2018; Tinetti et al. 2018; Venot et al. 2018 and Venot et al. in prep.). To fully benefit from these high-resolution data and accurately interpret those spectra, three-dimensional (3D) models that include a detailed description of both physical and chemical processes are needed. Such models have not been developed yet, because of the very high computational time they require. Indeed, full chemical networks used to model Hot Jupiters atmospheres (e.g. Moses et al. 2011; Venot et al. 2012) are very large, e.g. $\sim$2000 reactions including a hundred species. Resolving the corresponding large system of non-linear and highly coupled differential equations, in addition to the physical processes described in General Circulation Models (GCMs) would result in unreasonable computational times. To tackle this problem, most GCMs currently assume atmosphere at thermochemical equilibrium. Even if correct for very hot environments, this assumption might be wrong in cooler atmospheres and/or atmospheres subject to very strong transport.

Two attempts have been undertaken so far to couple dynamics and disequilibrium chemistry. First, Cooper & Showman (2006) have studied carbon chemistry in the atmosphere of HD 209458b thanks to a three-dimensional model. In this model, the repartition of carbon into carbon monoxide (CO) is determined by relaxing the mole fraction of CO to its equilibrium value and that into methane (CH4) is obtained through mass balance. This simple approach has shown that, at all latitudes, the mixing ratios of CO and CH4 are homogenised with longitude in the 1–1000 mbar pressure region. The authors argue that vertical mixing is more efficient than horizontal mixing in the atmosphere of HD 209458b. This method was later used by other groups. The same chemical relaxation was applied to the atmospheres of HD 209459b (Drummond et al. 2018a) and HD 189733b (Drummond et al. 2018b) using a different GCM, with the chemistry consistently coupled with the radiative transfer. They found also vertical mixing to be the dominant disequilibrium process over a large pressure range, albeit horizontal mixing has a non-negligible effect at low pressures and especially for CH4. Also, Mendonça et al. (2018) used another relaxation scheme (Tsai et al. 2018) to study the disequilibrium chemistry of CO, CO2, H2O, and CH4 in WASP-43b atmosphere. Contrary to the previous studies of Cooper & Showman (2006) and Drummond et al. (2018a, b), they found that horizontal quenching is the dominant process. This dissimilarity was attributed to differences in dynamical timescales, due to different model parameters (clear/cloudy atmosphere, surface gravity, rotation rate). As shown by Mayne et al. (2017), variations in numerical settings, such as numerical drag or boundaries conditions, could also be responsible for these different results.

An alternative method to study disequilibrium chemistry in Hot Jupiters has been undertaken by Agúndez et al. (2012, 2014). They have developed pseudo-2D photo-thermochemical models of HD 209458b and HD 189733b with full chemical kinetics. With this method, the 1D atmospheric column rotates along the equator to mimic horizontal mixing. The chemical composition is determined at steady-state by integrating a full chemical kinetics scheme (Venot et al. 2012). As Cooper & Showman (2006), they found an homogenisation of the chemical composition with longitude, but their conclusion on the relative importance of horizontal and vertical quenching phenomenon differs. Agúndez et al. (2012, 2014) have determined that vertical mixing and horizontal mixing are both at play in Hot Jupiter atmospheres. These various results obtained with different methods show the complexity of atmospheric modelling and the difficulty to draw a general behaviour. Finally, the best way to determine 3D maps of the atmospheric composition of exoplanets is to include a realistic and detailed description of chemical kinetics in 3D dynamical models.

While hot giant exoplanet atmospheric studies would undoubtedly benefit from coupled chemical and dynamical modelling in 3D, it seems to be more and more the case for Solar System giant planet interiors as well. Until now, only 1D models have been applied to the deep tropospheres of these planets to try to unveil their primordial composition (e.g. Visscher et al. 2010; Mousis et al. 2014; Cavalié et al. 2014, 2017). However, recent Juno/MWR observations of the meridional distribution of NH3 in Jupiter’s troposphere (Bolton et al. 2017) surprisingly show that NH3 is not well-mixed below its condensation level. This may then also be the case for other minor species, and 3D chemical and dynamical modelling of Jupiter’s troposphere may be required when attempting to use upper tropospheric abundance measurements to derive deep elemental composition. This is probably also valid for the other Solar System giant planets, in which the temporal and meridional variability in the appearance of convective storms and clouds (e.g. Fischer et al. 2011; de Pater et al. 2015; Irwin et al. 2016) may be the observable outcome of deep tropospheric inhomogeneities.

In this view, we have developed a reduced chemical scheme from the full chemical network published by Venot et al. (2012), called C0-C2 scheme. We have chosen as a template the atmosphere of the warm Neptune GJ 436b, where thermochemistry and disequilibrium chemistry are both at play, and with solar abundances. This reduced scheme has a size comparable to the ones already included in GCMs developed for terrestrial planets/satellites of the Solar System, that is to say 40 species and 284 reactions for Titan (Lebonnois et al. 2001), 12 species and 42 reactions for Mars (Lefèvre et al. 2004), 34 species and 121 reactions for Venus (Stolzenbach et al. 2014; Stolzenbach 2016; Gilli et al. 2017). Therefore, we presume that our reduced chemical scheme, containing 30 species en 181 reversible reactions (i.e. 362 reactions in total), should be implementable in 3D models. We explain in Sect.2 the method we have used to develop the chemical scheme and to determine its level of accuracy. Then, we present in Sect. 3 the reduced chemical scheme we have obtained and compare the atmospheric composition of GJ 436b obtained with this scheme with that obtained with the full original scheme. In Sect. 4, we study the range of validity of this scheme, i.e. the range of atmospheric conditions for which the chemical composition obtained with the reduced scheme is similar (or not) to the one obtained with the full chemical scheme for a set of species of interest, i.e. the observable ones. We present our conclusions in Sect. 5.

2 Method

2.1 Reduction of the chemical scheme

We have used the kinetic model of Venot et al. (2012) to perform the chemical scheme reduction presented in this paper. The initial full chemical scheme contains only neutral chemical reactions and we have excluded photodissociations. There are several methods that exist to reduce a chemical scheme, like the one proposed by Dobrijevic et al. (2011) for Solar System Giant Planet stratospheres. However such methods are not adapted to cases with reversible reactions. To reduce our chemical scheme, we have used the package ANSYS Chemkin-Pro Reaction Workbench (2017).

We have used the method DRGEP (Directed Relation Graph with Error Propagation), followed by a Sensitivity Analysis (SA). Initially, the Directed Relation Graph (DRG) method has been developed by Lu & Law (2005) and extended to the DRGEP method by Pepiot-Desjardins & Pitsch (2008) followed by Liang et al. (2009). The mathematical explanation of these methods is beyond the scope of this paper, but the detailed algorithms of these mechanism reduction methods were described in Lebedev et al. (2013). To summarize, from a list of species of interest (target species) defined by the user, in the DRGEP method, the significance of all the non-target species for description of the target species is determined by analysing the rates of transformations between them. Each non-target species is thus characterised by an importance index regarding to each target species. The maximum of these values (amount equals to the number of target species) define the overall importance of each non-target species in the chemical scheme. According to the level of accuracy/completeness desired by the user, all species with an overall importance index above the threshold value are included in the reduced scheme. In the SA method, the maximum value of the concentration sensitivity to each reaction is calculated and compared to a threshold value. The reactions with all sensitivity coefficients smaller than the threshold value are excluded from the mechanism.

We can mention that these approaches are commonly used in the field of combustion. Recently, Qiu et al. (2016) have used these approaches to transform an n-decane/$\alpha$-methylnaphthalene/polycyclic aromatic hydrocarbon (PAH) kinetic mechanism involving 108 species and 846 reactions into a reduced model including 56 species and 236 reactions. The newly developed mechanism was validated by experimental data in fundamental reactors, including mole fraction of key species in the ethylene premixed flame and jet stirred reactor, as well as ignition delay time of pure and mixed fuel. In our study, these methods identify and select the key reactions that govern the chemical composition for a given thermal profile and given elemental abundances. We have first employed the DRGEP methodology to identify and eliminate unimportant species and associated reactions. We then implemented a reaction sensitivity analysis module to eliminate less important reactions. The use of DRGEP before the sensitivity analysis is required owing to the fact that SA is very time-consuming.

We have applied these two methods to our nominal atmospheric structure of GJ 436b: a thermal profile of GJ 436b (Fig. 1) determined with the code ATMO (Tremblin et al. 2015). For the opacities, we used the parameters of the benchmark study of Baudino et al. (2017) and assumed a composition at thermochemical equilibrium with solar elemental abundances (but with a reduction of 20% of oxygen due to sequestration). In our kinetic model, we used the same elemental abundances and assumed a vertical mixing $K_{zz}$ = 108cm2s*-1*. With the reduced scheme, we aim at being able to reproduce the abundances of a series of species of interest : H2O, CH4, CO, CO2, NH3, and HCN. These species are both the most abundant ones and all have strong infrared spectral signature. Note that we could add acetylene (C2H2) to this list, as this species is potentially important for C-rich atmospheres (e.g. Moses et al. 2013b; Venot et al. 2015). We chose instead to focus on species that have already been observed in exoplanets, which is not the case of acetylene currently, albeit it is present in Jupiter’s atmosphere (Nixon et al. 2010). As we will see in Sect. 4, the absence of this species in our reduced scheme limits the range of validity of the latter to atmospheres in which C2H2 is not very abundant: atmospheres with both a high C/O ratio (larger than 1) and a high temperature (T $\gtrsim$1000 K) will not be perfectly represented with the reduced chemical scheme; C-rich atmospheres with a lower temperature can be studied with the reduced scheme as acetylene is not a major component of these atmospheres (Venot et al. 2015). The methodology used in (DRGEP+SA) stops the reduction process when the error on targeted parameters goes beyond the user-specified tolerance level during removal of species from the master mechanism. For the set of species cited above, we have chosen a value of relative tolerance of 10%, globally representative of the uncertainty on the predicted abundance determined by our uncertainty propagation study (see Sect. 3.2). Our goal is to generate a reduced mechanism able to reproduce the prediction of concentrations given by the original mechanism, within the error bars. We want to obtain a reduced mechanism as small as possible and valid on a very large range of operating conditions, allowing to study various kind of atmospheres (Fig. 1). The original mechanism includes 957 reversible reactions, involving 105 neutral species (molecule or radical). The reduced mechanism includes 181 reversible reactions involving 30 different species: H2, H, H2O, CO, N2, CH4, CH3, NH3, NH2, HCN, OH, CO2, H2CO, HCO, 3CH2, 1CH2, O(3P), CH3OH, CH2OH, CH3O, H2CN, HNCO, HOCN, CN, NCO, NH, NNH, N2H2, N2H3, and He (the latter only being involved as a third-body).

In terms of temperature and pressure, the reduced chemical scheme has the same range of validity of our full chemical scheme (i.e. [300-2500] K and [0.01–100] bar). The temporal evolution of abundances with time are very similar with both chemical schemes and they are both able to reproduce thermochemical equilibrium (Fig. 2). Finally, we obtained a reduction factor of about 3 for the number of species and about 5 for the number of reactions, which leaves hope for a possible coupling of such mechanism in a 3D model. For a given atmosphere, using the reduced scheme instead of the full chemical scheme in our 1D kinetic model permit $\sim$30$\times$ faster runs. This reduced scheme is available in the KIDA Database (Wakelam et al. 2012)111kida.obs.u-bordeaux1.fr.

2.2 Uncertainty propagation model

To go further than a simple comparison of the vertical abundances given by the kinetic model when using the full and the reduced chemical schemes, we performed an uncertainty propagation study with the full scheme to have an objective reference. Indeed, the precision of the kinetic model is related to the uncertainties in the rate coefficients of the chemical scheme. Uncertainties of all reactions combine to each other and propagate non-linearly in the kinetic model through the system of differential equations, resulting from the continuity equations (for details on continuity equations, see e.g. Moses et al. 2011; Venot et al. 2012; Drummond et al. 2016). To study this uncertainty propagation, we use the method described in Dobrijevic & Parisot (1998); Dobrijevic et al. (2003), applied to a chemical scheme whose reactions are reversible. The rate constant of each forward chemical reaction was pseudo-randomly estimated according to its own temperature-dependent uncertainty factors $F(T)$ following the data evaluation and the numerical method described in Hébrard et al. (2015). The corresponding reverse rate constant was determined using the reaction equilibrium constant. This way, thermochemical equilibrium is maintained in the deep layers of the atmosphere, or in absence of disequilibrium processes. The sets of perturbed rates are used to run the kinetic model and determine ‘alternative’ atmospheric compositions for our model of GJ 436b. They represent the uncertainty on the predicted abundances. We have performed 1000 Monte-Carlo simulations allowing us to have statistically significant results. All simulations are run for $t$=109s, which is the time needed to reach steady-state for the atmosphere of GJ 436b. The uncertainty factors for chemical reactions included in our initial C0-C2 scheme have been published in Hébrard et al. (2015). Note that we do not consider any uncertainty that is potentially associated with the NASA polynomial coefficients (McBride et al. 1993), used to calculate the equilibrium constant. To our knowledge, no uncertainty on these coefficients have been evaluated so far (Goldsmith et al. 2012). Such ambitious work of bibliographic census and analysis would be necessary for a complete evaluation of the uncertainty on the predicted abundances. This is beyond the scope of this study that aims at developing a reduced scheme reproducing the equilibrium state determined with a more complete one and nominal values of NASA polynomial coefficients.

3 Results

3.1 GJ 436b model

As can be seen on Figs. 3, 4, and 5, the agreement between the atmospheric compositions obtained with the two chemical schemes is very good. Table 1 gathers, for each selected species, the maximum difference of abundances obtained using the two chemical networks and at which pressure level this value is reached. We focus our study on the [0.1–1000] mbar region only, which is the region probed by infrared observations. The species that presents the largest difference is CO2, with 0.02% at 300 mbar. The other species present differences always lower than this value.

3.2 Uncertainty propagation model

The 1000 Monte-Carlo runs performed for the atmospheric model of GJ 436b allow us to determine the distribution of abundance profiles of each species. It can be seen on Figs. 3, 4, and 5 that these distributions are highly species- and altitude-dependent. The plots on the left column of these figures represent the vertical abundances profiles, whereas the plots on the right column represent the histogram of abundances at 1 mbar. We chose this pressure level because it is located in the quenched zone and the abundance is representative of the [0.1–1000] mbar region probed by infrared observations. For all species, thermochemical equilibrium is maintained in the deep atmosphere. The quenching area (proper to each species) remains the same but the exact quenching pressure level is slightly modified, leading to variations above this point, at lower pressures. For H2O, CH4, and NH3, the width of the distribution of abundance is very small, whereas this width is more important for species such as CO, CO2, and HCN. It is interesting to note that some important species, although not in the list of species that must be well described by the reduced scheme, are however well reproduced by the reduced scheme, e.g. H, N2, CH3.

We plot on each histogram the standard deviation ($\sigma$) of the distribution. For each species of the reduced scheme, the reduced chemical scheme gives an abundance well within the 1-$\sigma$ interval around the nominal abundance obtained with the full C0-C2 chemical scheme.

4 Range of validity of the reduced chemical scheme

The reduced chemical scheme has been designed originally for the atmosphere of GJ 436b, with solar elemental abundances. However, it turns out that this chemical scheme has a very large range of validity. We have used it to model GJ 436b with different atmospheric compositions in terms of metallicities and C/O ratios (Sect. 4.1). We have also modelled planets that have been studied with the C0-C2 scheme in the past: in Sect. 4.3 we show results for the hot Jupiters HD 209458b and HD 189733b (Venot et al. 2012). In Sect. 4.4, we show the results for the deep atmosphere of Uranus and Neptune (Cavalié et al. 2017) and in Sect. 4.5, we present the results obtained for SD 1110-like brown dwarf, modelled for the first time with the C0-C2 scheme. Finally, we explore a range where the reduced scheme is not valid anymore, by modelling HD 209458b with high C/O ratios (3 and 6 times solar) in Sect. 4.6. All the thermal profiles are presented in Fig. 1. It can be viewed that the range of temperatures thus scanned is large: $\sim$[300–2500] K. For each case, the kinetic code is run until steady-state of the atmosphere is reached (i.e. 107, 108, or 109s depending on the planet) and chemical composition obtained with the two schemes are compared.

4.1 GJ 436b with different elemental abundances

We modelled GJ 436b with the same thermal profile than in Sect. 3.1, but with different elemental abundances. On one hand, keeping the C/O/N ratios solar, we increase the metallicity by a factor of 10 and 100. In these cases, as can be seen on Fig. 6, the abundances of all species, except H, globally increase compared to the nominal case (Figs. 3, 4, 5). Note that an increase of the metallicity has an effect on the thermal profile (e.g. Charnay et al. 2015), that we didn’t take into account here as our goal is only to compare the outputs of the two chemical schemes. On the other hand, keeping O/H and N/H ratios equal to the solar values, we modify only the C/H ratio by a factor of 0.5 and 3 relatively to the solar value. One can observe on Fig. 6 that a decrease of the C/H ratio, leads to a decrease of the abundance of CH4, CO, CO2, and HCN. Inversely, an increase of C/H ratio leads to an increase of the abundance of these carbonaceous species. Methane becomes more abundant than water in this case of carbon-rich atmosphere. These evolutions of abundances with the metallicity and the elemental C/O ratio have been shown also by e.g. Moses et al. 2013a; Venot et al. 2014, 2015; Tsai et al. 2017.

The reduced chemical scheme describes very well the atmospheric composition of the four aforementioned models. The results obtained with the full chemical scheme and the reduced one are not distinguishable on log-log plots (Fig. 6). On Table 2, we report the maximum and mean differences of abundances between the two chemical schemes for these four models, within the pressure range of interest [0.1–1000] mbar. For each species, we observe that the maximal departures are reached for the high metallicity case. For this case, the maximal difference is 0.6% and is reached by CH4 at 600 mbar. For the other cases and species, deviations are always lower than this value.

We have also tested much higher metallicities for the atmosphere of GJ 436b: 500, 1000, and 10 000 $\times$ solar metallicities (not shown here). For all of these models, we have found a very good agreement between the reduced and the full chemical scheme, with maximal difference of CH4 between the two schemes of 1, 2, and 11 % for the cases 500, 1000, and 10 000 $\times$ solar metallicities, respectively (in the pressure range [0.1–1000] mbar). The other species present variations less than 0.5%.

4.2 Another warm Neptune: GJ 1214b

We have modelled another warm Neptune, GJ 1214b, with the average dayside profile calculated with another atmospheric model, the Generic LMDZ GCM (Charnay et al. 2015). We have assumed an atmospheric metallicity of 100 $\times$ solar metallicity (Lodders 2010) and used the corresponding parametrization of vertical mixing recommended by Charnay et al. (2015): $K_{zz}=3\times 10^{7}\times P^{-0.4}$cm2s*-1*, with $P$ in bar. Even if the thermal profile of GJ 1214b is close to that of GJ 436b (see Fig. 1), the difference of temperature ($\sim$100 K) in the quenching area (104 mbar) leads to important differences in chemical composition between the two planets while having the same elemental abundances. Indeed, with a metallicity 100 $\times$ solar, the transition of C-bearing species between CO and CH4 occurs at these temperature and pressure values ($\sim$1100 K and $\sim$10 bar). Whereas in GJ 436b’s atmosphere, CO is quenched when more abundant than CH4, in GJ 1214b, it is the opposite. Our model of GJ 1214b’s atmosphere, in complement of that of GJ 436b with a metallicity 100 $\times$ solar, shows that the reduced chemical scheme is accurate in these two different atmospheric configurations (see Fig. 7). In the pressure range of interest [0.1–1000] mbar, departures between the abundances obtained with the full and the reduced schemes are lower than 1% for all species for which the reduced scheme is designed (see Table 3). For this planet, CO and HCN are the two species that show the strongest variations in abundance between the two schemes. Note that these differences of temperature and thus of chemical compositions result more from differences in the atmospheric models (ATMO and Generic LMDZ GCM) used to calculate the thermal profiles (e.g. internal temperature, opacities, treatment of alkali) rather than differences between the two planets which are quite similar in terms of irradiation. It has been shown by Baudino et al. (2017) that different assumptions lead to important variations in the results of atmospheric models. Our results show that the chemical composition is consequently also sensitive to these hypotheses.

4.3 Hot Jupiters: HD 209458b and HD 189733b

We have tested the validity of the reduced chemical scheme on the famous hot Jupiters HD 209458b and HD 189733b. We have used the thermal profiles and eddy diffusion coefficients profiles of Moses et al. (2011) that were also used in Venot et al. (2012) with the full chemical scheme. We assume solar elemental abundances (Lodders 2010) with a reduction of 20% of oxygen, due to sequestration in the deep layers. The agreement between the results obtained with the two schemes is excellent (see Fig. 8). For the pressure range probed by observations, we listed in Table 3 the maximum difference of abundances between the two chemical schemes for these two models. Contrary to the other models presented earlier, we observe for HD 209458b some relatively important deviations between the two schemes for pressures lower than 4 mbar for HCN (60%), CH4 (600%), and NH3 (2000%). However, these deviations are acceptable as 1) they are at the limit of the pressure range probed by observations and 2) they concern species with a low abundance (¡10*-6*). Then, at even higher altitudes (from 10*-2* mbar), we observe deviations for CO2 and H2O (H and H2 also but they are not in the list of species for which the reduced scheme is designed). However, these pressure levels are not probed by observations. For HD 189733b, the agreement between the two chemical schemes is higher. The variations of abundances for species of interest are lower than 1% in the [0.1–1000] mbar region. The largest variation is due to HCN at 100 mbar.

4.4 Deep atmospheres of Uranus and Neptune

We have applied our kinetic model to the deep atmospheres of Uranus and Neptune, as in Cavalié et al. (2017), using the ‘3-layer’ thermal profiles calculated with the prescription of Leconte et al. (2017). Following Cavalié et al. (2017), we have set the eddy diffusion coefficient $K_{zz}$ to 108 cm2.s*-1*. Concerning the elemental abundances, we have assumed the following enrichment compared to the solar abundances (Lodders 2010): For both planets, we have assumed a N/H ratio of 0.5 times the solar value. For Uranus, we have assumed C/H and O/H ratios of 75 and 160 times the solar value, respectively. For Neptune, we have assumed C/H and O/H ratios of 45 and 480 times the solar value, respectively. These values correspond to the nominal cases of Cavalié et al. (2017), allowing to reproduce the observational constraints available at 2 bar for CO and CH4: y${}_{CH_{4}}$=0.04 for both planets and yCO=2.0$\times$10*-7* at Neptune and yCO¡2.1$\times$10*-9* at Uranus.

The tropospheric compositions determined with the full and reduced schemes are presented Fig. 9. Chemical abundances obtained with both schemes are very close. We report in Table 4 the maximum variations of abundances between the two chemical schemes for the models of Uranus and Neptune in the pressure range [106–107] mbar. For the giant planets, we choose this range because it includes the quenching level. In the methodology we apply to the study of these planets, this level is the most important and decisive for the conclusions drawn concerning the elemental composition of the atmospheres (Cavalié et al. 2014, 2017). For the two planets, the variations of abundances between the two chemical scheme are low ($<$ 10%) in the pressure range of interest. We remark that the maximum values are found for Neptune, which is not very surprising. Indeed, this atmosphere is the most enriched of this study (i.e. O/H and C/H are, respectively, 480 and 40 $\times$ the solar values), and thus the elemental conditions are the most distant from that of our nominal GJ 436b model used to develop the reduced chemical scheme.

4.5 Atmosphere of a SD 1110-like brown dwarf

We have finally applied our kinetic model to study the chemical composition of a brown dwarf’s atmosphere. With ATMO (Tremblin et al. 2016; Drummond et al. 2016), we have determined the thermal profile of a brown dwarf located at the L/T transition, with an effective temperature of 1000 K and log $g$=4.5. At the transition, L-type brown dwarfs are CO-dominated and become CH4-dominated as they cool down towards the T-dwarf sequence. Non-equilibrium chemistry at the CO/CH4 transition is thus an important process to correctly predict the transition as a function of effective temperature in the cooling sequence of brown dwarfs (Saumon et al. 2003). The model used for this study is similar to the T5.5 brown dwarf SD 1110 (see Stephens et al. 2009). The thermal structure has been calculated assuming thermochemical equilibrium with solar elemental abundances (Lodders 2010). The abundance of condensates are explicitly included in this calculation, which corresponds to a reduction of $\sim$20% of oxygen. For the kinetic model, we thus applied the same oxygen reduction and we have taken a constant vertical mixing of 108cm2s*-1*. Fig. 10 shows that the chemical composition obtained with the reduced chemical scheme is similar to that obtained with the full scheme. For brown dwarfs, the pressure range that can be probed by infrared observations is slightly deeper than for exoplanets. Thus, we have focused our study here to the [102–104] mbar range (see Morley et al. 2014). In this range, the largest variation between abundances obtained with the two schemes is of 3% and is due to HCN (at 2 bar) and CO2 (at 100 mbar). H2O and CH4 present variations of about the same order of magnitude (2% at 9 bar) (see Table 4).

4.6 Out of validity: C-rich hot Jupiter

As acetylene is not included in our reduced chemical scheme, the study of hot atmospheres rich in carbon might be problematic. To assess the extent of this error, we modelled HD 209458b, with our two chemical schemes, assuming two carbon-rich compositions: 3 $\times$ and 6 $\times$ the solar C/O ratio (and then removing 20% of oxygen). One can see on Fig. 11 that deviations between the abundances obtained with the two schemes appear for pressures lower than $\sim$10 mbar. These differences are due to the important mixing ratio of C2H2; species that is included only in the full chemical scheme. The high amount of carbon carried by this species in the C0-C2 scheme is distributed in the other carboneous species in the reduced scheme, mainly CH4, HCN, and CH3. For species of interest, the variations of abundances obtained with the two schemes reach very important values, up to 8$\times$105 % (see Table 5). Such differences of abundances can have consequences on the synthetic spectra computed using these results, especially because of the absence of C2H2 as a source of opacity when using the reduced scheme’s chemical composition. These results make us recommend not to use the reduced scheme to study hot carbon-rich atmosphere, or at least to be very cautious with the chemical composition predicted in the upper atmosphere (P¡10 mbar).

5 Conclusion

We have developed a reduced chemical scheme for exoplanetary atmospheres from the C0-C2 chemical network of Venot et al. (2012). This reduced scheme has been designed to reproduce the abundances of H2O, CH4, CO, CO2, NH3, and HCN and is available in the KIDA Database (Wakelam et al. 2012). For these species of interest, we have validated this new scheme through a comparison with the abundances predicted by the full scheme. In many different cases (metallicities, C/O ratios, thermal profiles) the abundances are very similar, with only some percent of deviations (often less than 1%) in the pressure range probed by infrared observations. The only limitation of this scheme concerns hot C-rich atmospheres. As the reduced chemical scheme does not include C2H2, the use of this scheme to study this kind of atmospheres (i.e. T$\gtrsim$ 1000K and C/O ¿ 1) is not recommended. For our model of GJ436b, we have also performed an uncertainty propagation study on the full chemical scheme. The abundances obtained with the reduced scheme are included in the 1-$\sigma$ interval around the C0-C2 nominal abundance.

An even more drastic reduction of the scheme is possible, decreasing the desired level of accuracy in our methodology. However, such a reduction would lead to a chemical scheme giving deviations of abundances with the full scheme larger than what we obtained with the scheme presented in this paper. Also, the range of validity in terms of temperature of this very reduced scheme might be narrowed down, which could be problematic when used in 3D models of planets with a large day/night temperature gradient, such as the hot Jupiters WASP-43b, WASP-18b, WASP-103b, WASP-12b (Parmentier & Crossfield 2017).

The advantage of this reduced network is that it runs faster than the full chemical scheme from which it is extracted ($\sim$30 times faster). This gain in computational time is due to the reduction in the size of the system of differential equations (continuity equations) that must be solved to determine the steady state of the atmosphere. But the most important application we expect for this scheme is to be incorporated in 3D models. This would be an important step forward in view of a good interpretation of the future observations that will be provided by the next generation of telescopes (JWST, ARIEL) and for understanding the global chemical composition of exoplanets, brown dwarfs and the deep composition of solar system giant planets. For these latters, the recent Juno/MWR data (Bolton et al. 2017), presenting an unexpected NH3 distribution below the condensation level, prove that robust 3D models accounting for chemistry and dynamics are highly needed to decipher their composition. These 3D models of the tropospheres of the solar system giant planets will be crucial to help select the latitude of any future entry probe (Arridge et al. 2014; Mousis et al. 2014, 2016, 2018).

Acknowledgements.

The authors thanks the referee for his careful reading of the manuscript which allowed its improvement. O.V. and T.C. thanks the CNRS/INSU Programme National de Planétologie (PNP) for funding support. PT acknowledges support from the European Research Council (grant no. 757858 – ATMO). BD and EH acknowledge support from the STFC Consolidated Grant ST/R000395/1.