NIR-Driven Moist Upper Atmospheres of Synchronously Rotating Temperate Terrestrial Exoplanets

TL;DR

This study uses 3D climate modeling to explore how water vapor distribution in the upper atmospheres of tidally locked Earth-sized exoplanets affects the detectability of water signatures in transmission spectra, influenced by stellar flux and spectral type.

Contribution

It provides new insights into the vertical water vapor transport driven by radiative heating and cloud effects in synchronously rotating exoplanets, enhancing understanding of H₂O detectability.

Findings

Water vapor in upper atmospheres increases with incident flux, but more gradually than 1D models suggest.

Upper atmospheric water vapor correlates with near-infrared incident flux, affecting spectral signatures.

Higher irradiation can amplify H₂O spectral features, easing observational detection.

Abstract

H$_2$O is a key molecule in characterizing atmospheres of temperate terrestrial planets, and observations of transmission spectra are expected to play a primary role in detecting its signatures in the near future. The detectability of H$_2$O absorption features in transmission spectra depends on the abundance of water vapor in the upper part of the atmosphere. We study the three-dimensional distribution of atmospheric H$_2$O for synchronously rotating Earth-sized aquaplanets using the general circulation model (GCM) ROCKE-3D, and examine the effects of total incident flux and stellar spectral type. We observe a more gentle increase of the water vapor mixing ratio in response to increased incident flux than one-dimensional models suggest, in qualitative agreement with the climate-stabilizing effect of clouds around the substellar point previously observed in GCMs applied to synchronously rotating planets. However, the water vapor mixing ratio in the upper atmosphere starts to increase while the surface temperature is still moderate. This is explained by the circulation in the upper atmosphere being driven by the radiative heating due to absorption by water vapor and cloud particles, causing efficient vertical transport of water vapor. Consistently, the water vapor mixing ratio is found to be well-correlated with the near-infrared portion of the incident flux. We also simulate transmission spectra based on the GCM outputs, and show that for the more highly irradiated planets, the H$_2$O signatures may be strengthened by a factor of a few, loosening the observational demands for a H$_2$O detection.

Tables

Table 1. Table 1 Stellar Properties

Star	Type	$T_{\mathrm{eff}}$	$L$ [$L_{\odot }$]	$M$ [$M_{\odot }$]	$R$ [$R_{\odot }$]	reference	source of spectrum
Sun	G2V	5772 K	1	1	1	-	Lean et al. (2005), SPARC/SOLARIS111http://solarisheppa.geomar.de/ccmi
HD 22049	K2V	5039 K	0.32	0.82	0.74	Baines & Armstrong (2012)	Segura et al. (2003)
Kepler-186	M1V	3755 K	0.055	0.544	0.523	Torres et al. (2015)	Allard et al. (2012)
GJ 876	M4V	3129 K	0.0122	0.37	0.3761	von Braun et al. (2014)	Domagal-Goldman et al. (2014)

Table 2. Table 2 Orbital period, to which the rotation period is equated, in the case of S X = 1 subscript 𝑆 𝑋 1 S_{X}=1 .

Star	$T_{\mathrm{eff}}$	Period [days] ($S_X=1$)
Sun	5772 K	365.3
HD 22049	5039 K	171.6
Kepler 186	3755 K	56.2
GJ 876	3129 K	22.0

Table 3. Table 3 The short-wave bands.

#	lower limit [µm]	upper limit [µm]
1	0.200	0.385
2	0.385	0.500
3	0.500	0.690
4	0.690	0.870
5	0.870	0.900
6	0.900	1.08
7	1.08	1.12
8	1.12	1.16
9	1.16	1.20
10	1.20	1.30
11	1.30	1.34
12	1.34	1.42
13	1.42	1.46
14	1.46	1.52
15	1.52	1.56
16	1.56	1.62
17	1.62	1.68
18	1.68	1.80
19	1.80	1.94
20	1.94	2.00
21	2.00	2.14
22	2.14	2.50
23	2.50	2.65
24	2.65	2.85
25	2.85	3.15
26	3.15	3.60
27	3.60	4.10
28	4.10	4.60
29	4.60	10.0

Equations

$$\Delta F(\lambda )=\frac{(R_{\mathrm{p}}+h_{\mathrm{eff}}(\lambda ){)}^2}{R_{\star }^2}.$$

$$k=2\pi r_{\mathrm{cld}}^2\cdot n=3\rho x_{\mathrm{cld}}/(2{\rho }_{H2O}{r}_{\mathrm{cld}}),$$

$$\frac{T}{\theta }\left(\frac{\partial \theta }{\partial t}+(\mathbf{v}\cdot \nabla )\theta \right)=\frac{Q_r+Q_l}{\rho c_p},$$

$$\frac{T}{\theta }\frac{\partial \theta }{\partial t}=-\frac{T}{\theta }\left(u\frac{\partial \theta }{\partial x}+v\frac{\partial \theta }{\partial y}\right)-\frac{T}{\theta }\left(\omega \frac{\partial \theta }{\partial p}\right)+\frac{Q_r}{\rho c_p}+\frac{Q_l}{\rho c_p},$$

Full text

1. Introduction

One of the primary interests in future observations of Earth-sized planets is the presence of H2O, both as a tracer of habitable conditions and as a clue to the origin and evolutional pathways of the planets. Its spectral signatures may be targeted in the near future through transmission spectroscopy using JWST or next-generation ground-based telescopes. Although H2O signatures have been detected in the atmospheres of hot Jupiters (e.g. Tinetti et al., 2007; Sing et al., 2016), detecting molecular signatures, including H2O, on temperate terrestrial planets is exceedingly challenging (Cowan et al., 2015) because of the small planetary radius and the small scale height (due to the lower temperature and presumably larger mean molecular weight). Additionally, while H2O is a dominant opacity source of Earth’s atmosphere, its signatures in the transmission spectra of temperate planets are not strong due to the “cold trap”, i.e., water vapor evaporated from the surface is transported upward, but condenses as air rises and cools, and most of it precipitates, leaving the stratosphere—where transmission spectroscopy typically probes—fairly dry. Several studies modeling transmission spectra of the Earth found modest signatures of water vapor (e.g. Ehrenreich et al., 2006; Kaltenegger & Traub, 2009; Bétrémieux & Kaltenegger, 2013; Misra et al., 2014). However, the efficiency of the cold trap depends on the atmospheric properties in general. Thus, modeling the water vapor mixing ratio profiles is an issue of interest in the context of evaluating the observability of atmospheric signatures of transiting terrestrial planets.

The mixing ratio of water vapor in the stratosphere is also closely related to planetary water loss and thus planetary habitability. Kasting et al. (1993) considered the maximum stratospheric water vapor mixing ratio that would allow for a planet to sustain an amount of water equivalent to the Earth’s oceans over the Earth’s age, $\sim 3\times 10^{-3}$, as a criterion for the inner edge of the habitable zone (the threshold for the moist-greenhouse stage). Motivated mainly by this “habitability” aspect, atmospheric models with different levels of complexity have examined stratospheric water vapor, or equivalently, the efficiency of the cold trap, with varying planetary properties.

Kasting et al. (1993) and Kopparapu et al. (2013) used a simple one-dimensional (1D) model with a saturated troposphere and an isothermal stratosphere to examine the response of the stratospheric water vapor mixing ratio to increased instellation (see also Kasting et al., 2015). They found a transition from an Earth-like dry upper atmosphere to a water-dominated atmosphere, driven by the increase of the surface temperature. The above-mentioned critical mixing ratio ($\sim 3\times 10^{-3}$) is reached when the surface temperature is $\sim$340 K. Kodama et al. (2015) also showed a similarly rapid transition to a water-dominated atmosphere as a function of incident flux using a 1D radiative-convective model.

Wordsworth & Pierrehumbert (2013, 2014) argued that the degree to which a condensible gas (e.g., water) is transported from the surface is represented by $\mathcal{M}=m_{v}p_{v}L/m_{n}p_{n}c_{p}T$ where $L$, $c_{p}$, and $T$ are the specific latent heat of the condensing gas, the specific heat capacity at constant pressure of the non-condensing gas, and surface temperature, respectively, while $m_{\{v,n\}}$ and $p_{\{v,n\}}$ are molar masses and the partial pressures of the condensing ($v$) or non-condensing ($n$) gases. A moist upper atmosphere occurs when $\mathcal{M}>1$. This indicates that water vapor profiles also depend on the background atmospheric properties, in addition to incident flux.

Atmospheric profiles of template terrestrial planets have also been studied using photochemical models combined with 1D radiative-convective models (Segura et al., 2003, 2005; Rauer et al., 2011; Hedelt et al., 2013; Rugheimer et al., 2013, 2015). Rauer et al. (2011) and Rugheimer et al. (2013, 2015) observed an enhanced water mixing ratio in the stratosphere for planets around late-type stars, which was attributed to the higher surface temperatures and photochemical production from CH4.

Inevitably, 1D models rely on several assumptions that cannot be treated within the model, including the global transport of heat and water vapor, as well as the effect of clouds. Various three-dimensional (3D) general circulation models (GCMs) revisited the atmospheric structure of terrestrial planets and have demonstrated the importance of accounting for the 3D circulation and climate heterogeneity (e.g. Ishiwatari et al., 2002; Abe et al., 2011; Leconte et al., 2013a, b; Wolf & Toon, 2014, 2015). Among others, Yang et al. (2013) pointed out that for synchronously rotating planets, the convection around the substellar point produces optically thick clouds that increase the planetary albedo and contribute to keep the surface temperature moderate for a much wider range of instellation (see also Yang et al., 2014; Way et al., 2015; Kopparapu et al., 2016). This situation is of particular relevance to the planets in the habitable zones of late-type stars, as they are likely to be synchronously rotating due to tidal locking (Dole, 1964; Kasting et al., 1993), although this need not always occur (Goldreich & Peale, 1966; Leconte et al., 2015). They primarily studied the highest instellation that allows for the model to converge as a proxy for the inner edge. While they also presented the stratospheric water vapor amounts, the discussion was limited.

Motivated by these works on the stabilized climate of synchronously rotating planets, we are specifically interested in the water vapor transport into the upper atmosphere, not only because it affects the water loss rate (and thus the conventional criterion for habitable conditions), but also because it limits our ability to detect H2O signatures of planets with surface liquid water, through transmission spectroscopy. Indeed, planets around later-type stars, which are more likely to be synchronously rotating, are best suited for future transit observations for several reasons, including the higher transit probability, larger signals in transmission spectra, and the abundance of late-type stars in the nearby universe.

In order to obtain better insights into the processes that control the water vapor mixing ratio in the upper atmosphere of synchronously rotating terrestrial planets with surface water, we conduct a series of experiments using the ROCKE-3D GCM, which calculates the atmospheric dynamic, thermodynamic, and radiative fields up to $0.1$ mbar ($\sim$65 km in the case of the Earth), sufficiently resolving the stratosphere. We specifically study the dependence of the water vapor mixing ratio on the incident flux (total incident flux and the spectral type), fixing planetary parameters (radius, gravity, atmospheric composition, and atmospheric pressure). We also present synthetic transmission spectra based on GCM outputs in order to evaluate the detectability of H2O signatures with future transmission spectroscopy.

The paper is organized as follows: Section 2 describes the setup of our GCM experiments and calculations of transmission spectra. Section 3 discusses our results of water vapor profiles as a function of the incident spectrum, and shows the modeled transmission spectra. The sensitivities of our results to the ocean model, the CO2 abundance, the spin period, and the water vapor continuum absorption are studied in Section 4. Finally, Section 5 summarizes our findings.

2. Model

2.1. GCM experiments

We use the ROCKE-3D GCM (Way et al., 2017) to obtain 3D atmospheric structures of terrestrial planets with surface water under varying irradiation, focusing on the H2O profile in the upper atmosphere. ROCKE-3D is a generalization of the ModelE2 GCM (Schmidt et al., 2014), which has been developed for the Earth at NASA’s Goddard Institute for Space Studies. The details of the ROCKE-3D GCM are presented in Way et al. (2017).

In the following, we first describe the specific configurations in our experiments to clarify our targets. Then, we also discuss in some detail the convection/cloud parameterizations and the radiation parameterization used in our GCM, as these are of particular importance to our results.

2.1.1 Our Configurations

In our experiments, the radius and the surface gravity of the planet are set to Earth’s value.

We only consider aquaplanets, i.e., planets whose surfaces are wholly covered with water. We use a thermodynamic ocean with 50 m depth for the fiducial runs for simplicity. However, in Section 4.1, we also present the results of runs with a fully dynamic ocean, to show that our conclusions are not very sensitive to ocean heat transport.

Our atmosphere has a mean surface pressure of 1 bar, and is composed mainly of N2, with only 1 ppm CO2, as in Kopparapu et al. (2016). In Section 4.2, we will briefly discuss the dependence of our results on the mixing ratio of CO2. We did not include the effects of aerosols and photochemistry.

The model was configured to have a horizontal grid with 4 degree resolution in latitude and 5 degree resolution in longitude, and 40 vertical layers covering up to 0.14 mbar (corresponding to $\sim 65$ km in the case of the Earth).

In our experiments, we mainly changed two parameters: the total incident flux and the spectral type of the star. In the following, we represent the incident flux normalized by the solar constant by $S_{X}$. We consider four types of stars: Sun (G2V), HD 22049 (K2V), Kepler-186 (M1V), and GJ 876 (M4V). The spectra of HD 22049 and GJ 876 are obtained from Segura et al. (2003) and Domagal-Goldman et al. (2014), respectively, while the spectrum of Kepler-186 is a modeled spectrum from BT-Settl (Allard et al., 2012) based on the stellar parameters reported in Table S1 of Quintana et al. (2014). The physical parameters of these stars, as well as the references, are summarized in Table 1, and the spectra are displayed in Figure 1.

Throughout the paper, we exclusively consider synchronously rotating planets due to their relevance to future transmission observations of close-in systems, as discussed in Section 1. The orbital period (identical to the spin period) is changed according to Kepler’s 3rd law, to be consistent with the total incident flux the planet receives and the stellar mass. Thus, the incident flux and the spin period are changed simultaneously for the fiducial set of runs. The orbital/spin periods corresponding to $S_{X}=1$ are shown in Table 2. We also isolate the effect of the spin period in Section 4.3.

Our fiducial models with a thermodynamic ocean typically equilibrate in $\sim$ 30 Earth years in model time. After the steady-state is reached, we average the 3D atmospheric profiles over an integer number of orbits that is approximately equal to 10 Earth years. All of the results presented in this paper are based on the averaged data created this way.

2.1.2 Cloud Scheme

In the ROCKE-3D GCM, Convection is triggered when air lifted from one model layer to the next higher layer becomes buoyant, according to a definition that includes both temperature and the molecular weight of H2O relative to dry air in assessing whether the lifted air is less dense than the environment. The cumulus parameterization uses a mass-flux closure determined by the mass required to establish neutral buoyancy at a cloud base over a specified adjustment time. The total mass flux is divided into two components (“plumes”) that entrain at different rates (Del Genio et al., 2007). The plume rises until the updraft speed decreases to zero. Entrainment is proportional to parcel buoyancy and inversely proportional to the square of updraft speed, following Gregory (2001). Mass is detrained at all levels above the level of neutral buoyancy. A downdraft forms at any level when an equal mixture of cloud and environment is negatively buoyant. Updraft convective condensate is partitioned into fractions that precipitate, detrain, and are carried upward by assuming a Marshall-Palmer particle size distribution and empirical size-fall speed relationships and calculating the fraction of the size distribution whose fall speed exceeds, is comparable to, or is less than the diagnosed updraft speed (Del Genio et al., 2005).

Detrained condensate is treated by the stratiform cloud microphysics (an updated version of Del Genio et al., 1996). Stratiform cloud fraction is diagnosed as a function of relative humidity and stability, while condensed water evolves using a prognostic equation that includes representation of all relevant microphysical source and sink processes. Partitioning between the liquid and ice phase is specified as a function of temperature in the mixed phase temperature range, and then modified by a parameterization of diffusional growth of ice when snow falls into a supercooled liquid water layer.

2.1.3 Radiation Scheme

For the radiation scheme, we adopt the Suite of Community RAdiative Transfer codes based on Edwards and Slingo111https://code.metoffice.gov.uk/trac/socrates (SOCRATES, Edwards, 1996; Edwards & Slingo, 1996) as described in Way et al. (2017). This radiation scheme uses the two-stream approximation of Zdunkowski & Korb (1985) with a diffusivity of $D=1.66$ for the long-wave (thermal) component, and of Zdunkowski et al. (1980) with $D=2$ for the short-wave (stellar) component. Opacities are treated using the correlated-$k$ method (Goody et al., 1989; Lacis & Oinas, 1991), with $k$-terms derived using exponential sum fitting of transmissions (Wiscombe & Evans, 1977) based on HITRAN2012 (Rothman et al., 2013). H2O continuum absorption is treated using the CAVIAR continuum (Ptashnik et al., 2011). Both H2O and CO2 lines are cut off at 25 cm*-1*, and the pedestal is removed from the H2O lines as required by the water vapor continuum parameterization. Overlapping gaseous absorption is treated using equivalent extinction (Edwards, 1996; Amundsen et al., 2017).

The optical properties of liquid clouds were derived using Mie theory, with the effective radius diagnosed from the cloud water content assuming a cloud droplet number concentration of 60 cm*-3* for liquid clouds. The optical properties of ice crystals are described in Edwards et al. (2007), and are based on the representation of ice aggregates introduced by Baran et al. (2001), with the effective dimension calculated using a parameterization based on temperature. Convective and stratiform clouds are treated separately, and if both are present in the same grid box they are treated as randomly overlapping to calculate the total cloud fraction within the grid box. Vertical cloud overlap is treated using the maximum-random overlap assumption, where clouds in adjacent layers overlap maximally, while clouds separated by one or more layers overlap randomly.

The configuration of the long-wave radiation scheme adopted here is the same as that used by the UK Met Office for global atmosphere configuration 7.0 (GA7.0; D. Walters et al. 2017, in preparation). For short-wave radiation, we use 29 bands, as indicated by the gray vertical lines in Figure  1 and listed in Table 3, in order to ensure accuracy of the short-wave absorption. The wavelengths are weighted internally in each band by each stellar spectrum when deriving $k$-coefficients and cloud optical properties to further improve accuracy. We check the accuracy of the short-wave radiation scheme by comparing the flux and heating rate obtained with $29$ bands to higher-resolution calculations, and find a good match; see the Appendix A for more details.

2.2. Simulation of transmission spectra

We also simulate transmission spectra based on the outputs of our GCM experiments, setting the modeled planet in a transiting geometry, in order to quantify the absorption features of water vapor. We follow the procedure described in Appendix B of Way et al. (2017), which is briefly summarized below.

First, we obtain the vertical profiles of temperature, water vapor, and cloud particles along the terminator (i.e., dayside-nightside boundary) by linearly interpolating the profiles of the GCM gridpoints adjacent to the terminator. Then, we trace the optical paths of the transmitted ray according to the profile of the refraction index, from the observer toward the stellar disk, in the same way as described in van der Werf (2008), following Misra et al. (2014). After that, the spectral opacity at different points along the path is calculated considering Rayleigh scattering, the absorption by atmospheric molecules, and absorption and scattering by cloud particles, excluding the effects of multiple scattering. The spectral signatures of the planetary atmosphere are discussed in terms of the transit depth $\Delta F(\lambda)$ where $\lambda$ is the wavelength, as well as the effective height, $h_{\rm eff}$, which is given as a function of planet radius, $R_{\rm p}$, and stellar radius, $R_{\star}$, by

[TABLE]

To compute the absorption coefficients of atmospheric molecules, we use line-by-line opacity tables based on the HITRAN2012 database (Rothman et al., 2013), prepared separately from the opacities used in the GCM.

Including the optical effects of clouds is tricky. Strictly speaking, one needs to know the size distribution and the spatial distribution of the cloud particles at the time scale of the planetary transit. In this paper, when calculating transmission spectra, we consider the effect of clouds only in an approximate manner with four assumptions. First, we assume that the radius of the cloud particles $r_{\rm cld}$ is $r_{\rm cld}=30~{}{\rm\mu m}$ (regardless of their phase), similar to the typical values of ice cloud particles on the Earth as well as in our GCM simulations. Second, we assume that the scattering cross section is given by that of the geometric optics limit, i.e. $2\pi r_{\rm cld}^{2}$. Third, we assume that cloud particles are uniformly distributed within each grid cell. Fourthly, we compute the transmitted light based on the optical properties averaged over 10 Earth years, rather than compute it based on the instantaneous optical properties before taking the temporal average. The third and fourth points above imply that the effects of clouds we find are the largest among the possible realizations given the particle size distribution, because of the nonlinear dependence of transmittance on the optical thickness. Under these assumptions, the absorption coefficient $k$ of the clouds is

[TABLE]

where $n$ is the number density of cloud particles, $\rho$ is the density of the ambient atmosphere, $x_{\rm cld}$ is the cloud water content provided by the GCM experiments, and $\rho_{\rm H2O}$ is the density of liquid/solid water, which is assumed to be $1{\rm g/cm}^{3}$. The uncertainties in the particle radius, scattering cross section, and the cloud water content could change the absorption coefficient by a factor of a few.

3. Results

3.1. General Trend in 3D Atmospheric Structures

Before examining the dependence of the water vapor mixing ratio in the upper atmosphere, we discuss the 3D structures of temperature and water vapor, as we need to know at which altitude the mixing ratio of water vapor becomes relatively uniform and hence can be regarded as representative of the humidity in the upper atmosphere.

Figure 2 shows the structures of temperature and water vapor mixing ratio along the equator for the planets around GJ 876 (M4V) and around the Sun (G2V), respectively, for varying incident flux normalized by the solar constant ($S_{X}$). In all panels the substellar point corresponds to $0^{\circ}$ longitude. In general, the temperature and the water vapor mixing ratio near the surface strongly depend on the distance from the substellar point, but are fairly uniform in the upper part of the atmosphere. Above about 1 mbar, which is located between 35 and 50 km altitude depending on the atmospheric temperature profiles, the water vapor mixing ratio becomes horizontally uniform within a few tens of percent. Therefore, in the following, we denote the globally averaged water vapor mixing ratio at 1.39 mbar by x_{\text{H{}_{2}O}}, and use it as representative of the water vapor abundance in the upper atmosphere.

We also observe optically thick clouds around the substellar point and a resulting increase of the planetary albedo as a function of total incident flux (not shown in figures), consistent with Yang et al. (2013, 2014), Kopparapu et al. (2016), and Way et al. (2016). As shown in Figure 2, the maximum surface temperature does not exceed 320 K even in the highly irradiated cases in our experiments where the water vapor mixing ratio near the model top is as high as $10^{-3}$; we will discuss this in the subsequent sections.

3.2. Response of Water Vapor Mixing Ratio

to Increased Incident Flux

In this section, we examine in more detail how the water vapor mixing ratio in the upper atmosphere depends on the incident flux (both the total flux and the spectral distribution).

The left panel of Figure 3 presents x_{\text{H{}{2}O}} versus $S_{X}$. Overall, we find an increase of x_{\text{H{}{2}O}} as a function of $S_{X}$. However, the slope is much more gradual than suggested by 1D models, in particular those of Kasting et al. (1993) and Wordsworth & Pierrehumbert (2013) with a small CO2 mixing ratio, which typically predict a rapid transition from an Earth-like dry stratosphere to an H2O-dominated atmosphere. Our result that the water vapor mixing ratio responds more slowly to the incident flux is qualitatively consistent with Yang et al. (2013), who studied synchronously rotating planets using CAM3 with an emphasis on the stability of the climate. Note that we cannot directly compare to their results, as our assumptions do not match perfectly; specifically, Yang et al. (2013) assumed a planet with a larger radius ($R=2R_{\earth}$) and higher gravity ($g=1.4g_{\earth}$), and a 60-day orbital period. The dependence on these other planetary parameters will be investigated in a future paper.

Because these are synchronously rotating planets where the surface temperature is stabilized by the optically thick clouds around the substellar point, one might think that this stabilized surface temperature would explain the slow response of water vapor to the increased total incident flux. To see the effect of the surface temperature, we have plotted x_{\text{H{}_{2}O}} as a function of surface temperature around the substellar point, where convection is strongest (Figure 3, middle panel). It is then clearly seen that the water vapor mixing ratio starts to rise at a fairly modest surface temperature and increases even more dramatically as a function of surface temperature than 1D models suggest. This indicates that the increasing surface temperature cannot by itself explain the increased stratospheric water vapor mixing ratios.

Instead, we find that the near-infrared (NIR) component of the incident flux appears to be a good predictor of stratospheric humidity across the full range of stellar spectral types (Figure 3, the right panel). Here, we integrated the star spectra between 0.9 and 3 µm (the interval is somewhat arbitrary) to obtain the NIR component of the incident flux. We will examine the reason for this in the next subsection.

3.3. Development of Large-scale Circulation

in the Upper Atmosphere

In Figure 4, we show the vertical profiles of temperature, water vapor mixing ratio, short-wave heating rate, and upward velocity, averaged over the region within a distance $0.2R_{\earth}$ from the substellar point. The variability of the profiles in the individual gridboxes that contribute to the average is negligible compared to the changes due to varying $S_{X}$. Rising motion exists at most altitudes in the substellar region. In the troposphere, where moist convection occurs, large-scale vertical motion is relatively insensitive to $S_{X}$. Above the convective layers, though, large-scale upward motion increases strongly with incident flux, concurrently with the increase in water vapor mixing ratio and heating rate in the upper atmosphere.

Here in the upper atmosphere we see the interplay among the radiative heating, water abundance, and the vertical air transport: when the incident flux increases, the radiative heating due to the absorption of short-wave radiation by water vapor and potential cloud particles, mainly in the NIR regime, also increases. The air then reacts by rising and adiabatically cooling. The induced upward motion of the atmosphere transports water vapor from the moist lower atmosphere to the upper atmosphere, lowering the vertical gradient of water vapor and increasing the moisture at high altitude, further increasing the short-wave heating.

To be more specific, the evolution of the profile satisfies the conservation equation for the potential temperature:

[TABLE]

or

[TABLE]

where $p$ is the pressure, $Q_{r}$ and $Q_{l}$ are the radiative and latent heating rates, $c_{p}$ is the specific heat of the atmosphere at constant pressure, $\rho$ is atmospheric density, $\{u,v,\omega\}$ are the conventional velocity components toward the east ($x-$axis), north ($y-$axis), and upward ($z-$axis; expressed here in pressure coordinates as $\omega\equiv dp/dt$), respectively. Given that the horizontal advection (the first term on the right side) has a negligible contribution since the horizontal gradients are small (Figure 2), and that the latent heating (the last term on the right side) is also negligible above the altitude to which convection penetrates, the increase of the radiative heating in the upper atmosphere must be balanced by the adiabatic cooling by vertical advection. The contribution of each term of the right side of equation (4) in the steady state is shown in Figure 5, for the substellar regions of low and high-irradiation cases with GJ 876 and the Sun. It is clearly seen that the net radiative heating in high-irradiation cases is compensated by the cooling by the vertical advection.

Figure 6 puts this upward motion in the substellar regions into a global context by presenting the velocity fields in the altitude-longitude plane along the equator for the low and high-irradiation cases for GJ 876 and the Sun. Here, the vertical components of the vectors are multiplied by 1000 to make them more visible, because large-scale vertical motions are about three orders of magnitude weaker than horizontal motions. The development of the large-scale circulation above the troposphere is clearly seen.

Altogether, Figures 4, 5 and 6 suggest that the stratospheric dynamical response to increased stellar radiation is similar for stars of different spectral type, but for a warmer star the onset of strong upward motion occurs at a much larger incident flux because less of the instellation is in the NIR where the radiation is more efficiently absorbed by water vapor and clouds. This is the reason for the behavior seen in the right panel of Figure 3.

We also note that the contribution of radiative heating by cloud particles above convective layers is not negligible. In order to see the radiative effect of clouds above convective layers, we created models that are comparable to the fiducial models, but where the influence of cloud particles above the 237 mbar level ($\sim 11$ km for the runs shown below) has been turned off in the radiation calculations (i.e., cloud particles in the upper atmosphere are transparent to radiation). Figure 7 displays the comparison between such “transparent high clouds” models (dashed black, dashed gray) and fiducial models (red, pink), for $S_{X}=1.0$ and $1.2$ with the incident spectrum of GJ 876. While these high clouds have a minimal effect on the surface temperature and atmospheric profiles in the convection layers, it is seen that radiative heating above the convection layer is substantially suppressed without the radiative effect of clouds, leading to a colder, and drier upper atmosphere. Note that we still see heating attributed to water vapor in the “transparent high clouds” models; for example, the water vapor mixing ratio and short-wave heating at 20-30 km altitude increase considerably from $S_{X}$ = 1.0 to $S_{X}$ = 1.2, even in the absence of cloud radiative effects, but the sensitivity to $S_{X}$ is weaker at higher altitudes. We also ran equivalent simulations for other types of stars and found that the overall effects of the “transparent high clouds” assumption are very similar.

In summary, atmospheric vertical motion and the radiative heating by water vapor and cloud particles interact with each other through positive feedbacks, and thus a modest increase of the incident radiative flux, $S_{X}$, can result in an order-of-magnitude or more increase of the water vapor mixing ratio in the upper atmosphere, while the surface temperature increases only modestly due to the stabilization by clouds. This picture is verified by a good correlation between the water vapor mixing ratio at 1.39 mbar and the NIR portion of the incident flux as shown in the right panel of Figure 3, insensitive to the stellar spectral types. Because the NIR flux is an external parameter that can be found from direct observations of the host star, it would be a useful guide to predict the high-altitude water vapor mixing ratio for a planet with surface liquid water.

3.4. Transmission Spectra

In this subsection, we show the dependence of the transmission spectra on the incident flux, taking the modeled planets around one of our M-type stars, GJ 876 as examples. Figure 8 shows the transmission spectra for the $S_{X}=0.6$, $1.0$ and $1.2$ cases, where the water vapor mixing ratio at 1.39 mbar is roughly $10^{-6}$, $10^{-4}$ and $3\cdot 10^{-3}$, respectively. The corresponding vertical profiles in the substellar region are shown in Figure 4, and note that the water vapor mixing ratio is even higher at medium altitudes (10-30 km). Although refraction due to the planetary atmosphere is included, its effect on the transmission spectra of planets around late-type stars like GJ 876 is minor (Bétrémieux & Kaltenegger, 2014; Misra et al., 2014).

Without the opacity due to clouds (colored dashed lines in Figure 8), the effective altitude of the water vapor absorption features increases from $\sim$15 km to $\sim$40 km when $S_{X}$ increases from $0.6$ to $1.2$. This level of effective altitude of water signatures is comparable to the strongest features (CO2) in the modeled visible-to-NIR transmission spectrum of the Earth (solid black line) , and substantially larger than the H2O signatures of the Earth. Note that the feature around 2.7µm is the overlap between H2O and CO2, and in our simulation cases with 1 ppm of CO2 the feature predominantly comes from H2O, while in the case of the Earth it mostly comes from CO2. Thus, for the highly irradiated cases shown here, H2O signatures will be a natural target when trying to observe the strongest atmospheric signatures of planets similar to the Earth.

However, we also find clouds in the upper atmosphere around the terminator, and the inclusion of clouds substantially raises the baseline of the transmission spectra. (Note that the high clouds discussed in the previous sections are those around the substellar point, while what matter here for the transmission spectra are those at the terminator.) These clouds at the terminator presumably form from water vapor that is transported from the substellar point and cools radiatively or adiabatically. Although these tenuous clouds at high altitude would be optically thin vertically, the grazing geometry of transmission observations cause them to significantly affect the transmission spectra due to the long optical path length. After including the effect of clouds, however, we still find the overall increase of water vapor signatures as the incident flux increases. Compared to the standard Earth model, or the low-irradiation case (water vapor mixing ratio $\sim 10^{-6}$) with the same assumptions for the cloud optical properties, the signal of a high-irradiation case could be larger by a factor of a few.

We can evaluate the detectability of these features for idealized cases where the only source of observational uncertainty is photon noise. For an Earth-sized planet 10 pc away, assuming a 6.5 m telescope (having the James Webb Space Telescope in mind) with 40 % throughput and a wavelength resolution $\mathcal{R}=$20, absorption features with 20 km differential effective altitude in this wavelength range (e.g., the feature around 2.7 µm for the cases with clouds) would become detectable (5$\sigma$) after $\sim 100$ hr of integration with an M4 star (assuming blackbody radiation with $R_{\star}=0.4R_{\odot}$, T=3500 K, similar to GJ 876), or after $\sim 15$ hr with an M8 star (with $R_{\star}=0.1R_{\odot}$, T=2500 K, similar to TRAPPIST-1), respectively.

4. Sensitivity Experiments

In this section we examine the sensitivity of our results to some of our assumptions. In particular, we discuss the effect of the assumption of a simple thermodynamic ocean (Section 4.1), the effect of the background CO2 mixing ratio (Section 4.2), the effect of the spin period (Section 4.3), and the effect of the uncertainty in the water vapor continuum absorption (Section 4.4).

4.1. Effect of Ocean Model

So far, we have used a thermodynamic ocean model with zero ocean heat transport and 50 m depth. Indeed, this simplified treatment has been regularly used in previous 3D atmospheric modeling studies of exoplanets. A more realistic approach is to include the dynamics of the ocean. Therefore, we check whether our results are sensitive to the prescription of the ocean by comparing our fiducial results and those obtained using a dynamic ocean module.

For this purpose, we ran the models for GJ 876 assuming a 900 m deep dynamic ocean. The models with the dynamic ocean typically equilibrated after a hundred to several hundred Earth years. In the same way as the fiducial models, we averaged the atmospheric profiles over approximately 10 Earth years after the model reached a steady state.

Figure 9 shows x_{\text{H{}{2}O}} for GJ 876 with varying $S_{X}$ both for our fiducial runs and those using the dynamic ocean module. There are minor differences between the two, with the maximum difference being less than an order of magnitude. The models with a dynamic ocean tend to have somewhat less x_{\text{H{}{2}O}} than the models with zero ocean heat transport. This appears to be related to the fact that ocean heat transport from the dayside to the nightside lowers the substellar surface temperature and weakens convection there, and thus there is less transport of water vapor to the upper atmosphere.

Overall, however, our fiducial models with a thermodynamic ocean give a reasonable order-of-magnitude estimate of water vapor mixing ratio in the upper atmosphere.

4.2. Effect of CO2

Our fiducial models assumed 1 ppm CO2, which was an arbitrary choice. Although a full investigation of the dependence on the CO2 mixing ratio is beyond the scope of this paper, here we briefly mention how much the results are changed if we assume an Earth-like CO2 mixing ratio. The results of the models with 300 ppm of CO2 are plotted in Figure 9. As shown, this level of variation of CO2 has a very small effect on the water vapor profile. Larger amounts of CO2 up to that of a CO2-dominated atmosphere could have a more substantial effect on the 3D profiles (Wordsworth & Pierrehumbert, 2013).

4.3. Effect of Spin Period

In our fiducial set of experiments, we changed the planetary spin period together with the stellar flux for a given star, because we consider synchronously rotating planets and the orbital period (which is equal to the spin period) must change to be consistent with the total incident flux and the host star mass according to Kepler’s 3rd law. Thus, strictly speaking, the change in the water vapor mixing ratio presented in Figure 3 includes the effect of both the varying incident flux and the varying spin period. Meanwhile, it is well-known that the spin period has first order effects on the atmospheric circulation and hence can potentially affect x_{\text{H{}_{2}O}}.

In order to isolate the effects of rotation, we performed experiments where we varied the spin period while keeping the incident flux fixed at $S_{X}=1$ with the GJ 876 spectrum, but allowing the planet to remain in synchronous rotation. The effect of varying spin period on x_{\text{H{}_{2}O}} is presented in Figure 10. It is considerably smaller than the dependence on $S_{X}$ seen in Figure 3, i.e., the increase in stratospheric water vapor as the planet is moved closer to its star is mostly an intrinsic sensitivity to the stellar flux rather than an effect of the changing planet rotation as the orbital period changes.

The spin effect is non-negligible, however; we find that x_{\text{H{}_{2}O}} becomes smaller as the spin period becomes very short, i.e., increasing planet spin partly offsets the effect of increasing stellar flux. Our interpretation of this tendency is as follows. At slow rotation (long spin period) the day-night circulation dominates the atmosphere. As the spin period becomes shorter, however, the atmosphere develops high-altitude zonal winds, as demonstrated in previous papers (Merlis & Schneider, 2010; Edson et al., 2011; Kopparapu et al., 2016). The super-rotating zonal wind seen in the equatorial region in the shorter period runs (not shown) efficiently transports water vapor from the dayside to the nightside, limiting the cloudiness in the substellar region, and resulting in a higher surface temperature (Kopparapu et al., 2016). Consistently, in the runs with a shorter spin period, the convection at the substellar point is deeper, and the temperature at the top of the convection where it moistens the atmosphere is colder, and thus the convection is less able to moisten the atmosphere there. Moreover, the water vapor transported upward by convection is more efficiently transported to the nightside by the zonal wind, as noted above; the water vapor mixing ratio on the nightside increases by about an order of magnitude when the rotation period decreases from 365 Earth days to 4 Earth days. These processes are likely to lead to the smaller upper atmosphere water vapor mixing ratio observed in our shorter shorter period runs.

4.4. Water Vapor Continuum

We have shown that the radiative heating caused by absorption of the incident NIR radiation by water vapor is one of the key processes governing the upper atmosphere humidity, and the water vapor continuum absorption is partly responsible for this. There is significant uncertainty in the strength of both the self- and foreign-broadened water vapor continua at NIR wavelengths (Campargue et al., 2016; Shine et al., 2016). Our fiducial models used the CAVIAR water vapor continuum (Ptashnik et al., 2011), which from recent laboratory measurements appears to overestimate particularly the self-broadened water vapor continuum (Shine et al., 2016). Another popular continuum model, MT_CKD (Mlawer et al., 2012), or its predecessor CKD (Clough et al., 1989), appears to be more in line with, or a bit below, the recent measurements of the self-broadened continuum. Measurements of the foreign-broadened continuum are sparse, but more consistent with CAVIAR as MT_CKD (and CKD) appears to underestimate the strength of this continuum (Shine et al., 2016).

In order to evaluate the impact of the water vapor continuum model on our results, we have also run simulations with the CKD 2.4 water vapor continuum. As CAVIAR and CKD 2.4 approximately bracket the available measurements, this should give a good indication of the sensitivity of our results to the adopted continuum model.

Figure 11 compares x_{\text{H{}_{2}O}} with the CAVIAR model (our fiducial model) and those with the CKD 2.4 model, indicating that the effect of the water vapor continuum is very small. We also compare in Figure 12 the vertical profiles of temperature and water vapor in the substellar region for the most highly irradiated (i.e., most moist) runs with the GJ 876 and Solar spectra. The resulting profiles overall closely match.

When the atmospheric structure is fixed, the difference in short-wave fluxes due to the different water vapor continua is observed predominantly in the lower atmosphere ($<$ 10 km), as expected since the strength of continuum absorption is $\propto P^{2}$, where $P$ is the pressure. Using the $P$-$T$ and specific humidity profiles at the substellar point from our fiducial run with the spectrum of GJ 876 and $S_{X}=1.2$, we find that the net short-wave clear-sky flux at the surface is larger with the CKD continuum compared to the CAVIAR continuum by about 7%, with the net surface clear-sky short-wave flux with the CAVIAR model being 838 W/m2. However, letting the system evolve, we find that such a difference is adjusted by small changes in temperature and humidity profiles and cloud patterns, leaving the basic atmospheric structure similar.

Based on these results, we are confident that the correlation between the humidity in the upper atmosphere and the amount of NIR incident radiation shown in Figure 3 is not sensitive to the choice of water vapor continuum model.

5. Summary

We have investigated the response of atmospheric water vapor profiles of synchronously rotating aquaplanets to changes in the incident stellar flux. We find a quasi-exponential increase of water vapor mixing ratios in the upper atmosphere as a function of total incident flux, resulting in a variation of several orders of magnitude as the total incident flux increases by a factor of two to three, but this is much more gradual than what is predicted by 1D models.

The increase of water vapor mixing ratio in response to larger incident flux is not primarily caused by the increased surface temperature as previous 1D models suggest, due to the relatively stable surface temperature at the substellar point. Instead, it is induced by the effect of the increased NIR portion of the incident flux on the high-altitude circulation. The enhanced radiative heating due to the absorption by water vapor and cloud particles drives a large-scale circulation of the air at altitudes above those to which convection penetrates, taking over the role of transporting water vapor from the moist lower atmosphere to the upper atmosphere, and further increasing the heating rate in the upper atmosphere. As a result, the moisture in the upper atmosphere increases dramatically, while the surface below it remains temperate. This picture is supported by the good correlation between the water vapor mixing ratio at 1.39 mbar and the NIR portion of the incident radiation, regardless of the spectral type of the host star.

Simulating transmission spectra based on the atmospheric profiles found in GCM experiments, we demonstrate that the increase of the water vapor mixing ratio at 1.39 mbar from $10^{-6}$ to $3\times 10^{-3}$ leads to an increase of the H2O absorption depth by a factor of a few, making the effective altitude of the H2O absorption features as high as 40 km. We also find the presence of high clouds at the terminator, which potentially increases the baseline effective altitude and thus weakens the absorption features. Our calculation may overestimate the effects of clouds on the transmission spectra due to the spatially and temporally averaged cloud properties of individual GCM grid cells. Additional uncertainties come from other assumptions we made. A more realistic inclusion of clouds in our simulated transmission spectra will be the scope of future work.

While our fiducial model assumed a thermodynamic ocean with zero ocean heat transport, comparable experiments with a dynamic ocean are also performed to demonstrate that this approximation of an ocean has only a modest effect on the water vapor mixing ratio, the size of which is within about an order of magnitude or less. Varying the CO2 mixing ratio from 1 ppm to 300 ppm has a negligible effect. We also confirmed that changing the planet spin period only has a secondary effect on the water vapor mixing ratio, resulting in an order-of-magnitude variation at most. Our results are not sensitive to the choice of water vapor continuum model within the known uncertainty.

In this paper, we limit our discussion on the water vapor profiles to the effects of the total incident flux and spectral energy distribution (i.e., the spectral type of the host star). These parameters are among the limited number of properties that can be constrained relatively easily by observations. In principle, however, water vapor distributions also depend on geophysical parameters including the surface gravity, atmospheric composition, atmospheric pressure (e.g., Turbet et al., 2016), and the planetary water abundance, as well as photochemistry, which depends on the detailed atmospheric profiles. The dependence of the water vapor mixing ratio on these planetary properties will be the focus of future study.

We thank Michael J. Way and Nancy Y. Kiang for providing the boundary conditions for ROCKE-3D and for the helpful discussions. The support from Maxwell Kelley and Igor Aleinov was critical in resolving the issues identified in the process of generalizing ROCKE-3D for the exotic configurations considered in this paper. We thank James Manners for providing access to SOCRATES and advice on coupling it to ROCKE-3D. We also thank Philip Von Paris for the Kepler-186 spectrum. Y. F. is deeply grateful to David S. Spiegel for the insightful discussions on the simulation suite for transmission spectra. Y. F.’s research was supported by an appointment to the NASA Postdoctoral Program at the NASA Goddard Institute for Space Studies, administered by Universities Space Research Association under contract with NASA, and the Grant-in-Aid No.15K17605 by the Japan Society for the Promotion of Science. We also acknowledge support from the NASA Astrobiology Program through the Nexus for Exoplanet System Science.

Appendix A Validation of radiation scheme

used in ROCKE-3D

In this paper we applied the ROCKE-3D GCM to planets around stars with different spectral types. In doing so, special caution is needed regarding the accuracy of the short-wave radiative transfer calculation, because the baseline model of ROCKE-3D was developed for the Earth, and the radiative transfer calculation was optimized for the solar spectrum. We therefore modified the bands used by the short-wave radiation calculation in SOCRATES in order to improve accuracy for low-mass stars while at the same time maintaining the accuracy for Sun-like stars. Our adopted short-wave bands in this study are given in Table 3.

In order to check the accuracy of the short-wave radiative transfer calculation used in the GCM, we compare it to a more accurate higher-resolution (280 bands) calculation using the same atmospheric profile. We use the atmospheric profile at the substellar point of a planet around GJ 876 with $S_{X}=1.2$ (see Figure 4), but without clouds. This profile was chosen because of its extremely high water vapor mixing ratio, which is very different from that of the Earth. Assuming this profile, we calculated short-wave radiative fluxes with 29 bands (used in the GCM) and with 280 bands, with both the Sun and GJ 876 as the incident spectra.

We show in Figure 13 the short-wave (stellar) downward fluxes, upward fluxes and heating rates we obtained. Differences between the 29- and 280-band configurations are small, $<$ 5 W/m2 for fluxes and $<$ 0.5 K/day for heating rates except at pressures $\ll$ 1 mbar, throughout the atmosphere. We therefore conclude that the radiation scheme adopted here is sufficiently accurate.