A CHEMICAL KINETICS NETWORK FOR LIGHTNING AND LIFE IN PLANETARY ATMOSPHERES

Two-body neutral–neutral and ion–neutral reactions follow the basic scheme:

$$\begin{eqnarray}&&{\rm{A}}+{\rm{B}}\to {\rm{Y}}+{\rm{Z}},\;\mathrm{and}\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{\rm{A}}}^{+}\;+\;{\rm{B}}\to {{\rm{Y}}}^{+}\;+\;{\rm{Z}}.\end{eqnarray} \tag{ 4 }$$

The rate constants for these reactions are approximated by the Kooij equation (Kooij 1893):

$$\begin{eqnarray}&&{k}_{2}=\alpha {\left(\displaystyle \frac{T}{300{\rm{K}}}\right)}^{\beta }\;{e}^{-\gamma /T},\end{eqnarray} \tag{ 5 }$$

where T [K] is the gas temperature³ , k₂ [cm³ s⁻¹] is the rate constant, and α [cm³ s⁻¹], β, and γ are constants characterizing the reaction. All of these reactions are reversed in our network and we use the rate coefficients for the best characterized direction for each reaction, which is typically the exothermic direction. For neutral–neutral reactions, even when exothermic, there is often a sizable barrier to reaction, allowing certain elements to be locked into non-equilibrium configurations at low temperatures effectively for eternity, because the barrier to the lower energy state is too large to be overcome in the current environment.

Ion–neutral reactions do not typically have barriers in the exothermic direction, and in many cases the rate constants are altogether temperature-independent, closely approximating the Langevin approximation. A notable exception are charge-exchange reactions,

$$\begin{eqnarray}&&{{\rm{A}}}^{+}\;+\;{\rm{B}}\to {{\rm{B}}}^{+}\;+\;{\rm{A}},\end{eqnarray} \tag{ 6 }$$

which, due to the differences in energy between ionic and neutral ground states, often contains barriers on the order of a few × 100 K.

The rate constants for the forward reactions are given in Appendix A with the label "2n," reactions 577–1352. These reactions are reversed following the scheme described in Appendix B. The ion–neutral reactions are also reversed, and are listed in Appendix A with "2i," reactions 1353–2569.

Reactions that involve a third body occur primarily in the two forms:

$$\begin{eqnarray}&&{\rm{A}}+{\rm{M}}\to {\rm{Y}}+{\rm{Z}}+{\rm{M}},\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{\rm{A}}+{\rm{B}}+{\rm{M}}\to {\rm{Z}}+{\rm{M}},\end{eqnarray} \tag{ 8 }$$

where M represents any third body. Decomposition reactions are well studied at high temperatures, being important for various combustion processes. Just as in Section 2.1, we choose the reactions best characterized, which in this case often involve endothermic reactions. The rate coefficients for the majority of these reactions follow the Lindemann form (Lindemann et al. 1922). In this form, we first determine the rate constants in the low-pressure (k₀ [cm⁶ s⁻¹]) and high-pressure $({k}_{\infty }\;$[cm³ s⁻¹]) limits:

$$\begin{eqnarray}&&{k}_{0}={\alpha }_{0}{\left(\displaystyle \frac{T}{300{\rm{K}}}\right)}^{{\beta }_{0}}\;{e}^{-{\gamma }_{0}/T},\end{eqnarray} \tag{ 9 }$$

$$\begin{eqnarray}&&{k}_{\infty }={\alpha }_{\infty }{\left(\displaystyle \frac{T}{300{\rm{K}}}\right)}^{{\beta }_{\infty }}\;{e}^{-{\gamma }_{\infty }/T}.\end{eqnarray} \tag{ 10 }$$

These are combined with the number density of the neutral third species, [M] [cm⁻³] to determine the reduced pressure, ${p}_{r}={k}_{0}[{\rm{M}}]/{k}_{\infty }$, and this can then be utilized to set the pressure-dependent effective "two-body" rate:

$$\begin{eqnarray}&&{k}_{2}=\displaystyle \frac{{k}_{\infty }{p}_{r}}{1+{p}_{r}}.\end{eqnarray} \tag{ 11 }$$

Sometimes this expression is multiplied by a dimensionless function $F(p,T)$ to more accurately approximate the transition between the low-pressure and high-pressure limits, and this provides the Troe form (Troe 1983). The coefficients for the Troe form are not explicitly given.

We favor using the rate constants for three-body combination reactions, and reversing these reactions to determine the rate of thermal decomposition. In many cases, however, the rate constants are unavailable. When we have only the rate coefficients for the decomposition reactions, we add an additional 500 K barrier to both the decomposition and three-body combination rate constants. This barrier is added in order to limit runaway three-body reactions that can result from reversing decomposition reactions at low temperatures.

Additionally, we incorporate a small number of radiative association reactions, of the form:

$$\begin{eqnarray}&&{\rm{A}}+{\rm{B}}\to {\rm{Z}}\;+\;\gamma ,\end{eqnarray} \tag{ 12 }$$

where γ is the radiated photon that carries the excess energy from the association. We appropriate the Kooij form for this reaction, as with two-body neutral–neutral reactions, in order to determine the rate constant ${k}_{\mathrm{ra}}$ [cm³ s⁻¹]. We then apply this rate constant, along with the rate constant for the corresponding three-body reaction, to the adduct form of the overall rate constant (Hébrard et al. 2013, their Equation (B.2)):

$$\begin{eqnarray}&&k=\displaystyle \frac{({k}_{0}[{\rm{M}}]F+{k}_{r}){k}_{\infty }}{{k}_{0}[{\rm{M}}]+{k}_{\infty }},\end{eqnarray} \tag{ 13 }$$

where the function F is from the Troe form of the transition from high to low pressure.

The rate constants for the forward reactions are given in Appendix A with the labels "2d" for the neutral species and "3i" for ion–neutral species. These reactions are reversed in the manner described by Appendix B. Reactions 1–420 are reactions of this type, for which each odd number reaction gives the low-pressure rate constant k₀ [cm⁶ s⁻¹] and each even number reaction gives the high-pressure rate constant ${k}_{\infty }$ [cm³ s⁻¹]. Reactions labeled "ra" are radiative association reactions, numbered 2974–2980.

A special set of three-body reactions are thermal ionization and three-body recombination reactions, which proceed by the pair of equations (analogous to Equations (7) and (8)):

$$\begin{eqnarray}&&{\rm{A}}+{\rm{M}}\to {{\rm{Z}}}^{+}+{e}^{-}+{\rm{M}},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{{\rm{A}}}^{+}+{e}^{-}+{\rm{M}}\to {\rm{Z}}+{\rm{M}}.\end{eqnarray} \tag{ 15 }$$

For which we again use published rates wherever possible for the ionization reactions, (Equation (14)), but in many cases here use the simple approximation:

$$\begin{eqnarray}&&{k}_{0}={\left(\displaystyle \frac{8\pi {e}^{8}}{{m}_{e}{k}_{B}T}\right)}^{1/2}{e}^{-I/{k}_{B}T},\end{eqnarray} \tag{ 16 }$$

where $e=4.9032\times {10}^{-10}$ esu is the elementary charge, ${k}_{B}=1.38065\times {10}^{-16}$ erg K⁻¹ is the Boltzmann constant, ${m}_{e}=9.1084\times {10}^{-28}$ g is the mass of the electron, and I is the ionization energy (here in units of erg) which we determine from the change in the Gibbs free energy for the reaction. ${k}_{\infty }$ is then estimated from k₀.

Three-body recombination and ionization reactions have been well studied, and in many cases have well characterized rate constants. Here we treat the three-body recombinations as the reverse reactions for the collisional ionization reactions, but the studied rate coefficients for these reactions generally have a temperature dependence of ${T}^{-4.5}$, at least for $T\gt 1$ K (Hahn 1997). This creates a problem for reversibility. Using these rates will not allow us to reproduce chemical equilibrium for plasmas and this is largely because we are not properly treating the time-dependent plasma conditions in which these rates are often measured. Many of these rate constants may accurately describe the time to achieve an equilibrium electron density in a regime where a strong ionizing source has recently been removed from the environment.

With this in mind, we instead set the recombination rate constants such that, when dissociative recombination reactions are disabled, the Saha equation is upheld.

These reactions and rate coefficients are also given in Appendix A. The ionization reactions are labeled "ti" and numbered 421–576. As with Section 2.2, the odd reactions are k₀ [cm⁶ s⁻¹] and the even numbers are ${k}_{\infty }$ [cm³ s⁻¹].

Finally, we include a series of dissociative recombination reactions, which take the form:

$$\begin{eqnarray}&&{{\rm{A}}}^{+}+{e}^{-}\to {\rm{Y}}+{\rm{Z}}.\end{eqnarray} \tag{ 17 }$$

These have rate constants parameterized in the form of Equation (5). The reverse reactions can in principle be calculated, and their rate constants could be calculated straightforwardly using the same principles used for the three-body reactions. This would effectively be analogous to the rates of three-body recombination for any third body, and we do not find that reversing these reactions changes the results much. When we compare with chemical equilibrium, however, we disable these reactions. The dissociative recombination reactions are taken only from the OSU 09 2010 high-temperature network (Harada et al. 2010), and shown in Appendix A, numbered 2777–2973, and labeled "dr."

Photochemistry is considered for the species H, H⁻, He, C, C(¹D), C(¹S), N, O, O(¹D), O(¹S), H⁻, C₂, CH, CN, CO, H₂, N₂, NO, O₂, OH, CO₂, H₂O, HO₂, HCN, NH₂, NO₂, O₃, C₂H₂, H₂CO, H₂O₂, NH₃, NO₃, CH₄, HCOOH, HNO₃, N₂O₃, C₂H₄, C₂H₆, CH₃CHO, C₄H₂, C₄H₄, Na, K, and HCl. The photoionization and photodissociation cross sections are taken almost entirely from PhIDRates ⁴ (Huebner & Carpenter 1979; Huebner et al. 1992; Huebner & Mukherjee 2015), with the exception of C₄H₂, C₄H₄, and N₂O₃, the cross sections of which are taken from the MPI-Mainz UV/VIS Spectral Atlas⁵ (Keller-Rudek et al. 2013).

We divide the cross sections between 200 bins each $\approx 50$ Å wide. A comparison between our binned cross sections and the raw cross sections from PhIDRates is plotted for an example reaction (Figure 1). The cross sections, both in the database and here are of the form $\sigma (\lambda )$ with σ in units cm² and wavelength in units of Å. The resolution for the UV cross sections is fairly low, and cannot encapsulate the fine structure of the UV emission lines or the UV cross sections. This is especially important when treating ionospheres of gas giants, since, e.g., the fine structure in the H₂ bands leave small spectral windows through which photons can penetrate and effectively ionize deeper in the atmosphere. Such a low-resolution spectrum will effectively close these windows and underestimate the ion production in the ionosphere (Kim & Fox 1994; Kim et al. 2014). High resolution is also a important for capturing where the UV flux and cross sections both peak; a low-resolution cross section can, in this case, underestimate the destruction rate of the species with this resonant photochemical cross section. As can be seen below, these issues do not significantly affect the comparisons of this model for HD 209458b, Jupiter, or Earth. For photoionization deep in the atmosphere, where high resolution is essential, the network itself need not be modified. The transport of UV photons line by line would need to be calculated.

Zoom In Zoom Out Reset image size

The tabulated chemical cross sections are combined with $F(\lambda ,z)$ [photons cm⁻² s⁻¹ Å⁻¹], the radiant flux density onto a unit sphere (hereafter called the actinic flux) located at atmospheric height, z [cm], to determine the photochemical rate constants,

$$\begin{eqnarray}&&{k}_{\mathrm{ph},i}(z)={\tau }_{f}{\displaystyle \int }_{1\;\mathring{\rm{A}} }^{{10}^{4}\;\mathring{\rm{A}} }{\sigma }_{i}\;F(\lambda ,z)\;d\lambda ,\end{eqnarray} \tag{ 18 }$$

where i is indexed over the molecules listed above, for which photochemistry is considered. ${\tau }_{f}$ is a dimensionless parameter representing the fraction of time (over a period much longer than the longest characteristic timescale for the atmosphere) the particular atmospheric region is irradiated; for tidally locked planets, ${\tau }_{f}=1$ (dayside) or 0 (nightside), the diurnal average for a rotating planet is ${\tau }_{f}=1/2$. The photoionization and photodissociation reactions are listed in Appendix A, reactions numbered 2570–2693, and labeled "pi" for photoionization reactions and "pd" for photodissociation reactions.

Cosmic-ray ionization and dissociation is parameterized by ζ (Rimmer & Helling 2013), to treat both direct ionization by galactic cosmic rays and ionization by secondary particles produced in air showers. The cosmic-ray ionization rate depends on the chemical species in question, since different species will have different chemical cross sections for the photons produced by cosmic rays, and this is accounted for by multiplying $\zeta (z)$ by a constant ${\kappa }_{\mathrm{CR},i}$ such that:

$$\begin{eqnarray}&&{k}_{\mathrm{CR},i}(z)={\kappa }_{\mathrm{CR},i}\zeta (z).\end{eqnarray} \tag{ 19 }$$

We treat low-energy cosmic rays ($E\lt 1$ GeV) for these objects as though they have been significantly shielded by the astrospheres of the host stars, and therefore set the fitting parameters for the incident cosmic-ray flux to $\alpha =0.1$ and $\gamma =-1.3$ in the equation for the flux of cosmic-ray particles:

$$\begin{eqnarray}j(E)=\left\{\begin{array}{ll}j({E}_{1}){\left(\displaystyle \frac{p(E)}{p({E}_{1})}\right)}^{\gamma }, & \mathrm{if}\;E\gt {E}_{2}\\ j({E}_{1}){\left(\displaystyle \frac{p({E}_{2})}{p({E}_{1})}\right)}^{\gamma }\;{\left(\displaystyle \frac{p(E)}{p({E}_{2})}\right)}^{\alpha }, & \mathrm{if}\;{E}_{\mathrm{cut}}\lt E\lt {E}_{2}\\ 0, & \mathrm{if}\;E\lt {E}_{\mathrm{cut}}\end{array}\right.\end{eqnarray} \tag{ 20 }$$

where $p(E)=\tfrac{1}{c}\sqrt{{E}^{2}\quad +\quad 2{{EE}}_{0}}$, ${E}_{0}=9.38\times {10}^{8}\;\mathrm{eV}$, ${E}_{1}={10}^{9}\;\mathrm{eV}$, and ${E}_{2}=2\times {10}^{8}\;\mathrm{eV}$, and the flux at E₁ is set to $j({E}_{1})=0.22$ cm⁻² s⁻¹ sr⁻¹ (GeV/nucleon)⁻¹. All of these parameters except α are observationally well constrained (Indriolo et al. 2009). For a demonstration of how α affects the cosmic-ray spectrum, and a discussion of the Monte Carlo transport we use for cosmic rays of energy $\lt 1\;\mathrm{GeV}$, see Rimmer et al. (2012) and Rimmer & Helling (2013). For ionization rate by cosmic rays of energy $\gt 1\;\mathrm{GeV}$, ${Q}_{\mathrm{HECR}}$ [cm⁻³ s⁻¹], we use the analytical method of Velinov & Mateev (2008).

Cosmic-ray reactions are listed in Appendix A, numbered 2694–2776, and labeled "cr."

At sufficiently high temperatures and pressures, a gas should rapidly settle into chemical equilibrium. An important test for a chemical network is that its steady state solution converges to the chemical equilibrium solution. To perform this test of our network, we solve the chemical kinetics at a single ($T,p$) point, using the rate constants from the Stand2015 network, disabling the cosmic-ray reactions, photochemistry, and dissociative recombination. We compute a time-dependent solution of the equation

$$\begin{eqnarray}&&\displaystyle \frac{{{dn}}_{i}}{{dt}}={P}_{i}-{L}_{i}.\end{eqnarray} \tag{ 21 }$$

We solve this equation for T = 1000 K and p = 1 bar, with solar abundances from Asplund et al. (2009). We compare our results to chemical equilibrium calculations using the Burcat polynomials (Burcat & Ruscic 2005), and plot our comparisons in Figure 2 and find excellent agreement. This agreement is not surprising; we have used the same thermochemical data to reverse our reactions, and only include reversed reactions in this test, so once the system achieves steady state, computationally achievable at this pressure and temperature, the chemistry has effectively settled into equilibrium.

Zoom In Zoom Out Reset image size

We also compare our electron number density to the electron number density achieved using the Saha equation, this time at a pressure of 10⁻⁴ bar and over a range of temperatures from 1000 to 10,000 K. This comparison is plotted in Figure 3. The comparison is virtually perfect when $T\gtrsim 2000$ K, unsurprising given the way the three-body recombination reactions are calculated (see Section 2.3). At ∼1000 K, our results diverge from the Saha equation. This is because the integrator does not reliably calculate mixing ratios below $\sim {10}^{-30}$. Indeed, at this stage, the electron number density achieves $\sim {10}^{-300}$ cm⁻³ while the H⁺ number density rests at $\sim {10}^{-60}$ cm⁻³, producing significant charge balance errors. These large errors in the charge balance fluctuate, and only appear when the ionization fraction is $\lesssim {10}^{-30}$, at which point ion–neutral chemistry is inconsequential.

Zoom In Zoom Out Reset image size

The coupled 1D continuity equations describing the time-dependent vertical atmospheric chemistry are

$$\begin{eqnarray}&&\displaystyle \frac{\partial {n}_{i}}{\partial t}={P}_{i}-{L}_{i}-\displaystyle \frac{\partial {{\rm{\Phi }}}_{i}}{\partial z},\end{eqnarray} \tag{ 22 }$$

where n_i [cm⁻³] is the number density of species i, and $i=1,\ldots ,{N}_{s}$, and N_s is the total number of species. P_i [cm⁻³ s⁻¹] is the rate of production and L_i [cm⁻³ s⁻¹] is the rate of loss of species i. The rightmost term is the vertical change in flux ${{\rm{\Phi }}}_{i}$ [cm⁻² s⁻¹] and represents the flux due to both eddy (K [cm² s⁻¹]) and molecular diffusion (D [cm² s⁻¹]),respectively, related as (Banks & Kockarts 1973, their Equation (15.14)),

$$\begin{eqnarray}{{\rm{\Phi }}}_{i} & = & -K\left[\displaystyle \frac{\partial {n}_{i}}{\partial z}+{n}_{i}\left(\displaystyle \frac{1}{{H}_{0}}+\displaystyle \frac{1}{T}\displaystyle \frac{{dT}}{{dz}}\right)\right]\\ & & -D\left[\displaystyle \frac{\partial {n}_{i}}{\partial z}+{n}_{i}\left(\displaystyle \frac{1}{{H}_{i}}+\displaystyle \frac{1+{\alpha }_{T}}{T}\displaystyle \frac{{dT}}{{dz}}\right)\right],\end{eqnarray} \tag{ 23 }$$

where H₀[cm] is the pressure scale height of the atmosphere at z [cm], H_i[cm] is the molecular scale height of the atmosphere for species i, and ${\alpha }_{T}$ is the thermal diffusion factor (Banks & Kockarts 1973; Yung & Demore 1999; Zahnle et al. 2006; Hu et al. 2012). For molecular diffusion coefficients, we adopt the Chapman–Enskog theory (Enskog 1917; Chapman & Cowling 1991). Eddy diffusion coefficients are either determined empirically, as with Earth and Jupiter, or are derived from global circulation models, as is the case for HD 209458b.

In Equation (23), the terms dealing with eddy diffusion and molecular diffusion are separated out, clarifying the four regions that Equations (22) and  (23) describe. (1) Deep within the atmosphere, where pressures and temperatures are sufficiently large, the thermochemistry dominates, and the equation simplifies to Equation (21). The atmospheric chemical composition converges to chemical equilibrium or at least to some stable quasi-equilibrium. (2) Higher in the atmosphere, the eddy diffusion may dominate, and the species are quenched. Their abundance is mixed evenly over a wide range of the atmosphere at timescales shorter than the chemical timescales. (3) Above this region, molecular diffusion may dominate, and at that point, species lighter than the mean molecular mass of the atmospheric gas will rise up, and species heaver than the mean molecular mass will settle down, and the chemistry will largely be determined by the individual scale heights of the atmospheric constituents. (4) Non-equilibrium processes, such as photochemistry or cosmic-ray chemistry, may create a fourth region, the composition of which is determined by irreversible chemical reactions.

Since the purpose of this paper is to introduce a new chemical kinetics network for lightning and prebiotic processes, our focus is not on the atmospheric dynamics (for this, see Lee et al. 2015). We therefore apply a simple approximation to Equation (22), inspired by Alam & Lin (2008). We first cast Equation (22) in a Lagrangian formulation, and consider Eddy diffusion to be moving small parcels of the gas vertically. We follow a single parcel as it moves up from the lower boundary of the temperature profile, and then returns down again. In reality, the parcel would be jostled in all three dimensions as it makes a complex journey up to the top of the atmosphere, but 1D transport models are unable to capture this effect in full.

The differential diffusion of molecules into and out of the parcel requires a different approach. The discrete formulas used by Hu et al. (2012, their Equation (9)) in the Lagrangian frame are

$$\begin{eqnarray}\displaystyle \frac{\partial {n}_{i,j}}{\partial t} & = & {P}_{i,j}-{L}_{i,j}{n}_{i,j}-{d}_{j+1/2}\displaystyle \frac{{n}_{\mathrm{gas},j+1/2}}{{n}_{\mathrm{gas},j+1}}{n}_{i,j+1}\\ & & -\left({d}_{j+1/2}\displaystyle \frac{{n}_{\mathrm{gas},j+1/2}}{{n}_{\mathrm{gas},j}}-{d}_{j-1/2}\displaystyle \frac{{n}_{\mathrm{gas},j-1/2}}{{n}_{\mathrm{gas},j}}\right){n}_{i,j}\\ & & +{d}_{j-1/2}\displaystyle \frac{{n}_{\mathrm{gas},j-1/2}}{{n}_{\mathrm{gas},j-1}}{n}_{i,j-1}.\end{eqnarray} \tag{ 24 }$$

Here, j represents the parcel being followed, $j-1$ the parcel directly beneath j, $j+1$ the parcel above j, and $j\pm 1/2$ an arithmetic average between j and $j\pm 1$. n without any i subscript represents ${n}_{\mathrm{gas}}$ at the relevant parcel, and

$$\begin{eqnarray}&&{d}_{j\pm 1/2}=\displaystyle \frac{{D}_{j\pm 1/2}}{2{({\rm{\Delta }}z)}^{2}}\left[\displaystyle \frac{(\bar{m}-{m}_{i})g{\rm{\Delta }}z}{{k}_{B}{T}_{j\pm 1/2}}-\displaystyle \frac{{\alpha }_{T}}{{T}_{j\pm 1/2}}({T}_{j\pm 1}-{T}_{j})\right].\end{eqnarray} \tag{ 25 }$$

$\bar{m}\;$[g] denotes the mean molecular mass of the atmosphere at z and m_i[g] the mass of species i.

Both the third and last terms on the right-hand side of Equation (24) do not depend on n_i and can therefore be treated as source terms, P_i. The fourth term can be treated as a term in L_i, such that molecules "destroyed" by this reaction are "banked," ${\rm{A}}\to \mathrm{BA}$. The "banked" molecules re-enter the parcel at a rate determined by the third and last terms on the right-hand side of the equation, thus conserving mass throughout the parcel's travels. Violations of this conservation do not appear here, but can be accounted for via further reactions, settling, condensation and evaporation, outgassing and escape, discussed in Appendix C. Although it is straightforward to handle atmospheric escape with this method, we do not do so for any of the test cases in this paper.

Equation (24) is solved within Argo in the same numerical manner as Nahoon (Wakelam et al. 2012) by the implicit time-dependent Gear method as incorporated by the Livermore Solver for Ordinary Differential Equations (DLSODE; Gear 1971; Brown & Hindmarsh 1989).

Once the fluid parcel has completed the atmospheric profile, the solar XUV actinic flux from 1 to 10000 Å as a function of depth, z [cm],⁷ and wavelength λ [Å] is calculated. We consider both the direct and approximate diffusive actinic flux. The local height-dependent actinic flux is calculated without any iteration on the local temperature. The cross sections for various photochemical reactions (Section 2.4) are multiplied by each vertical step ${({\rm{\Delta }}z)}_{j}$ [cm], where ${({\rm{\Delta }}z)}_{j}$ is the size of the step at height z_j. The total optical depth as a function of the wavelength takes the form

$$\begin{eqnarray}&&\tau (\lambda ,z)={{\rm{\Sigma }}}_{j}[{({\rm{\Delta }}z)}_{j}{{\rm{\Sigma }}}_{i}{\sigma }_{i}{n}_{{ij}}]+{\tau }_{s},\end{eqnarray} \tag{ 26 }$$

where i is summed over all species for which photoabsorption is considered (see Section 2.4 for a list of these species). ${\tau }_{s}$ is the optical depth due to Rayleigh scattering, and the actinic flux as a function of depth is defined as (Hu et al. 2012)

$$\begin{eqnarray}&&F(\lambda ,z)=F(\lambda ,{z}_{\mathrm{top}}){e}^{-\tau (\lambda ,z)/{\mu }_{0}}+{F}_{\mathrm{diff}},\end{eqnarray} \tag{ 27 }$$

where ${\mu }_{0}=\mathrm{cos}\theta$, where θ is the stellar zenith angle; we set ${\mu }_{0}=1/2$ for all calculations within this paper (see Hu et al. 2012, their Figure 7). ${F}_{\mathrm{diff}}$ denotes the actinic flux of the diffusive radiation, determined using the δ-Eddington two-stream method (Toon et al. 1989). Once the actinic flux is calculated, the photochemical rates are determined as in Section 2.4. Once the depth-dependent flux, $F(z,\lambda )$ [cm⁻² s⁻¹ Å⁻¹], is determined for all layers, the parcel's path through the atmospheric profile is repeated, now accounting for the photochemistry. The cosmic-ray ionization rate, $\zeta (z)$ [s⁻¹], is likewise calculated in a depth-dependent manner following Rimmer & Helling (2013) and incorporated into the chemistry (Section 2.4).

A new depth-dependent composition is constructed, then applied to Equation (26) to solve again for $F(z,\lambda )$. The value of $\zeta (z)$ does not change significantly between iterations. This process is repeated until the results converge; i.e., until the profile from the previous global calculation (transport + depth-dependent flux) agrees to within 1% the profile from the current global calculation. The number of repetitions depends on the parameters, but is typically between 5 and 12 global iterations. This iterative process is represented as a flow chart in Figure 5.

Zoom In Zoom Out Reset image size

This method is both simple and functional, requiring relatively little computational resources. It is also straightforward to adapt to diverse chemical environments, since it does not require the selection of "fast" and "slow" chemistry to ease computational speed. These strengths do not come without a cost: the simplistic dynamics does not transition as smoothly from the eddy diffusion regime to the molecular diffusion regime as the Eulerian formulation, and can result in steep changes over a handful of height steps.

In order to benchmark the Stand2015 chemical network in different planetary atmospheres, we test the molecular diffusion within Argo. We consider a 1D isothermal gas under a constant surface gravity, $g={10}^{3}\;\mathrm{cm}$ s⁻², with temperature T = 300 K, at hydrostatic equilibrium. The gas is initially composed of carbon and hydrogen atoms, each with a mixing ratio of ${X}_{0}({\rm{C}})=n({\rm{C}})/{n}_{\mathrm{gas}}=0.5$ and ${X}_{0}({\rm{O}})=0.5$ throughout. All chemistry is disabled. It is expected that the heavier species, carbon, will settle into the atmosphere, and the lighter species, hydrogen, will rise up, until they stratify. The analytic solution to this system is well known. The mixing ratio should be determined by the scale heights of the individual species such that, for the carbon abundance,

$$\begin{eqnarray}&&X({\rm{C}})=\displaystyle \frac{{X}_{0}({\rm{C}}){e}^{-z/{H}_{{\rm{C}}}}}{{X}_{0}({\rm{H}}){e}^{-z/{H}_{{\rm{H}}}}+{X}_{0}({\rm{C}}){e}^{-z/{H}_{{\rm{C}}}}},\end{eqnarray} \tag{ 28 }$$

where $X({\rm{C}})$ is the final steady state carbon mixing ratio, and ${H}_{{\rm{H}}}$ [cm] and ${H}_{{\rm{C}}}$ [cm] are the atmospheric scale heights for the hydrogen and carbon.

The code is run until steady state is achieved, when the carbon in the very upper atmosphere diffuses into the lower atmosphere. The steady state mixing ratio, as a function of height is compared the analytic mixing ratio, Equation (28), in Figure 6. The comparison is reasonable through the extent of the atmosphere.

Zoom In Zoom Out Reset image size

Below, we compare the results of our chemical kinetics to other results for HD 209458b and also for Jupiter and the Earth. Each of these objects has different boundary conditions and parameters. These conditions and parameters include the temperature profile of the object's atmosphere, the eddy diffusion profile, the elemental abundances, the initial composition at the lower boundary of the atmospheric profile, and the unattenuated UV flux. For HD 209458b, the conditions at the lower boundary of the atmospheric profile rapidly develop from the prescribed initial conditions toward chemical equilibrium. For Jupiter and the early Earth, the composition at the lower boundary is stable over the dynamical timescale (${{dn}}_{i}(z=0)/{dt}\approx 0$), and so the initial composition effectively acts as a lower boundary condition. The assumed elemental abundances and initial conditions at the lower boundary of the atmospheric profile are given in Table 1.

We take HD 209458b to have solar elemental abundances throughout its atmosphere, and set the initial conditions at the lower boundary of the atmosphere to be entirely atomic. The initial composition hardly matters here, since the composition quickly settles to chemical equilibrium at such a high temperature and pressure. The temperature profile and eddy diffusion profile for HD 209458b are both taken from Moses et al. (2011) so we can directly compare results.

Since HD 209458 is a G0 star, we use the solar UV flux. The unattenuated solar UV flux at 1 au is obtained from the SORCE data (Rottman et al. 2006) for 1–350 Å and 1150–10000 Å with data from PhIDRates for the 350–1150 Å range. The binned flux we use is plotted in Figure 7. This flux is adapted to HD 209458b by multiplying the solar UV flux by a factor of ${({d}_{\oplus }/{d}_{p})}^{2}$, where d_⊕ [au] is the distance from the Earth to the Sun and ${d}_{p}\approx 0.047\;\mathrm{au}$ is the approximate distance between HD 209458b and its host star. This may not be the most accurate approximation to the UV behavior of HD 209458, since it might have quite different activity from our Sun (Tu et al. 2015).

Zoom In Zoom Out Reset image size

For Jupiter, we use the temperature and eddy profiles from Moses et al. (2005). For consistency, we set the initial conditions at the lower boundary of Jupiter's atmosphere to be the same as Moses et al. (2005); see Table 1. The solar UV spectrum at 1 au is used for Jupiter, although multiplied by a factor of ${({d}_{\oplus }/{d}_{J})}^{-2}$, where ${d}_{J}\approx 4.5\;\mathrm{au}$ is the square of the distance between the sun and Jupiter.

For the present-day Earth, we use the measured surface mixing ratios from the US Standard Atmosphere 1976 (see Table 1) and the temperature profile from Hedin (1987, 1991), Figure 13. We use the present-day solar flux at 1 au as our incident UV flux.

We use the same chemical lower boundary conditions as from Kasting (1993) for the atmosphere of the early Earth (Table 1). The temperature profile for the early Earth is assumed to be the same as that of the present Earth (Hedin 1987, 1991), Figure 13. The UV field used for this model is that of the young Sun calculated using the scaling relationships of Ribas et al. (2005) for wavelengths between 1 and 1200 Å and the UV field of the solar analogue ${\kappa }^{1}$ Cet above 1200 Å (Ribas et al. 2010).

HD 209458b was first observed by Henry et al. (2000), and is one of a growing number of Hot Jupiters to have a measured spectrum, via transit (e.g., Queloz et al. 2000), and also in emission (e.g., Knutson et al. 2008). Various molecular species have been tentatively identified in the spectrum, such as TiO (Désert et al. 2008), water (Madhusudhan & Seager 2009; Swain et al. 2009; Beaulieu et al. 2010), CO, CO₂, and methane features (Madhusudhan & Seager 2009; Swain et al. 2009). HD 209458b has been extensively modeled with retrieval modeling (Madhusudhan & Seager 2009) and with hydrodynamic global circulation models (Showman et al. 2008). This planet has also been a popular target for non-equilibrium chemistry models such as those of Liang et al. (2003), Zahnle et al. (2009), Moses et al. (2011), Venot et al. (2012), Agúndez et al. (2014), and Lavvas et al. (2014).

We have chosen the atmosphere of HD 209458b as one candidate for benchmarking our results because it is well characterized and has been the subject of several non-equilibrium chemistry models, and it has a very high temperature even among Hot Jupiters. An additional benefit to HD 209458b is its suspected temperature inversion (Knutson et al. 2008, although this is debated; see also Schwarz et al. 2015), which allows us to test our chemistry at very high temperatures both at both high and low pressures. The thermal profile of HD 209458b from Moses et al. (2011) is shown in Figure 8. The local gas-phase temperature $T\gt 2000$ K both when $p\gt 100\;\mathrm{bar}$ and when the gas-phase pressure, $p\lt {10}^{-4}$ bar. This is a wide parameter space relevant for ion–neutral chemistry initiated via thermal ionization.

Zoom In Zoom Out Reset image size

We compare our results to the predictions of two different chemical kinetics models. (1) We compare our results to the results of Moses et al. (2011) with the ion–neutral chemistry disabled. (2) We compare the ionic abundances for our most abundant ions to the results of Lavvas et al. (2014). Also in this case, we disable cosmic-ray chemistry in order to draw a better comparison to the ion–neutral chemistry.

We compare our network and transport model to Moses et al. (2011) by examining the volume mixing ratios of major neutral species: H, H₂, He (hydrogen/helium chemistry), OH, H₂O, O, and O₂ (oxygen/water chemistry), N₂ and NH₃ (nitrogen chemistry), and CO, CH₄, and CO₂ (carbon chemistry). See Figure 9. These species were chosen because they are abundant and, in the case of H₂ and N₂, play an important role in the non-equilibrium chemistry. N₂ provides the reservoir for the transition between ${{\rm{N}}}_{2}\quad \rightleftharpoons$ NH₃. Other species were chosen because they contribute to features observed in transit spectroscopy (e.g., CO₂). The molecules CO and H₂O do both. Helium was chosen because its mixing ratio is not significantly affected by the chemistry. It changes with pressure due to molecular diffusion, and so it provides a useful comparison between our dynamical calculations and those of Moses et al. (2011).

Zoom In Zoom Out Reset image size

The transition of carbon between CO and CH₄, and nitrogen between N₂ and NH₃ is very sensitive to non-equilibrium chemistry, as ${\mathrm{CH}}_{4}\;\approx \;\mathrm{CO}$ when $p\sim 100\;\mathrm{bar}$ and $T\sim 2000$ K. As the pressure decreases rapidly while the temperature remains relatively high ($T\gt 1000$ K), the thermochemical equilibrium ratio for CH₄/CO plummets, approaching 10⁻⁷ at 0.1 bar in the HD 209458b atmosphere. The time it takes the carbon to meander from CH₄ to CO, however, becomes significantly longer than the relevant dynamical timescales (for HD 209458b, this timescale is prescribed by the eddy diffusion coefficient; see Bilger et al. 2013), and the CH₄ and CO abundances are quenched. The same sort of process governs the transition of nitrogen from N₂ to NH₃.

The pathways for both CH${}_{4}\quad \rightleftharpoons$ CO and ${{\rm{N}}}_{2}\quad \rightleftharpoons$ NH₃ interconversions are not well understood. In both cases, the paths competing with one another are often circuitous, and tend to be regulated by one of several reactions encountered along the journey, a slow rate-limiting step (Moses 2014). The timescale of the transition between species is almost entirely set by the rate by which that single reaction proceeds. As discussed in Section 2, rate coefficients can be frustratingly uncertain, with different estimations sometimes varying by more than an order of magnitude. For example, compare the rate experimental and theoretical rate constants for C₂H₆$\;\to$ CH₃ + CH₃ (Yang et al. 2009 and Kiefer et al. 2005, respectively). The path that one believes regulates these central transitions can be very different depending on what rate coefficients are used.

An illustrative example is the reaction ${\mathrm{CH}}_{3}+{{\rm{H}}}_{2}{\rm{O}}\to {\mathrm{CH}}_{3}\mathrm{OH}+{\rm{H}}$. Hidaka et al. (1989) has determined the rate for ${\mathrm{CH}}_{3}\;+\;{{\rm{H}}}_{2}{\rm{O}}\;\to \;\mathrm{Products}$, Reaction (29), proceeds with a barrier of $\approx 2670$ K (see Visscher et al. 2010 for a discussion on this reaction). With reasonable assumptions of the branching ratios for this reaction, namely that the branching ratios do not change much with temperature, one would set the same barrier to ${\mathrm{CH}}_{3}\;+\;{{\rm{H}}}_{2}{\rm{O}}\;\to \;{\mathrm{CH}}_{3}\mathrm{OH}\;+\;{\rm{H}}$, as done by Venot et al. (2012). However, Moses et al. (2011) carried out quantum chemical calculations for this reaction using MOLPRO and estimate a barrier for this particular branch of $\approx 10380$ K, much larger than the activation energies of the other branches. With the smaller barrier, the path carbon takes from CH₄ to CO proceeds as

$$\begin{eqnarray}{{\rm{H}}}_{2}+{\rm{M}} & \to & {\rm{H}}+{\rm{H}}+{\rm{M}}\\ {\mathrm{CH}}_{4}+{\rm{H}} & \to & {\mathrm{CH}}_{3}+{{\rm{H}}}_{2}\\ {\mathrm{CH}}_{3}\;+\;{{\rm{H}}}_{2}{\rm{O}} & \to & {\mathrm{CH}}_{3}\mathrm{OH}+{\rm{H}}\\ {\mathrm{CH}}_{3}\mathrm{OH}+{\rm{H}} & \to & {\mathrm{CH}}_{2}\mathrm{OH}+{{\rm{H}}}_{2}\\ {\mathrm{CH}}_{2}\mathrm{OH}+{\rm{M}} & \to & {{\rm{H}}}_{2}\mathrm{CO}+{{\rm{H}}}_{2}+{\rm{M}}\\ {{\rm{H}}}_{2}\mathrm{CO}+{\rm{H}} & \to & \mathrm{HCO}+{{\rm{H}}}_{2}\\ \mathrm{HCO}+{\rm{H}} & \to & \mathrm{CO}+{{\rm{H}}}_{2}\\ \mathrm{HCO}+{\rm{M}} & \to & \mathrm{CO}+{\rm{H}}+{\rm{M}}\end{eqnarray} \tag{ 29 }$$

$$\begin{eqnarray} &&{\mathrm{CH}}_{4}+{{\rm{H}}}_{2}{\rm{O}}\to \mathrm{CO}+3{{\rm{H}}}_{2}.\end{eqnarray} \tag{ 30 }$$

We adopt the rates of Moses et al. (2011) for this pathway, as well as the smaller rate coefficient for the three-body reaction H₂O + CH₂ + M $\to$ CH₃OH. An examination of our results would reveal that, as with Venot et al. (2012), the transition of carbon from CH₄ to CO is much more efficient than with Moses et al. (2011). We have examined the rates at which reactions proceed in our network and find another formation pathway:

$$\begin{eqnarray}{{\rm{H}}}_{2}+\mathrm{OH} & \leftrightarrow & {{\rm{H}}}_{2}{\rm{O}}+{\rm{H}}\\ \mathrm{OH}+{\rm{O}} & \to & {{\rm{O}}}_{2}+{\rm{H}}\\ {\mathrm{CH}}_{4}+{\rm{H}} & \to & {\mathrm{CH}}_{3}+{{\rm{H}}}_{2}\\ {\mathrm{CH}}_{3}+{\rm{H}} & \to & {\mathrm{CH}}_{2}+{{\rm{H}}}_{2}\\ {\mathrm{CH}}_{2}+{{\rm{O}}}_{2} & \to & \mathrm{COOH}+{\rm{H}}\\ \mathrm{COOH}\;+\;{{\rm{H}}}_{2}{\rm{O}} & \to & {\mathrm{CH}}_{2}{{\rm{O}}}_{2}+\mathrm{OH}\\ {\mathrm{CH}}_{2}{{\rm{O}}}_{2}+{\rm{M}} & \to & {\mathrm{CO}}_{2}+{{\rm{H}}}_{2}+{\rm{M}}\\ {\mathrm{CO}}_{2}+{\rm{H}} & \to & \mathrm{CO}+\mathrm{OH}\end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray} &&{\mathrm{CH}}_{4}+{\rm{O}}\to \mathrm{CO}+2{{\rm{H}}}_{2},\end{eqnarray} \tag{ 32 }$$

The atomic oxygen arises from thermal dissociation of OH or photodissociation of H₂O followed by diffusion downward. This pathway is critically dependent on Reaction (31). To our knowledge, the three-body rate coefficient for this reaction has not been determined. This reaction has instead appeared in our network as the reverse reaction of ${\mathrm{CH}}_{2}{{\rm{O}}}_{2}\;+\;\mathrm{OH}\;\to \;\mathrm{COOH}\;+\;{{\rm{H}}}_{2}{\rm{O}}$, for which we use an estimate based on reaction energetics (Mansergas & Anglada 2006). This pathway is highly uncertain, and removing it makes up the majority of the difference between our results and those of Moses et al. (2011) for methane between 1–10⁻⁴ bar. We suspect further differences owe to our different thermochemical constants and the use of slightly different solar abundances.

The path of nitrogen from NH₃ to N₂ is considerably more uncertain. The path is believed to roughly follow from NH₃ to NH via hydrogen abstraction, which will in turn react with another NH_X species to form N₂H_Y. This species will be destroyed either by reacting with hydrogen or via thermal decomposition, to form N₂. The reactions N₂H${}_{X+2}\;\to$ NH₂ + NH_X involve large uncertainties, which result in variations of the NH₃ quenched abundance by an order of magnitude. We find, similar to Moses et al. (2011), that

$$\begin{eqnarray}{{\rm{H}}}_{2}+{\rm{M}} & \to & {\rm{H}}+{\rm{H}}+{\rm{M}}\\ {\mathrm{NH}}_{3}+{\rm{H}} & \to & {\mathrm{NH}}_{2}+{{\rm{H}}}_{2}\\ {\mathrm{NH}}_{2}+{\rm{H}} & \to & \mathrm{NH}+{{\rm{H}}}_{2}\\ {\mathrm{NH}}_{2}+\mathrm{NH} & \to & {{\rm{N}}}_{2}{{\rm{H}}}_{2}+{\rm{H}}\\ {{\rm{N}}}_{2}{{\rm{H}}}_{2}+{\rm{H}} & \to & \mathrm{NNH}+{{\rm{H}}}_{2}\\ \mathrm{NNH}+{\rm{M}} & \to & {{\rm{N}}}_{2}+{\rm{H}}+{\rm{M}}\end{eqnarray} \tag{ 33 }$$

$$\begin{eqnarray} &&2{\mathrm{NH}}_{3}\to {{\rm{N}}}_{2}+3{{\rm{H}}}_{2}\end{eqnarray} \tag{ 34 }$$

with Reaction (33) as the rate-limiting step. The profile we have for NH₃ deviates considerably from the results of Moses et al. (2011), but this is in large part due to a difference in the nitrogen thermochemistry and initial abundances at high pressures propagating up through the atmosphere. Figure 9 shows that our quenching height is, in both cases, higher than for Moses et al. (2011), suggesting that the nitrogen in NH₃ migrates to N₂ more slowly in our network, even overtaking Moses et al. (2011) at $\sim {10}^{-4}$ bar, but that we start with less NH₃ than Moses et al. (2011). The increase in NH₃ abundance at $\sim 5\times {10}^{-6}$ bar is due to a formation path for NH₃ in Moses et al. (2011) that is less efficient in our network.

We conclude this section with a brief discussion of the most neutral ions, in comparison with Lavvas et al. (2014). We have plotted the most abundant ions in Figure 10. Note that for this paper, ${n}_{\mathrm{gas}}$ is a sum of all neutral gas particles, cations, ions, and electrons, so the mixing ratio of ions cannot increase above unity. This plot allows a direct comparison to Lavvas et al. (2014, their Figures 5 and 6). In our model, K⁺ is the most abundant ion deep within the atmosphere, followed by Mg⁺ and Fe⁺. Lavvas et al. (2014) does not consider these species, but they do not seem to have very much affect on the abundances of other ions deep within the atmosphere. When the pressure delves to 10⁻² bar, K⁺ deviates considerably between our results and those of Lavvas et al. (2014). This is likely due to the inclusion of several other ions in our model that become dominant charge carriers at this height, including several complex hydrocarbon ions, of the form ${{\rm{C}}}_{n}{{\rm{H}}}_{m}^{+}$. This indicates that ion–neutral chemistry can be significantly influenced by the variety of ions and neutral species under consideration. This will be especially true for the potassium chemistry. Our network contains a small number of potassium-bearing species. Including new species and reactions could significantly affect the degree of ionization. It will be interesting to discover how an expanded potassium and sodium chemistry affects the overall ion–neutral chemistry and the resulting abundances of trace species.

Zoom In Zoom Out Reset image size

Between 10⁻³ and 10⁻⁴ bar, Na⁺ overtakes K⁺ as the dominant positive charge carrier, and remains so until $\sim {10}^{-7}$ bar. This transition, the ratios between the ions, and the abundances of the ions are nearly identical between our model and that of Lavvas et al. (2014). Within the thermosphere of HD 209458b, there are some small discrepancies between our model and Lavvas et al. (2014) for He⁺, and quite large discrepancies for C⁺ which we suggest are owing to the non-Alkali photochemistry that Lavvas et al. (2014) include, but that we have not included here.

The atmosphere of Jupiter is divided into three regions: (1) the troposphere, where the gas-phase temperature T decreases with atmospheric height, (2) the stratosphere, where T is roughly constant with increasing height, and (3) the thermosphere, where T increases with height. In this section, we consider the chemical composition of Jupiter's stratosphere. The stratosphere of Jupiter is rich in hydrocarbons, owing to its large gas-phase C/O ratio, because the majority of the oxygen is locked in water ice and then gravitationally settles to below the tropopause. This is predicted to lead to a ${\rm{C}}/{\rm{O}}\sim 2\times {10}^{6}$ (Moses et al. 2005) in the absence of external sources of H₂O and CO₂ (Feuchtgruber et al. 1997; Moses et al. 2000a, 2000b), such as Shoemaker–Levy 9 (Cavalié et al. 2012). Jupiter's stratosphere provides an extreme example of how surface deposition can radically affect the C/O ratio, an effect more recently predicted for exoplanets and brown dwarfs (Bilger et al. 2013; Helling et al. 2014). The high C/O ratio, in combination with the large abundance of hydrogen (H₂ and CH₄ are the two most abundant volatiles in the stratosphere and lower thermosphere), means that the stratosphere of Jupiter is strongly reducing (Strobel 1983).

Fouchet et al. (2000) have observed ethane and acetylene in Jupiter's stratosphere. Ethylene has also been observed by Bézard et al. (2001). The stratospheric chemistry of Jupiter has been modeled by several groups, including Gladstone et al. (1996) and Moses et al. (2005). We adopt the lower boundary conditions and temperature profile that Moses et al. (2005) used and model the carbon–oxygen chemistry in the stratosphere of Jupiter, ignoring the nitrogen chemistry (most of the nitrogen will be locked in NH₃ ice). Boundary conditions are discussed in Section 4.1.

Our lower boundary is set to be identical to Moses et al. (2005). These boundary conditions are somewhat artificial; the carbon budget is controlled by the photochemistry and the dynamics. There is no effective destruction pathway for the stable hydrocarbons, but the timescale for their formation is often competing with the dynamical timescales. In the thermosphere, ∼10⁻⁷–10⁻⁸ bar, these hydrocarbons are lost through photodissociation and photoionization as well as molecular diffusion. At the base, the chemistry is halted once the dynamical timescale is reached, effectively treating the bottom boundary as an open boundary through which the hydrocarbons would continue to diffuse. In reality, the complex hydrocarbons are carried into Jupiter's deep atmosphere, where the high temperatures and pressures dissociate these hydrocarbons, and force the carbon budget to return to chemical equilibrium values: CH₄ with trace amounts of CO and other species. Visscher et al. (2010, their Figure 6) demonstrate how the carbon budget is set deep within Jupiter's atmosphere; we do not model this region.

With these reactions removed from the network, we ran the network using the temperature and K_zz profiles from Moses et al. (2005), shown in Figure 11. Comparisons between our results and a representative set of observations for the depth-dependent mixing ratios, for the species CH₄, C₂H₂, C₂H₄, C₂H₆, and C₄H₂, are shown in Figure 12. The observations for CH₄ are taken from Drossart et al. (1999) and Yelle et al. (1996), C₂H₂ observations are from Fouchet et al. (2000), Moses et al. (2005), and Kim et al. (2010), C₂H₄ observations are from Romani et al. (2008) and Moses et al. (2005), C₂H₆ observations are from Fouchet et al. (2000), Moses et al. (2005), Yelle et al. (2001) and Kim et al. (2010), and the C₄H₂ observations are from Fouchet et al. (2000) and Moses et al. (2005). We also incorporate observations for C₂H₂, C₂H₄, and C₂H₆ from Gladstone et al. (1996) and references therein.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Many of the published observations do not include error bars in atmospheric pressure. Additionally, there may seasonal in the pressure–temperature structure and the location of the homopause, which adds uncertainty to our predictions as a function of pressure. To account for these sources of uncertainty, we place error bars for the pressure at a factor of two above and below the published observations when errors in pressure were not given. These errors in pressure are of the same order as observations where errors in pressure are given. We do not compare our results for oxygen-bearing species, because the abundances of these species are expected to be greatly enhanced in the stratosphere by the addition of an external source of oxygen, such as Shoemaker–Levy 9.

The differences between our results and those of other models arise primarily because of different photochemistries and different rate constants, especially for the re-formation of methane after its photodissociation,

$$\begin{eqnarray}&&{\mathrm{CH}}_{3}+{{\rm{H}}}_{2}\to {\mathrm{CH}}_{4}+{\rm{H}},\;\mathrm{and}\end{eqnarray} \tag{ 35 }$$

$$\begin{eqnarray}&&{\mathrm{CH}}_{3}+{\rm{H}}+{\rm{M}}\to {\mathrm{CH}}_{4}+{\rm{M}}.\end{eqnarray} \tag{ 36 }$$

Differences between Jovian photochemical models can result in very large discrepancies between stratospheric abundances of complex hydrocarbons. The differences between Gladstone et al. (1996) and Moses et al. (2005) span several orders of magnitude in some cases (see Moses et al. 2005, their Figure 14).

Both ethane and acetylene agree reasonably well between our model and the observations, and the results for C₄H₂ lie more than a factor of five below the observational upper limits. Our predictions for the location of the methane homopause do not agree very well with observations. We use the eddy diffusion coefficient from Model C in Moses et al. (2005), and either this or the use of the Chapman–Enskog diffusion coefficient for Methane may be the source of the discrepancy. Our results are similar to the Model C results of Moses et al. (2005, their Figure 14). The molecule with the largest discrepancy between the two models is ethylene (C₂H₄), with the largest discrepancy between our predictions and the 1 millibar observations (ignoring the observation from Gladstone et al. 1996 that predicts a mixing ratio of $\sim {10}^{-8}$). In our model, the primary path of formation for ethylene follows from the photodissociation of ethane (Reaction 2679 in the network),

$$\begin{eqnarray}&&{{\rm{C}}}_{2}{{\rm{H}}}_{6}\;+\;\gamma \to {{\rm{C}}}_{2}{{\rm{H}}}_{4}+{{\rm{H}}}_{2},\end{eqnarray} \tag{ 37 }$$

and ethane is formed from CH₄ following paths to formation like this one:

$$\begin{eqnarray}2({\mathrm{CH}}_{4}\;+\;\gamma & \to & {}^{1}{\mathrm{CH}}_{2}\;+\;{{\rm{H}}}_{2}),\\ 2{(}^{1}{\mathrm{CH}}_{2}\;+\;{{\rm{H}}}_{2} & \to & {\mathrm{CH}}_{3}+{\rm{H}}),\\ {\mathrm{CH}}_{3}+{\mathrm{CH}}_{3}+{\rm{M}} & \to & {{\rm{C}}}_{2}{{\rm{H}}}_{6}+{\rm{M}};\end{eqnarray}$$

$$\begin{eqnarray} &&2{\mathrm{CH}}_{4}\;+\;2\gamma \to {{\rm{C}}}_{2}{{\rm{H}}}_{6}+2\;{\rm{H}}.\end{eqnarray} \tag{ 38 }$$

These differences may be resolved by a more careful accounting of pressure-dependent branching ratios, such as those of

$$\begin{eqnarray}&&{\rm{H}}+{{\rm{C}}}_{2}{{\rm{H}}}_{5}\to {\mathrm{CH}}_{3}+{\mathrm{CH}}_{3}\end{eqnarray} \tag{ 39 }$$

from Loison et al. (2015). We use the Kooij form for these reactions (Section 2.1), which does not account for the effect that pressure has on the rate constant.

Ion–neutral chemistry also makes a contribution, via the formation of C₂H₄ from the reaction

$$\begin{eqnarray}{\mathrm{CH}}_{5}^{+}\;+\;{{\rm{C}}}_{2}{{\rm{H}}}_{2} & \to & {{\rm{C}}}_{2}{{\rm{H}}}_{3}^{+}\;+\;{\mathrm{CH}}_{4}\\ {{\rm{C}}}_{2}{{\rm{H}}}_{3}^{+}+{e}^{-}+{\rm{M}} & \to & {{\rm{C}}}_{2}{{\rm{H}}}_{3}+{\rm{M}}\\ {{\rm{C}}}_{2}{{\rm{H}}}_{3}+{\mathrm{CH}}_{4} & \to & {{\rm{C}}}_{2}{{\rm{H}}}_{4}+{\mathrm{CH}}_{3},\end{eqnarray} \tag{ 40 }$$

and CH₅⁺ forms from a series of reactions starting with the photoionization of CH₃ and then a series of hydrogen abstractions, ${\mathrm{CH}}_{x}^{+}\;+\;{{\rm{H}}}_{2}\to {\mathrm{CH}}_{x+1}^{+}+{\rm{H}}$. It should be emphasized that this is not the primary formation pathway for ethylene, but it is an important path of formation in our chemistry and makes some contribution to the mixing ratios at 1 millibar.

Finally, there is a large discrepancy for CO, but this is not due to differences in the chemistry. Rather, this results from Moses et al. (2005) injecting CO, CO₂, and H₂O into Jupiter's stratosphere. The inclusion of this external source of oxygen-bearing species is justified by a number of data model comparisons mentioned at the beginning of this section. We neglected to include these external sources, and therefore oxygen-bearing species, especially H₂O and CO₂ (not shown), fail to agree with observations. Our results therefore suggest that some external source of oxygen-bearing species is necessary to explain the H₂O and CO₂ observations in Jupiter's stratosphere.

The Earth's atmosphere is well studied, and the profiles of trace species are well constrained, and the formation and destruction of these species is controlled by photochemistry and deposition. Comparing our results to the present-day Earth atmosphere therefore provides a comprehensive test of our chemical network (Section 4.4.1). Additionally, the connection between lightning-driven and NO_x chemistry⁸ has been extensively studied with experiments, observations, and models, and provides a useful regime in which to compare the results of Stand2015 applied to a lightning shock model (Section 4.4.2). It is important to find out what our model predicts in habitable environments before the onset of life, and so we apply our model to the early Earth (Section 4.4.3).

4.4.1. Present-day Earth Atmosphere

The best understood planetary atmosphere, in terms of both models and observations, is the atmosphere of the present-day Earth. Earth's atmosphere has been studied in situ, with the use of countless balloon experiments used to measure various trace elements, and remotely, with satellite measurements. Models of Earth's atmosphere range from simple to complex, both dynamically (1D diffusion to 3D global circulation models) and chemically (from treating only oxygen and hydrogen chemistry to modeling the transport and chemistry of chlorofluorocarbons and biological aerosols). Seinfeld & Pandis (2006) provide a useful introduction and review to the subject.

Zoom In Zoom Out Reset image size

Our interest is in validating our photochemical network to the present-day Earth, and not in coupling Earth's geochemistry to its atmospheric chemistry. We therefore make some simplifying assumptions when we set our boundary conditions. We compare our model to the contemporary Earth by setting the lower boundary conditions, temperature profile, and external UV field as given in Section 4.1. We present these comparisons for O₃, CH₄, and N₂O (Figure 14), NO and NO₂ (Figure 15), and OH and H₂O (Figure 16).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The data for O₃, CH₄, and N₂O is taken from the globally averaged mixing ratios from Massie & Hunten (1981). Following Hu et al. (2012), we apply error bars spanning an order of magnitude in mixing ratio to reflect the temporal and spatial variations. Our model fits the measured CH₄ to within the error bars throughout the atmosphere. The O₃ predicted by the model deviates from the data with errors at 15 km, and the N₂O deviates from the data with errors between 40 and 55 km. This may be due to an overestimation of the optical depth. If more UV photons in the model penetrated through to ∼10 km, the O₃ mixing ratios would be enhanced at 15 km, and the N₂O mixing ratios would be destroyed more efficiently deeper in the atmosphere.

The data for NO and NO₂ is taken from balloon observations at 35 deg N in 1993 (Sen et al. 1998), and here also we apply error bars spanning an order of magnitude to reflect spatial and temporal variations. As with Hu et al. (2012), we seem to overpredict the abundance of NO in the upper atmosphere (30–40 km). We find that this overprediction is due to Reaction 1300 in the network:

$$\begin{eqnarray}&&{{\rm{N}}}_{2}{\rm{O}}\;+\;{\rm{O}}{(}^{1}{\rm{D}})\to \mathrm{NO}+\mathrm{NO};\\ &&k=7.25\times {10}^{-11}\;{\mathrm{cm}}^{3}\;{{\rm{s}}}^{-1}.\end{eqnarray} \tag{ 41 }$$

We use the rate suggested by the JPL Chemical Kinetics and Photochemical Data for Use in Atmospheric Studies (Sander et al. 2011). If the rate constant for this reaction is decreased by a factor somewhere between 2 and 10, we come into much better agreement at 30–40 km, and worse agreement between 20–30 km (see Figure 15).

Finally, the data from OH and H₂O was taken from balloon measurements at various latitudes and heights in 2005 (Kovalenko et al. 2007). We plot each individual data point without error bars in order to represent the observed variations; changes at other points of the globe or at other times of the year or day may lead to more significant variations in the abundances. The H₂O predictions are within a factor of five of the observed water abundance, and our OH predictions lie within the measurements, indicating that the model correctly reproduces the water and OH mixing ratios.

4.4.2. Lightning Shock Model and NO_x Chemistry

It is also useful to the model's NO_x lightning-driven chemistry in the present-day atmosphere. For this purpose, we apply a simple shock model in order to explore the formation of NO_x species due to lightning at a single small region in the atmosphere. We employ the temperature and pressure calculations of Orville (1968, his Figures 1 and 3) and the timescaled results of Jebens et al. (1992, their Figures 2 and 3), fitting these to an exponential function. We use the following functions of temperature and pressure:

$$\begin{eqnarray}&&T(t)=300\;{\rm{K}}+(29800.0\;{\rm{K}})\;{e}^{-t/(55.56\mu {\rm{s}})};\end{eqnarray} \tag{ 42 }$$

$$\begin{eqnarray}&&P(t)=1.0\;\mathrm{bar}+(7.0\;\mathrm{bar})\;{e}^{-t/(5.88\mu {\rm{s}})}.\end{eqnarray} \tag{ 43 }$$

We start with present-day atmospheric chemistry at the base of the troposphere, except without the N₂O, NO, and NO₂ species, and with T = 300 K and p = 1 bar. The shock occurs at 1 ns, and is allowed to evolve until 0.1 ms. At this point the calculation is terminated, and another calculation initiated using for its initial conditions the final conditions of the shock model, except with temperature and pressure returned to 300 K and 1 bar, respectively. This model is run until 10⁴ s and results are shown in Figure 17.

Zoom In Zoom Out Reset image size

We find that the NO_x species are formed in our model thermally by the Zel'dovich mechanism (Zel'dovich & Raizer 1996):

$$\begin{eqnarray}&&{{\rm{O}}}_{2}+{\rm{M}}\rightleftharpoons {\rm{O}}+{\rm{O}}+{\rm{M}},\end{eqnarray} \tag{ 44 }$$

$$\begin{eqnarray}&&{{\rm{N}}}_{2}+{\rm{M}}\rightleftharpoons {\rm{N}}+{\rm{N}}+{\rm{M}},\end{eqnarray} \tag{ 45 }$$

$$\begin{eqnarray}&&{\rm{O}}+{{\rm{N}}}_{2}\to \mathrm{NO}+{\rm{N}},\end{eqnarray} \tag{ 46 }$$

$$\begin{eqnarray}&&{\rm{N}}+{{\rm{O}}}_{2}\to \mathrm{NO}+{\rm{O}},\end{eqnarray} \tag{ 47 }$$

$$\begin{eqnarray} &&{{\rm{O}}}_{2}+{{\rm{N}}}_{2}\to 2\mathrm{NO}.\end{eqnarray} \tag{ 48 }$$

We compare our NO yield to the lightning discharge experiments performed by Navarro-González et al. (2001). We use for our NO mixing ratio the values found before the end of the shock (10⁻⁴ s in Figure 17), between 10⁻² and 10⁻³, to (Navarro-González et al. 2001, their Equation (4)). We find that

$$\begin{eqnarray}P(\mathrm{NO}) & \approx & (2.4\times {10}^{22}\;{\rm{K}}/{\rm{J}})\displaystyle \frac{X(\mathrm{NO})}{{T}_{f}}\\ & \approx & 2-20\times {10}^{16}\;\mathrm{molecules}\;{{\rm{J}}}^{-1},\end{eqnarray} \tag{ 49 }$$

where T_f [K] is the "freeze-out" temperature after which the NO mixing ratio does not change appreciably over the timescale of the experiment, which we set to 1000 K (the approximate temperature of our model at $t\approx {10}^{-4}$ s). This is consistent with the production of NO in the "hot core" region of the experiment. This is also roughly consistent with the literature values for NO production of 10¹⁷ molecules J⁻¹ (Borucki et al. 1984; Price et al. 1997).

This is an order of magnitude comparison between the code and lightning experiments and models, and for a more complete comparison will need to be applied to a model atmosphere, where diffusion and photochemistry together will further process the NO_x species. We plan to do this in a future paper.

4.4.3. The Early Earth

The presence of life and the evolution of the Sun both have radically altered Earth's atmospheric chemistry. Oparin (1957) and Miller & Urey (1959) thought that the atmosphere of the early Earth ⁹ was largely reducing, dominated by methane, ammonia, and molecular hydrogen. Kasting (1993) made a compelling case that prebiotic formation of hydrogen would be too slow to allow for much molecular hydrogen in the atmosphere of the early Earth. Furthermore, a major constituent in the early Earth atmosphere needs to be a strong greenhouse gas, in order to compensate for the cooler young Sun. The atmospheric chemistry of the early Earth is difficult to determine, and a severe lack of data results in many possible early Earth chemistries. As an illustrative example, Tian et al. (2005) argue that hydrogen escape was less efficient during the first 1 Gyr as was previously thought.¹⁰ If Tian et al. (2005) are correct, then Earth's early atmospheric composition could have been reducing.

We present a model of the atmosphere of the early Earth, using the same lower boundary conditions as shown in Kasting (1993, his Figure 1), and a temperature profile for the present Earth (Hedin 1987, 1991)¹¹ , shown in Figure 13. The lower boundary conditions used for the early Earth are given in Section 4.1. We treat outgassing using the deposition method (Appendix C).

We compare our results to those of Kasting (1993, see his Figure 1). Our results are presented in Figure 18. The results compare reasonably well for CO and O₂, but not for H₂O and O. The CO abundance begins to increase at 30 km, 10 km higher than for Kasting (1993), and achieves a mixing ratio of $\approx 5\times {10}^{-3}$ at 60 km, which is within a factor of 2 of Kasting (1993). The O₂ likewise begins to rise above a mixing ratio of 10⁻⁶ 10 km higher in the atmosphere, and also achieves a mixing ratio of $\approx 2.5\times {10}^{-3}$, again within a factor of 2 of Kasting (1993). The water vapor profile is quite different, however. Instead of falling below a mixing ratio of 10⁻⁶ at 10 km, the H₂O mixing ratio in our model levels out at 5 × 10⁻⁴, increasing slightly at ∼50 km before plummeting. Also, the oxygen mixing ratio only reaches $\approx 3\times {10}^{-6}$, approximately two orders of magnitude below the mixing ratio predicted by Kasting (1993). These differences may be due to the different young solar UV fields assumed between ourselves and Kasting (1993), but we suspect that the differences are more likely due either to differences in the water condensation or the temperature profiles used. This seems especially likely for atomic oxygen, which is primarily destroyed by the reaction

$$\begin{eqnarray}&&{\rm{O}}+{{\rm{H}}}_{2}{\rm{O}}\to \mathrm{OH}+\mathrm{OH},\end{eqnarray} \tag{ 50 }$$

in spite of the sizeable 7640 K barrier. When the water vapor drops off at ∼55 km, this destruction route becomes unviable, and the atomic oxygen mixing ratio rapidly increases. 

Zoom In Zoom Out Reset image size