© ESO 2020

1 Introduction

Understanding how the early atmosphere of Earth emerged into habitable conditions is a central question not only for addressing our own origins, but also for interpreting the fascinating new data in exoplanetary science (for recent reviews, see, e.g., Jontof-Hutter 2019; Madhusudhan 2019). Magma oceans (MOs) consisting of hot and molten silicates in the mantle exist during the planetary accretion phase (Elkins-Tanton 2012) and facilitate the formation of an atmosphere by outgassing of volatile species such as H₂O, CO₂, H₂, CH₄, and CO. Thus, the MO period likely represents a key juncture on the pathway to habitability because atmosphere-interior couplings (e.g., Elkins-Tanton 2008; Lebrun et al. 2013; Hamano et al. 2013; Schaefer et al. 2016) established during this phase could have set the stage for subsequent atmospheric evolution. The duration of the magma ocean phase is potentially affected by the absence or presence of an atmosphere and the interior dynamics of the mantle (e.g., Lebrun et al. 2013; Nikolaou et al. 2019) that affect the thermal spectral evolution of the planet (e.g., Hamano et al. 2015; Katyal et al. 2019). Other uncertainties such as the initial volatile inventory of the mantle, for example, the amount of H₂O stored in the mantle (Genda et al. 2017; Meech & Raymond 2019), and compositional constraints such as the mantle oxygen fugacity (Hirschmann 2012; Schaefer & Fegley 2017) are also important quantities that might critically affect the oxidation state of outgassed volatiles and hence the subsequent climate and compositional evolution of the atmosphere.

The budget and origin of volatiles in the mantle of Earth is one of the main questions in geo-science (see, e.g., Meech & Raymond 2019). Recent studies suggest quite diverse scenarios that range from a wet, late accretion to a dry accretion to which water was added mostly during the “late veneer”, which is proposed to have occurred roughly around 80–130 Myr after the isolation of protoplanetary nebula at around 4.5 Gyr before present (Albarede et al. 2013). Abundances of highly siderophile elements (HSE) and other volatile species such as H₂O suggest thatvolatile-rich material was added to Earth after its core formed (Albarède 2009). Measurements of the abundances of water and carbon in the bulk silicate Earth (BSE) show that the late veneer might indeed have supplied 20–100% of the budget of the total hydrogen (H) and carbon (C) in the BSE (Wang & Becker 2013). Some recent studies do not support the late veneer as the main source of water on Earth (Dauphas 2017; Fischer-Gödde & Kleine 2017; Hopp & Kleine 2018) and argue that water was likely added to Earth already during the waning stages of accretion and thus was already present during the late veneer. Interestingly, the reaction of terrestrial water and Fe from the late veneer produced much hydrogen as a byproduct (Zahnle et al. 2020). There is increasing isotopic evidence indicating that the late veneer might also have supplied a significant amount of reducing material onto the surface of Earth and into the atmosphere (Genda et al. 2017; Greenwood et al. 2018; Zahnle et al. 2020).

The impact-degassing studies by Schaefer & Fegley (2007, 2010), Hashimoto et al. (2007) and Fegley & Schaefer (2012) also in general favored the outgassing of reduced gases upon impacts by certain types of reduced meteoritic materials such as ordinary and enstatite chondrites. Furthermore, changes in chondrites delivery rate are likely to affect the redox state of the mantle and hence the atmospheric composition and amount, as suggested recently by Schaefer & Fegley (2017). Their work also showed that the variation in oxygen fugacities for mixtures of primitive meteoritic material could lead to atmospheres ranging from highly oxidizing to highly reducing ones (also see Hirschmann 2012).

The redox state of the minerals and melts is related to the oxygen fugacity fO₂ of the system that is equivalent to the partial pressure of the gas or the availability of oxygen (Gaillard et al. 2015). This quantity is poorly constrained for the early Earth (Frost & McCammon 2008). The current understanding is that before the formation of the Earth core ~4.56 Gyr before present (Stevenson 1983), the lower fO₂ of the uppermantle was approximately three log units below the reducing iron-wüstite (IW) buffer, that is, IW-3 and it evolved to a higher fugacity value that was approximately equivalent to the current mantle of Earth, which resembles the quartz-fayalite-magnetite (QFM) buffer, that is, IW+3.5 (Wood et al. 1990; O’Neill 1991). Direct geological evidence regarding the redox state of the mantle is sparse, and it is not clear when the Earth may have been oxidized during the Hadean (Trail et al. 2011). However, the geological record mainly suggests that the upper mantle was initially reduced and became progressively oxidized between 4.6 and 3.9 Gyr before present (Kasting et al. 1993; Delano 2001; Schaefer & Fegley 2017), with a possible later further increase in redox state by ≈ 1.3 at the end of the Archean (Aulbach & Stagno 2016).

Previous interior modeling studies investigating atmospheres at and around the time of the magma ocean (e.g., Lebrun et al. 2013; Hier-Majumder & Hirschmann 2017; Nikolaou et al. 2019) have considered outgassing of volatiles such as H₂O and CO₂, which means that they considered the mantle to be constantly oxidized. These studies suggested outgassed pressures of H₂O and CO₂ ranging from a few tens to hundreds of bar depending upon uncertainties, for instance, in the initial volatile content, the timing of the MO, and internal properties of the mantle. Recent modeling studies have now started to investigate the effect of the redox state on the outgassing. Pahlevan et al. (2019) investigated the effect of fO₂ on the outgassing and atmosphere losses and provided evidence for an early oxidation of silicate Earth using the D/H of the oceans. Most recently, Ortenzi et al. (2020) calculated the expected outgassing rates for the reduced and oxidized rocky planets and focused on exploring the observational constraints for the atmospheric (and interior) redox state of exoplanets.

Hydrodynamic escape of H₂ during this time (~4.5 Gyr before present) was likely not sufficient to desiccate the planet (Kasting et al. 1993), so that conditions remained wet after the Moon-forming impact and crystallization of the final MO (Lammer et al. 2018). Hydrogen would therefore have been a sufficiently major component of an accretionary steam-based atmosphere during the MO phase. However, hydrogen evolution around the time of the MO is poorly understood because of the uncertainties in accretion, outgassing, ingassing, and escape (see, e.g., Tian 2015). Several studies have considered the energy-limited hydrodynamic escape of hydrogen from the early atmosphere (Archean) of Earth (Tian et al. 2005; Kuramoto et al. 2013; Zahnle et al. 2019) and calculated somewhat differing results based on including or excluding the effect of H₂ diffusive flux on the gas-density profile and processes such as radiative cooling and thermal conduction that these escape models treat, for example. Johnstone et al. (2019) have also discussed uncertainties in the energy-limited mass-loss formula when applied to the early atmospheres. For rocky planets around M dwarfs, large amounts of water could be lost depending on the uncertainties in the stellar luminosity, for instance (e.g., Tian & Ida 2015; Schaefer et al. 2016).

We apply a coupled suite of interior and atmospheric models to investigate outgassing and escape during the MO period and study the effect of varying the mantle fugacity on the atmospheric evolution. We also calculate theoretical atmosphericspectra that serve as a link for observations with early Earth-sized planets. In Sect. 2 we present the volatile speciation model and the atmospheric model. Section 3 presents the scenarios weadopted. Section 4 provides results of the effect of the redox state of the mantle on the infrared emission and transmission spectra. The interplay between outgassing and atmospheric loss of H₂ is also presented in this section. In Sect. 5 we discuss our findings. Finally, we provide a conclusion of our work in Sect. 6.

2 Methods and models

2.1 C-H-O based volatile speciation model

The outgassing and volatile chemical speciation were simulated following the equilibrium and mass balance method (French 1966; Holloway 1981; Fegley 2013; Gaillard & Scaillet 2014; Schaefer & Fegley 2017) that has recently been presented in detail by Ortenzi et al. (2020). We calculated the outgassed composition of the volatiles considering a broad range of temperatures, pressures, and redox states. The four common petrological buffers used in the literature and their mantle oxidation states are $$2\text{Fe}+{\text{SiO}}_2+{\text{O}}_2\rightleftharpoons {\text{Fe}}_2{\text{SiO}}_4,\text{QIF}$$(1) $$2\text{Fe}+{\text{O}}_2\rightleftharpoons 2\text{FeO},\text{IW}$$(2) $$3{\text{Fe}}_2{\text{SiO}}_4+{\text{O}}_2\rightleftharpoons 2{\text{Fe}}_3{\text{O}}_4+3{\text{SiO}}_2,\text{QFM}$$(3) $$\text{Ni}+1/2{\text{O}}_2\rightleftharpoons \text{NiO}.\text{NiNiO}$$(4)

Quartz-iron-fayalite (QIF) and IW buffers represent reduced conditions, and quartz-fayalite-magnetite (QFM) and nickel-nickel-oxide (NiNiO) describe the oxidizing redox states. The oxidation state (oxygen fugacity) for these buffers was calculated following the parameterization from Holloway et al. (1992), $${\log }_{10}{\text{fO}}_2=A-B/T+C(P-1)/T+Z,$$(5)

where the pressure (P) is in bars, the temperature (T) in is Kelvin, and the parameters A, B and C are defined in Table 1. Z is a positive or negative number denoting the deviation from the fugacity values with respect to the buffers, as stated in the table. We mainly used the IW buffer and chose Z to be −4 (highly reducing), 0 (reducing), and 4 (oxidizing) in order to investigate the effect of a range of redox mantle states on outgassing.

By simulating the oxidation state, that is, the oxygen fugacity fO₂ of the system, we are able to simulate the gas chemical speciation via the following equilibria: $$\text{CO}+1/2{\text{O}}_2\rightleftharpoons {\text{CO}}_2,\text{and}$$(6) $$2{\text{H}}_2+{\text{O}}_2\rightleftharpoons 2{\text{H}}_2\text{O}.$$(7)

To calculate the ratio between the carbon species (Eq. (6)), we considered the equilibrium constant K₁ for the equilibrium as $$K_1=\exp \left(\frac{-{\Delta }_r{G}_1^0}{RT}\right)=\frac{X_{{\text{CO}}_2}}{X_{\text{CO}}}\frac{1}{{\text{fO}}_2{}^{1/2}},$$(8)

where R is the universal gas constant (8.314 J K⁻¹ mol⁻¹), T is the temperature of the outgassed material in Kelvin, and ${\Delta }_r{G}_1^0$ is the Gibbs free energy of the reaction in Eq. (6), where $${\Delta }_r{G}_1^0={\Delta }_f{G}_{{\text{CO}}_2}^0-{\Delta }_f{G}_{\text{CO}}^0.$$(9)

Substituting ${\Delta }_r{G}_1^0$ and the calculated fO₂ in Eq. (6) gives us all the necessary parameters needed to calculate the carbon species. For Eq. (7), the Gibbs free energy of reaction ${\Delta }_r{G}_2^0$ is related only to the Gibbs free energy of formation of water (${\Delta }_f{G}_{{\text{H}}_2\text{O}}^0$), $${\Delta }_r{G}_2^0=2{\Delta }_f{G}_{{\text{H}}_2\text{O}}^0.$$(10)

The values of Δ_fG⁰ for the different species were calculated following the example of Fegley (2013) and compared to the literature (Chase 1999).

Similarly, the abundances of H₂ and H₂O Eq. (7) are related to the equilibrium constant K₂ and the fugacity fO₂ via $${\left(\frac{X_{{\text{H}}_2\text{O}}}{X_{{\text{H}}_2}}\right)}^2=K_2{\text{fO}}_2,$$(11)

where $X_{{\text{H}}_2\text{O}}$ and $X_{{\text{H}}_2}$ are the mole fractions of H₂O and H₂ and are related to the partial pressure of each of the species. The rate of outgassing of H₂O can therefore be related to the rate of outgassing of H₂ as $$r_{{\text{H}}_2\text{O}}=R_1{r}_{{\text{H}}_2},$$(12)

where $R_1={(K_2{\text{fO}}_2)}^{0.5}$, $r_{{\text{H}}_2}$, and $r_{{\text{H}}_2\text{O}}$ are the outgassing rates of H₂ and H₂O in units of m⁻² s⁻¹, respectively.

Fig. 1 Schematic showing coupled interior-atmosphere exchange via outgassing of reduced species such as H2 and CO (orange) for a reduced mantle (e.g., IW buffer) and oxidized species such as H2O and CO2 (blue) for an oxidized mantle (e.g., QFM buffer) as the magma ocean solidifies. The H2 in the atmosphere is lost to space via escape.

2.2 Coupled interior-atmospheric evolution

A schematic for interior-atmosphere coupling is presented in Fig. 1. We started by taking the variation in surface pressure p_s, surface temperature T_s, and volatile abundances $f_{{\text{H}}_2\text{O}}^{\text{init}}$ and $f_{{\text{CO}}_2}^{\text{init}}$ from the interior model output of Nikolaou et al. (2019), that is, a mantle temperature of 4000 K and an initial mantle composition of $X_{{\text{H}}_2\text{O}}=0.05$ weight percent (wt.%) (550 ppm) and $X_{{\text{CO}}_2}=0.01$ wt.% (130 ppm) as the input. We then applied the volatile speciation model from Ortenzi et al. (2020) as described in Sect. 2.1 to obtain the final outgassed molar abundance of the species as a function of oxygen fugacity fO₂ relative to the given mineral buffer (IW), as shown in Fig. 2. This figure shows that for a reduced mantle (between IW and IW-4) and increasing input H/C from top to bottom, we obtain atmospheres that are rich in CO, CO+H₂ mixtures, and H₂. On the other hand, for a more oxidized mantle with typically higher fugacity values (between IW and IW+4) and increasing input H/C from top to bottom (Fig. 2), we obtain an atmosphere that is rich in CO₂, H₂O +CO₂, and H₂O. The bottom panel of Fig. 2 corresponds to an initially assumed 100% H₂O which results in a almost pure H₂ atmosphere for a reduced mantle and a pure H₂O atmosphere for an oxidized mantle.

The initial outgassed volatiles H₂O and CO₂ taken from Nikolaou et al. (2019) with a mean molar mass μ_v react with the melt and result in a different composition of the species outgassed to the atmosphere as calculated from the speciation model (Sect. 2.1). When the number of moles of the species C-H-O is kept constant in the speciation model, a different atmosphere with a new mean molar mass μ_atm results, which has a different pressure at the bottom (surface) of the atmosphere, which we now call p_boA. This pressure is calculated based on the volatile mass balance (Bower et al. 2019), $$p_{\text{boA}}=p_{\text{s}}\left(\frac{{\mu }_{\text{atm}}}{{\mu }_v}\right),$$(13)

where p_s is the initial (surface) pressure of the volatiles calculated as a mono-gas atmosphere.

Starting from the mass of the outgassed species H₂O given by M = P*A∕g, we obtained the following equation for the initial H₂O outgassing rate from the interior. $$R_{{\text{H}}_2\text{O}}^{\text{init}}=\frac{1}{A}\frac{\text{d}N}{\text{d}t}=\frac{N_{\text{A}}1000}{g{m}_{{\text{H}}_2\text{O}}}\left(\frac{\text{d}P}{\text{d}t}\right)\text{molecules}{\text{m}}^{\text{-2}}{\text{s}}^{\text{-1}},$$(14)

where dN is the number of H₂O molecules outgassed in a time interval dt, N_A is Avogadro’s number (6.022 × 10²³), A is the area of the Earth surface (5.1 × 10¹⁴  m²), $m_{{\text{H}}_2\text{O}}$ is the molecular weight (18 g mol⁻¹), g = 9.81 m s⁻² is the surface gravity, and d p is the difference in outgassed pressure at a time interval dt taken from coupled model result of Katyal et al. (2019) and Nikolaou et al. (2019).

To obtain the outgassing rate of H₂, we used Eqs. (12) and (14). To account for the redox variation of the mantle and the volatile mass balance, we assumed that the initial rate of H₂O outgassing $R_{{\text{H}}_2\text{O}}^{\text{init}}$ in Eq. (14) is the total outgassing rate of combined (H₂ +H₂O). Therefore, $R_{{\text{H}}_2\text{O}}^{\text{init}}=r_{{\text{H}}_2}+r_{{\text{H}}_2\text{O}}$. Inserting Eq. (12) here, we obtain the rate of H₂ outgassing as $$r_{{\text{H}}_2}=\frac{R_{{\text{H}}_2\text{O}}^{\text{init}}}{1+R_1}\text{molecules}{\text{m}}^{\text{-2}}{\text{s}}^{\text{-1}}.$$(15)

Fig. 2 Volatile chemical speciation in terms of outgassed mole fraction vs. oxygen fugacity of the mantle. The x-axis indicates the oxygen fugacity range in logarithmic units relative to the IW buffer. For reference (in log units), buffers QIF = IW-1, QFM = IW+3.8, and NiNiO = IW+4.2 approximately (Wood et al. 1990). From top to bottom: BOM, active MO phase, EOM and a steam-atmosphere phase that are characterized by different initial H/C between hydrogen and carbon. The corresponding input values of the initial volume mixing ratio of volatiles, surface pressures Ps, and temperatures Ts (Nikolaou et al. 2019) for the speciation model are indicated at the top of each panel.

2.3 Convective lapse rate

The convection takes into account the heat transport by adiabatic expansion of a mixture of a nonideal condensable gas H₂O (because the critical point of water at T_c = 647 K and p_c = 220 bar is close to the T–p range considered in this study), denoted with subscript c, and ideal, noncondensable gases (mixture of CO₂, CO, and H₂) denoted with the subscript d.

For a saturated water vapor atmosphere, the pressure exerted by the condensable is given by the Clausius-Clapeyron equation, $$\frac{\text{d}{p}_{\text{s}}}{\text{d}T}=\frac{L_{\text{c}}(T)}{R{T}^2}{p}_{\text{s}},$$(16)

where p_s is the saturation vapor pressure and L_c(T) is the latent heat of vaporization of the condensable, that is, H₂O in this case. Inserting this expression for d lnp_c∕d lnT into the formula given in Pierrehumbert (2010) leads to the pseudoadiabatic slope as used in this study, which is given by (see also Ding & Pierrehumbert 2016) $$\frac{\text{d}\ln P}{\text{d}\ln T}=\frac{p_{\text{s}}}{P}\frac{L_{\text{c}}(T)}{R_{\text{c}}T}+\frac{p_{\text{d}}}{P}\frac{c_{\text{pd}}}{R_{\text{d}}}\frac{1+\left(\frac{c_{\text{pc}}}{c_{\text{pd}}}+\left(\frac{L}{R_{c}T}-1\right)\frac{L}{c_{\text{pd}}T}\right)r_{\text{sat}}}{1+\frac{L}{R_{\text{d}}T}{r}_{\text{sat}}},$$(17)

where c_pd and c_pc are the specific heat capacities of the noncondensable and condensable species, respectively. r_sat is the saturation mass mixing ratio, given as (Catling & Kasting 2017) $$r_{\text{sat}}=\varepsilon f_{\text{sat}}({\text{H}}_2\text{O})=\varepsilon \frac{p_{\text{sat}}}{P}.$$(18)

Here, ɛ is the ratio of mass of condensable and noncondensable species, P is the total pressure, and f_sat is the volume mixing ratio of condensable at saturation. Equation (17) reduces to Eq. (A.13) of Kasting (1988) at r_sat → 0. In this case, the atmosphere is dry and the temperature profile is obtained by the dry adiabatic lapse rate given by (from Eq. (17)) $$\frac{\text{d}\ln p_{\text{d}}}{\text{d}\ln T}=\frac{c_{\text{pd}}(T)}{R_{\text{d}}},$$(19)

where R_d is the gas constant for the noncondensable or the dry component of the mixture and is obtained by Eq. (A.7). When more than one noncondensable species is present, the mean molecular weight of the mixture of noncondensable is accounted for.

At temperatures that are high enough and for a condensable species reservoir much larger than the atmospheric pressure exerted by the dry gases such that r_sat ≫ 1, P → p_sat and p_d ∕P → 0, Eq. (17) reduces to the Clausius-Clapeyron relation (Pierrehumbert 2010), which is valid for an ideal gas equation. The temperature dependence of the latent heat of vaporization L_c (T) and saturation vapor pressure p_s of water cannot be neglected. For this, we used one of the most accurate formulations for calculating the saturation vapor pressure of water known as the Tetens formula (Tetens 1930). This formula provides a very good approximation for the saturation vapor pressure of water with an error lower than ~1% (Huang 2018) in the temperature range 0–100 °C and is given as follows: $$p_{\text{s}}(T)=\{\begin{array}{ll}{p}_{\text{s}}^{\text{ref}}\exp \left(\frac{17.625T}{T+238.3}\right)\text{},\hfill & \text{if }T\ge 0\hfill \\ p_{\text{s}}^{\text{ref}}\exp \left(\frac{21.875T}{T+265.5}\right)\text{},\hfill & \text{otherwise}.\hfill \end{array}$$(20)

Here $p_s^{\text{ref}}=610.78\text{Pa}$ is the reference water vapor pressure at the triple point, and T is given in °C and equal toT − 273 for T in units of Kelvin.

Next the term c_pc in Eq. (17) is the temperature-dependent specific heat capacity of water adapted from Wagner & Pruß (2002) and applicable to high-temperatures conditions given by (also see Katyal et al. 2019) $$\frac{c_{\text{pc}}(T)}{R_{\text{c}}}=1+n_3+{\sum }_{i=4}^8{n}_i\frac{{({\gamma }_i\tau )}^2\exp (-{\gamma }_i\tau )}{{[1-\exp (-{\gamma }_i\tau )]}^2},$$(21)

where τ = T_c∕T and T_c is the critical temperature (647 K). The values of the coefficients n_i and γ_i are taken from Table 6.1 of Wagner & Pruß (2002). The total heat capacity of the dry mixture c_pd is calculated as a linear combination of individual heat capacities of the noncondensable species weighted by their volume mixing ratios. The total c_p is thus the sum of specific heat capacity of moist and dry component, weighted by the respective volume mixing ratios, and is given by $$c_{\text{p}}(T)=\frac{{\sum }_i{x}_{\text{di}}{c}_{\text{pdi}}(T)+{\sum }_i{x}_{\text{vi}}{c}_{\text{pci}}(T)}{{\sum }_i{x}_{\text{di}}+{\sum }_i{x}_{\text{ci}}},$$(22)

which is valid for temperatures in the range 200–2000 K. The denominator in this equation is equal to 1 for our case. The specific heat capacity of dry gases such as CO₂, CO, and H₂ is derived from the Shomate equation (Parks & Shomate 1940) in units of kJ kg⁻¹ K, $$c_{\text{p}}=A+BT+C{T}^2+D{T}^3+E/T^2,$$(23)

where the coefficients A, B, C, D, and E are taken from Chase (1998) and T = Temperature∕1000 K for the temperature range valid between 298 and 6000 K. For temperatures below 298 K, the specific heat capacities of CO, CO₂, and H₂ are taken from analytical expressions based on a least-squares fit of data from Lide (2000). The specific heat capacity as a function of temperature in the range 200–3000 K for the four gases we considered is shown in Fig. 3.

In order to obtain the temperature profiles, we followed the same procedure as in Katyal et al. (2019), that is, for surface temperatures T_s > 647 K, we used the dry adiabatic lapse rate (Eq. (19)) by assuming a fixed surface temperature and pressure at the bottom of the atmosphere (Eq. (13)). The temperature profile follows the moist adiabatic lapse rate (second term in Eq. (17)) when the dry adiabatic lapse rate intersects with the moist adiabatic lapse rate. The related work of Katyal et al. (2019) calculated temperature profiles with only one gas species, that is, H₂O, in the atmosphere. We have additionally included CO₂, CO, and H₂ calculated from the volatile speciation model described in Sect. 2.1. Our model accounts for the convective processes in the atmosphere up to an altitude corresponding to the top of the atmosphere pressure p_ToA = 0.1 Pa (similar to Pluriel et al. 2019). The atmospheric layers are vertically spaced in log pressure coordinates. When the temperatures in the dry/moist adiabat fall below 200 K, the temperature structure in our calculations follows an isothermal profile with a constant tropospheric temperature of 200 K.

Fig. 3 Specific heat capacities as a function of temperature for volatile species H2O, CO, and CO2 on the left x-axis and H2 on the right x-axis as obtained from the Shomate equation described in the main text.

2.4 Radiative transfer code

The generic atmospheric radiation line-by-line infrared code (GARLIC; Schreier et al. 2014) was used to compute thermal emission and transmission spectra from10 to 30 000 cm⁻¹ with p, T and the composition of the atmosphere as the input. A detailed description of the radiative transfer code is also available in Katyal et al. (2019). For the verification and validation of GARLIC, see Schreier et al. (2018a,b). The absorption coefficients were calculated line-by-line (lbl) for H₂O, CO₂, H₂, and CO from the HITRAN2016 database (Gordon et al. 2017). Table 2 shows the list of various molecules and sources for the continuum available and relevant for this study.

Rayleigh scattering by molecule H₂O was obtained using the formalism described in Murphy (1977) (also see Scheucher et al. 2020). For other molecules such as CO₂, CO, and H₂, we have used the Rayleigh cross sections using the formalisms of Sneep & Ubachs (2005) and Shardanand & Rao (1977) across the whole spectrum in units of cm² molecules⁻¹ expressed as $${\sigma }_{\text{Rayleigh},\text{i}}(\lambda )={\sigma }_{0,i}{\left(\frac{{\lambda }_{0,i}}{\lambda }\right)}^{\alpha }.$$(24)

Here, α = 4, σ_0,i and λ_0,i are the reference cross sections and wavelengths for the various molecules used in this study that are listed in Table 3.

To mimic optically thick clouds in our simulations (Sect. 5.4), we assumed clouds were present throughout the atmosphere (Kaltenegger & Traub 2009) with the uppermost cloud deck coinciding with the lowermost T = 200 K altitude. The cloud-scattering cross section is expressed in a similar formalism as Eq. (24), but with α = 0. This therefore acts as a scaling term with no wavelength dependence. For reference, the terrestrial aerosol-scattering cross section is about 10⁻²⁷ cm² (Moran et al. 2018). That study (Moran et al. 2018) ruled out haze-scattering cross sections smaller than 9 × 10⁻²⁵ cm² to 3σ and 3 × 10⁻²³ cm² to 1σ. Therefore we assumed a cloud cross section ${\sigma }_{\text{cloud}}\equiv {\sigma }_{0_i}=6\times {10}^{-25}{\text{cm}}^2$, that is, 6 × 10² times larger than the Earth-scattering cross section in our simulations as presented in Sect. 5.4.

The aerosol optical depth integrated along the path is then given by $${\tau }_{\text{aerosol}}={\sigma }_{\text{aerosol}}\cdot N_z,$$(25)

where N_z is the column density of the molecular species in units of molecules cm⁻² (see also Kaltenegger & Traub 2009; Yan et al. 2015).

As a first step, we performed convective-radiative calculations without coupled chemistry similar to previous modeling studies of the magma ocean period. For future work, we plan to include the effect of chemical equilibrium in the lower atmospheric layers by applying the self-consistent convection-climate-photochemistry column model (Scheucher et al. 2020; Wunderlich et al. 2020) to the magma ocean period in order to study the effect upon transmission and thermal emission spectra and atmospheric escape, for example.

2.5 Transmission spectra

The volatile chemical speciation model along with the atmospheric model described in Sects. 2.1–2.4applicable to the early Earth magma ocean phase can be used to characterize exoplanets in the magma ocean phase (see Fig. 12 of Nikolaou et al. 2019), for instance, to detect the spectral features and obtain the vertical extent of the atmosphere observationally. We used a well-established diagnostic tool, transmission spectroscopy, to probe the molecular composition of exoplanetary atmospheres. It provides the effective height of the atmosphere as a function of wavelength (see, e.g., Benneke & Seager 2012).

In GARLIC, the theoretical transmission spectrum is expressed as the effective height of the atmosphere given as follows (Kaltenegger & Traub 2009; Schreier et al. 2018b; Katyal et al. 2019): $$h(\lambda )={\displaystyle {\int }_0^{\infty }A}(z,\lambda )\text{d}z={\sum }_i{A}_i(\lambda )\Delta h_i,$$(26)

where $A_i=1-{\mathcal{T}}_i$ is the absorption along the ith incident ray with a transmission ${\mathcal{T}}_i$ that traverses through the (exo)planetary atmosphere tangentially at a height h_i and continues to travel to the observer at Earth.

The atmospheric transit depth, t_atm is given by (Wunderlich et al. 2020) $$t_{\text{atm}}(\lambda )=\frac{{(R_{\text{p}}+h(\lambda ))}^2}{R_{\text{s}}^2}-\frac{R_{\text{p}}^2}{R_{\text{s}}^2},$$(27)

where R_p is the planetary radius, and R_s is the stellar radius.

2.6 Atmospheric escape

For the atmospheric escape of H₂, we considered hydrodynamic fluxes of escaping H₂, possibly diffusion-limited at the homopause ~100 km, in the presence of a static background atmosphere, similar to the approach of Kuramoto et al. (2013), Zahnle & Catling (2017) and Zahnle et al. (2019). At higher levels above the homopause, where the escaping H₂ gas is accelerated as a result of the absorption of stellar high-energy photons, hydrogen exists mainly in atomic form. For the quantitativeimplementation of the escape fluxes, we relied on Zahnle et al. (2019), who solved the hydrodynamic, possibly diffusion-limited, problem in the atmosphere. In particular, we used their Eq. (3), which parameterizes the transition from energy-limited escape to diffusion-limited escape. Diffusion-limited escape becomes effective when the atmosphere is more strongly irradiated and the escaping hydrogen cannot be replenished from below the homopause.

The diffusion-limited flux ϕ_dl generally provides an upper limit to the escape of lighter species such as hydrogen in a background heavier and stationary gas. The escape rate for H₂ can be written as (Hunten 1973) $${\varphi }_{\text{dl}}=b_{aj}{f}_{\text{tot}}({\text{H}}_2)\left(\frac{1}{H_a}-\frac{1}{H_{{\text{H}}_2}}\right),$$(28)

where b_aj is the binary diffusion coefficient between the background heavier gas a such as CO₂, CO and H₂O and the escaping gas H₂. Here, f_tot(H₂) is the hydrogen volume mixing ratio (VMR) at the homopause. H_a and $H_{{\text{H}}_2}$ represent the unperturbed scale heights (= kT∕mg) of the heavier gas a and the escaping gas H₂ at the homopause, respectively. We consider H₂ as the dominant H-bearing species, that is, f_tot ≡ f(H₂). The binary diffusion coefficients b_aj for H₂ in CO and CO₂ are roughly similar and taken to be 3 × 10²¹ m⁻²  s⁻¹ (Marrero & Mason 1972). We validated the escape rate obtained for the current Earth using Eq. (28). The total hydrogen VMR was taken to be f_tot(H₂) = 7.15 ppm with the background-dominant gas as N₂. The diffusion coefficient between the H₂ and N₂, b_aj of 1.7 × 10²¹ m⁻²  s⁻¹ was taken from Hunten (1973). We verified that the diffusion-limited escape of H₂ is ~ 1.5 × 10¹²  H₂  molecules m⁻²  s⁻¹, as also reported by Hunten & Donahue (1976). In terms of the loss of H-atoms to space, this is equivalent to a mass-loss rate of ~3 kg s⁻¹.

The energy-limited flux was obtained as the ratio of stellar XUV energy incident upon the planet and the energy required to lift a given mass out of the Earth potential well and into the space. The mass-loss rate of hydrogen using the energy-limited formula (Watson et al. 1981) for hydrodynamic escape (also see Koskinen et al. 2014; Hamano et al. 2015) is given by $${\dot{M}}_{\text{el}}=\frac{\pi \varepsilon R_{\text{p}}{R}_{\text{XUV}}^2{F}_{\text{XUV}}(t)}{G{M}_{\text{p}}}\text{g}{\text{s}}^{-1},$$(29)

where R_XUV is the radial distance at which the XUV energy from the star is deposited, and satisfies R_XUV > R_p, and ɛ is the heating efficiency <1. Upon assuming R_XUV ≡ R_p (Zhang 2020), this equation can be written as $${\dot{M}}_{\text{el}}=\frac{\pi \varepsilon R_{\text{p}}^3{F}_{\text{XUV}}(t)}{G{M}_p}\text{g}{\text{s}}^{-1}.$$(30)

Here, R_p is the Earth radius (6.4 × 10⁸ cm), M is the mass of Earth (5.4 × 10²⁷ gm), G is the gravitational constant (6.67 × 10⁻⁸ cm³ g⁻¹ s⁻²), F_XUV(t) is the time-evolving flux from the host star at XUV wavelengths at 1 AU obtained by the relation $F_{\text{XUV}}(t)=5{(4.5/t)}^{1.24}$ in units of erg cm⁻² s⁻¹, where t is the age of the Sun in billion years (Zahnle et al. 2019), and ɛ is taken to be 0.5 to be consistent with Zahnle et al. (2019). Owing to the fast temporal evolution in the XUV output of the young Sun (related to whether it was a fast, moderate, or slow rotator), the S value defined as S(t) ≡ F_XUV∕F_XUV⊙ can range between 10 and 100 during the Hadean (Tu et al. 2015; Lammer et al. 2018). The lifetime of H₂ atmospheres on early Earth is calculated as the ratio of mass-loss rate (Eq. (30)) to the total mass of H₂ in grams, $m_{{\text{H}}_2}$, in the atmosphere as $$t~\left[\frac{G}{\varepsilon \pi }\right]\left(\frac{M_p}{R_p^3}\right)\frac{m_{{\text{H}}_2}}{F_{\text{XUV}}}.$$(31)

Zahnle et al. (2019) have shown that the diffusion-limited escape applies to conditions of higher levels of irradiation, whereas energy-limited escape may be a better approximation for lower levels of irradiation. Building upon their findings, we assumed that the flux is given by the expression $${\varphi }_{{\text{H}}_2}=\frac{1\times {10}^{12}{f}_{{\text{H}}_2}S}{\sqrt{1+0.006S^2}}\text{molecules}{\text{cm}}^{-2}{\text{s}}^{-1},$$(32)

which is a good fit to their Fig. 5 if the VMR of H₂ is lower than 0.2.

The hydrogen mass-loss dM (in grams) due to escape at a particular time step dt is approximated as follows: $$\text{d}M=4\pi R_{\text{p}}^2{\varphi }_{{\text{H}}_2}{m}_{{\text{H}}_2}\text{d}t.$$(33)

Here, $m_{{\text{H}}_2}$ is the mass of H₂ in grams per molecule, and ${\varphi }_{{\text{H}}_2}$ is the escaping flux from Eq. (28) or Eq. (32).

The actual radius contribution in Eqs. (30) and (33) comes from the sum of the radius of the planet and height of the exobase, that is, (R_p +H_e). The height of the exobase, H_e, is given by the altitude level where the Lyman-edge opacity of hydrogen becomes 1. For the scenario presented in Sect. 4.6, we obtain that (R_p +H_e) ~ R_p (calculation not shown here).

3 Volatile outgassing scenarios

We used three main scenarios in this work. Scenario 1 is for a surface temperature of 3300 K and surface pressure of 76.7 bar, that is, at the beginning of the magma ocean (BOM) period, as shown in Fig. 2 (top) and Table 4. Scenario 2 is similar to scenario 1, but for a surface temperature of 1650 K and pressure of 395 bar, that is, at the end of the magma ocean (EOM) period, as shown in Fig. 2 and Table 5. The conditions where the MO phase ends are T_s = 1650 K at a solidification timescale of ~1 Myr according tothe coupled interior-atmospheric model results of Katyal et al. (2019) and Nikolaou et al. (2019). Previously, a coupled atmospheric-interior model by Salvador et al. (2017), which was applied around the MO phase, delivered P-T conditions with a similar H₂O-CO₂ composition as this study.

Scenario 3 is for representative cases taken from the study of Pluriel et al. (2019), as shown in Table 6 referring to an H₂O dominated and a CO₂ dominated atmosphere. Scenarios 1 and 2 cover a representative range of P, T, and mantle oxidation state during the MO period. The buffers in the three scenarios were chosen to represent the mantle redox state, ranging from reducing to oxidizing.

Our choice of three scenarios (with varying initial volatile abundances of H₂O and CO₂) and varying mantle fugacity, hence leading to different atmospheric compositions, is motivated by previous works in the literature that also explored the effect of composition by varying key species such as H₂, CO, CO₂, H₂O for magma ocean studies (Katyal et al. 2019; Nikolaou et al. 2019; Bower et al. 2019), Venus-type planets with CO₂ dominated atmosphere (Hamano et al. 2013; Pluriel et al. 2019), pure steam-dominated atmosphere (Katyal et al. 2019; Schaefer et al. 2016; Hamano et al. 2015), and the effect of adding H₂ to the atmosphere of early Mars (Ramirez et al. 2014).

4 Results

4.1 Temperature profiles

Figures 4 and 5 show the pressure-temperature (p–T) and altitude-temperature (z–T) profiles obtained using the convection lapse-rate formulation as described in Sect. 2.3 for scenarios 1 and 2 as shown in Tables 4 and 5, respectively. In Fig. 4 scenario 1.1 (95% initial CO₂, 5% initial H₂O, top left panel, p–T), the results suggest a dry adiabat in the lower unsaturated troposphere that is related to the high surface temperature. The green curve (oxidized buffer case) adiabat is steeper in this region than the blue and orange curves because it is CO₂-dominated, leading to alow heat capacity (Fig. 3) and accordingly to a steeper calculated lapse rate (Eq. (19)). This steep slope results in the green curve that intersects the saturated vapor curve (dashed red line) at lower pressures than for the blue and orange curves, which have lower CO₂ but more CO and H₂O (Table 5) owing to higher heat capacities (Fig. 3) and therefore shallower dry lapse rates (Eq. (19)).

In Fig. 4 for z–T (top right panel), the green curve features the lowest atmospheric height (geometric thickness) due to its larger atmospheric molecular weight and therefore smaller scale height than the other, lighter atmospheres shown by the orange and blue curves. Above the dry adiabat regime (e.g., above about 320 km for the orange line), the slope steepens corresponding to the wet adiabat as temperatures are low enough to allow for condensation. In the uppermost atmosphere, a fixed iso-profile temperature (T = 200 K) (e.g., occurring above about 380 km for the orange line) is imposed in our model for the radiative regime, thereby following other studies (Lupu et al. 2014; Pluriel et al. 2019; Katyal et al. 2019).

In the middle and lower panels of Fig. 4, the initial H₂O (CO₂) amounts are increased (decreased) (see scenarios 1.2 and 1.3 in Table 4), which leads to light, hydrogen-dominated atmospheres for the reduced buffer cases. Because H₂ has a heat capacity between that of H₂O and CO₂ or CO (Fig. 3), the green slope in the middle and lower left panels is now closer to the blue and orange slopes compared with the upper panel because the difference in heat capacities is now smaller. In the middle and lower right panels, results suggest extended scale heights for highly reduced atmospheres with dominant H₂ (e.g., scenario 1.3, IW-4) because its molecular weight is lower than that of oxidized species such as H₂O or CO₂ in the atmospheres (e.g., scenario 1.1 and 1.2, IW+4) and subsequently the change in altitude range of the figure is notable.

Figure 5 illustrates the EOM scenarios (Table 5). Compared with the earlier BOM case (Fig. 4) with a surface temperature of 3300 K, the surface temperature for the EOM scenario has cooled to 1650 K and the surface pressure has increased from 76.7 bar (BOM) to 395 bar (EOM). The cooler temperatures generally favor the formation of chemical species with weaker bond dissociation energies, that is, H₂ (4.52 eV) over H₂O (5.1 eV) andCO₂ (5.51 eV) over CO (11.1 eV) (where eV denotes electron Volt), which is consistent with the changes produced by the chemical speciation model (compare Table 5 with Table 4). Larger changes of up to a few dozen percent in the species mole fraction occur for scenario 2.2 compared with scenario 1.2, in which the relative amount of C-H-O (0.75:0.25 for H₂O:CO₂) is the most similar. Figure 5 (EOM) shows modest changes in comparison to Fig. 4 (BOM), which are associated with the change in p, T and corresponding changes in speciation and heat capacity, hence in the slope of the adiabats. The altitude decrease is mainly associated with the decrease in temperature.

Figures 4 and 5 show that (1) the transition between dry and moist adiabat (slope change) always appears to be located on the saturation curve P_sat and it effectively occurs at lower temperature than the saturation temperature (for similar pressure); and (2) for $f_{{\text{H}}_2\text{O}}$ ≥ 0.1 in the atmosphere, the moist adiabat always coincides with the water vapor saturation curve P_sat.

Figure 6 shows the temperature-pressure profiles obtained for scenarios 3.1 (upper panel) and 3.2 (lower panel) with H₂O and CO₂ dominated atmospheres of Table 6, respectively. In the case of scenario 3.1 with its initially H₂O-dominated atmosphere (Fig. 6 upper panel), for the reducing buffer IW case, the gradient of the dry adiabatic lapse rate is slightly reduced (orange line) due to chemical speciation of H₂ and CO as compared to the “no speciation” case (blue line). However, for an oxidizing buffer IW+4 (green line), the dry adiabatic lapse rate is similar to the case for “no speciation” because the volume mixing ratio of outgassed species is similar to the initial volatile outgassing. In Fig. 6, scenario 3.2 (lower panel), the chemical speciation of the initially assumed CO₂-dominated atmosphere (see Table 6), leads to a CO-rich (IW buffer) and CO₂-rich (IW+4 buffer) atmosphere. Thus, the moist adiabatic lapse rate changes as compared to the upper panel results and shows a deviation from pure water vapor saturation curve (dotted red line). This is because of a low water vapor volume mixing ratio $f_{{\text{H}}_2\text{O}}=0.02$ and a higher CO₂ volume mixing ratio f_CO2 = 0.98, leading to an unsaturated troposphere. Moreover, the vertical range of the wet adiabat regime (scenario 3.2) is smaller than in scenario 3.1 (upper panel), and the radiative regime (T = 200 K isothermal profile) is evident at 40 Pa, for instance, for the green curve (lower panel) but occurs near the top of the atmosphere, that is, around 0.1 Pa, for scenario 3.1. The differing behavior in scenario 3.2 compared with scenario 3.1 arises mainly from the differentchemical speciation (see panel legends), which affects the heat capacity and therefore the adiabatic gradients.

Fig. 4 Left: pressure-temperature profiles for the BOM scenarios 1.1 (upper), 1.2 (middle), and 1.3 (lower panel) for three different buffers IW-4, IW, and IW+4. The dashed red line represents the saturation water vapor curve. Pressures at the bottom of atmosphere values are taken from Table 4. Right: altitude-temperature profiles for scenarios 1.1 (upper), 1.2 (middle), and 1.3 (lower panel). We note the different altitude ranges for the panels on the right.

Fig. 5 Same as Fig. 4, but for p, T conditions and scenarios as summarized in Table 5, representing the EOM.

4.2 Thermal spectral emission

Using the temperature profiles as obtained in Figs. 4, 5 and composition of various species for various scenarios of Tables 4 and 5 as input to the line-by-line (lbl) radiative transfer code GARLIC, we obtain the thermal spectral emission from the atmosphere, also known as the outgoing longwave radiation (OLR). Figure 7 shows the thermal emission spectra of scenario 1.1(upper), 1.2 (middle), and 1.3 (lower) panel of Table 4 for the BOM. The spectral features that arise due to the presence of reduced species such as CO versus the oxidized species such as H₂O and CO₂ are easily distinguishable in the three panels. In the top panel of Fig. 7 (scenario 1.1), a prominent CO absorption feature is seen at 2.3 and 4.6 μm for the reduced buffer case IW-4 (blue curve), whereas distinct CO₂ features at 2.1, 4.3 and 15 μm are seen for the oxidized buffer case IW+4 (green curve). When the input H₂O outgassing value is increased (scenario 1.2), stronger H₂O features start to appear (orange and blue curves) along with the CO₂ absorption features for the oxidized buffer (green curve). Finally, when the input H₂O outgassing is increased to 100% (scenario 1.3), we obtain thermal spectra that are dominated by water vapor with distinct water features between 1 and 2.6 μm and at 6.2 μm for the oxidized case (green curve). The reduced atmosphere with IW-4 buffer produces an almost pure H₂ atmosphere with a similar spectrum as H₂O but with a higher emitted flux.

Figure 8 shows the thermal emission spectra of scenario 2.1 (upper), scenario 2.2 (middle), and scenario 2.3 (lower) of Table 5 for the EOM. One notable difference between Figs. 7 and 8 is a significant reduction in the emitted flux especially at the smaller wavelengths, which is mainly related to the lower surface temperature (i.e., 1650 K) for the latter. The effect of pressure at the bottom of the atmosphere on the thermal spectra is also seen and is discussed in detail in Sect. 4.3. The effect of the large H₂ continuum at 5–7 μm caused by self-collision of the H₂ molecules could reveal molecular H₂ as a main constituent of an atmosphere not known a priori.

Figure 9 compares the thermal emission spectra for reduced atmospheres overlying a reducing buffer (IW-4) for the various scenarios 2.1, 2.2, and 2.3 of Table 5. In this figure, the prominent absorption features are seen for CO at 2.3 and 4.6 μm, while CO₂ stands out at 2.1, 4.3, and15 μm. A small amount of CO₂ is always present even for this very reducing buffer calculated by the mass-equilibrium method (Sect. 2.1). However, the strength and width of CO₂ absorption features are seen to decrease as the volume fraction of initial CO₂ decreases from 0.95 to 0.0 (see the zoomed 15 μm feature in the right inset of Fig. 9). Another notable feature is visible in scenario 2.3 with pure H₂O (dissolved in the melt), which results in a H₂ -dominated atmosphere and displays a large continuum absorption in the 5–7 μm region caused by collision-induced absorption.

Figure 10 is same as Fig. 9, but now for the oxidizing atmosphere cases (IW+4) for scenarios 2.1, 2.2, and2.3 of Table 5. A reduction in the emitted flux is seen compared with the reduced atmosphere case (Fig. 9). This is because of enhanced absorption by higher concentrations of greenhouse gases such as CO₂ and H₂O that are present in the oxidized atmosphere. Thus, as the initial H₂O volume mixing ratio (VMR) increases (from scenario 2.1 through 2.3), a further reduction in the emitted flux is obtained. The most prominent absorption features are seen for CO₂ at 2.1, 4.3, and 15 μm and H₂O at 1.0, 1.2, 1.4, 1.9, and 2.6 μm, respectively.As a result of the thick overlying steam atmosphere, the atmospheric window (8–10 μm) becomes optically thick and the 6.2 μm H₂O feature is not prominent (green curves in scenarios 2.2 and 2.3). The inset of this figure shows the zoomed-in 4.3 μm CO₂ feature. The CO feature at 4.6 μm is not seen here as compared to the reduced case in Fig. 9. The inset on the far right shows a zoomed-in view of the 15 μm CO₂ absorption band as a function of CO₂ VMR and indicates the change in band depths when the input CO₂ VMR is reduced from 0.95 to 0.25 (red to blue curve) and is then absent for the green curve with no CO₂. In summary, the main result is a diminishing of the CO band at 4.6 μm and the appearance of distinct H₂O bands from ~1–2.6 μm. The H₂O feature at 6.2 μm is less prominent because of the thick steam atmosphere.

Fig. 6 Temperature profiles p − T for scenario 3.1 (upper panel) with an initially H2O-dominated atmosphere and scenario 3.2 (lower panel) with an initially CO2-dominated atmosphere (Table 6). The molar fraction of each species is shown in the legend. “No speciation” indicates that the values were input directly by the user and were not calculated by the C-O-H speciation model.

4.3 Dependence of OLR upon mantle fugacity

The thermal emission spectra were then averaged over the entire wavelength range to obtain the outgoing longwave radiation (OLR) as shown in Fig. 11 (right panels). The variation in pressure at the bottom of atmosphere p_boA (using Eq. (13)) for three different buffers is shown in Fig. 11 (left panels). The details of the various cases for which we investigated the effect of speciation on the p_boA and subsequently on the OLR for the BOM and EOM are shown in Tables B.1 and B.2 respectively.

For the BOM (Fig. 11; top right panel) with T_s = 3300 K, the highest OLR is obtained for the most reduced atmosphere (blue curve; IW-4 buffer) because the absorption bands are weaker and the atmospheric pressures (top left panel) are lower than the IW buffer (orange curve), wherein a mix of reduced and oxidized species are present in the atmosphere (see Table B.1) that contribute far more effectively to block radiation. It is interesting to note that even for the most oxidizing atmosphere (green curve; IW+4 buffer), less radiation is blocked by the atmosphere, which is mainly dominated by H₂O, than for the IW buffer.

At lower temperatures (Fig. 11; bottom right panel) with T_s = 1650 K (EOM), theemitted radiation is far lower than in the high temperature case (T_s = 3300 K) because for the latter, the radiation is emitted in the infrared and visible wavelengths, which leads to a stronger outgoing longwave radiation (OLR). Moreover, for the low temperature case (T_s = 1650 K), the OLR reaches a limit of ~280 W m⁻² (Nakajima et al. 1992; Goldblatt et al. 2013; Marcq et al. 2017) for water-dominated atmospheres (green curve; IW+4 buffer) at ~1800 K (also see Katyal et al. 2019). The calculated OLR does not drop below the OLR limit because of the impeding effect of absorption of radiation by H₂O even if present in moderate quantities. This is an important result that suggests that for a specific surface temperature, atmospheres with lower surface pressures (~4–200 bar) can lead to strong outgoing emission of radiation and therefore to a more effective cooling of the MO than for atmospheres with a high surface pressure (>200 bar). Studying the effect of pressure on OLR in this manner is analogous to the effect of a threshold temperature as shown in Marcq et al. (2017). We also refer to Fig. 10 of Nikolaou et al. (2019), where the OLR for pure steam atmospheres is plotted as a function of surface temperature and surface pressures. According to their Fig. 10 and Katyal et al. (2019), higher OLR values are obtained for atmospheres with lower surface pressure (for a similar surface temperature). This means that the outgoing radiation for lower pressure cases may reach the OLR limit, but at temperatures lower than ~ 1800 K as previously reported (Kopparapu et al. 2013; Katyal et al. 2019).

In Fig. 11, cases 1 and 2 are mainly dominated by CO₂ and/or CO (see Tables B.1 and B.2), which show only few absorption bands and thus display the strongest OLR. Cases 3 and 4, which are composed of a mixture H₂O+CO₂+CO+H₂ (see Tables B.1 and B.2), result in a minimum OLR where the spectra shows absorption lines of almost all the species present. Cases 5 and 6 are dominated by absorption through H₂O and/or H₂, which leads to a stronger OLR than the mixed-species cases 3 and 4, where absorption of all species strongly reduces the OLR for the high-temperature cases. For the low-temperature cases (lower right panel), the effect of the absorption by H₂O dominates because the H₂O band and continuum absorption cover the entire infrared wavelength regime.

Fig. 7 Comparison of thermal emission spectra for a reduced vs. an oxidized mantle during the BOM presented for scenario 1.1 (upper panel), 1.2 (middle panel), and 1.3 (lower panel) of Table 4. The output of the volatile speciation model is shown in the legend. Values written above the legend show the initial input values to the speciation model. Reference blackbody curves for four different temperatures are plotted. Absorption features ofkey chemical species are indicated. The spectra shown here have been binned to a resolution of λ∕Δλ = 1000.

Fig. 8 Same as Fig. 7, but for scenarios of Table 5. The blackbody curves for three different surface temperatures bracketing the emission spectra are plotted for reference.

4.4 Early H₂ atmosphere scenario

While hydrogen is prone to escape to space via strong incoming extreme UV (EUV) during the Hadean (more details in Sect. 4.6), it is also an important gas in enhancing the greenhouse effect in the atmosphere. The efficiency of heating the atmosphere, however, depends upon the surface partial pressure of hydrogen.

We explored a case similar to scenario 2.3 of Table 4 and shown previously in Fig. 9 with a pure H₂ atmosphere, but with a varying hydrogen surface pressure $p_{{\text{H}}_2}$ and a fixed surface temperature T_s = 2800 K. We obtained the thermal emission of radiation by varying the surface pressure $p_{{\text{H}}_2}$ in the range 0.6–185 bar, as shown in Fig. C.1 (upper panel). We observe a large number of undulations around 2 μm in the obtained thermal spectra, especially at higher surface pressures. These undulations are caused by the effects of collision-induced absorption (CIA) by H₂. The lower panel of Fig. C.1 shows a decrease in the OLR with the (surface) partial pressure of H₂ attributed to the higher absorption of radiation at higher pressures, as also shown in Fig. 11 (left panels). We therefore conclude that several bars of H₂ in the atmosphere could significantly affect the climate in the cases studied, leading to additional greenhouse warming (also see Ramirez et al. 2014). A similar trend in the OLR versus $p_{{\text{H}}_2}$ as shown here was suggested by Pierrehumbert & Gaidos (2011).

Fig. 9 Thermal spectral emission for a reducing atmosphere overlying a reduced mantle (IW-4) for scenarios 2.1, 2.2, and 2.3 of Table 5. The legend shows the mole fraction of the outgassed species obtained from the speciation model. Prominent features of absorbing species are seen at wavelengths marked in the figure. Reference blackbody curves for three different surface temperatures are plotted. The insets show the zoomed-in thermal spectra with the most prominent absorption features of CO2, CO and H2. The spectra have been binned to a resolution of λ∕Δλ = 1000.

Fig. 10 Same as for Fig. 9, but for an oxidizing atmosphere overlying an oxidized mantle (IW+4). Notable features are the emergence of prominent H2O features between 1 and 2 μm (scenario 2.3) and the disappearance of CO feature at 4.6 μm as compared to Fig. 9 (scenario 2.1 and 2.2).

4.5 Transmission spectra

We obtained the transmission spectra for the BOM and the EOM using the lbl radiative transfer code GARLIC and the temperature profiles (Sect. 4.1) and concentration as an input. The resultant spectra are shown as the wavelength-dependent effective height of the atmosphere (see Eq. (26) in Fig. 12). The spectra were binned to a resolution of R = λ∕Δλ = 1000 over the wavelength range shown.

As shown in Fig. 12, the effective height of the atmosphere increases to one order of magnitude for the reduced atmospheres (IW-4) from scenarios 1.1 to 1.3 (and 2.1 to 2.3) for both BOM (left panel) and EOM (right panel) scenarios. This is related to the increasing atmospheric scale height for the lighter H₂ –dominated atmospheres. Similarly, for the oxidized atmospheres (IW+4), the effective height increases from a dense and thick CO₂-dominated atmosphere (scenarios 1.1 and 2.1) to lighter H₂O-dominated atmospheres(scenarios 1.2, 1.3, 2.2, and 2.3). In the visible-optical range, the slope of the so-called Rayleigh-scattering extinction feature is somewhat steeper for the reducing atmosphere (blue line) because of the different scale heights. For scenarios 1.3 and 2.3, H₂O absorption features are seen ubiquitouslyfrom roughly 0.9 to 2.6 μm. For a high CO₂ initial volatile outgassing (scenario 1.1 and 2.1), CO absorption features are seen at 2.3, 4.6 μm for IW-4 buffer, and CO₂ absorption features are seen as shown for2.1, 4.3, and 15 μm for the IW+4 buffer. For the EOM scenarios (right panels), the cooler surface temperature (1650 K) leads to a decrease in scale height (H = kT∕mg) compared with the BOM (left panels) scenarios with higher surface temperature (3300 K).

A decrease in the depth of the absorption bands for CO (IW-4) and CO₂ (IW+4) is seen at 4.3 and 15 μm, respectively, which is related to a decrease in the input CO₂ from scenarios 1.1 to 1.2 at T_s = 1650 K (right panel). A similar trend is seenfor scenarios at T_s = 1650 K (right panel). For the oxidizing atmospheres (IW+4) in scenarios 1.3 and 2.3, the most prominent H₂O features between 1− 2 μm are seen, but the 6.2 μm H₂O feature is suppressed. This is likely related to the overlying optically thick atmosphere of H₂O (blanketing effect) with a very high surface pressure (~395bar), wherein the radiation will be dominated by the emission from the top of the atmosphere. This hinders chemical species in the lower atmosphere from being detected in the upper regions sampled by transmission spectroscopy. It is therefore difficult to probe dense atmospheres with high optical depths.

Fig. 11 Left: pressure at the bottom of the atmosphere pboA for various cases referring to different initial volatile abundances and three different mantle buffers: reduced (IW-4), IW, and oxidized (IW+4). The details of speciation for various cases on the x-axis are shown in Table B.1 (Ts = 3300 K) and Table B.2 (Ts = 1650 K) in the Appendix B. The horizontal dashed black line indicates the total surface pressure of the initial outgassed volatile for the various cases. Right: outgoing longwave radiation (OLR) for various cases corresponding to the atmospheric surface pressure shown in the left panels and for three different buffers.

4.6 Atmospheric evolution and escape over the MO period

Table 7 shows the range of surface partial pressure of H₂ for a highly reduced (IW-4) and an oxidized mantle (IW+4) over a period of MO solidification ~ 1  Myr. For the reducing-buffer case (IW-4), over the period of MO, the surface pressure of atmospheric hydrogen increases from 2.3 to 102 bar, while for the oxidized case, it increases from 0.02 to 2.7 bar. To investigate the evolution of H₂ along with other species in the atmosphere, we split the MO period into ten time steps where the change in species abundance from scenario 1.1 (BOM) to scenario 2.2 (EOM) in Table 7 is set to change with a linear time step of 0.1 Myr (only shown for H₂ in the table). The resulting changes in species abundances across the MO can be attributed to outgassing and chemical speciation. As a first step, the input data for the various scenarios were taken from Nikolaou et al. (2019) and Katyal et al. (2019), but in a future study, we plan to perform time-dependent, coupled interior-atmosphere evolutionary scenarios over the magma ocean phase applying a newly developed climate-chemistry atmospheric model (Scheucher et al. 2020; Wunderlich et al. 2020).

We now consider the evolution of atmospheric species over the MO period. The mole fraction f of the resultant outgassed species i is translated into the mass fraction via $w_i=f_i{M}_i/\overline{M}$, where M_i is the molecular weight of the species i and $\overline{M}$ is the mean molar mass of the atmosphere. The atmospheric mass of the individual species is obtained by multiplying the calculated mass fraction w_i with the total mass of the atmosphere M = P*A∕g.

Figure 13 shows the resulting atmospheric mass evolution of a reducing atmosphere lying above a reduced mantle, that is, IW-4 (upper panel) and an oxidizing atmosphere lying above an oxidized mantle, that is, IW+4 (lower panel) during the MO period. The circles show the mass of the H₂ escaping via diffusion-limited escape according to Eq. (33) as described in Sect. 2.6, where the color of the circles indicates the total mass-loss rate (g s⁻¹). The mass-loss rate obtained from energy-limited escape formalism as in Eq. (30) is obtained to be 5.2 × 10⁷ g s⁻¹ for S = 10 and 5.2 × 10⁸ g s⁻¹ for S = 100. As described in Sect. 2.6, (R_p + H_e) ~ R_p for this scenario (see Table 7). The mass loss of H₂ (in grams) is then calculated and shown as the yellow shaded region within these S ranges in Fig. 13 (both upper and lower panels). In general, the atmospheric H₂ mass loss (in grams) is seen to increase for both the escape processes because the H₂ abundance increases throughout the MO evolution timescale. The residual H₂ in the atmosphere after accounting for escape is seen to overlie the total H₂ content of the atmosphere, implying that a significant amount of H₂ remains in the atmosphere. We also show a shaded region obtained from the parameterized Eq. (32) of Zahnle et al. (2019), where the lower and upper boundaries of the shaded region indicate S = 10 and S = 100, respectively, and $f_{{\text{H}}_2}$ values are taken over the evolution timescale (see Table 7).

For the oxidized-mantle (IW+4) case (Fig. 13 lower panel), the H₂ abundance is very low (see Table 7), causing this to be the bottleneck for H₂ to diffuse through heavy CO₂ or H₂O up to the homopause. Moreover, the plentiful amount of XUV energy available here (with a high loss rate) can in fact lead to the complete removal of hydrogen. For the oxidizing atmosphere lying above an oxidized mantle, that is, IW+4, escape of H₂ is therefore not energy limited but diffusion limited. Finally, we conclude that for both the redox states of the mantle, the outgassing of H₂ into the atmosphere dominates the escape of H₂ (Kuramoto et al. 2013; Ramirez et al. 2014).

As an illustration, Fig. 14 compares the escape rates for scenario 2.2 (EOM) calculated for H₂ with our estimated H₂ outgassing rate obtained using Eq. (14) as described in Sect. 2.1.

The outgassing rates obtained by us during the magma ocean period in the Hadean (3.8–4.0 Gyr before present) for a variable redox state of the mantle are higher than in the Archean (2–4.0 Ga), that is, 9 × 10¹⁴ molecules m⁻² s⁻¹ (Tian et al. 2005; Kuramoto et al. 2013). First, one of our central results is that outgassing (red solid line) proceeds faster than diffusion-limited escape (closed red circles) by a factor of ×10 and energy-limited escape depending upon the uncertainty in the XUV radiation as evident from the S values displayed in the figure. Second, for an atmosphere with a high hydrogen abundance, that is, VMR(H₂) > 0.2(−4 < IW ≤ 0.5), the energy-limited escape is more effective, whereas for an atmosphere with lower hydrogen abundances, that is, VMR(H₂) < 0.2 (0.5 < IW < 4), the diffusion-limited escape is more effective, as shown in Fig. 14.

Fig. 12 Transmission spectra calculated by GARLIC shown as atmospheric effective height vs. wavelength for the BOM scenarios (left panels) and EOM scenarios (right panels) for three different redox buffers IW-4, IW, and IW+4.The mole fraction of final outgassed species in the atmosphere for various scenarios is the same as provided inTables 4 and 5. The effective height of a reduced atmosphere (with low molecular weight) for the case of reduced mantle IW-4 (blue curve; scenarios 1.2, 1.3, 2.2, and 2.3) is lower than an atmosphere with larger molecular weight (blue curve; scenarios 1.1 and 2.1) and oxidizing atmospheres (green curve; all scenarios). For more details, see the text. The spectra have been binned to a resolution of λ∕Δλ = 1000.

5 Discussion

We have studied the effect of the redox state of the mantle and volatile outgassing of major outgassed species H₂O and CO₂ from the interior and their interaction with the melt resulting in outgassing of newer species at the surface. The effect on the atmospheric profiles, thermal infrared emission, transmission, and escape of H₂ is discussed below.

5.1 Mantle redox and atmospheric temperature

Our study has suggested an interesting coupling between interior outgassing, volatile speciation due to the mantle redox state, atmospheric composition, and atmospheric temperature that affects atmospheric extent, mixing, and the location of the dry and moist adiabat. During the BOM, more oxidized atmospheres switch from dry to wet adiabat at lower pressures (Fig. 4, green lines) than reduced atmospheres due to differences in their mean molecular weight and heat capacities, which leads to a lower scale height. The pressure at which the moist adiabat is reached affects the onset of saturation, thereby affecting the microphysical processes such as condensation and ultimately ocean formation. Detailed microphysical processes are not included in our study and need further investigation with more detailed models. On the other hand, the more reduced atmospheres have larger scale heights, and the switch to the moist adiabat lies at higher altitudes (Fig. 4, blue curve). Precipitation initiated at high altitudes has to sediment over a greater distance to rain out, which could favor evaporation in the lower, warmer atmospheric layers, and which might therefore be less efficient at condensing to form the oceans. For the EOM case (Fig. 5), although the atmosphere cools (from 3300 to 1650 K at the surface) and the surface pressure increases by a factor of 5, the temperature structures for the EOM are similar to those for the BOM (Fig. 4). This is because the pressure increase (which increases atmospheric thickness) is offset by the cooling effect (which decreases atmospheric thickness). In summary, the above processes likely affected the evolution of climate on the early Earth by affecting key processes such as cooling to space and greenhouse heating, but also the onset of condensation and accordingly, the formation of terrestrial ocean.

5.2 Thermal emission and transmission spectra

Our results suggest that the thermal emission of radiation is affected by key parameters such as surface temperature, surface pressure, mantle redox state, and initial volatile content of the two main species H₂O and CO₂ that are outgassed during the magma ocean phase. The most prominent absorption features of oxidized versus reduced species in the atmosphere for the three different mantle oxidation states IW-4, IW, and IW+4 are shown in Figs. 7 and 8. For the reducing mantle cases, the volatile speciation model would result in a mixture of both oxidized and reduced species in the atmosphere depending upon the initial volatile contents $f_{{\text{H}}_2\text{O}}^{\text{init}}$ and $f_{{\text{CO}}_2}^{\text{init}}$ (see Tables 4, 5, B.1, and B.2). As shown in Fig. 11 (right panels), the OLR for the most reduced atmosphere case overlying a reduced mantle (IW-4) consists of H₂, CO or CO+H₂ (depending upon the case) and results in the highest OLR because these molecules have only a few absorption bands. On the other hand, we obtain lower OLR for other moderately reduced buffer, that is, IW, resulting in a mixed atmosphere of CO₂+CO+H₂+H₂O with the largest number of absorption bands covering the entire spectral range. Finally, the lowest OLR is obtained for an oxidizing atmosphere overlying an oxidizing buffer (IW+4), resulting in CO₂, H₂O or H₂O+CO₂ atmospheres(depending upon the case).

As a result, a planet with a reduced mantle buffer (IW-4) and low water content (cases 1 and 2) is expected to cool faster than a planet with a high water content in the atmosphere (cases 4, 5, and 6). For a moderately reduced mantle buffer (IW), the planet will cool down slowly as the OLR is the lowest for it. For a planet with an oxidized mantle buffer (IW+4), the cooling timescales will be slower than those of the IW-4, but faster than the IW case.

The initial outgassed water content along with the surface temperature and the redox state of the mantle are important factors driving planetary cooling. Figure 11 shows for lower temperature (T_s = 1650 K) and oxidized mantle (IW+4 buffer) that the OLR limit ~280 W m⁻² (Goldblatt et al. 2013; Marcq et al. 2017)occurs only for H₂O -rich atmospheres (cases 5 and 6). However, this is not the case for a high temperature (T_s = 3300 K) and an oxidized mantle (IW+4) because at such high temperatures, the emission occurs in visible wavelengths along with infrared, and H₂O is not able to block the radiation in the visible as effectively as in infrared. Radiation is thus not able to reach the OLR limit.

Oxidized atmospheres containing heavier species such as CO₂ result in a dense atmosphere with small effective heights for transmission spectra (Fig. 12, scenarios 1.1 and 2.1) as compared to a reduced and lighter atmosphere with H₂, resulting in extended atmospheres with a large effective scale height, as shown in scenarios 1.3 and 2.3 of Fig. 12. The slope of the Rayleigh-scattering feature in the ultraviolet and visible wavelength range could provide information on the bulk atmospheric composition. However, when aerosols, thick hazes, or clouds are considered, this slope changes considerably because of the change in the wavelength dependence of the scattering (see Eq. (25)), which changes the atmospheric optical depth and effective height of the atmosphere (also see Wunderlich et al. 2020).

Fig. 13 Evolution/escape of atmospheric species in terms of the atmospheric mass (left axis) for the duration of the magma ocean (x-axis) for a reduced mantle case (upper panel) and an oxidized mantle case (lower panel). The colored solid lines refer to atmospheric masses integrated over the whole planet that change due to outgassingand chemical speciation (see text). The dotted brown line shows the residual mass of H2 in the atmosphere after accounting for outgassing and DL escape. The shaded region in cyan shows the mass loss of H2 in the range for S = 10 to S = 100 (XUV relative to the modern value) using the Zahnle et al. (2019) best-fit formalism (see text). The shaded region in yellow shows the energy-limited loss of H2 for the lower range S = 10 to upper range S = 100. The dotted blue horizontal line shows the hydrogen content of one (modern) Earth ocean (MHO = 1.6 × 1023 g). The filled circles indicate the escaping mass of H2 obtained using DL. The colors indicate the respective mass-loss rates (right axis).

Fig. 14 Hydrogen escape and outgassing rates as a function of redox state of the mantle for scenario 2.2 (EOM). The solid red line shows the H2 outgassing rate at the end of magma ocean. The closed red circles show the diffusion-limited escape of H2 for the EOM. The shaded region in yellow shows the loss rate due to energy-limited escape within the range S = 10 (lower boundary) to S = 100 (upper boundary) i.e. XUV relative to the modern value. The vertical dotted line marks the region separating a high H2 VMR to low H2 VMR and illustrates the most efficient escape process occurring. The present Earth H2 degassing rate with an uncertainty range (Ramirez et al. 2014) is marked for reference.

5.3 Outgassing and escape: loss timescales

Our results suggest very different regimes in terms of outgassing of various species and escape of H₂ depending on the mantle redox state during the MO period. For the scenarios that consider the reducing and oxidizing state of the mantle, the outgassing rate of H₂ is greater than the escape rate, and there is enough interior outgassing of H₂ to form an atmosphere. A similar result was suggested by Ramirez et al. (2014) for early Earth and Mars.

Because the VMR of H₂ is high (e.g., in a reducing atmosphere), the escape of H₂ most likely occurs by energy-driven, that is, energy-limited, escape. On the other hand, if less hydrogen is present (e.g., in an oxidizing atmosphere), the diffusion of hydrogen through a heavy gas to reach the upper atmospheric regions and escape is limited. Therefore, the atmosphere would probably enter a DL configuration. Thus, depending upon the abundance of H₂ in the atmosphere, our results suggest a switch from energy-limited (EL) escape (high VMR) of H₂ to DL escape (low VMR), thus complying with Zahnle et al. (2019).

For a reduced early atmosphere with 102 bar of H₂ by the end of MO, we estimate the timescale of H₂ atmosphere removal to be ~10 Myr using both DL and EL (S = 10) escape. On the other hand, because the H₂ mass-loss rate and the surface partial pressure of surface hydrogen are lower (2.7 bar) for an oxidized atmosphere, the mass-loss timescale is estimated to be ~16 Myr assuming DL escape. Similar mass-loss timescales (~10⁶–10⁷ yr) for surface pressures of 10 to 100 bar are obtained by hydrodynamic escape modeling of Pahlevan et al. (2019). The results for total hydrogen mass-loss obtained via two different parameterized escape models (this work) and hydrodynamic model (Pahlevan et al. 2019) are therefore approximately similar. Moreover, as we showed in Fig. 13, the mass loss obtained using the energy-limited (hydrodynamic) approach following Zahnle et al. (2019) (cyan shaded region) approximates the diffusion-limited mass loss at higher levels of irradiation, S = 100 (also see Lammer et al. 2018). Furthermore, it is worth stating that several works (e.g., Lammer et al. 2011) have noted the potential importance of EUV in driving energy-limited escape on early Earth and speculated that this quantity is not known to within a factor 100 or more times the modern-day value depending upon whether the early Sun was a fast or slow rotator.

5.4 Clouds

Clouds could be rather common features in extrasolar atmospheres, and they could have a potentially strong effect upon the atmosphericspectra, climate, and photochemistry. Spectral features can be significantly reduced depending upon the extent (layer location andthickness) and the microphysical properties (size, shape, number density distribution, and refractive index) of the clouds occurring over a wavelength range that is sensitive to the cloud diameter, for example. The general effect of clouds upon atmosphericspectra and retrieval by the James Webb Space Telescope (JWST) was discussed by Mai & Line (2019) and Fauchez et al. (2019).

In magma ocean worlds, thick water cloud layers could form (see, e.g., Lebrun et al. 2013; Marcq et al. 2017) as the outgassed steam atmospheres cool. These cloud layers could help prevent the loss of planetary infrared radiation to space, and if thick water clouds were to extend over the region sampled by transmission spectroscopy, for instance, they might significantly reduce the strength of the spectral features detected.

In Fig. 15 we show the atmospheric transit depth calculated using Eq. (27) for planet Earth orbiting an M-type star (0.4 solar radii) in the top panel and a G-type star (bottom panel) for two different redox states of the mantle, IW+4 and IW. A clear increase in transit depth is visible for the IW buffer with 50% H₂O and 0.5% H₂ because the molecular mass of the atmosphere is lower than in the pure H₂O case (IW+4 buffer). Results also suggest a substantial increase in the transit depth when observing via an M-type star because its stellar radius is smaller, thus favoring the possibility of detecting atmospheric species. We also compared the cloud-free scenarios with the cloudy simulations. With a cloud deck (see Sect. 2.4), the spectral features are weakened (depending upon the cloud cross section), as shown in Fig. 15. The results of the comparison between cloud-free and cloudy scenarios (α = 0 and enhanced scattering cross section as compared to Rayleigh) agree well with Turbet et al. (2019). The transit depths in the optically thick cloudy scenarios can be considered as maximum values (because thick clouds effectively block shortwave radiation) compared with the smaller transit depths of the clear atmosphere cases without clouds. A sensitivity study by Moran et al. (2018) that involved increasing the haze-scattering cross section also suggested a weakening in spectral features when a global layer of Rayleigh-scattering haze was added to the TRAPPIST-1 d, e, and f planets. A more involved sensitivity study including the geometry of the cloud, such as its cross section (Kitzmann et al. 2011a,b), requires the use of detailed cloud microphysics (Zsom et al. 2012) that is applicable for Earth-like exoplanets. We did not explore this here and aim to investigate it in future related studies.

Fig. 15 Atmospheric transit depth of a cloud-free planetary atmosphere around an M star (top) and G star (bottom) at the EOM phase for an atmosphere simulated from IW+4 mantle buffer (red) and IW mantle buffer (blue) is compared with the transit depth for cloudy scenarios (α = 0 and enhanced cross-section as compared to the Rayleigh scattering).

5.5 Application to exoplanets

Several hundred rocky exoplanets classified as super-Earths with mass ranging between M = 1−8M_E and radius R = 1−1.8R_E are being discovered by various current space missions such as TESS¹ and will be discovered by future missions such as PLATO², JWST³, and CHEOPS⁴. Atmospheric characterization of these planets via spectroscopic analysis can benefit greatly from interior geochemical models that consider outgassing during the magma ocean phase (Bower et al. 2019) and volcanic outgassing after the magma ocean phase (Ortenzi et al. 2020). Ortenzi et al. (2020) showed that the redox state of the mantle affects the atmospheric thickness and overall evolution, thus constraining the interior chemical state of rocky exoplanets. Puffy outgassed atmospheres for hot, molten Earth-like planets often lead to an effect known as “radius inflation” that depends upon the initial volatile inventory of H₂O and CO₂ (Bower et al. 2019) and on the irradiation that they receive from the star (Turbet et al. 2019).

Furthermore, as shown in Fig. 15, the redox state of the mantle and initial inventory of volatiles can induce strong variations in the transit depths for planets orbiting G and M stars. This would translate into changes in the planetary radius to be observed by missions such as TESS and PLATO (Rauer et al. 2014). However, the detectability of molecular spectral features may be severely affected by the presence of clouds (see e.g. Fauchez et al. 2019).

Kaltenegger et al. (2007) modeled the observable spectra of Earth-like planets about two billion years ago to a present-day atmosphere. They reported a number of constraints, such as the resolution of the telescopes for detecting the molecular species in the atmosphere. It is therefore important to characterize the atmosphere of Earth through its geological history and that our study present useful insights on this. Recent interior-atmosphere modeling studies have shown that pure steam-based atmospheres of highly irradiated exoplanets orbiting M dwarfs such as GJ 1214 b (Schaefer et al. 2016) and LHS 3844 b are subject to strong XUV of the star and could lose the outgassed atmosphere and become a bare rock (Kreidberg et al. 2019). See also Kite & Barnett (2020) for a detailed study of atmospheric loss and subsequent revival relevant for exoplanets. In all these models, the water is dissociated into hydrogen and oxygen, which escape to space, and some remaining oxygen is dissolved into the magma ocean or remains in the atmosphere (e.g., Schaefer et al. 2016).

We have provided important testable predictions to inform future space missions such as the JWST, Ariel⁵ (Tinetti et al. 2018; Edwards et al. 2019), LIFE⁶ (Quanz et al. 2019), and ELT⁷ (Maiolino et al. 2013) that will characterize planetary atmospheres in terms of detecting species such as CO, CO₂, and H₂O in spectral observations of Earth-sized planets.

Planets in close-in orbits around cooler stars are favored targets in exoplanetary science because of their improved contrast ratios, higher transit frequency, and increased geometric transit probability (Scalo et al. 2007; Shields et al. 2016), for example. Furthermore, the detectability of atmospheric spectral signals from MO exoplanets (Ito et al. 2015) and thick atmospheres with H₂O and CO₂ with a high sensitivity and large spectral coverage of the JWST (NIRSpec, MIRI) instruments (Ducrot et al. 2020) would bring new possibilities to constrain such atmospheres.

6 Conclusions

The main conclusions of this study are listed below.

• 1.

  The redox state of the mantle can affect not only atmospheric composition, but also the vertical temperature structure and therefore the mixing and the hydrological cycle. These interior-atmospheric couplings could be important for steam-condensation timescales and accordingly, for the formation of the Earth oceans.

• 2.

  Reduced atmospheres emit thermal radiation more strongly, which leads to faster cooling than in oxidized atmospheres because the latter feature strong absorbers such as H₂O.

• 3.

  Reduced atmospheres (CO and H₂) have much larger spectral features in transmission than oxidized atmospheres (H₂O and CO₂) because the scale heights are larger.

• 4.

  Thick optical clouds and hazes can mute or weaken the spectral signatures of molecules in the transmission spectra and enhance the transit depth because of the cloud cover. It therefore becomes difficult to probe cloudy atmospheres. On the other hand, clouds absorb much of the radiation and reduce the outgoing longwave radiation by 1–2 orders of magnitude (Marcq et al. 2017), leading to longer magma ocean cooling timescales.

• 5.

  The interplay between outgassing and escape suggests that outgassing of hydrogen proceeds faster than the atmospheric escape of H₂ and that enough outgassed hydrogen is available to form an atmosphere over the magma ocean duration of ~1 Myr. The timescale for total H₂ mass loss, however, is estimated to be within some dozen million years. Furthermore, a more sophisticated hydrodynamical code to study the loss of H₂ and the heavy atmospheric species H₂O and CO₂ is desirable and planned for future studies.

Acknowledgements

We thank the anonymous reviewer for their valuable comments that greatly enhanced the quality of the manuscript. N.K. and G.O. acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Project-ID 263649064 – TRR 170. This is TRR 170 Publication No. 117. M.G. acknowledges funding from DFG – Project-ID GO 2610/1-1. F.S. acknowledges DFG project SCHR 1125/3-1.

Appendix A Gas constant for a mixture

When a mixture of two gases a and b, with partial pressure p_a and p_b, exists at the same volume V and temperature T, the mixture behaves like a perfect gas, and according to Dalton’s law of partial pressures, $$(p_a+p_b)V=(m_a{R}_a+m_b{R}_b)T.$$(A.1)

Because the mixture behaves like a perfect gas, for m = nM moles of a gas, the total pressure is given by $$\text{PV}=m{R}_{\text{mix}}T,$$(A.2)

where R_mix is the gas constant for the mixture. From Eqs. (A.1) and (A.2), mR_mix = m_aR_a + m_bR_b. Hence, $$R_{\text{mix}}=\frac{m_a{R}_a+m_b{R}_b}{m}.$$(A.3)

In terms of the molar volume fraction, $$R_{\text{mix}}=\frac{n_a{M}_a{R}_a+n_b{M}_b{R}_b}{{\sum }_i{n}_i{M}_i}.$$(A.4)

Because $MR=\overline{R}$, where $\overline{R}$ is the universal gas constant (=8.31415 J kg⁻¹ K⁻¹), this equation can be written as $$R_{\text{mix}}=\frac{n_a\overline{R}+n_b\overline{R}}{n_a{M}_a+n_b{M}_b}.$$(A.5)

For a single-component gas a, $R=\overline{R}/M_a$. For a two-component gas mixture, Eq. (A.5) can be written as $$R_{\text{mix}}=\overline{R}\left(\frac{n_a+n_b}{\mu }\right),$$(A.6)

where μ = n_aM_a + n_bM_b is the mean molecular weight of the mixture. For an ith-component mixture of gases, $$R_{\text{mix}}=\overline{R}\left(\frac{{\displaystyle \sum n_i}}{{\displaystyle \sum n_i}{M}_i}\right).$$(A.7)

Appendix B Volatile speciation model output

Appendix C OLR of H₂ atmospheres

Fig. C.1 Top: thermal emission spectra for an H2-dominated atmosphere overlying a highly reduced mantle (IW-4) for a varying surface partial pressureof hydrogen as shown in the legend. Bottom: OLR plotted against p${}_{{\text{H}}_2}$ for a fixed surface temperature of Ts = 2800 K showing a decrease in OLR with increase in H2 surface partial pressure.

References

All Tables

All Figures

	Fig. 1 Schematic showing coupled interior-atmosphere exchange via outgassing of reduced species such as H2 and CO (orange) for a reduced mantle (e.g., IW buffer) and oxidized species such as H2O and CO2 (blue) for an oxidized mantle (e.g., QFM buffer) as the magma ocean solidifies. The H2 in the atmosphere is lost to space via escape.
In the text

Fig. 2 Volatile chemical speciation in terms of outgassed mole fraction vs. oxygen fugacity of the mantle. The x-axis indicates the oxygen fugacity range in logarithmic units relative to the IW buffer. For reference (in log units), buffers QIF = IW-1, QFM = IW+3.8, and NiNiO = IW+4.2 approximately (Wood et al. 1990). From top to bottom: BOM, active MO phase, EOM and a steam-atmosphere phase that are characterized by different initial H/C between hydrogen and carbon. The corresponding input values of the initial volume mixing ratio of volatiles, surface pressures Ps, and temperatures Ts (Nikolaou et al. 2019) for the speciation model are indicated at the top of each panel.

In the text

	Fig. 3 Specific heat capacities as a function of temperature for volatile species H2O, CO, and CO2 on the left x-axis and H2 on the right x-axis as obtained from the Shomate equation described in the main text.
In the text

Fig. 4 Left: pressure-temperature profiles for the BOM scenarios 1.1 (upper), 1.2 (middle), and 1.3 (lower panel) for three different buffers IW-4, IW, and IW+4. The dashed red line represents the saturation water vapor curve. Pressures at the bottom of atmosphere values are taken from Table 4. Right: altitude-temperature profiles for scenarios 1.1 (upper), 1.2 (middle), and 1.3 (lower panel). We note the different altitude ranges for the panels on the right.

In the text

	Fig. 5 Same as Fig. 4, but for p, T conditions and scenarios as summarized in Table 5, representing the EOM.
In the text

	Fig. 6 Temperature profiles p − T for scenario 3.1 (upper panel) with an initially H2O-dominated atmosphere and scenario 3.2 (lower panel) with an initially CO2-dominated atmosphere (Table 6). The molar fraction of each species is shown in the legend. “No speciation” indicates that the values were input directly by the user and were not calculated by the C-O-H speciation model.
In the text

Fig. 7 Comparison of thermal emission spectra for a reduced vs. an oxidized mantle during the BOM presented for scenario 1.1 (upper panel), 1.2 (middle panel), and 1.3 (lower panel) of Table 4. The output of the volatile speciation model is shown in the legend. Values written above the legend show the initial input values to the speciation model. Reference blackbody curves for four different temperatures are plotted. Absorption features ofkey chemical species are indicated. The spectra shown here have been binned to a resolution of λ∕Δλ = 1000.

In the text

	Fig. 8 Same as Fig. 7, but for scenarios of Table 5. The blackbody curves for three different surface temperatures bracketing the emission spectra are plotted for reference.
In the text

Fig. 9 Thermal spectral emission for a reducing atmosphere overlying a reduced mantle (IW-4) for scenarios 2.1, 2.2, and 2.3 of Table 5. The legend shows the mole fraction of the outgassed species obtained from the speciation model. Prominent features of absorbing species are seen at wavelengths marked in the figure. Reference blackbody curves for three different surface temperatures are plotted. The insets show the zoomed-in thermal spectra with the most prominent absorption features of CO2, CO and H2. The spectra have been binned to a resolution of λ∕Δλ = 1000.

In the text

	Fig. 10 Same as for Fig. 9, but for an oxidizing atmosphere overlying an oxidized mantle (IW+4). Notable features are the emergence of prominent H2O features between 1 and 2 μm (scenario 2.3) and the disappearance of CO feature at 4.6 μm as compared to Fig. 9 (scenario 2.1 and 2.2).
In the text

Fig. 11 Left: pressure at the bottom of the atmosphere pboA for various cases referring to different initial volatile abundances and three different mantle buffers: reduced (IW-4), IW, and oxidized (IW+4). The details of speciation for various cases on the x-axis are shown in Table B.1 (Ts = 3300 K) and Table B.2 (Ts = 1650 K) in the Appendix B. The horizontal dashed black line indicates the total surface pressure of the initial outgassed volatile for the various cases. Right: outgoing longwave radiation (OLR) for various cases corresponding to the atmospheric surface pressure shown in the left panels and for three different buffers.

In the text

Fig. 12 Transmission spectra calculated by GARLIC shown as atmospheric effective height vs. wavelength for the BOM scenarios (left panels) and EOM scenarios (right panels) for three different redox buffers IW-4, IW, and IW+4.The mole fraction of final outgassed species in the atmosphere for various scenarios is the same as provided inTables 4 and 5. The effective height of a reduced atmosphere (with low molecular weight) for the case of reduced mantle IW-4 (blue curve; scenarios 1.2, 1.3, 2.2, and 2.3) is lower than an atmosphere with larger molecular weight (blue curve; scenarios 1.1 and 2.1) and oxidizing atmospheres (green curve; all scenarios). For more details, see the text. The spectra have been binned to a resolution of λ∕Δλ = 1000.

In the text

Fig. 13 Evolution/escape of atmospheric species in terms of the atmospheric mass (left axis) for the duration of the magma ocean (x-axis) for a reduced mantle case (upper panel) and an oxidized mantle case (lower panel). The colored solid lines refer to atmospheric masses integrated over the whole planet that change due to outgassingand chemical speciation (see text). The dotted brown line shows the residual mass of H2 in the atmosphere after accounting for outgassing and DL escape. The shaded region in cyan shows the mass loss of H2 in the range for S = 10 to S = 100 (XUV relative to the modern value) using the Zahnle et al. (2019) best-fit formalism (see text). The shaded region in yellow shows the energy-limited loss of H2 for the lower range S = 10 to upper range S = 100. The dotted blue horizontal line shows the hydrogen content of one (modern) Earth ocean (MHO = 1.6 × 1023 g). The filled circles indicate the escaping mass of H2 obtained using DL. The colors indicate the respective mass-loss rates (right axis).

In the text

Fig. 14 Hydrogen escape and outgassing rates as a function of redox state of the mantle for scenario 2.2 (EOM). The solid red line shows the H2 outgassing rate at the end of magma ocean. The closed red circles show the diffusion-limited escape of H2 for the EOM. The shaded region in yellow shows the loss rate due to energy-limited escape within the range S = 10 (lower boundary) to S = 100 (upper boundary) i.e. XUV relative to the modern value. The vertical dotted line marks the region separating a high H2 VMR to low H2 VMR and illustrates the most efficient escape process occurring. The present Earth H2 degassing rate with an uncertainty range (Ramirez et al. 2014) is marked for reference.

In the text

	Fig. 15 Atmospheric transit depth of a cloud-free planetary atmosphere around an M star (top) and G star (bottom) at the EOM phase for an atmosphere simulated from IW+4 mantle buffer (red) and IW mantle buffer (blue) is compared with the transit depth for cloudy scenarios (α = 0 and enhanced cross-section as compared to the Rayleigh scattering).
In the text

	Fig. C.1 Top: thermal emission spectra for an H2-dominated atmosphere overlying a highly reduced mantle (IW-4) for a varying surface partial pressureof hydrogen as shown in the legend. Bottom: OLR plotted against p${}_{{\text{H}}_2}$ for a fixed surface temperature of Ts = 2800 K showing a decrease in OLR with increase in H2 surface partial pressure.
In the text