Carbonate-silicate cycle predictions of Earth-like planetary climates and testing the habitable zone concept

Abstract

In the conventional habitable zone (HZ) concept, a CO₂-H₂O greenhouse maintains surface liquid water. Through the water-mediated carbonate-silicate weathering cycle, atmospheric CO₂ partial pressure (pCO₂) responds to changes in surface temperature, stabilizing the climate over geologic timescales. We show that this weathering feedback ought to produce a log-linear relationship between pCO₂ and incident flux on Earth-like planets in the HZ. However, this trend has scatter because geophysical and physicochemical parameters can vary, such as land area for weathering and CO₂ outgassing fluxes. Using a coupled climate and carbonate-silicate weathering model, we quantify the likely scatter in pCO₂ with orbital distance throughout the HZ. From this dispersion, we predict a two-dimensional relationship between incident flux and pCO₂ in the HZ and show that it could be detected from at least 83 (2σ) Earth-like exoplanet observations. If fewer Earth-like exoplanets are observed, testing the HZ hypothesis from this relationship could be difficult.

Introduction

Newton first alluded to the concept of a stellar habitable zone (HZ) in his 1687 Principia¹ by noting that Earth’s liquid water would vaporize or freeze at the orbits of Mercury and Saturn, respectively². Later, Whewell noted that “the Earth’s orbit is in the temperate zone of the Solar System”³. Since then, the definition of the stellar HZ has been refined, reaching its modern incarnation based on climate models^4,5.

Current HZ calculations⁶ find that around a Sun-like star, an Earth-like planet could remain habitable between 0.97 and 1.70 AU. The inner edge of the HZ is set by loss of surface water and the outer edge is set by the maximum greenhouse of a CO₂ atmosphere where extensive CO₂ condensation and increased Rayleigh scattering prevent any further greenhouse warming from CO₂ (refs. ^6,7). This definition of the HZ only considers H₂O and CO₂ as greenhouse gases, so Earth-like planets warmed by other greenhouse gases (e.g., H₂ or CH₄) could remain habitable at bigger orbital distances^5,8,9. However, CH₄-rich atmospheres in the HZ may not be possible without life to generate substantial CH₄ (refs. ^10,11). In addition, more complex climate models have shown the HZ might extend to smaller orbital distances, perhaps interior to Venus’ orbit, with appropriate planetary conditions^12,13,14,15.

Residing within the HZ does not guarantee habitable surface conditions. Crucially, greenhouse gas abundances (and planetary albedo) must conspire to produce clement surface conditions. For example, by most estimates, Mars resides within the Sun’s HZ but is not habitable because there is insufficient greenhouse warming from CO₂, in part because of the lack of volcanic outgassing of CO₂. Thus, considering the planetary processes that control atmospheric CO₂ abundances on Earth-like planets in the HZ is necessary to constrain planetary habitability.

The prevailing hypothesis is that CO₂ levels are controlled by a weathering thermostat¹⁶. This can explain why Earth has maintained a clement surface throughout its history despite the  ~30% brightening of the Sun over the past  ~4.5 Gyr^{17,18,19,20,21}. The changing luminosity of the Sun with time is similar to moving a planet through the HZ, and so the same CO₂ weathering process responsible for maintaining habitability on the Earth through time, the carbonate–silicate weathering cycle, may similarly stabilize planetary climates within the HZ.

In the carbonate–silicate cycle, atmospheric CO₂ dissolves in water and weathers silicates on both the continents and seafloor, which releases cations and anions^{16,22,23,24,25,26,27}. On the continents, the weathering products, including dissolved SiO₂, HCO${\,}_{3}^{-}$, and Ca⁺⁺, wash into the oceans where the HCO${\,}_{3}^{-}$ combine with cations like Ca⁺⁺ to create CaCO₃, which precipitates out of solution. The net process converts atmospheric CO₂ into marine carbonate minerals (i.e., CaCO₃). Also, seafloor weathering occurs when seawater releases Ca⁺⁺ ions from the seafloor basalt and CaCO₃ precipitates in pores and veins. Subsequently, the carbonates within sediments and altered seafloor can be subducted.

Carbon returns to the atmosphere from outgassing. If CO₂ outgassing increases above the steady-state outgassing rate, a planet’s surface temperature rises. This leads to increased rainfall and continental weathering as well as potentially warmer deep-sea temperatures and more seafloor weathering^21,24,28. Increased weathering draws down atmospheric CO₂ and stabilizes the climate over   ~10⁶-year timescales on habitable, Earth-like planets²⁹.

One- and three-dimensional (1D, 3D) climate calculations of HZ limits^4,6,14 assume that a carbonate–silicate weathering cycle is functioning but do not explicitly include it. The assumed presence of the carbonate–silicate cycle would predict that atmospheric CO₂ of Earth-like planets increases with orbital distance in the HZ^4,6,29. In particular, future telescopic observations, e.g., NASA’s Habitable Exoplanet Imaging Mission (HabEx)³⁰ and Large Ultraviolet Optical Infrared Surveyor (LUVOIR)³¹, could search for the CO₂ trend to test the HZ hypothesis^32,33,34. Previous studies^29,35 have suggested the carbonate–silicate weathering cycle could alter predictions of pCO₂ in the HZ, but it is important to know the exact relationship we are looking for. Also, while an increase of pCO₂ with orbital distance in the HZ may be true if all Earth-like exoplanets have the exact same properties as the modern Earth, the trend becomes less certain if HZ planetary characteristics deviate from those of the modern Earth. There could be considerable variability in atmospheric CO₂ throughout the HZ, perhaps even enough to obscure a monotonic trend with orbital distance.

Here, we show that uncertain physicochemical and geophysical parameters in the carbonate–silicate weathering cycle²⁶ cause scatter in pCO₂ with orbital distance. We then demonstrate that future telescopes must observe at least 83 (2σ) HZ planetary atmospheres to confidently detect our predicted relationship between atmospheric CO₂ and orbital distance, and confirm the HZ hypothesis.

Results

Stable pCO₂ abundances from our numerical model

We use a coupled climate and carbonate–silicate weathering model (see Methods, subsection “Numerical carbonate–silicate cycle modeling”) to explore pCO₂ on Earth-like planets in the HZ. The model considers numerous planetary properties, listed in Table 1, and their effect on the carbonate–silicate weathering cycle to calculate a planet’s steady-state pCO₂ and surface temperature. If the globally averaged, steady-state surface temperature is below 248 K, we assume the planet is completely frozen and uninhabitable at the surface, as shown by 3D climate models³⁶. Similarly, we assume planets are uninhabitable beyond 355 K, above which surface water would be rapidly lost to space³⁷ (see Methods, subsection “Numerical carbonate–silicate cycle modeling” for additional details on these assumed temperature constraints).

Table 1 Parameter ranges for our numerical model.

We randomly generated 1050 habitable, stable, Earth-like exoplanet climates using uniform distributions of the model parameters in Table 1. A total of 1200 random, initial parameter combinations were considered but we eliminated those that resulted in planets that froze completely or were too hot to retain their surface oceans. As colored dots, Fig. 1 shows habitable, steady-state solutions.

Fig. 1: The expected distribution of stable, Earth-like exoplanet climates from our habitable zone weathering model.

The horizontal axis shows incident flux, S, normalized to the solar constant (S_⊕) and the corresponding orbital distance in Astronomical Units (AU) above the plot. The vertical axis shows the atmospheric CO₂ partial pressure (pCO₂) in bar. Each point represents a climate in steady state. The black labeled contours show the mean global surface temperature for the given pCO₂ and incident flux. The white region below the 248 K contour is where our model assumption of a liquid ocean is no longer plausible so no planets are shown in that region. Above the 355 K contour, Earth-like planets are too hot to retain their liquid oceans for billions of years. Similar to the frozen planets, such hot planets are not considered habitable. Modern Earth and Mars are shown by black squares. The blue histogram at the bottom of the figure shows the number of stable planets in each incident flux bin. The color of each simulated planet shows the relative point density in the plot at that location. The color was calculated using a kernel-density estimate with Gaussian kernels and rescaled from 0 to 1. A color value of 0 represents the lowest relative point density, 1 the highest. The log-linear line of best fit between pCO₂ and S is shown in red. The slope of the red, best fit line is 3.92 ± 0.24 (95%) with units −${\mathrm{log}\,}_{10}$(pCO₂ [bar])/(S/S_⊕). Our model predicts that atmospheric CO₂ should increase with orbital distance in the HZ.

Our model predicts that atmospheric CO₂ abundances should broadly increase and narrow in their spread with orbital distance in the HZ (Fig. 1), consistent with other models of CO₂ in the HZ^29,38. As justified next in section “Habitable zone climate theory revisited”, the scatter is about a nominal linear trend between incident flux, S, and log(pCO₂), which is different from a non-linear trend in models that assume a constant surface temperature in the HZ from negative feedbacks^32,34 but do not actually model the carbonate–silicate feedbacks. If future missions are to test the HZ concept by searching for a trend between incident flux, S, and pCO₂ (refs. ^32,33,34), they could search for the fundamental S–pCO₂ relationship shown in Fig. 1.

Below, we show that a log-linear relationship between pCO₂ and S may be the default in the HZ if Earth-like carbonate–silicate weathering is ubiquitous on habitable planets. In fact, the trend is elucidated by combining climate theory with carbonate–silicate cycle theory in what follows.

Habitable zone climate theory revisited

A conventional assumption is that the carbonate–silicate weathering cycle will approximately maintain a stable, temperate surface temperature for an Earth-like planet moved about in the HZ^6,7,39 or even a constant temperature^32,34. Thus, if we moved the modern Earth outward in the HZ, the smaller incident flux would initially cause the planet to cool. The cooler temperature would lower the CO₂ weathering rate causing CO₂ to accumulate in the atmosphere until the temperature returned to its nominal value of 289 K. Figure 2 shows this scenario with the dotted blue 289 K contour, which gives the pCO₂ value required to maintain a 289-K surface temperature for the modern Earth as it moves about the HZ. The line was calculated from a radiative-convective climate model described in the Methods below, subsection “Habitable zone 1D climate model” (see Eq. (8)).

Fig. 2: The relationship between incident flux and atmospheric CO₂ for Earth-like planets regulated by a carbonate–silicate weathering cycle.

The horizontal axis shows incident flux, S, normalized to the solar constant (S_⊕) and the corresponding orbital distance in Astronomical Units (AU) above the plot. The vertical axis shows the atmospheric CO₂ partial pressure (pCO₂) in bar. The dotted blue curve labeled 289 K shows the pCO₂ value required to maintain a 289 K surface temperature for the given incident flux, S. The conventional assumption of CO₂ in the HZ stipulates that pCO₂ will adjust to maintain a temperate or even constant surface temperature. Under this assumption, moving the modern Earth (labeled black square) outward in the HZ would have the planet approximately follow the dotted blue 289 K contour. The colored points and gray curves show the modern Earth moving outward in the HZ with a functioning carbonate–silicate weathering cycle, calculated from Eq. (6). We consider two temperature and pCO₂ dependencies for continental weathering in this plot. The strong temperature dependence contour (labeled Strong T-dep.), uses a temperature and pCO₂-dependent weathering factor of αT_e = 2.3, which implies a strong temperature feedback on continental weathering compared to the pCO₂ feedback (see Eq. (7)). The moderate temperature dependence contour (labeled Moderate T-dep.), uses a temperature and pCO₂-dependent weathering factor of αT_e = 7.5. These two values for αT_e result in two different paths the Earth can take as it moves outward in the HZ. The planet color shows the mean surface temperature. Log-linear fits to the colored points of the Strong T-dep. and Moderate T-dep. contours have r² values of 0.959 and 0.999, respectively. Thus, even for a strong temperature dependence of continental weathering, our coupled climate and weathering model predicts an approximately log-linear relationship between incident flux and pCO₂ on Earth-like planets in the HZ.

The constant, 289 K surface temperature contour in Fig. 2 is a non-linear relationship between incident flux, S, and log(pCO₂) but it does not consider the temperature and pCO₂ feedbacks inherent to the carbonate–silicate weathering cycle. We demonstrate that if these feedbacks are taken into account, surface temperature declines with orbital distance, as mentioned in previous work²⁹, and the relationship between S and log(pCO₂) is actually approximately linear for Earth-like planets in the HZ.

If Bond albedo is fixed, the surface temperature, T_s, for an Earth-like planet in steady state varies approximately linearly with incident flux, S^5,40,41. This linear relationship between T_s and S arises from energy balance and from water vapor feedback and can be expressed as

$${F}_{{\rm{SOL}}}={F}_{{\rm{OLR}}}=\left(\frac{1-{A}_{{\rm{B}}}}{4}\right)S=a+b{T}_{{\rm{s}}},$$

(1)

where F_SOL is the incoming solar radiation flux, F_OLR the outgoing long-wavelength radiation flux, A_B the Bond albedo, and a and b are empirical constants. From satellite measurements of the modern Earth and radiative calculations, for T_s in K, the empirical constants in Eq. (1) are approximately a = −370 W m⁻² and b = 2.2 W m⁻² K⁻¹ (ref. ⁴¹).

Solving for T_s in Eq. (1), the surface temperature is given by

$${T}_{{\rm{s}}}=\left(\frac{1-{A}_{{\rm{B}}}}{4b}\right)S-\frac{a}{b}.$$

(2)

Under the conventional assumption that the HZ is regulated by a CO₂–H₂O greenhouse effect where H₂O concentrations respond to changes in pCO₂, the temperature offset in Eq. (2),  −a/b, is a function of pCO₂. Thus, surface temperature, as a function of S and pCO₂, is given by

$${T}_{{\rm{s}}}\left(S,{{\rm{pCO}}}_{2}\right)=\left(\frac{1-{A}_{{\rm{B}}}}{4b}\right)S+f\left({{\rm{pCO}}}_{2}\right),$$

(3)

where $f\left({{\rm{pCO}}}_{2}\right)$ is a function that depends on pCO₂. For the modern Earth at 1 AU, f(pCO₂) = −a/b. For pCO₂ ≤ 0.1 bar, the CO₂ greenhouse effect is logarithmic in pCO₂, i.e., $f({{\rm{pCO}}}_{2})\propto \mathrm{ln}\,({{\rm{pCO}}}_{2})$^42,43. Above   ~0.1 bar, weaker CO₂ absorption features become important and f(pCO₂) deviates from $\propto \mathrm{ln}\,({{\rm{pCO}}}_{2})$^43,44.

As pCO₂ increases for an Earth-like planet moved outward in the HZ, the surface temperature will follow Eq. (3) while the rate at which CO₂ is removed from the atmosphere will adjust according to the carbonate–silicate weathering feedback. Quantitatively, the pCO₂- and T_s-dependent flux of CO₂ removal due to the continental weathering flux, F_w (in mol CO₂ per unit time) is described by

$${F}_{{\rm{w}}}=\rho {\left(\frac{{{\rm{pCO}}}_{2}}{{{\rm{pCO}}}_{2}^{{\rm{mod}}}}\right)}^{\alpha }\exp \left(\frac{{T}_{{\rm{s}}}\left(S,{{\rm{pCO}}}_{2}\right)-{T}_{{\rm{s}}}^{{\rm{mod}}}}{{T}_{{\rm{e}}}}\right),$$

(4)

where ρ is a constant determined by the continental weathering properties of the modern Earth, α a dimensionless constant between 0.1 and 0.5 and regulates the pCO₂ dependence of continental silicate weathering, T_e a constant between 10 K and 40 K and represents the e-folding temperature dependence of continental weathering. The range for T_e has been empirically constrained for the surface temperatures relevant to habitable, Earth-like planets from lab measurements and Phanerozoic geologic constraints^16,26. Finally, pCO${\,}_{2}^{{\rm{mod}}}=288\times 1{0}^{-6}$ bar and ${T}_{{\rm{s}}}^{{\rm{mod}}}=289$ K are the modern Earth’s preindustrial pCO₂ and surface temperature, respectively²¹.

The range for α on the Earth was determined empirically from geologic constraints over the past 100 Myr²⁶. We assume that this derived range for α applies to the Earth through time^21,45,46 and the Earth-like exoplanets modeled here that have a carbonate–silicate cycle. However, better proxy data for the ancient Earth or observing the carbonate–silicate cycle on habitable exoplanets^32,34 may be necessary to understand if the assumed range for α applies more generally to habitable planets.

In Eq. (4), we assume seafloor weathering is negligible, which is a reasonable approximation for the modern Earth²¹, and illustrative for our purposes of deriving a simple, analytic relationship between S and pCO₂. Here, we seek to predict the behavior of the modern Earth in the HZ to gain intuitive understanding, whereas in our numerical model we consider a broad range of properties for Earth-like planets on which seafloor weathering may be important.

The modern Earth, and all Earth-like planets considered in this work, are assumed to be in steady state, in which the flux of CO₂ from volcanic outgassing is equal to the rate of CO₂ removal from weathering, F_w. If we assume a HZ planet with CO₂ outgassing the same as the modern Earth’s, F_w remains constant despite changes in S and pCO₂. Setting ${T}_{{\rm{s}}}\left(S,{{\rm{pCO}}}_{2}\right)={T}_{{\rm{s}}}^{{\rm{mod}}}$ and ${{\rm{pCO}}}_{2}={{\rm{pCO}}}_{2}^{{\rm{mod}}}$ for the modern Earth, from Eq. (4), we see that F_w = ρ and

$$1={\left(\frac{{{\rm{pCO}}}_{2}}{{{\rm{pCO}}}_{2}^{{\rm{mod}}}}\right)}^{\alpha }\exp \left(\frac{{T}_{{\rm{s}}}\left(S,{{\rm{pCO}}}_{2}\right)-{T}_{{\rm{s}}}^{{\rm{mod}}}}{{T}_{{\rm{e}}}}\right).$$

(5)

Equation (5) must hold for a modern Earth within the HZ. If it did not, F_w would not balance CO₂ outgassing, which would result in either complete CO₂ removal, or CO₂ accumulation without bound.

Solving for ${T}_{{\rm{s}}}\left(S,{{\rm{pCO}}}_{2}\right)$ in Eq. (5), we find

$${T}_{{\rm{s}}}\left(S,{{\rm{pCO}}}_{2}\right)={T}_{{\rm{s}}}^{{\rm{mod}}}-\alpha {T}_{{\rm{e}}}\mathrm{ln}\,\left(\frac{{{\rm{pCO}}}_{2}}{{{\rm{pCO}}}_{2}^{{\rm{mod}}}}\right).$$

(6)

Equating Eq. (6) to Eq. (3) and rearranging gives

$$\left(\frac{1-{A}_{{\rm{B}}}}{4b}\right)S={T}_{{\rm{s}}}^{{\rm{mod}}}-\left[\alpha {T}_{{\rm{e}}}{\mathrm{ln}}\,\left(\frac{{{\rm{pCO}}}_{2}}{{{\rm{pCO}}}_{2}^{{\rm{mod}}}}\right)+f\left({{\rm{pCO}}}_{2}\right)\right].$$

(7)

If $f({{\rm{pCO}}}_{2})\propto \mathrm{ln}\,({{\rm{pCO}}}_{2})$, which is the case for pCO₂ ≤ 0.1 bar^43,44, then $S\propto -\mathrm{ln}\,({{\rm{pCO}}}_{2})$. However, even if f(pCO₂) deviates from log-linearity with pCO₂, S will become increasingly log-linear with pCO₂ as αT_e increases. In Eq. (7), increasing αT_e will cause the $\mathrm{ln}\,({{\rm{pCO}}}_{2})$ term to dominate the relationship between S and pCO₂. Intuitively, increasing αT_e decreases the temperature dependence of continental weathering relative to its pCO₂ dependence. Note that bigger T_e reduces the temperature dependence of continental weathering while bigger α increases the pCO₂ dependence of continental weathering (Eq. (4)).

In addition to predicting a linear relationship between ${\mathrm{log}}\left({\mathrm{pCO}}_{2}\right)$ and S, the carbonate–silicate cycle implies that moving an Earth-like planet outward in the HZ will cause ${T}_{{\rm{s}}}\left(S,{{\rm{pCO}}}_{2}\right)$ to decrease. For increasing orbital distance, pCO₂ must increase for ${T}_{{\rm{s}}}\left(S,{{\rm{pCO}}}_{2}\right)$ to increase in the HZ. From Eq. (5), pCO₂ will be greater than pCO${\,}_{2}^{{\rm{mod}}}$ in such cases so ${T}_{{\rm{s}}}\left(S,{{\rm{pCO}}}_{2}\right)$ must be less than ${T}_{{\rm{s}}}^{{\rm{mod}}}$. This decrease in ${T}_{{\rm{s}}}\left(S,{{\rm{pCO}}}_{2}\right)$ as S decreases is shown in Fig. 2. Physically, the power law dependence of weathering on pCO₂ can balance volcanic outgassing at lower surface temperatures in the outer HZ.

 Figure 2 shows the approximately log-linear relationship between pCO₂ and S for the modern Earth moved outward in the HZ. The gray lines and colored circles in Fig. 2 show the expected pCO₂ value for the given incident flux S, calculated from Eq. (6). For each S value in Fig. 2, Eq. (6) was solved for pCO₂ by using Eq. (8), the polynomial fit for surface temperature based on a 1D climate model (described in the Methods, subsection “Habitable zone 1D climate model”), assuming values of αT_e.

The value of αT_e affects the slope of the relationship between S and pCO₂ due to the carbonate–silicate weathering cycle, shown in Fig. 2. From above, the ranges for α and T_e are 0.1 ≤ α ≤ 0.5 and 10 ≤ T_e ≤ 40 (ref. ²¹), so 1 ≤ αT_e ≤ 20. If we consider uniform distributions of α and T_e, then roughly 95% of αT_e values will be greater than 2.3. If α = 0.23 and T_e = 10 K then αT_e =  2.3, which is used for the Strong T-dep. curve in Fig. 2. The mean of each parameter, α = 0.3 and T_e = 25 K gives αT_e = 7.5, which corresponds to the Moderate T-dep. curve in Fig. 2. For αT_e ≤ 2.3 the colored points and gray curves become increasingly similar to the constant 289 K surface temperature contour in Fig. 2. However, for uniform distributions of α and T_e,   ~95% of αT_e values are greater than 2.3, so an approximately log-linear relationship between S and $\mathrm{log}\,({{\rm{pCO}}}_{2})$ is the default expectation for Earth-like planets in the HZ.

Observational uncertainty for pCO₂ in the HZ

In the log-linear fit shown as the solid red line in Fig. 1, which is the expected relationship between pCO₂ and S that we have derived above, the r²-value is 0.49. Thus, about half the variance in log(pCO₂) is described by changes in incident flux. The slope is 3.92 ± 0.24 (95%) with units −${\mathrm{log}\,}_{10}$(pCO₂ [bar])/(S/S_⊕), so our model predicts a trend of increasing atmospheric CO₂ with orbital distance, which future missions might detect^32,33,34. However, there is sufficient spread in our simulated planets that confirming the HZ hypothesis from such a trend may be difficult.

This difficulty is readily seen if we consider a log-uniform distribution for pCO₂ on Earth-like planets in the HZ. If we randomly generate 1050 such planets, where 10⁻⁴ ≤ pCO₂ ≤ 10 bar is sampled log-uniformly, 0.35 ≤ S ≤ 1.05 is sampled uniformly, and impose the same constraints on surface temperature for habitability as in Fig. 1, then the log-linear line of best fit through the uniform planet data has a slope of 3.76 ± 0.465 (95%) with units −${\mathrm{log}\,}_{10}$(pCO₂ [bar])/(S/S_⊕). Thus, measuring just the log-linear trend between pCO₂ and S in the HZ is unlikely to test the HZ hypothesis as this measurement cannot confidently detect the presence of the carbonate–silicate weathering cycle—it is indistinguishable from that of randomly distributed pCO₂ between the surface temperature limits for habitability.

The inability to differentiate between the log-linear trends for weathering-mediated and random pCO₂ vs S in the HZ is due to the assumed surface temperature constraints we impose in our model (between 248 and 355 K for planets in the HZ, see Methods, subsection “Numerical carbonate–silicate cycle modeling”). Such temperature constraints are necessary as the carbonate–silicate weathering cycle can only operate when water, as liquid, is present at the planetary surface. Even without the carbonate–silicate weathering cycle, a minimum surface temperature for habitable planets, which must exist, will result in increasing pCO₂ with orbital distance, as shown by the constant temperature contours in Fig. 1.

To test the HZ hypothesis, we propose searching for the two-dimensional (2D) distribution of planets in the S-pCO₂ phase space that arises from the carbonate–silicate weathering cycle. This S-log(pCO₂) relationship is shown by the point density in Fig. 1, where the distribution of habitable, stable planets is not log-uniformly distributed over pCO₂. Rather, around the best-fit line, there is an abundance of planets in the outer HZ at high pCO₂, a dearth of low pCO₂ planets between   ~0.9 and   ~0.7 S/S_⊕, and few high-pCO₂ planets throughout the HZ compared to the log-uniform pCO₂ case. These differences are expected features of the carbonate–silicate weathering cycle due to the temperature- and pCO₂-dependent nature of the weathering feedback. Recall from section “Habitable zone climate theory revisited”  that, as S decreases, the lowered temperature will reduce weathering causing pCO₂ to increase. This results in the lack of low-pCO₂ planets in the middle of the HZ and the high abundance of habitable planets in the outer HZ (purple shaded region in Fig. 1). Similarly, for large pCO₂, the temperature is warmer and pCO₂ higher than that of modern Earth so the carbonate–silicate weathering cycle acts to lower pCO₂, which reduces the number of high-pCO₂ planets throughout the HZ relative to the outer HZ.

To detect the prevalence of the carbonate–silicate weathering cycle and test the validity of the HZ concept, future observations should measure the 2D S-pCO₂ distribution of habitable, Earth-like exoplanets. This measured distribution can be compared to the distribution of habitable planets we predict in Fig. 1 to determine if Earth-like planets in the HZ are consistent with the S-pCO₂ predictions of the carbonate–silicate weathering cycle.

A test of the 2D phase space of S and pCO₂ in the HZ is shown in Fig. 3, which was produced using a 2D Kolmogorov–Smirnov (KS) test. The 2D KS test compares the statistical similarity of a sample distribution to a reference distribution^47,48,49. For Fig. 3, the reference distribution was comprised of 500 randomly generated planets from the log-uniform distribution for pCO₂ described above (10⁻⁴ ≤ pCO₂ ≤ 10 bar, 0.35 ≤ S ≤ 1.05S/S_⊕, and surface temperature between 248 and 355 K). The sample distribution was generated by randomly selecting a number of planets from Fig. 1 equal to the number of observed exoplanets. For a given number of observed exoplanets in Fig. 3, the horizontal axis, we ran the KS test 10,000 times then calculated the mean and standard deviation from those runs, shown by the gray contour and shaded region. This resampling is necessary as the 2D KS test is a nonparametric approximation that two data sets come from the same underlying population⁴⁸. We note that below  ~20 observed planets and for probabilities above   ~0.1, the 2D KS test used here can be unreliable⁴⁹. These limitations do not invalidate the analysis shown in Fig. 3, as we want to know, with 95% confidence, that a log-uniform pCO₂ distribution can be ruled out if real exoplanets follow the distribution shown in Fig. 1, which corresponds to the gray line and shaded contour dipping below the 0.05 probability value, shown by the horizontal black line, at 83 observations in Fig. 3.

Fig. 3: The probability observed exoplanets will accidentally match a log-uniform distribution for pCO₂ in the HZ if the true pCO₂ distribution is regulated by the carbonate–silicate weathering cycle, as shown in Fig. 1.

This probability is shown on the vertical axis. The horizontal axis shows the number of observed Earth-like exoplanets. The solid gray curve and corresponding shaded gray region show the mean probability and 2σ uncertainty, respectively, that the observed planets, sampled from the planets shown in Fig. 1, match a log-uniform pCO₂ distribution in the HZ. This curve is calculated from 10,000 two-dimensional Kolmogorov–Smirnov tests (see Results, subsection “Observational uncertainty for pCO₂ in the HZ”). The solid, horizontal black line highlights the 5% probability line. At the bottom of the figure, the labeled points and error bars show the number of Earth-like exoplanets the next generation of proposed space telescopes are expected to observe (telescope data in Table 2). The vertical scaling of the telescope points is arbitrary, only the horizontal position and extent of the 1σ error bars is significant. To rule out a log-uniform pCO₂ distribution with 95% confidence, future telescopes would need to observe at least 83 Earth-like planets.

Thus, confidently detecting the carbonate–silicate weathering cycle will require many exoplanet observations, as shown in Fig. 3. Proposed NASA telescopes, HabEx and LUVOIR, are expected to observe between 3 and 115 Earth-like exoplanets^30,31 (see Table 2). The ranges for each mission concept are shown by the colored circles with error bars in Fig. 3. Only the nominal capability of LUVOIR-A, the variant of the proposed LUVOIR space telescope with a primary mirror diameter of 15 m, would provide sufficient Earth-like exoplanet detections to confidently discriminate between a log-uniform pCO₂ distribution in the HZ and a pCO₂ distribution regulated by the carbonate–silicate weathering cycle. A caveat is that this calculation does not consider the instrument uncertainty in derived pCO₂ measurements for each telescope or that other processes not considered in our model may alter pCO₂ in the HZ, as discussed below.

Table 2 The number of expected Earth-like exoplanets observed by each platform from the HabEx³⁰ and LUVOIR³¹ final reports.

Discussion

Our model assumes that the full variation and uncertainty in Earth’s carbon cycle parameters through time (Table 1) are representative of habitable Earth-like exoplanets generally. This assumption is a reasonable first-order approximation as the bulk composition and geochemistry of rocky exoplanets appear similar to Earth’s⁵⁰. However, the validity of this assumption likely depends on the parameter in question. For example, it is probably reasonable to expect habitable exoplanets to have a wide range of land fractions and outgassing fluxes, but it is unclear whether there is as much natural variability in the temperature dependence of silicate weathering. An improved mechanistic understanding of weathering on Earth^51,52 might reduce these uncertainties.

Other weathering feedbacks have been proposed to operate on the Earth through time, such as reverse weathering⁵³. In reverse weathering, cations and dissolved silica released from silicate weathering are sequestered into clay minerals rather than carbonates so that CO₂ remains in the atmosphere, warming the climate and reducing ocean pH. Reverse weathering is thought to be strongly pH dependent and as ocean pH decreases, reverse weathering turns off, acting as a climate stabilization mechanism similar to the carbonate–silicate cycle. The importance of reverse weathering is so poorly constrained through Earth’s history⁴⁵ that it does not make sense to consider it in our model. However, with future constraints from geology and lab measurements, reverse weathering might alter the stable CO₂ abundances of our modeled atmospheres shown in Fig. 1.

At both the inner and outer edges of the HZ, our model assumes that abundant liquid water exists at the planetary surface because, without a liquid surface ocean, the carbonate–silicate weathering cycle ceases and CO₂ cannot be sequestered after outgassing. Beyond these temperature bounds, other processes must regulate pCO₂. This is a caveat to consider in future observations. As we see from Fig. 1, Mars has low atmospheric CO₂ and low incident flux. Frozen exoplanets similar to Mars, populating the white area under the 248 K contour in Fig. 1, could exist in exoplanet surveys. Similarly, planets devoid of surface water, such as Venus, might exist at high pCO₂ within the HZ. If future observations detect such planets without confirming the existence of a liquid surface or surface temperature, it could introduce additional uncertainties in any relation between orbital distance and atmospheric CO₂. Detecting a surface ocean, one of the most important surface features to confirm when searching for biosignatures and habitability^54,55,56, is also important to interpret trends of CO₂ in the HZ.

Because we only consider variations on an Earth-like planet, our model predictions may underestimate the inherent variability in habitable exoplanetary conditions. Planets very different from the modern Earth, such as waterworlds without a carbonate–silicate weathering cycle⁵⁷ or CH₄-rich worlds^58,59, could introduce additional uncertainty in an observed relationship between S and pCO₂ in the HZ. Despite such uncertainties, future missions should measure the relationship between S and pCO₂ in the HZ, or possibly a sharp transition in pCO₂ at the inner edge of the HZ due to loss of surface water and subsequent shutoff of surface weathering^38,60. A more complex model than presented here is necessary to predict such a jump in pCO₂ at the inner edge of the HZ. However, if the carbonate–silicate weathering cycle is indeed ubiquitous, as is typically assumed in HZ calculations, then the relationship between incident flux and pCO₂ may follow the S-pCO₂ relationship predicted in Fig. 1. If no such relationship is observed, then the carbonate–silicate weathering cycle may have limited influence on planetary habitability and the limits of the conventional HZ could need revision. Alternatively, the HZ hypothesis could be incorrect and the long-term climate of HZ planets could be set by phenomena beyond those considered here.

A previous version of this work was published as part of a Ph.D. thesis⁶¹.

Methods

Habitable zone 1D climate model

We use the Virtual Planetary Laboratory (VPL) 1D radiative-convective climate model^5,62 to generate surface temperatures for an Earth-like planet at various pCO₂ and incident fluxes. We consider incident fluxes between 1.05S_⊕ and 0.35S_⊕, the HZ limits for a Sun-like star⁶, and atmospheric CO₂ partial pressures between 10⁻⁶ and 10 bar. We assume the atmosphere is comprised of CO₂ and H₂O. If the CO₂ partial pressure is below 1 bar, we set the initial atmospheric pressure to 1 bar and add N₂ to the atmosphere such that the total surface pressure is 1 bar. We fix the stratospheric water vapor concentration to the modern Earth value and follow the Manabe–Wetherald relative humidity distribution in the troposphere with empirical constraints based on the modern Earth⁶³.

We fit the surface temperature output, T_s in K, from the climate model with a fourth-order polynomial in $\mathrm{ln}\,({{\rm{pCO}}}_{2})$ and normalized stellar flux, as follows:

$$\begin{array}{l}{T}_{{\rm{s}}}\left(S,{{\rm{pCO}}}_{2}\right)=4.809-222.0X-68.44{X}^{2}-6.737{X}^{3}-0.206{X}^{4}\\ +1414XY+446.4{X}^{2}Y+44.41{X}^{3}Y+1.364{X}^{4}Y\\ -2964X{Y}^{2}-978.4{X}^{2}{Y}^{2}-98.86{X}^{3}{Y}^{2}-3.059{X}^{4}{Y}^{2}\\ +2655X{Y}^{3}+907.5{X}^{2}{Y}^{3}+92.87{X}^{3}{Y}^{3}+2.892{X}^{4}{Y}^{3}\\ -868.4X{Y}^{4}-304.6{X}^{2}{Y}^{4}-31.48{X}^{3}{Y}^{4}-0.985{X}^{4}{Y}^{4}\\ +1045Y-1496{Y}^{2}+1064{Y}^{3}-281.1{Y}^{4}.\end{array}$$

(8)

Here, CO₂ partial pressure pCO₂ is in bar, $X=\mathrm{ln}\,({{\rm{pCO}}}_{2})$, and Y = S/S_⊕ is the incident flux, S, normalized to the solar constant, S_⊕. Figure 4 shows the agreement between the 1D climate model and the polynomial fit used in this work.

Fig. 4: The relative error between our fourth-order polynomial fit and the full 1D radiative-convective climate model.

Our polynomial fit is valid between 1.05S_⊕ and 0.35S_⊕, where S_⊕ is the solar constant. The polynomial fit is valid for atmospheric CO₂ abundances between 10⁻⁶ and 10 bar. The surface temperatures predicted by the polynomial fit reproduce the results of the 1D climate model. The maximum error in predicted surface temperature between the polynomial fit and the 1D climate model is  ~3%.

Numerical carbonate–silicate cycle modeling

To calculate the steady-state pCO₂ in the atmospheres of Earth-like planets in the HZ, we use a weathering model that describes pCO₂ on the Earth through time^21,26. We summarize the model below and highlight how the model in this work differs from previous implementations^21,26. These previous implementations provide a comprehensive explanation and justification of the model parameterizations, and empirical and theoretical basis. The model, as a Python script, is available in the Supplementary Data and contains a complete description of the model equations and parameters (see the file weathering_model.py).

The weathering model balances the flux of outgassed CO₂ against the loss of carbon due to continental and seafloor weathering, which result in precipitation of carbonates in the ocean and seafloor pore space. Quantitatively, for time t, this is described by time-dependent equations where we normalize to the mass of the ocean, M_o (nominally, an Earth ocean, 1.35 × 10²¹ kg):

$$ \frac{{\mathrm{d}}C}{{\mathrm{d}}t} =\frac{{F}_{{\rm{out}}}\,+\,{F}_{{\rm{carb}}}\,-\,{P}_{{\rm{o}}}\,-\,{P}_{{\rm{p}}}}{{M}_{{\rm{o}}}}\\ \frac{{\mathrm{d}}A}{{\mathrm{d}}t} =2\times \frac{{F}_{{\rm{sil}}}\,+\,{F}_{{\rm{carb}}}\,+\,{F}_{{\rm{diss}}}\,-\,{P}_{{\rm{o}}}\,-\,{P}_{{\rm{p}}}}{{M}_{{\rm{o}}}}.$$

(9)

Here, C is the non-organic carbon content of the atmosphere–ocean system in mol C kg⁻¹, and A is the carbonate alkalinity in mol equivalents (mol eq). Carbonate alkalinity (henceforth alkalinity) is the charge-weighted sum of the mol liter⁻¹ concentration of bicarbonate and carbonate anions, [HCO${\,}_{3}^{-}$] + 2[CO${\,}_{3}^{2-}$]. F_out is the global CO₂ outgassing flux, F_carb and F_sil are the continental carbonate and silicate weathering fluxes, F_diss is the rate of seafloor basalt dissolution, and P_p and P_o are the pore and ocean precipitation fluxes. The fluxes on the right-hand side of Eq. (9) (F_out, F_carb, P_o, P_p, F_sil, F_carb, and F_diss) are given in mol C yr⁻¹ for dC/dt and in mol eq yr⁻¹ for dA/dt.

The alkalinity that enters the ocean from weathering will balance a +2 charge cation (e.g., Ca⁺⁺), which is why a factor of 2 enters in the definition of dA/dt in Eq. (9). Hence, geochemists often think of alkalinity in terms of the balance of cations produced in weathering, principally Ca⁺⁺. This reasoning arises because the weighted sum of carbonate and bicarbonate concentrations must balance the charge of conservative cations minus conservative anions $\left({\mathrm{i}}.{\mathrm{e}}., \, 2[{\mathrm{Ca}}{^{++}}]\,+\,2[{\mathrm{Mg}}{^{++}}]\,+\,{\mathrm{Na}}^{+}+\ldots{\hbox{-}}[{\mathrm{Cl}}^{-}]{\hbox{-}}{\ldots}{\!}\right)$, ignoring minor contributions from weak acid anions and water dissociation. Weathering releases cations and carbon speciation adjusts to ensure charge balance, so that the cation release is effectively equivalent to carbonate alkalinity.

To improve the rate of model convergence and range of model inputs over which Eq. (9) converges, we do not consider the seafloor pore space and atmosphere–ocean as separate systems. This differs from previous versions of the model²¹, which considered the atmosphere–ocean and pore space independently. Rather, we approximate the atmosphere–ocean and pore space as a single entity in Eq. (9). This simplification does not appreciably change the model output for atmospheric CO₂ because we run the model to steady state in all cases, where the atmosphere–ocean and pore space reach approximate equilibrium. In the next section, we present additional details on our model implementation and discuss the agreement between our no-pore model and the original, two-box model²¹.

A second modification is the range of incident stellar fluxes over which the model can be run. Previously, the model described here was used to study the Earth through time²¹ and thus only considered solar fluxes between S_⊕ (the modern solar constant) and early Earth’s 0.7S_⊕ (S_⊕ = 1360 W m⁻²). We extend that range to include the entire conservative HZ of a Sun-like star, roughly 1.05S_⊕ to 0.35S_⊕⁶. We use Eq. (8), the fourth-order polynomial fit to a 1D climate model, to calculate surface temperatures throughout the HZ. The Bond albedo of the planet is calculated dynamically by the climate model and thus included implicitly in our polynomial fit.

With the coupled climate and weathering model, we generate steady-state, Earth-like planets by randomly sampling plausible initial model inputs. The ranges for each parameter we consider are representative of the Earth through time²¹ and shown in Table 1. These ranges represent very broad uncertainties of the carbonate–silicate cycle on the Earth through time and so are appropriate for Earth-like planets. We conservatively assume a uniform distribution for each parameter range shown in Table 1.

We parameterize the internal heat of an Earth-like planet conservatively using the planet’s age, ranging 0–10 Gyr, which is the approximate habitable lifetime of an Earth-like planet around a Sun-like star⁶⁴. The equation for planetary heat relative to the modern Earth, Q, is given by

$$Q={\left(1-\frac{4.5-\tau }{4.5}\right)}^{-{n}_{{\rm{out}}}},$$

(10)

where τ is the age of the planet in Gyr, and n_out is the scaling exponent for internal heat, with a range given in Table 1.

The parameter ranges shown in Table 1 represent the uncertainty of the carbonate–silicate weathering cycle on the Earth through time²¹. Implicit in our assumed parameter ranges is that continental land fraction, f_land, and biological weathering fraction, f_bio, have increased from 0 when the Earth formed to 1 on the modern Earth. Similarly, the relative internal heat, Q, is assumed to be large when the Earth is young and unity for the modern Earth. Therefore, on the modern Earth, where f_land = 1, f_bio = 1, and Q = 1, the weathering rate is maximized and outgassing rate is relatively small (see Methods, subsection “Validity of carbon cycle parameterizations to exoplanets” for a discussion on the importance of these three parameters in our model). This is seen in Fig. 1, where the modern Earth appears near the lower bound for predicted pCO₂ in the HZ. If the continents on an exoplanet were more easily weathered or outgassing much lower than on the modern Earth, such exoplanets could have pCO₂ values well below the modern Earth value shown in Fig. 1. We do not consider such exoplanets in this model, so the results presented here are only applicable to planets similar to the Earth through time.

Our model assumes that each simulated planet is habitable, i.e., it has a stable, liquid surface ocean, a necessity for the carbonate–silicate cycle to operate. For a mean surface temperature below 248 K, Earth-like planets are likely completely frozen³⁶, which we use as a lower temperature bound in the model. While 248 K is below the freezing point of water, it is a global mean surface temperature and 3D models show that the range 248–273 K for this parameter does not preclude the existence of a liquid ocean belt near the equator. At the other temperature extreme, a hot, Earth-like planet can rapidly lose its surface oceans due to high atmospheric water vapor concentrations that are photolyzed and subsequently lost to space. This upper temperature bound on habitability occurs at   ~355 K³⁷. Above 355 K, Earth-like planets are unlikely to remain habitable for more than   ~1 Gyr³⁷ and cannot operate a carbonate–silicate cycle over geologic timescales. We use these two temperature bounds, 248 K and 355 K, as the limits for habitability in our model. Any modeled planet with a final surface temperature outside these limits is uninhabitable and removed from our results.

We limit HZ planets to those with pCO₂ below 10 bar. For most Earth-like planets in the HZ, 10 bar of CO₂ results in planets with surface temperatures well above 355 K, which are not habitable on long time scales. If we impose a fixed stratospheric water vapor concentration in the 1D climate model and modify the tropospheric water vapor concentration based on empirical data from the modern Earth, we enable the 1D climate model to accurately model habitable, Earth-like planets through much of the HZ. But in the outer HZ, with more than  ~10 bar of CO₂, this assumption overestimates atmospheric water vapor concentrations and leads to artificially warm planets, so we reject such cases. Above   ~10 bar of CO₂ in the outer HZ, assuming a saturated troposphere for water vapor, increasing atmospheric CO₂ may not lead to additional warming⁶. Rather, the surface cools in such scenarios because additional CO₂ leads to increased Rayleigh scattering and no additional warming. Because Earth-like planets in the outer HZ would be frozen and uninhabitable even with CO₂ partial pressures above  ~10 bar, we impose a 10-bar limit for CO₂ in the outer HZ. This limit agrees with previous CO₂ limitations in coupled climate and weathering models²⁹.

Combined ocean and pore space model justification

The carbon cycle model used in this work was previously derived as a two-box model²¹, where the atmosphere–ocean and the seafloor pore space were separated. In this work, we combine the ocean–atmosphere and the pore space into a single unit. This modification can be implemented in the original model²¹ by assuming that the pH of the pore space is the same as the pH of the ocean, and assuming that the alkalinity and carbon content of the ocean and pore space are the same. The dissolution and precipitation fluxes can then be calculated without treating the ocean–atmosphere and the pore space as different systems. This modification allows the model to converge quicker over a wider range of parameter combinations.

To validate our combined model, we ran the modern Earth through both the original, two-box model²¹ and our modified model at ten different incident fluxes between S_⊕ and 0.7S_⊕. The average error in predicted CO₂ values between our model and the two-box model was 2.8%, with a minimum error of 2.3%, and a maximum error of 3.6%. Given the large uncertainties in model inputs (Table 1), the few percent error introduced by our simplified model is unimportant.

For each parameter combination in our simplified model, we start with the modern Earth then impose a step change for each model parameter. We then run the simulation for 10 Gyr or until the system reaches steady state. We consider the model to have reached steady state when extrapolation of the rate of change of pCO₂ for 1 Gyr changes pCO₂ by <1%. Typically, the model converges within a few Myr to a few tens of Myr. Rarely (2 of the 1200 planets simulated in this work), parameter combinations will not reach steady state after 10 Gyr. Simulations with combinations of exceptionally high outgassing rates and low CO₂ weathering rates can enter a regime were atmospheric CO₂ builds without bound, never converging. Such model results are beyond the range of validity of our model.

Validity of carbon cycle parameterizations to exoplanets

The parameterization of weathering in our model has been empirically validated for the modern Earth^16,21,26. The exponential temperature dependence of continental weathering is a reasonable approximation that agrees with field and lab measurements¹⁶ and can reproduce the climate results of more complex models^21,26. Similarly, the power-law parameterization for the pCO₂ dependence of continental weathering agrees with data from the modern Earth²⁶ and can even be approximately derived from equilibrium chemistry arguments for an Earth-like exoplanet³⁸. The bulk geochemistry of rocky exoplanets may be similar to Earth’s⁵⁰, so we expect our weathering parameterization to reasonably approximate Earth-like planets in the HZ. However, uncertainties in how the carbonate–silicate weathering cycle regulates climate on Earth persist²¹, so the predicted variations in pCO₂ in our model may not capture the true variability of pCO₂ in the HZ. Below, we show that our broad parameterization of the carbonate–silicate weathering cycle may encompass the plausible range of pCO₂ for the Earth through time, but improved understanding the carbonate–silicate weathering cycle may be necessary to know if such variations are indeed representative of the Earth through time and applicable to Earth-like planets generally.

The rate of weathering depends strongly on the intrinsic features of a planet, such as the CO₂ outgassing rate and the properties of its continents. Changes in continental uplift rate, lithology, and configuration are parameterized in our model through the f_land and f_bio terms. The parameters f_land and f_bio linearly scale the weathering flux and could analogously be considered a continental weatherability scaling factor. For the ranges of f_land and f_bio considered in our model (see Table 1), changes in the continental weatherability alone can generate pCO₂ values spanning  ~4 orders of magnitude. This broad parameterization likely encompasses pCO₂ perturbations due to continental weatherability changes caused by large volcanic eruptions or changes in continental configuration. Indeed, the largest, constrained change in pCO₂ due to such events on Earth may be closer to  ~1 order of magnitude, coeval with the eruption of the Siberian Traps⁶⁵.

The importance of continental weatherability (f_land and f_bio) on pCO₂, relative to other parameters, is shown in Fig. 5. Figure 5 was generated by sampling uniform distributions for each model parameter shown in Table 1 across its listed range. When one parameter was varied, all other parameters were held constant at their modern Earth value, which we define as: ${F}_{{\rm{mod}}}^{{\rm{out}}}=6$ Tmol C yr⁻¹, n = 1.75, x = 1, T_e = 25 K, α = 0.3, ξ = 0.3, f_land = 1, S_thick = 1, ${F}_{{\rm{carb}}}^{{\rm{mod}}}=10$ Tmol C yr⁻¹, f_bio = 1, ${a}_{{\rm{grad}}}=1.075$, γ = 0.2, β = 0.1, m = 1.5, E_bas = 90 kJ mol⁻¹, and Q = 1. Note that we incorporate n_out and τ from Table 1 into Q, the internal heat (see Eq. (10)), which is the parameter of interest. We show two different values for S in Fig. 5, S = 1.0S_⊕ in the left panel and S = 0.5S_⊕ in the right panel. For both values of S, Fig. 5 shows that variations in f_land and f_bio alone can alter pCO₂ by orders of magnitude.

Fig. 5: The spread in steady-state pCO₂ from varying a single model parameter.

Each parameter in Table 1 is shown on the vertical axis (note that n_out and τ are incorporated into Q, see Eq. (10)). The left panel shows an incident flux of S = 1.0S_⊕. The right panel shows an incident flux of S = 0.5S_⊕. For each parameter, we held all other parameters constant at the modern Earth value (see text) and randomly sampled 100 values for the parameter in question from uniform distributions of the ranges given in Table 1. The horizontal extent of the gray shaded region shows the range of possible pCO₂ values when all other parameters are fixed. The thickness of each gray shaded region shows the relative abundance of steady-state planets at that pCO₂. The thickest regions show maximum relative abundance, the thinnest regions show the minimum relative abundance. No surface temperature limits on habitability were imposed for the simulated planets. At low pCO₂, three parameters (f_land, f_bio, and Q) dominate the spread in pCO₂. At higher pCO₂, the temperature- and pCO₂-dependence of continental silicate weathering (T_e and α) are also important.

The internal heat of the planet, Q, plays a similarly important role in setting pCO₂. The rate of CO₂ outgassing is determined by Q and our broad parameterization of Q allows pCO₂ to vary by orders of magnitude throughout the HZ, as shown in Fig. 5.

The rate of CO₂ outgassing and continental weatherability drive the majority of the spread in pCO₂ shown in Fig. 1. This is readily seen in Fig. 6, which shows the results of 300 random parameter combinations from uniform distributions of the parameters in Table 1 except for Q, f_land, and f_bio, which were all fixed to 1. Of the 300 parameter combinations, 235 remained above 248 K and are shown in Fig. 6. Comparing Fig. 6 to Fig. 1, it is readily apparent that the broad uncertainty in pCO₂ from our results is due to variations in intrinsic planetary properties (Q, f_land, and f_bio) rather than uncertainties in the tuning parameters of our carbon cycle parameterization.

Fig. 6: The expected distribution of stable, steady-state pCO₂ on Earth-like planets if Q, f_land, and f_bio are fixed to 1, i.e., modern Earth values.

Except for fixing Q = f_land = f_bio = 1, this figure is generated identically to Fig. 1. The spread in pCO₂ in the outer HZ is due to the temperature- and pCO₂-dependence of continental weathering (T_e and α). This is expected from Eq. (4), which shows that T_e and α will be increasingly influential as pCO₂ and surface temperature deviate from the modern Earth values, as discussed in the Results, subsection “Habitable zone climate theory revisited”. Without changes in Q, f_land, and f_bio, there is little spread in pCO₂ in the inner HZ.

The outgassing rate and continental properties of habitable exoplanets remain unknown. Thus, our broad parameterization of those terms, which align with possible conditions on Earth throughout its history, are a reasonable approximation. If an Earth-like, carbonate–silicate weathering cycle is common on habitable planets, then these parameters may largely determine pCO₂ on such planets and generate a range for pCO₂ at a given orbital distance similar to that shown in Fig. 1.