PHOTOCHEMISTRY IN TERRESTRIAL EXOPLANET ATMOSPHERES. I. PHOTOCHEMISTRY MODEL AND BENCHMARK CASES

The coupled one-dimensional continuity–transport equation that governs the chemical composition is

$$\begin{equation} \frac{\partial n}{\partial t}=P-nL-\frac{\partial \Phi }{\partial z}, \end{equation} \tag{ 1 }$$

where n is the number density of a certain species (cm⁻³), z is the altitude, P is the production rate of the species (cm⁻³ s⁻¹), L is the loss rate of the species (s⁻¹), and Φ is the upward vertical transport flux of the species (cm⁻² s⁻¹). The flux can be approximated by the eddy diffusion together with the molecular diffusion, viz.,

$$\begin{equation} \Phi =-KN\frac{\partial f}{\partial z} - DN\frac{\partial f}{\partial z} +Dn\bigg (\frac{1}{H_0}-\frac{1}{H}-\frac{\alpha _{\rm T}}{T}\frac{dT}{dz}\bigg), \end{equation} \tag{ 2 }$$

where K is the eddy diffusion coefficient (cm² s⁻¹), D is the molecular diffusion coefficient (cm² s⁻¹), N is the total number density of the atmosphere, f ≡ n/N is the mixing ratio of the species, H₀ is the mean scale height, H is the molecular scale height, T is the temperature (K), and α_T is the thermal diffusion factor (Banks & Kockarts 1973; Levine 1985; Krasnopolsky 1993; Yung & DeMore 1999; Bauer & Lammer 2004; Zahnle et al. 2006). The first term of Equation (2) represents eddy diffusion, and the last two terms represent molecular diffusion. The vertical transport flux is written in terms of the derivatives of the mixing ratio rather than in terms of the number density because this form is simpler and more straightforward to implement in numerical schemes. We show that Equation (2) is equivalent to the standard vertical diffusion equation in Appendix A.

The molecular diffusion coefficient D and the thermal diffusion factor α_T can be determined from the gas kinetic theory, but the eddy diffusion coefficient is a more speculative parameter and must be estimated empirically. For molecular diffusion, we use the following expressions in cm² s⁻¹ for H and H₂ in an N₂-based atmosphere

$$\begin{eqnarray} D({\rm H, N_2}) & = & \frac{4.87 \times 10^{17} \times T^{0.698}}{N},\nonumber \\ D({\rm H_2, N_2}) & = & \frac{2.80 \times 10^{17} \times T^{0.740}}{N} ; \end{eqnarray} \tag{ 3 }$$

in a CO₂-based atmosphere

$$\begin{eqnarray} {\rm D}({\rm H, CO_2}) & = & \frac{3.87 \times 10^{17} \times T^{0.750}}{N},\nonumber \\ {\rm D}({\rm H_2, CO_2}) & = & \frac{2.15 \times 10^{17} \times T^{0.750}}{N} ; \end{eqnarray} \tag{ 4 }$$

and in an H₂-based atmosphere

$$\begin{equation} {\rm D}({\rm H, H_2}) = \frac{8.16 \times 10^{17} \times T^{0.728}}{N}. \end{equation} \tag{ 5 }$$

In Equations (3)–(5), T has units of K and N has units of cm⁻³. The functional form of Equations (3)–(5) is derived from the gas kinetic theory and the coefficients are obtained by fitting to experimental data (Marrero & Mason 1972; Banks & Kockarts 1973). The expressions for molecular diffusion coefficients are valid for a wide range of temperatures, up to 2000 K, except for D(H₂, CO₂) that is only valid at low temperatures, i.e., T < 550 K (Marrero & Mason 1972). D(H, CO₂) is assumed to be 1.8 times larger than D(H₂, CO₂) (Zahnle et al. 2008). The thermal diffusion factor α_T is taken as a constant for H and H₂ as α_T = −0.38 (Banks & Kockarts 1973). The negative sign of α_T corresponds to the fact that the lightest molecules tend to diffuse to the warmest region.

The eddy diffusion coefficient is the major uncertainty in the one-dimensional photochemistry model. For Earth's atmosphere, the eddy diffusion coefficient can be derived from the number density profile of trace gases (e.g., Massie & Hunten 1981). The eddy diffusion coefficient of Earth's atmosphere is characterized by the convective troposphere and the non-convective stratosphere, which is not necessarily applicable for the exoplanets. The eddy diffusion coefficient may be typically parameterized as

$$\begin{eqnarray} && K=K_T \quad (z<z_T),\nonumber \\ && K={\rm min}\bigg (K_H, \quad K_T\bigg (\frac{N(z_T)}{N}\bigg)^{1/2}\bigg)\quad (z>z_T), \end{eqnarray} \tag{ 6 }$$

where K_T, K_H and z_T are independent parameters satisfying K_H > K_T. This formula is adapted from Yung & DeMore (1999) with the requirement of continuity in K. Our code can either import eddy diffusion coefficient from a file or specify eddy diffusion coefficient according to Equation (6).

The goal of the photochemistry model is to obtain a steady-state solution of each species, or a set of n(z), which makes the left-hand side of Equation (1) vanish and satisfies the boundary conditions. Assuming N_x species in the model and N_l equally stratified layers of the atmosphere, we transform the continuity Equation (1) into a discrete form as

$$\begin{equation} \frac{\partial n_i}{\partial t}=P_i-n_iL_i-\frac{\Phi _{i+1/2}-\Phi _{i-1/2}}{\Delta z}, \end{equation} \tag{ 7 }$$

where the subscript i denotes physical quantities in the ith layer, and the subscripts i + 1/2 and i − 1/2 mean that the flux is defined at the upper and lower boundary of each layer. It is physically correct to define the flux term at the boundary of each layer, which is also crucial for numerically preserving hydrostatic equilibrium for an atmospheric transport scheme that uses number density as independent variables. According to Equation (2),

$$\begin{eqnarray} \Phi _{i+1/2} & = & {-}(K_{i+1/2}+D_{i+1/2}) N_{i+1/2} \frac{f_{i+1}-f_i}{\Delta z} + D_{i+1/2}\frac{N_{i+1/2}}{2}\nonumber \\ & &\times\bigg [\frac{(m_a-m)g}{k_{\rm B}T_{i+1/2}}-\frac{\alpha _{\rm T}}{T_{i+1/2}}\frac{T_{i+1}-T_{i}}{\Delta z}\bigg ] (f_{i+1}+f_{i}), \end{eqnarray} \tag{ 8 }$$

where m_a is the mean molecular mass of the atmosphere, m is the molecular mass of the species, g is the gravitational acceleration, and k_B is the Bolzmann constant. We have approximated f_{i + 1/2} by (f_{i + 1} + f_i)/2 in Equation (2).

The combination of Equations (7) and (8) gives the following second-ordered centered discrete differential equation to be solved numerically,

$$\begin{eqnarray} \frac{\partial n_i}{\partial t} & = & P_i-n_iL_i + \bigg (k_{i+1/2}\frac{N_{i+1/2}}{N_{i+1}}-d_{i+1/2}\frac{N_{i+1/2}}{N_{i+1}}\bigg)n_{i+1} \nonumber \\ & & - \bigg (k_{i+1/2} \frac{N_{i+1/2}}{N_i} +d_{i+1/2}\frac{N_{i+1/2}}{N_i}\nonumber\\ && + k_{i-1/2} \frac{N_{i-1/2}}{N_i} - d_{i-1/2}\frac{N_{i-1/2}}{N_i}\bigg)n_i \nonumber \\ & & + \bigg (k_{i-1/2} \frac{N_{i-1/2}}{N_{i-1}}+d_{i-1/2} \frac{N_{i-1/2}}{N_{i-1}}\bigg)n_{i-1}, \end{eqnarray} \tag{ 9 }$$

where

$$\begin{eqnarray*} k_{i+1/2} & = & \frac{K_{i+1/2}+D_{i+1/2}}{\Delta z^2},\nonumber \\ d_{i+1/2} & = & \frac{D_{i+1/2}}{2\Delta z^2} \bigg [\frac{(m_a-m)g\Delta z}{k_{\rm B}T_{i+1/2}}-\frac{\alpha _{\rm T}}{T_{i+1/2}}(T_{i+1}-T_i)\bigg ]. \nonumber \end{eqnarray*}$$

In principle, ∂n/∂t = 0 is equivalent to a set of N_xN_l nonlinear algebraic equations. The set of nonlinear algebraic equations may be solved numerically by Newton–Raphson methods (Press et al. 1992). However, in practice, we find that Newton–Raphson methods require an initial guess to be in the vicinity of the solution. Instead, we treat the problem as a time-stepping problem by evolving the system according to Equation (9) to the steady state.

Implicit numerical methods are implemented to solve Equation (9). Due to the orders-of-magnitude differences in chemical loss timescales, the system is numerically stiff. We use the inverse-Euler method for the time stepping, as in most previous photochemistry models (e.g., Kasting et al. 1985; Nair et al. 1994; Zahnle et al. 2008). For each time step, we need to invert a matrix of dimension N_xN_l. As seen in Equation (9), the variation in any layer only depends on the number density in that layer and adjacent layers. As a result, the matrix is by nature block tridiagonal with a block dimension of N_x, which is solved efficiently by the Thomas algorithm (Conte & DeBoor 1972). The time step is self-adjusted in a way that the code updates the time step after each iteration according to the variation timescale of the whole chemical-transport system, i.e., the minimum variation timescale of each species at each altitude. As the system converges, larger and larger time steps are chosen.

The criteria of convergence to the steady-state solution is that the variation timescale of each species at each altitude is larger than the diffusion timescale of whole atmosphere, and the fluxes of gain and loss balance out for all species. The gain fluxes include surface emission, chemical production, and condensation; whereas the loss fluxes include chemical loss, dry deposition, wet deposition, atmospheric escape, and condensation. We require that the ratio between the net global flux and the column-integrated number density for all long-lived species is small compared to the diffusion timescale. This condition is very important in the determination of convergence. If any possible long-term trends of major species can be detected, we run the code for an extended time period to test the convergence rigorously.

The current photochemistry model can compute concentrations of 111 molecules or aerosols made of C, H, O, N, and S elements. These species are CO₂, H₂, O₂, N₂, H₂O, O, O(¹D), O₃, H, OH, HO₂, H₂O₂, N, NH, NH₂, NH₃, N₂H₂, N₂H₃, N₂H₄, N₂O, NO, NO₂, NO₃, N₂O₅, HNO, HNO₂, HNO₃, HNO₄, C, CO, CH₄, CH, CH₂, CH¹₂, CH₃, CH₂O, CHO, CH₃O, CH₃O₂, CHO₂, CH₂O₂, CH₄O, CH₄O₂, C₂, C₂H, C₂H₂, C₂H₃, C₂H₄, C₂H₅, C₂H₆, C₂HO, C₂H₂O, C₂H₃O, C₂H₄O, C₂H₅O, HCN, CN, CNO, HCNO, S, S₂, S₃, S₄, S₈, SO, SO₂, SO¹₂, SO³₂, SO₃, H₂S, HS, HSO, HSO₂, HSO₃, H₂SO₄, OCS, CS, CH₃S, CH₄S, CH₃NO₂, CH₃ONO₂, CH₅N, C₂H₂N, C₂H₅N, C₃H₂, C₃H₃, CH₃C₂H, CH₂CCH₂, C₃H₅, C₃H₆, C₃H₇, C₃H₈, C₄H, C₄H₂, C₄H₃, C₄H₄, C₄H₅, 1-C₄H₆, 1,2-C₄H₆, 1,3-C₄H₆, C₄H₈, C₄H₉, C₄H₁₀, C₆H, C₆H₂, C₆H₃, C₆H₆, C₈H₂, H₂SO₄ aerosols, S₈ aerosols, and organic hazes.

An important feature of the photochemistry code is the flexibility to choose a subset of molecules and aerosols to study a particular problem. Since the computation time for each time step scales with N³_x, decreasing the number of chemical species in the main loop greatly reduces the computation time. As is common practice in the effort to reduce the stiffness of the system and improve the numerical stability, "fast" species with relatively short chemical loss timescales are computed directly from the diagnostic equation, namely,

$$\begin{equation} n=\frac{P}{L}, \end{equation} \tag{ 10 }$$

which implies that photochemical equilibrium can be achieved within each time step. The choice of fast species varies with the atmospheric composition and should be considered on a case-by-case basis. The choices of "fast" species were usually hard-wired in previous photochemistry codes, which fundamentally limited their application to certain specific types of atmospheres (e.g., to Earth and to Mars). Our photochemistry code offers the flexibility to adjust the choice of species in photochemical equilibrium, which yields the capacity of treating atmospheres having very different oxidation states with the same code. We evaluate P and L in Equation (10) at each time step using the number density of all species from the previous time step and determine the number density of fast species according to Equation (10) after the current time step linearly. We neglect nonlinearity due to multiple species in photochemical equilibrium; as a result, we do not allow strongly inter-dependent species to be concurrently considered as in photochemical equilibrium, for example for S₂, S₃, and S₄. Our approach yields stable convergence to the steady-state solution from virtually any initial test solutions. We always verify the mass balance of the steady-state solution.

We compiled a comprehensive list of chemical and photochemical reactions from the literature. The production and loss rates in Equation (7) are provided from all chemical and photochemical reactions that produce or consume the relevant molecule. In the generic model, we included 645 bimolecular reactions, 85 ter-molecular reactions, and 93 thermal dissociation reactions. We included the thermal dissociation reactions to make the photochemistry model potentially adaptable to simulate hot planets. We mainly used the updated reaction rates from the online NIST database⁴ and refer to the JPL publication (Sander et al. 2011) for the recommended rate when multiple measurements are presented in the NIST database. For ter-molecular reactions, we used the complete formula suggested by the JPL publication (Sander et al. 2011) when the low-pressure-limiting rate and the high-pressure-limiting rate are available. For all reaction rates, we use the Arrhenius formula to account for the dependence on temperature when the activation energy data are available; otherwise, we adopted the value of experimental measurements, which are usually under ∼298 K. We performed careful comparisons of our reaction lists and corresponding reaction rates with those used by Nair et al. (1994), Pavlov & Kasting (2002), and Zahnle et al. (2008). Most of our reaction rates are the same as those used by other codes. For a dozen reactions, we find updated reaction rates listed in either the JPL publication or the NIST database updated after the publication of the cited photochemistry models. Also, for the general purpose of our photochemistry model, we include the reactions that are very slow at low temperatures (i.e., 200–400 K) but may be become important at higher temperatures >1000 K. In addition, a number of reactions that lack laboratory-measured rates may be important for low-temperature (i.e., 200–400 K) applications. We have adopted the rates for unmeasured sulfur reactions from Turco et al. (1982), Kasting (1990), and Moses et al. (2002), and the rates for reactions of C > 2 hydrocarbons from Yung & DeMore (1999). All chemical reactions and their reaction rates are tabulated in Table 1.

Chemical reaction rates may sensitively depend on temperature. The reported reaction rates are usually valid in certain specific temperature ranges. For example, reaction rates recommended by the JPL publication (Sander et al. 2011) are the reactions relevant to Earth's atmosphere, in general valid within the small range of 200–300 K. Hence it may be problematic to extend the rate expressions to high-temperature cases, such as hot planets. At the high-temperature end, the valid temperature range for some of the reaction rates measured in the combustion chemistry is a few thousand K. These reactions may be important for modeling hot Jupiters (e.g., Moses et al. 2011), but they should not be included in the computation of low-temperature cases, e.g., for Earth-like planets. We annotate the low-temperature reactions and the high-temperature reactions in Table 1. Moreover, most chemical reactions involving free radicals have very small activation energy, and thus weak temperature dependencies. Usually the reactions rate of free radicals are measured at room temperature and can be used in photochemistry models under different temperatures (Sander et al. 2011). Common atmospheric free radicals are: O(¹D), OH, H, N, NH₂, N₂H₃, C, CH, CH₂, CH¹₂, CH₃, CHO, CHO₂, CH₃O, CH₃O₂, C₂, C₂H, C₂H₃, C₂H₅, C₂HO, C₂H₃O, C₂H₅O, CN, CNO, S, S₂, SO₃, HS, HSO, HSO₂, HSO₃, CH₃S, and CS. For reactions with these molecules, if no activation energy is reported in the literature, we assume that the activation energy is negligible and their kinetic rates do not depend on temperature.

Molecules in the atmosphere absorb ultraviolet (UV) and visible light from the host star. If the absorbed photon carries enough energy, the molecule may be photodissociated to form free radicals. The photodissociation rate is proportional to the number density of photons with UV and visible wavelengths at each altitude. For direct stellar radiation, the optical depth τ includes the contribution from molecular absorption τ_a, Rayleigh scattering τ_r, and aerosol particle extinction τ_m. For the multiple-scattered (diffusive) radiation, we use the δ-Eddington 2-stream method implemented based on the formulation of Toon et al. (1989). The actinic flux of the diffusive radiation is F_diff = 2(F⁺ + F⁻) where F⁺ and F⁻ are the diffusive flux in the upward and downward direction. The photolysis flux at a certain altitude includes both the direct radiation and the diffusive radiation, i.e.,

$$\begin{equation} F(\lambda, z)=F_0(\lambda)\exp {[-\tau (\lambda, z)/\mu _0]} + F_{\rm diff}, \end{equation} \tag{ 11 }$$

where F₀ is the radiation flux at the top of the atmosphere where τ = 0, and μ₀ is the angle of the path of sunlight. By default we assume the zenith angle of the star to be 573 (see Appendix B for justification). The photodissociation rate is then

$$\begin{equation} J(z)=\frac{1}{2}\int q(\lambda)\sigma _a(\lambda)L(\lambda, z)d\lambda, \end{equation} \tag{ 12 }$$

where σ_a is the absorption cross section, q(λ) is the quantum yield that is defined as the ratio between the yield of certain photodissociation products and the number of photons absorbed, and L(λ, z) is the actinic flux with units of quanta cm⁻² s⁻¹ nm⁻¹.⁵ The 1/2 factor is included to account for diurnal variation (e.g., Zahnle et al. 2008), which is not used to model the dayside of a tidally lock planet. With the general-purpose model, we treat 70 photodissociation reactions. Different branches resulting from the photodissociation of one molecule are treated as different photodissociation reactions. The photodissociation reactions, the sources of data for cross sections and quantum yields, and the rates on the top of Earth's atmosphere are tabulated in Table 2.

For the UV and visible cross sections and the quantum yields, we use the recommended values from the JPL publication (Sander et al. 2011) when available. We also use the cross sections from the MPI-Mainz-UV-VIS Spectral Atlas of Gaseous Molecules⁶ when the JPL recommended values are not available or incomplete. There are a number of molecules of atmospheric importance that lack UV and visible cross sections and quantum yields, and we have estimated photolysis rates for them in the following way. The photodissociation timescale of S₂ has been measured to be ∼250 s at 1 AU from the Sun (De Almeida & Singh 1986), from which we estimate the photolysis rate of S₂. For other molecules that have no UV or visible-wavelength cross sections available in the literature, we assume their photolysis rate to be the same as another molecule that has similar structure. For example, the photolysis rate of HNO is assumed to be the same as HNO₂ (Zahnle et al. 2008); the photolysis rate of HSO is assumed to be the same as HO₂ (Pavlov & Kasting 2002); and the photolysis rate of N₂H₂ is assumed to be the same as N₂H₄.

Temperature dependencies of photolysis cross sections and quantum yields are considered. Notably, at 200 K, compared with room temperature, N₂O, N₂O₅, HNO₃, OCS, and CO₂ have smaller UV cross sections, leading to photolysis rates more than 10% lower; whereas NO₃ has a larger UV cross section, leading to a photolysis rate more than 10% higher. SO₂ has complex band structures in its UV spectrum and the cross sections depend on temperature as well. We take into account any other temperature dependencies reported for a gas at temperatures outside 290–300 K (see Table 2 for notes on temperature dependencies). In most cases, cross sections are also measured at lower temperature such as 200 K, primarily for Earth investigations (Sander et al. 2011). We use linear interpolation to simulate cross sections between 200 and 300 K, and do not extrapolate beyond this temperature range. It is worthwhile noting that UV cross sections of almost all gases at temperatures significantly higher than room temperature are unknown, and temperature dependencies of cross sections are not considered in previous high-temperature photochemistry models of hot exoplanets (e.g., Moses et al. 2011; Miller-Ricci Kempton et al. 2012). This may be a plausible simplification for gases without significant band structure in UV, but for gases such as H₂, CO₂, and SO₂, whose band structures are sensitive to temperature, extrapolation of cross sections to high temperature might induce significant errors. One needs to be cautious about the uncertainty of photolysis rate at temperatures much higher than room temperature.

Rayleigh scattering from atmosphere molecules introduces additional optical depth, particularly important for attenuation of UV radiation. The optical depth due to Rayleigh scattering is

$$\begin{equation} \tau _r(\lambda, z)=\int _z^{\infty }N(z^{\prime })\sigma _r(\lambda) dz^{\prime }, \end{equation} \tag{ 13 }$$

in which the Rayleigh scattering cross section is (e.g., Liou 2002)

$$\begin{equation} \sigma _r(\lambda) = C_r \frac{8\pi ^3(m_r(\lambda)^2-1)^2}{3\lambda ^4N_s^2}, \end{equation} \tag{ 14 }$$

where m_r is the real part of the refractive index of the molecule, C_r is a corrective factor to account for the anisotropy of the molecule, and N_s is the number density at the standard condition (1 atm, 273.15 K). The refractive index depends on the main constituent in the atmosphere. The refractive index of Earth's atmosphere is from Seinfeld & Pandis (2006), the refractive index of H₂ is from Dalgarno & Williams (1962), the refractive index of N₂ is from Cox (2000), the refractive index of CO₂ is from Old et al. (1971), and refractive indices of CO and CH₄ are given in Sneep & Ubachs (2005). In principle, the correction factor C_r depends on the molecule and the wavelength, but C_r is usually within a few percent with respect to the unity, except for C_r ∼ 1.14 for CO₂ (Sneep & Ubachs 2005). In the following we assume C_r = 1.061, the value for Earth's atmosphere at ∼200 nm (Liou 2002).

For the spectrum of a solar-type star (G2V), we use the Air Mass Zero (AM0) reference spectrum produced by the American Society for Testing and Materials.⁷ The AM0 spectrum covers a wavelength range from 119.5 nm to 10 μm. For the extreme-UV spectrum, we use the average quiet-Sun emission provided by Curdt et al. (2004).

Microphysical processes involved in the formation of atmospheric aerosols are nucleation, condensational growth, and coagulation (Toon & Farlow 1981; Seinfeld & Pandis 2006). It is beyond the purpose of this work to simulate these microphysical processes in detail. For photochemically produced aerosols, the competition between coagulation and sedimentation mainly determines the particle size distribution. A complete treatment of the atmospheric aerosols involves solving the steady-state size distribution function (e.g., Seinfeld & Pandis 2006).

In our photochemical model, we simplify the problem by assuming the particle radius to be a free parameter. In Earth's atmosphere, aerosols formed in the atmosphere often have a lognormal size distribution around 0.1–1 μm (Seinfeld & Pandis 2006). We define a particle radius parameter, r_p, to be the surface area average radius. This parameterization allows us to separate the complexity of aerosol formation from the photochemistry model, as well as to explore how the particle size affects the overall chemical composition.

We compute the production and loss rate of a molecule in the condensed phase (i.e., aerosols) based on its condensation timescale. When the molecule becomes supersaturated at a certain altitude, condensation can happen and aerosols form. The condensation/evaporation timescale is given by Hamill et al. (1977) and Toon & Farlow (1981) as

$$\begin{equation} \frac{1}{t_c}=\frac{m}{4\rho _{\rm p}}\bigg (\frac{8k_{\rm B}T}{\pi m}\bigg)^{1/2}\frac{n_{\rm g}-n_{\rm v}}{r_{\rm p}}, \end{equation} \tag{ 15 }$$

where t_c is the condensation/evaporation timescale, m is the mass of molecule, ρ_p is the particle density, k_B is the Boltzmann constant, T is the atmospheric temperature, n_g is the number density of the corresponding gas, and n_v is the saturated vapor number density at the corresponding pressure. The formula (15) is suitable for both condensation and evaporation processes. When n_g > n_v, the gas phase is saturated, so the condensation happens and t_c > 0. When n_g < n_v, the gas phase is unsaturated, so the evaporation is possible and t_c < 0. The production or loss rate of the molecule in the condensed phase is

$$\begin{equation} P=\frac{n_g}{t_c}, \quad \quad L=\frac{1}{t_c}. \end{equation} \tag{ 16 }$$

We include gravitational settling in the mass flux term of the continuity–transport equation (Equation (1)) for aerosol particles in addition to eddy mixing. The additional gravitational downward flux of the aerosol particle is

$$\begin{equation} \Phi _{\rm F} = -v_{\rm F} n_c, \end{equation} \tag{ 17 }$$

where v_F is the settling velocity of the particle in the atmosphere. The settling velocity is reached when the gravitational force is balanced by the gas drag. For aerosols with diameter of order of 1 μm, the settling velocity is reached within 10⁻⁵ s in Earth's atmosphere (Seinfeld & Pandis 2006). Therefore, we assume the falling velocity to be the settling velocity. The settling velocity can be derived from the Stokes' law (Seinfeld & Pandis 2006) as

$$\begin{equation} v_{\rm F} = \frac{2}{9}\frac{r_{\rm p}^2\rho _{\rm p}gC_{c}}{\mu }, \end{equation} \tag{ 18 }$$

where g is the gravitational acceleration, μ is the viscosity of the atmosphere, and C_c is the slip correction factor related to the mean free path (λ) of the atmosphere as

$$\begin{equation} C_c = 1+\frac{\lambda }{r_{\rm p}} \bigg [1.257+0.4\exp \bigg ({-}\frac{1.1r_{\rm p}}{\lambda }\bigg)\bigg ]. \end{equation} \tag{ 19 }$$

The treatment of photochemically produced aerosols described above is applicable to any molecules that could reach saturation as a result of photochemical production. In Earth's atmosphere, the photochemical aerosols include sulfate aerosols (H₂SO₄), sulfur aerosols (S₈), organic hazes, nitric acid aerosols (HNO₃), and hydrochloric acid aerosols (HCl). For now we have implemented sulfate aerosols (H₂SO₄) and sulfur aerosols (S₈). These aerosols and water vapor are the common condensable materials at habitable temperatures that commonly exist in planetary atmospheres (e.g., Kasting et al. 1989; Pavlov et al. 2000; Seinfeld & Pandis 2006). The required data for including an aerosol species in the photochemistry model is the saturation vapor pressure. Saturation vapor pressure of H₂SO₄ is taken as recommended by Seinfeld & Pandis (2006) for atmospheric modeling, with a validity range of 150–360 K. Saturation pressure of S₈ is taken as the total sulfur saturation pressure against liquid sulfur at T > 392 K and solid (monoclinic) sulfur at T < 392 K tabulated and expressed by Lyons (2008). In addition to aerosols, we use Equations (15) and (16) to compute the process of condensation of water vapor in the atmosphere. We do not consider evaporation of condensed water in the atmosphere because water droplets may grow by aggregation and rapidly precipitate out (Seinfeld & Pandis 2006). Saturation pressure of water is taken as that against ice at temperature lower than 273.16 K (Murphy & Koop 2005), and that against liquid water at temperature higher than 273.16 K (Seinfeld & Pandis 2006).

We treat the optical effect of aerosols by considering both scattering and absorption. Assuming a homogeneous sphere, the Mie theory (see Van de Hulst 1981 for a detail description) computes extinction cross section (σ_ext), single scattering albedo (w_s), and asymmetric factor (g_asym), based on the following parameters: the refractive index of the material (m_r + m_ii), r_p, and the wavelength. In this paper, we use the refractive index of S₈ aerosols from Tian et al. (2010) for the UV and visible wavelengths and from Sasson et al. (1985) for infrared (IR) wavelengths. We use the refractive index of H₂SO₄ aerosols (assumed to be the same as 75% sulfuric acid solution) from Palmer & Williams (1975) for UV to IR wavelengths, and Jones (1976) for far IR wavelengths.

The atmospheric chemical composition of a terrestrial exoplanet is ultimately determined by boundary conditions. The upper boundary conditions describe the atmospheric escape. The lower boundary conditions describe the surface emission, the deposition of molecules and aerosols to the surface, or the presence of a large surface reservoir of certain molecule (e.g., H₂O). The boundary conditions need to be properly provided to capture the physics of an exoplanet atmosphere.

The upper boundary conditions describe the atmospheric escape of an terrestrial exoplanet. The escape rates of exoplanet atmospheres are fairly uncertain depending on stellar soft X-ray and UV luminosity, exosphere chemistry, existence of magnitude fields, etc (e.g., Yelle et al. 2008; Tian 2009). We thus provide the following options of specifying escape rates in the photochemistry code:

• Type 1: $\Phi _{N_l+1/2}=0$, or no escape;
• Type 2: $\Phi _{N_l+1/2}=n_{N_l} V_{\rm lim}$, where V_lim is the diffusion-limited escape velocity;
• Type 3: $\Phi _{N_l+1/2}$ is an assigned nonzero value.

Here we use the same notation for flux as in Equations (7) and (8), such that the upper boundary condition replaces the flux at the upper boundary of the layer N_l in Equation (8). For the Type 2 upper boundary condition (atmospheric escape), the diffusion-limited velocity (V_lim) is

$$\begin{equation} V_{\rm lim}=D_{N_l+1/2}\bigg (\frac{1}{H_0}-\frac{1}{H}\bigg), \end{equation} \tag{ 20 }$$

where $D_{N_l+1/2}$ is the molecular diffusion coefficient evaluated at the top of the atmosphere. The diffusion-limited flux is the highest escape flux of an atmosphere in hydrostatic equilibrium (Hunten 1974). For Mars, the Jeans escape of H₂ reaches the diffusion-limited flux when the exobase temperature is above 400 K (Zahnle et al. 2008). We use the Type 2 upper boundary condition for the escape of H and H₂ in our model, and generally we assume no escape for all other species, i.e., the Type 1 upper boundary condition. The Type 3 upper boundary condition may be used when processes above the neutral atmosphere are important. For example, an influx of atomic N can represent the photodissociation of N₂ in the upper atmosphere. The Type 3 boundary condition may also be used when hydrodynamic escape has to be considered.

The lower boundary conditions describe the interaction between the atmosphere and the surface, which includes surface emission and surface deposition. The three types of lower boundary conditions are:

• Type 1: n₁ is assigned;
• Type 2: Φ_1/2 = −n₁V_DEP where V_DEP is molecule-specific dry deposition velocity;
• Type 3: Φ_1/2 is assigned.

Again we use the same notation for flux as in Equations (7) and (8), and Φ_1/2 replaces the flux at the lower boundary of the layer 1 in Equation (8).

The Type 1 lower boundary condition presents a large reservoir at the surface and for this paper we use this condition for water vapor to simulate the effect of a surface with oceans. This approach is equivalent to setting the relative humidity at the surface to be a constant. Note that specifying the Type 1 lower boundary condition means decreasing the number of free variables; n₁ is fixed as the lower boundary condition and is no longer considered as a variable in the main computation loop.

The Type 2 lower boundary condition specifies the dry deposition velocity, which is a key parameter that determines the chemical composition of the atmosphere and is a major unknown. The deposition velocity depends on both the dynamical properties of the lower atmosphere and the chemistry of the planet's surface. With the number density at the bottom layer (n₁) computed from the photochemistry model, the interaction between the bottom layer of atmosphere and the surface consists of two steps: first, the molecular transport across a thin stagnant layer of air adjacent to the surface, called the quasi-laminar sublayer; second, the uptake at the surface (Seinfeld & Pandis 2006). A parameterization of dry deposition velocity involving these two steps is described in Appendix C.

Physically the dry deposition requires a sink at the surface; for a gas without effective surface sink, the dry deposition velocity should be zero. In particular, the surface deposition of a number of gases, including CO, H₂, CH₄, and NH₃, is primarily removed by microorganisms on Earth (e.g., Kharecha et al. 2005; Seinfeld & Pandis 2006; J. F. Kasting 2012, private communication), and the canonical values for their dry deposition velocities are not applicable for presumably abiotic planets. Take carbon monoxide for an example. If a planet has no ocean, there is no known reaction that can consume CO at the surface, and therefore the dry deposition velocity of CO for an abiotic desiccated planet should be zero (like Mars). If the surface has an ocean, the dissolved CO may be naturally and slowly converted into acetates by OH in sea water, or biologically converted into acetates at a much faster rate. The CO deposition velocity has been estimated to be 1.2 × 10⁴ cm s⁻¹ for the most efficient biological removal and 10⁻⁸–10⁻⁹ cm s⁻¹ on an abiotic ocean planet (Kasting 1990; Kharecha et al. 2005). For another example, it is common to use a fairly large dry deposition velocity for SO₂ (∼1 cm s⁻¹) to study the sulfur cycle in Earth's marine atmosphere (e.g., Toon et al. 1987). In contrast, the dry deposition velocity is assumed to be zero, or reduced by an arbitrary factor of up to 1000, to mimic a putative Mars ocean that is believed to be saturated with dissolved SO₂ and other sulfur species (e.g., Halevy et al. 2007; Tian et al. 2010). In general, the dry deposition velocity depends on a broad context of planetary geochemistry, notably on the surface mineralogy, the acidity of ocean, and the surface pressure and temperature. Therefore, it is critical for exoplanet exploration to understand the interaction between atmospheric chemistry composition of terrestrial exoplanets and the surface deposition.

In addition to the dry deposition at the surface, we also include the wet deposition throughout the atmosphere as a removal process for soluble species. We use the parameterization of Giorgi & Chameides (1985) as

$$\begin{equation} k_R(z)=f_R\times \frac{n_{\rm H_2O}(z) k_{\rm H_2O}(z)}{55A_V[L\times 10^{-9}+(H^{\prime }RT(z))^{-1}]}, \end{equation} \tag{ 21 }$$

where f_R is a reduction factor adjustable in the model, $k_{\rm H_2O}$ is the precipitation rate taken to be 2 × 10⁻⁶ s⁻¹, A_V is the Avogadro's number, L is the liquid water content taken to be 1 g m⁻³ in the convective layer near the surface, and H' is the effective Henry's Law constant measuring the solubility of the molecule in the unit of mol dm⁻³ atm⁻¹. The effective Henry's law constant may differ from the standard Henry's law constant when taking into account dissociation in the aqueous phase (Giorgi & Chameides 1985; Seinfeld & Pandis 2006). In the model we use the effective Henry's Law constant published in Giorgi & Chameides (1985) as well as the standard Henry's law constants from the NIST Chemistry Webbook⁸. Since the parameterization of Equation (21) is primarily for modeling Earth's atmosphere, the specific reduction factor f_R should be applied when it is reasonable to believe the hydrological cycle is reduced on an exoplanet (e.g., Tian et al. 2010).

The Type 3 lower boundary condition has assigned flux that represents the surface emission. Note that the dry deposition (Type 2) and assigned flux (Type 3) lower boundary conditions can be used at the same time, but assigned number density (Type 1) boundary condition overrules other lower boundary conditions. For example, SO₂ may be deposited to the surface at a rate proportional to the number density of the bottom layer, and also be emitted from the surface to the atmosphere at an assigned flux.

We finish this section with a definition of the so-called "redox balance," the requirement of which arises when the fixed mixing ratio lower boundary condition (Type 1) is used in the photochemistry model. Following Kasting & Brown (1998), Zahnle et al. (2006), and Segura et al. (2007), we define H₂O, N₂, CO₂, and SO₂ as redox neutral, and assign the redox number (${\cal R}$) of any H–O–C–N–S molecule as the number of hydrogen in excess, i.e.,

$$\begin{equation} {\cal R}({\rm H_{a_1}O_{a_2}N_{a_3}C_{a_4}S_{a_5}}) = a_1 - 2a_2 + 4a_4 + 4a_5. \end{equation} \tag{ 22 }$$

Note that our definition of the redox number differs from that of Zahnle et al. (2006) and Segura et al. (2007) in the way that we count the number of H and they count the number of H₂. With our definition, for example, the redox number of H₂S is 6, and the redox number of H₂SO₄ is −2. The redox balance says that the total redox influx to the atmosphere (i.e., surface emission) should be balanced by the total redox outflux from the atmosphere (i.e., atmospheric escape, dry deposition, and wet deposition), otherwise the atmosphere is being oxidized or reduced. The redox balance is equivalent to the conservation of total budget of hydrogen in the atmosphere, and is equivalent to the conservation of the total number of electrons in the atmosphere. One might consider the redox balance to be redundant. The redox balance is redundant if the lower and upper boundary conditions for all species are specified in fluxes (and not in mixing ratios). In many applications, however, it is useful to use the Type 1 lower boundary condition, which specifies mixing ratio. In that case, the mass conservation would not necessarily guarantee the redox balance. After all, an imbalance in the redox budget for a steady-state solution indicates either bugs in chemistry kinetics or unphysical boundary conditions. We explicitly check the redox balance for all of our model outputs.

As a test of transport-related schemes of the photochemical model, we first consider a transport-only case. In theory, for species whose removal timescales are significantly larger than their transport timescales, such as CO₂ in Earth's atmosphere, their mixing ratios do not change with altitude, i.e., they are well mixed.

We turn off the chemical network deliberately and compute a model using Earth's temperature profile to test the eddy diffusion transport, molecular diffusion transport, condensation, and rainout schemes of our chemical-transport model. A valid transport model should preserve hydrostatic equilibrium, predict well mixed mixing ratios for long-lived species, and predict a mixing ratio gradient for species that is rapidly removed in the atmosphere. We have verified that the transport scheme, although written in terms of number density rather than mixing ratio, preserves hydrostatic equilibrium. The fact that the code maintains hydrostatic equilibrium indicates that the transport scheme is numerically correct. We verify that long-lived species are well mixed, such as CO₂ throughout the atmosphere, as shown in Figure 1. We assign a mixing ratio for CO₂ at the surface and the code is able to predict a well-mixed vertical profile. We also verify that for H₂ molecular diffusion tends to increase the mixing ratio when the molecular diffusion coefficient is comparable with the eddy diffusion coefficient (i.e., the homopause; see Figure 1). With a diffusion-limited escape, the effect of molecular diffusion and the effect of escape on the H₂ vertical profile near the homopause largely cancel out (Figure 1), consistent with the definition of diffusion-limited escape. We also show the behavior of H₂S, a poorly mixed species, to compare with well-mixed CO₂. In this simplified model the only way to remove H₂S is photolysis, which requires photons with wavelengths shorter than 260 nm. Most of the photons in this wavelength range are effectively shielded by the major species O₂, so that removal rate of H₂S by photolysis near the surface is smaller than that in the stratosphere (see Figure 1).

Zoom In Zoom Out Reset image size

Water vapor content in a rocky exoplanet's atmosphere is controlled by the water reservoir at the planetary surface, vertical transport, and condensation. In our model, the mixing ratio of water vapor at the surface is assigned according to appropriate temperature-dependent relative humidity. Water vapor can be transported up into the atmosphere by eddy diffusion, and as the temperature decreases the atmosphere may become supersaturated in water as the altitude increases. For an Earth-like planet atmosphere, it is appropriate to assume a water vapor mixing ratio of 0.01 at the surface (which corresponds to a relative humidity of about 60%). The condensation scheme becomes effective when water vapor is supersaturated, which keeps the water vapor profile along the saturation profile. Such an approach to computing the atmosphere water vapor content is commonly adopted by other previously described photochemistry models of terrestrial planets (e.g., Nair et al. 1994; Yung & Demore 1999; Zahnle et al. 2008).

The one-dimensional transport model tends to saturate the tropopause and overpredict the amount of water vapor in the stratosphere. In fact, when using the US Standard Atmosphere 1976 as the temperature profile, the tropopause temperature is 217 K, much warmer than the required "cold trap" temperature (∼200 K) in order to maintain a dry stratosphere (water vapor mixing ratio in a few ppm). Common one-dimensional transport-condensation schemes may overpredict the amount of water vapor above the cold trap by a few orders of magnitude, and therefore overpredict the strength of HO_x cycle in the stratosphere. This is a well-known problem in modeling Earth's atmosphere, mainly due to spatial and temporal variability of tropopause temperatures (K. Emmanuel 2012, private communication). Tropospheric convection is most effective in transporting water vapor in tropics, where the tropopause temperature is the lowest and the number density of water vapor at the tropopause is the lowest compared with other latitudes. It is therefore likely that the tropopause is highly unsaturated as a global average. We use the temperature profile of the equatorial region in January (see Figure 1), and limit the water vapor saturation ratio to within 20%. In this way we reproduce the water vapor profile of the US Standard Atmosphere 1976, which has a dry stratosphere (Figure 1). The classic Manabe–Wetherald relative humidity profile for one-dimensional photochemistry models of Earth's atmosphere also has 20% relative humidity at the tropopause (Manabe & Wetherald 1967).

As a result of water vapor condensation and rain, soluble gases are removed from the atmosphere by wet deposition (rainout). We expect to see a decreasing slope in mixing ratio of the species being rained out, which depends on the rainout rate, itself dependent on the solubility of the species. For example, SO₂ is much more soluble than H₂S, so that SO₂ is more effectively rained out in the troposphere and cannot accumulate there (see Figure 1).

We compute a photochemical model to simulate the present-day Earth, and compare our results with globally averaged measurements. The most important photochemical process in Earth's atmosphere is formation of the ozone layer, which leads to the temperature inversion in the stratosphere (e.g., Seinfeld & Pandis 2006). On the other hand, the tropospheric chemistry is mostly controlled by the hydrological cycle and surface emission of CH₄, NO_x (i.e., NO and NO₂), and SO₂, and involves coupled processes of photochemistry and deposition. The reproduction of current Earth atmospheric composition is therefore a comprehensive validation of various aspects of the photochemistry model.

We start with a nominal 80% N₂ and 20% O₂ composition, having well-mixed CO₂ of 350 ppm. We assume the water vapor mixing ratio at the surface is 0.01 (i.e., Type 1 lower boundary condition), and let the water vapor be transported in the atmosphere and condense near the tropopause. The temperature profiles are assumed to be the 1976 US Standard Atmosphere (subsequently referred to as "Model A") or a reference tropical atmosphere in January of 1986 from COSPAR⁹ International Reference Atmosphere (CIRA, subsequently referred to as "Model B;"), and the eddy diffusion coefficient is adopted from Massie & Hunten (1981), who have derived it from vertical profiles of trace gases. We model the atmosphere from 0 to 86 km altitude as 43 equally spaced layers. The incoming solar radiation is cut off for wavelengths shorter than 100 nm to account for the thermosphere absorption. The photolysis rates at the top layer of the atmosphere are tabulated in Table 2. Key trace species are assigned emission rates from the surface according to typical global values (Seinfeld & Pandis 2006). Emission rates and dry deposition velocities of these species are tabulated in Table 3.

Our models correctly simulate the photochemistry in Earth's stratosphere. First, the formation of the ozone layer is correctly predicted. The modeled vertical profiles and amounts of ozone are consistent with the globally averaged measurements, as shown in Figure 2. Second, the nitrogen oxide cycle and the HO_x cycle are properly simulated. In Earth's stratosphere, chemical species are rapidly converted to each other within the NO_x group and the HO_x group, and both cycles lead to catalytic destruction of ozone (e.g., Seinfeld & Pandis 2006). We reproduce the vertical profiles of NO_x and HO_x in the stratosphere (see Figure 2), qualitatively consistent with mid-latitude measurements. It appears that the models tend to overpredict the amounts of NO in the lower stratosphere. Third, for relatively long lived species in the stratospheres like N₂O and CH₄, our models produce correct altitude gradient of mixing ratios.

Zoom In Zoom Out Reset image size

The key element of the tropospheric chemistry is the production of the hydroxyl radical OH, because OH is the major removal pathway for most species emitted from the surface. It has been established that OH in Earth's troposphere is produced by reactions between O(¹D) and H₂O, and O(¹D) is produced by photolysis of O₃ (Seinfeld & Pandis 2006). In turn, the main source of tropospheric ozone is the NO_x cycle, with non-neglible effects of HO_x cycle and hydrocarbon chemistry as well. We find the surface ozone mixing ratio is 19 ppb and 22 ppb for Model 1976 and Model 1986 respectively, and the surface OH number density is 2.0 × 10⁶ cm⁻³ for Model 1976 and 2.2 × 10⁶ cm⁻³ for Model 1986, very close to the commonly adopted OH number density of 1.0 × 10⁶ cm⁻³. We also confirm that the steady-state mixing ratios of major trace gases near the surface are consistent with ground measurements (see Table 3). The main removal mechanism of CO, CH₄, and NH₃ is through reactions of OH, and N₂O is long-lived in the troposphere and is transported up to the stratosphere. Considering the significant spatial and temporal variability of the amount of OH, we conclude that our photochemistry model correctly computes the chemistry in Earth's troposphere.

Our photochemistry model correctly treats the sulfur chemistry in Earth's atmosphere (see Table 3). All sulfur-bearing emission, if not deposited, is oxidized in multiple steps in the troposphere and eventually converted into sulfate. Sulfite (S^{4 +}) and sulfate (S^{6 +}) are soluble in water and effectively removed from the atmosphere by rainout. Our photochemistry models simulate these processes, and find steady-state mixing ratios of sulfur-bearing species (e.g., H₂S and SO₂) consistent with the ground measurements. The models also predict the saturation of sulfate and the formation of sulfate aerosols as expected. The modeled mixing ratio of sulfate is slightly larger than observations, which might be related to the fact that we do not consider the hydration of H₂SO₄, nor the sulfate-facilitated cloud formation in our models.

In summary, we validate our photochemistry model by reproducing the chemical composition of current Earth's atmosphere (see Table 3 and Figure 2). The photochemistry model is successful in predicting the formation of the ozone layer, treating key chemical cycles in both the stratosphere and the troposphere, computing oxidation of hydrocarbon, ammonia, and sulfur-bearing species in the troposphere, and transporting long-lived species from the troposphere to the stratosphere. These aspects involve all chemical kinetics, photolysis, and transport processes, which not only verify that our photochemistry model is suitable for applications in oxidizing atmospheres, but also allows our model to be applied to other atmospheric scenarios.

We validate our photochemistry model by simulating the atmosphere of present-day Mars. The current atmosphere of Mars is a thin CO₂-based atmosphere and its bulk chemical composition is known from extensive ground-based observations (see Krasnopolsky 2006 and references therein). We validate our photochemistry code by reproducing the observed mixing ratios of CO and O₃, and comparing with previous results regarding Martian atmosphere photochemistry. Key parameters and results of the present-Mars atmosphere model are shown in Table 4 and Figure 3. We emphasize that it is important to use temperature-dependent UV cross sections for modeling the photolysis and UV penetration in the Martian environment (see Anbar et al. 1993 for an error analysis of temperature dependence of CO₂ in the Martian atmosphere).

Zoom In Zoom Out Reset image size

The most important aspect of Martian atmospheric photochemistry is the stabilization of the CO₂ atmosphere by H, OH, and HO₂ (commonly referred to as odd hydrogens; Nair et al. 1994; Krasnopolsky 2006; Zahnle et al. 2008). Odd hydrogens, as catalysts, facilitate the recombination of photochemically produced CO and O₂ and maintain the CO₂ dominance of the Mars' atmosphere (e.g., Nair et al. 1994; Krasnopolsky 2006; Zahnle et al. 2008). A pure CO₂ atmosphere is unstable against photolysis, because CO₂ can only be restored with a slow three-body reaction CO + O + M $\longrightarrow$ CO₂ + M (e.g., Yung & DeMore 1999). Odd hydrogen species, including H, OH, and HO₂, are produced by photolysis of water vapor; and trace amounts of odd hydrogens can effectively stabilize the CO₂ dominant atmosphere (e.g., Nair et al. 1994; Krasnopolsky 2006; Zahnle et al. 2008). Nonetheless, one-dimensional photochemistry models tend to overaddress the problem, predicting amounts of CO, O₂ and O₃ several times smaller than the observed values (e.g., Nair et al. 1994; Krasnopolsky 2006). By considering a slow dry deposition of H₂O₂ and O₃ to the surface, Zahnle et al. (2008) are able to predict the amount of CO, O₂, and H₂ that match observations. We follow the assumptions of Zahnle et al. (2008), and confirm the finding of appropriate amounts of O₂, CO and H₂ (see Table 4 and Figure 4). We therefore reproduce the photochemical stability of the current Martian atmosphere.

The timescale of CH₄ removal can also be used for model validation. Rapid variation (over several years) of the amount of CH₄ in Mars' atmosphere has been reported (e.g., Mumma et al. 2009), but the modeling of coupled general circulation and gas-phase chemistry find no known gas-phase chemistry path that allow such a rapid removal (e.g., Lefèvre & Forget 2009). Based on our fiducial model of Mars' atmosphere, the loss timescale for CH₄ is computed to be about 240 years, within the same order of the magnitude of Lefèvre & Forget (2009).

In summary, we validate our photochemistry model by reproducing the atmospheric compositions of current Earth and Mars. We find that the model gives consistent results compared to observations and previous photochemistry models. All physical and chemical processes including photolysis, chemical reactions, transports, condensation, and deposition are rigorously tested in these examples.

We here list our main findings on general chemistry properties of atmospheres on habitable terrestrial exoplanets.

Atomic hydrogen (H) is a more abundant reactive radical than hydroxyl radical (OH) in anoxic atmospheres. Atomic hydrogen is mainly produced by water vapor photodissociation; and in anoxic atmospheres the main ways to remove atomic hydrogen is its recombination and reaction with CO. This is in contrast to oxygen-rich atmospheres (e.g., current Earth's atmosphere) in which H is quickly consumed by O₂. As a result, removal of a gas by H is likely to be an important removal path for trace gases in an anoxic atmosphere. Atomic oxygen is the most abundant reactive radical in CO₂-dominated atmospheres. Due to the photochemical origin of the reactive species including H, OH, and O, their abundances in the atmosphere around a quiet M dwarf are two orders of magnitude lower than their abundances around a Sun-like star.

Dry deposition velocities of long-lived compounds, notably major volcanic carbon compounds including methane, carbon monoxide, and carbon dioxide, have significant effects on the atmospheric oxidation states. The specific choice of dry deposition velocities for emitted gases and their major photochemical byproducts in the atmosphere is critical to determining the atmospheric composition on terrestrial exoplanets.

Volcanic carbon compounds (i.e., CH₄ and CO₂) are chemically long-lived and tend to be well mixed in terrestrial exoplanet atmospheres, whereas volcanic sulfur compounds (i.e., H₂S and SO₂) are short-lived. CH₄ and CO₂ have chemical lifetime longer than 10,000 years in all three benchmark atmospheres ranging from reducing to oxidizing, implying that a relatively small volcanic input can result in a high steady-state mixing ratio. The chemical lifetime CO, another possible volcanic carbon compound, ranges from 0.1 to 700 years depending on the OH abundance in the atmosphere.

We find abiotic O₂ and O₃, photochemically produced from CO₂ photolysis, build up in the 1 bar CO₂-dominated atmosphere if volcanic emission rates of reducing gases (i.e., H₂ and CH₄) are more than one order of magnitude lower than current Earth's volcanic rates. Abiotic O₂ can be a false positive for detecting oxygenic photosynthesis, but the combination of O₂/O₃ and reducing gases remains a rigorous biosignature.

4.2.1. H₂-dominated Atmospheres

The main reactive agent in the H₂-dominated atmosphere is atomic hydrogen (H). The abundance of atomic hydrogen is five orders of magnitudes higher than that of hydroxyl radical (Figure 4). The source of both H and OH is water vapor photodissociation. In the H₂-dominated atmosphere, most of the OH molecules produced from water vapor photodissociation react with H₂ to reform H₂O and produce H through the reactions

$$\begin{eqnarray*} && \hspace{46pt}{\rm H_2O + h\nu \rightarrow H + OH,}\hspace{8pc} ({\rm C1})\\ && \hspace{46pt} {\rm OH + H_2 \rightarrow H + H_2O}.\hspace{7.9pc}({\rm C2}) \end{eqnarray*}$$

Therefore the abundance of atomic hydrogen is much higher than that of OH. We note that H production via water photodissociation is much more efficient than the direct photodissociation of H₂, which requires radiation in wavelengths less than 85 nm. As water vapor is the primary source of H and OH in anoxic atmospheres, the amounts of H and OH depends on the mixing ratio of water vapor above the cold trap, which is in turn sensitively controlled by the cold trap temperature. Water vapor mixing ratio spans 3 orders of magnitudes for the cold trap temperature ranging in 160–200 K; consequently, the number densities of H and OH in the atmosphere can easily vary by one order of magnitude, depending on the cold trap temperature. The removal of atomic hydrogen is mainly by recombination to H₂, which can be more efficient with the presence of CO, via

$$\begin{eqnarray*} &&\hspace{3pc}{\rm H} + {\rm CO} + {\rm M} \rightarrow {\rm CHO} + {\rm M},\hspace{7pc}\ ({\rm C3})\\ &&\hspace{16pt}{\rm H} + {\rm CHO} \rightarrow {\rm CO} + {\rm H}_2.\hspace{8.5pc}\ ({\rm C4}) \end{eqnarray*}$$

As a result of low OH abundance in the H₂-dominated atmosphere, both CH₄ and CO are long-lived and therefore well mixed. CH₄ emitted from the surface is slowly oxidized in the atmosphere into methanol, which gives methane a chemical lifetime of 8 × 10⁴ years. CO is produced by CO₂ photodissociation in our benchmark model; it is emitted by volcanoes on Earth at a much lower rate than CO₂, but it can presumably be the main carbon-bearing gas produced by volcanoes if the upper mantle is more reduced than the current Earth's (e.g., Holland 1984). We find that CO is long-lived in the H₂ atmosphere as well.

The lack of efficient atmospheric sink of CH₄ and CO in the H₂ atmosphere implies that surface sink, if any, is the major sink for these two carbon compounds. CH₄ and CO have zero or very small dry deposition velocities on an abiotic planet, so they can build up to significant amounts in the H₂ atmosphere (Figure 4). If a nonzero deposition velocity is adopted for CH₄ and CO, their steady-state mixing ratios will be much lower than the benchmark model. For example, using a deposition velocity for CH₄ that is standard on Earth (∼0.01 cm s⁻¹) results in a mixing ratio of less than 1 ppb for CH₄, compared with a mixing ratio of 6 ppm in the benchmark model. The dry deposition velocity is indeed the controlling factor for the steady-state abundance of the long-lived carbon compounds. In comparison, the emitted sulfur compounds (H₂S and SO₂) are short-lived. In the H₂ atmosphere, sulfur emission from the surface is readily converted to elemental sulfur aerosols (S₈) in the atmosphere.

Interestingly, CO₂ is fairly long-lived and well mixed in the H₂ atmosphere, meaning that the reduction of CO₂ in the H₂-dominated atmosphere is not efficient (Figure 4). Only 1/3 of emitted CO₂ is reduced in the atmosphere, and the rest is deposited to the surface by dry and wet deposition. It is the balance between surface emission and surface deposition that sets the steady-state mixing ratio of CO₂.

4.2.2. N₂-dominated Atmospheres

Both reducing radicals (i.e., H) and oxidizing radicals (i.e., O and OH) are relatively abundant in the N₂ atmosphere compared with the H₂ atmosphere (Figure 4). Like in H₂-dominated atmospheres, H abundance is orders of magnitude larger than OH abundance, because most of OH molecules from water photolysis react with H₂ to reform H₂O and produce H. The molecular hydrogen that consumes OH and boosts H is emitted volcanically, so a general N₂ atmosphere can be more oxidizing (i.e., having lower H and higher OH) if there is a lower volcanic H₂ emission and mixing ratio than the benchmark model.

An important feature of the N₂-dominated atmosphere is that both H and OH are relatively abundant near the surface. Comparing the N₂-dominated atmosphere with the H₂-dominated atmosphere (Figure 4), we find that both the OH and the H number densities are higher near the surface due to the lower H₂ number density. The relatively high OH number density leads to relatively fast removal of CO by

$$\begin{equation*} \hspace{5pc}{\rm CO} + {\rm OH}\rightarrow {\rm CO}_2 + {\rm H},\hspace{7pc} ({\rm C5}) \end{equation*}$$

as shown in Figure 4. With a low CO abundance the recombination of H via reactions (C3) and (C4) is inefficient. Counterintuitively, a high OH number density helps preserve H in this specific case. This example shows the complexity and the nonlinearity of an atmospheric chemical network. The feature of simultaneous high OH and H abundances near the surface is sensitive to the specification of surface hydrogen emission and eddy diffusion coefficients (see Section 4.3 for relevant rationale).

The chemical lifetimes of CH₄ and CO mainly depend on the amount of OH. In the N₂ atmosphere, CH₄ is well mixed because its chemical lifetime is long. CH₄ is photodissociated and oxidized slowly in the atmosphere into methanol (CH₄O) with a chemical lifetime of ∼6 × 10⁴ years. The photolysis of methane is a secondary source of atomic hydrogen, which concentrates at the pressure level of 10 Pa. Interestingly, methane photolysis causes the apparent trough of the O₂ mixing ratio profile at ∼10 Pa. It is the shielding of the UV radiation that dissociates methane by methane itself that determines this pressure level (see Appendix D for an analytical formula for assessing the pressure level at which the photolysis of a certain gas is important).

CO₂ is actively photodissociated into CO and O in the upper atmosphere, but most of the CO produced is efficiently converted back to CO₂ by reacting with OH (reaction (C5)). In equilibrium, the net chemical removal of CO₂ is minimal, so the CO₂ mixing ratio is set by the balance between emission and deposition. The steady-state amount of CO₂ can be directly estimated by the mass balance between emission and deposition, viz.,

$$\begin{equation} F = f_1 = \frac{n_1}{N_1} = \frac{\Phi _{1/2}}{V_{\rm DEP}N_1}, \end{equation} \tag{ 23 }$$

where the overall mixing ratio (F) is the same as the near-surface mixing ratio (f₁), and Φ_1/2 is the surface emission rate. For CO₂ in the N₂-dominated atmosphere, using Equation (23) with V_DEP(CO₂) = 1 × 10⁻⁴ cm s⁻¹, we find that the steady-state mixing ratio is 1.3 × 10⁻⁴, consistent with the full photochemical model (see Table 6 and Figure 4). The steady-state mixing ratio of a long-lived volcanic gas that is primarily removed by surface deposition is inversely proportional to its dry deposition velocity.

Sulfur-bearing gases emitted from the surface are effectively converted into elemental sulfur and sulfuric acids. Elemental sulfur, mainly in the condensed phase (i.e., aerosols), is the major sulfur-bearing species in the steady state because of relatively high H₂ mixing ratio in the benchmark model. In separate numerical simulations we find that sulfuric acid aerosols may outnumber elemental sulfur aerosols when the H₂ emission is reduced by more than 1 order of magnitude and the atmosphere is more oxidizing than the benchmark model.

4.2.3. CO₂-dominated Atmospheres

CO₂ photodissociation produces CO, O, and O₂. Atomic oxygen is the most abundant reactive radicals in the CO₂-dominated atmosphere (Figure 4), because H is readily removed by O₂ and OH is readily removed by CO. As a result of low H and OH, CO is long-lived in the atmosphere and can build up to very high mixing ratios in the atmosphere depending on its dry deposition velocity; and CH₄ emitted from the surface is also long-lived in the atmosphere with a chemical lifetime of ∼6 × 10⁴ years. The steady-state abundance of CO and CH₄ is therefore controlled by their dry deposition velocities. Most of the emitted SO₂ is deposited to the surface, because there are very few OH or O radicals near the surface. A small fraction of the SO₂ is transported upwards and converted into sulfuric acid aerosols in the radiative layer.

Based on our benchmark model of CO₂-dominated atmospheres, we now revisit whether photochemically produced O₂ can cause false positive for detecting oxygenic photosynthesis. Oxygen and ozone are the most studied biosignature gases for terrestrial exoplanet characterization (e.g., Owen 1980; Angel et al. 1986; Léger et al. 1993, 1996; Beichman et al. 1999). One of the main concerns of using O₂/O₃ as biosignature gases is that O₂ may be produced abiotically from photodissociation of CO₂. A number of authors have studied the abiotic production of oxygen in terrestrial atmospheres, either for understanding prebiotic Earth's atmosphere (e.g., Walker 1977; Kasting et al. 1979; Kasting & Catling 2003), or for assessing whether abiotic oxygen can be a false positive for detecting photosynthesis on habitable exoplanets (Selsis et al. 2002; Segura et al. 2007). Selsis et al. (2002) found that photochemically produced oxygen can build up in CO₂-dominated atmospheres without any surface emission or deposition. The results of Selsis et al. (2002) was challenged by Segura et al. (2007), who had additionally considered the surface emission of reducing gases including H₂ and CH₄, and concluded that abiotic oxygen from CO₂ photodissociation is not likely to build up in the atmosphere on a planet having an active hydrological cycle.

We find that the steady-state number density of O₂ and O₃ in the CO₂-dominated atmosphere is mainly controlled by the surface emission of reducing gases such as H₂ and CH₄, and without surface emission of reducing gas photochemically produced O₂ can build up in a 1 bar CO₂-dominated atmosphere. In addition to the benchmark model, we have simulated CO₂-dominated atmospheres with relatively low and zero emission rates of H₂ and CH₄ (Table 7 and Figure 6). We see that the O₂ mixing ratio near the surface increases dramatically in 1 bar CO₂ atmospheres when the emission of reducing gases decreases. O₂ is virtually nonexistent at the surface for the Earth-like emission rates of H₂ and CH₄, but O₂ mixing ratio can be as high as 10⁻³ if no H₂ or CH₄ is emitted (Figure 6). In particular, if no H₂ or CH₄ is emitted, the O₃ column integrated number density can reach one third of the present-day Earth's atmospheric levels (Table 7), which constitutes a potential false positive.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 7. Mixing Ratios of O₂ and O₃ and Redox Budget for CO₂-dominated Atmospheres on Rocky Exoplanets having Different Surface Emission of Reducing Gases

Chemical	CO2-dominated Atmospheres
species	Φ1/2(H2) = 3 × 1010 cm s−1	Φ1/2(H2) = 3 × 109 cm s−1	No H2 Emission
Column-averaged mixing ratio
O2	6.4E−7	3.8E−6	1.3E−3
O3	7.0E−11	3.7E−10	1.3E−7
Redox balance
Atmospheric escape
H	−7.0E+8	−2.4E+8	−1.2E+6
H2	−5.9E+10	−5.8E+9	−2.0E+6
Surface emission
H2	6.0E+10	6.0E+9	0
CH4	2.4E+9	2.4E+9	0
H2S	1.8E+9	1.8E+9	1.8E+9
Dry and wet deposition
O3	0	6.3E+4	1.1E+10
HO2	3.1E+4	7.8E+7	1.2E+8
H2O2	1.7E+5	1.5E+9	2.7E+9
CO	−3.7E+9	−5.6E+9	−1.5E+10
CH2O	−1.3E+6	−7.2E+4	0
Organic Haze	−2.0E+6	−8.1E+1	0
H2S	−4.9E+8	−2.5E+8	−3.2E+8
H2SO4	1.2E+8	6.3E+7	2.9E+7
Balance	4.9E+3	−9.5E+2	−1.3E+3

Notes. The redox number for each species is defined according to Equation (22). In the redox balance, all values have unit of molecule cm⁻² s⁻¹, and defined as positive for hydrogen flux into the atmosphere. The redox budget for our atmosphere models is balanced, meaning that the atmosphere is not becoming more oxidized or reduced.

Download table as:  ASCIITypeset image

In the 1 bar CO₂-dominated atmosphere O₃ can potentially build up to a false-positive level even on a planet with active hydrological cycle. Segura et al. (2007) have based their conclusion on simulations of 20% CO₂ 1 bar atmospheres with and without emission of H₂ and CH₄ and simulations of 2 bar CO₂ atmospheres with emission of H₂ and CH₄. We have been able to reproduce all results of Segura et al. (2007) quantitatively to within a factor of two. Where we differ from Segura et al. (2007) is that we successfully simulated high CO₂ 1 bar atmospheres with minimal volcanic reducing gas emission (Figure 6). This is a parameter space that Segura et al. (2007) did not cover, but we find that this is the parameter space for high abiotic O₂.

We here provide rationale for our specification of the atmospheric temperature profile, the eddy diffusion coefficients, and the dry deposition velocities.

First, the surface temperatures for the three scenarios are assumed to be 288 K, and the semi-major axis of the planet is adjusted according to appropriate amounts of greenhouse effect. The semi-major axes around a Sun-like star for H₂-, N₂-, and CO₂-dominated atmospheres are found to be 1.6 AU, 1.0 AU, and 1.3 AU, respectively. We have compared the planetary thermal emission flux and the incidence stellar flux to determine the semi-major axis, a similar procedure as Kasting et al. (1993) taking into account CO₂, H₂O, and CH₄ absorption, H₂ collision-induced absorption, 50% cloud coverage, and a Bond albedo of 30%. The temperature profiles are assumed to follow the appropriate dry adiabatic lapse rate (i.e., the convective layer) until 160 K (H₂ atmosphere), 200 K (N₂ atmosphere), and 175 K (CO₂ atmosphere) and to be constant above (i.e., the radiative layer). We simulate the atmosphere up to the altitude of about 10 scale heights or the pressure level of 0.1 Pa. The adopted temperature profiles are consistent with significant greenhouse effects in the convective layer and no additional heating above the convective layer for habitable exoplanets. The results discussed above do not change significantly if these temperature profiles are changed by several tens of K.

Second, we have used eddy diffusion coefficients empirically determined on Earth and scaled the values to account for different mean molecular masses. The eddy diffusion coefficient at a certain pressure level is assumed to be that of current Earth's atmosphere at the same pressure level (see Figure 1), scaled by 6.3, 1.0, and 0.68 for the H₂-, N₂-, and CO₂-dominated atmospheres, respectively, to account for different dominant molecules. We have roughly scaled the eddy diffusion coefficient assuming K∝H₀ where H₀ is the atmospheric scale height. The reasoning of the scaling is as follows. According to the mixing length theory, K∝lw, where l is the typical mixing length and w is the mean vertical velocity. The mixing length is a fraction of the pressure scale height (e.g., Smith 1998). The mean vertical velocity is related to the vertical convective energy flux, as F∝pw, where p is the pressure (e.g., Lindzen 1990). For a certain planet F should have the same order of magnitude for different atmosphere compositions. As we apply the scaling from pressure surface to pressure surface, we have roughly K∝H₀. The scaling is an approximation and we only intend to provide a consistent description of eddy diffusion for atmospheres with very different mean molecular mass. The eddy diffusion coefficients for the three scenarios are also consistent with their temperature profiles, featuring minima near the tropopause. The general results discussed above are not sensitive to the variation of eddy diffusion coefficients by an order of magnitude.

Third, the typical deposition velocity on Earth is sometimes not directly applicable for terrestrial exoplanets, because the major surface sink of a number of gases (notably, H₂, CH₄, CO) on Earth is actually microorganisms. In this paper we focus on the scenarios assuming no biotic contribution in neither surface emission nor surface deposition. For H₂ and CH₄, a sensible dry deposition velocity without biotic surface sink is zero (J. Kasting, 2012, private communication). For CO₂, the dry deposition velocity should match with the timescale of weathering and carbonate formation, which is 10,000 years (Archer 2010). For CO, the rate limiting step for converting CO into bicarbonate has been proposed to the hydration of CO in the ocean, which corresponds to a deposition velocity of 1.0 × 10⁻⁹ to 10⁻⁸ cm s⁻¹ (Kharecha et al. 2005). The adopted or assumed values of dry deposition velocities are tabulated in Table 5.

For completeness, we now comment on our assumptions for nitrogen chemistry and organic haze. We do not track the bulk nitrogen cycle, instead assuming 10% N₂ for H₂- and CO₂-dominated atmospheres. Our model does not treat abiotic nitrogen fixation by lightning (see Kasting & Walker 1981; Zahnle 1986; Tian et al. 2011 for specific analysis), so N₂ is considered as inert in our models. We include the formation of elemental sulfur aerosols and sulfuric acid aerosols in our models; but we do not include the less-understood hydrocarbon chemical network for hydrocarbon molecules that have more than two carbon atoms or the formation of hydrocarbon haze. Organic haze may be formed in anoxic atmospheres based on methane photolysis and hydrocarbon polymerization. The number density of hydrocarbons that have more than four carbon atoms (usually considered to be condensable) is ∼4 orders of magnitude smaller than the number density of C₂H₆ at the steady state (e.g., Allen et al. 1980; Yung et al. 1984; Pavlov et al. 2001). To account for the loss of carbon due to possible organic haze formation and deposition, we apply an ad hoc dry deposition velocity of 1.0 × 10⁻⁵ cm s⁻¹ for C₂H₆. This small velocity results from a scaling based on typical sub-micron particle deposition velocity (∼0.2 cm s⁻¹, Sehmel 1980) and the number density ratio of C₄ hydrocarbons and C₂H₆ (e.g., Yung et al. 1984). An accurate treatment for organic haze is not the purpose of this paper, because the efficiency of haze formation depends on oxidation states of atmospheres, stellar UV radiation, and a number of less-understood reaction rates. Nonetheless, for all simulated scenarios, the carbon loss due to haze formation is less than 1% of the methane emission flux; therefore we do not expect that our simplification of organic haze formation would impact our results. We will fully explore the formation of organic haze and its implication on atmospheric chemistry in a separate paper (R. Hu et al. 2013, in preparation).