Toward RNA Life on Early Earth: From Atmospheric HCN to Biomolecule Production in Warm Little Ponds

In Figure 1, we illustrate our early Earth atmospheric models by focusing on the main sources and sinks for the key molecular species relevant to controlling and determining HCN chemistry. As an overall principle, H₂ and CO₂ are the main species that determine whether the environment is reducing or oxidizing. It is the balance between these two species that determines the concentration of CH₄, the main precursor to HCN.

Zoom In Zoom Out Reset image size

In the top left panel, impact degassing produces H₂. The degassing rate at each epoch is calculated by combining equilibrium H₂ production rates from enstatite chondrite impactors via the reaction Fe + H₂O → FeO + H₂ (Zahnle et al. 2020) with the bombardment rate on early Earth based on mathematical fits to the observed lunar cratering record (Chyba 1990; Pearce et al. 2017). The main sinks for H₂ include UV photodissociation, hydrodynamic escape to space, and disequilibrium chemistry.

In the top right panel, the main source of CO₂ on early Earth is volcanic outgassing (Zahnle et al. 2020). We use a constant Earth-like volcanic CO₂ outgassing rate in all our models (Hu et al. 2012). The main sinks for CO₂ are photodissociation in the upper atmosphere and rain-out in the lower atmosphere.

We utilize a potentially abundant source of CH₄ to the Hadean atmosphere from serpentinization and Fischer–Tropsch-type (FTT) synthesis. This begins with water-dependent processes in hydrothermal systems wherein Fe- and Mg-rich ultramafic rocks (e.g., olivine) in mid-ocean ridges and forearc systems produce H₂. Then, H₂ reacts with the aqueous CO₂ in these environments in the presence of mineral catalysts to produce CH₄ (Holm et al. 2015; Guzmán-Marmolejo et al. 2013). Abiotic methane production has been observed in hydrothermal systems (Fiebig et al. 2007; Bradley & Summons 2010); however, experiments of FTT synthesis from olivine typically produce very low yields (Krissansen-Totton et al. 2018). Naturally occurring catalysts such as awaruite or chromite greatly speed up FTT synthesis (Bradley 2016); however, the availability of these catalysts in primordial hydrothermal systems is still somewhat uncertain. Equilibrium models suggest this hydrothermal process can sustain ∼2–2.5 parts per million (ppm) of CH₄ in the early atmosphere (Guzmán-Marmolejo et al. 2013). We use the calculated CH₄ outgassing rate from these models (Guzmán-Marmolejo et al. 2013). The main sinks for methane are UV photodissociation and disequilibrium chemistry.

Finally, the main source of HCN is photodissociation or lightning dissociation of species such as N₂ and CH₄ followed by radical chemistry. In the case of lightning, this radical chemistry takes place at high temperature in the lightning channel. The main sinks for HCN are UV photodissociation and rain-out. For further details on the source and sink rates, see Appendix A.

We model the early Earth atmosphere during its reducing phase at 4.4 bya for calculated habitable surface temperatures of 78°C (Model A) and 51°C (Model C), as well as its oxidizing phase at 4.0 bya for calculated surface temperatures of 51°C (Model B) and 27°C (Model D). These particular surface temperatures result from the radiative transfer calculations for our chosen starting atmospheric compositions. These models differ in atmospheric composition, solar luminosity, UV irradiation intensity, HCN and radical production from lightning, and impact bombardment rate (see Table 1 for model details).

Table 1. Summary of the Four Early Earth Atmospheric Models in This Work

Model	Description	Date (bya)	Ps (bar) a	Ts (°C) b	Molar Composition	Surface Flux (cm–2 s–1)	Lightning (cm–2 s–1) c
A	Early Hadean	4.4	1.5	78	H2: 90%	H2: 2.3 × 1011	H: 4.1 × 106
	(reducing)				N2: 10%	CO2: 3.0 × 1011	CO: 6.3 × 103
					CH4: 2 ppm	CH4: 6.8 × 108	OH: 4.4 × 103
					H2O: Figure A1	H2O: 2.0 × 109	NO: 5.0 × 101
							3O: 4.4 × 100
							4N: 7.2 × 10−1
							$\mathrm{HCN}$: 4.4 × 10−1
							O2: 2.2 × 10−2
							CH3: 1.6 × 10−3
B	Late Hadean	4.0	2	51	CO2: 90%	H2: 2.3 × 1010	CO: 1.8 × 105
	(oxidizing)				N2: 10%	CO2: 3.0 × 1011	O2:8.5 × 104
					CH4: 10 ppm	CH4: 6.8 × 108	NO: 7.4 × 103
					H2O: Figure A1	H2O: 2.0 × 109	OH: 3.5 × 103
							3O: 5.5 × 102
							H: 9.7 × 101
							4N: 3.2 × 10−3
							$\mathrm{HCN}$: 2.9 × 10−6
C	Early Hadean	4.4	1.13	51	H2: 90%	H2: 2.3 × 1011	H: 1.1 × 105
	(reducing)				N2: 10%	CO2: 3.0 × 1011	OH: 1.3 × 102
					CH4: 1 ppm	CH4: 6.8 × 108	CO: 8.1 × 101
					H2O: Figure A1	H2O: 2.0 × 109	NO: 1.6 × 100
							3O: 1.6 × 10−1
							4N: 2.2 × 10−2
							$\mathrm{HCN}$: 4.1 × 10−3
							O2: 7.4 × 10−4
							CH3: 1.3 × 10−5
D	Late Hadean	4.0	2	27	CO2: 90%	H2: 2.3 × 1010	CO: 1.8 × 105
	(oxidizing)				N2: 10%	CO2: 3.0 × 1011	O2: 8.5 × 104
					CH4: 1.5 ppm	CH4: 6.8 × 108	NO: 7.4 × 103
					H2O: Figure A1	H2O: 2.0 × 109	OH: 3.5 × 103
							3O: 5.5 × 102
							H: 9.7 × 101
							4N: 3.2 × 10−3
							$\mathrm{HCN}$: 2.9 × 10−6

Notes. Initial compositions for each model are chosen to (a) align with typical assumptions for reducing (H₂-dominant) or oxidizing (CO₂-dominant) conditions at the chosen epoch, and (b) yield calculated P–T profiles with habitable surface conditions (i.e., 0 °C ≤ T_s ≤ 100 °C) (see Figure A2 for P–T profiles).

^a P_s : surface pressure. ^b T_s : surface temperature. ^{c 3}O and ⁴N represent the triplet and quartet spin states for oxygen and nitrogen atoms, respectively. Our chemical network also includes singlet oxygen (¹O) and doublet nitrogen (²N).

Download table as:  ASCIITypeset image

ChemKM is a spatially 1D chemical kinetic model for disequilibrium atmospheric chemistry calculations that makes use of the Double-precision Livermore Solver for Ordinary Differential Equations (DLSODE) from the ODEPACK collection (Stepleman 1983). The error for this solver is controlled by the relative error tolerance and the absolute error tolerance, which are set to 10⁻⁵ and 10⁻⁹⁹ respectively to ensure numerical stability.

ChemKM computes atmospheric chemistry based on input chemical network, P–T profile, eddy diffusion (mixing) profile, top-of-atmosphere (TOA) solar radiation, wavelength-dependent photochemical reactions, and influxes and outfluxes of species at the surface and TOA. ChemKM has been benchmarked with several other chemical kinetic codes, ⁷ and has been used in the past to simulate the atmosphere of Titan (Pearce et al. 2020b), as well as those of cold, hot, and ultrahot gaseous exoplanets (Molaverdikhani et al. 2019a, 2020a). Atmospheric rain-out was newly developed in ChemKM for this work. Simulations were run on the cluster for approximately one week. In this time, our models reached 1–50 million years of simulated time.

CRAHCN-O is a consistent reduced atmospheric hybrid chemical network now containing 259 two- and three-body reactions for the production of HCN and H₂CO in atmospheres dominated by any of H₂, CO₂, N₂, CH₄, and H₂O. We introduce 28 new reactions to CRAHCN-O in this work, in order to avoid the atmospheric buildup of species that previously had no reaction sinks (see Tables A1 and A2 for details). We have tested an oxygenless version of this network (CRAHCN) by modeling HCN production in Titan's atmosphere, and our computed HCN profile agreed very well with the Cassini observations (Pearce et al. 2020b).

petitRADTRANS is a 1D radiative transfer code based on the correlated-k method for gas absorption and the Guillot temperature model (Guillot 2010; Mollière et al. 2019). It is typically used to model exoplanet atmospheres to obtain transmission and thermal spectra (e.g., Molaverdikhani et al. 2019a; Mollière et al. 2020; Wang et al. 2020). We build upon its existing functionality to calculate P–T profiles self-consistently with tropospheric water vapor in sequence with the Arden Buck equation (Buck 1981).

P–T structure calculations are performed using petitRADTRANS (Mollière et al. 2019), and atmospheric chemistry is calculated using ChemKM (Molaverdikhani et al. 2019a, 2020a) coupled with an updated version of the CRAHCN-O chemical network (Pearce et al. 2020a, 2020b). We are the first to calculate composition-dependent P–T profiles for modeling HCN chemistry in the early Earth atmosphere. Past models have used a general habitable P–T profile or estimated the surface temperature using an analytic equation for a moist adiabat (Zahnle 1986; Tian et al. 2011; Zahnle et al. 2020). The code used to compute P–T profiles is available at https://gitlab.com/mauricemolli/petitRADTRANS.

For complete details on our atmospheric models, see Appendix A.

We follow the thermodynamic treatment from Chameides & Walker (1981) for the lightning production of HCN and other species on early Earth. Based on our initial atmospheric compositions, we calculate the equilibrium abundances of $\mathrm{HCN}$, H₂, N₂, H₂O, CO₂, CH₄, O₂, NO, OH, H, ⁴N, CO, ³O, and CH₃, in a 1 cm² lightning channel extending through the lowest layer in our atmospheres at a freeze-out temperature T_F = 2000 K. We use thermochemical data from the JANAF tables (Stull & Prophet 1971), and the ChemApp Software library for Gibbs free energy minimization (distributed by GTT Technologies, http://gtt.mch.rwth-aachen.de/gtt-web/).

We then use the resultant mixing ratios to calculate the influx of each of these species into the lowest layer of our atmospheres. These species were chosen because they are dominant equilibrium products in the early Earth lightning models by Chameides & Walker (1981). The freeze-out temperature (T_F ) was chosen to most accurately model the HCN produced in a lightning strike, as this is the key species of interest in this paper. Although freeze-out temperatures typically range from 1000–5000 K across species, one freeze-out temperature must be chosen to conserve the elemental abundances in the lightning strike. A nonequilibrium approach was also considered; however, an extensive high-temperature (up to 30,000 K) chemical kinetic network would be required and is perhaps unnecessary given the equilibrium timescale of <1 μs above 10,000 K compared to the ∼10 μs timescale of eddy diffusion and the ∼100 ms cooling time of a lightning channel (Hill et al. 1980).

The main assumption of our models is that the surface of the Earth maintained habitability (i.e., 0°C < T < 100°C), which is key for the presence of WLPs and the origin of life. We begin with assumed reducing and oxidizing atmospheric compositions for the early and late Hadean, respectively, and calculate the initial P–T profiles and tropospheric water vapor based on these compositions. We adjust both initial methane concentration and surface pressure to obtain calculated temperature profiles that fall within the habitable range. We smooth the initial water profiles from our calculations to 1% at the surface, and include ocean–atmosphere coupling by imposing an ocean evaporation rate of 2 × 10⁹ cm⁻² s⁻¹ to maintain a water mixing ratio of ∼0.1%–1% at the surface.

In Table 1, we summarize the four early Earth atmospheric models in our study. We model two epochs that vary in atmospheric composition, solar luminosity, UV irradiation intensity, and asteroid bombardment rate. These models correspond to the early Hadean at 4.4 bya (billion years ago) and the late Hadean at 4.0 bya. We compute two habitable P–T profiles for each model by slightly adjusting the methane content and/or surface pressure. The luminosity, UV intensity, and asteroid bombardment rate at each epoch are based on stellar evolution models (Heller et al. 2020; Baraffe et al. 2015), observations of solar analogs (Ribas et al. 2005), and the lunar cratering record (Chyba 1990), respectively.

Our atmospheric models are coupled (via rain-out) with the sources and sinks WLP model we first developed in Pearce et al. (2017). Biomolecule abundances are described by first-order linear differential equations and are solved numerically. The evolving concentrations of nucleobases, ribose, formaldehyde, and 2-aminooxazole in our WLP models are driven by the rate of incoming HCN from rain-out, and biomolecule losses due to UV dissociation, seepage, and hydrolysis. Given that experimental reaction rates are fast (≲ a few days), we apply experimental reaction yields to our HCN pond concentrations in order to estimate the pond concentrations of formaldehyde, nucleobases, ribose, and 2-aminooxazole. However, we recognize that at the lower pond concentrations here, these reactions could take much longer. The code used to compute biomolecule concentrations in ponds is available at https://github.com/bennski/Wet_Dry_Cycling_Pond_Model.

We utilize the experimental result that UV irradiation of liquid water produces solvated electrons, enabling a chemical pathway from HCN to formaldehyde (H₂CO) (Yi et al. 2020) in ponds. H₂CO can also enter ponds directly from the atmosphere (Pinto et al. 1980). Aqueous solutions containing HCN and H₂CO can produce nucleobases (i.e., adenine, guanine, cytosine, uracil, thymine) (Oró 1961; Larowe & Regnier 2008; Ferus et al. 2019), which are the base-pairing components of RNA and DNA, as well as ribose (Butlerow 1861; Breslow 1959), which binds with phosphate to make up the RNA backbone.

Furthermore, irradiated and wet–dry cycled or flowing solutions of HCN in the presence of phosphorus and dissolved salts enable the production of 2-aminooxazole, a key intermediate in the Powner–Sutherland pathway for producing the pyrimidine building blocks of RNA (cytidine and uridine monophosphate) (Powner et al. 2009; Ritson et al. 2018; Yi et al. 2020). Finally, Becker et al. (2018) recently presented a pathway to RNA nucleosides that involves the wet–dry cycling of solutions containing HCN and other atmospheric precursors. Such prebiotic chemistry experiments and models are based on the assumption that species such as HCN and H₂CO would be present and concentrated in WLPs on early Earth.

We considered including formose-like reactions occurring from more complex aldehydes (e.g., glycolaldehyde, glyceraldehyde) (Cassone et al. 2018); however, given the greater complexity of these species, it is expected that they would be produced in the atmosphere in much lower concentrations than HCN and H₂CO. However, it would be valuable to explore reactions that produce these more complex aldehydes in planetary atmospheres so that these sugar precursors can be included in a future model.

Lastly, we understand that straightforwardly applying experimental biomolecule production yields from HCN to obtain estimates of biomolecule influxes in our WLP model has issues given the high reactant concentrations and ideal conditions of each experimental setup; however, given that there are no chemical kinetic rate coefficients for the aqueous production of nucleobases, ribose, and 2-aminooxazole, we utilize experimental yields in our model to obtain reasonable first-order estimates. We emphasize that these biomolecule concentrations should be understood as upper bounds in the absence of any concentrating mechanism beyond evaporation.

In Table A6 we summarize the sources and sinks of our pond models. See Appendix A.2 and Pearce et al. (2017) for complete details regarding these models.

We have found that ponds that are roughly 1 m in radius and depth and are cylindrical in shape are an optimal fiducial estimate for subsequent RNA polymer synthesis by wet–dry cycles (Pearce et al. 2017). We use the "intermediate" hot early Earth environment from Pearce et al. (2017), which is based on the seasonal sinusoidal precipitation rates in Indonesia (Reichle et al. 2011; Berghuijs & Woods 2016). Precipitation coupled with evaporation and seepage produces a natural wet–dry cycle within the pond that has a ∼6 month wet phase followed by a ∼6 month dry phase. Various pond environments were explored in Pearce et al. (2017), and were found to produce similar results in terms of peak nucleobase concentrations.

In Figures 2(A), (B), (C), and (D), we plot the temporal evolution of atmospheric HCN and H₂CO in these four early Earth models. To give some context for what would be considered high atmospheric HCN abundances, Cassini observed HCN mixing ratios in the heavily reducing atmosphere of Titan to be ∼10⁻⁷–10⁻⁶ near the surface (150–300 km) and ∼10⁻⁴–10⁻² in the upper atmosphere (700–1050 km) (Magee et al. 2009; Vinatier et al. 2010; Adriani et al. 2011; Koskinen et al. 2011).

Zoom In Zoom Out Reset image size

In the early Hadean (reducing) model A, HCN mixing ratios decrease from 10⁻¹¹ at the surface to 10⁻¹⁴ from t = 100 yr to 400 yr. Then, HCN slowly builds back up again to 2.7 × 10⁻¹¹ over the next 50 million years. Moving up in altitude, we see the main region of HCN production from UV radiation at ∼500–600 km, which provides an HCN abundance of ∼10⁻⁷–10⁻⁶ at these altitudes after 10 million years. We isolate the predominant energy source for HCN production in Figures 2(E) and (F) by turning off UV photochemistry and lightning chemistry, respectively, for our early Hadean (reducing) model A.

In the late Hadean (oxidizing) model B, HCN mixing ratios decrease at the surface from 10⁻¹⁴ to 10⁻¹⁶ from t = 100 yr to ∼40,000 yr. Then, HCN slowly builds back up to 10⁻¹⁵ over the next 2 million years. UV production of HCN produces a peak abundance of ∼10⁻⁸ at ∼50 km at 2 million years.

One of the most important results is that neither oxidizing model produces nearly as much HCN as our reducing models. Our calculations reveal that HCN production near the surface is about four orders of magnitude more favorable in reducing conditions than it is in oxidizing conditions. It is worth noting that the oxidizing conditions at 4.0 bya have reduced gases such as H₂ in the ppb–ppm range. Even more oxidizing conditions could be present after 4.0 bya with the declining H₂ impact degassing rate allowing oxygen species to become more dominant.

We learned that UV photochemistry is crucial for atmospheric HCN production on early Earth (Figures 2(E) and (F)). A major result of our simulations is that without UV photochemistry, HCN would be over eight orders of magnitude less abundant at the surface during the early Hadean. Even with an increased lightning flash density representative of storms occurring during volcanic eruptions (1 × 10⁴ flashes km⁻² yr⁻¹, Hodosán et al. 2016), we see no enhancement in HCN concentration beyond the photochemical result (see Figure B4).

Formaldehyde, given these atmospheric models, would need to come from elsewhere, such as UV-driven, aqueous chemistry in WLPs. In comparison to HCN, H₂CO is much less abundant at the surface in our early Hadean (reducing) models. H₂CO builds up from ∼10⁻²³ to ∼10⁻¹⁶ over 50 million years. The mixing ratio for H₂CO is at its highest value of 10⁻⁹–10⁻⁸ in the middle of the atmosphere of these models after 50 million years. In the late Hadean (oxidizing) model B, H₂CO increases at the surface from ∼10⁻¹⁶ to ∼10⁻¹³ in ∼1 million years. We did not explore atmospheric production of H₂CO further given its considerably low abundances in all models.

The temporal evolution of the dominant atmospheric species at the lowest (surface) layer in the atmosphere is shown in Figure 3. The atmospheric rain-out rates for HCN and H₂CO from this layer provides the influx rates into the WLPs.

Zoom In Zoom Out Reset image size

One striking result of our atmospheric chemistry calculations is that the HCN and CH₄ mixing ratios are correlated after the first few hundred years in all our models. The relevant abundance curves are shown in bold in Figure 3 to emphasize this point, as the HCN and CH₄ abundances trace one another closely from t = 300 yr onward.

In Figure B2, we plot the molar ratio of HCN to CH₄ over time and find an average value of ∼(1–4) × 10⁻⁷ for the early Hadean (reducing) models and (3–5) × 10⁻⁶ for the late Hadean (oxidizing) models.

CH₄ and $\mathrm{HCN}$ surface abundances also follow the trend of the H₂ surface abundances in the late Hadean (oxidizing) models B and D. This is because H₂ drives the evolution of reducing and oxidizing conditions in our models. We see that CH₄ abundances are stable in reducing (high H₂) conditions, and that CH₄ abundances are depleted in oxidizing (high O₂ and OH) conditions.

There is an anticorrelation between $\mathrm{HCN}$ and O₂ in the late Hadean (oxidizing) models. This is because a high O₂ environment leads to the oxidation of carbon into species such as H₂CO rather than the reduction of carbon into HCN. This is also why we also see a correlation between O₂ and H₂CO in the late Hadean (oxidizing) models.

These correlations and anticorrelations are consistent with the atmospheric observations we see for Titan and present-day Earth, respectively. Titan's atmosphere is abundant in HCN (10⁻⁷–10⁻²) (Vinatier et al. 2010; Adriani et al. 2011; Koskinen et al. 2011; Magee et al. 2009) due to the high abundances of reducing gases such as CH₄ (∼5.7%) and H₂ (∼0.1%) and low concentrations of oxidizing gases such as CO₂ (10–20 ppb) and H₂O (0.5–8 ppb) (Catling 2015; Hörst 2017). On the other hand, HCN on Earth today is present in low abundances ∼10⁻¹⁰ (Cicerone & Zellner 1983) because of the high abundance of oxidizing gases in our atmosphere such as O₂ (21%), H₂O (0%–3%), and CO₂ (∼400 ppm) and modest ∼1 ppm levels of CH₄.

We tested the hypothesis presented in several exo-atmosphere studies that the C/O ratio plays a central role in controlling their chemical composition (Madhusudhan 2012; Mollière et al. 2015; Rimmer & Rugheimer 2019; Molaverdikhani et al. 2019a, 2019b, 2020b). In our early Hadean (reducing) models, the C/O ratio in the lowest atmospheric layer increases from 10⁻³ to 10⁻² in the first few thousand years, and then to 1 after ∼30 million years. In the late Hadean (oxidizing) models, the C/O ratio decreases from 0.5 to ∼10⁻² in the first few thousand years, and then increases to 0.02–0.03 over the next 1–2 million years. Evidently the C/O ratios in our models are set mainly by the increase and decrease of CO₂ and CO at the surface, and in the late Hadean (oxidizing) models, also by the fluctuation of O₂. Our models show that the abundances of key biomolecule precursor species CH₄ and HCN are not strongly dependent on the C/O ratio.

These results are not entirely surprising. In a study of the effects of C/O ratio on atmospheric HCN abundance, Rimmer & Rugheimer (2019) found that for a given C/O ratio, the presence of methane leads to considerably more HCN. Considering these results, we suggest that CH₄ abundance is a better fundamental driver for HCN production than the C/O ratio on its own.

In Figure 4(A), we display the concentrations of adenine in our model WLP from aqueous production for different HCN rain-out (influx) rates from our early Hadean (reducing) and late Hadean (oxidizing) atmospheric models (see Figure B1 for rain-out rates). Adenine concentrations are displayed as shaded regions to cover the range of experimental yields of adenine production from HCN. We also display for comparison the concentrations of adenine from meteoritic and interplanetary dust particle (IDP) delivery calculated using the same source/sink pond models in Pearce et al. (2017).

Zoom In Zoom Out Reset image size

Given the range of experimental conditions, there are additional uncertainties beyond the shaded regions that are not displayed here. Aqueously produced biomolecule concentrations could be lower than the abundances calculated here due to (1) the low nanomolar-range HCN pond concentrations compared to typically molar-range experimental concentrations, (2) differences in the radiation field used in experiments for photolytic H₂CO production versus the solar radiation incident on WLPs, and (3) uncertainties in the hydrolysis and photodestruction rates for each biomolecule. Conversely, biomolecule concentrations could be higher than our calculated values due to various concentration mechanisms including adsorption to mineral surfaces, and sequestration into mineral gels or amphiphilic matrices (Damer & Deamer 2020; Deamer 2017). Given these uncertainties, we present the following biomolecule concentrations as an upper bound in the absence of concentration mechanisms.

Adenine concentrations from aqueous production peak at 7.3 nM (1 ppb) and 0.2 pM (0.03 pptr) for our early Hadean (reducing) and late Hadean (oxidizing) models, respectively. The early Hadean (reducing) adenine concentrations are approximately three orders of magnitude smaller than the peak adenine concentration from meteoritic delivery of 10.6 μM (1.43 ppm); however, the adenine concentration from aqueous production can be sustained for more than 100 million years rather than days.

Adenine concentrations from delivery by IDPs are the most dilute in our WLP models, peaking at ∼10⁻¹⁰ μM (Pearce et al.2017).

In Figure 4(B), we plot the pond concentrations of HCN and H₂CO from atmospheric rain-out, as well as the concentrations of nucleobases, ribose, H₂CO, and 2-aminooxazole from aqueous HCN-based production. HCN concentrations peak at 0.04 μM and reduce to approximately 3 nM when the water level in the pond is highest.

We learn from these models that formaldehyde in WLPs likely did not come directly from the atmosphere during the early Hadean. H₂CO concentrations from aqueous photolytic production peak at 1.5 nM, which is 3–4 orders of magnitude higher than the maximum H₂CO concentration from atmospheric rain-out (4.1 × 10⁻⁵ nM). On the other hand, for our late Hadean (oxidizing) model B, the H₂CO concentration from rain-out is 0.34 nM, which is two orders of magnitude higher than the H₂CO concentration from photolytic production. Therefore, in oxidizing conditions H₂CO in WLPs may also come directly from the atmosphere.

Our model solves another main limitation of meteorites as a source of prebiotic nucleobases, in that cytosine and thymine are not present in meteorites (Callahan et al. 2011; Pearce & Pudritz 2015, 2016). During the early Hadean, guanine, cytosine, uracil, and thymine concentrations peak at 8.2, 1.4, 0.7, and 0.5 nM, respectively. 2-aminooxazole and ribose concentrations peak at 0.045 and 0.015 nM, respectively. Meteorites would have provided brief (∼a few weeks) enhancements in nucleobase concentrations up to three orders of magnitude above the base abundance from in situ production.

In Figure 5 and Table B1, we summarize the peak concentrations of key biomolecules and biomolecule precursors in our model WLPs for our early Hadean (reducing) and late Hadean (oxidizing) models A and B. We also include the early Hadean (reducing) model A with seepage turned off, representing a scenario where biomolecules are not exposed to this sink due to blockage of rock pores by amphiphilic multilamellar matrices or mineral gels, or adsorption onto mineral surfaces. Concentrations of biomolecules from the early Hadean (reducing) model A are approximately four orders of magnitude higher than concentrations from the late Hadean (oxidizing) model B. When seepage is turned off, biomolecule concentrations increase in model A to ∼1–400 nM. It is interesting to note that there are order-of-magnitude differences in the concentrations of purine nucleobases (adenine and guanine) versus the pyrimidine nucleobases (cytosine, uracil, and thymine). Does this have implications for the frequencies of various bases in the first replicating RNA polymers?

Zoom In Zoom Out Reset image size

We have developed multiple new methods that have greatly enhanced the capabilities of nonequilibrium calculations of atmospheric HCN on early Earth. These include the calculation of composition-dependent pressure–temperature (P–T) profiles using a radiative transfer code, the inclusion of lightning chemistry, and the time-dependent influx of H₂, CO₂, and CH₄ from impact degassing, volcanism, and oceanic geochemistry, respectively.

Past nonequilibrium atmospheric models for the Archean (∼3.8–2.5 bya) have computed HCN production for a range of CH₄, H₂, and CO₂ abundances using a commonly suggested input P–T profile for early Earth with surface temperatures of 273–288 K (Zahnle 1986; Tian et al. 2011). These models imposed CH₄ abundances in the 10–1000 ppm range, either did not fix H₂ or imposed H₂ abundances in the 0.01%–1% range, and imposed CO₂ abundances in the 0.04%–3% range. Our strategy is different, as we begin with initial reducing (4.4 bya) or oxidizing (4.0 bya) conditions that are thought to represent each epoch, and calculate the composition-dependent input P–T profiles using a radiative transfer code. As a result, our surface temperatures are 27°C–78°C hotter than these past Archean models. Then, we allow the concentrations of all species to evolve over time based on their source and sink rates at each epoch.

For example, CO₂ reaches a steady surface abundance of ∼0.05%–0.1% in our models based on the source rates from volcanic outgassing and association chemistry balanced with the sink rates from atmospheric rain-out and photodissociation. Our steady-state abundances are similar to the lower end of CO₂ abundances from past Archean models. Similarly, our end-of-simulation H₂ abundances (18%–19% and ppb–ppm for reducing and oxidizing models, respectively) are balanced by the source rate from epoch-dependent impact bombardment, and the sink rates from hydrodynamic escape, photodissociation, and chemistry. The steady-state H₂ abundances in our reducing models are at least an order of magnitude higher than those calculated or imposed in the Archean models, while the steady-state H₂ abundances in our oxidizing models are at least two orders of magnitude lower than the values calculated or imposed in those models.

CH₄ abundances are based on source rates from oceanic sources balanced with sink rates from photodissociation and chemistry. Surface CH₄ abundances reach end-of-simulation values that are 0–4 orders of magnitude lower than the range of CH₄ concentrations imposed in the past Archean models. Given the correlation between HCN production and CH₄ abundance, our calculated surface HCN concentrations tend to be lower than those calculated in these Archean models. The Archean models compute HCN mixing ratios near the surface to be ∼10⁻¹²–10⁻⁷. The lower end of this range is similar to the surface HCN abundances in our reducing 4.4 bya models of ∼10⁻¹¹, but three orders of magnitude higher than the HCN mixing ratios in our oxidizing 4.0 bya models of ∼10⁻¹⁵ after ∼1 Myr. This three-orders-of-magnitude discrepancy is due to the fact that our 4.0 bya models are more oxidizing than the previous Archean models. In other words, our 4.0 bya models contain higher concentrations of oxygen species such as O₂ and H₂O, and lower concentrations of reducing gases such as H₂ and CH₄.

Zahnle et al. (2020) modeled nonequilibrium chemistry following a large-body impact with input P–T profiles based on a simple analytic equation for a moist adiabat. These models began with concentrations of H₂, CO, CO₂, CH₄, and NH₃ that result from various impactors equilibrating with different mineral redox buffers. They found that impactors of at least the size of Vesta (525 km) are required to sustain the high temperatures required for rapid methane production. These models produced post-impact atmospheric CH₄ abundances of ∼3%, which resulted in similarly high (∼a few %) HCN abundances for a few million years after impact. These post-impact CH₄ abundances are four orders of magnitude higher than the ppm-range steady-state CH₄ abundances from our reducing models. The analytic equation Zahnle et al. (2020) used for obtaining a habitable surface and P–T profile for their nonequilibrium chemistry models does not consider the strength of various opacity sources such as H₂–H₂ collisionally induced absorption (CIA) and CH₄. Zahnle et al. (2020) compute surface temperatures of ∼320 K for a post-Vesta-impact atmosphere of 3.9 bar of H₂ and 0.17 bar of CH₄. We are unable to obtain habitable surfaces when modeling H₂ atmospheres >2 bar using the equations of radiative transfer. In our radiative transfer models, we find H₂–H₂ CIA produces a strong greenhouse effect above ∼1.13 bar of H₂. The surface temperature of our early Hadean (reducing) model A, which has 1.5 bar of H₂, 2 ppm of CH₄, and ppm-range H₂O, is reaching the limits of habitability at 78°C. This suggests that the resultant high atmospheric pressures of H₂ and CH₄ from the large-body impacts modeled by Zahnle et al. (2020) would produce atmospheric temperatures too high for WLPs to exist. More research needs to be done to understand whether a large-body impact could provide high HCN rain-out rates to WLPs after a substantial amount of H₂ escapes from the upper atmosphere, allowing the surface to cool to habitable temperatures.

Atmospheric models of nitrogen-rich rocky exoplanets that use C/O ratio as an adjustable parameter produce HCN mixing ratios of 10⁻⁸–10⁻⁷ for atmospheric C/O ratios near 0.5, and HCN mixing ratios of ∼10⁻³ for C/O ratios >1.5 (Rimmer & Rugheimer 2019). We do not use C/O as an adjustable parameter, because we find that the balance of outgassing and losses of species such as CO₂, H₂O, and CH₄ in our models leads to surface C/O ratios that vary from ∼0.001 to 1 over the course of the simulations. These C/O ratios are generally lower than those explored by Rimmer & Rugheimer (2019).

Pinto et al. (1980) modeled the chemical kinetics of formaldehyde production in a primitive Earth atmosphere composed of N₂, H₂O, CO₂, CO, and H₂. They employed a highly reduced network of eight photochemical reactions and 31 neutral reactions between 13 species. Their computed formaldehyde rain-out rates (∼10¹¹ mol yr⁻¹) are approximately three orders of magnitude higher than our maximum formaldehyde rain-out rates for Model B (∼10⁸ mol yr⁻¹). Part of this discrepancy is due to the higher efficiency of their H₂CO rain-out. They used a rain-out (scavenging) coefficient (cm³ s⁻¹) that is based on H₂CO rain-out on Earth today, whereas we compute rain-out using a deposition velocity (cm s⁻¹) that is calculated from a thin-film parameterization model (Wagner et al. 2002) and is used in more recent atmospheric models that follow the treatment of deposition velocity (Ranjan et al. 2020; Hu et al. 2012). The treatment used in Pinto et al. (1980) produces an H₂CO rain-out rate that is about a factor of 20 greater than our models would produce if we had the same surface H₂CO density. Considering these differences, we may be underestimating H₂CO rain-out by up to a factor of ∼20. However, the main source of discrepancy is because their atmosphere is more oxidizing than all of our models. They do not include CH₄ or other reduced carbon species (e.g., C₂H₂, C₂H₆) in their network, which, if given initial abundances, would change the oxidation state of the atmosphere and decrease the production of oxidized carbon species such as H₂CO.

Miller–Urey experiments have shown that reducing conditions are more favorable than oxidizing conditions for biomolecule production (Schlesinger & Miller 1983; Cleaves et al. 2008). For example, Schlesinger & Miller (1983) found a difference of ∼3–4 orders of magnitude in amino acid yields when switching from reducing (H₂-dominant) to oxidizing (CO₂-dominant) experimental conditions. This is consistent with our results, where we have shown that the early reducing phase of the Hadean eon at 4.4 bya produces atmospheric HCN and RNA building blocks ∼4 orders of magnitude higher in concentration than during the late oxidizing phase at 4.0 bya. Some Miller–Urey experiments have shown reasonable success in amino acid production in CO₂/N₂ conditions when buffering the solution with calcium carbonate (Cleaves et al. 2008); however, this is likely demonstrating catalytic effects that increase amino acid production yields from low HCN concentrations.

Of critical importance is how high pond concentrations of biomolecules such as nucleobases, ribose, and 2-aminooxazole would need to be in order for nucleotide synthesis and subsequent RNA polymerization to occur. Published laboratory experiments that react nucleobases to produce nucleosides and nucleotides typically use 400 μM–100 mM concentrations of nucleobases and obtain 0.01%–74% yields (Ponnamperuma et al. 1963; Fuller et al. 1972; Saladino et al. 2017; Nam et al. 2018). The lowest end of this experimental concentration range is three orders of magnitude higher than the maximum nucleobase concentrations from our early Hadean (4.4 bya) model without seepage (∼400 nM). Mechanisms to increase nucleobase concentrations in WLPs would likely be necessary given that experiments at nanomolar concentrations to produce nucleotides and RNA have not been reported, and achieving high yields requires optimal conditions. Such mechanisms include adsorption to mineral surfaces, and sequestration into amphiphilic multilamellar matrices or mineral gels (Damer & Deamer 2020; Deamer 2017).

The porosity of WLPs is difficult to assess. While something is known about the porosity of volcanic basalts (Hamouda et al.2014), pond environments are likely to be messier affairs. Various kinds of debris and films would be expected to accumulate at the bottom of ponds, such as amphiphilic multilamellar matrices or mineral gels. These would clog the pores, reducing losses due to seepage. It is difficult to estimate how large this effect would be. We can, however, estimate an upper limit to abundances of biomolecules by computing models for which seepage is set to zero. Turning off seepage in our models provides an increase in nucleobase concentrations by 1–3 orders of magnitude after 500–10,000 yr (see Figure B3). In this scenario, hydrolysis takes over as the the main concentration-limiting sink.

Ribose and 2-aminooxazole in our models are the most dilute of the RNA building blocks, with concentrations in the sub-nanomolar range. Again, laboratory experiments use much higher concentrations than this for nucleotide synthesis, typically in the millimolar to molar range (Ponnamperuma et al. 1963; Fuller et al. 1972; Nam et al. 2018; Powner et al. 2009). Saladino et al. (2017) used only 8 μM ribose to produce nucleosides in a neat formamide solution; however, because this experiment was performed with a different solvent, it is unclear whether similar results could occur in aqueous solution. Our results encourage additional experimental work to determine the levels that are sufficient for nucleotide synthesis in realistic prebiotic conditions.

Self-consistent disequilibrium atmospheric simulations are carried out iteratively using the consistent reduced atmospheric hybrid chemical network oxygen extension (CRAHCN-O) (Pearce et al. 2020a, 2020b) coupled with a 1D chemical kinetic model (ChemKM) (Molaverdikhani et al. 2019a, 2020a), with input pressure–temperature (P–T) structures calculated using petitRADTRANS (Mollière et al. 2019).

CRAHCN-O now contains 259 one-, two-, and three-body reactions, whose rate coefficients are gathered from experiments when available (∼40%), and are otherwise calculated using accurate, consistent, theoretical quantum methods (∼60%). Approximately 93 of the reactions in CRAHCN-O were missing from the literature prior to their discovery in Pearce et al. (2019, 2020a, 2020b). This network can be used to calculate HCN and H₂CO chemistry in atmospheres characterized by any of N₂, CO₂, CH₄, H₂O, and H₂. We added 28 new reactions to CRAHCN-O after our preliminary simulations showed the artificial buildup of species that previously had no reaction sinks. We list the new two-body and three-body reactions in Tables A1 and A2, respectively.

Table A1. New Two-body Reactions Added to CRAHCN-O for Our Early Earth Atmospheric Models, and Their Experimental Arrhenius Coefficients

Reaction Equation	α	β	γ	Source(s)
HCCO + NO → HCNO + CO	1.4 × 10−11	0	−320	Carl et al. (2000)
$\mathrm{HCCO}+\mathrm{NO}\to \mathrm{HCN}+{\mathrm{CO}}_{2}$	6.1 × 10−12	−0.72	−200	Carl et al. (2000)
HCCO + 3O → CO + CO + H	1.6 × 10−10	0	0	Baulch et al. (1992)
HCCO + H → CO + 3CH2	2.1 × 10−10	0	0	Glass et al. (2000); Frank et al. (1988)
NCO + O2 → CO2 + NO	1.3 × 10−12	0	0	Schacke et al. (1974)
NCO + NO → N2 + CO2	1.6 × 10−11	0	0	Cooper et al. (1993); Cooper & Hershberger (1992)
NCO + 3O → NO + CO	6.4 × 10−11	−1.14	0	Becker et al. (2000)
NCO + H → NH + CO	2.2 × 10−11	0	0	Becker et al. (2000)
HO2 + 3O → O2 + OH	5.4 × 10−11	0	0	Baulch et al. (1992)
HO2 + OH → H2O + O2	4.8 × 10−11	0	−250	Baulch et al. (1992)
HO2 + H → H2O + 3O	5.0 × 10−11	0	866	Baulch et al. (1992)
HO2 + H → O2 + H2	7.1 × 10−11	0	710	Baulch et al. (1992)
HO2 + H → OH + OH	2.8 × 10−10	0	440	Baulch et al. (1992)
O2 + HCO → CO + HO2	8.5 × 10−11	0	850	Tsang & Hampson (1986)
O2 + C2H → HCCO + 3O	1.0 × 10−12	0	0	Tsang & Hampson (1986)
O2 + C2H → HCO + CO	4.0 × 10−12	0	0	Tsang & Hampson (1986)
${{\rm{O}}}_{2}+\mathrm{CN}\to \mathrm{NCO}+{}^{3}{\rm{O}}$	1.1 × 10−11	0	−205	Baulch et al. (1992)
O2 + 3N → NO + 3O	4.5 × 10−12	1.0	3720	Baulch et al. (1994)
O2 + CH → OH + CO	5.0 × 10−11	0	0	Lichtin et al. (1984, 1983)
				Duncanson & Guillory (1983); Messing et al. (1979)
O2 + C → CO + 3O	3.0 × 10−11	0	0	Geppert et al. (2000); Dorthe et al. (1991)
				Becker et al. (1988); Husain & Young (1975)
				Husain & Kirsch (1971); Braun et al. (1969)
				Martinotti et al. (1968)
NO + 4N → N2 + 3O	3.1 × 10−11	0	0	Atkinson et al. (1989)
NO + 2N → N2 + 3O	6.0 × 10−11	0	0	Herron (1999)
$\mathrm{NO}+{\rm{C}}\to \mathrm{CN}+{}^{3}{\rm{O}}$	2.5 × 10−11	0	0	Baulch et al. (1994)
C2H + CH4 → C2H2 + CH3	3.0 × 10−12	0	250	Tsang & Hampson (1986)
C2H + 3O → OH + CO	1.7 × 10−11	0	0	Baulch et al. (1992)

Note. These are the most efficient sink reactions for species that would otherwise erroneously build up over timescales of tens to hundreds of millions of years. The Arrhenius expression is $k(T)=\alpha {\left(T/300\,{\rm{K}}\right)}^{\beta }{e}^{-\gamma /T}$.

Download table as:  ASCIITypeset image

All values are experimental, except for the low- and high-pressure limit rate coefficients for C₂H + H + M → C₂H₂ + M and the high-pressure limit rate coefficient for C + H₂ +M → ³CH₂ + M, as there were no experimental values available. For these reactions, we calculate the rate coefficients using the same validated theoretical and computational quantum methods developed in Pearce et al. (2020a, 2020b) for the other three-body reactions in CRAHCN-O.

The ChemKM code takes as input: (a) a chemical network, (b) an atmospheric P–T structure, (c) an eddy diffusion profile to characterize turbulent mixing, (d) the solar radiation spectrum at the top-of-atmosphere (TOA), (e) wavelength-dependent photochemical reactions, and (f) incoming and outgoing molecular fluxes from the surface and TOA (i.e., impact degassing, volcanic outgassing, rain-out, hydrogen escape, and chemical production from lightning). ChemKM uses the plane-parallel two-stream approximation to calculate radiative transfer, and includes both photoabsorption and Rayleigh scattering. The pressure profiles remain static throughout the simulations; therefore we must assume that pressure has reached equilibrium with the influx and outflux of atmospheric gases. This assumption is valid for our models, as our H₂ impact degassing rates never exceed the H₂ escape rates.

P–T structures for early Earth models are calculated using the petitRADTRANS software package (Mollière et al. 2019). petitRADTRANS is a 1D radiative transfer code that uses the Guillot analytic temperature model (Guillot 2010) and correlated-k opacity tables to solve for atmospheric temperatures of planets with no surface boundary condition. Visible opacities are calculated using the Planck mean, and infrared opacities are calculated using the Rosseland mean (Parmentier & Guillot 2014). Our models are all 100 layers, from surface pressures of ∼1–2 bar to TOA pressures of 1 × 10⁻⁸ bar. We implement cloudless and hazeless models given the large uncertainties of these parameters for the early atmosphere. We note that the lack of biogenic cloud condensation nuclei would have led to a shorter lifetime for optically thick convective clouds during the Hadean (Rosing et al. 2010), reducing their contribution to both the greenhouse and albedo when compared with those of modern Earth (Charnay et al. 2020).

In Figure A1 we display the initial water vapor profiles in the troposphere of our four early Earth models. We make an incremental improvement over the standard water vapor profile of Manabe & Wetherald (1967) for Earth's atmosphere by calculating tropospheric water abundances in two steps. In the first step, we iterate the Arden Buck equation (Buck 1981), which is dependent on temperature, and the P–T calculations from petitRADTRANS, which are dependent on water composition. To avoid a runaway greenhouse due to water vapor feedback, we parameterize the strength of water vapor feedback by decreasing the relative humidity (RH) by β = 6% for every 1°C of warming (Held & Soden 2000). In step two, we smooth out these profiles to avoid numerical instabilities in ChemKM. We use the calculated water vapor at the tropopause and linearly increase the mixing ratio moving downwards in altitude to reach a typical water vapor abundance of 1% for wet rocky planets (Hu et al. 2012). The tropopause is chosen to be 0.14 bar, similar to the present-day Earth (Robinson & Catling 2014).

Zoom In Zoom Out Reset image size

The three main inputs for our P–T structure calculations are (a) input composition, (b) equilibrium temperature (T_eq), and (c) internal temperature (T_int). Initial guesses for input compositions were selected to represent reducing (H₂/N₂-dominant) or oxidizing (CO₂/N₂-dominant) phases of the early (4.4 bya) and late (4.0 bya) Hadean eon. Surface pressure and methane abundance are adjusted from 1 to 2 bar and 1 to 10 ppm, respectively, to maintain habitable surface temperatures (i.e., 0 °C ≤ T_s ≤ 100 °C) (see Table 1). Equilibrium temperatures are calculated using the equation

$$\begin{eqnarray}&&{T}_{\mathrm{eq}}={\left(\displaystyle \frac{(1-A){L}_{\odot }}{16\pi \sigma {a}^{2}}\right)}^{1/4},\end{eqnarray} \tag{ A1 }$$

where T_eq is equilibrium temperature, A is albedo, L_⊙ is solar luminosity, σ is the Stefan–Boltzmann constant, and a is the semimajor axis of the planet.

Luminosities for the Sun at 4.4 bya (0.705 L_⊙) and 4.0 bya (0.728 L_⊙) are obtained from a precomputed stellar evolution model of a Sun-like star (Heller et al. 2020; Baraffe et al. 2015). The Hadean Earth would have been mostly covered in water; therefore, albedo is taken to be 0.06, which is consistent with a cloudless water world (Roesch et al. 2002).

Internal heat flow is taken to be three times the present value, which is compatible with thermal modeling of the Hadean (Sleep 2010). Using the Stefan–Boltzmann law, this results in internal temperatures of T_int = 43.3 K.

In Figure A2, we display the P–T profiles for our four Hadean models.

Zoom In Zoom Out Reset image size

In Figure A3, we display the eddy diffusion profile used for all our early Earth models. This is the standard profile for early Earth and analogous exoplanets (Ranjan et al. 2020; Arney et al. 2016; Tian et al. 2011; Zahnle 1986).

Zoom In Zoom Out Reset image size

The TOA radiation for our models is based on the solar mean (Thuillier et al. 2004) with a solar zenith angle of 50°; however, we increase the UV flux to simulate the increased activity of the young Sun. In Table A3, we display the multiplicative factors used in our models for each UV wavelength interval. Values are based on observations of young solar analogs (ages ∼0.1–7 Gyr) (Ribas et al. 2005).

Our 33 photochemical reactions mostly match those from the Titan models by Hébrard et al. (2012); however, we update the H₂O absorption cross sections with the recent near-UV experimental measurements from Ranjan et al. (2020), we remove erroneous CO₂ absorption below 202 nm (Ranjan et al. 2020), and we add photochemistry for O₂ and HO₂ following treatments in Hu et al. (2012) and CH₃OH using experimental cross sections from Lange et al. (2020) and Burton et al. (1992).

We couple planetary surface processes to our atmospheric models by adding influxes of species to the lowest layer of our atmospheres. These include H₂ impact degassing, CO₂ outgassing from volcanoes, CH₄ outgassing from hydrothermal systems, H₂O evaporation from the ocean, and chemical production (e.g., HCN, CO, ³O, H) due to lightning.

Equilibrium chemistry calculations performed by Zahnle et al. (2020) for enstatite chondrite impactors suggest that H₂ degassing via the reaction Fe + H₂O → FeO + H₂ scales linearly with impactor mass, at a rate of ∼10⁻²¹ mol H₂ cm⁻² g⁻¹ impactor. Mathematical fits to the lunar cratering record provide us with an estimate of the rate of impactors on early Earth at a given epoch. In Figure A4, we display our model bombardment rates, which lie between the minimum and maximum bombardment fits (Pearce et al. 2017; Chyba 1990). These bombardment rates are 1.2 × 10²⁵ g Gyr⁻¹ and 1.2 × 10²⁴ g Gyr⁻¹ at 4.4 bya and 4.0 bya, respectively. Multiplying the H₂ degassing abundance per unit mass by the mass delivery rates at 4.4 bya and 4.0 bya gives us H₂ impact degassing rates of 2.3 × 10¹¹ cm⁻² s⁻¹ and 2.3 × 10¹⁰ cm⁻² s⁻¹ for the early and late Hadean models, respectively.

Zoom In Zoom Out Reset image size

We model H and H₂ escape using the approximation developed by Zahnle et al. (2019, 2020), which blends energy-limited and diffusion-limited escape. The equation is

$$\begin{eqnarray}&&{\left(\displaystyle \frac{{{dN}}_{{{\rm{H}}}_{2}}}{{dt}}\right)}_{{\rm{esc}}}=-\displaystyle \frac{{AS}}{\sqrt{1+{B}^{2}{S}^{2}}}\displaystyle \frac{{N}_{{{\rm{H}}}_{2}}}{{{\rm{\Sigma }}}_{j}{N}_{j}}\,({{\rm{cm}}}^{-2}\,{{\rm{s}}}^{-1})\end{eqnarray} \tag{ A2 }$$

where A = 2 × 10¹² cm⁻² s⁻¹, B² = 0.006, S is the extreme-UV and far-UV irradiation relative to the modern Sun (i.e., 30 and 9 for 4.4 and 4.0 bya, respectively), ${N}_{{{\rm{H}}}_{2}}$ is the number of H₂ molecules, and N_j is the number of molecules of species j.

The limits of H₂ degassing to assume a static pressure profile in equilibrium with H₂ escape would be ∼2.3 × 10¹³ and 1.4 × 10¹³ cm⁻² s⁻¹ at 4.4 and 4.0 bya, respectively.

We use a constant CO₂ outgassing rate of 3.0 × 10¹¹ cm⁻² s⁻¹ in all our models, which is consistent with Earth-like volcanic outgassing used in other atmospheric models (Hu et al. 2012).

Guzmán-Marmolejo et al. (2013) modeled the production of CH₄ in hydrothermal systems for an Earth-like planet, and calculated a production rate of 6.8 × 10⁸ cm⁻² s⁻¹. Guzmán-Marmolejo et al. (2013) also modeled H₂ production in hydrothermal systems; however, the rates of H₂ produced in hydrothermal environments are orders of magnitude lower than the rates of H₂ production from impact degassing.

We do not include loss of CH₄ due to haze production, as our CH₄/CO₂ ratio never exceeds 0.1 (which is a common identifier for haze production) (Trainer et al. 2006).

In Table A4, we list the deposition velocities for the species that are rained out of the lowest layer in our atmospheric models.

Table A4. Deposition Velocities for the Chemical Species in Our Early Earth Models

Species	Deposition (cm s−1)	Source
CO2	1 × 10−4	Archer (2010)
CH3OH	0.1	Wohlfahrt et al. (2015)
O2	1 × 10−8	This work
CO	1 × 10−8	Kharecha et al. (2005)
H2O2	0.5	Hauglustaine et al. (1994)
C2H6	1 × 10−5	Hu et al. (2012)
HO2	1	Ranjan et al. (2020)
H2CO	0.1	Wagner et al. (2002)
HCO	0.1	Ranjan et al. (2020)
$\mathrm{HCN}$	7 × 10−3	Tian et al. (2011)
OH	1	Ranjan et al. (2020)
3O	1	Ranjan et al. (2020)
H	1	Ranjan et al. (2020)

Download table as:  ASCIITypeset image

We apply a CO₂ deposition velocity of 1.0 × 10⁻⁴ cm s⁻¹, which is estimated in Hu et al. (2012) to produce a CO₂ lifetime consistent with the lifetime of silicate weathering on Earth. We use a CO deposition velocity of 1 × 10⁻⁸ cm s⁻¹, calculated from a two-box model (Kharecha et al. 2005). We also use a deposition velocity of 1 × 10⁻⁸ cm s⁻¹ for O₂, given its similar solubility and diffusivity to CO (Harman et al. 2015). We use the standard HCN deposition velocity of 7.0 × 10⁻³ cm s⁻¹ that is used in other early Earth models (Tian et al. 2011; Zahnle 1986). Additional deposition velocities are chosen to be consistent with other atmospheric models for rocky exoplanets (Hu et al. 2012; Ranjan et al. 2020). Major species not listed in this table (e.g., H₂ and CH₄) are not very soluble in water, therefore we do not include rain-out for these species (Sander 2020).

A.1.1. Lightning

Lightning chemistry in the context of the origin of life was first developed experimentally in the 1950s (Miller 1953). The fundamental Miller–Urey experiment involves sending an electric discharge through a combination of reduced gases to trigger dissociation. The radicals produced in this process then react to form biomolecule precursors such as HCN and H₂CO (Miller 1957b; Bada 2016). These precursors condense into a reservoir, where aqueous chemistry produces biomolecules such as amino acids (Miller 1953; Bada 2013) and nucleobases (Ferus et al. 2017b).

Present-day Earth has an average global lightning flash density of ∼2 flashes km⁻² yr⁻¹ (Hodosán et al. 2016). However, above just the oceans, this average density drops to 0.3–0.6 flashes km⁻² yr⁻¹. Given the smaller coverage of continental crust above sea water during the Hadean, we set the global lightning flash density for our models to 1 flash km⁻² yr⁻¹; however, we also explore an average lightning flash density measured locally during volcanic eruptions on Earth today (∼10⁴ flashes km⁻² yr⁻¹).

We considered both nonequilibrium and equilibrium approaches for modeling lightning chemistry. For our nonequilibrium approach, we integrated the production of key radicals for the first ∼40 μs of a lightning strike using the pressure and temperature evolution from Ardaseva et al. (2017). However, this approach had accuracy issues as a complete high-temperature reaction network is required to accurately calculate the chemical evolution within a cooling lightning channel. This approach also led to some nonsensical results such as HCN production that is independent of lightning flash density.

Therefore, we use an equilibrium approach for modeling the lightning production of HCN and other species based on the lightning chemistry models for HCN and NO production by Chameides & Walker (1981). Lightning channels heat up to a point (∼30,000 K) (Orville 1968) where the equilibrium timescale is less than 1 μs (Hill et al. 1980). In fact, ab initio molecular dynamics simulations of electric discharges suggest the timescale for lightning chemistry may be of the order of picoseconds (Cassone et al. 2018). This is fast compared to the hundred millisecond cooling timescale of a lightning channel, as well as the 10 μs eddy diffusion timescale in the lowest layer in our atmospheres.

Chemical abundances rapidly reach equilibrium while the lightning channel is above the freeze-out temperature (T_F ). The freeze-out temperature is the temperature at which the concentration of a species can still be described by its equilibrium value. Beyond this point, there is not enough time at a given temperature for equilibrium to be reached, and thus the concentrations are frozen into the gas for the remainder of the cooling of the lightning channel. Reaction rate coefficients that break down a species are used to roughly determine the freeze-out temperature, e.g., HCN + MCN + H + M. Typical freeze-out temperatures range from 1000 to 5000 K. The freeze-out temperature for HCN is ∼2000–2500 K for lightning strikes similar in energy to Earth today (10⁵ J m⁻¹) (Chameides & Walker 1981). Other species such as NO have higher freeze-out temperatures near 3000–3500 K for similar lightning discharge energies, but can also be ∼2000 K for the highest discharge energies (10¹⁵ J m⁻¹). We adopt T_F = 2000 K for our equilibrium calculations to estimate the mixing ratios for HCN and 13 other dominant equilibrium products in the early Earth lightning models by Chameides & Walker (1981).

Equilibrium calculations are performed using the thermochemical data from the JANAF tables (Stull & Prophet 1971), and the ChemApp software library (distributed by GTT Technologies, http://gtt.mch.rwth-aachen.de/gtt-web/).

In Table A5, we display the equilibrium mixing ratios based on the initial abundances in our four early Earth models for a freeze-out temperature of T_F = 2000 K.

Table A5. Equilibrium Abundances (Molar Mixing Ratios) from Lightning Chemistry Occurring in Our Four Early Earth Models

Species	Model A	Model B	Model C	Model D
$\mathrm{HCN}$	1.4 × 10−10	1.8 × 10−13	5.1 × 10−11	1.8 × 10−13
H2	8.9 × 10−1	8.7 × 10−1	8.9 × 10−1	8.7 × 10−1
H	1.3 × 10−3	6.0 × 10−6	1.4 × 10−3	6.0 × 10−6
N2	9.9 × 10−2	9.8 × 10−2	9.9 × 10−2	9.8 × 10−2
4N	2.3 × 10−10	2.0 × 10−10	2.7 × 10−10	2.0 × 10−10
H2O	9.9 × 10−3	9.7 × 10−3	9.9 × 10−3	9.7 × 10−3
3O	1.4 × 10−9	3.4 × 10−5	1.9 × 10−9	3.4 × 10−5
OH	1.4 × 10−6	2.2 × 10−4	1.6 × 10−6	2.2 × 10−4
NO	1.6 × 10−8	4.6 × 10−4	1.9 × 10−8	4.6 × 10−4
O2	6.9 × 10−12	5.3 × 10−3	9.1 × 10−12	5.3 × 10−3
CH4	1.4 × 10−11	4.2 × 1021	4.0 × 10−12	4.2 × 10−21
CH3	5.0 × 10−13	2.3 × 10−20	1.6 × 10−13	2.0 × 10−20
CO	2.0 × 10−6	1.1 × 10−2	9.9 × 10−7	1.1 × 10−2
CO2	4.8 × 10−9	8.7 × 10−1	2.4 × 10−9	8.7 × 10−1

Note. Thermodynamic simulations are based on initial concentrations in Table 1 for a freeze-out temperature of T_F = 2000 K.

Download table as:  ASCIITypeset image

For our fiducial models, we introduce HCN and other species to the bottom layer of the atmosphere at a rate that corresponds to 1 flash km⁻² yr⁻¹. The species influx rates from lightning chemistry are calculated along a lightning channel extending through the first layer of the atmosphere using the following equation:

$$\begin{eqnarray}&&\displaystyle \frac{d[{\rm{M}}]}{{dA}\,{dt}}=\displaystyle \frac{{n}_{M}P{\rm{\Delta }}H\dot{f}{\sigma }_{l}}{{k}_{{\rm{B}}}T\gamma }{\left(\displaystyle \frac{1}{100}\right)}^{3}{\left(\displaystyle \frac{1}{10,000}\right)}^{2},\end{eqnarray} \tag{ A3 }$$

where $d[{\rm{M}}]/{dA}\,{dt}$ is the molar concentration of species M produced per cm² per second, n_M is the molar mixing ratio of species M produced in the lightning strike (cm⁻³), ΔH is the height of the lowest atmospheric layer (cm), $\dot{f}$ is the lightning flash density in flashes km⁻² yr⁻¹, σ_l is the cross section of the lightning channel (∼1 cm²), γ = 3600 × 24 × 365.25 s yr⁻¹, and the remainder is unit conversion.

In Table A6, we display the sources and sinks for nucleobases, ribose, and 2-aminooxazole in our warm little pond models. All biomolecule reaction yields are based on HCN. For cytosine, uracil, thymine, ribose, and 2-aminooxazole, which require formaldehyde as a reactant, we use a formaldehyde yield from HCN of 3.6%, which is three times the glyceronitrile yield from radiolytic aqueous HCN experiments performed by Yi et al. (2020). We assume, to first order, that the UV radiation incident at the pond surface allowed this aqueous photolytic reaction to proceed with similar yields to the Yi et al. (2020) radiolytic experiments.

We consider two yields for adenine production from HCN in our pond solutions based on experiments: a lower yield of 0.5% based on aqueous reactions of HCN (Oró & Kimball 1961), and an upper yield of 18% based on HCN reactions with more ideal conditions for forming adenine (e.g., solutions containing NH₃ and ammonium formate (Hill & Orgel 2002; Wakamatsu et al. 1966).

For guanine and uracil we consider lower yields of 0.0067% and 0.0017%, respectively, based on experiments of frozen ammonium cyanide solutions (Miyakawa et al. 2002). We use theoretical yields for the upper bounds. Guanine has a theoretical yield of 20% based on the theoretical HCN-based reaction equation $5{\rm{HCN}}+{{\rm{H}}}_{2}{\rm{O}}\to {\rm{guanine}}+{{\rm{H}}}_{2}$ (Larowe & Regnier 2008). Uracil has a theoretical yield of 50%; however, this is based on H₂CO as a limiting reagent ($2{\rm{HCN}}+2{{\rm{H}}}_{2}{\rm{CO}}\to {\rm{uracil}}+{{\rm{H}}}_{2}$ (Larowe & Regnier 2008). Experiments of the Kiliani–Fischer synthesis of glyceronitrile produce H₂CO as an intermediate from aqueous solutions of HCN (Yi et al. 2020). Yields of glyceronitirile production are 1.2%; however, theory suggests three intermediate HCN molecules are involved in this reaction. Considering this, we apply a yield of 3.6% for H₂CO production from HCN. For uracil, this results in an upper yield of 50% × 3.6% = 1.8%.

Given the lack of experiments producing cytosine and thymine from aqueous HCN, we only consider the theoretical upper yields for these species based on the reaction equations $3{\rm{HCN}}+{{\rm{H}}}_{2}{\rm{CO}}\to {\rm{cytosine}}$ and $2{\rm{HCN}}+3{{\rm{H}}}_{2}{\rm{CO}}\to {\rm{thymine}}+{{\rm{H}}}_{2}{\rm{O}}$ (Larowe & Regnier 2008). Again, since H₂CO is the limiting reagent for these theoretical reactions, we apply the yield of 3.6% for H₂CO production from HCN, resulting in yields of 3.6% and 1.2% for cytosine and thymine production from HCN, respectively.

For 2-aminooxazole, we consider a yield of 0.11% based on radiolytic experiments of aqueous solutions of HCN (Yi et al. 2020).

Finally, experiments of ribose synthesis from H₂CO have identified ribose as a product, but yields remain uncertain. Shapiro (1988) suggests 1% as an upper bound, therefore we consider this yield and also apply the yield of 3.6% for H₂CO production from HCN to obtain an overall yield of 0.036%.

The sink rates for biomolecule photodestruction, seepage, and hydrolysis are chosen to match those modeled in Pearce et al. (2017) so that we can directly compare our results with their calculations of biomolecule concentrations in WLPs from meteorites and interplanetary dust particles. Photodestruction and seepage rates are consistent among all biomolecules given the lack of experimental data for each individual biomolecule in our study. Comparison studies of biomolecule photodegradation show differences of about an order of magnitude in photostability between biomolecules (Poch et al. 2015; Todd et al. 2019). Therefore, we expect approximately an order of magnitude additional uncertainty in our biomolecule concentrations during the dry phase, when irradiation is turned on.

In Figure B1, we display the rain-out rates for $\mathrm{HCN}$, H₂CO, and CO₂ as a function of time. These water-soluble species, and a few others, are removed from the lowest layer of our atmospheric models at each time step. The HCN rain-out rate for the early Hadean (reducing) model A is 1.3 × 10⁻⁷ kg m⁻² yr⁻¹ at 50 million years. This is about five orders of magnitude higher than the HCN rain-out rates for the late Hadean (oxidizing) models of ∼10⁻¹² (same units). The HCN rain-out rate for model C is approximately a factor of 5 lower than that for model A at 10 million years. We use the 50 million year and 800,000 year HCN rain-out rates from models A and B, respectively, for influx into our subsequent warm little pond models.

Zoom In Zoom Out Reset image size

The H₂CO rain-out rate for the early Hadean (reducing) model A is 6.0 × 10⁻¹² kg m⁻² yr⁻¹ at 50 million years. This is ∼2 orders of magnitude lower than the H₂CO rain-out rates for the late Hadean (oxidizing) models of ∼10⁻⁹ kg m⁻² yr⁻¹ at 1 million years.

In Figure B2, we display the molar abundance ratio of HCN/CH₄ in the lowest atmospheric layer from t = 500 yr onwards. For the early Hadean (reducing) models A and C, the average HCN/CH₄ ratio during this period is ∼(1–4) × 10⁻⁷. For the late Hadean (oxidizing) models B and D, the average HCN/CH₄ ratios are (3–5) × 10⁻⁶. The HCN/CH₄ ratios are constant at 10⁻⁸ for models A and C from t = 500 to 1 million years, but vary by up to three orders of magnitude beyond this point. Multiple oxygen species begin to fluctuate in abundance after ∼1 million years, including H₂O, CO, and OH, which affects $\mathrm{HCN}$ production but not the CH₄ abundance. Similarly, fluctuations in species such as O₂ and H₂ in models B and D after ∼600,000 yr cause the HCN/CH₄ ratios to fluctuate beyond this point.

Zoom In Zoom Out Reset image size

In Figure B3, we turn off seepage in our model WLPs from Figure 4, and display the results from t = 0 to 10,000 yr. This model represents a scenario where the rock pores at the base of the WLP are blocked by, e.g., amphiphilic multilamellar matrices or mineral gels (Deamer 2017; Damer & Deamer 2020). In the absence of seepage, hydrolysis takes over as the main biomolecule sink. We see species with high hydrolysis rates such as cytosine, adenine, and guanine reach a steady state in just a few hundred years. Hydrolysis rates for uracil and thymine are lower, allowing them to build in concentration to near micromolar abundances in 10,000 yr. 2-aminooxazole and ribose are not given hydrolysis rates in our models, and thus their concentrations here are to be considered maxima in the absence of aqueous chemical sinks. We summarize the no-seepage model biomolecule concentrations at t = 10,000 yr in Table B1.

Zoom In Zoom Out Reset image size

In Table B1, we display the peak concentrations of biomolecules in our model WLPs for our early Hadean (reducing), late Hadean (oxidizing), and no-seepage models.

In Figure B4, we plot the molar HCN mixing ratio as a function of altitude for our early Hadean (reducing) model A using an increased lightning flash density 10⁴ times greater than our fiducial rate. This lightning flash density represents an average value measured near volcanic eruptions on Earth today (Hodosán et al. 2016).

Zoom In Zoom Out Reset image size

We see no changes to the HCN profile when using an enhanced lightning flash density. This suggests that volcanically active regions of early Earth would not have higher HCN abundances than the global average from photochemistry.