Habitability of Exoplanet Waterworlds

Simulations of rocky-planet formation yield many rocky small-radius (R_pl < 1.6 R_⊕) planets that have planet water mass fractions f_W = 10–1000  × f_W,⊕ (Raymond et al. 2004, 2007; Ciesla et al. 2015; Mulders et al. 2015; Zain et al. 2018). Earth's f_W (∼10⁻⁴) appears to be the result of chance (e.g., Raymond et al. 2007; Schönbächler et al. 2010; Lichtenberg et al. 2016; Morbidelli et al. 2016) (Section 5.2). Indeed, simulations suggest f_W ≈ f_W,⊕ may be uncommon on planets in general (e.g., Tian & Ida 2015).

Among previous studies of waterworlds (e.g., Kuchner 2003; Léger et al. 2004; Selsis et al. 2007; Fu et al. 2010; Abbot et al. 2012; Levi et al. 2017; Unterborn et al. 2018), the closest in intent to our own are those of Kitzmann et al. (2015) and Noack et al. (2016). Kitzmann et al. (2015) consider carbonate-system equilibria and find a "sweet-spot" in total planet C similar to our own, but they do not consider geological processes. Noack et al. (2016) emphasize the development of high-pressure H₂O-ice (HP-ice) layers that may isolate the liquid ocean from the nutrients supplied by silicate-rock leaching. Noack et al. also confirm the result of Kite et al. (2009) that deep oceans suppress C exchange between the convecting mantle and the water ocean. We conservatively do not count planets with HP-ice as examples of habitability (see Section 6.4 for a discussion). This can set a T_surf-dependent upper limit on the depth of a habitable ocean (provided T_surf ≲ 375 K; Figure 3 below), and we consider only <8 GPa oceans. We go beyond Noack et al. (2016) by considering the effect of C on climate, by using an N-body code to calculate how volatile delivery and giant impacts "load the dice" by regulating the fraction of planets that can have cycle-independent planetary habitability and—most importantly—by tracking T_surf and p_CO2 for 10 Ga on the habitable worlds we model.

Most of the physical and chemical processes we discuss have been investigated previously in an Earth context. The novel aspect of our paper is that we apply these ideas to exoplanet waterworlds (Section 2).

In a CO₂–H₂O atmosphere, high surface temperatures will cause high H₂O mixing ratios in the stratosphere: this is referred to as a moist greenhouse state. In the moist greenhouse, H₂O molecules are photolyzed by λ ≲ 240 nm photons, producing H₂. This H₂ is subsequently lost in a hydrodynamic outflow of the upper atmosphere, an outflow that is underpinned by absorption of ≲100 nm photons (Kasting 1988). Loss of H₂ from H₂O tends to dry out the planet, and so moist-greenhouse worlds are outside the "conservative HZ" limits of Kasting et al. (1993).

However, for an initial water endowment of >50 Earth oceans in the HZ of an FGK star, the moist greenhouse cannot cause complete ocean loss (Kasting 1988; Lammer et al. 2009; Luger & Barnes 2015; Zahnle & Catling 2017) (Figure 1). For example, integrating the XUV-flux estimate of Lammer et al. (2009) for 1 au and a G-type star gives a loss of no more than 20 Earth oceans over 4.5 Ga. This is small compared to the water available on waterworlds. Moreover, the XUV-limited water loss is an upper bound, as other processes likely limit the water loss to a lower value (e.g., Kulikov et al. 2006; Wordsworth & Pierrehumbert 2013b; Tian 2015; Owen & Alvarez 2016; Bourrier et al. 2017).

Zoom In Zoom Out Reset image size

In principle, oceans can dry out if ocean water reacts with rocks to form hydrated minerals (Mustard 2018). In practice, for waterworlds that have solid rock interiors, atmosphere-interior exchange of H₂O is not important for setting the surface water inventory. That is because solid Earth mantle rocks have an H-storage capacity of only 1–10 Earth oceans (e.g., Hirschmann 2016); this is a small fraction of the total water available to a waterworld (Tikoo & Elkins-Tanton 2017; for a contrary view, see Marty 2012). For such worlds, the ocean a planet has after cooling from the last giant impact is the ocean that planet will keep.

Planetary habitability is modulated by CO₂ (Figure 2). Very high partial pressures of atmospheric CO₂ (p_CO2) lead to temperatures too high for life. However, intermediate p_CO2 can stave off global surface ice cover (up to a low-insolation limit where CO₂ must condense) and thus extend habitability (Figure 2) (Kopparapu et al. 2013). p_CO2 would have varied widely if Earth had always lacked a negative feedback on p_CO2 (Kasting & Catling 2003). On Earth, >90% of C is stored in rocks, and C cycles between the atmosphere+ocean+biosphere and rocks every ∼300 kyr in the modern era (Knoll et al. 2012); this is geologically rapid. As a result of this rapid geochemical cycling, a small initial imbalance between the rate of C release from rocks (volcanic/metamorphic outgassing) and C uptake into rocks (weathering uptake) could lead to global surface ice cover (a snowball) or alternatively lead to a moist greenhouse. An initial moist greenhouse with XUV-limited ocean loss would lead to loss of Earth's ocean water (Figures 1 and 21). Shallow-ocean rocky planets are vulnerable, unless geochemical cycles provide a negative feedback (e.g., carbonate-silicate weathering feedback) that can maintain habitability.⁷

Zoom In Zoom Out Reset image size

For worlds with oceans a few times deeper than the Earth's ocean—enough to drown the land, but not as deep as the oceans we model in this paper—it is more difficult to equalize C uptake into rocks and C outgassing (Abbot et al. 2012; Foley 2015). Even though C uptake by weathering of seafloor rocks might make up some of this imbalance (Coogan et al. 2016; Coogan & Gillis 2018), it is less likely that a small initial imbalance between C release and C uptake could be restored via a negative feedback on worlds with oceans a few times deeper than the Earth's.

Although waterworlds are buffered against water loss (Section 2.1), the surface temperature of that water—icy, heat-sterilized, or somewhere in between?—depends on p_CO2 (Figure 2). p_CO2 on waterworlds is set by atmosphere–ocean partitioning of C. The ocean is a major reservoir of C that equilibrates with the atmosphere on <10⁸ yr timescales. For example, Earth's pre-industrial partition of C is 60 parts in the ocean for one part in air, with an atmosphere–ocean exchange time of ∼10³ yr.

Zoom In Zoom Out Reset image size

Therefore, p_CO2 depends on C abundance, planet water mass fraction (≡ dilution), and ocean pH. The last matters because of the strong effect of H⁺ activity on carbonate equilibria (Figure 4). pH ≈ −log₁₀[H⁺], where the square brackets here denote molarity. (Despite the similar notation, pH has nothing to do with the partial pressure of H₂.) p_CO2 in the atmosphere is in equilibrium with dissolved CO₂ in the ocean:

$$\begin{eqnarray}&&{{\rm{p}}}_{\mathrm{CO}2}=\displaystyle \frac{[{\mathrm{CO}}_{2(\mathrm{aq})}]}{{k}_{H}}\end{eqnarray} \tag{ 1 }$$

where the square brackets here denote concentration, and k_H is a solubility constant. However, p_CO2 can be much less than expected by dividing the total inorganic C (i.e., [C]) in the ocean by k_H. Specifically, if the pH is high enough to form HCO₃⁻ or CO${}_{3}^{2-}$ at the expense of CO_2(aq) (Figure 4),

$$\begin{eqnarray}&&({\mathrm{CO}}_{2}(\mathrm{aq})+{{\rm{H}}}_{2}{\rm{O}}\rightleftharpoons {{\rm{H}}}_{2}{\mathrm{CO}}_{3})\rightleftharpoons {{\rm{H}}}^{+}\\ &&\quad +\,{\mathrm{HCO}}_{3}^{2-}\rightleftharpoons 2{{\rm{H}}}^{+}+{\mathrm{CO}}_{3}^{2-}\end{eqnarray} \tag{ 2 }$$

then p_CO2 can be very low.

Zoom In Zoom Out Reset image size

In this paper, we will refer to the sum of the electrically neutral species $({\mathrm{CO}}_{2}(\mathrm{aq})$ and ${{\rm{H}}}_{2}{\mathrm{CO}}_{3})$ as CO_2T; in Earth seawater, the concentration of ${{\rm{H}}}_{2}{\mathrm{CO}}_{3}$ is ≲0.3% that of ${\mathrm{CO}}_{2}(\mathrm{aq})$.

When ocean pH is high (H⁺ concentration is low), then by Le Chatelier's principle carbonic acid gives up its H⁺. Thus, at high pH, C is hosted in the ocean as CO${}_{3}^{2-}$ and HCO₃⁻, ${\mathrm{CO}}_{2(\mathrm{aq})}$ is low, and by Equation (1), p_CO2 is low. In this case, the fraction of total C in the atmosphere is small (Figure 4). This describes the modern Earth (Zeebe & Wolf-Gladrow 2001; Ridgwell & Zeebe 2005; Butler 1982; Zeebe 2012). In other words, high pH effectively sequesters CO₂ from the atmosphere by driving the carbonate equilibria to the right. This raises the ratio $[{\mathrm{CO}}_{3}^{2-}$]/[${\mathrm{CO}}_{2(\mathrm{aq})}]$ (Equation (2)), and so decreases the p_CO2 in equilibrium with the ocean (Equation (1)). For a fixed total atmosphere+ocean C inventory, p_CO2 must fall. Conversely, when ocean pH is low (H⁺ concentration is high), then by Le Chatelier's principle, C in the ocean exists mostly as CO_2T, and the fraction of total C in the atmosphere is large.

What sets pH? Ocean pH rises when dissolved rock (e.g., Ca²⁺, Na⁺) is added to the ocean. That is because charge balance is maintained almost exactly in habitable oceans. Adding a mole of Na²⁺ to an ocean relieves one mole of H⁺ from their duty of maintaining charge balance, so they revert to H₂O. Loss of H⁺ raises pH, and sucks C out of the atmosphere (Equation (1), Figure 4). This effect of dissolved rock on p_CO2 does not require carbonate minerals to form. However, C can exist mostly as solid CaCO₃ if [Ca] is high (or even as solid Na₂CO₃ minerals, if [Na] is extremely high). The relationship between pH and p_CO2 and pH is an equilibrium, and so it equally true to think of p_CO2 as driving pH.

Assuming constant atmosphere+ocean C content, when T rises ocean pH will fall. This fall is due to changes in k_H, the analogous equilibrium constants in Equation (2), and increasing self-dissociation of water.

A general model for p_CO2(t) would be complicated (Holland 1984; Hayes & Waldbauer 2006; Krissansen-Totton et al. 2018). Complexity occurs because p_CO2 is coupled to pH, pH is affected by rock–water reactions, and the controls on rock–water reactions are complex and evolve with time. Fortunately, for waterworlds, the problem is simpler.

The waterworld climate evolution problem is simplified by five factors. Each factor tends to reduce the number of fluxes that we have to model. The cumulative effect of these five factors is to uncouple climate evolution from ocean–mantle geochemical cycling (Figure 6). The most important factor is that C exchange between the convecting mantle and the water ocean is limited on waterworlds—Factor #4. This is because seafloor pressure on waterworlds curtails adiabatic decompression melting, which is the dominant mechanism of volcanism on rocky planets (Kite et al. 2009; Figure 5; Appendix A). The five factors, which we collectively term the "waterworld approximation," are as follows.

• 1.  
  H₂ is negligible.
• 2.  
  Water loss to space is small.
• 3.  
  H₂O in the atmosphere+ocean greatly outweighs H₂O species in the silicate mantle.
• 4.  
  C exchange between the deep mantle and the water ocean shuts down within O(10⁸) yr after the last giant impact.
• 5.  
  No land.

We explain each factor below.

• 1.  
  H₂/He is negligible. H₂ blankets of the appropriate thickness to prolong habitability within the HZ of Kopparapu et al. (2013) require fine-tuning to form. They also swiftly escape to space (e.g., Wordsworth 2012; Owen & Mohanty 2016; Odert et al. 2018; but see Ramirez & Kaltenegger 2017). While H₂-rich, low-density planets likely exist in the HZ (Rogers 2015; van Eylen et al. 2018), we expect that they will have uninhabitably hot rock-volatile interfaces ("surfaces") (Rogers 2012) and we do not model them. We do not track H₂ outgassing by Fe+H₂O $\to \,$ FeO+H₂. In the context of modern planet formation models, this reaction requires mixing of Fe with the volatile envelope during collisions between planets and planetary embryos. However, planetary embryos contain metal cores of >500 km diameter, which merge quickly (Dahl & Stevenson 2010; Jacobson et al. 2017, but see also Genda et al. 2017). We assume that the mantle is sufficiently oxidized that any volcanically outgassed C is in the form of CO₂, and that H₂ outgassing is minor. Therefore, we consider an H₂-free, CO₂+H₂O(± N₂) atmosphere.
• 2.  
  Water loss to space is small, relative to the initial water complement, for waterworlds orbiting FGK stars. This is explained in Section 2.1.
• 3.  
  H₂O in the atmosphere+ocean outweighs H₂O in the silicate mantle. Upon magma-ocean crystallization shortly after the last giant impact, the silicate-magma ocean exsolves H₂O (Elkins-Tanton 2011). The H₂O concentration that is retained in the now-frozen magma (=silicate rock) is limited by the upper mantle silicate minerals' H storage capacity.⁸ This storage capacity is small relative to the total water on waterworlds (Hirschmann 2006; Cowan & Abbot 2014, but see Marty 2012 for a contrary view). As a result, the water in the ocean greatly exceeds the water in the silicate mantle, so any subsequent cycling (e.g., Schaefer & Sasselov 2015; Komacek & Abbot 2016; Korenaga et al. 2017) will have little (fractional) effect on ocean depth. Hypothetical rock-hydration feedbacks (e.g., Kasting & Holm 1992), if they are real, would also be limited by the small mantle storage capacity for H. Moreover, rock-hydration feedbacks do not reduce ocean depth on worlds that have oceans much deeper than on Earth (Korenaga et al. 2017). For seafloor pressure >10 GPa, the H₂O storage capacity of mantle rock just beneath the stagnant lid jumps to ∼0.5–1.0 wt% (Ohtani 2005). Even for such high seafloor pressures, however, the dominant H reservoir is still the ocean; this is because the H₂O storage capacity of average mantle rock is still ≪0.5 wt%. Moreover, in this paper we only consider <8 GPa seafloor pressures. To sum up, after the magma ocean has crystallized, ocean depth can be considered as a constant.
• 4.  
  C exchange between the convecting mantle and the water ocean shuts down within O(10⁸) yr after the last giant impact. To shut down C cycling between a planet's mantle and ocean, it is necessary to shut down volcanic outgassing. For this to occur, it is sufficient to shut down volcanism. (Other modes of outgassing such as geysers are ultimately driven by volcanism.) Volcanism is curtailed by seafloor pressure (Figure 5; Kite et al. 2009). The seafloor pressure threshold needed to curtail volcanism depends on how the mantle convects, and is smaller for stagnant-lid convection (∼1–3 GPa). We analyze this effect in Appendix A. For simplicity, we impose an abrupt 1 GPa cutoff on volcanism. In reality some volcanism probably continues above this seafloor pressure, but (we assume) not enough to buffer atmosphere+ocean C. Because volcanism is intimately connected to plate tectonics, it is plausible that with little or no volcanism there is little or no tectonic resurfacing (Sleep 2000a). Without tectonic resurfacing, the planet's mantle convects the same way as the mantles of Mars, Mercury, and the Moon: stagnant-lid worlds (O'Rourke & Korenaga 2012; Noack et al. 2017; Dorn et al. 2018). In stagnant-lid mode, the mantle's convecting interior is cooled from above by conduction across a stable, very viscous boundary layer (the stagnant lid) (Stevenson 2003). We assume stagnant-lid volcanism in this paper, with little or no role for plate tectonics. Plate tectonics is considered in Section 6.5.After volcanic outgassing has shut down, C is not released from the mantle into the ocean; neither is C subducted into the deep mantle, because no-volcanism stagnant-lid tectonics lacks subduction. Therefore, we treat the atmosphere+ocean+(shallow lithosphere) as a closed system with respect to C after 10⁸ yr (Figure 6).Even without sustained volcanism, water and rock still react on waterworlds. Most importantly, volcanism will occur ≲10⁸ yr after the last giant impact, for example as the dregs of the magma ocean extrude. This volcanism is key to setting ocean chemistry, and is discussed in Section 3. Even after ∼10⁸ yr, slow alteration of porous seafloor rocks can continue, due to diffusion and/or changes in seafloor T that shift mineral equilibria. This slow process, which we do not model, supplies chemical energy and nutrients (Section 6.4). Finally, carbonate precipitation and dissolution at (or just below) the seafloor may occur.
• 5.  
  No land. Mountain-peak height is set by tectonic stress in competition with gravity-driven crustal spreading (Melosh 2011). The tallest peaks in the solar system are ∼20 km above reference level. Such peaks are drowned on waterworlds, so there is no land. With no land, there is no continental weathering.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Now, for a given planet mass (M_pl), we can write

$$\begin{eqnarray}&&{p}_{\mathrm{CO}2}={p}_{\mathrm{CO}2}({f}_{W,t={t}_{0}},{c}_{C,t={t}_{0}},{[X]}_{t={t}_{0}},{L}_{* }(t),t)\end{eqnarray} \tag{ 3 }$$

where f_W is planet water mass fraction, L_* is stellar luminosity, ${[X]}_{t=0}$ is initial ocean chemistry (including mineral phases in equilibrium with the ocean, if any), and t is time in Gyr. (All terms are defined in Table 1.) The start time t₀ = 0.1 Gyr; our simulations have 0.1 Gyr timesteps (Figure 6). Continental weathering flux does not appear in Equation (3), because we assert "no land." Ocean-interior exchange flux does not appear in Equation (3), because we assert "no post-0.1 Ga C exchange between the deep mantle and the water ocean." (Figure 6). The only time dependence in Equation (3) is warming as the star brightens. If we can solve Equation (3) to get p_CO2(t), we can find the corresponding surface temperature (from Figure 2) and thus find the duration of habitable surface water as the star brightens.

In this subsection, we justify (1) treating the atmosphere+ocean CO₂ mass fraction (c_C) as a free parameter, varying between planets from 10⁻⁵ to 10⁻³ of planet mass, such that c_C/f_W varies from 10⁻⁶ to 10⁻¹; and (2) using salinities O(0.1–1) mol kg⁻¹. We define c_C to be the mass of C in the atmosphere+ocean after the magma ocean has crystallized but before any carbonate formation (excluding C in the metal core or silicate mantle), divided by the mass of the planet, and multiplied by 44/12 to get the CO₂-equivalent mass fraction. In this paper we do not consider compositions that are very C-rich, i.e., comet-like. Comets made only a minor contribution to Earth's volatile budget (Marty et al. 2016), but exoplanets might have comet-like compositions. These require a more sophisticated treatment including clathrates and other C-rich phases (Bollengier et al. 2013; Levi et al. 2014, 2017; Marounina & Rogers 2017; Ramirez & Levi 2018), which we do not attempt in this paper. Busy readers may skip to Section 3.2.

Seawater composition on stagnant-lid waterworlds is set by (1) delivery of CO₂, and (2) supply of cations via water–rock reaction in the 10⁸ yr after the last giant impact.⁹

• (1)  
  C delivery. C dissolves into liquid Fe-metal in the immediate aftermath of giant impacts, when liquid silicate and liquid Fe-metal equilibrate (Kuramoto & Matsui 1996; Hirschmann 2012; Dasgupta 2013). C is so iron-loving that, even if 2% of a 1 M_⊕ planet is composed of C (equivalent to 7 wt% CO₂!), only <100 ppm escapes the core (Figure 2 in Dasgupta 2013). The lure of the Fe-metal core for C was so strong that most of the C that existed in the Earth prior to the Moon-forming impact is believed to have been trapped in the core during the magma ocean era associated with giant impacts. Therefore, the C in our bodies (and in Earth's ocean) is thought to represent C delivered after the last giant impact (Bergin et al. 2015). C's behavior under giant impact conditions contrasts with that of H (Kuramoto & Matsui 1996). H does not partition as strongly into liquid Fe-metal under giant impact conditions as does C. Therefore, in contrast to C, H now in Earth's oceans could have been delivered to the Earth prior to the last giant impact. Indeed, the H in Earth's oceans is thought to have arrived during the main stage of planet assembly (O'Brien et al. 2014, 2018; Rubie et al. 2015; see Albarède 2009 for an opposing view). Giant impacts are stochastic, so the ratio of CO₂-equivalent C mass in the atmosphere+ocean to the planet water content (c_C/f_W) is expected to vary stochastically (Hirschmann 2016) (Section 6.5). Another source of variation in c_C/f_W is the varying extent (which we do not track) of the reaction Fe+H₂O $\to \,$ FeO+H₂. Given these uncertainties, we vary c_C from 10⁻⁵ to 10⁻³, such that c_C/f_W varies from 10⁻⁶ to ∼10⁻¹. c_C = 10⁻⁵ corresponds to an Earth/Venus-like C abundance, and we find that c_C = 10⁻³ all but ensures a p_CO2 > 40 bar atmosphere, which we deem uninhabitable.
• (2)  
  Supply of cations via water–rock reaction. Rock with the composition of Earth's oceanic crust (CaO ≈ MgO ≈ FeO ≈ 10 wt%; Gale et al. 2013) can form {Ca, Mg, Fe}CO₃ minerals, and potentially trap up to 70 (bar CO₂)/(km crust). To what extent is this potential realized? The fluid composition resulting from water–rock reactions <10⁸ yr after the last giant impact can be estimated (a) by analogy to similar worlds, (b) with reaction-transport models, and (c) by geophysical reasoning about the mass of rock available for water–rock reactions.

  • (a)  
    Previous work on extraterrestrial oceans—such as Europa and Enceladus—suggests salinities O(0.1–1) mol kg⁻¹. For example, Glein et al. (2015) and Hand & Chyba (2007) estimate salinity using spacecraft data, and McKinnon & Zolensky (2003), Zolotov (2007), and Zolotov & Kargel (2009) estimate salinity using models.
  • (b)  
    We use the aqueous geochemistry thermodynamics solver CHIM-XPT (Reed 1998) to explore reactions between water and rock (Appendix D). Our rock input for these runs has the composition of Earth's oceanic crust (Gale et al. 2013). Rocks release cations into the ocean that—even when not forming carbonates—can control ocean chemistry. Rock-sourced cations (e.g., Na⁺ and Ca²⁺) substitute in the ocean charge balance for H⁺. H⁺ removal raises pH, which in turn raises the fraction of dissolved inorganic carbon that is stored as HCO₃ and CO${}_{3}^{2-}$ rather than as dissolved CO₂ (Figure 4). As dissolved CO₂ decreases, p_CO2 decreases in proportion (by Henry's law). Two cations, Na⁺ and Ca²⁺, are used in our waterworld evolution model. Other cations contribute little to the charge balance of outlet fluid in our CHIM-XPT runs (Appendix D), and so are ignored. CHIM-XPT runs (for T ≤ 600 K) indicate that Mg and Fe likely form silicates, not carbonates, and so Mg and Fe do not draw down CO₂. Non-participation of {Mg, Fe} cuts C uptake by a factor of ∼3. CaCO₃ forms, and buffers [C] in outlet fluid to ≲0.1 mol kg⁻¹. Complete use of Ca to form CaCO₃, with only minor MgCO₃/FeCO₃, matches data for carbonated Archean seafloor basalts (Nakamura & Kato 2004; Shibuya et al. 2012, 2013).
  • (c)  
    Geophysical arguments about the mass of rock available to react with fluid may be divided into: (i) eruption history arguments, and (ii) pervasiveness arguments. We assume that rock–water reactions occur at T_seafloor < 650 K (for reasons given in Appendix C).

    • (i)  
      Eruption history. Alteration by water is much easier for extrusive rocks (lavas) than for intrusive rocks (sills and dikes). That is because intrusives have ≳10³-fold lower permeability, which denies entry to altering fluids (Alt 1995; Fisher 1998). Intrusive:extrusive ratios (I/E) on Earth are scattered around an average of ∼5:1 (6.5:1 for oceanic crust) (White et al. 2006; White & Klein 2014). Extrusives less than 0.5 km below the seafloor—corresponding to the highest permeabilities in oceanic crust—are the site of most of the alteration within the oceanic crust on Earth. Altered crust on an exoplanet could be much thicker than 0.5 km if the extrusive pile were thicker. A thicker extrusive pile is likely for stagnant-lid mode. In this mode, each lava erupts at the surface, cools and is rapidly altered, and then is buried by later lavas. Extrusives whose alteration leaves them completely leached of an element will give a waterworld ocean with an aqueous molality (denoted by square brackets) of

      $$\begin{eqnarray}&&\displaystyle \frac{23}{1000}[\mathrm{Na}]=\mathrm{Na}^{\prime} \displaystyle \frac{{z}_{\mathrm{cr}}{\rho }_{\mathrm{cr}}}{{z}_{{oc}}{\rho }_{w}}{({\rm{I}}/{\rm{E}})}^{-1}\,\end{eqnarray} \tag{ 4 }$$

      where z_oc is ocean thickness, z_cr is crust thickness, Na' = 0.028 is the Na mass concentration in the crust (Gale et al. 2013), ρ_cr = 3 × 10³ kg m⁻³ is crust density, and ρ_w ∼ 1300 kg m⁻³ is ocean density (this is the density of water at 1.8 GPa and 373 K; Abramson & Brown 2004). For z_oc = 100 km, z_cr = 50 km, and I/E = 5, this gives [Na] = 0.3 mol kg⁻¹.
    • (ii)  
      Pervasiveness. Earth's modern oceanic crust is altered mainly near veins and cracks (Alt 1995). However, more pervasive alteration of seafloor extrusive rocks took place under C-rich, perhaps warmer climates ∼3.4 Ga on Earth (Nakamura & Kato 2004). Moreover, Earth's ocean pH is thought to have been moderate for 4 Ga (pH = 7.5 ± 1.5; Halevy & Bachan 2017; Krissansen-Totton et al. 2018), and either higher pH or lower pH would be more corrosive. Finally, water is highly corrosive as a supercritical fluid, and seafloor T and P on waterworlds can reach supercritical conditions (Figure 3). Therefore, a wide range of alteration—including pervasive alteration—is possible on other worlds (Vance et al. 2007, 2016; Kelemen & Hirth 2012; Neveu et al. 2015).

Taken together, (2a)–(2c) suggest that, for waterworlds,

• 1.  
  both Na and Ca, but not Fe and Mg, participate in CO₂ drawdown if rocks interact with water; and
• 2.  
  the reasonable range of rock-layer thicknesses for complete/pervasive interaction with water is 0.1 km (Earth-like) to 50 km. The biggest value corresponds to an unusually low I/E = 1:1, pervasive alteration of extrusives, and a crust thickness of 100 km, which is large relative to models and relative to solar system data (Taylor & McLennan 2009; O'Rourke & Korenaga 2012; Plesa et al. 2014; Tosi et al. 2017). The largest values might be thought of as including weathering of impactors and meteoritic dust. For all but the very largest impacts, impact-target materials (Sleep & Zahnle 2001) will not contribute cations, because the ocean will shield the seafloor from impacts. The only ejecta will be water and impactor materials (Marinova et al. 2008; Nimmo et al. 2008).

The corresponding cation abundances (for a 100 km deep ocean, and neglecting Ca-mineral formation) are 0.003–1.4 mol kg⁻¹ Na, and 0.008–4 mol kg⁻¹ Ca. Lumping Na and Ca as moles of equivalent charge (Eq), we obtain 0.02–9 Eq/kg.

p_CO2(t) is set by sea-surface equilibration between the well-mixed ocean and the atmosphere (Figure 6). This equilibrium is rapid compared to the time steps in our model.

To find pH(T_surf) and p_CO2(T_surf), for each of a wide range of prescribed seawater cation content and dissolved inorganic C content, we use CHIM-XPT. We found (batch) equilibria in the system H₂O–HCO₃⁻–Ca²⁺–Na⁺, including formation of calcite (CaCO₃) and portlandite (Ca(OH)₂).

We consider the full range of geophysically plausible cation contents. In Section 4, we initially describe results for three cation cases: (1) [Na] ≈ 0.5 mol kg⁻¹, [Ca] ≈ 0–higher-pH, but lacks minerals; (2) [Ca] ≈ [Na] = 0–a lower-pH case; (3) [Ca] = 0.25 mol kg⁻¹, [Na] ∼ 0–higher-pH. Many of the runs in case (3) form calcite (or, for pH > 9, portlandite). We do not track the fate of Ca-mineral grains, which store most of the Ca in our higher-pH equilibrium model output. Depending on ocean depth and grain size, minerals might stay suspended in the water column, sink and redissolve, or pile up on the seafloor. Whatever the fate of individual grains, equilibrium with calcite buffers shallow-ocean (and thus atmosphere) CO₂ content for ocean chemistries that form calcite.

To find CO₂ solubilities, we used tables of CO₂ solubility in pure H₂O, fit to experimental data. Specifically, we interpolated and extrapolated in Table 3 of Carroll et al. (1991) and Table 3 of Duan & Sun (2003), also using the constraint that CO₂ solubility = 0 when p_atm = p_H2O. For p_H2O, we used the formulation in Appendix B of Duan & Sun (2003). CO₂ solubility is reduced by salts. For a 1 mol kg⁻¹ solution at 10 bar and 303 K, salting-out reduces solubility by 18% (Duan & Sun 2003). However, we do not include the direct effect of Na⁺ or Ca²⁺ on CO₂ solubility in our model, reasoning that the pH effect is more important and that the range of C considered in our model is so large that the effect of the ions on CO₂ solubility is less important. Using the CO₂ solubilities and the carbonate system equilibria from the CHIM-XPT calculations, we compute p_CO2(T_surf, [Na], [Ca], [C]). Equipped with p_CO2(T_surf, [Na], [Ca], [C]), we next calculate (for a given planet surface gravity, and fixing atmospheric molecular mass = 44 g/mole) p_CO2(T_surf, [Na], [Ca], c_C). Here, c_C is the atmosphere+ocean column C abundance and can (in principle) greatly exceed the ocean C-storage capability. To ameliorate interpolation artifacts, we smooth out each p_CO2(T_surf) curve using the MATLAB smooth function with a 30°C full-width bandpass. (Appendix E has more details about how ocean-atmosphere equilibration is calculated.)

Next, we interpolate the ${T}_{\mathrm{surf}}({L}_{* },{p}_{\mathrm{CO}2})$ results of a 1D radiative–convective climate model (Wordsworth & Pierrehumbert 2013b). Their study neglects clouds, assumes a fixed relative humidity of 1.0, and adjusts the surface albedo to a value (0.23) that reproduces present-day Earth temperatures with present-day CO₂ levels. (Comparison to the results of models that include clouds, and omit the albedo contribution of bright continents, indicates that this assumption will not affect the qualitative trends presented in the current paper; see Section 6.5.) We extrapolate the log-linear fit to p_CO2 = 3 × 10⁻⁶ bar. T_surf always increases with L_*. T_surf usually increases with p_CO2 as well, except near the outer edge of the HZ (p_CO2 > 8 bar in Figure 7). Near the outer edge of the HZ, adding CO₂ leads to cooling because the Rayleigh scattering of incoming shortwave radiation from extra CO₂ exceeds the greenhouse warming from extra CO₂. Where Wordsworth & Pierrehumbert (2013b) find multiple stable equilibria in their 1D atmospheric radiative-convective model, we use the warmer of their two solutions, on the assumption that impacts intermittently allow the atmosphere+(wave-mixed ocean) (which have relatively low thermal inertia) to jump onto the warm branch. Once the atmosphere+(wave-mixed ocean) have reached the warm branch, the deep ocean will gradually adjust to the new, higher temperature. To find ocean–atmosphere equilibria, we find {T_surf, p_CO2} combinations that satisfy L_*(t), [Na], and [Ca] (Figure 7), as explained in Section 4.1. Although the ocean T varies with depth within the ocean (Figure 3), the CO₂ solubility increases with depth, and so the ocean surface temperature T_surf is the temperature that matters for the purpose of finding ocean–atmosphere equilibria. (Appendix E gives more details, and shows the results of a sensitivity test using the 1D radiative–convective climate model of Ramirez et al. 2014.)

Zoom In Zoom Out Reset image size

At each timestep, we check for ice. Ice at the ocean surface (ice I) will have steady-state thickness <20 km (due to geothermal heat). We assume initially ice-covered surfaces have an albedo that declines due to dark exogenic contaminants (meteoritic dust) on a timescale ≪10 Gyr, back to albedo = 0.23. This value is used for convenience, as it is the surface+clouds albedo used by Wordsworth & Pierrehumbert (2013b). This lowering of initially high albedo is reasonable because, in our model, planets develop surface ice not at all, or very early. Early-developed ice cover will sweep up dark planet-formation debris (Löhne et al. 2008). The low albedo of debris is one of the reasons that all icy worlds in the solar system have dark surfaces except terrains which have been tectonically resurfaced and/or created after the stage of high impact flux in the early solar system. Lower albedo allows stellar luminosity increase to cause ice-covered surfaces to melt (Section 6.5). Melting with these assumptions always leads to a habitable state (this need not be true if CO₂ outgassing is allowed; Turbet et al. 2017; Yang et al. 2017). If ice cover develops later in planetary history, then initially high albedos can stay high, but this is very uncommon in our model. In addition to checking for surface ice, we also evaluate whether high-pressure ice phases (ice VI, ice VII) are stable at any depth between the ocean surface and the seafloor along a pure-liquid-H₂O adiabat (Appendix B, Figure 3). We do not track the climate or pH in detail when HP-ice is present (see Levi & Sasselov 2018), and in particular we do not redistribute CO₂ between clathrates, the ocean, and the atmosphere. This is because we conservatively assume that if HP-ice is present, the world is not habitable.

Charting planet evolution: With the waterworld approximation, the cations available to the ocean are constant with time after 0.1 Gyr. Thereafter, T_surf and the total amount of C in the atmosphere+ocean (including C in Ca-minerals, if any) are the sole controls on p_CO2. p_CO2 usually increases with T_surf, due to the low solubility of CO₂ in warm water, and related T-dependent shifts in the carbonate system equilibria (Figure 7). p_CO2 = p_CO2(T_surf) is a single-valued function for any cation and total-C concentration (Figure 7). T_surf = T_surf(p_CO2) greenhouse curves (green lines) can be computed using a climate model given semimajor axis a and L_*(t). To find the p_CO2 and T_surf for a given time, we locate the green line corresponding to that star age, then find its lowest-T_surf stable intersection with the appropriate p_CO2(T_surf) line (Figure 7).

As L_* increases with time (Bahcall et al. 2001), T_surf = T_surf(p_CO2, L_*) increases for a given p_CO2. This is shown by the rightward drift of the green lines (greenhouse curves) in Figure 7. The green lines are closely spaced at p_CO2 ∼ 0.2–20 bars, corresponding to a small T_surf rise for a given increase in L_* (Figure 2). This prolongs the duration of habitable surface water (longer τ_hab). By contrast the green lines are widely spaced at low p_CO2, so a small increase in L_* corresponds to a rapid T_surf rise (Figure 2). The T_surf increase is usually more than would occur without ocean chemistry (Kitzmann et al. 2015), because the greenhouse gas CO₂ is usually exsolved with increasing T_surf. For a given total atmosphere+ocean C inventory (e.g., 1.1 × 10⁶ kg m⁻² CO₂-equivalent), the p_CO2 is much higher (for a given T_surf) for the cation-poor case relative to the cation-rich cases. This is due to the higher pH of the cation-rich cases (Figure 4). However, p_CO2 increases more steeply with T_surf for the cation-rich cases.

For most waterworld-evolution tracks, each time step is independent of the previous ones. Therefore, the colored portions of the waterworld evolution tracks (ice-free oceans) in Figures 8–10 should be valid whether or not HP-ice is present earlier in that world's history (gray portions of the lines).

Zoom In Zoom Out Reset image size

Stability and instability: To see that T_surf–p_CO2 equilibria are usually stable, consider the point around 50°C, 1 bar. Suppose we add some energy to the system, perturbing T_surf. This drives CO₂ from the ocean to the atmosphere—moving along the blue line—which will cause more warming, leading to a positive feedback on the initial perturbation. However, at fixed L_*, the value of p_CO2 that is in equilibrium with the increased surface temperature is even higher than the increased CO₂. Therefore, the climate is stable.

For P > 1 bar, exsolving CO₂ can be a negative (stabilizing) feedback on climate evolution. For example, consider the point at 40°C, 10 bar in Figure 7. An increase in L_* corresponding to 1 Gyr of stellar evolution at constant p_CO2 would cause a T_surf increase of ∼40 K. However, the climate is constrained to evolve along a curve of constant seawater cation abundance. Therefore, CO₂ is exsolved from the ocean. Cooling from Rayleigh scattering per unit of added CO_2(g) exceeds greenhouse warming per unit of added CO_2(g), so adding CO_2(g) cools the planet (Ramirez et al. 2014). For the purple dashed line, the T_surf increase is only <10 K; a negative feedback. This negative feedback is not a geochemical cycle.

CO₂'s solubility in water (for fixed p_CO2) increases for T_surf > (110–160)°C. This feedback is mildly stabilizing, but we neglect $\partial {p}_{\mathrm{CO}2,\mathrm{GH}}/\partial T$ < 0 in the following discussion as these cases are rarely important.

For a {p_CO2, T_surf} point to undergo runaway warming in response to small T_surf increases, two conditions must be satisfied. (1) $\partial {p}_{\mathrm{CO}2,\mathrm{OC}}/\partial T$ > $\partial {p}_{\mathrm{CO}2,\mathrm{GH}}/\partial T$, where the subscript _OC refers to ocean chemistry curves, and the subscript _GH refers to greenhouse curves. Otherwise, the positive feedback has finite gain, and there is no runaway. (2) $\partial {p}_{\mathrm{CO}2,\mathrm{GH}}/\partial T$ > 0. If (2) is not satisfied, then there is a stabilizing, negative feedback.

Exsolution-driven climate instabilities (runaway feedbacks) are possible in the model. For example, consider the red thick-dashed line corresponding to Na = 0.25 mol kg⁻¹, c_C = 3.5 × 10⁶ kg m⁻² in Figures 8 and 9. At 2.0 Gyr and 1.1 au, two stable states are possible for these inputs: {37°C, ∼2 bar} and {75°C, ∼10 bar} (Figure 7). As L_* increases, both stable states get warmer. Then, at 2.3 Gyr, the GH curve no longer intersects the thick red dashed line for T_surf < 60°C. The cool branch has disappeared; as a result, the climate suddenly warms. (Emergence of cold solutions as solar luminosity increases is much less common in our model.) Because of interpolation artifacts and our very approximate treatment of CO₂ solubility, one should be careful about reading out quantitative gradients for specific data points from the plots. The rate of change associated with these instabilities is unlikely to risk whole-ocean microbial extinction (in our opinion), provided that both the initial climate and final climate are habitable. The model-predicted instabilities are discussed further in Appendix F.

Zoom In Zoom Out Reset image size

Duration of habitable surface water (τ_hab): In some cases, habitability can persist for many Gyr (Figures 8 and 9). As L_* rises, either T_surf eventually exceeds 450 K, or CO₂ exceeds ∼40 bar, either of which we deem sufficient to extinguish surface habitability (for a discussion, see Section 6.4). Different cation cases with the same CO₂-equivalent atmosphere+ocean C content (c_C) have different p_CO2(t) trajectories (different colors of the same thickness in Figure 9). This is because, for the same total amount of C, the fraction of C stored in the ocean differs between cation cases. However, the trajectories for different cation cases bunch together at high c_C (thick lines of all colors in Figure 9). That is because, for high total atmosphere+ocean C inventories, the ocean is saturated and most of the C is in the atmosphere. When (C in atmosphere) ≫ (C in ocean), C_atm ≈ C_total, independent of how much C is in the ocean. The colors also bunch together for very low c_C. For such worlds, the p_CO2 is too low to affect climate until T_surf > 100°C. As a result, such worlds have H₂O as the only greenhouse gas. Then, very low c_C worlds deglaciate at the same time (for a given planet orbital semimajor axis), trek quickly across the narrow "neck" at the bottom of the diagram on Figure 2, and then suffer the runaway greenhouse.

In Figure 8, cases with high CO₂-equivalent atmosphere+ocean C content have longer durations with habitable surface water (τ_hab). The longest-lived worlds maintain p_CO2 in the 0.2–20 bar range for Gyr (Figure 9).

Worlds with planet water mass fraction = 0.1 (Figure 10) show greater persistence of HP-ice. This is because f_W = 0.1 seafloor pressure is 12 GPa (for M_pl = 1 M_⊕), and these high pressures favor HP-ice (Figure 3). For a given CO₂-equivalent atmosphere+ocean C content, the climate is colder, because the bigger ocean dilutes the CO₂. For high-f_w worlds, because the ocean is the dominant reservoir of C, small fractional changes in the storage of C in the ocean can lead to large fractional changes in p_CO2.

Zoom In Zoom Out Reset image size

The initial pH of the habitable surface water ranges from 3 to >11. In our model, as atmosphere p_CO2 rises during the interval of habitable surface water, ocean pH falls (Zeebe 2012), typically by <1 pH unit. For {p_CO2, T_surf} values that allow habitable surface water, pH is higher for the "0.25 mol kg⁻¹ Ca" case and highest of all for the "0.5 mol kg⁻¹ Na" case.

In summary, ocean chemistry is key to waterworld climate. This is because of the following effects. (1) Ocean pH sets the ocean-atmosphere partitioning of C. This in turn sets the sensitivity of p_CO2 to increases in total C abundance (the Revelle factor). The Revelle factor corresponds to the differences in spacing between the y-axis positions of lines of the same color (cation abundance) but different thickness (c_C) in Figures 7–10. (2) Carbonate-system equilibria are T-dependent. This causes the upwards slope of the colored lines in Figure 7–10.

Figure 11 shows, for f_W = 0.01, M_pl  =  1 M_⊕waterworlds, the onset time, shutoff time, and duration of habitable surface water. In each case, nonzero onset times (thin solid lines) correspond to the meltback time for an early-established ice lid on the ocean. Shutoff times (thin dash–dot lines) correspond to either the runaway greenhouse or to excessive (>40 bar) p_CO2. With increasing CO₂-equivalent atmosphere+ocean C content, for all panels in Figure 11, τ_hab (thick solid lines in Figure 11) rises and then falls over an interval of ∼ 10⁶ kg m⁻², i.e., O(10²) bars of CO₂-equivalent. This exceeds the O(10) bar p_CO2 range for optimum habitability in Figure 2. That is because around the habitability optimum, as C is added to the system, most of the extra C is partitioned into the ocean and does not contribute to p_CO2 (Archer et al. 2009).

Zoom In Zoom Out Reset image size

Too much C makes the planet uninhabitably hot. Figure 13, which is a simplistic cartoon of C partitioning, shows why. Starting with a small CO₂-equivalent atmosphere+ocean C content, p_CO2 is initially low. Adding C will eventually overwhelm rock-sourced cations. The ocean then becomes acidic. Once the ocean is acidic, extra C spills into the atmosphere and p_CO2 becomes high. Supposing the ocean storage capacity for C to be fixed at ∼0.5 mol kg⁻¹ C (in reality capacity will rise with ${p}_{\mathrm{CO}2};$ Equation (1)), and for pH low enough that all aqueous C is stored as CO_2T, then ocean CO₂-equivalent content is 20 g kg⁻¹—for a 100 km deep ocean, 3 × 10⁶ kg m⁻². If extra C goes into the atmosphere, then an additional 4 × 10⁵ kg m⁻² leads to an uninhabitably hot surface. Although this "bucket model" (over-)simplifies carbonate system chemistry, it explains the ∼10⁶ kg m⁻² upper limit for habitability for the Na-free, Ca-free case (middle panel in Figure 11). The C threshold for planet sterilization is raised if cations are available (left panel and right panel of Figure 11), because they neutralize the carbonic acid. For example, adding 0.5 mol kg⁻¹ of positive charges (0.25 mol kg⁻¹ of Ca²⁺, or 0.5 mol kg⁻¹ of Na⁺) means that an extra 0.5 mol kg⁻¹ of C must be added to get nontrivial CO₂ in the atmosphere. For a 100 km deep ocean, this is 200 bar = 2 × 10⁶ kg m⁻². This quantity corresponds to the vertical offset between the habitability optimum in the middle panel, versus the habitability optimum in the two side panels, of Figure 11. This quantity of 2 × 10⁶ kg m⁻² also corresponds to the vertical offset between the upper limit for habitability in the middle panel, versus the upper limit for habitability in the two side panels, of Figure 11. Neutralization by cations also explains why τ_hab is brief (and uniform) for low C content when cation content is high. For these worlds, p_CO2 is so low (because pH is so high) that these planets effectively have H₂O+N₂ atmospheres. Examples of evolutionary tracks for such worlds include the nearly vertical tracks to the right of Figure 8. Pure-H₂O atmospheres lead to planets having habitable surface water for <1 Gyr if they orbit G-stars, according to the 1D model of Goldblatt (2015).

τ_hab is longest for semimajor axes a ∼1.1 au. For a < 1.1 au habitable surface water is available from t = 0 at the inner edge of the HZ, but the runaway greenhouse comes soon. For a ≫ 1.1 au, the surface is frozen until near the end of the model run (and the end of the main sequence, which we decree to occur at 10 Gyr for the G-type host star, arrives swiftly). For a ∼ 1.5 au, only a narrow range of high CO₂ permits habitability, corresponding to the wing-shape extending to the top right on all panels. CO₂ condensation (Turbet et al. 2017) might clip this wing of high-semimajor-axis habitability. Only around 1.1 au does habitable surface water both begin early and end late. The 1.1 au maximum in τ_hab corresponds to the "bull's-eyes" in Figure 11.

On a 10×-deeper ocean (f_W = 0.1; Figure 12), the habitability optimum shifts to 10× larger values of C. This shift is because more C is needed to overwhelm the ocean sinks (Figure 13). The range of C that gives long durations of surface liquid water appears narrower on a log scale. That is because the range from "minor p_CO2 to 40 bars p_CO2" is the same (∼4 × 10⁵ kg m⁻² C for fixed cations) in absolute terms, so smaller in fractional terms. The upper limit for habitability has also moved higher, but only by a factor of two. This is because [C]-rich tracks that reach 40 bar p_CO2 at T < 100°C always have HP ice for f_W = 0.1, and so do not count as habitable (thick stubby gray lines that terminate within the bottom left corner of Figure 10). The strong dependence of the lifetime of habitability on semimajor axis for f_W = 0.01 is more subdued for f_W = 0.1. This is because almost all f_W = 0.1 worlds start with HP-ice due to the deeper ocean (Figure 3), and such worlds are not counted as habitable in Figure 12 until temperature has risen sufficiently to melt all the high-pressure ice.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

So far we have considered three specific cation-abundance cases. Now we consider a much wider range of cation abundances.

Figure 14 shows how the duration of habitable surface water (τ_hab) depends on CO₂-equivalent atmosphere+ocean C content (c_C) and cation abundance (i.e., the extent of rock–water reaction) for semimajor axis = 1.1 au and planet water mass fraction = 0.01. For high c_C (black circles in Figure 14), C swamps all available oceanic and crustal sinks, and so piles up in the atmosphere (Figure 13). This leads to sterilizing surface temperatures. When cations released from rock dominate ocean chemistry, pH is high, and almost all C is sequestered as CO${}_{3}^{2-}$ ions, HCO₃⁻ ions, or CaCO₃ (Kempe & Degens 1985; Blättler et al. 2017). Therefore little C remains for the atmosphere, and so p_CO2 is minimal and has little role in setting planet climate. The planet then has habitable surface water for only <1 Gyr (dark blue disks at the top of Figure 14). For cation-poor oceans (lower part of plot), τ_hab is longest for the [C] that gives p_CO2 ∼ 1 bar at habitable temperatures. This is as expected, because the range of L_* that permits habitable surface water is widest for p_CO2 ∼ 1 bar (Figure 2). Durations dwindle for lower [C] and vanish entirely for higher [C] (excessive p_CO2 and/or T_surf). When cations in the ocean balance [C] (black horizontal line), the maximum in habitable-surface-water duration occurs at larger C. We suppose that the altered rock layer has the composition of mid-ocean ridge basalt, that it is completely leached of the charge-equivalent of its Na and Ca content, and that the leached thickness is log-uniformly distributed from 0.1 to 50 km (Section 3.2). For the purposes of Figures 14 and 15 we assume that carbonates do not form; if they did, then a factor of ∼2 more C would be needed to neutralize the cations. Equilibrium with minerals (e.g., CaCO₃, NaAlSi₂O₆) might buffer ocean cation content to <0.5 mol kg⁻¹, which roughly corresponds to the middle red line in Figure 14. Most points above this line have short τ_hab. Therefore, our decision to allow cation contents > 0.5 mol kg⁻¹ is geochemically less realistic, but conservative for the purpose of estimating τ_hab on waterworlds.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

To explore habitable-surface-water duration, we do a parameter sweep. We consider a hypothetical ensemble of M_pl = M_⊕, planet water mass fraction = 0.01, semimajor axis = 1.1 au worlds (Figure 15). Suppose that log(c_C) is uniformly distributed (independent of the leached thickness) from −5 to −3 (0.3%–30% of water content; C/H > 1 is cosmochemically implausible) (Section 3.2). This prior likely overweights small values of c_C, because for small values of non-core planet C mass fraction (c_C), C can be stored mainly in the silicate mantle (Hirschmann 2016). For these cosmochemically and geophysically reasonable priors, the duration of habitable surface water is shown in Figure 15. The lower panel shows that C has a negative or neutral effect on τ_hab for slightly more than half of the parameter combinations tested. About a quarter of worlds never have habitable surface water, due to excessive initial p_CO2. Another quarter have p_CO2 so low, due to efficient rock leaching, high-pH oceans, and storage of C as CO${}_{3}^{2-}$ ion or HCO₃⁻ ion, that CO₂ has minor climate effect. As a result, the (short) duration of habitable surface water is about the same as for an N₂+H₂O atmosphere, and <1 Gyr. The remaining half of worlds have longer τ_hab due to CO₂. In ∼9% of cases, the waterworlds stay habitable for longer than the age of the Earth.

Figure 15 also shows the results for planet water mass fraction = 0.1. Durations of habitable surface water are shorter, because HP-ice can only be avoided for a narrow (steamy) range of temperatures. Reducing the maximum temperature for life from 450 to 400 K would reduce τ_hab to <1 Gyr on f_W = 0.1 worlds. For deeper oceans, fewer worlds are sterilized by too-high CO₂ values: for these worlds, the solution to pollution is dilution.

Results depend on the mean value of the ratio of CO₂-equivalent atmosphere+ocean C content to planet water content (i.e., c_C/f_W). If c_C often exceeds 10⁻³, then a greater fraction of waterworlds will have thick CO₂-rich atmospheres and be too hot for life.

Mean waterworld lifetime is ∼2 Gyr (area under the curves in the lowest panel of Figure 15). This is not very much less than the maximum planet lifetime for a = 1.1 au, M_pl = 1 M_⊕, which is 7 Gyr. After 7 Gyr, the runaway greenhouse occurs even on planets where geochemical cycles aid habitability (Caldeira & Kasting 1992; Wolf & Toon 2015). Therefore, cycle-independent planetary habitability is not much less effective at sustaining habitable surface water relative to a hypothetical geochemical cycle that maintains p_CO2 at the best value for habitability.

The effect of C on habitability in our model is positive (Figure 15, lowest panel). Our parameter sweep suggests that the median τ_hab is comparable in the random-allocation-of-C case relative to a C-starved case. Maximum τ_hab is greatly increased (Carter 1983). This gain outweighs the losses due to planets that form with too much C for habitable surface water. Furthermore, life requires C, so very low [C] in the ocean could hinder life's origin and persistence (Kah & Riding 2007).

We perform N-body simulations including simple gravitational interactions with a smooth, 1D gas disk with a surface density proportional to ${r}^{-3/2}$, but not including the back-reaction of disk perturbations on the planets. This results in a simple strong-migration scenario, similar to that in many other studies (e.g., Cossou et al. 2014; Ogihara et al. 2015; Izidoro et al. 2017; Sun et al. 2017; Raymond et al. 2018; Ronco & de Elía 2018). For the purposes of this study, the N-body simulations serve to provide illustrative migration and collisional histories. Therefore, details of the N-body simulations are not essential; some readers may choose to skip to Section 5.2.

We employ the REBOUND N-body code (Rein & Liu 2012) and its IAS15 integrator (Rein & Spiegel 2015), with additional user-defined forces to represent gravitational tidal drag from the disk (Zhou & Lin 2007). Each simulation is initialized with dozens of planetary embryos near or beyond the H₂O snowline, orbiting a solar-mass star. No gas giants exist at the start of our simulation, consistent with the wetter-than-solar system cases we model (Batygin & Laughlin 2015; Morbidelli et al. 2016). The following parameters are varied across simulations: the total mass in embryos (from 10 to 40 M_⊕); the outer edge of the disk (from 3 to 10 au); and a scale factor for the mass of the gas disk relative to a fiducial minimum mass solar nebula (from 2 to 16 times more massive than the minimum mass solar nebula).

For each set of model parameters, we generate multiple sets of initial conditions, resulting in a total of >145 N-body runs. Embryo initial masses are ∼1 Mars mass. The initial semimajor axes are drawn from a power-law distribution with a power-law index of −1.5. Orbits are initialized with modest eccentricities and low inclinations (relative to the plane of the gas disk). The initial orbital eccentricities and inclinations are drawn from Rayleigh distributions with a scale parameter of 0.1 for eccentricity, and 0.1 rad for inclination. In practice, the eccentricities and inclinations within a system rapidly damp to a level where damping is balanced by N-body excitation from neighbors. This choice of initial conditions helps to avoid a very brief phase of rapid collisions before the system is able to relax to a physically plausible set of inclinations and eccentricities. The three remaining angles (argument of pericenter, longitude of ascending node, mean anomaly) were drawn from uniform distributions ($0\leftrightarrow 2\pi$).

While the gas disk is present, orbital migration is parameterized following Equations (19) and (20) from Zhou & Lin (2007). The gravitational accelerations due to tides from the gas disk are proportional to: (a) the difference in the velocity of the planet (or embryo) from the circular velocity at the planet's current position (i.e., proxy for the gas velocity), (b) the planet–star mass ratio, (c) a 1/a² term that accounts for a ${a}^{-3/2}$ power-law dependence of the gas-disk surface density, and (d) a time-dependent factor that exhibits first slow dispersal, and then rapid dispersal due to photoevaporation. This prescription results in orbital migration, eccentricity damping, and inclination damping that have self-consistent timescales (Zhou & Lin 2007). For each of our simulations, the gas-disk mass initially undergoes exponential decay with a characteristic timescale of 2 Myr. Starting at 3 Myr, the gas-disk decay accelerates due to photoevaporation, so the exponential decay rate is decreased to 50,000 yr. After 4 Myr, the gas-disk mass has essentially dispersed, so it is set to zero and the simulation becomes a pure N-body integration.

Following gas disk dispersal, planets continue to scatter and collide with one another. We consider planets as having collided when they come within a radius corresponding to a density of 0.125× that of Mars. Collisions are assumed to be perfect mergers in the N-body code.

Of our simulations, 29% of the runs resulted in one or more planets in the HZ, based on the zero-age main sequence luminosity of the Sun (∼0.7× present solar luminosity). We set aside N-body simulations where no planets survived exterior to the HZ, as the results of those simulations may be impacted by edge effects (i.e., if embryos formed beyond the outer end of our initial set of embryos). In some cases, multiple planets in the HZ at the end of the simulation have semimajor axes differing by less than <0.05 au from one another. In these cases, we merge such planets, as they are expected to collide if we were able to extend each of the N-body simulations. We do not include the loss of water due to these (final) giant impacts. After applying these filters, we are left with ∼30 HZ planets.

Embryo H₂O content is set as follows. We represent the cumulative effects of the evolving H₂O-ice snowline on the water mass fraction of planetary embryos ${f}_{W,o}$ (Piso et al. 2015; Hartmann et al. 2017) by

$$\begin{eqnarray}{f}_{W,o} & = & {f}_{W,\max }(a-{a}_{\mathrm{sl}}/{a}_{\mathrm{sl}}),{a}_{\mathrm{sl}}\lt a\lt 2{a}_{\mathrm{sl}}\\ {f}_{W,o} & = & {f}_{W,\max },\ a\gt 2{a}_{\mathrm{sl}}\\ {f}_{W,o} & = & 0,a\lt {a}_{\mathrm{sl}}\end{eqnarray} \tag{ 5 }$$

where a_sl = 1.6 au is an anchor for the evolving snowline (Mulders et al. 2015). For each N-body simulation, we consider f_W,max = {50%, 20%, 5%, 1%}. The highest value, 50%, corresponds to theoretical expectations from condensation of solar-composition gas (Ciesla et al. 2015). Intermediate values match density data for Ceres (Thomas et al. 2005) and inferences based on Earth's oxidation state (Rubie et al. 2015; Monteux et al. 2018). Five percent corresponds to meteorite (carbonaceous chondrite) data. Why are the carbonaceous chondrites relatively dry? One explanation involves heating of planetesimals by short-lived radionuclides (SLRs; ²⁶Al and ⁶⁰Fe), leading to melting of water ice. Meltwater then oxidizes the rocks, and H₂ returns to the nebula (e.g., Rosenberg et al. 2001; Young 2001; Monteux et al. 2018). Because most protoplanetary disks probably had fewer SLRs than our own (Dwardakas et al. 2003; Gaidos et al. 2009; Gounelle 2015; Lichtenberg et al. 2016; see also Jura et al. 2013), this process would operate differently (or not at all) in other planetary systems. Another explanation for the dryness of carbonaceous chondrites relies on Jupiter's early formation (Morbidelli et al. 2016); gas giants of Jupiter's size and semimajor axis are not typical. The range of f_W,max we consider corresponds to processes occurring during growth from gas+dust to embryos. We neglect H₂O loss by XUV-stripping at the pre-embryo stage (Odert et al. 2018).

Growing planets shed water during giant impacts. For volatile-rich planets, at the embryo stage and larger, the most plausible mechanism for losing many wt% of water is giant impacts. The shift to giant-impact dominance can be seen by comparing Figures 14 and 15 of Schlichting et al. (2015), and extrapolating to the smallest volatile layer considered in this paper (10⁸ kg m⁻²). Such deep volatile layers greatly inhibit volatile loss by r < 50 km projectiles, and blunt the escape-to-space efficiency of larger projectiles. (For volatile-poor planets, by contrast, small impacts are extremely erosive; Schlichting et al. 2015.) Combined with the expectation that most of the impacting mass is in the largest projectiles, erosion by giant impact is the main water loss process for waterworlds.

H₂O is attrited by individually tracked giant impacts (Genda & Abe 2005; Marcus et al. 2010b; Stewart et al. 2014), using the relative-velocity vectors from the N-body code. The dependence on relative velocity of fractional ocean loss is taken from Stewart et al. (2014). Because the parameterization of Stewart et al. (2014) is for f_W ∼ 10⁻⁴, and ocean loss efficiency decreases with ocean thickness (Inamdar & Schlichting 2016), this is conservative in terms of forming waterworlds. More recent calculations find lower ocean-loss values (Burger et al. 2018). When embryos collide with other embryos, or with planets (Maindl et al. 2017; Golabek et al. 2018), Fe-metal cores are assumed to merge efficiently. We do not track protection of H₂O from giant impacts by dissolution within deep magma oceans formed during earlier giant impacts (Chachan & Stevenson 2018). This is also conservative in terms of forming waterworlds.

Planet assembly is relatively gentle in our model. Most planet-forming collisions occur while disk gas maintains low planet eccentricities and inclinations. As a result, $| {v}_{\infty }|$ is small, and because giant-impact loss is sensitive to $| {v}_{\infty }|$, these early collisions have individually minor effects on planet mass. When f_W ≥ 0.15, an individual giant impact does not remove a large percentage of the initial planet volatile endowment (Figure 16).

Zoom In Zoom Out Reset image size

Planet water mass fraction is anti-correlated with M_pl, due to the larger number of devolatilizing giant impacts needed to assemble larger planets. For a given M_pl, worlds assembled via collisions between equal-mass impactors have higher f_W than worlds that are built up piecemeal (one embryo added a time). This is because piecemeal assembly requires a greater number of giant impacts (Inamdar & Schlichting 2016). There is no obvious trend of f_W with a (nor M_pl with a) within the HZ.

The final planet water mass fractions are sensitive to f_W,max, as expected (Ciesla et al. 2015; Mulders et al. 2015). For a given f_W,max, seafloor pressures are bunched because both the weight per kg of water, and the number of water-removing giant impacts, increase with planet mass. None of the f_W,max = 0.05 or f_W,max = 0.01 planets have seafloor pressures >1 GPa. By contrast, 97% of the f_W,max = 0.2 worlds have 1–8 GPa seafloor pressure. Finally, of the f_W,max = 0.5 worlds, 7% have seafloor pressure <1 GPa, 3% have seafloor pressure ∼6 GPa, and 90% have seafloor pressure > 8 GPa. These relatively low seafloor pressures are the consequence of low surface gravity: M_pl averages 0.3 M_⊕ for this ensemble. To calculate the gravitational acceleration experienced by the ocean, we follow Valencia et al. (2007) and let ${R}_{\mathrm{seafloor}}/{R}_{\oplus }=({M}_{\mathrm{pl}}/{M}_{\oplus }){}^{0.27}$. If more massive planets had emerged in our simulations, then we would expect seafloor pressures to be higher for a given f_W,max.

We injected the {a, M_pl, f_W} tuples from the f_W,max = 0.2 planet assembly run into the waterworld evolution code, varying the CO₂-equivalent atmosphere+ocean C content and varying cation-leaching using the same procedure and the same limits as in Section 4. The results (Figure 17) show that habitable surface water durations > 1 Gyr occur for ∼20% of the simulated parameter combinations. τ_hab > 4.5 Gyr occurs for ∼5% of the simulated parameter combinations. Only around a third of trials never have habitable surface water. Never-habitable outcomes can correspond to always-icy surfaces (especially near the outer edge of the HZ), or to p_CO2 ≳ 40 bars. Never-habitable outcomes become less likely as planet water mass fraction increases, mainly because of dilution. These results are illustrative only because mean (and median) M_pl for this ensemble is 0.3 M_⊕. The lower gravity of 0.3 M_⊕ worlds relative to our 1 M_⊕ reference will move the inner edge of the HZ further out than is assumed in Figure 17 (Kopparapu et al. 2014).

Zoom In Zoom Out Reset image size

Waterworlds cannot be distinguished from bare rocks by measuring radius and density. Even if radius and density measurements were perfect, a small increase in core mass fraction could offset a 100 km increase in ocean depth (Rogers & Seager 2010). Indeed, scatter in core mass fraction is expected for planets assembled by giant impacts (Marcus et al. 2010a). The best-determined planet radius is that of Kepler-93b; 1.481 ± 0.019 R_⊕ (Ballard et al. 2014), an error of 120 km in radius. Yet a 120 km deep ocean is enough to be in the cycle-independent planetary regime discussed in this paper, and a 0 km deep ocean is a bare rock. Small uncertainties in star radius, or in silicate-mantle composition (e.g., Dorn et al. 2015; Stamenković & Seager 2016) can also affect retrievals. Therefore, all small-radius HZ planets are potentially waterworlds, and this will remain the case for many years (Dorn et al. 2017; Simpson 2017; Unterborn et al. 2018).

The main climate feedback in our code is the decrease in CO₂ solubility with increasing temperature, which is usually destabilizing (Section 4.1). We now speculate on possible feedbacks, not included in our code, that might stabilize climate.

• 1.  
  When carbonates dissolve, CO₂ is consumed (Zeebe 2012). This is because (for circumneutral pH) one mole of Ca²⁺ charge-balances two moles of HCO₃⁻, whereas CaCO₃ pairs one mole of Ca with only one mole of C (Figure 4). Carbonates dissolve with rising T (Dolejš & Manning 2010), so the corresponding drawdown of CO₂ is a climate-stabilizing feedback. This feedback is curtailed on Earth by the build-up of a lag of insoluble sediment, for example continent-derived silt and clay. Insoluble sediments would be in short supply on waterworlds, and so the carbonate-dissolution feedback could be stronger on waterworlds.
• 2.  
  Pervasive aqueous alteration lowers the density of the upper crust. Low density of altered layers may cause rising magma to spread out beneath the altered layer, failing to reach the surface as an extrusive lava, and instead crystallizing in the subsurface as an intrusion. Intrusions have low permeability and so are resistant to alteration (in contrast to extrusive lavas that have high permeability). This negative feedback on aqueous alteration might limit cation supply.
• 3.  
  When high-pressure ice forms, it excludes C from the crystal structure. This raises the C concentration in the ocean, and thus p_CO2. This negative feedback on cooling will have a strength that depends on the extent to which C is taken up by clathrates (Levi et al. 2017). Moreover, according to Levi et al., the solubility of CO₂ in water in equilibrium with CO₂ clathrate decreases with decreasing temperature. This solubility behavior is an additional negative feedback.
• 4.  
  High C concentration lowers pH. This favors rock leaching, which is a negative feedback on p_CO2 (Schoonen & Smirnov 2016). But because the water/rock ratio is effectively a free parameter over many orders of magnitude, the effect of water/rock mass ratio is more important in our model (Figure 18). Moreover, high-pH NaOH solutions corrode silicate glass, a positive feedback.
• 5.  
  It has been suggested that carbonate in seafloor basalt is an important C sink on Earth, and has a temperature-dependent formation rate, and so might contribute to Earth's climate stability (Coogan et al. 2016).

These feedbacks might extend τ_hab beyond the already long durations calculated in Section 4.

Zoom In Zoom Out Reset image size

The number of planet-years of habitable surface water enabled by cycle-independent planetary habitability on waterworlds (T_nocycle,ww) can be estimated as follows:

$$\begin{eqnarray}&&{T}_{\mathrm{nocycle},\mathrm{ww}}=\displaystyle \sum _{i=1}^{{\eta }_{\oplus }{N}_{* }}{\tau }_{\mathrm{hab},i}({f}_{W,i},{c}_{C,i},{M}_{{pl},i},\ \ldots )\end{eqnarray} \tag{ 6 }$$

where N_* is the number of Sun-like stars in the Galaxy, and ${\eta }_{\oplus }$ (≈0.1) is the number of HZ small-radius exoplanets per Sun-like star. Given a τ_hab look-up table (e.g., Figure 15), Equation (6) can be evaluated (e.g., Figure 17) given estimates (from plate formation models constrained by data) of the distribution of f_W, c_C, M_pl, and so on. How does T_nocycle,ww compare to the number of planet-years of habitable surface water maintained by geochemical cycles (T_cycle)? We have

$$\begin{eqnarray}&&\displaystyle \frac{{T}_{\mathrm{cycle}}}{{N}_{* }{\eta }_{\oplus }}\sim \,{x}_{\mathrm{IC}}\,{x}_{\mathrm{evo}}\,5\,\mathrm{Gyr}.\end{eqnarray} \tag{ 7 }$$

Here, x_IC is the fraction of small-radius HZ exoplanets with suitable initial conditions for habitability maintained by geochemical cycles, x_evo accounts for planets where geochemical cycles initially maintain habitability but subsequently fail to do so, and 5 Gyr is a typical maximum duration of habitable surface water (before the runaway greenhouse) for a planet initialized at random log(a) within the HZ of a Sun-like star (Figure 2).

Despite the familiarity of habitability maintained by geochemical cycles, it is possible that T_cycle ≪ T_nocycle,ww. This is because x_IC and x_evo are both very uncertain. This uncertainty is because predicting geochemical cycles involving habitable environments using basic models is difficult. This difficulty can be traced to the low temperatures of habitable environments. These low temperatures bring kinetic factors—such as grain size, permeability, catalysts, and rock exposure mechanisms—to the fore (White & Brantley 1995). It is also difficult to forensically reconstruct the processes that buffer Earth's p_CO2 (Broecker & Sanyal 1998; Beerling 2008). Negative feedbacks involving geochemical cycles probably buffer the post-0.4 Ga atmospheric partial pressure of CO₂ (p_CO2) in Earth's atmosphere (e.g., Kump et al. 2000; Zeebe & Caldeira 2008; Stolper et al. 2016). The effectiveness and underpinnings of these feedbacks on Earth are poorly understood (Edmond & Huh 2003; Maher & Chamberlain 2014; Galy et al. 2015; Krissansen-Totton & Catling 2017). Therefore, it is hard to say how often these feedbacks break down or are overwhelmed due to the absence of land, supply limitation/sluggish tectonics, inhibition of carbonate formation due to low pH or to SO₂, contingencies of tectonic or biological evolution, and so on (Kopp et al. 2005; Bullock & Moore 2007; Halevy & Schrag 2009; Abbot et al. 2012; Edwards & Ehlmann 2015; Foley 2015; Flament et al. 2016; Genda 2016; Galbraith & Eggleston 2017; Foley & Smye 2018). Without a mechanistic understanding of planet-scale climate-stabilizing feedbacks on Earth, we cannot say whether these feedbacks are common among planets, or very uncommon (Lacki 2016).

Is 450 K a reasonable upper T limit for habitability? A temperature of 450 K exceeds the highest temperature for which life has been observed to proliferate (395 K, Takai et al. 2008). However, no fundamental barriers have been identified to life at >400 K (COEL 2007). Theory tentatively suggests 423–453 K as the upper T limit for life in general (Bains et al. 2015). If we lower the limit for life to 400 K, then our model worlds with planet water mass fraction (f_W) = 0.1 will almost always have HP-ice at times when their surface temperature is habitable (Figure 10). It is conceivable that Earth's T_surf was >350 K at >3 Ga based on isotopic data (e.g., Tartèse et al. 2017, but see also Blake et al. 2010) and on the T tolerance of resurrected ancestral proteins (Akanuma et al. 2013; Garcia et al. 2017).

Is 40 bar a reasonable upper p_CO2 limit for habitability? We use this upper limit mainly because, above this value, fine-tuning of L_* is increasingly required to avoid uninhabitably high T_surf (Figures 2, 19). Therefore, 40 bar is a conservatively low upper limit. Above 70 bar and 300 K, CO₂ is a supercritical fluid. Supercritical CO₂ is an excellent organic solvent and bactericidal agent. However, within the ocean, life might persist. We consider near-sea-surface habitability in this paper, but note that life proliferates at 1 GPa (Sharma et al. 2002).

Extreme pH values destabilize simple biological molecules, but life on Earth solves these problems (COEL 2007). The pH requirements for the origin of life might be more restrictive.

Life as we know it requires nutrients other than C, e.g., P. Yet it is unclear to what extent life requires resupply of such nutrients (Moore et al. 2017). To the contrary, given an initial nutrient budget and a sustained source of photosynthetic energy, a biosphere might be sustained by heterotrophy, i.e., biomass recycling. Biomass recycling is easier for waterworlds that lack land and ocean-mantle geochemical cycling, because biomass cannot be buried by siliciclastic sediment, nor subducted by tectonics. If biomass recycling is incomplete (Rothman & Forney 2007), then to stave off productivity decline will require nutrient resupply. Nutrient resupply will be slower (relative to Earth) on waterworlds that lack ocean–mantle geochemical cycling, because land is absent (so sub-aerial weathering cannot occur), and volcanism is absent or minor. Moreover, the larger water volume will dilute nutrients. In our model, nutrients are supplied by water–rock reactions at the seafloor early in planet history, with limited nutrient resupply thereafter (by slow diffusive seafloor leaching, or minor ongoing volcanism). Resupply might be shut down entirely by HP-ice. However, if HP-ice both convects and is melted at its base (as on Ganymede?), then nutrient supply to the ocean may continue (e.g., Kalousová et al. 2018; Kalousová & Sotin 2018). Alternatively, nutrients can be supplied to habitable surface water by bolides and interplanetary dust, just as Earth's seafloor is nourished by occasional whale-falls. On a planet with habitable surface water, chemosynthesis (which requires an ongoing supply of rock) is greatly exceeded as a potential source of energy by photosynthesis. If this occurs, then (given the likelihood of nutrient recycling), it has not clear how potentially observable parameters would scale with geological nutrient flux. Therefore, it is hard to make a testable prediction about the nutrient budgets of water-rich exoplanets that have habitable surface water.

Low-temperature hydrothermal vents, which are implied by our models, are a frequently proposed site for the origin of life (e.g., Russell et al. 2010).

Astrophysics. It would be difficult to gather enough water to form a waterworld if HZ planets form without a major contribution of material from beyond the nebula snowline.

If Fe-metal and atmophiles do not mix during giant impacts, then it is difficult to trap C in the core (Rubie et al. 2015). This would lead to an atmosphere+ocean that inherits the C/H anticipated for nebular materials, which is more C-rich than considered here. The reaction Fe⁰ + H₂O $\to \,{\mathrm{Fe}}^{2+}{\rm{O}}$ + H₂, with subsequent escape-to-space of H, could also raise C/H.

Our results are sensitive to the assumption that initially ice-coated worlds gather enough meteoritic debris during the late stages of planet accretion to lower albedo. If we are wrong and albedos of initially ice-coated worlds all stay so high that surface ice does not melt, then the only habitable waterworlds will be those that are unfrozen for zero-age main sequence stellar luminosity (i.e., a ≲ 1.1 au). For a ≲ 1.1 au, the runaway greenhouse occurs relatively early in stellar main sequence evolution, so if surface ice does not melt, then most waterworlds with habitable surface water will orbit stars younger than the Sun. Because many HZ planets have a > 1.1 au, then if surface ice does not melt that would lower the number of waterworlds with habitable surface water by a factor of >2.

We assume that our worlds lack H₂, so it is reasonably self-consistent to assume waterworlds have CO₂ (and not CH₄) atmospheres. If some H₂ remains, then it could fuel early hydrodynamic escape of H₂O, and potentially aid Fischer–Tropsch type synthesis of organics (Gaidos et al. 2005). We assume a fixed p_N2, ∼0.8 bar, but variation in the abundance of N volatiles may be important, especially if planets draw material from a ≳ 10 au (Atreya et al. 2010; Marounina et al. 2018). N₂ affects T_surf by pressure broadening, Rayleigh scattering, etc. (Goldblatt et al. 2013), and NH₃ lowers the freezing-point. On Earth, the air–ocean partition of N₂ is ∼10³:1 (Ward 2012); a very massive ocean becomes a major N₂ reservoir.

We neglect tidal heating (Henning et al. 2009). There is no evidence that tidal heating exceeds radiogenic heating for any planet in the HZ of an FGK star older than 100 Myr.

Geophysics. In order to make our calculations simpler, we understate the pressure needed to completely suppress C exchange between the convecting mantle and the water ocean. This pressure could be closer to 3 GPa than 1 GPa (Noack et al. 2016).

We assume stagnant-lid tectonics. We suspect stagnant lid is the default mode of planetary tectonics, in part because stagnant-lid worlds are common in the solar system whereas plate tectonics is unique, and in part because of the difficulty of explaining plate tectonics (Korenaga 2013). Also, plate tectonics might require volcanism (Sleep 2000a), and volcanism is curtailed on waterworlds (Kite et al. 2009). What if plate tectonics nevertheless operates on waterworlds (Noack & Breuer 2014)? In that case, relative to stagnant-lid mode the lithosphere might store less C. That is because plate tectonics involves a thin carbonated layer that is continually subducted and decarbonated (Kelemen & Manning 2015), whereas stagnant-lid mode slowly builds-up a potentially thick layer of thermally stable carbonated rocks (Pollack et al. 1987; Foley & Smye 2018). For a planet with plate tectonic resurfacing (or heat-pipe tectonics), the ocean may interact with a large (time-integrated) mass of rock (Sleep & Zahnle 2001; Valencia et al. 2018). Increased rock supply makes it more likely that ocean chemistry would reach equilibrium with minerals (Sillen 1967). If plate tectonics operated without volcanism, then the seafloor would be composed of mantle rocks. These, when aqueously alterated, produce copious H₂ (Klein et al. 2013). Copious H₂ production would also occur if the crust extruded at the end of magma ocean crystallization was very olivine-rich (Sleep et al. 2004). The alteration of mantle rocks to form serpentine would lead to a Ca²⁺-rich ocean, with high pH (Frost & Beard 2007). In addition, H₂ production would affect atmospheric composition and climate (e.g., Wordsworth & Pierrehumbert 2013a).

Geochemistry and climate. We use results from the 1D climate model of Wordsworth & Pierrehumbert (2013b) in this study (e.g., Figure 2), which uses HITRAN H₂O absorption coefficients and a fixed relative humidity of 1.0. Subsequent work by Ramirez et al. (2014) uses improved (HITEMP) H₂O absorption coefficients and an adjustable relative humidity. The results of a sensitivity test using the Ramirez et al. (2014) work are shown in Appendix E. Recent 3D GCM results validate the trends shown in Figure 2 (e.g., Leconte et al. 2013; Wolf & Toon 2015; Popp et al. 2016; Kopparapu et al. 2017; Wolf et al. 2017). However, the amplitudes, location, and even existence of the instabilities discussed in Section 4 may depend on the specifics of the climate model, including how relative humidity and cloud cover respond to changing T_surf (e.g., Yang et al. 2016).

Our treatment of both carbonate formation and CO₂ solubility in cation-rich fluids in Section 4 is crude and this could shift the optimum CO₂-equivalent atmosphere+ocean C content for habitability by a factor of ∼2.

We neglect Cl. We do so because Cl is geochemically much less abundant than Na and C (Grotzinger & Kasting 1993; Sleep et al. 2001; Sharp & Draper 2013; Clay et al. 2017). Because we neglect Cl, our model understates the acidity of low-C-content worlds. Therefore, the vertical stripe of "optimal c_C" may be drawn at too high a value of CO₂-equivalent atmosphere+ocean C content. The H₂O/Cl ratio depends on the mechanism for water loss, if any. Production of H₂ on planetesimals by oxidation of Fe⁰ and Fe²⁺, followed by escape-to-space of H₂, gets rid of water but not Cl. Mechanical ejection of water to space during impacts gets rid of both.

Although basalt with an elemental composition similar to that of Earth's oceanic crust—as used in Appendix D—is ubiquitous in the solar system (Taylor & McLennan 2009), the silicate Earth is depleted in Na relative to meteorites that are believed to represent planetary building blocks (the chondrites) (Palme & O'Neill 2014). The chondrites are themselves depleted in Na relative to equilibrium condensation calculations. A waterworld formed from low-T condensates, or formed from chondrites but without Na depletion, would have Na:Al ∼1. Peralkaline magma would be common, stabilizing Na silicates. Na-silicates dissolve readily, releasing aqueous Na. A Na-rich ocean would have high pH and low p_CO2. With Na-carbonates exposed at the surface, the dwarf planet Ceres may be a local example of a Na-rich ocean world (Carrozzo et al. 2018).

Large-amplitude stochastic differences in planet composition can result from a collision (versus a near miss) with an individual planetary embryo.

Our modeled water–rock reactions, which use freshwater as the input fluid, raise pH due to release of (for example) Na⁺ and Ca²⁺ (Appendix D) (MacLeod et al. 1994). Na is released from seafloor rocks (Staudigel & Hart 1983) in our model, consistent with experiments (Gysi & Stefánsson 2012) and Na-leaching of preserved 3.4 Ga seafloor crust (Nakamura & Kato 2004). Moreover, on waterworld seafloors, jadeite is the stable Na-hosting phase, and (at least for T > 400°C) the jadeite–albite transition is a maximum in [Na] (Galvez et al. 2016). Na release into waterworld oceans is consistent with Na uptake at Earth hydrothermal vents (Seyfried 1987; Rosenbauer et al. 1988), because Earth seawater is Na-rich (0.5 mol kg⁻¹), which favors Na uptake in hydrothermal systems. In summary, waterworld ocean Na concentrations O(0.1–1) mol kg⁻¹ are geologically reasonable.

Our analysis suggests the following science opportunities for waterworld research.

• 1.  
  More detailed modeling of water–rock reaction throughout waterworld evolution, using newly available constraints (e.g., Pan et al. 2013).
• 2.  
  New measurements (via numerical and/or lab experiments) of solubilities, etc., for high-pressure, low-temperature conditions relevant to waterworlds.
• 3.  
  Map out rocky-planet atmospheres and atmosphere–ocean equilibria in volatile-abundance space, taking into improving constraints on the galactic range of H/Cl, H/N, and H/C (Bergin et al. 2015). Our results are sensitive to the mean value of the ratio (c_C/f_W) of CO₂-equivalent atmosphere+ocean C content to planet water mass fraction.
• 4.  
  Seek to confirm or reject hints of life proliferating at T > 400 K on Earth (Schrenk et al. 2003). If life at T > 400 K is confirmed, this would underline the habitability of worlds with planet water mass fraction (f_W) = 0.1, regardless of whether or not HP-ice is sufficient to prevent habitability (Figure 10).
• 5.  
  Investigate the long-term stability of rock–water mixtures on waterworlds. During giant impacts, rock and water are fully miscible. On planets with f_W ≪ 0.01, the mixture cools quickly, and the rock settles out, because rock is insoluble in cool water. A rock–water-mixture lifetime cannot exceed the conductive cooling timescale,

  $$\begin{eqnarray}&&{\tau }_{\mathrm{unmix}}\lesssim \displaystyle \frac{\rho {c}_{p,f}{\rm{\Delta }}Z{\rm{\Delta }}{z}_{\mathrm{bl}}}{{(2.32)}^{2}{k}_{\mathrm{cond}}}\sim 1\,\mathrm{Gyr}{\left(\displaystyle \frac{{\rm{\Delta }}Z}{200\mathrm{km}}\right)}^{2}\end{eqnarray} \tag{ 8 }$$

  where ρ is characteristic fluid density (∼1100 kg m⁻³), c_p,f is fluid heat capacity (∼4 × 10³ J kg⁻¹ K⁻¹), ${\rm{\Delta }}Z$ is ocean depth, ${\rm{\Delta }}{z}_{\mathrm{bl}}$ is the thickness of the conductive boundary layer, k_cond is characteristic fluid thermal conductivity (1 W m⁻¹ K⁻¹), and we have made the upper-limit assumption ${\rm{\Delta }}Z\,=\,{\rm{\Delta }}{z}_{\mathrm{bl}}$ (Turcotte & Schubert 2011). On worlds with large planet water mass fraction, a parcel of the dense, hot, rock+water layer is stable to rapid convective swapping with the overlying cool rock-poor aqueous fluid. Convective inhibition slows mixture cooling and allows super-adiabatic temperature gradients (persistence of hot layers at depth). Thus, an ocean at habitable temperature might be underlain not by a solid seafloor but instead by a hot rock–water mixture. The true lifetime will be set by double-diffusive convection (Radko 2013; Friedson & Gonzales 2017; Moll et al. 2017), taking account of the exotic behavior of aqueous fluids at P > 10 GPa (e.g., Redmer et al. 2011; Tian & Stanley 2015), and could be much less than the Equation (8) upper limit. It is not obvious whether such a system would be more or less habitable than a fluid ocean underlain by a solid crust.
• 6.  
  Combine eH modeling (Schaefer et al. 2016; Wordsworth et al. 2018) with pH modeling (this study).
• 7.  
  Evaluate the climate feedbacks proposed in Section 6.2. Explore the bestiary of exsolution-driven climate feedbacks suggested by our simple model (Section 4, Figure 7, Kite et al. 2011). Determine if any apply to Earth history.