CO2 condensation is a serious limit to the deglaciation of Earth-like planets

TL;DR

This study uses a 3D climate model to show that CO2 condensation on Earth-like planets can prevent deglaciation by trapping CO2 in polar caps and under water ice, especially beyond 1.27 AU from a Sun-like star.

Contribution

It demonstrates the limiting effect of CO2 condensation on the ability of Earth-like planets to escape from glaciation, highlighting the importance of polar CO2 ice and water ice cover in climate stability.

Findings

Planets beyond 1.27 AU may never deglaciate due to CO2 condensation.

Polar CO2 ice caps can trap significant amounts of atmospheric CO2.

Water ice cover can sequester CO2, further inhibiting deglaciation.

Abstract

It is widely believed that the carbonate-silicate cycle is the main agent to trigger deglaciations by CO$_2$ greenhouse warming on Earth and on Earth-like planets when they get in frozen state. Here we use a 3D Global Climate Model to simulate the ability of frozen planets to escape from glaciation by accumulating enough gaseous CO$_2$. We find that Earth-like planets orbiting a Sun-like star may never be able to escape from glaciation if their orbital distance is greater than $\sim$ 1.27 AU (Flux $<$ 847 W m$^{-2}$ or 62$\%$ of the Solar constant), because CO$_2$ would condense at the poles forming permanent CO$_2$ ice caps. This limits the amount of CO$_2$ in the atmosphere and thus its greenhouse effect. The amount of CO$_2$ that can be trapped in the polar caps depends on the efficiency of CO$_2$ ice to flow laterally as well as its graviational stability relative to subsurface water ice. The flow of CO$_2$ ice from poles to equator is mostly controlled by the bottom temperature, and hence by the internal heat flux. We find that a frozen Earth-like planet located at 1.30 AU of a Sun-like star could store as much as 1.5/4.5/15 bars of dry ice at the poles, for internal heat fluxes of 100/30/10 mW m$^{-2}$. But these amounts are lower limits. For planets with a significant water ice cover, we show that CO$_2$ ice deposits should be gravitationnally unstable. They get buried beneath the water ice cover in geologically short timescales of $\sim$~10$^4$ yrs, mainly controlled by the viscosity of water ice. CO$_2$ would be permanently sequestered underneath the water ice cover, in the form of CO$_2$ liquids, CO$_2$ clathrate hydrates and/or dissolved in subglacial water reservoirs (if any). This would considerably increase the amount of CO$_2$ trapped and further reduce the probability of deglaciation.

Tables

Table 1. Table 1: Physical Parameterizations used for GCM calculations.

Explored parameters	Values
CO2 partial pressures (in bar)	[0.01, 0.1, 0.4, 1.0, 1.5, 2.0, 3.0]
i) Stellar flux (in W m-2)	[1130, 1033, 949, 874, 808, 750, 697, 650]
ii) Flux compared to Earth (S${}_{\text{eff}}$)	[0.827, 0.756, 0.695, 0.640, 0.592, 0.549, 0.510, 0.476]
iii) Equivalent distance (in AU)	[1.10, 1.15, 1.20, 1.25, 1.30, 1.35, 1.40, 1.45]
Obliquity	0, 23.5∘
CO2 ice clouds	radiatively active, inactive

Table 2. Table 2: Physical Parameterizations used for the GCM calculations.

Physical parameters	Values
H2O ice Albedo	0.6
H2O ice thermal conductivity	2.5 W m${}^{-1}$K-1
H2O ice emissivity	1.0
Ground H2O ice thermal inertia	2000 J m-2 s-1/2 K-1
CO2 ice Albedo	0.6
CO2 ice thermal conductivity	0.5 W m${}^{-1}$K-1
CO2 ice emissivity	0.9
CO2 ice density	1.5
Surface Topography	Flat
Surface roughness coefficient	0.01 m

Equations

$$h_{\text{max}}=\frac{{\lambda }_{{\text{CO}}_2}(T_{\text{melt}}-T_{\text{surf}})}{F_{\text{geo}}},$$

$$\frac{1}{r}\frac{\partial}{\partial r}\Bigg{(}r~{}\frac{2A_{0}{(\rho g)}^{n}}{n+2}~{}h^{n+2}~{}{\mid\frac{\partial h}{\partial r}\mid}^{n-1}~{}\frac{\partial h}{\partial r}\Bigg{)}~{}+~{}a(r)~{}=~{}0,$$

$$h(r)=h_{\text{max}}\text{ for }r⩽R_{\text{cap,1}},$$

$$h(r)~{}=~{}h_{\text{max}}\Bigg{(}\frac{\int_{r}^{R_{\text{cap,2}}}{\Big{(}-\frac{1}{r_{2}}~{}\int_{R_{\text{cap,1}}}^{r_{2}}r_{1}~{}a(r_{1})~{}dr_{1}\Big{)}}^{\frac{1}{n}}dr_{2}}{\int_{R_{\text{cap,1}}}^{R_{\text{cap,2}}}{\Big{(}-\frac{1}{r_{2}}~{}\int_{R_{\text{cap,1}}}^{r_{2}}r_{1}~{}a(r_{1})~{}dr_{1}\Big{)}}^{\frac{1}{n}}dr_{2}}\Bigg{)}^{\frac{n}{2n+2}}\text{ for }r>R_{\text{cap,1}},$$

$${\tau }_{\text{RT}}=\frac{13\eta }{\Delta \rho gb},$$

$$U_s=\frac{2}{9}\Delta \rho g(r^2/\eta )$$

$$h_{\text{max}}=\frac{{\lambda }_{{\text{CO}}_2}(T_{\text{melt}}-T_{\text{surf}})}{F_{\text{geo}}},$$

$$T_{\text{melt}}=T_{\text{ref}}{e}^{\frac{(\frac{1}{{\rho }_{\text{liq}}}-\frac{1}{{\rho }_{\text{sol}}})}{L_{\text{fus}}}(g{\rho }_{\text{sol}}{h}_{\text{max}}+P_{{\text{CO}}_2}+P_{{\text{N}}_2}-P_{\text{ref}})},$$

$$T_{\text{surf}}=\frac{1}{\frac{1}{T_{\text{ref}}}-\frac{R}{L_{\text{sub}}{M}_{{\text{CO}}_2}}\ln (\frac{P_{{\text{CO}}_2}}{P_{\text{ref}}})},$$

$$\begin{cases}h_{\text{max}}~{}=~{}\frac{\lambda_{\text{CO}_{2}}}{F_{\text{geo}}}~{}\Big{(}\frac{1}{\frac{1}{T_{\text{ref}}}-\frac{R}{L_{\text{sub}}~{}M_{\text{CO}_{2}}}~{}\ln(\frac{P_{\text{surf}}}{P_{\text{ref}}})}\Big{)}\\ -\frac{L_{\text{fus}}}{(\frac{1}{\rho_{\text{liq}}}-\frac{1}{\rho_{\text{sol}}})~{}g\rho_{\text{CO}_{2}}}\times~{}W\Big{(}\frac{-(\frac{1}{\rho_{\text{liq}}}-\frac{1}{\rho_{\text{sol}}})~{}g\rho_{\text{CO}_{2}}~{}\lambda_{\text{CO}_{2}}~{}T_{\text{ref}}~{}e^{\frac{(\frac{1}{\rho_{\text{liq}}}-\frac{1}{\rho_{\text{sol}}})}{L_{\text{fus}}}~{}(P_{\text{surf}}-P_{\text{ref}})}}{L_{\text{fus}}~{}F_{\text{geo}}}\\ \times~{}e^{\frac{(\frac{1}{\rho_{\text{liq}}}-\frac{1}{\rho_{\text{sol}}})~{}g\rho_{\text{CO}_{2}}\frac{\lambda_{\text{CO}_{2}}}{F_{\text{geo}}}}{L_{\text{fus}}}~{}\Big{(}\frac{1}{\frac{1}{T_{\text{ref}}}-\frac{R}{L_{\text{sub}}~{}M_{\text{CO}_{2}}}~{}\ln(\frac{P_{\text{surf}}}{P_{\text{ref}}})}\Big{)}}\Big{)}\\ \end{cases}$$

$$M_{\text{tot}}=M_{\text{atm}}+M_{\text{surf}}=\frac{S{P}_{{\text{CO}}_2}}{g}+A{h}_{\text{mean}}{\rho }_{{\text{CO}}_2},$$

$${\rho }_{{\text{CO}}_2}A{h}_{\text{mean}}[P+\delta P]>{\rho }_{{\text{CO}}_2}A{h}_{\text{mean}}^{\text{max}}[P+\delta P],$$

$$M_{\text{tot}}-\frac{S[P+\delta P]}{g}>{\rho }_{{\text{CO}}_2}A{h}_{\text{mean}}^{\text{max}}[P+\delta P],$$

$$M_{\text{tot}}-\frac{SP}{g}-\frac{S\delta P}{g}>{\rho }_{{\text{CO}}_2}A{h}_{\text{mean}}^{\text{max}}[P]+\delta P{\rho }_{{\text{CO}}_2}\frac{d(A{h}_{\text{mean}}^{\text{max}})}{d{P}_{{\text{CO}}_2}}[P].$$

$$\Big{(}M_{\text{tot}}-\frac{S~{}P}{g}-\rho_{\text{CO}_{2}}~{}Ah_{\text{mean}}[P]\Big{)}-\frac{S~{}\delta P}{g}>\delta P~{}\rho_{\text{CO}_{2}}~{}\frac{d(Ah_{\text{mean}}^{{\text{max}}})}{dP_{\text{CO}_{2}}}[P].$$

$$A\frac{d{h}_{\text{mean}}^{\text{max}}}{d{P}_{{\text{CO}}_2}}[P]+h_{\text{mean}}^{\text{max}}\frac{dA}{d{P}_{{\text{CO}}_2}}[P]<-\frac{S}{{\rho }_{{\text{CO}}_2}g}.$$

$$\frac{d{h}_{\text{mean}}^{\text{max}}}{d{P}_{{\text{CO}}_2}}<-\frac{S}{gA{\rho }_{{\text{CO}}_2}}=\frac{-10}{9.81\times 1.5\times 10^3}=-7\times 10^{-4}\text{m}{\text{Pa}}^{-1}.$$

Full text

1 Introduction

There are geological evidences that Earth had to face multiple episodes of global or quasi-global glaciation, at the end of the Archean 2.45-2.22 Gya and during the Neoproterozoic era, 710 and 650 Mya (Evans et al.,, 1997; Hoffman et al.,, 1998; Hoffman and Schrag,, 2002).

It is widely believed that the carbonate$-$silicate cycle (Walker et al.,, 1981; Kasting et al.,, 1993; Kump et al.,, 2000) was the main agent to trigger deglaciations on ancient Earth. In particular, on a completely frozen planet, also called ’Hard Snowball’, the weathering of CO2 is stopped ; continuous volcanic outgassing builds up atmospheric CO2 that warm up the climate until liquid water is produced in the equatorial regions.

By extension, this mechanism may be crucial to stabilize the climate of Earth-like exoplanets. It is even central for the definition of the classical Habitable Zone (Kasting et al.,, 1993; Forget and Pierrehumbert,, 1997; Kopparapu et al.,, 2013), which assumes that planets can build up CO2 atmospheres (as massive as wanted) suitable for the stability of surface liquid water.

Planets near the outer edge of the Habitable Zone may suffer, as Earth did, episodes of complete glaciation. Recent work by Menou, (2015) has even shown that, in the case of planets lacking land vascular plants, temperate climates may not be stable, because at high CO2 partial pressure, the increased weathering rate does not allow temperate stable solutions. This effect, which is enhanced for planets weakly irradiated by their star, may systematically drive Earth-like exoplanets toward episodes of glaciation.

Can Earth-like planets in the Habitable Zone of their star always escape from episodes of glaciation ? In this work, we use 3D Global Climate Model (GCM) simulations of snowball planets to study the ability of the increased CO2 greenhouse effect resulting from volcanic outgassing to drive them out of glaciation. As the CO2 outgassed by volcanoes accumulates in the atmosphere, the temperature of condensation of CO2 can exceed the surface temperature of the poles. This leads to the trapping of extra outgassed CO2, forming permanent CO2 polar ice caps, and thus seriously limits the efficiency of the carbonate$-$silicate cycle. This possibility has already been suggested in Pierrehumbert, (2005) and Pierrehumbert et al., (2011), but has never been explored yet. We propose in this paper a detailed study of this scenario.

2 Method

We use in this paper the 3-Dimensions LMD Generic Global Climate Model to study the deglaciation of Earth-like planets orbiting in a circular orbit around a Sun-like star (Sun spectrum and luminosity) in response to increasing amounts of atmospheric CO2 (from 0.01 to 3 bars), and for orbital distances ranging from 1.10 to 1.45 Astronomical Units (AU). Detailed information on the model can be found in Appendix A.

Our simulations were designed to represent completely frozen Earth-like planet characteristics. These include the radius (6370 km), the gravity field (9.81 m s*-2*), the obliquity (23.5∘) and the rotation speed (7.28$\times$10*-5* rad s*-1*). The roles of obliquity, planetary mass and rotation rate on the ability of planets to escape glaciation episodes are discussed in sections 3.2 and 3.4. Eventually, most of the simulations were performed for a uniformly flat topography. The effect of topography is discussed in section 3.4.5.

All the simulations performed in this study were forced initially in a cold and dry state, assuming:

A uniform and complete ice cover. 2. 2.

Atmospheric and surface temperatures arbitrarily fixed to 230 Kelvins everywhere. 3. 3.

No water vapor, no clouds.

Depending on CO2 partial pressure, obliquity, or parameterization of clouds, the simulations evolve in different steady state climate regimes that are discussed in the next section.

3 Glaciation escape limited by CO2 atmospheric collapse

3.1 Reference case

We performed 40 simulations of Earth-like planets (as described in section 2 – with a 23.5*∘* obliquity and radiatively inactive CO2 ice clouds) starting from an initially cold (T = 230 K everywhere) and global glaciation state for multiple irradiation fluxes and CO2 partial pressures (see table 1). We fix the N2 partial pressure to be constant and equal to 1 bar, in order to be consistent with the case of Earth. All simulations were run long enough to reach equilibrium111up to 30 Earth years for the thickest atmospheric configurations..

We find that, depending on the CO2 partial pressure and the stellar flux, these simulations can evolve in three different climate regimes:

The greenhouse effect of CO2222which is enhanced by the pressure broadening of N2. is sufficient to raise the surface temperatures in equatorial regions above the melting temperature of water ice333273 K, here.. In this case, the positive feedback due to the decrease of the surface albedo (from 0.6 to 0.07) and the greenhouse effect of water vapor drive the planet out of glaciation. 2. 2.

The greenhouse effect of CO2 is too weak to trigger a deglaciation. The planet stays in a snowball state. 3. 3.

The greenhouse effect of CO2 is too weak to raise the surface temperatures of the poles above the condensation temperature of CO2. In this case, atmospheric CO2 collapses and the planet is locked in a global glaciated state, with two permanent CO2 ice polar caps.

Figure 1 shows that, for a given solar flux (808 W m*-2* here), the choice of the initial CO2 partial pressure can either drive the planet out of glaciation (pCO2 = 3 bars) or cause a permanent collapse of the atmosphere at the poles (pCO2 = 1 bar).

These results are summarized in Figure 2b, for multiple stellar fluxes and CO2 partial pressures. Each region of the diagram denotes a steady state climate regime reached by the planet:

In red, the planet is partially or totally deglaciated. 2. 2.

In white, the planet is entirely frozen. 3. 3.

In blue, the planet is entirely frozen and gaseous CO2 has permanently collapsed at the poles. This situation actually occurs when winter CO2 frost formation rate exceeds seasonal summer sublimation.

We now put these results in context of an active carbonate-silicate cycle. A planet that enters in a global glaciation state must have initially a low CO2 atmospheric content that places it in the lower part of the diagram (in Figure 2b). From this point, CO2 is outgassed by volcanoes and accumulates in the atmosphere. As the CO2 partial pressure increases, the planet moves up in the diagram until:

the planet reaches the red zone first. The planet is able to escape glaciation. 2. 2.

the planet reaches the blue zone first. All the extra CO2 possibly outgassed by volcanoes condenses at the poles and the planet is locked in a perpetual glaciation state.

For an (Earth-like) 23.5*∘* obliquity, we find this limit to occur for planets located at more than 1.27 AU from a Sun-like star (Flux $<$ 847 W m*-2*, S${}_{\text{eff}}$444The effective flux S${}_{\text{eff}}$ is defined as the ratio between the incoming stellar flux on the planet and that on the Earth. $<$ 0.62), with an optimal CO2 partial pressure around 1 bar (see Figure 2b).

However, if a planet starts with a large enough CO2 atmospheric content so that it lies above the blue zone, then the CO2 should resist atmospheric collapse and the deglaciation becomes possible.

3.2 Null obliquity case

For a 0*∘* obliquity planet, because insolation in the polar regions is lowered, the equator-to-pole temperature contrast is amplified. Consequently, we find that permanent condensation of CO2 at the poles happens for much lower equivalent distances (d $>$ 1.13 AU, Flux $<$ 1070 W m*-2*, S${}_{\text{eff}}$ $<$ 0.78 ; see Figure 2a).

3.3 Effect of CO2 ice clouds

During polar nights, the temperatures at the poles are very low. CO2 condenses at the surface but also in the atmosphere, forming CO2 ice clouds that may have a significant warming effect (Forget and Pierrehumbert,, 1997).

We performed multiple simulations for planets at 23.5*∘* obliquity including radiatively active CO2 ice clouds. For this, we used a constant cloud condensation nuclei [CCN] of 105 kg*-1* designed to maximize the greenhouse effect of CO2 ice clouds (Forget et al.,, 2013). This idealized configuration corresponds to an endmember description of the radiative warming effect of CO2 ice clouds. In particular, Kitzmann, (2016) has shown recently, using a more refined radiative transfer scheme, that a configuration like the one used in this section probably overestimates the warming effect of CO2 ice clouds.

Our simulations (example in Figure 3) show that the CO2 ice clouds distribution evolves significantly with seasons. In winter, CO2 ice clouds form a quasi-complete polar cover and have a powerful greenhouse effect that limits the condensation of CO2 at the poles. In summer, CO2 clouds dissipate as insolation increases. As a result, CO2 clouds have a very limited impact on the albedo and therefore do not contribute much to prevent polar CO2 ice summer sublimation.

For these two reasons, CO2 ice clouds seriously limit the collapse of CO2 at the poles. Figure 2c shows that, with radiatively active CO2 clouds, CO2 condensation would occur for equivalent distance $>$ 1.40 AU (Flux $<$ 697 W m*-2*, S${}_{\text{eff}}$ $<$ 0.51).

3.4 Sensitivity study

3.4.1 The albedo of ices

The choice of water ice albedo can severely affect our results. For albedo higher than 0.67, we find that CO2 condensation can occur on Earth-like planets (23.5∘ obliquity, radiatively inactive CO2 ice clouds) for equivalent distance as low as 1.15 AU (Flux $\sim$ 1033 W m*-2*, S${}_{\text{eff}}$ $\sim$ 0.76), corresponding to the reduced luminosity of the Sun 3.8 Gya calculated from standard solar evolution models (Gough,, 1981).

We note that water ice albedo is considerably reduced around cool stars (Joshi and Haberle,, 2012), making the scenario of CO2 polar condensation less efficient. However, planets orbiting in the Habitable Zone of M-dwarfs are also subject to tidal locking. The temperature on the nightside of a synchronous planet can be extremely low, favoring the CO2 condensation. This possibility is explored and discussed in details in Turbet et al., 2017b .

We also explored the effect of CO2 ice albedo that could potentially be very high (see section A.3, Kieffer et al., 2000) and we found it to be much less important.

Eventually, we acknowledge that a real Snowball might have regions of open continent not covered by water ice, which would lower albedo. These could act as sources of dust, and volcanic aerosols could also darken the surface (Abbot and Pierrehumbert,, 2010).

3.4.2 The rotation rate

The rotation rate of Earth has evolved in time from 1.2$\times$10*-4* (4 Gya) to 7.3$\times$10*-5* rad s*-1* (now) due to the effect of the Moon’s tidal friction (Walker and Zahnle,, 1986; Zahnle and Walker,, 1987). More generally, Earth-like planets could harbor a wide range of rotation rates.

In this study, we find that rotation rate is a very important parameter for CO2 polar condensation. A high rotation rate reduces the latitudinal transport and thus increases the equator-to-pole surface temperature gradient (Kaspi and Showman,, 2015), favoring CO2 condensation at the poles. We find that, for a rotation rate $\Omega$ = 4 $\Omega_{\text{Earth}}$ (23.5∘ obliquity, radiatively inactive CO2 ice clouds), CO2 can collapse as early as 1.15 AU (Flux $\sim$ 1033 W m*-2*, S${}_{\text{eff}}$ $\sim$ 0.76).

To go further, we investigated the case of archean Earth, 3.8 Gya (S${}_{\text{eff}}$ $\sim$ 0.76), through simulations of a completely frozen planet, assuming a 23.5*∘* obliquity, a rotation period of 14 hours and a water ice albedo set to 0.6. We find that CO2 can condense at the poles seasonally, but never permanently.

3.4.3 The planetary mass

Massive planets also have a large radius that may be responsible for the weakening of the poleward heat transport, favoring the collapse of CO2. For example, assuming an internal composition similar to the Earth, a 5 M${}_{\text{Earth}}$ planet would have a radius around 1.5 R${}_{\text{Earth}}$ or 104 km (Seager et al.,, 2007). We performed GCM simulations for a 5 M${}_{\text{Earth}}$ planet, for a CO2 partial pressure of 1 bar and a N2 partial pressure of 1 bar. In these conditions, we find that the collapse of CO2 can occur for equivalent distances as low as 1.18 AU (Flux $<$ 981 W m*-2*, S${}_{\text{eff}}$ $<$ 0.72).

However, because accretion of volatiles should increase with planetary mass, the atmospheric mass should also increase with it (see Kopparapu et al., (2014) section 2 for more details). Therefore, we expect a massive planet to have a greater content of background gas (e.g. N2). For a 5 M${}_{\text{Earth}}$ planet, assuming that P${}_{\text{surface}}$ scales as R${}_{\text{planet}}$2.4. (equation 3 in Kopparapu et al., (2014)), we roughly get a partial pressure of N2 equal to 3.3 bars 555We remind the reader that volatile delivery to a planet is stochastic in nature and may be a weak function of planetary mass. By way of comparison, Venus (at $\sim$ 0.8 M${}_{\text{Earth}}$) has $\sim$ 3 bars of N2 in its atmosphere.. Using GCM simulations, we find that in this more realistic case, CO2 collapses for equivalent distances greater than 1.30 AU (Flux $<$ 808 W m*-2*, S${}_{\text{eff}}$ $<$ 0.59).

3.4.4 The atmospheric composition: the role of N2 partial pressure

The result of the previous section shows that N2 partial pressure can have a significant impact on CO2 collapse. As pN2 increases, the atmosphere thickens and (1) the poleward heat transport increases and (2) the greenhouse effect of CO2 increases because of pressure broadening of CO2 by N2. Both effects tend to prevent CO2 condensation.

Quantitatively, our simulations show that, for an Earth-like planet with a CO2 partial pressure of 1 bar (with 23.5∘ obliquity, radiatively inactive CO2 ice clouds), doubling pN2 from 1 to 2 bars leads to the condensation of CO2 for equivalent orbital distances greater or equal to 1.35 AU (+0.08 AU difference).

3.4.5 The topography

To investigate the possible effects of topography on CO2 polar condensation, we ran several simulations in which we emplaced a mountain range similar in size to the Himalaya (5000 km altitude plateau, 106 km2) at several latitudes (0∘N, 30∘N, 45∘N, 60∘N, 90∘N). We find that CO2 permanent condensation can occur on top of the mountain range for distances lower than 1.27 AU, only for latitudes greater than 45∘N. For the 60∘N latitude case, CO2 condensation starts as low as 1.20 AU (Flux $\sim$ 949 W m*-2*, S${}_{\text{eff}}$ $\sim$ 0.69).

Nonetheless, the total amount of CO2 that could condense and be trapped above a mountain range is rather low, because it is limited by the area of the mountain range.

4 How much CO2 ice can be trapped ?

In this section, we investigate two processes that should control the amount of CO2 that can be trapped on the surface or subsurface: 1) flows of polar CO2 ice to lower latitudes and sublimation; and 2) burial of CO2 ice beneath the water ice cover due to higher density.

4.1 Maximum size of CO2 ice glaciers

At first sight, the main limit of the trapping of CO2 as ice instead of greenhouse gas is the size of the solid CO2 polar reservoirs. When CO2 is outgassed by volcanoes, the atmospheric pressure stays constant but the size of the CO2 polar caps grows. At some point, the CO2 ice caps form glaciers that can flow efficiently toward equatorial regions, CO2 ice being much less viscous than water ice. In the process, CO2 ice can be sublimated and reinjected in the atmosphere.

In this section, we give estimates of the maximum amount of CO2 that can be trapped in steady state CO2 ice polar caps, for a planet with a 23.5*∘* obliquity, radiatively inactive CO2 ice clouds, located at a distance of 1.30 AU (F = 808 W m*-2*) and for both CO2 and N2 partial pressures set to 1 bar.

4.1.1 CO2 ice caps radial extent

First, we get from our GCM simulations the radial extent of the two permanent CO2 ice polar caps, defined for a positive annual mean CO2 condensation rate. We note the corresponding radius of the ice cap $R_{\text{cap,}1}$. Then, in GCM simulations, we artificially extend (in the radial direction) the CO2 ice caps to take into account the glacier flow, until the net globally averaged annual condensation$/$sublimation rate vanishes. In practice, we used the following algorithm:

We run the simulation until globally averaged annual mean rate of condensation/sublimation is constant. 2. 2.

If this rate is (roughly) null, CO2 ice caps are in a dynamical equilibrium. We stop here. 3. 3.

If this rate is positive, we artificially increase the size of the CO2 ice caps of one GCM grid (in the latitudinal direction, and for each longitude) by emplacing arbitrarily a large enough amount of CO2 ice. We go back to step 1.

This method gives us a good estimate of the steady state maximal extent of the CO2 ice polar caps. We note the corresponding radius of the ice cap $R_{\text{cap,}2}$.

With and without glacier flow, we obtain polar cap radii $R_{\text{cap,}1}$ = 2070 km and $R_{\text{cap,}2}$ = 2800 km. This corresponds to a latitudinal extent down to 71.5*∘N and 65∘*N 666By symmetry, results are identical in the southern hemisphere., respectively. Figure 4 shows zonal condensation/sublimation annual mean rates along the North Pole CO2 glacier.

In total, the two CO2 ice caps have an area of 27 millions of km2 (respectively 49 millions of km2 if considering the glacier flow). This corresponds to 5$\%$ (respectively 10$\%$) of the total area of the (Earth-like) planet.

4.1.2 CO2 ice cap inner region: a thickness limited by basal melting

In the region within the CO2 ice caps with positive annual mean condensation rates, the CO2 ice maximum thickness should be mainly limited by basal melting induced by the geothermal heat flux, noted $F_{\text{geo}}$.

Assuming that the temperature inside the CO2 glacier rises linearly with depth, with a lapse rate fixed by the geothermal heat flux (conductive regime), the maximum CO2 ice thickness $h_{\text{max}}$ is given by:

[TABLE]

where $\lambda_{\text{CO}_{2}}$ is the thermal conductivity of CO2 ice, $T_{\text{surf}}$ is the temperature at the top of the glacier and $T_{\text{melt}}$ is the melting temperature of CO2 ice at the base of the glacier. As the latter is a function of the pressure below the CO2 glacier, it implicitly depends on $h_{\text{max}}$ so that the above equation must be solved numerically, once the local surface temperature is known (see Appendix B for details on the calculations). The resulting maximum CO2 ice cap thickness $h_{\text{max}}$ is plotted on Figure 5 as a function of the internal heat flux and the CO2 partial pressure. For geothermal heat fluxes of 100/30/10 mW m*-2* (red stars in Figure 5), basal melting condition gives CO2 ice maximum thicknesses of 120/360/1200 m.

The exact latitude at which the glacier maximum thickness is determined by either basal melting condition or by glacier flow is difficult to determine. Moreover, in the region of transition between these two regimes, a basal liquid CO2 flow could persist777Such basal flow could carry heat away from basal melting regions and thus also reduce the heat flux conducted through the CO2 ice caps (Cuffey and Paterson,, 2010; Menou,, 2013). and therefore speed up the glacier flow.

For simplicity, we fixed the thickness of the steady-state CO2 glacier in the regions with positive annual mean condensation rate (regions of latitude $>$ 71.5*∘*N ; see Figure 4) at the constant maximum thickness h${}_{\text{max}}$ derived from the basal melting condition. This maximum thickness serves as a boundary condition for the calculation of the glacier flow in regions with positive annual mean sublimation rate (see next section).

4.1.3 CO2 ice cap outer region: a thickness limited by glacier flow and sublimation

To get estimates of CO2 glacier profiles in the regions where sublimation dominates over condensation (latitude $<$ 71.5*∘*N, see Figure 4), we use a simple isothermal glacier model 888This is a simple approach that does not physically describe how the accumulating CO2 ice spreads to feed the ablation zone. We acknowledge that a more physically complete model should be developed in the future. (assuming no slip at the base of the glacier), following a similar approach than Menou, (2013). The steady-state equation satisfied by a thin isothermal CO2 glacier with a flat base and a no slip boundary condition at its bottom is (Fowler,, 2011; Paterson,, 2010)999This is the classic solution for a glacier shape.:

[TABLE]

with $h(r)$ the thickness of the CO2 ice cap, assumed (for simplicity) radially symmetric, A0 the flow rate constant (in Pa s*-n*), n the power-law creep exponent, g the surface gravity, $\rho$ the density of CO2 ice and a(r) the annual mean condensation rate assumed (also for simplicity) to be a function of r only, and derived from our GCM simulations (see Figure 4).

First, assuming $\frac{\partial h}{\partial r}~{}=~{}0$ for $r~{}=~{}R_{\text{cap,1}}$, we integrate equation 2 once. Then, assuming $h~{}>~{}0$, $\frac{\partial h}{\partial r}~{}<~{}0$ and $a(r)<0$, we separate the variables $h$ and $r$, and we integrate a second time, assuming $h$($R_{\text{cap,2}}$) = [math] at the edge of the glacier. Eventually, assuming that $h$($R_{\text{cap,1}}$) = $h_{\text{max}}$, with $h_{\text{max}}$ the maximum thickness of the glacier calculated using the basal melting condition, we can constrain the flow rate A0 and get the following solutions for the thickness $h$ of the CO2 ice glacier:

[TABLE]

[TABLE]

Figure 6 presents 9 possible steady-state CO2 ice glacier profiles, for 3 different geothermal fluxes and 3 different flow laws (n = 1, 4.5 and 7) chosen consistent with previous experimental studies (Durham et al.,, 1999) and works about the martian CO2 southern polar cap (Nye et al.,, 2000).

4.1.4 Calculation of the CO2 maximum total reservoir

After integration of the steady state CO2 ice radial profiles, we find that the maximum amount of CO2 ice that can be stored in the two polar caps is approximately equal to 1.5/4.5/15 bars101010These quantities are not affected much by the choice of flow law (see Figure 6). for geothermal heat fluxes of 100/30/10 mW m*-2*, respectively. On Earth, the current geothermal heat flux is $\sim$ 80 mW m*-2* but can varies a lot, typically from 20 mW m*-2* up to 400 mW m*-2*, depending on the regions (Davies,, 2013). Approximately 60$\%$ (respectively 40$\%$) of CO2 ice is trapped in region of $r<R_{\text{cap,1}}$ (respectively $r>R_{\text{cap,1}}$), meaning that the accumulation zone and the ablation zone contribute significantly and equally to the total size of the CO2 ice reservoir.

It is now possible to calculate, for a given rate of CO2 volcanic outgassing, the time required for a hard snowball planet to fill the CO2 ice polar reservoirs and thus potentially escape glaciation. Assuming a present-day Earth CO2 volcanic outgassing rate of $\sim$ 7 teramoles/year or 60 bars/Gy (Brantley and Koepenick,, 1995; Jarrard,, 2003), we get durations of glaciation periods of $\sim$ 42/93/270 My for internal heat fluxes of 100/30/10 mW m*-2*. In fact, because CO2 volcanic outgassing rate is dependant on the internal heat flux, the duration of the glaciation phases could be drastically increased for planets with low geothermal heat fluxes (typically lower than 80 mW m*-2* as on present-day Earth). This suggests that planets near the outer edge of the habitable zone may need to be quite volcanically active in order to remain habitable.

When (and if) CO2 ice caps become full, the system suddenly reaches a runaway phase, ending up with the complete vaporization of CO2 polar ice caps and thus a hot atmosphere. This may have implications in the calculation and description of ’limit-cycles’ (Menou,, 2015; Haqq-Misra et al.,, 2016; Abbot,, 2016) and more generally on the habitability of distant Earth-like planets. Detailed calculations of the stability of CO2 ice caps when full can be found in Appendix C.

It is also conceivable that catastrophic events such as meteoritic impacts could trigger deglaciation.

4.2 Gravitational stability and CO2 sequestration

For planets covered by several hundred meters of water ice111111We remind here that the Earth Global Equivalent Layer of surface water is $\sim$ 3km., owing to the high density of CO2 ice ($\sim$1500 kg m*-3*) or liquids ($\sim$1170 kg m*-3*) compared to water ice ($\sim$930 kg m*-3* at 200-220 K), the CO2 deposits should rapidly become gravitationally unstable. This would result in a burial of the CO2 deposits at the base of the H2O ice layer. Such a density contrast should initiate Rayleigh-Taylor instabilities at a timescale that can be estimated at first order assuming two isoviscous layers using the following relationship (Turcotte and Schubert,, 2001):

[TABLE]

with $\eta$ the viscosity of the more viscous layer, $\Delta\rho$ the density contrast between the two layers and $b$ the characteristic size of the domain. Depending whether CO2 is liquid or solid at the interface with the H2O layer, the density contrast would range between 240 and 570 kg m*-3*. For a CO2 ice deposit thickness of about 100 m, this density contrast results in a stress at the interface between the two layers of the order of 1 MPa. For such a stress level and temperature ranging between 200 and 220 K, the viscosity of the CO2 ice layer is estimated between 1012-1013 Pa s, based on available experimental data (Durham et al.,, 1999). At the same temperature and stress conditions, water ice has a viscosity $\leq$ 1016 Pa s, for grain size lower than 1 mm, based on experimental data (Durham et al.,, 2001; Durham and Stern,, 2001; Goldsby and Kohlstedt,, 2001). Therefore, this would be the layer controlling the Rayleigh-Taylor timescale. Assuming that a thickness of 100 m of CO2 deposit is representative of the characteristics size of the domain, the R-T timescale is 7$\times$103 yrs $\times$ ($\eta$/1016 Pa s), which is geologically short. In the case of basal melting of the CO2 layer, we should keep in mind that the generated CO2 liquids will be unstable relative to the underneath water ice layer, but also relative to the overlying CO2 ice layer as it has a density intermediate between those of CO2 ice and water ice. When CO2 melting occurs at the base of the CO2 layers, a competition is expected between downward drainage through the water ice layer (similar to crevasse hydrofracturing in Earth’s glaciers; Krawczynski et al., 2009) and liquid injection in the overlying CO2 ice layer and refreezing (similar to dike propagation). Assessing in details such complex processes are out of scope of the present paper. Nevertheless, based on these considerations, we can safely argue that the accumulation of a large volume of liquid CO2 at the base of the CO2 ice layer is unlikely.

Once gravitationally destabilized, the CO2 ice deposit would sink at the base of the water ice shell at a rate that is determined mostly by the viscosity of water ice and the size of the CO2 ice diapir. The time required for a CO2 ice diapir to cross the water ice layer can be estimated using the Stokes velocity, the terminal velocity of a sphere falling through a constant viscosity medium (Ziethe and Spohn,, 2007):

[TABLE]

For a diapir radius $r$ of 100 m (comparable to thickness of the CO2 deposit) and a viscosity of water ice of 1015-1016 Pa s, this leads to a velocity of 0.04-0.4 m$/$yr. As temperature increases as a function of depth, the viscosity of water ice is expected to decrease with depth, resulting in an acceleration of the diapir fall. A 100-m diapir of CO2 ice would thus not need more than $\sim$ 10${}^{4}~{}$yrs to reach the bottom of a 1-km thick water ice layer.

Even if the CO2-H2O interface is beneath the melting temperature of CO2 ice, melting may occur during the diapir fall, which will reduce the density contrast and hence the fall velocity. However, as mentioned above, melting may promote fracturing of the water ice medium and rapid drainage of the CO2 liquids due to its high density relative to the surrounding water ice. The rate of melting should also depend on the efficiency of heat transfer between the diapir and the surrounding ice. As the thermal diffusive timescale ($\sim$r${{}^{2}/\kappa}$) for a 100-m-size diapir of CO2 is of the order of 103 yrs and is thus $\leq$ to the expected diapir fall timescale, melting during the descent would be efficient. Detailed modelling would be required though to determine the exact volume of generated melt during descent.

CO2 ice should completely melt and equilibrate thermally with the surrounding H2O media once stabilized at the bottom of the water ice shell. The temperature and pressure conditions at the bottom of the water ice layer depend on its thickness and on the geothermal flow. For geothermal heat flux between 50 and 100 mW -2, typical of the Earth, the melting of water ice would be reached for depth ranging between 1.5 and 3 km, corresponding to pressure of about 14 and 28 MPa. This pressure range is well above the saturation vapor pressure (3.5 MPa at 273 K; Lide, 2004), so that CO2 is highly stable in the liquid phase. Destabilizing the CO2 liquids into gas phase would require the 273 K-isotherm at a depth of only 375 m, corresponding to abnormally high geothermal heat flux of 400 mW -2, 5 times larger than the average heat flux on Earth. Even if the density of liquid CO2 decreases with increasing temperature as it equilibrates with the surrounding water ice media, it remains always denser than water ice (Span and Wagner,, 1996), and therefore should always accumulate at the bottom of the ice shell. At T = 273 K and pressure between 14 and 28 MPa, CO2 liquids have a density very close to that of liquid water (994 and 1048 kg m*-3*, respectively, using the equation of state of Span and Wagner, (1996)), so that CO2 coexists with H2O at the ice-water interface.

Pressure and temperature conditions expected at the bottom of the ice layer are in the stability field of CO2 clathrate hydrate (Sloan,, 1998; Longhi,, 2005), therefore CO2 should rapidly interact with H2O molecules to form clathrate hydrate. Clathrate hydrates are non-stoichiometric compounds consisting of hydrogen-bonded H2O molecules forming cage-like structures in which guest gas molecules, such as CO2, can be trapped (Sloan,, 1998). Once formed, these clathrates are very stable and can be dissociated only if the temperature is raised about 5-10 K above the melting point of water ice. The storage of CO2 in the form of clathrate should be particularly efficient in the case of warm bottom conditions as CO2 liquids and liquid water coexist. In the cold bottom condition with no water ice melting, clathration is still possible but would be less efficient due to the physical separation between the two components. As CO2 clathrate hydrate has a density of about 1150 kg m*-3* (assuming full cage occupancy, Sloan, (1998)), they would rapidly sink at the bottom of the water liquid layer, ensuring an almost complete clathration of CO2. Part of the CO2 should dissolve in the liquid water during the clathrate sinking. The relative proportion of CO2 trapped in the form of clathrate hydrate or dissolved in the water layer would depend on the volume of CO2 that is buried at the base of the ice shell and on the volume (thickness) of the water layer.

In summary, as long as the water ice shell exceeds a few hundred meters, CO2 should remain sequestered either in the form of CO2 liquids, in the form of CO2 clathrate hydrate or dissolved in the subglacial water layer. Release of gaseous CO2 may occur only due to an increase of surface insulation and increase of geothermal flux resulting in a significant thinning and breaking-up of the water ice shell. The total amount of CO2 that can be stored in the H2O layer (by any of these three processes) depends on the total abundance of H2O of the planet as well as CO2/H2O ratio. The CO2 inventory on Earth, mostly in the form of carbonate rocks on the continents, is of the order of 102 bars (Walker,, 1985). If the total CO2 inventory exceeds the total amount of CO2 that can be stored in the water layer, then the planet should be able at some point to escape from global glaciation, extra CO2 being likely returned to the atmosphere through cryovolcanism.

Evaluating the maximum amount of CO2 that can be trapped underneath the water ice cover require however a detailed description of the H2O layer structure as well as thermodynamic models predicting the partitionning of CO2 between the different phases.

5 Conclusions

We highlight in this paper a new scenario that would prevent distant Earth-like planets orbiting Sun-like stars from escaping episodes of complete glaciation. When a terrestrial planet reaches a Snowball state, CO2 weathering ceases and CO2 can accumulate in the atmosphere because of volcanic outgassing. As CO2 builds up in the atmosphere, the temperature of condensation of CO2 can exceed the surface temperature of the poles, which leads to the trapping of all extra CO2 possibly outgassed.

Using LMD Global Climate Model simulations designed for Earth-like fully glaciated planets, we show that this mechanism can work from orbital distances as low as 1.27 AU (Flux $\sim$ 847 W m*-2*, S${}_{\text{eff}}$ $\sim$ 0.62). By comparison, the most recent estimate of the outer edge of the Habitable Zone (e.g. the maximum greenhouse limit) in similar conditions is 1.67 AU (Flux $\sim$ 490 W m*-2*, S${}_{\text{eff}}$ $\sim$ 0.36) (Kopparapu et al.,, 2013).

This limit can occur at even lower distances for planets with 1) a low obliquity, 2) a low N2 partial pressure, 3) a high rotation rate and 4) a high value of the water ice albedo (we chose in this work the conservative value of 0.6).

Conversely, this limit can occur at significantly higher distances when taking into account the radiative effect of CO2 ice clouds. In winter, CO2 ice clouds that form at the cold poles scatter back thermal infrared radiation of the surface. In summer, clouds are dissipated and therefore have almost no impact on the bond albedo. For these two reasons, CO2 ice clouds can have a powerful warming effect at the poles, limiting CO2 collapse at orbital distances greater than 1.40 AU (Flux $<$ 697 W m*-2*, S${}_{\text{eff}}$ $<$ 0.51).

For each possible configuration, the amount of CO2 that can be trapped in polar CO2 ice caps depends on the efficiency of CO2 ice to flow laterally as well as its graviational stability relative to subsurface water ice. The flow of CO2 ice from polar to low latitudes regions is mostly controlled by its basal temperature (due to both the conductivity and melting temperature of CO2 ice being low), and hence by the internal heat flux of the planet. We find for example that a frozen Earth-like planet located at 1.30 AU of a Sun-like star could store as much as 1.5/4.5/15 bars of CO2 ice at the poles, for internal heat fluxes of 100/30/10 mW m*-2*.

But these amounts are in fact lower limits. CO2 ice being denser than water ice ($\sim$ 1.7 $\times$ the volumetric mass of water ice), we find that CO2 ice deposits should be gravitationally unstable and get buried beneath the water ice cover in a geologically short timescale of $\sim$ 104 yrs, mainly controlled by the viscosity of water ice. If water ice cover exceeds about 300 meters (e.g. 10$\%$ of the Earth hydrosphere), then CO2 should be permanently sequestered underneath the water ice cover, in the form of CO2 liquids, CO2 clathrate hydrates and/or dissolved in subglacial water reservoirs (if any). This would considerably increase the amount of CO2 trapped and further reduce the probability of deglaciation.

6 acknowledgements

M.T. thanks Emmanuel Marcq for fruitful discussions related to this work during the legendary LMD PLANETO meetings. We acknowledge the feedbacks from Jim Kasting, William B. McKinnon, Dorian Abbot and an anonymous referee, which improved both the style and content of the manuscript. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No. 679030/WHIPLASH).

Appendix A The LMD Generic Global Climate Model

A.1 Model core

This model originally derives from the LMDz 3-dimensions Earth Global Climate Model (Hourdin et al.,, 2006), which solves the primitive equations of geophysical fluid dynamics using a finite difference dynamical core on an Arakawa C grid.

The same model has been used to study many different planetary atmospheres of low-irradiated planets, including Archean Earth (Charnay et al.,, 2013), past climates of Mars (Forget et al.,, 2013; Wordsworth et al.,, 2013; Turbet et al., 2017a, ), or exoplanets (Wordsworth et al.,, 2011; Turbet et al.,, 2016; Turbet et al., 2017b, ).

The simulations presented in this paper were all performed at a horizonal resolution of 96 $\times$ 96 (e.g. 3.75*∘* $\times$ 1.875*∘*) in longitude $\times$ latitude. In the vertical direction, the model is composed of 15 distinct atmospheric layers that were designed using hybrid $\sigma$ coordinates (where $\sigma$ is the ratio between pressure and surface pressure).

To account for thermal conduction in the icy ground, we used a 18-layers thermal diffusion soil model. Mid-layers depths range from d0 $\sim$ 0.1 mm to d17 $\sim$ 18 m, following the power law dn = d0 $\times$ 2n with $n$ being the corresponding soil level, chosen to take into account both diurnal and seasonal thermal waves. We assumed thermal inertia of the ground ice to be constant and equal to 2000 J m*-2* s*-1/2* K*-1*.

We considered (in GCM simulations) the internal heat flux and/or the thermal heat flux F${}_{\text{ground}}$ conducted from an hypothetical underlying ocean to be zero121212This is a reliable assumption for thick enough (typically $>$ 100 m) water ice covers expected on hard Snowball Earth-like planets..

A.2 Radiative Transfer

The GCM includes a generalized radiative transfer for a variable gaseous atmospheric composition made of a mix of CO2, N2 and H2O (HITRAN 2012 database (Rothman et al.,, 2013)) using the correlated-k method (Fu and Liou,, 1992; Eymet et al.,, 2016)) suited for fast calculation. For this, we decomposed atmospheric temperatures, pressures, and water vapor mixing ratio into the following respective 7 x 8 x 8 grid, similar to Charnay et al., (2013): Temperatures = $\{$100,150, .. ,350,400$\}$ K ; Pressures = $\{$10*-6*,10*-5*, .. ,1,10$\}$ bars ; H2O Mixing Ratio = $\{$10*-7*,10*-6*, .. ,10*-2*,10${}^{-1}\}$  mol of H2O $/$ mol of air (H2O+CO2 here).

CO2 collision-induced absorptions (CIA) and dimer absorption (Wordsworth et al.,, 2010) ) were included in our calculations as in Forget et al., (2013) and Wordsworth et al., (2013), as well as N2-N2 collision-induced absorption (Richard et al.,, 2012) and its role on the pressure broadening. We also added H2O self and foreign continuums calculated with the CKD model (Clough et al.,, 1989), with H2O lines truncated at 25 cm*-1*.

For the computation, we used 32 spectral bands in the thermal infrared and 36 in the visible domain. 16 non-regularly spaced grid points were used for the g-space integration, where g is the cumulated distribution function of the absorption data for each band. We used a two-stream scheme (Toon et al.,, 1989) to take into account radiative effects of aerosols (CO2 ice and H2O clouds) and Rayleigh scattering (mostly by N2 and CO2 molecules), using the method of Hansen and Travis, (1974).

We note that, for the calculation of absorption coefficients, we used only main isotopes: 12C16O2 and 1H216O and 12N2. It has been shown in similar conditions (Halevy et al.,, 2009) that the radiative effect of isotopic composition should be small. Nonetheless, we take this opportunity to encourage dedicated studies about the influence of isotopic composition on the radiative properties of (thick) atmospheres.

A.3 H2O and CO2 physical properties

Both CO2 and H2O cycles are included in the GCM used in this work.

1. Melting, freezing, condensation, evaporation, sublimation and precipitation of H2O physical processes are all included in the model. In the atmosphere, water vapor can condense into water ice particles clouds. At the surface, we fix the H2O ice albedo at 0.6131313This is the standard value used in the Snowball Model Intercomparison (SNOWMIP) project (Pierrehumbert et al.,, 2011). and we use an emissivity of 1.

2. In our model, CO2 can condense to form CO2 ice clouds and surface frost if the temperature drops below the saturation temperature. Atmospheric CO2 ice particles are sedimented and thus can accumulate at the surface. The CO2 ice layer formed at the surface can sublimate and recycle the CO2 in the atmosphere. The CO2 ice on the surface contributes to the surface albedo calculation: if the CO2 ice layer exceeds a threshold value of 1 mm thickness, then the local surface albedo is set immediately to the albedo of CO2 ice (0.6 in this work). On Mars, the albedo of CO2 ice can substantially vary (Kieffer et al.,, 2000) with insolation and presence of dust. Without dust, the albedo can become very high so that 0.6 is probably a lower estimate. For CO2 ice, we use an emissivity of 0.9. The radiative effect of CO2 ice clouds is discussed in details in section 3.3

Physical parameters used for both CO2 and H2O ices are summarized in table 2.

Appendix B Computation of maximal CO2 ice thickness before basal melting

Compared to water ice, the temperature dependence of CO2 ice thermal conductivity is rather flat in the 100-250 K range (Schmitt et al., 1997; Part I, Figure 4). We can thus estimate with a good approximation the maximum CO2 ice thickness $h_{\text{max}}$ (limited by basal melting) by:

[TABLE]

where $\lambda_{\text{CO}_{2}}$ is the thermal conductivity of CO2 ice, $T_{\text{surf}}$ is the temperature at the top of the glacier and $T_{\text{melt}}$ is the melting temperature of CO2 ice at the base of the glacier.

We derive then $T_{\text{melt}}$ using the Clausius-Clapeyron relation (CO2 solid-liquid equilibrium) at the base of the glacier:

[TABLE]

with $\rho_{\text{sol}}$ and $\rho_{\text{liq}}$ the volumetric mass of liquid and solid CO2, $L_{\text{fus}}$ the latent heat of fusion of CO2 ice, $P_{\text{ref}}$ and $T_{\text{ref}}$ the coordinates of the triple point of CO2, and $P_{\text{CO}_{2}}$ and $P_{\text{N}_{2}}$ the partial surface pressures of CO2 and N2, respectively. The pressure at the base of the glacier $P_{\text{melt}}$ is estimated from the equation $h_{\text{max}}$ = $\frac{(P_{\text{melt}}-P_{\text{surf}})}{g~{}\rho_{\text{sol}}}$, with $P_{\text{surf}}$ = $P_{\text{CO}_{2}}+P_{\text{N}_{2}}$. We remind that we choose pN2 = 1 bar for all our simulations.

We now derive the surface temperature $T_{\text{surf}}$ as a function of the surface pressure $P_{\text{surf}}$, using another Clausius-Clapeyron relation (CO2 solid-gas equilibrium):

[TABLE]

with $L_{\text{sub}}$ the latent heat of sublimation of CO2 ice and $M_{\text{CO}_{2}}$ the molar mass of CO2. It is in fact assumed here that the global mean temperature at the top of the CO2 polar ice caps is constant and equal to the temperature of condensation of CO2. This approximation should remain valid as soon as there is enough CO2 in the atmosphere (or CO2 stays a dominant atmospheric species).

From the set of equations 7-9, we get a new equation on $h_{\text{max}}$ of the form $h_{\text{max}}~{}=~{}e^{h_{\text{max}}}$, after several variable changes. We solve explicitely this equation using the lambert W function, and we obtain the following expression of the maximum possible CO2 ice cap thickness $h_{\text{max}}$:

[TABLE]

This equation is used to directly compute $h_{\text{max}}$ as a function of the internal heat flux and CO2 partial pressure (see Figure 5).

Note that we assumed a CO2 ice volumetric mass density of 1500 kg m*-3* and a CO2 ice thermal conductivity of 0.5 W m${}^{-1}~{}$K*-1* (Schmitt et al., 1997; Part I, Figure 4).

Appendix C Stability of CO2 ice caps

What happens once CO2 ice caps have reached their maximum size, as calculated in section 4.1.4 ? Using calculations, we show in this appendix that, depending on a few parameters, the CO2 ice caps when full may (or may not) be stable.

The total mass of CO2 available at the surface and in the atmosphere is:

[TABLE]

with $A$ the area of the CO2 ice caps, S the total area of the surface (i.e. 4$\pi$R${}_{\text{p}}^{2}$) and $h_{\text{mean}}$ the thickness of the CO2 ice caps averaged over $A$. When CO2 ice caps are full, $h_{\text{mean}}~{}=~{}h_{\text{mean}}^{{\text{max}}}$.

We assume now that the system – initially with full CO2 ice caps – is slightly perturbated from $P_{\text{CO}_{2}}~{}=~{}P$ to $P+\delta P$. By conservation of CO2 mass, the new mass of CO2 ice caps is $M_{\text{surf}}[P+\delta P]$ = $\rho_{\text{CO}_{2}}~{}Ah_{\text{mean}}[P+\delta P]$141414Hereafter, ’[’ and ’]’ are used to bracket variables..

If this quantity is lower than the maximum mass of CO2 ice caps, e.g. $\rho_{\text{CO}_{2}}~{}Ah_{\text{mean}}^{{\text{max}}}[P+\delta P]$, then the system is stable. Otherwise, the condition of unstability can be written formally:

[TABLE]

which can be rewritten, using mass conservation:

[TABLE]

and using Taylor series:

[TABLE]

Assuming that CO2 ice caps were initially full, this yields to $h_{\text{mean}}[P]~{}=~{}h_{\text{mean}}^{{\text{max}}}[P]$ and gives:

[TABLE]

For $\delta P>0$, we have then:

[TABLE]

We know first from GCM simulations that the area of the ice caps $A$ is not a monotonic function of $P_{\text{CO}_{2}}$ (see Figure 2). $A$ is in fact extremely rich in information, because it depends on a subtle combination of the greenhouse effect of CO2, the surface thermal emission, the condensation temperature of CO2 and the global heat atmospheric redistribution. As $P_{\text{CO}_{2}}$ grows, GCM simulations tell us that $\frac{dA}{dP_{\text{CO}_{2}}}$ should vary from positive to negative values. Fortunately, we also know from the various GCM simulations shown in Figure 2b (23.5*∘* obliquity, radiatively inactive CO2 clouds) that the maximum of $A$($P_{\text{CO}_{2}}$) should for example roughly lie at a CO2 partial pressure of 1 bar in this specific configuration.

It is possible to derive an analytical - complicated yet - expression of $\frac{dh_{\text{mean}}^{{\text{max}}}}{dP_{\text{CO}_{2}}}$ using equations 10 and 4. It tells us that, whatever the configuration, $\frac{dh_{\text{mean}}^{{\text{max}}}}{dP_{\text{CO}_{2}}}$ is always negative. As pCO2 increases, the temperature at the top of the CO2 ice caps also increases, which limits the thickness of the CO2 ice caps151515Equation 4 shows that the mean thickness of the CO2 ice caps is proportional to the maximum thickness calculated by basal melting. (see equation 1).

We illustrate now the previous calculations with the experiment extensively described in this work (1.30 AU, 23.5*∘* obliquity, radiatively inactive CO2 ice clouds, CO2 partial pressure of 1 bar). Based on GCM simulations, we assume that $\frac{dA}{dP_{\text{CO}_{2}}}|_{pCO_{2}=1bar}~{}\sim~{}0$ m2 Pa*-1* and we have therefore (with $\frac{A}{S}~{}\sim~{}0.1$; see section 4.1.1):

[TABLE]

Our calculations indicate that this condition is valid for an internal heat flux roughly lower than 80 mW m*-2* (case 1 in Figure 7). It means that, for an internal heat flux lower than 80 mW m*-2*: (1) CO2 ice caps when full would not be stable, and (2) the maximum amount of CO2 that can be trapped in the $\{$surface, atmosphere$\}$ system is attained for pCO2 lower or equal to 1 bar.

Conversely, for an internal heat flux higher than 80 mW m*-2* (case 2 in Figure 7), CO2 ice caps when full would be stable (unless $\frac{dA}{dP_{\text{CO}_{2}}}$ term becomes somehow significant).