Evolution and Spectral Response of a Steam Atmosphere for Early Earth with a Coupled Climate–Interior Model

The early atmospheric evolution of Earth and other terrestrial planets is likely to be characterized by the presence of one or several magma oceans (MOs) with a hot and largely molten mantle lying above the core (Elkins-Tanton 2012). An MO consists of molten silicates undergoing turbulent convection and is a natural outcome of core formation and giant energetic impacts during planetary accretion (Bonati et al. 2019) or via comets and meteorites (Abe & Matsui 1985; Canup 2004).

During the accretion of planetesimals and following giant impacts, a transient atmosphere due to impact devolatilization was present above the MO, but the planet did not acquire its initial water content through ingassing from such an atmosphere. Recent models of ingassing of nebular hydrogen into the MO show that only 1% of water entered the planet in this way (Wu et al. 2018). The current general consensus is that Earth acquired its water inventory largely via accretion of chondritic materials (Russell et al. 2017). Also, experiments of metal-silicate partitioning require the presence of a hydrous MO to account for the characteristics of certain trace elements of Earth (e.g., Righter & Drake 1999). At any rate, when Earth's MO started to solidify, it already had a significant amount of water to be outgassed, no matter how it was acquired.

The chemical composition of the atmosphere built as a result of outgassing from the interior depends on the compositions of different meteoritic classes of materials, e.g., carbonaceous, ordinary, or enstatite chondrites as discussed by Schaefer & Fegley (2007, 2010). Furthermore, it has been argued that a mantle rich in hydrogen would give rise to a reduced atmosphere mainly consisting of CO, CH₄, NH₃, and H₂, whereas a mantle more abundant in oxygen would give rise to an oxidized atmosphere mainly consisting of CO₂ and H₂O (Zahnle et al. 2010; Gaillard & Scaillet 2014). Recently, the oxygen fugacities for the different meteoritic materials (both reduced and oxidized) have been explored in greater detail by Schaefer & Fegley (2017). In the present paper, we investigate the atmospheric evolution by assuming an oxidized mantle whereby the important factor remains the planetary outgassing.

Atmospheric 1D column models focusing on the "greenhouse effect" with surface temperatures exceeding the critical temperature of water (647 K) have been employed by Kasting (1988) and Nakajima et al. (1992) that obtained an outgoing longwave radiation (OLR) limit of 280–300 W m⁻². Therefore, water photolysis and H escape are both significant processes that impact the atmospheric composition and mass as shown by Schaefer et al. (2016). H escape may lead to the buildup of O₂. The planet GJ 1132b, studied by Schaefer et al. (2016), could be trapped in a long-term MO, so water photolysis followed by H escape could be important processes. Goldblatt et al. (2013) studied the transition to runaway greenhouse and obtained a radiation limit of 282 W m⁻² for modern Earth. This limit, however, is obtained for T_s < 1800 K and breaks down when the surface becomes hot enough to radiate in the visible wavelength regime, and then the OLR is seen to rise.

Usually in these models, the atmosphere is assumed to be in global energy balance, i.e., OLR or the emitted radiation by the planet is equal to the absorbed shortwave incoming radiation from the host star. This condition, however, might not be satisfied for cases with a deep MO lying below a massive steam atmosphere formed as a result of an additional heat flux from the interior. Such scenarios require more detailed investigations (e.g., Hamano et al. 2013).

In the case of coupled atmospheric–interior models, the evolution of the atmosphere and thermal cooling of the MO are invariably linked to each other and are influenced by several assumptions in the numerical modeling. First of all, there is the assumption of a gray or nongray radiative transfer approach in the atmospheric models. For example, various studies linking the atmosphere–interior model have used a gray radiative transfer approach (Matsui & Abe 1986; Elkins-Tanton 2008; Lebrun et al. 2013), whereas Hamano et al. (2015) and Schaefer et al. (2016) have used a nongray approach. A second assumption regards the composition of the atmosphere. A radiative-convective 1D atmospheric model assuming pure water composition, and based on the gray approach of Nakajima et al. (1992), has been coupled with an interior model by Hamano et al. (2013). On the other hand, both H₂O and CO₂ have been assumed as the main radiative species by Matsui & Abe (1986), Elkins-Tanton (2008), and Lebrun et al. (2013). It has been shown by Elkins-Tanton (2011) that an MO with as low as 0.1 wt% (1000 ppm) of water can potentially outgas hundreds of bars of water. Moreover, the formation of water oceans on rocky planets (e.g., due to the collapse of a steam atmosphere) suggests that steam (H₂O) is the major volatile reservoir for most planetary surfaces and their mantles. Therefore, in this paper we use only water as the constituent of the atmosphere.

Finally, various parameterizations in the interior model also affect the thermal cooling of the MO. For example, the cooling of the MO might also get influenced by factors such as the initial volatile content in the molten mantle (Salvador et al. 2017). The influence of initial volatile water content (varied from 0.1 to up to 1000 Earth ocean masses M_EO) on the solidification timescale of an MO has been explored by Schaefer et al. (2016), with a focus on the atmospheric escape of H and O. Other factors affecting the delay of the MO are the viscosity of the fluid and the depth of the MO (Solomatov 2007; Lebrun et al. 2013), as well as assumptions in the melting and solubility curves, which are explored in more detail in Nikolaou et al. (2019).

During the MO solidification phase, the growth and evolution of an atmosphere also affect the planet's thermal spectra significantly. In order to study the spectral evolution of hot molten planets with up to 50 Earth oceans (M_EO) of water, a nongray radiative transfer model was developed by Hamano et al. (2015). A planet with bulk initial water content exceeding 1 wt% would outgas large volumes of volatile into the atmosphere. Their work describes the growth of a massive steam atmosphere during the MO solidification phase, and their calculations have indicated that the blanketing effect of a "water-dominated" atmosphere prevents the thermal radiation from escaping the planet and prolongs the MO solidification timescale to ∼4 Myr for type I planets and 100 Myr for type II planets (Hamano et al. 2013). The planets, designated as type I and type II by are located at a critical distance of 1 and 0.7 au from the host star, respectively.

Very recently, Bonati et al. (2019) have discussed the possible direct detection of MO planets in nearby young stellar associations with future facilities, for example, E-ELT⁵ and LIFE.⁶ As shown by them, the young and close stellar objects (Table 1 of Bonati et al. 2019) with highest frequencies of giant impacts (e.g., β Pictoris having an age of 23 Myr) leading to the formation of MOs have a higher probability for detection of MOs. For a steam atmosphere, only a few planets with a flux density of the order of 10⁻² to 10⁻³ (μJy) would be detectable in the β Pictoris stellar association with an observation time of 50 hr with E-ELT through direct imaging (Bonati et al. 2019). The PLATO⁷ mission, on the other hand, will provide accurate ages to a precision of 10% and radii to a precision of 2% for a large sample of planetary systems (Rauer et al. 2014), which will lead to better characterization of young MO planets. Additionally, accurate estimation of the brightness temperature of a planet and the contrast ratio between planet and star can determine whether the atmospheric signatures are detectable or not (Lupu et al. 2014; Hamano et al. 2015; Marcq et al. 2017). Such studies motivate the investigation of exoplanetary atmospheres through the emission and transmission spectra, which can provide information on the planet's subsequent evolution to habitable conditions.

During the Hadean period, effects such as photolysis of water molecules, the atmospheric escape of H due to the high stellar irradiation, and thermal dissociation of H₂O due to the high surface temperatures will lead to the formation of additional species. In the latter case, the atmosphere in the lower layers close to the surface would mainly consist of other H- and O-bearing species, e.g., H₂ and O₂. The buildup of abiotic O₂ under such conditions has also been discussed in detail by Kasting (1995) and Schaefer et al. (2016). In this paper, we investigate the possible abiotic buildup of O₂ and other O-bearing species due to thermal dissociation of water. We also show the contribution of these newly formed species to the transmission spectra as a function of wavelength.

Table 1 provides a comparison of different coupled interior–atmospheric models from the literature, along with the methods used for radiative transfer calculations in each of them.

With the growing need for high-resolution infrared (IR) and microwave spectra to compare with the observations of the planets, line-by-line (lbl) modeling of atmospheric radiative transfer is essential. Therefore, in this paper, we utilize the lbl radiative transfer model GARLIC (Schreier et al. 2014; discussed in detail in Section 2.3) to calculate the emission from a steam atmosphere, which is formed as a result of outgassing from Earth's interior. The main motivation of our study is to investigate the impact of a time-dependent outgassing of volatiles on the thermal and chemical evolution of a steam atmosphere during the MO phase of the early Earth. Our calculation of the atmospheric radiative flux for H₂O without the radiative effects of H₂O speciation into H₂ and O₂ is used by the companion paper to assess the changes in the surface temperature and to obtain the solidification timescale of the MO. Also, the companion paper discusses the effects of several parameters that affect the MO cooling process. The evolution of surface temperature as a function of time is presented in the companion paper, as well as this paper. Finally, we study the effects of thermal dissociation of water on the OLR and the transmission spectra.

A brief outline of the paper is as follows. Section 2 describes the numerical models and the framework. Section 3 shows results of the OLR as a function of a prescribed [T, p] grid, which are an input to the interior model. The results of both models are combined together and presented in Section 4. In Section 5.1, the results of the wavelength-dependent effective height of the atmosphere for the whole duration of MO are presented. Section 6 discusses the effects of thermal dissociation of water on the OLR and the effective height of the atmosphere. Finally, we provide a discussion in Section 7, followed by a summary and conclusions in Section 8.

2.1. Interior Model

Thermal cooling of the planetary interior and the outgassing of volatile species (H₂O here) are obtained with the use of a 1D parameterized interior model. The model consists of a spherically symmetric mantle divided into layers of 1 km thickness assumed for a global MO stage for an Earth-sized planet. To ensure global mantle melting, we assume an initial mantle potential temperature of 4000 K. However, this can vary with different mantle melting curves. For the melting curves we use, the mantle is molten by assuming a potential temperature (T_p) of 3000 K and above. We start our calculations at T_p = 4000 K to be consistent with former studies (Lebrun et al. 2013; Schaefer et al. 2016) that calculate solidification timescales of the MO. The vaporization of minerals at such high temperatures leading to the formation of "silicate clouds" is not included since the phase during which they were present was relatively very short, ∼1000 yr (Zahnle et al. 2010). However, the vaporization of silicate minerals is important for estimating the atmospheric type in terms of, e.g., temperature and the relative amounts of vaporized metal and mineral species as shown by Miguel et al. (2011), for example, via cloud formation. However, more extensive laboratory data for the equilibrium composition of volatilized silicate minerals in hot steam atmospheres would be desirable.

The hot MO emits a heat flux, F_MAG, from the interior via the planetary surface owing to turbulent convection. The temperature profile, i.e., the variation of surface temperature (T_s) with pressure (p) in the mantle, is assumed to be adiabatic. The melt fraction ϕ(r) in the interior is calculated using the equation (Solomatov 2007)

$$\begin{eqnarray}&&\phi (r)=\displaystyle \frac{{T}_{{\rm{s}}}(r)-{T}_{\mathrm{sol}}(r)}{{T}_{\mathrm{liq}}(r)-{T}_{\mathrm{sol}}(r)},\end{eqnarray} \tag{ 1 }$$

where T_sol and T_liq are the solidus and liquidus temperatures, respectively. The critical melt fraction is taken to be ϕ_c = 0.4 (Lebrun et al. 2013; Schaefer et al. 2016) and divides the mantle into two zones. If ϕ > ϕ_c, the melt has a liquid-like behavior, and if ϕ < ϕ_c, the melt exhibits a solid-like behavior, indicating the end of the MO phase and the simultaneous volatile outgassing from it. The critical value of the melt fraction separates the liquid-like magma from the solid-like mantle phase, and the limit between the two zones is known as the rheology front (RF). The mantle temperature corresponding to the RF at the surface, denoted as T_RF,0, is a function of specific melting curves, explained in more detail in the companion paper.

The outgassing of water as a function of temperature is calculated using the 1D interior model assuming phase equilibrium with the melt and using the solubility curve of vapor for a given water concentration in the silicate melt (Caroll & Holloway 1994). Comparing the interior temperature profile with the melting curves of the mantle, adopted by assuming a non-evolving KLB-1 peridotitic composition (Zhang & Herzberg 1994; Herzberg et al. 2000; Hirschmann 2000; Fiquet et al. 2010), yields the volume fraction of the mantle in the solid and the liquid phases (e.g., Nikolaou et al. 2019). The fractionation of the water vapor into the two respective reservoirs is then calculated, and the amount of water that is supersaturated in the melt is released into the atmosphere. Thus, the pressure of outgassed water vapor at each stage of the simulation is obtained. The MO stage comes to an end when the potential temperature reaches T_RF,0 at the surface, i.e., 40% melt at the surface. For the melting curves employed in this study, this happens at T_RF,0 ∼ 1650 K.

By the use of the same 1D interior model, the convective heat flux out of the MO is calculated.

The initial water inventory chosen by us is slightly more than 1 Earth ocean (EO = 1.39 × 10²¹ kg). It is equivalent to 550 ppm or 400 bars. The motivation for this was to investigate the extent by which these starting conditions could reproduce Earth's ocean reservoir of 300 bars. Using this water inventory, Nikolaou et al. (2019) have obtained the outgassing of water from the mantle using an interior model coupled with a gray atmospheric model and calculated the variation of surface temperature with the outgassed pressure, denoted as [T, p] henceforth and shown in Figure 1. As an assumption, we prescribe their model output as the surface temperature and pressure for our atmosphere model. The range of outgassed pressure during the MO solidification corresponds to ${P}_{{{\rm{H}}}_{2}{\rm{O}}}\sim [4,300]$ bars (see Figure 1), which is also the range of the input chosen for the calculation of atmospheric structure and flux of the OLR, F_OLR, as explained below.

Zoom In Zoom Out Reset image size

2.2. Atmospheric Structure

In this section, we derive the atmospheric structure for the prescribed surface temperature and pressure as shown in Figure 1. The atmosphere is assumed to be consisting of only H₂O and is divided into N = 66 atmospheric layers equally spaced on a logarithmic scale of pressure and calculated by fixing the pressure at the top of the atmosphere (TOA) to be 1 Pa. The temperature structure of the atmosphere is obtained by assuming a dry adiabatic lapse rate for the atmospheric layers whose temperature is above the critical temperature of water (647 K). Since our calculations begin with a low pressure (4 bars) of water at higher temperatures, therefore, an assumption of the ideal gas equation for the calculation of dry adiabatic lapse is reasonable for our purpose. However, we additionally consider the T-dependence of C_p over the whole range of temperatures and pressures. The dry adiabatic lapse rate is calculated from the surface to the top by integrating the equation $d\mathrm{ln}(T)/d\mathrm{ln}(P)\,=R/{C}_{p}(T)$, where R is the gas constant and C_p(T) is the specific heat capacity of water with a temperature dependence, which is obtained using the following equation taken from Wagner & Pruß (2002) (and also mentioned by Hamano et al. 2013):

$$\begin{eqnarray}&&\displaystyle \frac{{C}_{p}(T)}{R}=1+{n}_{3}+\displaystyle \sum _{i=4}^{8}{n}_{i}\displaystyle \frac{{\left({\gamma }_{i}\tau \right)}^{2}\exp (-{\gamma }_{i}\tau )}{{\left[1-\exp (-{\gamma }_{i}\tau )\right]}^{2}},\end{eqnarray} \tag{ 2 }$$

where τ = T_c/T and T_c is the critical temperature (647 K). The values of coefficients n_i and γ_i are taken from Table 6.1 of Wagner & Pruß (2002).

The remaining procedure for calculating the atmospheric structure is the same as that of Goldblatt et al. (2013). The tropopause temperature T₀ for high surface temperature runs results from the moist adiabatic lapse rate when the pressure at the TOA reaches 1 Pa. T₀ for our runs ranges from 352 to 204 K (for decreasing surface temperatures) and reaches 204 K at T_s = 2900 K. For departures from ideal gas behavior, a self-consistent derivation of moist adiabat requires the use of a non-ideal gas equation at high densities (Wordsworth & Pierrehumbert 2013b).

The temperature profiles for the three [T, p] pairs are shown in Figure 2. The top panel illustrates that for a high temperature (T_s = 4000 K) a dry troposphere is present with no moist layer. However, as the surface temperature decreases (T_s = 1800 and 800 K; middle and bottom panels, respectively), moist layers of thickness 100 and 200 km are seen, respectively. The moist upper layers are cooler than the lower layers; hence, they act as a cold trap, and any exchange of heat from the TOA to space becomes less efficient.

Zoom In Zoom Out Reset image size

2.3. Atmospheric Radiative Transfer

For a gaseous, nonscattering atmosphere assuming local thermodynamic equilibrium, the intensity (radiance) at the TOA is given by the integral of the Schwarzschild equation

$$\begin{eqnarray}&&I(\nu )={I}_{{\rm{b}}}(\nu ){e}^{-{\tau }_{{\rm{b}}}(\nu )}+{\int }_{0}^{{\tau }_{{\rm{b}}}(\nu )}B(\nu ,T(\tau ^{\prime} )){e}^{-\tau ^{\prime} }\ d\tau ^{\prime} ,\end{eqnarray} \tag{ 3 }$$

where I_b is the background radiation and B is Planck's function for a blackbody with temperature T.

The monochromatic transmission ${ \mathcal T }$ and optical depth τ (at position s = 0, equivalent to τ = 0) are functions of wavenumber ν and path length s and are described by Beer's law as follows:

$$\begin{eqnarray}&&{ \mathcal T }\left(\nu ,s\right)={e}^{-\tau (\nu ,s)}.\end{eqnarray} \tag{ 4 }$$

In order to solve for the transmission of radiation (Equation (4)), it is important to obtain the optical depth τ, which is written as

$$\begin{eqnarray}&&\tau (\nu ,s)={\int }_{0}^{s}\alpha (\nu ,s^{\prime} )\ {ds}^{\prime} ,\end{eqnarray} \tag{ 5 }$$

where α is the absorption coefficient. It is usually written as the sum of the absorption cross sections κ_m of the molecule, scaled by its number density n_m as follows:

$$\begin{eqnarray}&&\alpha (\nu ,s)=\displaystyle \sum _{m}{\kappa }_{m}\left(\nu ,p(s),T(s)\right){n}_{m}(s).\end{eqnarray} \tag{ 6 }$$

For high-resolution lbl models, the absorption cross section of molecule m is given by the superposition of many absorption lines l, each described by the product of a temperature-dependent line strength S_l and a normalized shape function g describing the broadening mechanism. The calculation for the absorption coefficient and the absorption cross section for the molecules is described in the next section.

Assuming the net energy balance of the incoming and outgoing radiation at the TOA, along with the cooling flux from the MO, one can write

$$\begin{eqnarray}&&{F}_{\mathrm{OLR}}-{F}_{\mathrm{SUN}}={F}_{\mathrm{MAG}}.\end{eqnarray} \tag{ 7 }$$

Here F_OLR is the flux of the OLR calculated using GARLIC at a viewing angle of 38° from zenith (See Section 3 and Appendix A) and F_SUN = (1 − α)S₀/4, where S₀ is the solar constant with 100% luminosity and a fixed surface albedo of a = 0.3 is used. For simplicity, we do not include the effect of clouds and Rayleigh scattering.

2.4. GARLIC

The "Generic Atmospheric Radiation Line-by-Line Infared Code," or GARLIC (Schreier et al. 2014), finds its applications in several fields such as exoplanet studies (Rauer et al. 2011; Hedelt et al. 2013; Vasquez et al. 2013b, 2013c) and for Venus-type atmospheres (Hedelt et al. 2011; Vasquez et al. 2013a). GARLIC has been recently used to study Earth's atmosphere (Hochstaffl et al. 2018; Xu et al. 2018) for retrieval purposes in the field of remote sensing. Moreover, GARLIC has been validated using the transmission spectroscopic studies involving ACE-FTS IR observations of Earth (Schreier et al. 2018b). Furthermore, an intercomparison of radiance obtained from the three lbl models: GARLIC, KOPRA, and ARTS, in 18 channels throughout the mid-IR is presented in Schreier et al. (2018a).

In this work, the line absorption coefficients for molecule H₂O and other species are calculated by GARLIC using the 2010 HITEMP line list, i.e., the high-temperature molecular spectroscopic database by Rothman et al. (2010), and a "CKD continuum" model (Clough et al. 1989) for H₂O. The latest spectral database, HITRAN 2016 (Gordon et al. 2017), is also used toward the end for comparison. A Voigt profile (Schreier et al. 2014), representing the combined effect of pressure broadening together with Doppler broadening, is used.

The high-resolution spectra in GARLIC are obtained by construction, i.e., the wavenumber grid spacing is adjusted to the line widths δν = γ/n, where n is by default 4 and γ is the half-width at half-maximum dependent on the pressure, temperature, and molecular mass. Since the Lorentz width is proportional to the pressure, lines become narrower in the upper atmosphere. In GARLIC, the final grid (i.e., for all layers and molecules combined) is defined by the smallest δν.

We have verified the GARLIC code for a pure steam atmosphere at a low temperature of T = 300 K with a background surface pressure of 1 bar with the published results of Goldblatt et al. (2013) and Schaefer et al. (2016). A spectrally averaged OLR of 282.4 W m⁻² is obtained, which is close to the result of Schaefer et al. (2016), i.e., 274.4 W m⁻² as seen in Figure 3. It is also in a good agreement with Goldblatt et al. (2013), who obtained 282 W m⁻² for modern Earth using the lbl code SMART (Meadows & Crisp 1996).

Zoom In Zoom Out Reset image size

At high surface temperatures and pressures and for a fixed tropopause temperature of 344 K, we compared the OLR values with the former studies using lbl codes. The comparison is presented in Table 2. The values are comparable for T_s = 1000 K, but for higher temperature T_s = 2000 K, there seems to be a large deviation. The difference at high surface temperatures is attributed to the choice of the water continuum and due to the fact that the continuum is poorly understood at high temperatures.

Table 2.  Comparison of OLR for High-temperature Runs Calculated Using GARLIC with the Former lbl Studies

T0	Ts	Ps	OLR
(K)	(K)	(bars)	(W m−2)
344	1000	100		795.9a	801.4b
344	1000	250	789c		784.54b
344	2000	250	969c		1400b

Notes.

^aSchaefer et al. (2016). ^bThis study. ^cGoldblatt et al. (2013).

Download table as:  ASCIITypeset image

2.5. Effective Height

The effective height of the atmosphere, which is a function of the wavenumber, describes the extent of the atmosphere and is directly measured from the transit observations. For calculating the effective height from the transmission spectra, we use the following equation:

$$\begin{eqnarray}&&h(\nu )={\int }_{0}^{\infty }\left(1-{ \mathcal T }(\nu ,{z}_{i})\right){{dz}}_{i},\end{eqnarray} \tag{ 8 }$$

where the integral is calculated over all the limb transmission spectra with tangent altitude z_i for every incident ray i passing through the atmosphere terminating at the TOA. For more details, see Schreier et al. (2018b).

2.6. Chemical Equilibrium with Applications (CEA)

At T_s > 2000 K, water is expected to undergo thermal dissociation, leading to a somewhat different composition of the atmosphere that may influence the OLR and thereby the MO lifetime. Therefore, we employ the NASA CEA (Chemical Equilibrium with Applications) model⁸ to calculate the chemical equilibrium compositions and properties of various species formed as a result of the thermal dissociation of water at such high temperatures. The CEA calculates the chemical equilibrium concentrations of species from any set of initial composition(s) for a given [T, p]. Some of the built-in applications of the program include calculation of theoretical rocket performance, Chapman–Jouguet detonation parameters, shock tube parameters, and combustion properties (McBride et al. 2002; Snyder 2016). The CEA program utilizes a Gibbs free energy minimization routine to obtain the equilibrium-state composition for a given [T, p]. The thermodynamic data for the species are obtained from the NASA CEA database (Snyder 2016).

The scale height of the atmosphere also depends on the composition of the atmosphere and is calculated using the same procedure as described in Section 2.2 but with an updated mean molecular weight ${\mu }^{* }={f}_{{{\rm{H}}}_{2}{\rm{O}}}{\mu }_{{{\rm{H}}}_{2}{\rm{O}}}+{f}_{{\rm{H}}}{\mu }_{{\rm{H}}}+{f}_{{{\rm{H}}}_{2}}{\mu }_{{{\rm{H}}}_{2}}$ + ${f}_{{\rm{O}}}{\mu }_{{\rm{O}}}\,+{f}_{{{\rm{O}}}_{2}}{\mu }_{{{\rm{O}}}_{2}}+{f}_{\mathrm{OH}}{\mu }_{\mathrm{OH}}$, where the additional species are formed as a result of thermal dissociation of water.

Figure 4 shows a comparison between the temperature profiles of water and the CEA composition for a surface temperature of 4000 K and a surface pressure of 4 bars. Only a minor difference in the scale height is observed while comparing CEA composition with water. Also, the effect of the CEA composition on the OLR is small and is discussed later in Section 6.

Zoom In Zoom Out Reset image size

2.7. Framework

The planetary thermal emission is obtained for 400 [T, p] points as described in Section 2.7. These are calculated using the lbl code GARLIC for a down-looking observer at a viewing angle of 38° with the zenith angle as also chosen by Segura et al. (2003), Rauer et al. (2011), and Grenfell et al. (2011). This angle is considered to be the "mean angle" or "equivalent latitude" at which the solar insolation incident upon Earth is 341 W m⁻², which is equal to one-quarter of the present solar constant S₀ = 1361 W m⁻². The solar radiation, F_SUN, in our model is taken to be one-quarter of S₀(1 − α), which justifies the OLR calculation at a 38° viewing angle.

The flux (integral of radiance over the solid angle) in the units of W m⁻² μm⁻¹ at each [T, p] as a function of wavelength (ranging from mid-UV to IR) is shown in Figure 5 for a 100% steam atmosphere and for four [T, p] pairs. The effective radiating temperature, defined as ${T}_{\mathrm{rad}}=\sqrt[4]{{F}_{\mathrm{OLR}}/\sigma }$, is obtained, where F_OLR is calculated by integrating the flux over 5000 points of wavenumber. Using Planck's law, the electromagnetic flux corresponding to T_rad at the radiative emitting (top) layer is shown in Figure 5 (dashed line).

Zoom In Zoom Out Reset image size

For a [T, p] pair with high temperature (T_s = 4000 K) and low pressure (P_s = 4 bars), as in the top panel of Figure 5, the peak of the maximum emission is at short wavelength (following Wien's displacement law). Hence, the radiation can escape through the window regions at short wavelengths. At lower surface temperatures (lower panels), the amount of water vapor in the atmosphere increases (because of increased outgassing), leading to a decrease in the emitted flux. Also, the peak of the emission moves toward longer wavelengths (following Wien's displacement law), where radiation is effectively absorbed by the atmosphere. Due to this, the spectroscopic window region becomes optically thick, e.g., at 1 μm. As a consequence, the radiation drops by ∼10 orders of magnitude from the top panel (4000 K, 4 bars) to the bottom panel (800 K, 300 bars).

A comparison between the effective radiative temperature and the surface temperature of early Earth for variable outgassing and a fixed surface pressure of 300 bars (e.g., Kasting 1988; Marcq 2012; Marcq et al. 2017) is shown in Figure 6.

Zoom In Zoom Out Reset image size

We notice here that the effective radiating temperature T_rad is less than the surface temperature T_s. This is because of the greenhouse effect of water. Since we take into account a variable outgassing ranging between 4 and 300 bars of water, the effective radiative temperatures in the range of 265 K < T_rad < 2300 K are higher as compared to the fixed surface pressure (300 bars), where the range is much smaller, i.e., 265 K < T_rad < 1600 K. For example, Marcq (2012) obtained an effective temperature of 350–400 K for T_s = 2350 K for a fixed water pressure of ∼300 bars, whereas we obtain an effective temperature of 600–700 K by considering variable pressure. In the latter case, an efficient cooling of the planet is expected. A "blanketing effect," i.e., obscuring of the surface due to steam lying in the atmosphere, begins at ∼1800 K. For T_s < 1800 K, the planet is effectively radiating at T_rad = 265 K and a fixed OLR of 282 W m⁻². This means that the TOA's thermal spectrum is bounded by the flux from an effective radiation temperature of 265 K in this region.

Using the OLR calculated for the 400 grid points of [T, p], a bilinearly interpolated grid is calculated (Nikolaou et al. 2019) to facilitate an effective coupling of the atmospheric model with the interior model within a continuous range of outgassed pressures [4, 300] bars as a function of T_s [800, 4000] K, based on the actual [T, p] output of the interior model. The output of the interpolation is shown in Figure 7. We then investigate the coupled interaction of the interior and atmosphere in Section 4.

Zoom In Zoom Out Reset image size

The calculation for the MO evolutionary model (described in detail in Nikolaou et al. 2019) begins at a very high surface temperature of 4000 K, which corresponds to a molten state of the full mantle. The time period of formation of Earth is assumed to be 100 Myr. Therefore, to account for the calculation of the onset of thermal cooling of the mantle, 100 Myr is subtracted from the time of evolution in our model.

By assuming a radiation balance at the TOA, an iterative loop (starting at an initial time step) with a convergence criteria is run until the following fluxes are equal, meaning that they numerically converge to the same value:

$$\begin{eqnarray}&&| [{F}_{\mathrm{OLR}}(T,p)-{F}_{\mathrm{SUN}}]-{F}_{\mathrm{MAG}}| \lt \epsilon ,\end{eqnarray} \tag{ 9 }$$

where  = 0.1 W m⁻² is equal to the uncertainty in the OLR values obtained by using GARLIC. Once the iterative loop in Equation (9) is satisfied at a particular time step, the surface temperature T_s is then obtained, and a new potential mantle temperature T_p is calculated at the next time step by using the MO evolutionary model (described in detail in the companion paper). If Equation (9) is not satisfied, the surface temperature is adjusted by a factor ΔT_s, the corresponding F_OLR is chosen from the interpolated grid (see Figure 7), and F_MAG is recalculated. The steps are repeated until the surface temperature reaches approximately 1650 K (i.e., T_RF,0), as mentioned before. Therefore, the calculations for the MO evolutionary model have been carried out until T_s ∼ 1650 K. Note that this value is different from the T_RF,0 of Hamano et al. (2013), i.e., ∼1370 K, due to different assumptions in their RF. We do not perform calculations for the case of a warming planet, i.e., when a negative left-hand side of Equation (9) is obtained. The case for the latter is designated as planet type II by Hamano et al. (2013).

The time evolution of T_s is shown in Figure 8. As seen in the top panel of Figure 8, initially both potential temperature T_p and surface temperature T_s start to cool efficiently, but at different rates, and reach T_s ∼ 2150 K after 0.25 × 10⁴ yr. Thereafter, cooling slows down because the atmosphere starts to limit the escape of heat flux from the interior. This is evident from the rise of outgassed steam pressure as seen toward the right of the top panel of the figure. It can be clearly seen that the outgassed steam pressure rises from 4 bars to about 300 bars during the evolutionary phase of the MO. Finally, the MO solidification is achieved at ∼1 Myr. This is similar to the time obtained by Hamano et al. (2013) when they set the surface solidification temperature of T_RF,0 ∼ 1700 K, close to our value of 1650 K.

Zoom In Zoom Out Reset image size

Furthermore, in the bottom panel of Figure 8, it is seen that F_OLR approaches the radiation limit of ∼282 W m⁻² only by the end of the MO phase. A reduction in the F_MAG (=F_OLR − F_SUN) with time indicates thermal cooling of the MO.

Figure 9 shows the spectral response of the thermal emission for different times elapsed during MO evolution. For smaller timescales in early evolutionary stages (i.e., red solid line), due to high surface temperature and less surface pressure of water vapor in the atmosphere, the thermal radiation is quite high and escapes to space efficiently. However, these thermal fluxes rapidly decrease with time owing to the decrease in the surface temperature and an increasing surface pressure at later times ∼10³–10⁴ yr. This leads to an increased optical thickness of the atmosphere (due to water outgassing) and a "closing" of the associated IR absorption window.

Zoom In Zoom Out Reset image size

As illustrated in this figure, thermal fluxes at later times ∼10⁴–10⁶ yr tend to become constant for λ > 10 μm, which implies that the steam atmosphere is lying in the radiation limit. This is also evident from the bottom panel of Figure 8, where F_OLR is seen to become constant at a value of 282 W m⁻² (OLR limit).

The largest contribution to the emission spectra of the hottest scenarios is from the shorter wavelengths (U band) because of the absence of strong water absorption lines. In this case, a detailed diagnostic measure of a planet's atmosphere can be obtained from its thermal emission as seen in the visible and near-IR wavelengths, i.e., M, L, and K_s bands (e.g., water absorption bands in the 1–8 μm region). With the decrease in T_s and growth of the steam atmosphere, the U band's contribution to the OLR decreases with time, whereas the M, L, and K_s bands initially decrease and then become constant with time, indicating that the OLR limit has been reached.

Notably, the spectral band feature at ∼6.3 μm in Figure 9, related to the absorption by water vapor in the atmosphere, is seen to broaden (because of pressure broadening) with the increase of the amount of outgassed steam, suggesting that the atmosphere is saturated with water.

5.1. Transmission Spectra

Transmission spectroscopy is a widely used tool to constrain the composition of planetary atmospheres. For exoplanetary studies, it has been recently used in this way to provide the first hints on bulk atmospheric composition and clouds of Earth-size planets, for example, Trappist-1b and c (Wit et al. 2016), water bands in Neptune-sized exoplanets HAT-P11b (Fraine et al. 2014) and HAT-P-26b (Wakeford et al. 2017), and the H₂-rich atmosphere of planet GJ 3470b (Ehrenreich et al. 2014).

Using our 1D atmospheres and the radiative transfer code, we calculated the transmission spectrum ranging from 0.34 to 30 μm. It is shown as the wavelength-dependent effective height in Figure 10 for the whole duration of the MO. The spectra are computed at a very high resolution, as explained in Section 2.4 using the HITEMP 2010 database, and then binned to a resolution of λ/Δλ = 1000. As the MO solidifies over time with the decrease in the surface temperature (see Figure 8), the effective height of the atmosphere is reduced. Note that the abundance of water molecules for a different elapsed time of the MO is the same. The depth of the atmosphere where the absorption by the water vapor takes place is found to be ∼100 km, in contrast to 10–20 km for modern Earth (Wunderlich et al. 2019).

Zoom In Zoom Out Reset image size

As a comparison, the results of the effective height using the latest database HITRAN 2016 (Gordon et al. 2017) are presented in Figure 11. The variation in the absorption features in the two cases is due to the fact that HITEMP 2010 consists of many more molecular bands and line transitions at high temperatures. HITRAN 2016 (Gordon et al. 2017), on the other hand, has an expanded database of water vapor as compared to the previous HITRAN versions. The difference in the two databases is appreciable at high temperatures and low pressures and looks the same at T_s < 2100 K as illustrated in Appendix B.

Zoom In Zoom Out Reset image size

For the early stages of the MO, it can be expected that, due to thermal dissociation of H₂O at T_s > 2000 K, additional species such as H, O, O₂, O, and OH are formed. The molar fraction of these species as a function of temperature is shown in Figure 12. For this example, we set P_s = 4 bars (corresponding to T_s = 4000 K as per the interior outgassing model). It can be seen that at T_s = 4000 K only 13% of H₂O is present, with the rest dissociated into H (32.9%), H₂ (17%), O (16%), OH (15%), O₂ (4.7%), and other trace elements. Moving to lower temperatures, the molar fraction of the other species reduces and water constitutes 100% at T_s ∼ 2000 K. Note that there is a peak in H₂ and O₂ fractions at temperatures around 3500–4000 K. This is potentially relevant for assessing O₂ false positives in atmospheric biosignature science.

Zoom In Zoom Out Reset image size

The molar fraction as a function of temperature obtained using the NASA CEA program is verified with another equilibrium chemistry codes, e.g., IVTANTHERMO (L. Schaefer, personal communication), which is also used in Schaefer & Fegley (2004) and recently in Schaefer & Fegley (2017).

The higher the atmosphere, the more chances there are to detect it through the pioneering transit observations. Hence, a planet in the early stages of MO would have perhaps more chances of detection through characterization of its composition.

Atmospheric concentration profiles (altitude versus molar fraction) for various species and temperature–pressure combinations are shown in Figure 13. The surface temperature and pressure are according to the outgassing rates of the interior, as mentioned on the top of each panel in this figure. It is to be noted that the initial concentration of H₂O is set by the outgassing and is therefore dependent on the surface temperature.

Zoom In Zoom Out Reset image size

Figure 14 shows a comparison between the emission thermal spectra or OLR as a function of wavelength for a 100% steam atmosphere and the equilibrium chemistry composition (obtained using CEA) for T_s = 4000 K and P_s = 4 bars. The OLR for the CEA composition is found to be lower by 5%–6% than the pure steam case. The species that dominate the CEA composition, namely, H₂ and H, might be responsible for lowering the OLR. However, the O-bearing species, namely, O₂, OH, and O, did not result in any change in the OLR value. The impact of CEA on the MO solidification time is hence expected to be negligible. More discussion related to absorption by H₂ is presented in Section 7.3.

Zoom In Zoom Out Reset image size

We also compare the effective height of an atmosphere constituting pure water versus the individual species of the CEA profile at T_s = 4000 K and P_s = 4 bars using the HITRAN 2016 spectral database as shown in Figure 15. The lowered fraction of H₂O in the CEA case results in a slightly higher effective height (cyan line) as compared with 100% abundance of H₂O (orange line) according to the equation H = kT/μg. Due to the thermal dissociation of water at higher temperatures, O- and H-bearing species contribute a significant amount to the effective height of the atmosphere. Different absorption behavior with altitude leads to spectral features arising at different altitudes as seen in Figure 15. One can see strong absorption features for the species, namely, for H₂O, H₂, and O₂, but the contribution of each species toward the transmission as a whole is not linearly additive. Due to this property, the effective height of the atmosphere for the CEA profile with all the species (red line) is seen to be overlapping with the CEA profile when considering only water as absorber in the lbl calculations in cyan. The difference in the spectra of H₂O of the CEA profile (cyan line) and spectra of pure H₂O atmospheres (orange line) is mainly due to the increase in the scale height of the atmosphere, so the overall effective height increases.

Zoom In Zoom Out Reset image size

One of the key objectives of the paper was to study the thermal spectral evolution of the early Earth, i.e., during the MO phase using an lbl atmospheric model. We have performed a systematic study of the spectral evolution for a steam-based atmosphere by including the time-dependent outgassing rates obtained from Nikolaou et al. (2019). We briefly summarize and discuss our results below.

7.1. Model Comparison

7.2. MO Solidification Times

7.3. H₂ as Greenhouse Gas

In this paper, we have investigated the evolution of the steam atmosphere during the MO period of the early Earth using the atmospheric lbl code GARLIC based on a non-gray approach. The evolution of the atmosphere is invariably linked to the thermal cooling of the MO. We assumed a time-dependent outgassing of water vapor, leading to a variable outgassed pressure ranging between 4 and 300 bars and surface temperatures between 800 and 4000 K. Using variable outgassed pressure instead of choosing a constant value, i.e., 270 bars, we determine higher estimates for the effective radiating temperature, as compared to the previous studies. In this paper, we have studied the atmospheric spectral evolution by calculating the OLR and the effective height of the atmosphere as a function of MO cooling timescales. The thermal dissociation of water at high temperatures and its effect on the OLR and effective height of the atmosphere are also presented. The main conclusions of our study are listed below.

(1) For relatively hot surfaces (T_s > 2500 K), a larger flux of radiation leaves the planet owing to a smaller pressure of water on the surface. This leads to a significant and rapid cooling of the MO. We find out that such atmospheres usually have high effective radiating temperatures (1000–2300 K) and hence are easier to characterize observationally.

(2) For relatively cool surface temperatures (T_s < 2500 K), a limited flux is able to escape to space, with moderate effective radiating temperatures (265–800 K), and a strong blanketing of the surface is seen. Such atmospheres are difficult to observe, probably due to the thick overlying layer of clouds.

(3) Due to the blanketing effect of water, only a limited amount of radiation can escape to space in order to cool the planet. This limiting value of the OLR, known as the OLR limit, is estimated to be 282 ± 2 W m⁻² using the lbl code GARLIC. This value is consistent with previous studies that reported a range between 280 and 310 W m⁻² (Abe & Matsui 1988; Goldblatt et al. 2013; Hamano et al. 2013; Kasting 1988; Kopparapu et al. 2013; Marcq et al. 2017). At this stage, the planet is in a runaway greenhouse regime, radiating at an effective temperature of 265 K. Moreover, we conclude that the difference between the surface temperature and the effective radiating temperature for a planet is a strong function of its surface pressure and is important for estimating the radiative heating due to the greenhouse gases present in the atmosphere.

(4) By coupling the results of the climate model (this work) and the interior model (companion paper), the solidification timescale for the MO is found to be 1 Myr at a solidification temperature of 1650 K, assuming a pure steam atmosphere derived from initial water reservoir slightly greater than one Earth ocean (400 bar). One of the main applications of our work would be to compare the results of the lbl model "GARLIC" (nongray) used in this study with the gray approach (Lebrun et al. 2013). This has been done by Nikolaou et al. (2019), who have found a delay in the MO cooling timescale by a few hundred thousand years by using results of our lbl (nongray) atmospheric model as compared to a gray model. Furthermore, they have also discussed several other factors responsible for the delaying of the MO timescale.

(5) The transmission spectra for the whole time duration of the MO for pure water atmospheres are presented as wavelength-dependent effective height yielding ∼100 km depth of the atmosphere, where the absorption by water molecule takes place. The effective height of the atmosphere decreases as the MO solidifies or the surface temperature reduces.

(6) Finally, we have included the thermal dissociation of H₂O at T_s > 2000 K calculated using chemical equilibrium analysis (CEA) and studied its effect on the OLR. We suggest an H₂-rich early Earth's atmosphere with ∼10%–30% volume fraction for surface temperature ranging between 3500 and 4000 K, respectively, i.e., during the initial evolutionary phase of the MO. It is responsible for lowering the OLR by 1%–6% as compared to the pure water case. The thermal dissociation also resulted in an abiotic formation of O₂. H₂ formed in the atmosphere may slowly diffuse through the atmosphere and hence be transported to the higher layers and escape. Photodissociation of H₂O, which is yet another important process for H₂ formation and escape in these atmospheres, has not been included in this study. For the CEA profile, we also studied the transmission spectra. In this spectra (Figure 15), the thermal dissociation of water into H- and O-bearing species led to an increase in the effective height of the atmosphere and can be clearly seen in the beginning of the MO (T_s = 4000 K and P_s = 4 bars), which decreases as the surface temperature cools down.

Our results provide important implications for observing hot molten and young planets having short MO phases, especially for planets with a low initial volatile inventory and with an evolving atmosphere. The presence of H₂O, H₂, and perhaps O₂ in early thermal spectra of the planet could suggest the evolutionary phase of the planet and also act as biomarkers for characterizing its atmosphere.

Moreover, by including CO₂ as another constituent of the atmosphere and studying its thermal dissociation, this work invokes further investigation related to chemistry in the upper atmosphere, leading to the possible formation of organics in the early Earth's atmosphere.

We thank the referee for useful comments that led to the improvement of the manuscript. N.K. acknowledges funding from the German Transregio Collaborative Research Centre "Late Accretion onto terrestrial Planets (LATP)" [TRR170, subproject C5]. M.G. acknowledges funding from DFG project GO 2610/1-1. A.N. and N.T. acknowledge funding from the Helmholtz Association (project VH-NG-1017). F.S. acknowledges DFG project SCHR 1125/3-1. N.T. also acknowledges support from the German Research Foundation (DFG) through the SPP 1833 "Building a habitable Earth" (project TO 704/2-1). N.K. thanks R. Wordsworth and L. Schaefer for providing OLR data for comparison in Figure 3 and Table 2. N.K. thanks S. Städt for discussions on observer's geometry setup for GARLIC.

In this appendix, we compare the OLR obtained using spectral databases HITRAN 2012 (Rothman et al. 2013) and HITEMP 2010 (Rothman et al. 2010). We also choose two different viewing angles, 0° and 38°, and show the comparison of OLR versus temperature in Figure 16. As can be seen, OLR is the same for the two databases up to T_s = 1800 K (Figure 16(a)) and different for higher temperatures (Figure 13(b)). We note that a lower OLR is obtained with the HITEMP 2010 database, as it includes more water lines at higher temperatures. The magnitude of the difference is shown in Figure 16(c), i.e., we find that the OLR difference between the two databases rises monotonically after T_s = 1800 K and the OLR for the HITRAN 2012 database is almost twice as high as that for the HITEMP 2010 database around T_s = 3000 K.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The comparison for the two databases is seen in Figure 17.