Publisher’s note: this comment was edited on 29 August 2024. The following text is not identical to the original comment, but the adjustments were minor without effect on the scientific meaning.

Review of « a simple physical model for glacial cycles » by Sergio Pérez-Montero.

1. General impressions

The contribution of the authors is a very welcome addition to the literature of models of glacial-interglacial cycles. The model that they introduce, called PACCO, features 0-d ice sheet dynamics along with a simple linear predictor of CO_2 concentration. It is argued that it is as an adequate device for  "testing the different hypotheses and isolating them from each other".

1.1 Mechanism attribution

How much we learn about physics and balance of mechanisms in a simple model is an interesting question. I would suggest the following argument. Consider a model with thousands of lines of codes, with the representation of highly detailed mechanisms based on physical principles and parameterisations. Suppose then that this model successfully reproduces one or several glacial cycles. For sure, such achievement required some tuning (the Saltzman 1990 argument), but one might argue that, in such a big model, the relative importance of different components or mechanisms in the resulting ice age dynamics is fairly robust. It does not mean that it is faithful to the truth, but that it is hard to change this balance by tuning. Now consider a low-order model with, say, 10 state variables such as PACCO.  One might suspect that different choices of parameterisations or even different assumptions may have a large impact on the relative role of different components. This is both the strength and the weakness of the low-order models. A strength, because it offers opportunities of experimentation and addresses sensitivities not apparent in the bigger model; a weakness, because as different assumptions might converge to similar mathematical expressions, interpretation is not as straightforward as it may look. This modelling ambiguity was I would argue perhaps most severe in Saltzman - Maash models (e.g. Saltzman and Maash 1990, hereafter SM90), where only one out of the three differential equations was non-linear (the one for the carbon cycle); the carbon cycle was therefore by design given the responsibility of causing the limit cycle.

For sure, PACCO has more physically explicit mechanisms than SM90, but I suspect caution must still be exercised when attributing a "role" to a component --- I will make a case below that perhaps the respective roles of isostasy, basal melt and snow/ice aging might be sensitive to questionable assumptions (through comments 1.2 and 1.3).

Before going further, I would suggest the authors to have a look and perhaps position themselves with respect to two articles which they did not cite, yet  I suspect that have been somehow influential --- because some of the comments I am making here have already been made there,  and also because of snow aging was already the attention of a controversy at the time:

          -        Gallée et al 1992 introduced what they called a "sectorially" averaged climate model, with different parameterisations affecting the heat balance of the ice sheet including snow aging (based on previous works by L.D.D Harvey --- also relevant for your section 2.7)

          -        Tarasov and Peltier 1997 reported a simulation of the latest termination with an energy balance model coupled to an ice sheet model. They commented the Gallée et al. 1992 paper as well as other, similar efforts, and specifically, they made the following point: "Although the addition of each of these plausible mechanisms [they quoted here the list of mechanisms] to the reduced model has led to improvements in their ability to acceptably simulate the ice age cycle, it has remained unclear as to whether difficulties encountered may not be due to an incorrect representation of the processes already included [.].

It could be unfair to take this quote out of context, but I believe the Tarasov and Pelter question is legitimate: if a model includes a number of processes (in PACCO: we have ice sheet dynamics, lithosphere adjustment, snow ageing, parameterization of ice discharge itself obtained depending on a parameterization of the bottom balance), each of these processes being necessarily highly idealised, what robust conclusions can be obtained about the relative importance of each of these processes in the emerging phenomenon which is, in this case, the amplitude, timing and shape of glacial interglacial cycles?

Let us be fair: the authors of PACCO do not go as far as quantifying the relative importance of each mechanism. However, if I read well, they attribute the asymmetry of cycles to the delayed isostatic response, and give a central role to ice aging for the timing of glacial-interglacial cycles. 

Verbitsky et al. 2018 (VCV18) includes many fewer processes of PACCO and so cannot possibly investigate the relative importance of these different mechanisms (as the authors of the present contribution correctly pointed out). For example, it ignores isostasy and snow ageing. For sure, not representing a mechanism does not imply that it plays no role. But nevertheless the question remains: how is it that VCV18 produces glacial interglacial cycles with decent amplitude, timing and shape with only ice flow dynamics, basal temperature feedback and a general, linear climate feedback, while PACCO needs isostasy and snow ageing ? One might argue, VCV18 implicitly represents the isostasy effect in its equations. VCV18 authors were aware of this, and this is the reason why they pitched their paper by emphasizing the role of adimensional numbers, specifically with their attempt to define a "V-number" measuring the ratio between positive and negative feedbacks.

1.2 Mechanism for the delayed feedback

When it comes to the "asymmetry" of cycles, I suspect the crucial element is the presence of a delayed feedback (see, e.g., Quinn et al. 2018 who coded a glacial-interglacial model with delay equations). What component generates delayed feedback ? True, isostasy does. But the thermal balance of the basal ice does, too.  On this point, there is one concern that has already been mentioned in the public discussion by M. Verbitsky. Equation 34 linking the conduction heat flux to the surface air temperature is disputable. The point of controversy here relates to propagation time needed by air temperature signal to reach the bottom layer, which involves an advection timescale.  This effectively introduces a delayed feedback which can be crucial for the generation of a high amplitude glacial cycle with it its asymmetrical shape.

1.3 Fixed-length assumption

Similary, the fixed ice sheet length (L) assumption may significantly distort the albedo feedback and all feedbacks associated to elevation (I think this assumption would affect all dynamical aspects of the model) I'd like to confirm here that I identified these concerns before actually reading Mikhail Verbitsky's comments, though it is not surprising that we both spotted this issue given our common modelling history, and, to be clear, I haven't inspected the discussion that followed the original community comment.

Comments 1.1, 1.2 and 1.3 do not by any means question the legitimacy of the present study but they come as a word of caution about what one should view as the most adequate device for "testing the different hypotheses and isolating them from each other". I would argue that for that aim the device must have the output and resolution needed to be confronted with detailed observations which are independent enough, thus providing independent constraints. Only in this case, I would argue, is mechanism attribution less ambiguous. Besides these concerns (heat flow parameterisation, fixed length, and caveat about mechanism attribution), most comments below can easily be addressed because they are of editorial nature.

Editorial Comments

1. Introduction

- Some comments about the introduction have already been made by Mikhail Verbitsky. Just adding here that Milankovitch did not "postulate" that  glacial - interglacial variability is caused by changes in the insolation. Milankovitch can be seen, as summed up by Berger (2021), as the father of (paleo-)climate modelling. He reasoned on the physics of the radiative balance and effectively computed the variations in temperature to obtain an "effective snow line" (expressed as a latitude) caused by the quasi-periodic changes in insolation. The "postulate" about the astronomical origin of ice ages has to be found earlier (specifically, speaking of the attribution of ice ages to summer insolation , J.J. Murphy 1876. See reviews by Berger 1988 and Berger 2021.)

- The phrase "conceptual model" is widely used, but it may be a bit ambiguous. Perhaps use the distinction between "phenomenological" (describing the observations in a way consistent with the theory but not derived from it --- Paillard 1998 would be phenomenological, for example) and "theoretical" models? In that sense PACCO is indeed a theoretical model with some empirical parameterisations; Verbitsky et al. 2018 also present their model as a theoretical construction. There is no sharp line between phenomenological and theoretical, but I suspect this distinction would help for the presentation of PACCO.

2a. Model description

- Would it be possible to have all parameters listed in one table by alphabetical order? It was at places a bit painful to find the definition and value of the different parameters while reading the text (e.g., after equation 16, or to that \dot s and k_{\dot s} are model parameters, but doesn't mention that T_ref is also listed in the parameter table. It took me a while to find T_thr which indeed is also in table 2.

- Synthetic insolation forcing: I found introducing P_e and then setting it to zero confusing. Why not simply explain that one major difference between the synthetic signal and the real one is the amplitude modulation of precession by eccentricity? (Which indeed can be relevant for resonance phenomena)

- line 300 : I did not understand what is meant by the "Dynamic nature of thermodynamic hypothesis".

2b. Result analysis

line 216 : replace "higher periodicity"  by longer periods.

I suspect (but this can be discussed) that the word "resonance" is used a bit abusively here. Resonance would refer to a disproportionate increase in amplitude for certain frequencies of the forcing (and then distinguishing linear from non-linear resonance, non-linear resonance occurring when the matching frequency depends on the amplitude of forcing). It is not phase locking either because the model does not have a self-sustained and oscillation. So I believe the right phrase here is just non-linear response with the output period being a multiple of that of the forcing.

line 220: "on the contrary, even a moderate increase in insolation induces a termination too easily". (If this is indeed what you meant; though it might be argued, high sensitivity to insolation is also what you need for the MIS12-11 transition).

Line 370: the sentence is a little bit unclear. Isn't it normal to have high sensitivity to sliding?

In the sentence "one possible mechanism could be related to the rigidity of the substratum on which ice is formed" : please clarify or rephrase: rigidity is not a mechanism.

Figure 13: For clarity just specify which h_geo was used for the I,k,m,o panels.

References:

Gallée H., J. P. Ypersele, T. Fichefet, I. Marsiat, C. Tricot and A. Berger (1992), Simulation of the last glacial cycle by a coupled, sectorially averaged climate-ice sheet model. Part II : Response to insolation and CO2 variation, Journal of Geophysical Research, (97) 15713-15740 doi:10.1029/92JD01256

Saltzman B. (1990), Three basic problems of paleoclimatic modeling: a personal perspective and review, Climate Dynamics, (5) 67-78 doi:10.1007/BF00207422

Saltzman B. and K. A. Maasch (1990), A first-order global model of late Cenozoic climate, Transactions of the Royal Society of Edinburgh Earth Sciences, (81) 315-325 doi:10.1017/S0263593300020824

Tarasov L. and R. W. Peltier (1997), Terminating the 100 kyr ice age cycle, Journal of Geophysical Research, (102) 21665-21693, doi:10.1029/97JD01766

Berger A. (1988), Milankovitch theory and climate, Reviews of Geophysics, (26) 624-657

Berger A. (2021), Milankovitch, the father of paleoclimate modeling, Climate of the Past, (17) 1727-1733 doi:10.5194/cp-17-1727-2021

Quinn C., J. Sieber, A. S. Heydt and T. M. Lenton (2018), The Mid-Pleistocene Transition induced by delayed feedback and bistability, Dynamics and Statistics of the Climate System, (None) None doi:10.1093/climsys/dzy005

Authors investigated Pleistocene glacial cycles that are known to be driven by the nonlinear response of the climate to solar forcing. Though complex models accurately simulate these cycles they are computationally expensive. Simpler conceptual models lack physical detail, though. The Physical Adimensional Climate Cryosphere mOdel (PACCO) used in this study aims to balance complexity and simplicity by focusing on the interaction between climate and Northern Hemisphere ice sheets. Authors show that PACCO effectively reproduces 100,000-year glacial cycles by incorporating ice-sheet dynamics, thermodynamics, and ice aging. The study reveals that ice aging and delayed isostatic response are key to matching geological records, making PACCO a valuable tool for studying glacial cycles.

Major comments on the paper:  

I have to say that I like how the paper was written. It goes from a simple model, where linear insolation forcing is introduced to the one-dimensional model as described by SIA equations. I do not agree with RC1/CC1’s description that it even might have 2 dimensions, as here H is not a dimension but merely a product of equations. I only miss here the motivation for all constants and order of complexity they introduced to the model. For example, in table 2, Length is 1000 km based on Oerlemans et al., 2008. But why is simply that taken and not given a reasoning? The same goes for many equations that are not fully motivated. Please note that I am a big fan of simple conceptual models as I believe that we can learn more from simple yet fast models compared to expensive and long-time running Stokes models. 

Then, from section 2.2 onwards, authors speak of periodicities at time scales of 60-80-100-120 kyears. I must disagree as, e.g., Figure 5 shows pronounced peeks around 20 and 40 kyears for all values of the sliding parameter, and around 60 kyears for lower values of the sliding parameter. Above that, there is no significant peak which might suggest that what we see is just a model artifact or an aliasing error. And this is additionally supported with the figure 5c. The same holds for figure 7 (ISOS simulation). When we go to figure 9 (RISOS simulation), the authors statement does not hold at all, because the peaks really depend on the sliding parameter and cannot be generalized. So here, for me, the result of AGING simulation is the only result that can claim that there is a periodicity between 80-120 kyears (the correct time really depends on the sliding parameter).

Another thing that should be changed in the following manuscript is explanations and captions of all figures. For example: we can see a change in dominant peaks which is logical as the authors introduce the model complexity. However, it is not explained why the energy of the peaks for the same sliding factor changes between the runs (see for example the difference in peaks at 60 kyears between Figure 5 and Figure 7. Also, here the vertical axes would really be beneficial to make a good comparison.

Additionally, the authors did explain why we see differences in periodicity between the sliding parameters, and why their statements are “true” for the lower values of the sliding parameter. There to say, for which glacier could we use these results to possibly predict the future evolution? 

Thus, to conclude: I think the paper is a valuable contribution to understand the glacier cycles and ice-sheet dynamics. However, I would not consider it for publication just yet. 1. The model and its equations need to be explained and introduced better as as CC1 says: “it looks like a verbalized computer code”. 2. The figures are not “self-standing” and need a better explanation of all lines, better comparison between experiments, and y-axes would be beneficial. 3. Better discussion is necessary as it is not obvious how this model can be used and what is a future of it.