Review of "Tracing the Snowball bifurcation of Aquaplanets through time reveals a fundamental shift in critical-state dynamics" by Georg Feulner, Mona Sofie Bukenberger and Stefan Petri

Reviewer: Aiko Voigt

The authors study the location of the Snowball Earth bifurcation in terms of atmospheric CO2 as a function of insolation in the range of 1361-1034 Wm-2. As the sun becomes stronger over time, the insolation range covers the time from today to 3600 Ma before present, meaning that the work studies the bifurcation as a function of time. The authors apply a model of intermediate complexity with a simplified atmosphere model in aquaplanet setup, which allows them to sweep through a broad range of insolation and CO2 values. Their two main findings are i) that for lower insolation values the critical CO2 decreases logarithmically as insolation increases but drops faster for higher insolation values, and ii) that the nature of critical states (defined as states before the runaway ice-albedo feedback sets in) is different between low and high insolation. For low insolation values, the critical ice edge is located in the midlatitudes (termed the Ferrel state by the authors), whereas for higher insolation values it is located in the subtropics (the Hadley state). The authors ascribe this difference in critical states to the meridional gradient in insolation and wind-driven sea-ice transport. The text is well written and the graphics are of high quality (except for two minor questions, see below).

My main criticisms is the following. From reading the text it seems the authors suggest that critical states with a sea ice cover around 50% or with a sea ice edge equatorward of 30 deg were not possible. Yet, there are several studies that have found such states. The conclusion of the change in critical state dynamics thus seems not as robust as described by the authors. I also found that in some cases the comparison with previous studies seems a bit lopsided. I elaborate this below as part of my main comments.

Overall, however, this is a well conducted and well presented paper that addresses a question that was so far not studied. I am confident the authors can address my concerns and recommend minor revisions.

Main comments:

1. In the conclusion section (L350ff) the authors argue that critical states with a sea ice cover around ~40% are not possible (the exact numbers are model dependent). The argument is made based on the Ferrel vs. Hadley states, and is allegedly supported by comparison to the work of Yang et al. However, when checking the figures in Yang et al. (2012a) I believe I found some inconsistencies with the authors' arguments. Specifically, Fig. 2 of Yang shows that there are stable states with a sea ice fraction of 50%, contradicting the statement that "...  there are no stable states with global sea-ice fractions between ∼ 40% and ∼ 60% for a present-day continental configuration." Probably even more severe, Fig. 16b of Yang et al. (2012a) shows that there is a stable state with 70% sea ice cover for 90% insolation. In my understanding such a state contradicts the Ferrel-Hadley-state argument of the authors. There might be other inconsistencies with the Yang et al results.

2. I am missing a discussion about the fact that critical states with sea ice margins quite close to the equator have been found in models, e.g., Voigt and Abbot (2012), Abbot et al. (2011, http://dx.doi.org/10.1029/2011JD015927) and Braun et al. (2022, https://doi.org/10.1038/s41561-022-00950-1). Overall, this makes me think that the changes in the critical state dynamics - although operating in the Climber model used here - are not as robust and fundamental as described by the authors.

Other comments:

L10 and L145: Is the change in the CO2-insolation function related to the change in the critical state dynamics? This is not clear to me.

L27: It is unclear to me what you mean by "for even lower solar luminosities". What does "even" refer to.

L80: Pierrehumbert et al., 2011 (doi:10.1146/annurev-earth-040809-152447) compared Snowball initiation in three AGCMs in aquaplanet setup (their Fig. 4). These models did not include ocean and sea ice dynamics, but used the same coordinated setup. Also, Hoerner et al, JAMES, 2022 (https://doi.org/10.1029/2021MS002734) used an aquaplanet setup to study the impact of sea ice thermodynamics on Snowball initiation. Maybe these are interesting references?

L101: Some more discussions on the atmosphere model, its limitation and the impacts of its limitations would be desirable. For example, are the Hadley and Ferrel cell boundaries fixed in time, or can they move with the seasonal cycle? How does this impact the P-E patterns and hence snow on sea ice and surface albedo? Do the authors think that this matters? This would also be helpful for the wind argument made around L262 in the result section.

L111: The agreement with the Liu et al (2013) work seems cherry picking and in my view is a weak argument. There are other studies for which the agreement would be much lower, as in fact can be seen from Fig. 1 of the paper.

Table 1: I would find it helpful if the S/S0 ratio could be included in the table, as the ratio is used in Figs. 1 and 2.

L140 and L193: The ~0ppm CO2 value for today's insolation is consistent with Voigt and Marotzke, 2010, who found that removing all CO2 would lead to a Snowball in the coupled ECHAM5/MPI-ESM model (using present-day continents).

L147ff: I do not understand what the authors mean by baseline warming from water vapor. I also wonder how clouds are treated in Climber.

L162: Voigt et al., 2011, Climate of the Past showed that moving continents to the tropics cools the climate and facilitates Snowball initation. This is in line with the argument made by the authors and maybe worth including.

L175: I agree with the statement that sea ice dynamics was found to facilitate Snowball initiation. Yet I do not agree that previous studies robustly found that simplified oceans make Snowball initiation more difficult. There are at least three counter examples. Poulsen and Jacob (2004, doi:10.1029/2004PA001056) stated that "The wind-driven ocean circulation transports heat to the sea-ice margin, stabilizing the sea-ice margin.". Rose (2015, https://doi.org/10.1002/2014JD022659) also found a stabilizing role of ocean heat transport. This relates to the argument made in L215 regarding the lack of a full ocean. Voigt and Abbot (2012, https://doi.org/10.5194/cp-8-2079-2012) show explicitly that setting ocean heat transport to zero makes Snowball initiation easier, and they argue that this is related to the subtropical wind-driven ocean cells (see their Figs. 12 and 13).

L180: The study of Pierrehumbert et al., 2011 (see above) tested for albedo values in 3 models, showing that ice albedo differences are key.

L198: I believe Lewis et al., 2003 used prescribed surface winds, because of which they could not make robust statements of the impact of sea ice dynamcics. See the discussion of the Lewis work in Voigt and Abbot (2012; page 3 left column).

L261: Is the fuzzy transition a result of seasonal averaging over fully ice covered grid boxes or does the model allow for partially ice covered boxes?

L274: I am wondering about the role of the wind-driven subtropical ocean cells below the Hadley cells. These cells should be represented by the ocean model and are expected to work towards Snowball initiation (see my comment regarding L175).

Fig. 1: I appreciate the very nice summary of previous modeling work in the figure. Some relevant studies seem to be missing, however. I suggest adding the results of Pierrehumbert et al. (2011), Voigt and Abbot (2012), Hoerner et al. (2022, https://doi.org/10.1029/2021MS002734) and Braun et al. (2022, https://doi.org/10.1038/s41561-022-00950-1). I apologize that these are all studies that I co-authored, I am listing them here since they are missing and I know of them. There might be additional relevant work.

Are Figs. 5 and 6 needed given the zonal symmetry and the zonal-mean plots in Fig. 7?

Fig. 8: I do not understand the meaning of the legend in panel a and the color coding of the lines.

Feulner et  al. traced the snowball bifurcation of an aquaplanets through the history of the Earth using a climate model of intermediate complexity. To my knowledge, this has never been done before. Importantly, they did get some very interesting results from such practice and found that the critical climate state was different under low and high solar constants. In particular, they found that once the sea-ice edge crosses ~40° latitude, it would march forward to the equator without any external forcing when the solar constant is low. While when the solar constant is high, this critical latitude is at ~30° latitude, i.e. the boundary of the Hadley cell. Therefore, I think the work provides new knowledge about the stability of the Earth’s climate to the society and definitely worth publication. However, there are still a few relatively small things that need to be clarified. Especially, the mechanism for the stability of the ‘Hadley state’ may be explained better.

The major reason that the ice edge cannot be stabilized at ~30° latitude when the solar constant is low, I think, is because the atmosphere+ocean heat transport across the 30° latitude exceeds the energy that can be received by the oceans equatorward. This will cause a continuous cooling of the tropical region and eventually allows the ice edge to march forward towards the equator. The rate of this cooling should be inversely proportional to the solar constant and is indeed well indicated by their Fig. 4. Therefore, I hope they can demonstrate this mechanism more clearly by showing explicitly the total meridional heat transport at 30°S and 30°N and the solar energy received by the oceans within 30° latitudes. These should be shown for the transient stage in one of the simulations, for example, at year 100 of the 1500 Ma simulation in Fig. 4.

If the point above can be confirmed, then the mechanism for stabilizing the ‘Hadley state’ may need to be modified (such as the last sentence of the abstract). It is stabilized by enhanced oceanic heat transport once the sea-ice edges approach the boundary of Hadley cell. The enhanced heat transport is expected to be due to stronger easterly winds as normally seen in other models and thus stronger poleward Ekman transport. This mechanism always works as clearly shown in Fig. 4 but it can stabilize the climate only momentarily when solar constant is low because the enhanced heat transport extracts more energy than the tropical ocean can receive; the heat content of the tropical ocean is drained out quickly. While for a high solar constant, a balance can be achieved easily (other feedback processes are naturally involved, especially the outgoing longwave radiation, latent heat flux etc. so that the tropical ocean will lose less energy in these ways) unless the CO₂ concentration is lowered further. This is likely also the major reason that the slope in Fig. 1 increases once such balance can be achieved.

Minor questions:

• How are the poles treated in the ocean module of this model since an aquaplanet is simulated?
• Is the boundary of the Ferrel cell also fixed in this model?
• Can the oceanic and atmospheric heat transport be calculated separately in this model?
• L319: ore -> or

Yonggang Liu