Evolved Climates and Observational Discriminants for the TRAPPIST-1 Planetary System

We use the Spectral Mapping and Atmospheric Radiative Transfer code (SMART), originally developed by David Crisp (Meadows & Crisp 1996; Crisp 1997), for both high-resolution spectra and within the climate model to compute radiative fluxes and heating rates. SMART is capable of producing accurate synthetic planetary spectra (e.g., Robinson et al. 2011). SMART uses the Discrete Ordinate Radiative Transfer code (DISORT; Stamnes et al. 1988, 2000) to numerically compute the radiation field. SMART can combine the radiative fluxes from overlapping stellar and thermal sources and includes the extinction in each layer from absorbing gases (due to visible-to-far-infrared vibrational–rotational transitions, collision-induced absorption, and UV predissociation bands), Rayleigh scattering, cloud and aerosol scattering, a wavelength-dependent surface albedo, and a stellar source spectrum. We also use SMART to produce transit transmission spectra, including refraction (Robinson 2017).

We use VPL Climate, a general-purpose, 1D radiative–convective equilibrium, terrestrial planet climate model (Meadows et al. 2018a; Robinson & Crisp 2018). This model uses SMART for rigorous spectrum-resolving (line-by-line) radiative transfer, coupled to a mixing-length convection parameterization. As input, VPL Climate requires an initial atmospheric state (consisting of a temperature profile and gas mixing ratios), physical parameters, and fundamental laboratory and experimental data.

We use a 1D coupled climate–photochemical model in this work, rather than a 3D GCM, to more rigorously and self-consistently explore both the radiative and photochemical effects of the TRAPPIST-1 stellar spectrum on atmospheric composition and climate for non-Earth-like terrestrial planet atmospheres. For Proxima Centauri b, Meadows et al. (2018a) performed a 1D climate comparison with the 3D GCM results of Turbet et al. (2016) and demonstrated that differences in globally averaged temperatures between 1D and 3D models were small. For N₂- and CO₂-dominated atmospheres, the difference was less than 20 K when the planet was synchronously rotating, and 5 K or less when the planet was in a 3:2 spin–orbit resonance state. In slower, synchronously rotating planets, enhanced day–night temperature contrasts (which cannot be spatially resolved by 1D models) may induce atmospheric collapse (Turbet et al. 2016, 2018). However, even if synchronously rotating, the TRAPPIST-1 planets, with orbital periods of 1.5–19 days, are in an orbital regime closer to rapid rotators, with the possibility of sufficient Coriolis force to drive zonal circulation and help prevent atmospheric collapse (Yang et al. 2014; Kopparapu et al. 2016). For the desiccated planets that we model here, we also assume dense (≥10 bar) atmospheres that will have smaller day–night contrast and are more likely to be stable against collapse (Hu & Yang 2014; Meadows et al. 2018a; Turbet et al. 2018). Although our habitable planet example for TRAPPIST-1 e has only a 1 bar atmosphere, Hu & Yang (2014) showed that day–night temperature contrasts for synchronously rotating planets with oceans were much smaller when a full dynamical ocean model, which strongly augments atmospheric heat transport, was used in a 3D GCM. With a dynamical ocean, even a relatively thin atmosphere was more stable against collapse.

The VPL Climate model components are summarized below and described in more detail in Appendix A. The inputs are described in Section 2.6.

To self-consistently model the effects of M dwarf host-star UV radiation on our modeled atmospheres, we couple VPL Climate to a publicly available atmospheric chemistry model.⁶ This code is based on Kasting et al. (1979) and Zahnle et al. (2006) and can simulate a range of planetary redox states, including high-oxygen atmospheres (Arney et al. 2016; Schwieterman et al. 2016; Meadows et al. 2018a). This photochemical model is described in detail in Meadows et al. (2018a). Briefly, it divides the nominal model atmosphere into plane-parallel layers and assumes hydrostatic equilibrium. A vertical transport scheme includes molecular and eddy diffusion. Boundary conditions can be set for all species at the top and bottom of the atmosphere, including initial gas mixing ratios, outgassing fluxes, and surface deposition velocities. We describe the specific boundary conditions for each type of atmosphere in Section 3.

Here, we updated the photochemistry model by extending the photochemically active wavelength range considered by the model to include Lyα for a wide range of species. Updates to the model are described in Appendix B, and validations for Earth and Venus are given in Appendix C.

The climate and chemistry models can be run independently or as a coupled model. When the climate model is run independently, the convergence criterion is met when all atmospheric layers are flux-balanced within 1 W m⁻² and the heating rate for each layer is less than 10⁻⁴ K per day. When run independently, the atmospheric chemistry model tracks changes in gas concentrations in each time step and computes the relative error. It adjusts the next time step based on the largest error. The model checks the time step length, and when it reaches 10¹⁷ s in less than ∼100 time steps, we consider the model to be converged.

To run coupled climate–photochemistry experiments, we pass the converged profile from VPL Climate (temperature and gas mixing profiles) to the chemistry code, which is run to convergence. The modified gas mixing ratios from the chemistry code are then passed back to VPL Climate to compute a new equilibrium temperature structure and condensable gas mixing profiles. These profiles are then passed back to the chemistry model, and this process is repeated until global convergence is achieved.

Plausible upper limits of water and oxygen inventories of postrunaway atmospheres are modeled using the atmospheric escape framework introduced by Luger et al. (2015) and Luger & Barnes (2015). Briefly, we compute the XUV-driven, energy-limited escape rate of hydrogen modified by the hydrodynamic drag of oxygen as in Hunten (1973), assuming photolysis is fast enough that it is not the rate-limiting factor. As in Luger & Barnes (2015), we assume surface sinks are inefficient at removing photolytically produced oxygen, yielding a strict upper bound on the amount of oxygen retained in the atmosphere. Unlike Luger & Barnes (2015), we model the decrease in the efficiency of the flow due to the increase in the mixing ratio of oxygen, which causes the flow to transition to diffusion-limited escape (Schaefer et al. 2016). We also relax the assumption of constant escape efficiency _xuv, instead modeling it as a function of the incident XUV flux as in Bolmont et al. (2017). For the TRAPPIST-1 planets, _xuv varies between about 0.01 (at high XUV flux) and 0.10 (at moderate XUV flux).

Our climate, photochemical, and atmospheric escape models require a number of inputs that tailor the model to the planetary atmosphere studied. Inputs include planetary properties (e.g., orbit and radius), atmospheric gas mixing ratios for each constituent in each layer, gas absorption properties (i.e., cross section data and line lists), thermodynamic data for condensable gases, particle optical properties for aerosols, wavelength-dependent surface albedo data, and stellar spectral energy distribution (SED), including an estimate of total XUV flux over time.

2.6.1. Planetary Properties and Model Atmospheres

2.6.2. Gas Absorption Data

2.6.3. Aerosol Optical Properties

The single-scattering optical properties of clouds and aerosols are defined in precomputed tables for a range of specific particle types. Each particle type is defined by its composition, its associated wavelength-dependent refractive indices, and its particle size and shape distributions. Given these properties, a single-scattering model is used to compute the wavelength-dependent extinction and scattering efficiencies (Q_ext and Q_abs) and the phase function moments (for Mie scattering) or the particle asymmetry parameter (g for some nonspherical particles) used in our radiative transfer calculations.

For the cloudy aqua planet, we specify Earth-like cirrus (water–ice) and stratocumulus (liquid water) clouds, each of cumulative optical depth τ = 5. For cirrus clouds, the optical properties  from B. Baum's Cirrus Optical Property Library (Baum et al. 2005)⁸ are used. The particles consist of a distribution of 45% solid columns, 35% plates, and 15% 3D bullet rosettes, spanning 2–9500 μm with a cross section weighted mean diameter of 100 μm. Stratocumulus cloud optical properties are based on refractive indices of water from Hale & Querry (1973) and calculated using Mie scattering, with a two-parameter gamma distribution (a = 5.3, b = 1.1) and mean particle radius 4.07 μm.

For Venus-like sulfuric acid aerosols, we use refractive indices for a range of specific acid concentrations between 25% and 100% from Palmer & Williams (1975). The sulfuric acid concentration is computed assuming vapor pressure equilibrium of the H₂O and H₂SO₄ gases with the condensed solution. The effective optical properties are then calculated by spline interpolation of the tables. Our photochemical model calculates the monodisperse aerosol radii at every atmospheric layer given the coagulation, sedimentation, and diffusion timescales (see, e.g., Pavlov et al. 2001) in phase equilibrium. For climate and spectral modeling, we convert the monodisperse distributions into log-normal distributions with the modal radii equal to the monodisperse particle radii and geometric standard deviation equal to 0.25, similar to the particle distribution findings for individual Venus aerosol modes as described in Crisp (1986). The differential optical depth profile used in the radiative transfer calculations for each planet is calculated from the log-normal distribution, assuming that the total mass in each layer is the same as in the monodisperse distribution. The optical depth dτ is calculated by

$$\begin{eqnarray}&&d\tau ={dz}{\int }_{{R}_{\min }}^{{R}_{\max }}\pi {r}^{2}{Q}_{\mathrm{ext}}(\lambda ,r)n(r){dr},\end{eqnarray} \tag{ 1 }$$

where Q_ext is the extinction efficiency at a reference wavelength and n is the number density distribution.

2.6.4. Thermodynamic Data

2.6.5. Surface Spectral Albedo

We use different surface compositions based on the type of climate we are simulating. Basalt was used for the volcanic Venus-like planets, and a desert surface (primarily kaolinite) was used for the post-ocean-loss O₂-dominated atmospheres. For the aqua planet, we used open ocean. Surface reflectance data are from the U.S.G.S. spectral library (Clark et al. 2007),¹⁰ except for the ocean longward of 2 μm, which is from the ASTER spectral library (Baldridge et al. 2009)¹¹  (see Figure 1). The integrated, energy-flux-weighted surface albedos for these surfaces for our TRAPPIST-1 stellar spectrum are 0.09 (basalt), 0.37 (desert), and 0.02 (ocean). For comparison, this albedo for snow used in the VPL Spectral Earth Model (Robinson et al. 2011) is 0.16.

Zoom In Zoom Out Reset image size

2.6.6. Stellar Spectra

Characterizing the climate and computing spectra require the stellar spectrum of TRAPPIST-1 (2MASS J23062928-0502285). The available data on this star in the UV to IR are limited to photometry and to measurements of Lyα from the Hubble Space Telescope (Bourrier et al. 2017). TRAPPIST-1 has an effective temperature of 2559 ± 50 K, [Fe/H] =0.04 ± 0.08, ${M}_{* }=0.0802\pm 0.0073\,{M}_{\odot }$, and ${R}_{* }=0.117\,\pm 0.0036\,{R}_{\odot }$ (Gillon et al. 2017). To model the panchromatic UV–NIR stellar spectrum necessary for line-by-line climate–photochemical modeling, we use the 2500 K, [Fe/H] = 0.0, log g = 5.0 spectral model of the PHOENIX v2.0 spectral database, which spans from 0.25 to 5.5 μm and has been shown to faithfully reproduce the visible–NIR spectra of M dwarfs (France et al. 2013). We compared this PHOENIX model to the photometric observations (Gillon et al. 2016) and found that the model SED, binned to available photometric bands, is within 3σ.

The intensity and wavelength dependence of the stellar UV spectrum is a critical factor affecting atmospheric chemistry and the subsequent atmospheric composition and climate. For example, Grenfell et al. (2007) noted changes in the atmospheres of planets around F-G-K stars as a function of orbital distance, including the impact of outgassed species on climate, while Grenfell et al. (2014) noted that UV changes within 200–350 nm have the largest effect on ozone levels for planets around late-type M dwarfs. In the absence of complete UV observations of TRAPPIST-1, we fit a UV spectrum compiled for Proxima Centauri (Meadows et al. 2018a). We scaled Proxima's UV spectrum to the Lyα flux ratio of TRAPPIST-1 to Proxima Centauri, a factor of one-sixth at 1 au (Bourrier et al. 2017). The UV flux was stitched to the PHOENIX SED at around 0.3 μm, where photospheric emission began to dominate (i.e., where the flux from the PHOENIX model began to rise above the UV). This method provides only a possible UV flux environment for the planets of TRAPPIST-1 and does not represent the star's actual activity (except the reconstructed Lyα). Given the variable activity of M dwarfs, UV observations of TRAPPIST-1 are needed to properly characterize its planets. Our compiled spectrum is shown in Figure 2 and is available for download online using the VPL Spectral Explorer.¹²

Zoom In Zoom Out Reset image size

2.6.7. XUV Flux

To determine plausible oxygen inventories for our O₂-dominated atmospheres, we computed the loss of hydrogen and oxygen during the pre-main-sequence phase of TRAPPIST-1 (see Figure 3). The maximum initial water inventory considered was 20 Earth oceans (by mass), the approximate maximum ocean loss for planets b, c, and d under our assumptions. Planets e, f, and g may each lose between three and six Earth oceans of water, while h may lose less than one ocean. This ocean loss results in the generation of dense O₂ atmospheres with maximum pressures (assuming no sinks) from 22 to 5000 bar.

Zoom In Zoom Out Reset image size

The amount of water each planet may have lost and the subsequent surface pressure of atmospheric oxygen vary due to instellation, atmospheric scale height, and planetary surface area. The degree of H₂O loss is a strong function of instellation and time spent interior to the HZ inner edge. For example, as the host star dims, within 50–400 Myr after formation, planets h, g, f, and e successively enter the habitable zone and cease to lose significant amounts of H₂O, while planets b, c, and d continue to lose H₂O (Figure 3). Although the orbits (and levels of instellation received) for these planets are well spaced in resonance chains (Luger et al. 2017b), they exhibit variations in H₂O loss not related to their instellation levels. Under diffusion-limited escape, we assume that H₂O loss is proportional to the inverse of scale height (Hunten 1973), which is in turn proportional to temperature and inversely proportional to gravity and mean molecular mass. The small variations in H₂O loss among b–d and e–h are in part due to the different scale heights for each planet.

In comparison to ocean loss, the resultant surface pressure of equivalent atmospheric oxygen generated by water photolysis does not show a similar increasing trend with instellation; the sequence we calculated for the TRAPPIST-1 system (from highest to lowest O₂ surface pressure) is c, d, b, e, g, f, h. Compared to planet b, c loses similar amounts of H₂O and is similar in size, but because c has a much larger mass and stronger surface gravity, it builds up a higher surface pressure of equivalent atmospheric O₂. That is, c can accrue a similar mass of remaining O₂, but the stronger surface gravity results in a larger surface pressure (since surface pressure is the product of surface gravity and the integrated atmospheric mass per unit area). The sequence of e, g, f, and h is also affected by each planet's respective gravitational acceleration, as opposed to an effect of radius; while f and g are similar in radius, their nominal surface gravities are very different, which results in a higher surface pressure of oxygen for g.

In the following subsections, we describe our input assumptions and resultant environmental states (atmospheric temperature profiles as a function of pressure, and vertical profiles of principal atmospheric constituent abundances) for each atmosphere. Globally averaged surface temperatures for all modeled climates are given in Figure 4, and O₃ column densities are given in Table 3.

Zoom In Zoom Out Reset image size

Table 3.  Ozone Column Densities [cm⁻²]

Planets	b	c	d	e	f	g	h
Aqua planet, clear, 1 bar				$3.1\times {10}^{18}$
Aqua planet, cloudy, 1 bar				$2.9\times {10}^{18}$
O2, desiccated, 10 bar	$2.5\times {10}^{18}$	$6.3\times {10}^{18}$	$1.1\times {10}^{19}$	$4.7\times {10}^{19}$	$2.2\times {10}^{22}$	$6.0\times {10}^{22}$	$1.3\times {10}^{23}$
O2, desiccated, 100 bar	$7.2\times {10}^{17}$	$2.4\times {10}^{18}$	$6.4\times {10}^{18}$	$8.8\times {10}^{18}$	$6.1\times {10}^{22}$	$2.7\times {10}^{23}$	$1.4\times {10}^{24}$
O2, outgassing, 10 bar	$1.6\times {10}^{18}$	$5.2\times {10}^{18}$	$9.9\times {10}^{18}$	$2.7\times {10}^{19}$	$8.5\times {10}^{20}$	$2.4\times {10}^{21}$	$4.0\times {10}^{21}$
O2, outgassing, 100 bar	$1.7\times {10}^{19}$	$2.1\times {10}^{20}$	$2.7\times {10}^{21}$	$2.4\times {10}^{22}$	$7.2\times {10}^{23}$	$6.6\times {10}^{24}$	$1.8\times {10}^{25}$
Venus-like, clear, 10 bar	$2.3\times {10}^{15}$	$2.1\times {10}^{15}$	$5.4\times {10}^{15}$	$1.0\times {10}^{15}$	$6.9\times {10}^{15}$	$8.6\times {10}^{15}$	$1.5\times {10}^{16}$
Venus-like, cloudy, 10 bar		$1.0\times {10}^{15}$	$1.1\times {10}^{16}$	$4.6\times {10}^{15}$	$7.6\times {10}^{15}$	$9.1\times {10}^{15}$	$1.5\times {10}^{16}$
Venus-like, clear, 92 bar	$1.8\times {10}^{15}$	$1.9\times {10}^{15}$	$5.1\times {10}^{15}$	$1.3\times {10}^{15}$	$6.7\times {10}^{15}$	$9.1\times {10}^{15}$	$1.6\times {10}^{16}$
Venus-like, cloudy, 92 bar		$1.0\times {10}^{15}$	$9.5\times {10}^{14}$	$4.2\times {10}^{15}$	$7.8\times {10}^{15}$	$9.8\times {10}^{15}$	$1.6\times {10}^{16}$

Note. Compare to global-average Earth ozone column depth of $8\times {10}^{18}$ cm⁻² (https://ozonewatch.gsfc.nasa.gov/).

Download table as:  ASCIITypeset image

3.2.1. O₂-dominated Atmospheres: Desiccated

In Figure 5, we present environmental states for completely desiccated, oxygen-dominated atmospheres for the seven TRAPPIST-1 planets, assuming 10 or 100 bar of surface pressure. These surface pressure values were chosen as representative of possible states, because although TRAPPIST-1's pre-main-sequence evolution may produce up to thousands of bar of O₂ in the most extreme case (Figure 3), much of this O₂ is likely to be reincorporated into the land surface, ocean, lithosphere, or mantle (e.g., Hamano et al. 2013; Schaefer et al. 2016; Wordsworth et al. 2018) or, depending on each planet's magnetic field, lost to nonhydrodynamic, top-of-atmosphere processes over time (e.g., Khodachenko et al. 2007; Lammer et al. 2007, 2011; Ribas et al. 2016), such as the "electric wind" on Venus (Collinson et al. 2016). While an abiotic O₂-dominated planetary atmosphere has yet to be discovered, and Venus does not retain significant quantities of atmospheric oxygen, it may be possible to test the existence of this hypothesized ocean-loss-generated atmosphere with JWST.

Zoom In Zoom Out Reset image size

We initialized these atmospheres with 95% O₂, 0.05 bar CO₂ (consistent with Meadows et al. 2018a), and the remainder N₂, each specified as a fixed surface mixing ratio, with no outgassing. Since this case assumes a planet that was severely depleted in volatiles in both its atmosphere and mantle, these atmospheres are considered completely desiccated (i.e., they contain no water, nor any other hydrogen-bearing species). The remaining composition is dominated by photochemically produced trace gases.

The composition of photochemically produced trace gases, principally O₃ and CO, increases monotonically as a function of irradiation from the parent star (i.e., as a function of semimajor axis). The equilibrium concentration of ozone is due to the competing effects of vertical transport and the Chapman cycle, whereby photolysis of O₂ produces atomic oxygen that can combine with O₂ to form O₃, which is directly lost by photolysis if atomic oxygen combines to form O₂. For the hotter planets (b–d), O₃ photolysis is efficient down to the surface from absorption in the broad Chappuis band. However, the colder planets receive less UV flux as a result of distance from the star, so the photolysis rates are slower. In the coldest atmospheres (f, g and h), eddy transport can overwhelm photolysis and reaction rates for O₃, resulting in well-mixed vertical profiles and saturated O₃ absorption bands. Furthermore, photolysis of O₂ drops faster with semimajor axis because its cross section is limited to the FUV and easily saturates compared to O₃, which includes a strong NUV band and the weak Chappuis band in visible wavelengths, which are not easily saturated.

The monotonic increase in CO with semimajor axis in these atmospheres is due to CO₂ photolysis and lack of recombination. In the terrestrial solar system atmospheres, the recombination to CO₂ is efficient due to photolytic water products (e.g., OH⁻). Since this environment is devoid of hydrogen, CO can only recombine with atomic oxygen under a slow, density-dependent three-body reaction (Gao et al. 2015).

These modeled environments exhibit surface temperatures (152–406 K, from h to b) that are similar to their equilibrium values (173–400 K; Gillon et al. 2017; Luger et al. 2017b), largely due to low greenhouse gas abundances. The 100 bar cases are cooler at the surface than the 10 bar cases, which indicates an antigreenhouse effect is occurring. The atmospheric structures for these completely desiccated planets each exhibit a stratospheric temperature inversion, like Earth, due to the radiative effects of the photochemically produced O₃ and O₂–O₂ collision-induced absorption (CIA), which is significant in these dense O₂ atmospheres. A sensitivity test removing each constituent to determine the radiative effect suggests that O₃ and O₂–O₂ contribute equally to this heating structure, but via different mechanisms: O₂–O₂ directly absorbs M dwarf NIR irradiation, while O₃ traps outgoing thermal radiation.

The pressure level at which the stratosphere peaks is different than on Earth due to the O₂–O₂ CIA. Earth's stratosphere is caused by O₃ NUV absorption and peaks around ∼0.001 bar, while the pressure level at which the peak in these O₂ atmospheres occurs is much closer to the surface, around 1–10 bar, due to the density dependence of CIA. This stratospheric absorption induces heating in the atmosphere, exceeding heating at the planetary surface, causing a surface temperature inversion similar to Earth and Mars at night (e.g., André & Mahrt 1982; Hinson & Wilson 2004) and inhibiting a troposphere from forming. The causes for these conditions on Earth and Mars at night are due instead to rapid cooling of the surface, rather than warming of the atmosphere.

3.2.2. O₂-dominated Atmospheres: Outgassing

Here we consider the case where a planet has lost its oceans but continues to outgas from a volatile-rich mantle, as is likely the case on Venus (e.g., Bullock & Grinspoon 2001). For other planetary systems, including TRAPPIST-1, this could also be due to delivery by impactors (Chyba 1990; Kral et al. 2018). For all seven planets, we specify Earth-like outgassing fluxes (F), which are relatively well known: 1.68 × 10¹¹ cm² s⁻¹ H₂O (Burton et al. 2000, ∼$10{F}_{{\mathrm{CO}}_{2}}$), 1.0 × 10⁸ cm² s⁻¹ H₂S, 1.35 × 10⁹ cm² s⁻¹ SO₂ (Carn et al. 2017), 6.88 × 10⁹ cm² s⁻¹ CO₂ (Williams et al. 1992, ∼$3.5{F}_{{\mathrm{SO}}_{2}}$), and 3.7 × 10⁵ cm² s⁻¹ OCS (Belviso et al. 1986, ∼${10}^{-4}{F}_{{\mathrm{CO}}_{2}}$). These fluxes represent volcanic outgassing, so they are distributed within one scale height of the surface (∼6–20 km). We specify minimal nonzero dry deposition rates of 10⁻⁸ cm s⁻¹ for all outgassed species to prevent a "runaway" state.

As shown in Figure 5, identical outgassing fluxes result in different levels of trace gases across the seven planets. Some species vary monotonically with irradiation (H₂O, O₃, CO₂), while others do not (CO, SO₂). These trends are modulated by surface gravity (planets b, d, and h have low surface gravity) and atmospheric temperature. There are several interacting chemical networks: carbon, ozone, and sulfur. First consider the carbon network: the primary carbon species are CO₂ and its photolytic product, CO. The CO profiles demonstrate the competing effects of photolysis and vertical transport. The hotter planets have a steep CO gradient, driven by strong photolysis, while the colder planets have well-mixed CO profiles, where vertical transport overcomes photolysis rates and the peak abundance drops. The CO₂ levels are not wholly a function of photolysis; for a given surface pressure, they have different atmospheric column masses (as discussed in Section 3.1 on atmospheric escape). Even though the hotter planets have larger mixing ratios of CO₂, because of hotter temperatures their atmospheres are less dense and have smaller number densities. Photolysis of CO₂ from identical outgassing fluxes results in lower CO₂ inventories (i.e., number densities) in the more irradiated atmospheres. These CO₂ results differ from the fully desiccated environments because of the Earth-like active surface fluxes we assumed here.

The ozone profiles are very similar to the desiccated O₂-dominated planets. The production of O₃ is the same, but the outgassing planets have more sinks, primarily water vapor. As a result, the outgassing cases overall have lower column densities of O₃. Since the outer planets can be very cold, they cannot maintain much water vapor (Figure 5), so they too can accumulate substantial quantities of O₃, particularly g and h.

The abundance and vertical distribution of species in the sulfur network (SO₂, SO₃, and H₂SO₄) are a balance between water vapor availability, temperature, and photolysis rates. The colder atmospheres are water-limited due to condensation and lower photolysis rates, so atmospheric SO₂ inventories are highest. The formation of SO₃ is from the combination of SO₂ with free oxygen atoms. Sulfuric acid is a primary sink for SO₃ and is formed by (Krasnopolsky 2012)

$$\begin{eqnarray}&&{\mathrm{SO}}_{3}+{{\rm{H}}}_{2}{\rm{O}}+{{\rm{H}}}_{2}{\rm{O}}\,\longrightarrow \,{{\rm{H}}}_{2}{\mathrm{SO}}_{4}+{{\rm{H}}}_{2}{\rm{O}}\end{eqnarray} \tag{ 2 }$$

and can either thermally decompose or condense and rain out. The inner planets receive sufficient irradiation to permit efficient conversion of SO₂ to SO₃ and H₂SO₄.

Volatiles in these atmospheres have a large effect on planetary climate and bridge the results of the O₂-desiccated and Venus-like atmospheres. The coldest planets (g and h) maintain temperature profiles similar to the desiccated cases, including near-surface temperature inversions, because they are too cold to maintain quantities of water vapor sufficient to heat the lower atmospheres. Conversely, the remaining planets generate sufficient water abundances to maintain a positive lapse rate. To prevent a full runaway greenhouse, we have limited the stratospheric water relative humidity to 0.1% for b and 1.0% for c, while the remaining planets are limited to 10%. Planets b, c, and d show the tendency to continue into a runaway greenhouse state, which is supported in our results by their  abundances of water vapor well mixed throughout the atmospheric column. Reducing the outgassing of H₂O in our model would more closely approximate the Venus-like H₂O abundance. Only TRAPPIST-1 e has a moderate surface temperature (271–314 K) in this outgassing case, which supports the possibility for it to outgas an ocean and recover habitability.

3.2.3. Venus-like Atmospheres

Here we model Venus-like environments, which may be common around M dwarfs (Kane et al. 2014). For this environmental state, we assume that the complete loss of oceans during the superluminous pre-main-sequence phase removed a potentially strong sink for soluble gases (e.g., CO₂ and SO₂). Subsequent loss or sequestration of oxygen over time (Schaefer et al. 2016), and initial and subsequent outgassing of volatiles, may have allowed a Venus-like high-CO₂ atmosphere to develop. Even though planet h may not have lost an entire Earth ocean (see Figure 3), it could have initially outgassed significant amounts of CO₂. We modeled Venus-like atmospheres for all seven planets (Figure 6)  using the 10 bar boundary conditions used to validate observations of Venus, following Krasnopolsky (2012; see Appendix C). Venus atmospheric chemistry is typically modeled in sections: by thermochemistry, which dominates the lower atmosphere below the cloud deck (e.g., Krasnopolsky 2013), and by photochemistry, which is responsible for the directly observable altitudes at and above the cloud deck (e.g., Krasnopolsky 2012; Zhang et al. 2012). We use photochemistry to model these atmospheres but do not include thermochemical reactions, which could be a useful improvement in future work for the planets with the highest surface temperatures.

Zoom In Zoom Out Reset image size

Our Venus atmospheres are initiated as either 10 or 92 bar of total pressure with Venus-like levels of CO₂ (96.5%) and N₂ (3.5%). For the 92 bar cases, we assumed a constant mixing ratio for every species from 10 bar to the surface. We initiated the photochemical model with the trace-gas lower boundary conditions required to reproduce observations of Venus following the modeling of Krasnopolsky (2012), consisting of fixed  mixing ratios at 10 bar of 30 ppm H₂O, 4.5 ppm H₂, 5.5 ppm NO, 28 ppm SO₂, 1 ppm OCS, and 400 ppm HCl. We calculate the photochemical–kinetic equilibrium profiles of these and other trace gases subject to these boundary conditions. These trace gases absorb at wavelengths that are complementary to CO₂ and can have a substantial effect on climate (Bullock & Grinspoon 2013; Lee et al. 2016).

For purposes of climate–photochemical modeling for observational discriminants, we separately model both clear-sky and cloudy Venus-like atmospheres, except for planet b, which does not condense sulfuric acid aerosols in our photochemistry model. Both clear-sky and cloudy atmospheres display a similar, Venus-like profile (dashed black line in Figure 6), despite being irradiated by a very different SED than the Sun. The temperature profiles are characterized by long adiabats from the surface each to a ∼0.05–1.0 bar tropopause. The surface temperatures are the hottest of the environments we simulated, ranging from 259–465 K (200–398 K cloudy) for planet h (outside the outer edge of the HZ) to 714–927 K for b, which receives approximately twice the stellar flux of Venus.

As in the O₂-outgassing cases, the Venus-like atmospheric photochemistry is dominated by the formation of sulfuric acid from the combination of water with SO₃, a photolytically derived product of SO₂. This sulfur chemistry is not driven as strongly as that of Venus due to lower MUV–NUV flux from TRAPPIST-1, so the rates of both SO₂ and SO₃ photolysis are much lower than for Venus. However, SO₂ is still susceptible to photolytic destruction, which results in the trend seen in Figure 6, that SO₂ is more abundant farther from the star. The abundance of SO₃ is a balance between its formation rate via SO₂ photolysis and the rate at which it combines with water to form H₂SO₄. In our simulations, water vapor decreases as a function of semimajor axis in both clear-sky and cloudy cases due to this condensation into H₂SO₄ aerosols, as well as into water and water-ice in the cooler planets (e, f, g, and h, although we do not model the water clouds that would form as a result of this process).

The prevalence of H₂SO₄ vapor depends on its chemical formation, where and whether it condenses, and its destruction via thermal decomposition in the lower atmosphere. In our simulations, H₂SO₄ formed in abundances similar to that on Venus (up to 17 ppm versus 8.5 ppm for Venus; Krasnopolsky 2015) for the hotter TRAPPIST-1 planets and drops in abundance with semimajor axis due to the reduced availability of H₂O at lower temperatures. The highest abundances of H₂SO₄ and SO₃ are seen in the stratospheres of the hotter planets (b and c), with peak vapor concentrations decreasing with semimajor axis and at lower altitudes. Similarly, the peak optical depths of the clouds drop with semimajor axis due to less effective formation processes. Since H₂SO₄ can persist closer to the surface in the colder atmospheres, the semimajor-axis-dependent total optical depth peaks for planet d, due to the competing effects of formation/condensation and thermal decomposition/rainout (see Figure 7). Notably, for planet b, the modeled atmosphere was too hot for condensation. Planet d exhibited the highest aerosol optical depth due to its low surface gravity, which allows aerosols to be more easily lofted in the atmosphere.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The sulfuric acid aerosols in our cloudy Venus atmospheres substantially reduced the surface temperatures. The net top of atmosphere stellar irradiance was reduced by ∼50%–60%, resulting in temperature ranges of 200–616 K (10 bar) and 398–779 K (92 bar), h to c. As a caveat, it is possible that sulfuric acid production in Venus's atmosphere may vary by orders of magnitude (Gao et al. 2014), and it is unknown how an alien M dwarf spectrum may affect the more uncertain aspects of Venus-like atmospheric chemistry, particularly of the unknown UV absorber, which we exclude.

3.2.4. Aqua Planet

Here we show simulated transit transmission and emission spectra for our evolved worlds to support observation planning with JWST for this planetary system. These are noiseless spectra generated using SMART, sampled at 1 cm⁻¹ and convolved with a 1 cm⁻¹ half-width at half-maximum slit function, which can be used as the model input for instrument and observation simulators that calculate noise sources and instrument sensitivity. These spectra can also be compared with known or anticipated instrument noise floors to determine zeroth-order signal detectability. In a subsequent paper (J. Lustig-Yaeger et al., in preparation), these spectra will be used to quantify the detectability of our simulated atmospheres with JWST and assess optimum observing techniques for identifying key planetary characteristics.

We have also produced direct-imaging reflectance spectra. Reflectance spectra of late-type M dwarfs may first be observed by large ground-based instruments, or later by HabEx/LUVOIR. The spectra we simulated, including those not presented here, are available online using the VPL Spectral Explorer¹³ or upon request.

3.3.1. Transit Transmission Spectra

In Figures 9–11, we show transit transmission spectra of our modeled TRAPPIST-1 planetary environments, covering nearly the entire wavelength range of potential JWST exoplanet observations (0.5–20 μm). These spectra are presented as "relative transit depth," that is, the signal originating from the atmosphere as compared to the solid planetary body. Note in this paper we do not account for stellar limb darkening, so "transit depth" here refers to (R_p/R_*)², where R_p is the radius of the solid body and R_* the radius of the star. To assess the signal originating from the atmosphere, we can expand the transit depth,

$$\begin{eqnarray}&&\displaystyle \frac{{dF}}{F}={\left(\displaystyle \frac{{R}_{p}+{R}_{a}}{{R}_{* }}\right)}^{2}={\left(\displaystyle \frac{{R}_{p}}{{R}_{* }}\right)}^{2}+\displaystyle \frac{2{R}_{p}{R}_{a}}{{R}_{* }^{2}}+{\left(\displaystyle \frac{{R}_{a}}{{R}_{* }}\right)}^{2},\end{eqnarray} \tag{ 3 }$$

and define the relative transit depth of the atmosphere, valid for atmospheres with small scale height:

$$\begin{eqnarray}&&\displaystyle \frac{{{dF}}_{a}}{F}\approx \displaystyle \frac{2{R}_{p}{R}_{a}}{{R}_{* }^{2}},\end{eqnarray} \tag{ 4 }$$

where R_a is the vertical extent of the atmosphere. Although the radius of the solid body and altitude of the atmosphere for exoplanets cannot be separately measured, this equation allows us to relate predicted variations in transit depth of the modeled spectra for different-sized bodies to the physical extent of the modeled atmospheres.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The seven planets exhibit a variety of features in transmission, the strengths of which depend on each planet's characteristics. TRAPPIST-1 b displays the strongest features, owing to its higher temperature and moderately low surface gravity, both of which increase the atmospheric scale height and therefore its transmission signal. Planet b shows molecular signals exceeding 200 ppm, which is well above the putative 20, 30, and 50 ppm noise floors for JWST NIRISS SOSS, NIRCam grism, and MIRI LRS, respectively (Greene et al. 2016). Planet d exhibits molecular signals nearly as strong as b, in large part because d has the lowest surface gravity of these planets. Due to geometry and the effects of refraction  (Misra et al. 2014; Meadows et al. 2018a), the surface and near-surface atmospheres of the 10–100 bar atmospheres cannot be probed with JWST transmission spectroscopy. All of the modeled atmospheres we present include CO₂, and despite a broad range of abundances, they all exhibit strong CO₂ absorption at 2.0, 2.8, 4.3, and 15 μm, though the hotter or higher abundance CO₂ atmospheres have additional CO₂ features, particularly shortward of 2 μm. In the following subsections, we present the spectra for each environment.

O₂-Dominated Atmospheres, Desiccated: Transit transmission spectra for the desiccated, high-O₂ atmospheres are presented in Figure 9. As in Schwieterman et al. (2016) and Meadows et al. (2018a), the characteristic spectral features of a high-O₂ atmosphere are strong O₃ and O₂–O₂ CIA bands together. The CIA bands are strongest at 1.06 and 1.27 μm. A mixture of O₂, O₃, and O₂–O₂ is present at 0.5–1.3 μm, and a mixture of CO₂, CO, and O₃ at 1.5–6 μm. The strongest O₃ features are the Chappuis band at ∼0.6 μm and the 9.6 μm feature, with weaker bands at 2.5, 3.3, 3.6, 4.75, and 5.6 μm. These weaker bands are generally only available here because of high levels of O₃ combined with relatively low levels of CO₂. Most of this wavelength range could be accessible to JWST NIRISS and NIRSpec instruments. Features of CO₂ and O₃ are also present at 6–20 μm, which could be probed by the JWST MRS instrument. In addition to the complete lack of any water features, the presence of CO features (2.35 μm and 4.6 μm) is characteristic of this severely desiccated atmosphere (e.g., Gao et al. 2015), though these will be very challenging to observe with JWST.

O₂-Dominated Atmospheres, Outgassing: The transmission spectra of the outgassing atmospheres are shown in Figure 9. The detectability of outgassed trace species varies by planet, particularly for H₂O and SO₂. Water features at 1.4, 1.9, 2.6, 3.0–3.3, and 6.3 μm rival the strength of CO₂ bands for the hottest planets (b, c, and d), due to stratospheric water vapor in excess of the moist greenhouse limit, indicative of planets experiencing water loss. The colder planets have insufficient H₂O to generate H₂SO₄, so they maintain higher levels of SO₂, which shows small features of less than 60 ppm at 4.0, 7.3, 8.7, and 19 μm. Much larger SO₂ fluxes or lower water vapor abundances would be required to generate SO₂ features above the putative noise floor of JWST (20–50 ppm; Greene et al. 2016).

Venus-like Atmospheres: Clear-sky and cloudy Venus-like transmission spectra are shown in Figure 10. With a clear sky, Venus-like atmospheres exhibit deep absorption features, primarily from CO₂. Water vapor is present in the clear-sky spectra at 0.9–2.0 μm and at 6.3 μm. Sulfuric acid aerosols truncate the minimum altitude probed by the transmission measurement, severely reducing the strength of the absorption features. These Venus-like planets may have absorption features in transmission approaching ∼90 ppm.

Aqua Planet: Clear-sky and cloudy aqua planet transmission spectra for TRAPPIST-1 e are shown in Figure 11. Due to the large angular size of TRAPPIST-1 as seen from the planet, the clear-sky atmosphere can be probed nearly to the surface, while the cloud deck  truncates the spectrum at 11 km. A number of molecular absorption bands are present, including O₃ (∼0.6, 4.7, 9.6 μm), O₂ (0.76, 1.27 μm), CH₄ weakly (3.3 and 7.7 μm), and H₂O (6.3 μm). Even with an abundance of only ∼290 ppm, CO₂ is the most observable molecule in these atmospheres (1.6, 2.0, 2.7, 4.3, 15 μm). The bulk atmospheric constituent, nitrogen, exhibits a 4.1 μm CIA band, which could provide a sensitive probe of its partial pressure (Schwieterman et al. 2015).

3.3.2. Emission Spectra

We present thermal emission spectra as normalized flux (Figures 12 and 13). Below ∼5 μm, stellar flux reflected from the surface or scattered in the atmosphere contributes significantly to the spectra.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

O₂-Dominated Atmospheres, Desiccated: The desiccated, O₂-dominated atmospheres (Figure 12) are distinctive in emission spectra, due to lower temperatures, lack of water bands, and their unusual temperature structures. The stratospheric temperature peak causes some bands to be seen in emission rather than absorption, depending on how cool the particular planet is. In the colder atmospheres, the O₃ band at 9.6 μm exhibits varying hot wing emission at 10.4 μm from the stratospheric temperature peak. In the hotter atmospheres with lower O₃ abundance, part of the 9.4 μm CO₂ "hot band" also emerges in emission at 9.0 μm. CO₂ varies between absorption and emission at ∼10–18 μm. As a potential insight into isotope fractionation for ${}^{18}{\rm{O}}{/}^{16}{\rm{O}}$, the fundamental asymmetric stretch ${v}_{1}=1\to 0$ transition of the ${}^{16}{{\rm{O}}}^{12}{{\rm{C}}}^{18}{\rm{O}}$ isotopologue at 7.3 and 7.9 μm (Rothman et al. 2013) appears in emission. The O₂ band at 6.4 μm exhibits strong emission, but it is very narrow and is present only because of the complete lack of water vapor. Carbon monoxide is not observable.

O₂-Dominated Atmospheres, Outgassing: The emission spectra of the O₂-dominated, outgassing atmospheres are given in Figure 12. The radiative–convective processes and resultant emission spectra are largely dominated by water vapor in the atmosphere. The warmer planets, with larger quantities of water vapor throughout the atmospheric column, exhibit massive H₂O absorption in the 6.3 μm band. Conversely, the colder, drier planets exhibit SO₂ absorption at 7.3, 8.8, and 19 μm. These outgassing atmospheres have much less ozone than the desiccated ones, due to destruction catalyzed by the presence of water vapor, so their O₃ features are much weaker, particularly for the warmer planets.

Venus-like Atmospheres: The Venus-like atmospheres (Figure 12) emit very little flux, except through narrow windows between strong absorption bands, similar to Venus. This suppression of emitted flux is due to the combined Venus greenhouse of primarily CO₂, H₂O, and SO₂. The emission windows in the hotter atmospheres are at 2.4, 3.5, and 5.5–7.0 μm. The cooler, clear-sky planets emit some flux around 10 μm, between CO₂ hot bands. In the clear-sky case, water vapor absorbs weakly in the 5.5–7.0 μm window, and SO₂ absorbs at 7–9 μm, but these are obscured in the cloudy case. The remaining emission is reduced by CO₂ absorption. The cloudy spectra exhibit more shallow absorption features and are more similar to a blackbody spectrum because the thermal emission is primarily from the cloud deck, though the colder planets also emit above the cloud deck from the 15 μm CO₂ band.

Aqua Planets: The potentially habitable aqua planet atmospheres (Figure 13) have emission spectra very similar to that of Earth (see Robinson et al. 2011). The near-surface absorption from the water vapor continuum is apparent in the clear-sky case, particularly at 10–13 μm, and the cloudy case emits from its global cloud deck at a colder temperature. Both cases exhibit absorption from H₂O (6.3 μm), O₃ (9.6 μm), and CO₂ (15 μm).

We conducted a study of H₂O loss and O₂ accumulation for the TRAPPIST-1 system using an update to the model of Luger & Barnes (2015), and we found that ocean loss affects the entire HZ. Even planets that are beyond the outer edge during the main sequence phase (e.g., TRAPPIST-1 h) have the potential to lose a significant ocean fraction and generate O₂. These results are broadly consistent with a study on ultracool dwarf stars conducted by Bolmont et al. (2017). Note that Bolmont et al. based their TRAPPIST-1 work on earlier system parameters (Gillon et al. 2016), which included only planets b, c, and a range of semimajor axes for d. Where we can compare, we predict up to an order of magnitude larger amounts of H loss and O₂ accumulation, likely due to our much older stellar age (5 Gyr, Burgasser & Mamajek 2017, versus 850 Myr, Bolmont et al. 2017) and assumptions for the planetary radii and masses, which we based on additional observations (Grimm et al. 2018).

Our model and calculation of the escape efficiency _xuv are consistent with the value inferred for low-mass planets in the Kepler population, considering their sizes and densities as a function of irradiation (Owen & Wu 2013). Another metric, the "cosmic shoreline," supports our use of thermally driven hydrodynamic escape. The "cosmic shoreline" is an empirically derived relationship between planetary irradiation and escape velocity that divides planets with atmospheres from airless rocks and appears to apply broadly to the solar system and the population of exoplanets with known densities (Zahnle & Catling 2017).

A thorough discussion of caveats for assessing H escape and O₂ accumulation was given by Bolmont et al. (2017). We do not include the effects of O₂ loss processes other than thermal escape, such as surface  sinks (Schaefer et al. 2016; Wordsworth et al. 2018) and nonthermal atmospheric loss (Collinson et al. 2016; Airapetian et al. 2017),  which are important and warrant further detailed study. At the location of Proxima Centauri b for an Earth twin, the calculated nonthermal escape rate may not currently exceed ;Earth-like outgassing replenishment rates, which may allow a secondary atmosphere to accumulate (Garcia-Sage et al. 2017). However, nonthermal escape could have contributed to net O₂ loss during the pre-main-sequence phase. Therefore, the O₂ generated by our evolution models represents maximum possible inventories, using nominal parameters for stellar and planetary properties. We reduced the O₂ abundance for our climate–photochemistry and spectral models in consideration of these loss mechanisms, which we did not include in our models. Given the importance of oxygen as a biosignature, the generation of abiotic, O₂-rich post-ocean-loss atmospheres should be considered when interpreting future observations of M dwarf terrestrial planets.

These diverse environments produce a wide range of habitable and uninhabitable surface temperatures (Figure 4), suggesting that the region around a star where Venus-like planets may be common evolutionary outcomes (Kane et al. 2014) may extend into and beyond the classical habitable zone for late-type M dwarfs. Even though some of these modeled environments do not support surface liquid water, if dense CO₂ atmospheres do form on terrestrial planets after an early runaway greenhouse phase, our results suggest that surface temperatures in the habitable range or hotter may occur beyond the classical maximum greenhouse limit of the HZ outer edge (where warming from additional CO₂ is offset by higher Rayleigh scattering in a temperate H₂O–CO₂ greenhouse; Kasting et al. 1993). Even TRAPPIST-1 h, beyond the HZ outer edge (${S}_{{\rm{h}}}=0.148{S}_{\oplus }$), could have a surface temperature as hot as 465 K if early atmospheric evolution resulted in a clear-sky Venus-like atmosphere, and 398 K with the expected sulfuric acid clouds. These high temperatures may be stable over long time periods, because known sinks for CO₂—such as deposition to an ocean and carbonate–silicate weathering (Walker et al. 1981)—will not be available, although a surface carbonate buffer may form and keep the atmosphere at an equilibrium CO₂ value that is less than 92 bar (Hashimoto et al. 1997). Figure 4 also shows that across the diversity of atmospheres considered, the instellation at the modeled TRAPPIST-1 e's distance from the star was most likely to maintain a habitable surface temperature. This suggests that the habitable zone concept is still useful as an initial means of narrowing down the search for habitable conditions on terrestrial planets, even when planetary atmospheres are not truly Earth-like.

For late-type M dwarf stars in particular, our work indicates that the maximum greenhouse limit may not apply, or may apply at a distance much farther from the star than previously calculated (Kopparapu et al. 2013). We have shown that the H₂O–CO₂ greenhouse used in traditional habitability calculations can be enhanced significantly by other greenhouse gases, including SO₂. While SO₂ is a greenhouse gas in desiccated atmospheres like that of Venus (e.g., Lee et al. 2016), it is less likely to be present in the atmosphere of a habitable planet with an ocean because it is highly soluble. If an initially desiccated planet beyond the HZ were to recover habitability via outgassing and maintain a warm temperature, other gases would be needed to support a strong greenhouse. Greenhouse gas abundances can also be affected by energetic particle events, which can drive chemistry that generates trace amounts of greenhouse gases (Airapetian et al. 2016). These events are likely to remove greenhouse gases, by the destruction of O₃ in Earth-like atmospheres by repeated flaring (Tilley et al. 2017) or by the destruction of CH₄, which may form hydrocarbon hazes if oxygen radicals are not present (Arney et al. 2017). In the absence of sufficient greenhouse gases, the coldest atmospheres (particularly the 10 bar cloudy, Venus-like h) may be subject to collapse, but most will be resistant, as the condensation temperature for O₂ is very low and day–night temperature contrasts are lower in thicker atmospheres (Turbet et al. 2018).

Our results indicate that the interplay between stellar irradiation, photochemistry, chemical kinetics, and condensation shaped our simulated planetary atmospheric compositions, with water vapor, ozone, sulfur species, and aerosols displaying the greatest variation. The decreases in H₂O (largely due to temperature) and increases in O₃ abundance (due to an interplay between photolysis rates, water abundance, and vertical transport) as a function of orbital distance are consistent with trends found by previous studies (e.g., Grenfell et al. 2007). The water vapor profile in the warmer atmospheres is representative of planets in a moist greenhouse state, with substantial water vapor well mixed throughout the atmosphere and with no cold trap. Consequently, these atmospheres may suffer severe water loss over gigayears, unless the escape rate is balanced by outgassing. Unless they begin with warm, Venus-like atmospheres, the outer planets are too cold to recover from freezing temperatures without a brightening star or larger inventories of greenhouse gases, which may eventually accumulate from outgassing or impactor delivery.

Ozone and CO both exhibit characteristic patterns across our modeled atmospheres. In those containing ≳0.2 bar O₂, at ∼200 Pa, O₃ is about 10⁻⁵ mol/mol in all of the simulated atmospheres, despite different levels of irradiation, H₂O, and O₂. Above this layer, the O₃ abundance is controlled by photolysis, while below, mixing and other temperature- and water-dependent loss mechanisms dominate. In desiccated atmospheres, CO abundance is controlled by the slow three-body CO₂ recombination reaction (Gao et al. 2015), while in atmospheres with water vapor, the abundance of CO is controlled by the availability of H-bearing recombination catalysts, which are more prevalent near the planetary surface, removing the CO at lower altitudes. Large amounts of CO in the presence of O₂/O₃ are a false positive discriminant (Schwieterman et al. 2016), indicating that oxygen is more likely to have been generated by photolysis of CO₂ (Gao et al. 2015) rather than photosynthesis.

Aerosols affect the planetary climates by reflecting incoming stellar radiation and either absorbing, emitting, or scattering thermal radiation. Sulfuric acid aerosols can reduce the surface temperatures of hot, partly desiccated planets, which caused the cloudy Venus-like atmospheres we modeled to drop by up to nearly 200 K (Figure 4) compared to the clear-sky cases. Sulfuric acid production also varied as a function of irradiation. The warmer planets more easily formed H₂SO₄ due to higher water vapor abundances and more rapid SO₂ photolysis. This production was offset at lower altitudes by thermal decomposition. TRAPPIST-1 b did not condense sulfuric acid in our model because atmospheric temperatures were too high. If any of the TRAPPIST-1 planets have Venus-like atmospheres, we expect that b will not form sulfuric acid clouds, but the other planets likely would. Although we did not model the radiative effects of aerosols for the O₂ outgassing environments, sulfuric acid condensation occurred in the photochemical model, so their climates could also be affected by these clouds. Water and water-ice clouds impacted the temperature structure of the aqua planet, but had little effect on surface temperature (Figure 8). This was due to competing effects as stellar radiation interacted with the cloud properties: the upper tropospheric cirrus clouds effectively scatter stellar radiation, but the absorption properties of water and water-ice in the NIR also make them effective absorbers of M dwarf peak stellar radiation, especially given the spectral energy distribution of TRAPPIST-1. Exoplanet aerosols for terrestrial planets are a rich field for future study.

The parameter space for possible planetary characteristics is large, so the assumptions regarding atmospheric surface pressure, surface material, atmospheric bulk composition, geological or biological surface fluxes, and atmospheric transport can have a large impact on the equilibrium atmospheric state. The differences between our O₂ desiccated and outgassing atmospheres demonstrate the effects of surface fluxes on the planetary environment. These factors are important when considering what atmospheres may be possible on a given exoplanet, and for comparing possible planetary states and among different models. Wolf (2017) found that TRAPPIST-1 e requires more CO₂ than Earth to sustain temperate conditions, but that work underestimated the reduction in bolometric albedo due to the late M dwarf SED (Turbet et al. 2018). It is possible that Turbet et al. (2018) also overestimated the surface albedos of ocean-dominated TRAPPIST-1 planets, as they calculated a mean water-ice albedo of 0.21. We used wavelength-dependent surface albedos in our climate calculations that depend on the planetary environment (Figure 1). We found that both pure ocean and melting ice surfaces had very low stellar flux-weighted average albedos (0.02 and 0.03 respectively), while a typical snow surface averaged only 0.16 around TRAPPIST-1. Wolf (2017) found that 0.1 bar of CO₂ with 1 bar N₂ is required for a surface temperature of ∼270 K, while modern Earth CO₂ levels of ∼400 ppm yield 240 K. Turbet et al. (2018) calculated 250 K under similar conditions, while we find that with or without clouds, only preindustrial Earth CO₂ abundance of 280 ppm is required in a 1 bar N₂/O₂ atmosphere on an ocean-covered planet to reach an Earth-like globally averaged surface temperature. It is crucial to note that surface albedo is one of the primary drivers of climate state (Godolt et al. 2016), particularly in an atmosphere transparent to its host star's radiation.

Here we compare spectra derived from the complete set of simulated atmospheres for TRAPPIST-1 b, c, and d, whose higher temperatures make them more easily observable, and e, which is a potentially habitable candidate, to identify observational characteristics that could be used to discriminate between these environments and their evolutionary histories. Spectral features that can be used to discriminate between the modeled environments are given in Table 4. Figure 14 shows spectra for the atmospheres we modeled for TRAPPIST-1 b, c, and d. The stronger carbon dioxide bands are prominent in all of our planetary spectra, regardless of whether CO₂ is a bulk constituent or a trace gas, and these CO₂ bands are therefore a good indicator of the presence of a terrestrial atmosphere, but are not as useful for discriminating among different environments. However, if the weaker CO₂ bands are also present, this does indicate a higher, more Venus-like abundance. As seen in planet d's transmission spectra in Figure 14, cloudy, Venus-like atmospheres are well discriminated by sulfuric acid absorption longward of 3 μm (though these clouds could also form in O₂-dominated atmospheres with outgassing).

Zoom In Zoom Out Reset image size

The massive O₂ atmospheres may reveal themselves in transmission via the presence of O₂–O₂ collision-induced absorption at 1.06 and 1.27 μm (Schwieterman et al. 2016), a wavelength region strongly backlit by the host star. Note that our photochemical model suggests that cloud formation in the outgassing atmosphere is likely due to sulfuric acid condensation (at ∼1 ppb concentration). However, this will not affect the visibility of the O₂–O₂ features since the aerosol scattering properties and geometry of the TRAPPIST-1 system permit transmission shortward of 2.5 μm for the innermost planets. The favorably compact geometry of the TRAPPIST-1 system allows transit transmission to probe relatively deep into the atmosphere, where O₂–O₂ absorption is strongest, producing transit signals as high as 70 ppm (versus ∼3 ppm for an early-type M dwarf; Schwieterman et al. 2016). However, as demonstrated in the inset in Figure 14, the broad O₂–O₂ bands at 1.06 and 1.27 μm overlap with CO₂ bands at 1.05, 1.13, 1.24, 1.29, and 1.32 μm in high CO₂ atmospheres. An atmosphere with high levels of both CO₂ and O₂ may be common if a planet that experienced complete ocean loss retained a large fraction of the CO₂ that was dissolved in the oceans or subsequently outgassed from the interior. These weak CO₂ features could confuse retrieval and interpretation of the O₂–O₂ features without sufficient resolution, signal-to-noise ratio, or corroboration from stronger CO₂ features elsewhere in the spectrum.

These O₂-dominated spectra also show strong features from O₃ and SO₂. A strong ozone Chappuis band around ∼0.6 μm may indicate desiccation, as it traces ozone deep in the atmosphere, where water vapor would normally destroy it (Meadows et al. 2018a). As clearly shown in Figure 14 in all three panels, the 9.6 μm O₃ feature is strong and broad for both O₂-rich atmospheres, but not present in the Venus-like spectra primarily due to the lack of O₂ in the atmosphere. Detection of SO₂ absorption in either transmission or emission would be indicative of an outgassing interior, and possibly an atmosphere and surface with low water abundance, due either to early desiccation or temperatures below freezing. The O₂ outgassing case is unique among the environments considered for the hotter planets in showing very strong water features, especially near 1.4 and 6.3 μm. This is because these planets have high levels of water vapor (up to 0.1% mol/mol) in their stratospheres, so they are in a moist greenhouse state, approaching a runaway greenhouse. These strong H₂O bands are indicative of a planet in the process of losing water vapor, rather than indicating habitable, ocean-bearing conditions.

The bottom panel of Figure 14 demonstrates the differences between the evolved environments in emission using TRAPPIST-1 c as an example of the hotter planets. The water band at 6.3 μm is the primary discriminant. At high temperatures and high abundance (e.g., the O₂ outgassing case), H₂O absorption can dominate over absorption from other interesting molecules at 5–9 μm, such as O₂, SO₂, and CO₂. For the Venus-like, clear-sky atmosphere, the low abundance of H₂O ≲ 1 ppm fails to suppress the outgoing radiation from the hot lower atmosphere, resulting in an emission peak at these wavelengths. The cloudy Venus case shows that clouds can adversely affect both outgoing flux and the strength of molecular features by truncating the atmosphere at higher, colder altitudes.

Figure 15 compares the spectra of the O₂- and CO₂-dominated atmospheres with the clear and cloudy aqua planets for TRAPPIST-1 e. The relative transit depths are naturally offset by the truncating effects of clouds, which allow transmission to probe only altitudes above the cloud deck, and the obscuring effect of atmospheric opacity for atmospheres of different total atmospheric pressure (i.e., 1 versus 10 bar). In the thermal infrared, clouds reduce the overall outgoing flux by emitting at colder temperatures. The O₂ desiccated and Venus-like atmospheres for TRAPPIST-1 e display features and corresponding observational discriminants similar to the hotter planets. However, the O₂ outgassing environment was not hot enough to maintain high levels of stratospheric H₂O, so it exhibits a more Earth-like, modestly absorbing 6.3 μm feature. The overall emission spectrum of that environment is very similar to that of aqua planets, but with additional CO₂ bands due to the higher levels of CO₂ we assumed.

Zoom In Zoom Out Reset image size

Water is not a good discriminant among some of these atmospheres, so discriminating the habitable aqua planet from the uninhabitable environments is not straightforward. The aqua planet has water vapor features in the near-infrared in transmission at 0.94, 1.15, 1.4, 1.9, 2.6, 3.7, and most strongly at 6.3 μm. Unfortunately, these water features have comparable signals among the clear-sky atmospheres (except the completely desiccated case), and among the cloudy cases, despite vastly different levels of tropospheric H₂O, due to the cold trapping of water in the aqua planet atmosphere.

The bulk composition of these atmospheres can potentially be distinguished by weak spectral effects of their dominant gases O₂, CO₂, and N₂. The Earth-like aqua planet atmosphere exhibits absorption from the collision-induced N₂–N₂ band at 4.1 μm, which is sensitive to the partial pressure, and could be diagnostic of atmospheric pressure if this challenging observation could be made (Schwieterman et al. 2015). At high abundances of CO₂ (as low as 0.5 bar here), the wings of the 4.3 μm CO₂ band overpower the N₂–N₂ CIA. The presence of the weak NIR–MIR CO₂ bands  (including the numerous bands shortward of 2 μm, the ${}^{16}{{\rm{O}}}^{12}{{\rm{C}}}^{18}{\rm{O}}$ isotopologue absorption bands at 3.6, 3.8, 4.0, 7.3, and 7.9 μm, and hot bands at 9.5 and 10.5 μm) is indicative of a warm, CO₂-rich atmosphere. As in Earth's atmosphere, O₂–O₂ CIA is weakly present in the aqua planet atmospheres, but the bands are considerably stronger at higher O₂ abundance, which would suggest that an abiotic, post-ocean-loss source for the oxygen was more likely.

Potentially habitable environments may be constrained by indicators of desiccation, such as O₃, SO₂, and CO, as these gases react with water vapor and its photolytic by-products. Both the modeled O₂-dominated and aqua planet atmospheres exhibit ozone absorption at 0.6 and 9.6 μm. However, a strong indication of desiccation is high ozone abundance, as constrained by the presence of a long wavelength tail of the Chappuis band (out to 1 μm), the 2.5 and 3.6 μm features, and additional weak bands at 9 and 13 μm, which arise in the O₂ cases, particularly the colder and desiccated environments. Retrievals of ozone abundances in conjunction with climate and photochemical modeling could help constrain the likely abundances of tropospheric water vapor. The detection of SO₂ would also suggest a desiccated atmosphere, because SO₂ easily combines with H₂O to form sulfuric acid. The presence of CO, if not masked by CO₂ or O₃, would indicate desiccation because H₂O photolytic by-products such as hydroxyl are efficient catalysts at recombining CO₂ from CO and free oxygen. While the absence of water vapor absorption does not necessarily indicate the lack of tropospheric water vapor or surface water, desiccation indicators from O₃, SO₂, and CO, even with weak water vapor features, would suggest a desiccated and likely uninhabitable environment.

Signs of atmospheric escape and ocean loss may be probed by observing isotopologue bands HDO and ${}^{18}{{\rm{O}}}^{12}{{\rm{C}}}^{16}{\rm{O}}$. Measuring the D/H ratio in H₂O and the ¹⁸O/¹⁶O ratio in CO₂ could be used to probe isotopic fractionation, which occurs during atmospheric escape and ocean loss because lighter elements are preferentially lost over heavy elements (e.g., Hunten 1982). D/H fractionation from atmospheric loss has been measured for both Venus (D/H ∼ 120 versus Earth; De Bergh et al. 1991) and Mars (D/H ∼ 4 versus Earth; Encrenaz et al. 2018). In particular, such D/H fractionation would indicate an environment that had the majority of its ocean reservoir depleted (Hunten 1982). Future work should assess oxygen fractionation due to massive ocean loss in an environment with some retained atmospheric oxygen.

Our results and those of others demonstrate that the atmospheres and surfaces of M dwarf terrestrials have likely been strongly shaped by several factors, including stellar evolution, atmosphere and ocean loss, photochemistry, and geological (and biological) surface fluxes (Segura et al. 2005; Grenfell et al. 2014; Rugheimer et al. 2015; Meadows & Barnes 2018; Meadows et al. 2018a). In particular, interior outgassing and photochemistry are key processes for modeling terrestrial exoplanets, as both drive disequilibria in the atmosphere. Our results show a strong effect from outgassing, such as from SO₂ and H₂O, which, together with their photochemical by-products and reactions with other species, shape the temperature inversions and lapse rates as well as produce radiative greenhouse (and antigreenhouse) effects. Although our environments do not explicitly contain biogenic fluxes, surface fluxes from life processes also drive disequilibria, and this is one of the principal means that will be used to search for life on exoplanets (Simoncini et al. 2013; Krissansen-Totton et al. 2016; Meadows 2017). Photochemistry is critically important in assessing planetary composition and in understanding the formation and destruction of key species, some of which may be outgassed or biogenic. As we have shown, this is especially true for planets at different orbital distances, including those orbiting M dwarfs, where the alien UV spectrum can drive chemistry that the Sun does not induce in the Earth's atmosphere. The species outgassed or generated from photochemistry may be directly observed, and their vertical distribution may affect other gas absorption features because of their effects on climate.

To illustrate the effects of modeling factors—including outgassing and photochemistry—on predicted spectra, we directly compare a sampling of our spectra for TRAPPIST-1 d and e with Morley et al. (2017), who used radiative transfer and thermochemical equilibrium models to conduct a broad assessment of spectroscopic observables for the TRAPPIST-1 planets. For planet d, we compare our radiative–convective/photochemical equilibrium 10 bar clear-sky and cloudy Venus-like atmospheres with the Venus-like atmospheres from Morley et al., identified by their albedos (A = 0.0, 0.3, and 0.7), which correspond to clear-sky to global cloud coverage (Figure 16). The spectra of clear-sky, Venus-like atmospheres modeled in both approaches are similar in both transmission and emission, except for stronger SO₂ features from the thermochemical model, which are reduced in our model (and Venus's true atmosphere) by H₂SO₄ formation. The SO₂ in the thermochemical model suppresses outgoing flux at wavelengths of 7–10 μm when compared to both our clear and cloudy atmospheres. Although the integrated flux of the A = 0.0 atmosphere is very close to that of our clear-sky Venus, the distribution of flux is different. Our absorption features are not as deep (e.g., 7–9 μm), which is likely due to our atmospheres having warmer stratospheres. Our atmospheres also show higher emission near 6 μm because weaker water vapor absorption permits escape of radiation from the hot lower atmosphere. This is likely due to different assumptions about the water vapor profile between the two models.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our spectra differ considerably from Morley et al. (2017) for cloudy Venus-like cases. Photochemical models produce results tailored to the spectrum of the host star, whereas thermochemistry is driven purely by the atmospheric temperature profile at a given orbital position and so is more generic. The cloudy Venus differs significantly from the clear-sky case in both transmission and emission due to the truncation of the atmospheric column by cloud scattering and absorption at higher, colder altitudes, which reduce the outgoing thermal flux and can raise the tangent height when observed in transmission. Although the effects of clouds on surface temperature can be approximated by changing the Bond albedo in the thermochemical models, this process does not take into account the radiative effect of the clouds on the temperature profile or on transmission and emission spectra. Sulfuric acid scattering and absorption in the cloudy Venus-like transmission spectrum, and emission from the cloud deck in direct imaging or secondary eclipse, are distinctive features of a Venus-like planet not captured in models without vertical cloud distributions.

As shown in Figure 17, the Earth-like atmospheres also show significant differences in composition between the two modeling approaches, primarily because we include the effects of water vapor clouds, photochemically produced gases such as O₃, and outgassed species such as CH₄. These characteristic features of Earth's spectrum are not seen in the spectra generated by Morley et al. (2017) using a thermochemical equilibrium model. Our photochemically consistent clear-sky aqua planet atmosphere also produces significantly stronger water vapor absorption bands (likely due to the inclusion of HDO isotopologue absorption at 3.6 and 9 μm), although the cloudy case is perhaps the more realistic model.

Several preliminary observational constraints on atmospheric composition have been obtained of planets in the TRAPPIST-1 system, to which we can compare our modeling results. Spectra taken with the Hubble Space Telescope (de Wit et al. 2016, 2018) ruled out low-molecular-weight atmospheres, although these spectra were also consistent with flat spectra within the error bars of ∼200–350 ppm. For TRAPPIST-1 b, Delrez et al. (2018) observed a signal difference of 208 ± 110 ppm between the 3.6 and 4.5 μm Spitzer bands, which they suggest is due to CO₂ absorption. These results, and the K2 photometry from the Gillon et al. (2017) discovery paper, are shown with error bars in Figure 18, along with the band-integrated transit depths for our modeled spectra. We have adjusted the assumed solid body radius in these spectra to fit the observations for the solid body and wavelength-dependent atmospheric absorbing radius. Our simulated spectra are consistent with these data within 2σ error, and the O₂-dominated atmospheres fit within 1σ. The 3.6 μm filter bandpass contains absorption by the CO₂ isotopologue bands at 3.6, 3.8, and 4.0 μm, plus the wings of the main isotopologue 4.3 μm band, whereas the 4.5 μm filter bandpass contains the 4.3 and 4.9 μm CO₂, 4.6 μm CO, and 4.75 μm O₃ bands. Rayleigh scattering and ozone absorption may contribute to the larger transit depth measured by K2. A different distribution of UV flux than assumed in this work could enhance ozone levels, further increasing the transit depth in the K2 band and the Spitzer 4.5 μm band.

Zoom In Zoom Out Reset image size

Our simulated atmospheres are consistent with the observed data, but do not agree with the atmospheric scale height and temperature derived by Delrez et al. (2018) for TRAPPIST-1 b. Delrez et al. assumed that the 208 ppm transit signal between the observed Spitzer bands was due to two scale heights from the CO₂ signal in an Earth-like atmosphere—with mean molecular mass 28 g mol⁻¹, CO₂ as the only absorbing gas, and excluding CH₄ and H₂O, which have transitions in the 3.6 μm band. However, their derived temperature of 1400 K and scale height of 52 km (1800 K and 100 km with water vapor included) are not consistent with the assumption of an Earth-like atmosphere. Our self-consistent modeled environments, which include a case with a greenhouse comparable to that of Venus, have emission temperatures that vary with wavelength (325–575 K; Figure 19) but overall stay close to the equilibrium temperature (∼400 K; Gillon et al. 2017), even though they have surface temperatures spanning 406–714 K.

Zoom In Zoom Out Reset image size

In Figure 20, we show that the scale heights at the effective transit altitude for our modeled atmospheres vary considerably as a function of wavelength and are within values typical of solar system terrestrials (i.e., 8–16 km for Venus and Earth; Meadows & Crisp 1996), spanning ∼9–11 scale heights (3–5 scale heights when convolved with the Spitzer bands), and are consistent with the Spitzer observations. Scale height inferences from the size of potential features are degenerate with temperature, species, abundances, aerosols, and opacities, and so are difficult to accurately assess without properly accounting for a range of atmospheric possibilities  within a retrieval framework (e.g., Line et al. 2013).

Zoom In Zoom Out Reset image size

If the existence of an exoplanetary M dwarf terrestrial atmosphere is confirmed, near-term observations  will provide constraints on the environmental effects of M dwarf stellar evolution on orbiting planets as a function of distance from the star and position in the habitable zone. Our modeling work provides hypotheses for planetary states as a result of evolution and composition that could be tested by these upcoming observations, which may constrain atmospheric temperatures and species abundances using the discriminants we have presented. In a subsequent paper, we will use the modeled planetary states and spectra with instrument noise models to predict observational requirements to distinguish the current composition and evolutionary histories of the TRAPPIST-1 planets.

VPL Climate is based on the line-by-line, multistream, multiple-scattering radiative transfer code SMART (described in Section 2).  VPL Climate solves the 1D surface-atmosphere thermodynamic energy equation as an initial value problem using time stepping. To step toward equilibrium in a model atmosphere without incurring the computational penalty of a full radiative transfer calculation at every time step, VPL Climate uses SMART to generate spectrally resolved Jacobians describing the response of the radiative fluxes to changes in temperature and, where applicable, gas mixing ratios. The Jacobians consist of layer-by-layer, line-by-line stellar and thermal source terms and their derivatives, and derivatives of layer reflectivity, transmissivity, and absorptivity. These quantities are used in a linear flux-adding approach to determine heating rates, as developed by Robinson & Crisp (2018). If, at any time step, the atmospheric properties are outside the linear range of the Jacobians, time stepping is suspended, and updated fluxes and layer radiative properties and their Jacobians are recomputed using SMART. High-resolution fluxes from SMART are spectrally binned to 10 cm⁻¹ (reflected stellar) or 1 cm⁻¹ (emitted thermal) intervals to calculate the stellar heating and thermal cooling in each model layer at each time step.

Spectrally dependent molecular absorption cross sections for rotational–vibrational transitions of gases used by SMART are evaluated using the line-by-line model LBLABC (Meadows & Crisp 1996). LBLABC employs nested spectral grids that fully resolve the narrow cores of individual absorption lines for every atmospheric layer and include the line contributions up to 1000 cm⁻¹ from the line center. Cross sections for CO₂ are computed using sub-Lorentzian wings with experimentally determined χ correction factors to simulate the effects of quantum-mechanical line mixing. Lines for H₂O are computed using a line shape model with super-Lorentzian wings to parameterize the far-wing quasi-continuum absorption due to the finite duration of molecule collisions. Opacities for O₂ are computed using super-Lorentzian line wings (Hirono & Nakazawa 1982) with an exponent of 1.958. Line parameters, collision-induced absorption, and UV–visible cross sections are obtained from a variety of sources (see Section 2.6.2).

Convective processes transport heat vertically in planetary atmospheres. Convection can result in both sensible heating (related to the energy required to change the temperature of an air parcel due to vertical mixing) and latent heating (related to the energy required to change the phase of a substance, e.g., due to condensation and evaporation). These processes can transport substantial amounts of heat from the surface to the upper troposphere, especially for planets with greenhouse gases, which can prevent the deep atmosphere (or surface) from radiating directly to space. For example, early 1D modeling results showed that convective processes reduce Earth's surface temperature relative to pure radiative equilibrium conditions by ∼50 K (Manabe & Strickler 1964). Convective processes are also responsible for transport of trace species in the atmosphere, substantially altering atmospheric chemistry (e.g., Fleming et al. 1999; Charnay et al. 2015).

In VPL Climate, mixing length theory is used to calculate convective heat fluxes and heating rates, as this approach is more physically rigorous and versatile than the convective adjustment typically used in 1D exoplanet models. Convective adjustment sets the troposphere temperature profile to a moist or dry adiabat defined by Earth's atmospheric conditions (e.g., Kopparapu et al. 2013; Godolt et al. 2016). In these regimes, convective adjustment does not allow for atmospheric lapse rates that exceed the adiabat, which occurs on planets with more tenuous atmospheres, like Mars (e.g., Hinson & Wilson 2004). Convective adjustment also stabilizes atmospheric layers perfectly at every time step, even if such instantaneous stability is unphysical (e.g., requires convective motions faster than the speed of sound). Furthermore—and especially problematic for coupling to chemistry models—convective adjustment provides no direct insight into the vertical transport of mass in the atmosphere. In contrast, mixing-length theory uses fundamental physical properties of the atmosphere and its constituents (e.g., air parcel temperature, density, specific heat capacities, and static stability) to estimate the vertical heat transport rate, expressed as a vertical "eddy diffusivity," K. This quantity is then used to calculate the resulting heat and mass fluxes. Mixing-length theory is therefore applicable to a wider variety of atmospheric temperature and composition regimes and can be used to model convection for a diverse range of exoplanet atmospheres.

The convective heating rate is related to the divergence of the convective energy flux:

$$\begin{eqnarray}&&{q}_{c}=-\displaystyle \frac{1}{\rho {c}_{{\rm{p}}}}\displaystyle \frac{\partial {F}_{{\rm{c}}}}{\partial z},\end{eqnarray} \tag{ 5 }$$

where z is altitude, c_p is the specific atmospheric heat capacity, and ρ is the air density. The convective energy flux F_c is given by

$$\begin{eqnarray}&&{F}_{{\rm{c}}}=-\rho {c}_{{\rm{p}}}{K}_{{\rm{h}}}\left(\displaystyle \frac{\partial T}{\partial z}+{{\rm{\Gamma }}}_{\mathrm{ad}}\right),\end{eqnarray} \tag{ 6 }$$

where K_h is the eddy diffusion coefficient for heat, and the adiabatic lapse rate is ${{\rm{\Gamma }}}_{\mathrm{ad}}=g/{c}_{{\rm{p}}}$, where g is the acceleration due to gravity. The atmosphere is unstable to convection when $\partial T/\partial z\gt -{{\rm{\Gamma }}}_{\mathrm{ad}}$. In this case, the eddy diffusion coefficient for heat K_h is given by

$$\begin{eqnarray}&&{K}_{{\rm{h}}}={l}^{2}{\left[-\displaystyle \frac{g}{T}\left(\displaystyle \frac{\partial T}{\partial z}+{{\rm{\Gamma }}}_{\mathrm{ad}}\right)\right]}^{1/2},\end{eqnarray} \tag{ 7 }$$

where l is the mixing length. The nominal mixing length is given by Blackadar (1962):

$$\begin{eqnarray}&&l=\displaystyle \frac{{kz}}{1+{kz}/{l}_{0}}.\end{eqnarray} \tag{ 8 }$$

Here, k is von Kármán's constant and l₀ is the mixing length in the free atmosphere, which we express as

$$\begin{eqnarray}&&{l}_{0}={f}_{{\rm{z}}}H,\end{eqnarray} \tag{ 9 }$$

where H is the pressure scale height, given by

$$\begin{eqnarray}&&H=\displaystyle \frac{{RT}}{\mu g}.\end{eqnarray} \tag{ 10 }$$

Here, R is the ideal gas constant, μ is the mean molar mass, and f_z is the proportionality constant.

The preceding convection parameterization was used by Robinson & Crisp (2018), with ${f}_{{\rm{z}}}=1.0$ and with K_h = 0 under stable conditions. However, turbulent eddies still provide some vertical mass transport in stable conditions, due to shear instability and gravity wave breaking (e.g., Hodges 1969; Kondo et al. 1978; Lindzen & Forbes 1983; Canuto et al. 2008). Based on measurements of Earth and Venus eddy diffusion rates (see Figure 21), we specify the following parameterization for the eddy diffusion coefficient for mass transport:

$$\begin{eqnarray}&&{K}_{{\rm{m}}}={K}_{{\rm{m}},0}{(\rho /{\rho }_{0})}^{-1/2}(1-{e}^{-P/{P}_{0}})+{K}_{{\rm{m}},1}{\rm{\Lambda }}+{K}_{{\rm{m}},2},\end{eqnarray} \tag{ 11 }$$

where Λ is a stability-related parameter that is nonzero only for the aqua planet. To determine Λ, we iterate the following from the surface upward, layer-by-layer:

$$\begin{eqnarray}{{\rm{\Lambda }}}_{k}=\left\{\begin{array}{ll}{\rm{\Gamma }}/{{\rm{\Gamma }}}_{\mathrm{ad}}, & \mathrm{if}\ {\rm{\Gamma }}/{{\rm{\Gamma }}}_{\mathrm{ad}}\gt {{\rm{\Lambda }}}_{0}\\ {{\rm{\Lambda }}}_{k+1}{e}^{-{z}_{k}/{z}_{k+1}}, & \mathrm{otherwise}\end{array}\right.,\end{eqnarray} \tag{ 12 }$$

where ${\rm{\Gamma }}=-\partial T/\partial z$. In the first term of Equation (11), the density dependence represents the correlation of eddy diffusion with breaking gravity waves in terrestrial atmospheres (Lindzen & Forbes 1983; Izakov 2001), and the exponential decay term allows the reduction in eddy diffusion as the gravity wave amplitude decreases to zero around pressure P₀. We do not calculate molecular diffusion in our climate model, but this is calculated in the photochemical model. The second term in Equation (11) is dependent on stability and includes additional eddy diffusion in neutral and unstable conditions, such as in Earth's troposphere. Equation (12) slowly reduces the contribution of stability in the transition between unstable and stable regions. The overall profile is then consistent with Earth examples in Massie & Hunten (1981) and Brasseur & Solomon (2006), and with the eddy diffusion profiles assumed in exoplanet photochemical-kinetics studies by other authors (e.g., Segura et al. 2007; Hu et al. 2012). The evolved environments we simulate have thick atmospheres and no oceans and therefore are more like Mars and Venus, whose eddy diffusivities are well modeled by breaking gravity waves alone (Izakov 2001), so in those cases Λ = 0. See Table 5 for nominal values.

Zoom In Zoom Out Reset image size

Table 5.  Adjustable Convection Model Parameters

Parameter	O2-dominated	Venus-like	Aqua planet	Description
${f}_{{\rm{z}}}$	0.01	0.01	0.01	Mixing length fraction
U [m s−1]	0.1	0.1	10	Surface wind speed
z0 [m]	0.005	0.005	0.0002	Surface roughness length
${\rho }_{0}$ [kg m−3]	1	1	1	Baseline density for eddy diffusion
P0 [Pa]	1	1	1	Pressure for breaking of gravity waves in eddy diffusion
${K}_{{\rm{m}},0}$ [m2 s−1]	0.5	1.0	0.5	Eddy scaling coefficient for mass, gravity waves
${K}_{{\rm{m}},1}$ [m2 s−1]	⋯	⋯	20	Eddy scaling coefficient for mass, stability
${K}_{{\rm{m}},1}$ [m2 s−1]	⋯	0.5	⋯	Constant minimal eddy diffusion
${{\rm{\Lambda }}}_{0}$	⋯	⋯	0.1	Minimum stability factor
${f}_{\min }$	10−4	10−4	10−4	Minimum stability scaling factor for heat
${f}_{\max }$	0.15	0.15	0.15	Maximum stability scaling factor for heat

Download table as:  ASCIITypeset image

In stable conditions, the eddy diffusion coefficient for heat is reduced compared to that for momentum (e.g., Kondo et al. 1978). This reduction occurs gradually as a function of stability; Kondo et al. (1978) provided a fit based on Richardson number (Ri). However, our mixing length code cannot compute Ri in the free atmosphere, because it does not include an explicit description of the vertical shear in horizontal wind, so we use Γ/Γ_ad as a representative stability parameter for this purpose and reduce the eddy diffusion coefficient for heat by

$$\begin{eqnarray}&&{K}_{{\rm{h}}}={f}_{{\rm{h}}}{K}_{{\rm{m}}}{e}^{-\alpha {\rm{\Gamma }}/{{\rm{\Gamma }}}_{\mathrm{ad}}},\end{eqnarray} \tag{ 13 }$$

where we set

$$\begin{eqnarray}&&\alpha \equiv \displaystyle \frac{1}{2}\mathrm{ln}\left(\displaystyle \frac{{f}_{\min }}{{f}_{\max }}\right),\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{f}_{{\rm{h}}}\equiv {f}_{\max }{e}^{\alpha }.\end{eqnarray} \tag{ 15 }$$

The minimum and maximum scalings from K_m to K_h are denoted by f_min and f_max, respectively, and given in Table 5.

Equations (5)–(11) are composed almost entirely of fundamental physical properties of the local atmospheric conditions: pressure, temperature, density, gravity, and specific heat capacity. Since the vertical transport of heat, mass, and momentum in terrestrial exoplanet atmospheres is generally unknown, we chose to use this physically motivated approach to convection, rather than convective adjustment to a prespecified lapse rate that can be different from the adiabatic lapse rate (e.g., Earth is ∼6.5 K km⁻¹, Mars ∼2.5 K km⁻¹, and Venus ∼8.0 K km⁻¹; Catling & Kasting 2017). The parameterizations we have chosen are able to reproduce eddy diffusion coefficients for mass consistent with measurements and kinetics studies of Earth and Venus (see Figure 21).

To simulate a potentially habitable, ocean-covered (aqua) planet, we implemented an option for latent heat exchange and kinematic mixing of a condensable gas. The layer-by-layer vertical mixing of a condensable species is calculated by solving the continuity-transport equation in 1D:

$$\begin{eqnarray}&&\displaystyle \frac{\partial r}{\partial t}=-\displaystyle \frac{1}{\rho }\displaystyle \frac{\partial {F}_{m}}{\partial z},\end{eqnarray} \tag{ 16 }$$

where the kinematic flux is given by

$$\begin{eqnarray}&&{F}_{{\rm{m}}}=-\rho {K}_{{\rm{m}}}\displaystyle \frac{\partial r}{\partial z}.\end{eqnarray} \tag{ 17 }$$

Here, K_m is the eddy diffusion coefficient for mass transport, and r is the mixing ratio of the condensable species. As a result of either radiative cooling or an upwelling condensable species, a gas may become supersaturated (i.e., when the partial pressure at a level exceeds the saturation vapor pressure). We assume sufficient cloud condensation nuclei such that any condensate above the saturation point condenses immediately (i.e., within a single time step). The saturated mass mixing ratio is determined by

$$\begin{eqnarray}&&{r}_{\mathrm{sat}}=\displaystyle \frac{\mu }{{\mu }_{\mathrm{atm}}}\displaystyle \frac{{P}_{\mathrm{sat}}}{P}.\end{eqnarray} \tag{ 18 }$$

At the globally averaged temperatures, it is unlikely that the globally averaged water mixing ratios will approach saturation, though there is nonetheless condensation, evaporation, and latent heat exchange. The simplest solution here is to use a different mixing ratio criterion for saturation. We compute condensation if relative humidity exceeds that maximum humidity using a modified form of Equation (2) from Manabe & Wetherald (1967) for Earth:

$$\begin{eqnarray}&&{rh}={{rh}}_{{\rm{s}}}\max \left(\displaystyle \frac{p/{p}_{{\rm{s}}}-0.02}{1-0.02},0.2\right),\end{eqnarray} \tag{ 19 }$$

where s refers to the surface value. The modification is our specification of a minimum value that depends on the surface humidity (Manabe & Wetherald 1967 made a similar adjustment to obtain Earth's stratospheric H₂O abundance). While Manabe & Wetherald (1967) use rh_s = 0.77, we find that rh_s = 0.70 matches our globally averaged Earth data best (see Appendix C).

To determine the heating rate due to condensation, we evaluate the conserved quantity, the specific moist enthalpy:

$$\begin{eqnarray}&&h=({c}_{{\rm{p}}}{r}_{{\rm{a}}}+{c}_{{\rm{p}},{\rm{v}}}{r}_{{\rm{v}}}+{c}_{{\rm{p}},{\rm{c}}}{r}_{{\rm{c}}})T+{{Lr}}_{{\rm{v}}},\end{eqnarray} \tag{ 20 }$$

where the subscript c refers to the condensed form, v refers to the vapor phase of the condensable, a refers to the vapor-phase bulk atmosphere, and L is the latent heat of evaporation. The expansion of the layer reduces the net latent heat, as a result of isobaric work (PdV), which per unit mass is equivalent to (R/μ) dT. Including this term with the change in specific moist enthalpy due to a change in condensate and the resultant temperature change, and then solving for the change in temperature of the air parcel, we find the net heating rate during time step Δt is

$$\begin{eqnarray}&&{q}_{\mathrm{lh}}=\displaystyle \frac{{\rm{\Delta }}T}{{\rm{\Delta }}t}=\displaystyle \frac{{\rm{\Delta }}{r}_{{\rm{c}}}}{{\rm{\Delta }}t}\displaystyle \frac{L-T({c}_{{\rm{p}},{\rm{c}}}-{c}_{{\rm{p}},{\rm{v}}})}{{c}_{{\rm{p}}}(1-{r}_{{\rm{v}}}-{r}_{{\rm{c}}})+{c}_{{\rm{p}},{\rm{v}}}{r}_{{\rm{v}}}+{c}_{{\rm{p}},{\rm{c}}}{r}_{{\rm{c}}}-R/\mu },\end{eqnarray} \tag{ 21 }$$

where the quantity being condensed is ${\rm{\Delta }}{r}_{{\rm{c}}}\equiv {r}_{{\rm{c}}}-{r}_{\mathrm{sat}}$, r is the mass mixing ratio, and c_p is the specific heat.

We compute the surface fluxes of sensible heat, latent heat, and condensable gas using similarity theory (e.g., Pierrehumbert 2011, pp. 395–409), given by

$$\begin{eqnarray}&&{F}_{{\rm{c}},\mathrm{sh}}=\rho {c}_{{\rm{p}}}{C}_{{\rm{D}}}U({T}_{{\rm{s}}}-{T}_{\mathrm{sbl}})\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}&&{F}_{{\rm{c}},\mathrm{lh}}=\rho {{LC}}_{{\rm{D}}}U({r}_{{\rm{s}}}-r)\end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray}&&{F}_{{\rm{m}}}=\rho {C}_{{\rm{D}}}U({r}_{{\rm{s}}}-r),\end{eqnarray} \tag{ 24 }$$

where U is the surface wind speed, T_s is the temperature of the surface, T_sbl is the temperature of the surface atmospheric boundary layer, r_s is the surface mixing ratio boundary condition (determined from a humidity assumption or fixed mixing ratio), and C_D is the coefficient of drag at the surface, which we calculate as in Pierrehumbert (2011):

$$\begin{eqnarray}&&{C}_{{\rm{D}}}={\left(\displaystyle \frac{k}{\mathrm{log}z/{z}_{0}}\right)}^{2},\end{eqnarray} \tag{ 25 }$$

where z₀ is the surface roughness height, which we specify as 0.005 m for the desiccated planets (corresponding to flatlands with no vegetation) and 0.0002 m for open ocean (Davenport et al. 2000; Jarraud 2008). As the atmospheres are much thicker than Earth, we assume surface wind speeds of 0.1 m s⁻¹, much less than on Earth but consistent with measurements of Venus's surface atmosphere (Ainsworth & Herman 1975; Counselman et al. 1979). For the aqua planet, we use 10 m s⁻¹, which is roughly consistent with Earth's global-average wind speed (7 m s⁻¹; Meissner et al. 2001).

We include the radiative effects of aerosols in cloudy cases for both the aqua planet and Venus-like atmospheres. Since this is a 1D climate model, we present end-member cases of global clear and cloudy cases, which can be used to assess the observational discriminants of these atmospheres. We have not yet implemented a cloud microphysics model in our climate model, so we specify altitude-dependent differential optical depths and aerosol optical properties from Earth or Venus, as appropriate. Aerosol parameters are described in Section 2.6.3.

We first approach convergence using four-stream radiative transfer for a single global-average-approximation solar zenith angle of ∼60°. The final convergence of VPL Climate is conducted using eight-stream radiative transfer with four-point Gaussian quadrature integration of the stellar radiance, which increases the accuracy of the convergence result (Cronin 2014; Hogan & Hirahara 2015). In future work, we will quantitatively assess the impact of solar radiance integration in our climate model. Kitzmann (2016) determined that increasing the number of radiative transfer streams from eight to 16 had no effect on the radiative fluxes.

In addition to the updates presented in the preceding appendix, we validated Earth using the boundary conditions given in Table 7. In particular, this includes fluxes (not fixed mixing ratios) of the important Earth constituents CO, CH₄, NO, H₂S, H₂SO₄, OCS, N₂O, CS₂, and dimethyl sulfide (C₂H₆S), along with vertically distributed fluxes of SO₂. Oxygen is fixed at the surface at 21% and nitrogen at 78%. Carbon dioxide is fixed at the surface to the current value of 400 ppm. The flux of NO is calibrated to best reproduce the effects of NO_x chemistry, as we do not include chlorine species or other pollutants not listed here. The water vapor profile is from Manabe & Wetherald (1967) using a fixed surface humidity of 70%, assuming the tropopause is at 14 km and that the minimum saturation value in the troposphere is 20% of the surface value. Vertical profiles of important species are shown in Figure 22 along with several sources for comparison.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Table 7.  Earth Validation Lower Boundary Conditions

Species	Boundary condition
O	${v}_{\mathrm{dep}}=1$
O2	r0 = 0.21
H2O	rh0 = 0.70
H	${v}_{\mathrm{dep}}=1$
OH	${v}_{\mathrm{dep}}=1$
HO2	${v}_{\mathrm{dep}}=1$
H2O2	${v}_{\mathrm{dep}}=0.2$
H2	${v}_{\mathrm{dep}}=2.4\times {10}^{-4}$
CO	${v}_{\mathrm{dep}}=0.03$, $F=3.7\times {10}^{11}$
HCO	${v}_{\mathrm{dep}}=1$
H2CO	${v}_{\mathrm{dep}}=0.2$
CH4	$F=1.6\times {10}^{11}$
NO	${v}_{\mathrm{dep}}=1.6\times {10}^{-2}$, $F=1.0\times {10}^{9}$
NO2	${v}_{\mathrm{dep}}=3.0\times {10}^{-3}$
HNO	${v}_{\mathrm{dep}}=1$
H2S	${v}_{\mathrm{dep}}=0.02$, $F=2.0\times {10}^{8}$
SO2	${v}_{\mathrm{dep}}=1.$, $F=9.0\times {10}^{9}$ dist. 14 km
H2SO4	${v}_{\mathrm{dep}}=1.$, $F=7.0\times {10}^{8}$
HSO	${v}_{\mathrm{dep}}=1$
OCS	${v}_{\mathrm{dep}}=0.01$, $F=1.5\times {10}^{7}$
HNO3	${v}_{\mathrm{dep}}=0.2$
N2O	$F=1.53\times {10}^{9}$
HO2NO2	${v}_{\mathrm{dep}}=0.2$
CO2	${r}_{0}=4.0\times {10}^{-4}$
CS2	$F=2.0\times {10}^{7}$
C2H6S(DMS)	$F=3.3\times {10}^{9}$
N2	r0 = 0.78

Notes. All unlisted species have a lower boundary condition ${v}_{\mathrm{dep}}=0$. Deposition velocities (${v}_{\mathrm{dep}}$) are given in cm s⁻¹, fluxes (F) are given in molecules cm⁻² s⁻¹, r₀ is a fixed surface mixing ratio, and rh₀ is the surface relative humidity.

Download table as:  ASCIITypeset image

This is the first time our photochemical model has been used for a Venus-like planet. This has required careful selection of boundary conditions, reactions, and the eddy diffusion profile. Specifically for Venus, few additional changes were required to the model itself. In addition to the model updates listed in Appendix B, we updated saturation vapor pressure data, added an additional lower boundary condition, compiled a new Venus input template, and conducted a validation against recent Venus literature.

We updated the saturation vapor pressure and temperature data for sulfuric acid solutions using finer intervals of 1% by fitting splines to published refraction and absorption data, which are available at sulfuric acid concentrations of 25, 38, 50, 75, 84.5, and 95.6% (Palmer & Williams 1975).

We drew heavily on recent Venus literature (primarily Krasnopolsky 2012; Zhang et al. 2012) to construct our Venus template, consisting of important species, their boundary conditions, and relevant reactions. The boundary conditions are listed in Table 8. As in previous work (e.g., Krasnopolsky 2012; Zhang et al. 2012), we do not extend our photochemical template for Venus to the surface. Instead, we set the lower boundary at 28 km (∼11.5 bar), as we do not include thermochemical reactions, which are typically treated as a separate regime from the photochemical region (e.g., Krasnopolsky 2013). Because the lower boundary is not the surface, the "deposition velocity" v depends on the scale height and eddy diffusion coefficient (Krasnopolsky 2012):

$$\begin{eqnarray}&&v=\displaystyle \frac{{K}_{m}}{2H},\end{eqnarray} \tag{ 26 }$$

so the model can calculate the associated downward flux. The reactions in our Venus template are primarily reproduced from those given in Krasnopolsky (2012).

Table 8.  Venus Validation Lower Boundary Conditions

Species	Boundary condition
O2	${v}_{\mathrm{dep}}=2v$
H2O	${r}_{0}=3.0\times {10}^{-5}$
H2	${r}_{0}=4.5\times {10}^{-9}$
NO	${r}_{0}=5.5\times {10}^{-9}$
SO2	${r}_{0}=2.8\times {10}^{-5}$
OCS	${r}_{0}=1.0\times {10}^{-6}$
HCl	${r}_{0}=4.0\times {10}^{-7}$
CO2	r0 = 0.965
N2	r0 = 0.035

Notes. All unlisted species have a lower boundary condition ${v}_{\mathrm{dep}}=v\,={K}_{m}/2H$. Boundary conditions are either deposition velocities (${v}_{\mathrm{dep}}$), given in cm s⁻¹, or a fixed surface mixing ratio r₀. The Venus lower boundary is at 28 km (∼11.5 bar).

Download table as:  ASCIITypeset image

The validation for Venus was primarily for major climate and observable species, particularly those involving Venus's sulfuric acid clouds (e.g., H₂O and SO₂; see Figure 23). We find good agreement with recent modeling by Krasnopolsky (2012), on which our reaction and boundary condition template is based. However, there are still disagreements between Venus photochemical modeling in the literature among the abundances retrieved from different observations, particularly observed oxygen abundances and in the features of the SO₂ profile. Future improvements to our Venus reactions can be made as those researchers who focus on Venus solve these outstanding issues.

Our photochemical model includes simplistic aerosol formation, including sulfuric acid aerosols, based on production/loss lifetimes (Pavlov et al. 2001). We found this model fails to produce the specific distribution of the multiple aerosol modes of Venus and is unable to generate aerosols as large as the mode 3 (3.85 μm) particles. Nonetheless, it does self-consistently generate aerosols at altitudes consistent with the cloud deck of Venus. We find that the transmission spectrum for our model Venus is roughly consistent with the cloud parameterization from Crisp (1986). Therefore, this provides an acceptable starting point for self-consistent climate–chemistry modeling of Venus-like exoplanets with aerosol formation.