the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 25 Jun 2024
A simple physical model for glacial cycles
Abstract. Glacial cycles are the norm in Pleistocene climate variability. Models of varying degree of complexity have been used to answer the question of what causes the nonlinear response of the climate system to the periodic forcing from the Sun. At one end of the spectrum of complexity are comprehensive models which aim to represent all involved processes in a realistic manner. However, their high computational cost precludes their use in the ultra-long simulations needed. At the other end are conceptual models which are computationally less demanding but which generally lack a physical basis. Most of them yield good results in terms of capturing the shape and patterns of glacial cycles as indicated by the geological record, thus making it very difficult to identify the underlying mechanisms. Here we present a conceptual model that aims to physically represent the interaction between the climate and the Northern Hemisphere ice sheets while eliminating spatial dimensions in some of the fundamental ice-sheet thermodynamic and dynamic equations. To this end, we describe the Physical Adimensional Climate Cryosphere mOdel (PACCO) from its simplest to its most complex configuration. We discuss separately the implications of different fundamental mechanisms such as ice-sheet dynamics and thermodynamics, glacial isostatic adjustment and ice-sheet albedo aging for our model. We conclude that ice-sheet dynamics and a delayed isostatic response are sufficient to produce resonance around periodicities of 100 kyr, although the forcing has a spectrum concentrated around lower values. In addition, ice-sheet thermodynamics and ice aging separately enhance the model nonlinearities to provide 100 kyr periodicities in good agreement with reconstructions. However, we found that it is easier to reproduce the late Pleistocene glacial cycles using the simpler process of ice aging. Overall, PACCO is a valuable tool for analyzing the different hypotheses present in the literature.
- Preprint
(2940 KB) - Metadata XML
- BibTeX
- EndNote
Status: final response (author comments only)
-
CC1: 'Comment on egusphere-2024-1842', Mikhail Verbitsky, 02 Jul 2024
Simple but physics-based models of ice ages have been introduced many years ago (e.g., Birchfield and Weertman, 1978; Chalikov and Verbitsky, 1984, etc.); nevertheless, phenomenological models of glacial rhythmicity are, as Michel Crucifix (2012) said, still “seductive”. Their obvious lack of underlying physics is compensated by introduction of artificial thresholds and other Boolean statements that, indeed, allow them to reproduce empirical data but do not add much to our understanding of the physical nature of glacial periods. Therefore, every attempt to introduce “A simple physical model for glacial cycles” should be welcomed.
Unfortunately, the presented paper has (I frankly hesitate to say that, but for the lack of a better word…) a flaw that needs to be addressed before one proceeds to the results.
You approximate ice discharge as q = v H/L, where L is a constant. Even intuitively, making the horizontal size of evolving Northern Hemisphere Pleistocene ice sheets to be a constant, is strange, but from the basics of ice physics it is simply incorrect: H and L are not independent, and H is proportional to L^(1/2) (e.g., Verbitsky, 1992, Bahr et al, 2015). Taking this into account may dramatically change the dynamical properties of your governing equation (1).
Few comments which are less significant, but may be you consider them helpful:
- Your model is not “adimensional” as you claim in your PACCO abbreviation. For example, your equations have dimensional H, measured in meters, and dimensional time, measured in seconds. You didn’t attempt to make your system adimensional. You simply do not resolve space.
- It would be very helpful, if you present your final equations all together in one place instead of forcing a reader like me to do this work and substitute let say tau_d into v_d into v into q into dH/dt. At this point, the model description looks like a verbalized computer code (that I suspect it is).
- Since you position your model as a physical model, it would be nice if you explain the basic physical nature of the governing equations. For example, I was able to figure out that your equation (1) is scaled kinematic condition on the ice sheet free upper surface, but I am not sure I can explain with the same confidence equation (29) especially when Q is not defined.
- Line 35: “However, most of them rely on very mathematical approaches and include artificial or imposed thresholds and trends (Paillard, 1998; Paillard and Parrenin, 2004; Gildor and Tziperman, 2001; Verbitsky et al., 2018; Ganopolski, 2024).” First, what “very mathematical” means? And, second, I am not aware about “artificial or imposed thresholds and trends” in Verbitsky et al (2018).
- And finally, a bit of funny thing. Your introduction begins with the sentence that makes a reader to believe that Paillard in 2001 and Ganopolski in 2024 introduced glacial-interglacial variability to the world. With all due respect to celebrated scientists, I would suggest someone like Agassiz to be mentioned there. It would also help to explain (next sentence) how Milankovitch was able to offer his theory before them.
Respectfully,
Mikhail Verbitsky
References
Bahr, D. B., Pfeffer, W. T., and Kaser, G.: A review of volume-area scaling of glaciers, Rev. Geophys., 53, 95–140, doi:10.1002/2014RG000470, 2015.
Birchfield, G.E. and Weertman, J., 1978. A note on the spectral response of a model continental ice sheet. Journal of Geophysical Research: Oceans, 83(C8), pp.4123-4125.
Chalikov, D.V. and Verbitsky, M.Y., 1984. A new Earth climate model. Nature, 308(5960), pp.609-612.
Crucifix, Michel. "Oscillators and relaxation phenomena in Pleistocene climate theory." Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 370, no. 1962 (2012): 1140-1165.
Verbitsky, M.Y. Equilibrium ice sheet scaling in climate modeling. Climate Dynamics 7, 105–110 (1992). https://doi.org/10.1007/BF00209611
Citation: https://doi.org/10.5194/egusphere-2024-1842-CC1 -
CC2: 'Additional considerations regarding ice-sheet thermodynamics', Mikhail Verbitsky, 05 Jul 2024
I think I owe you a more explicit explanation of my discomfort with Equation (29).
It is a common knowledge that for typical ice sheets the Peclet number (Pe) is of order of 10. This means that temperature advection dominates heat diffusion and an ice-flow trajectory has a near-constant temperature determined by its value on the upper free surface of the ice sheet (e.g., Grigorian et al, 1976, Morland, 1984, etc). The thickness of the basal boundary layer where, instead, the heat diffusion begins to dominate, is proportional to Pe^(-1/2)*H and is about 100 m. The timescale of the upper-surface temperature “delivery” to the basal boundary layer is the same as the timescale of ice growth. Your equation (29) seems to describe the heat balance of such basal boundary layer. Its thickness H_b = 10 m implies that you assume Pe to be even larger than 10. Nevertheless, the mechanism of delayed cold ice delivery to the bottom layer is absent in equation (29) and replacing it with conduction term (34) is very difficult to justify. Obviously, the absence of the vertical-temperature-advection timescale may have significant implications for the entire system dynamical properties.
Respectfully,
Mikhail Verbitsky
References
Grigoryan, S. S., M. S. Krass, and P. A. Shumskiy. "Mathematical model of a three-dimensional non-isothermal glacier." Journal of Glaciology 17, no. 77 (1976): 401-418.
Morland, L. W. "Thermomechanical balances of ice sheet flows." Geophysical & Astrophysical Fluid Dynamics 29, no. 1-4 (1984): 237-266.
Citation: https://doi.org/10.5194/egusphere-2024-1842-CC2 - AC2: 'Reply on CC2', Sergio Pérez Montero, 19 Jul 2024
-
AC1: 'Reply on CC1', Sergio Pérez Montero, 19 Jul 2024
Dear Mikhail Verbitsky,
Thank you for your comments. We respond to all your concerns (CC1 and CC2) in the attached pdf document. Here, we summarize the main changes that we consider to be implemented in the document:
- We will add a variable for the typical horizontal spatial scale of the ice sheet following your suggestion for the driving stress (Equation 10). However, we will keep as a constant when calculating ice discharge (Equation 3).
- We will clarify some terms and add some references that we consider relevant thanks to your comments.
- We will modify the thermal balance in the basal boundary layer equations.
Finally, we would like to emphasize again that even though these comments are useful and have changed a few things in the manuscript (requiring recalibration in some cases), our results are qualitatively very similar and thus the conclusions are robust.
We hope our answers clarify your concerns.
Best regards,
Sergio Pérez-Montero et al.
-
CC3: 'Reply on AC1', Mikhail Verbitsky, 20 Jul 2024
Dear Sergio Pérez-Montero, Dear Colleagues,
My response is enclosed as a pdf file.
Best regards,
Mikhail Verbitsky
-
AC3: 'Reply on CC3', Sergio Pérez Montero, 31 Jul 2024
Dear Mikhail Verbitsky,
Thanks for your comment, you can find the answer to this comment in the attached pdf.
Regards,Sergio Pérez-Montero et al.
-
CC4: 'Reply on AC3', Mikhail Verbitsky, 01 Aug 2024
Dear Sergio Pérez-Montero, Dear Colleagues,
I appreciate your response. I think we can finally arrive to a common vision. My final comment is just a suggestion.
Certainly, you have every right to add a bit of phenomenology in your otherwise physical model and parameterize the q-term as being proportional to vH because you observe this in the 3-D model results. However, to preserve physical “purity” of your q ∼ vH term, I would suggest the reasoning based on CC3.1 as it has been derived from the mass conservation law. Indeed,
q ∼ L_ocean* vH/L^2 (CC3.1)
If we suggest that L_ocean (ice sheet boundary exposure to the ocean) increases with increased ice sheet area, i.e. L_ocean ∼ S ∼ L^2, then CC3.1 becomes:
q∼ vH
This may be another substantiation of your parameterization.
I wish you successful publication,
Mikhail Verbitsky
Citation: https://doi.org/10.5194/egusphere-2024-1842-CC4
-
CC4: 'Reply on AC3', Mikhail Verbitsky, 01 Aug 2024
-
AC3: 'Reply on CC3', Sergio Pérez Montero, 31 Jul 2024
-
RC1: 'Comment on egusphere-2024-1842', Michel Crucifix, 03 Aug 2024
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.
- 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
- 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
Citation: https://doi.org/10.5194/egusphere-2024-1842-RC1 - AC4: 'Reply on RC1', Sergio Pérez Montero, 13 Nov 2024
-
RC2: 'Comment on egusphere-2024-1842', Anonymous Referee #2, 26 Aug 2024
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.
Citation: https://doi.org/10.5194/egusphere-2024-1842-RC2 - AC5: 'Reply on RC2', Sergio Pérez Montero, 13 Nov 2024
Video supplement
Representation of the Late Pleistocene using the AGING experiment from PACCO. Sergio Pérez-Montero https://github.com/sperezmont/Perez-Montero-etal_YYYY_ESD/blob/a172bb2654cf525d4e5096365489514b0f62641f/figures/pacco_animation_s1_a1.gif
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 539 | 198 | 39 | 776 | 14 | 16 |
- HTML: 539
- PDF: 198
- XML: 39
- Total: 776
- BibTeX: 14
- EndNote: 16
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1