Orbital and millennial-scale forcing of the Patagonian Ice Sheet throughout the Last Glacial Cycle

Domain: earth_science
Published: 2025-10-02
Editorial outcome: minor_revision

Abstract
During the Last Glacial Maximum (23,000 to 19,000 years ago), the Patagonian Ice Sheet covered the central chain of the Andes between 38°S and 55°S. The paleoclimatic evidence from Patagonia and New Zealand suggests that the maximum glacier expansion of the Southern Hemisphere mid-latitudes was desynchronized with the Northern Hemisphere glacial history. Here we present numerical simulations of the Patagonian ice sheet throughout the Last Glacial Cycle. Our analysis suggests that the Patagonian ice sheet had two main periods of advance, during the Marine Isotope Stage 4 and late Marine Isotope Stage 3, experiencing inter-millennial scale variability. We show that the Patagonian Ice Sheet long-term evolution can be attributed to changes in the integrated summer insolation, which combines the summer duration and insolation intensity and has an obliquity-like periodicity. We further suggest that this metric also modulated the behaviour of glaciers over the entire Southern Hemisphere mid-latitudes.

Introduction
During the last glacial period, particularly during the Last Glacial Maximum (LGM, ~23,000 to 19,000 years ago1), the Earth hosted extensive ice sheets beyond those in Antarctica and Greenland. Notably, North America, northern Europe, and Patagonia were covered by the North American Ice Sheet complex, the Eurasian Ice Sheet complex, and the Patagonian Ice Sheet (PIS), respectively2. As the waxing and waning of ice masses is largely controlled by variations in temperature and precipitation, past ice sheets provide important information on past climate changes.

The PIS extended between approximately 38°S and 55°S, and recent geochronological advances have refined our understanding of its extent and chronology3. As such, studies of the evolution of the PIS provide unique insights into the Southern Hemisphere mid-latitudes climate system, for which there is a notorious lack of proxy information4. At the LGM, the PIS covered an estimated area and volume comparable with the size of the former British-Irish Ice Sheet that covered the British Isles during the LGM and was part of the Eurasian Ice Sheet complex3,5,6. Notably, based on geochronological reconstructions, it has been proposed that the PIS reached its maximum extent at around 35ka, during the Marine Isotope Stage (MIS)3. Relatively stable conditions followed until 30ka when it started to retreat until a temporary stabilization at around 20ka before the final deglaciation3. Although the PIS extension over the last deglaciation has been relatively well constrained using moraine dating and glacial erratics, large uncertainty remains in both the pre-LGM evolution and the dynamical processes of the ice sheet related to millennial-scale variability, which are not resolved by geochronological reconstructions3,7–9.

The temporal evolution of the PIS, in particular its asynchronous history when compared to the ice masses in the Northern Hemisphere during the global LGM, is still a matter of debate4,10–12. In general, paleoclimate proxy reconstructions suggest an extended LGM in the Southern Hemisphere, where the onset of maximum glacial conditions commenced from about 35ka, even though this pattern is subject to spatial heterogeneity and millennial-scale variability13. The asynchronous timing of the Southern Hemisphere mid-latitude ice fields and glaciers, and particularly of the PIS, has been linked to the latitudinal migration and strength of the southern westerly winds10,14,15 as well as to the expansion of the Antarctic sea-ice, suggesting a dominant role of obliquity16. More recently, the equatorial migration and strengthening of the southern westerly winds, linked to global cooling and a steeper meridional temperature gradient, have been proposed as main drivers for the glacial expansion, as well as interhemispheric teleconnections4.

Simulating the spatio-temporal evolution of the PIS over the Last Glacial Cycle (LGC, ~120,000 years ago to present day) is notoriously difficult due to the lack of continuous climate forcing fields17. However, a recently published alkenone-based sea surface temperature (SST) record from marine sediment core MR16-09 PC03 provides us with an opportunity to reconstruct the Patagonian climate over the past 140ka as the core is situated about 150 km offshore the Taitao Peninsula in southern Chile/central Patagonia (Fig. 1), and has a temporal resolution high enough to resolve millennial-scale variability9.

Here we present results from transient model simulations of the PIS throughout the complete LGC using the numerical ice sheet model SICOPOLIS18 and making use of the MR16-09 PC03-based climate reconstructions. To this end, we employ a glacial index approach using two different alkenone-based temperature proxies UK37 and UK'37. The two indices differ in alkenone C37:4, which occurs mainly in cold regions. Both alkenone records show strong orbital-scale changes, suggesting glacial-interglacial sea surface temperature differences between 5 and 8°C (Fig. 1) and reproducing coherent chronologies with the glacial histories of New Zealand and South America, suggesting the presence of large-scale common mechanisms9,10,14,19. At the core site, the choice of the index is critical, in particular since millennial-scale variability associated with cold Antarctic stadials is well pronounced in the UK37 record, but almost lacking in the UK'37 record (Fig. 2)9. Therefore, the use of both proxy records in our study allows us to elucidate the role of millennial-scale Antarctic cold stadials in the dynamics of the PIS, even though Hagemann et al.9 argue that UK37 is preferable at the core location.

The goal of this study is to explore the spatio-temporal history of the PIS throughout the LGC, assessing the timing and dynamics of glacial advances and retreats. In particular, we investigate the role of millennial-scale climate variability in forcing and response and revisit the orbital forcing mechanisms of the PIS.

Results
The simulation forced by the UK37 record suggests an ice inception at around 115ka (Fig. 2b). During the MIS5, the ice conditions remained relatively stable, entering a quasi-steady state. By approximately 70ka, coinciding with the onset of MIS4, the model reconstructs a significant increase in the ice extent. The growth of the ice sheet during this period resembles the extent reached during the LGM. This growth is punctuated by a millennial-scale variability, associated with two prominent cold pulses, which facilitated ice sheet expansion, with a notable increase around 65ka (Fig. 2a). However, these pulses did not solely contribute to the ice sheet’s overall growth but instead marked a plateau of sustained glaciation through the late MIS4 to early MIS3, stabilizing near 40ka. During the late MIS3 and the MIS2, millennial-scale variability has then triggered cold pulses that have promoted the ice sheet growth, reaching the maximum ice extent during the LGC. During this period, the results suggest that the ice sheet has experienced waxing and waning rather than a smooth glacial history (Fig. 2b) as opposed to the earlier suggestions based on geological reconstructions3.

The simulation driven by UK'37 reveals a distinctive trajectory, with rapid ice sheet growth during MIS5, reaching ice extent comparable to the LGM. Notably, around 95ka, the model exhibits a marked and abrupt reduction in ice extent, followed by a rapid recovery at approximately 90ka. This decrease is associated with a warming signal in sea surface temperature, which is much less pronounced in the UK37 record (Fig. 2a). During the MIS4, the simulation performed with UK'37 showcases quasi-stable ice extent behavior, with only minor fluctuations, and a nearly absent millennial-scale variability as opposed to UK37-based simulation. The impact of the millennial-scale variability on the ice sheet geometry and thickness is mainly over the north- and south-western margins as well as over the eastern margin in the vicinity of MR16-09 PC03 (Supplementary Fig. 6). A modest growth occurs during the MIS4-MIS3 transition (Fig. 2b), but during MIS3 and MIS2, the ice sheet remains in a relatively stable condition, again indicating that UK'37 does not capture millennial-scale variability during the MIS2, suggesting minimal variations in the ice sheet geometry.

In summary, the simulation with the UK37 forcing suggests two main periods of glacier advance, one during the MIS4 and another during the late MIS3-early MIS2, with millennial-scale variability throughout the interval with the maximum ice extent. The modeled PIS is mainly drained by fast-flowing outlet glaciers to the east, while even larger velocities can be found in the western ice margin (Supplementary Figs. 3 and 5). The velocities of ice streams that terminate in the Pacific Ocean range over 1000 ma−1, whereas the velocities of land-terminating fast-flowing outlet glaciers at the eastern margins of the PIS are generally within 100 ma−1. The ice divide mainly extends north-south along the Patagonian Andes between. On the other hand, the simulation performed with UK'37 shows less variability, representing the LGM as a phase of stability with limited dynamical shifts in the ice sheet configuration (Fig. 3).

Methods
Icesheet model and climate forcing
We apply the open-source, three-dimensional, thermomechanical ice sheet model SICOPOLIS (SImulation COde for POLythermal Ice Sheets18) for the area between 80°W and 62°W and between 36°S and 58°S. The ice model setup was presented in a previous publication17 and is briefly summarized here. Supplementary Table 1 contains an overview of the applied model parameters.

The inception and evolution of the PIS in our model are driven by the surface mass balance (SMB), which is calculated as the difference between the applied fields of accumulated precipitation and surface ablation. The latter is computed using a positive-degree-day (PDD) model31, based on the given near-surface air temperature field. PDD parameters have been selected based on contemporary and paleo studies in the area17,26,32. Surface mass accumulation is assumed to depend on monthly precipitation and temperature fields. The transition between solid and liquid precipitation is linearly proportional to variations in air temperature with 0°C to 2°C limiting purely solid or purely liquid, respectively33. As the model domain surface evolves, discrepancies between the prescribed (fixed) topography used in the climate model snapshots and the dynamic one in SICOPOLIS are accounted for by implementing a near-surface air temperature lapse-rate correction of −6.5 K km−1. The precipitation changes by 7.3% for each degree Celsius of air temperature change34. Glacial isostatic adjustment of this bedrock produced by temporal variations in the ice mass load is accounted for through an elastic lithosphere-relaxing asthenosphere (ELRA) model35 using standard parameter values. To isolate the atmospheric dynamics and the ocean dynamics, the model is set as such that it calves as soon as it reaches the coast.

Aiming to explore the glacial history of the PIS before the global LGM, we perform transient simulations of the PIS throughout the full LGC. To generate a climate state at any given model time, LGM and PI conditions from MPI-ESM1-2-LR have been used36, representing peak glacial and interglacial conditions. Based on previous studies, MPI-ESM1-2-LR stands as the best-fitted paleoclimate model in Patagonia for equilibrium and transient studies17,26. These two snapshots are then subject to a weighted interpolation following a glacial index approach with

GI(t) = X(t) − XPI / XLGM − XPI (1)

where t is time and X is the SST derived from the offshore records37. The time-dependent climate-forcing fields are then given by

T(t) = TPI + GI(t)(TLGM − TPI) (2)

P(t) = PPI (1 − GI(t)) (1 − PLGM / PPI) (3)

where T and P represent the temperature and precipitation fields through time, respectively.

In the present study, the time-dependent weight GI is derived from the offshore sediment core MR16-09 PC03. The core covers the past 140ka, being retrieved at around 46°S, 77°W from a water depth of 3082m9. SST records were derived from two different alkenone-based proxies from the same core, UK37 and UK'37. The computation of GI uses the 19–23ka mean of SST to define peak glacial conditions. Likewise, peak interglacial conditions are defined near the PI period by averaging over the last 3ka. Each simulation is initiated from ice-free conditions.

Calculation of integrated summer energy
The integrated summer energy29 has been calculated using the daily mean insolation data from ref. 38. The integrated summer energy was calculated assuming ISE = ∑ i Wi ⋅ 86400), with Wi being the daily mean insolation and i equal to one when Wi is greater than or equal to a threshold τ and equal to zero otherwise29.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "Review of “Orbital and millennial-scale forcing of the Patagonian Ice Sheet throughout the Last Glacial Cycle”\nAndrés et al present a modeling study of the Patagonian Ice Sheet (PIS) through the last glacial cycle that explores why the\nPIS reached it maximum likely before the global LGM, and to a lesser extent, why this is out of phase with Northern\nHemisphere ice sheets. This study adds to an expansive body of literature over the last several decades that have all\nfocused on 1) when were glaciers/ice sheets in the Southern Hemisphere expanded, 2) when did they retreat, 3) what forced\nthese changes, and 4) how/why is this signature different from the Northern Hemisphere.\nI am not an ice-sheet modeler, but I am a Quaternary paleoclimatologist/geochronologist who thinks about these questions\nabove. As such, I found this paper quite interesting and it had me thinking about my own research and sifting through papers\nand furthermore, I think this manuscript is generally well written. I honestly cannot comment if this manuscript is “correct”, but\nthe methods seem sound, the discussion reasonable, and I do think it a worth addition to the broader discussion surrounding\nSouthern Hemisphere glaciation. I only have a few suggestions and questions:\n1) At first, I struggled a bit with “What is the main selling point here?”, which frankly is probably not a good sign for a\nmanuscript like this. I think in the end the main selling points are a) a scheme that uses integrated insolation as a driver of\nlong-term change provides the best data-model fit, and 2) more traditional treatments of insolation are unable to explain the\npattern of PIS glaciation. I encourage the authors to really emphasize integrated insolation. For example, its mentioned in\nthe abstract, but the first thing I read (and thought was the main takeaway) was that the PIS advanced during MIS/3…to\nwhich my response was “we’ve known that for awhile.”\n2) Selection of sediment core MR16-09 PC03. This seems a bit arbitrary, and I do wonder about how reproducible these\nresults are using other temp records. I strongly suspect this record was used because it’s the longest record from the region\nand spans the last glacial cycle (and probably the most recent record?). But it may be worth exploring how other records\naffect the results. ODP sites 1233 and 1234 slightly north of the record used here also have alkenone SST records. I believe\nat least one of those also spans the entirety of the last glacial. How easily can these be incorporated. At the bare minimum, I\ndo not think these records can be completely ignored….are they different/similar to MR16-09 PC03? Should other records\nbe plotted on a figure as well?\n3) Like #2, I find it hard to believe I just read an entire Southern Hemisphere glaciation paper without seeing a single ice-\ncore record plotted, let alone used. The authors routinely mention that “millennial-scale Antarctic cold stadials played a\ncritical role….” I think the authors nee to at the bare minimum show these on a figure somewhere. If I’m playing devil’s\nadvocate, what would happen if you simulated the last glacial cycle using a temp record from Antarctica, or at least WAIS\n(realizing that WAIS does not go back 140ka)? Has this been done before? Maybe what I’m also getting at here is how do\nthe authors simulations compare to previous efforts. Since this is Nature Communications I think you have a little more\nleeway with article length.\n4) I like how the authors made a decent attempt to bring in the geologic record. I think many modeling groups often don’t\neven bother, so this is a nice touch. My only comment on this aspect would be how the authors evaluate model-data success\nor mismatches. There are a lot of references to how the author’s efforts *broadly* do a good job at reproducing the geological\nrecord of PIS change, but it is a bit tough to gauge from the text and figures. I will say that having the 35-ka outline in Figure\n3 is very helpful although I was wondering a bit about the slight mismatch between model output and the geologic record at\n35 ka. Might be nice for the authors to comment on that. Also, why does the model output appear shifted east at 35 ka? Or\nmaybe another way to put it is that it undershoots ice on the west side but overshoots on the east side. I think this paper\nwould be stronger if the comparison between model output and the geological record were a bit more robust, maybe there is\nsome way to quantify this. The most obvious targets/benchmarks are MIS 4 (PIS at max near max as constrained by\ngeologic record) and MIS 2 where there are numerous dated moraines.", "reviewer_2": "This manuscript is one of the first to attempt to model the Patagonian Ice Field / Ice Sheet through the Last Glacial Period\nusing a new ocean core that provides a alkenone-based sea surface temperature record (MR16-09 PC03) from Hagemann\net al. (2024). The ability to use such a record is really interesting and fascinating, as such attempts have previously not been\nable to do so. I do commend the manuscript for attempting to do so and being at the forefront of modelling the Patagonian Ice\nField using such recent and new data.\nNoteworthy results\n- Modelling the Patagonian Ice Field through the Last Glacial Period itself is noteworthy\n- Multiple advances of the Ice Field both within MIS 4 and 3 that are varied due to inter-millennial variability\n- That, using the records provided, the simulations agree with PATICE (Davies et al., 2020), but disagree with a smooth\ndeglaciation but rather varied deglaciated with interannual fluctuations\n- Summer and winter isolation do not provide a clear explanation for ice variability, but more closely follow summer energy.\nDoes the work supports the conclusions?\nThe evidence provided in Figure 2, the model output, and the comparison with other records provides clear evidence for the\nconclusions presented within the paper. Would have been nice to have seen more of a comparison with more outlines of the\nPatagonian Ice Sheet, to see it through different timesteps, and if the model with the new SST record could match more then\njust the 35ka transition.\nFlaws in data analysis:\nWhile the output itself of the ice extent and thickness are valuable in and of themselves, I think more could bet said from your\nmodel which you have negated to do so. It would be really interesting and useful to understand what is happening wit the\nvelocity fields, how the ice is evacuating itself through the ice sheet. Does this change as the ice sheet involves through the\nLast Glacial Period? Is the thickness realistic, and does it match what other models have done over Patagonia for the LGM\nfor example? While it may capture 35ka well, does it actually match other time periods well?\nMethodology:\nThis modelling is a continuation of a modelling studies from the same authors over Patagonia (Castillo-Llarena et al., 2024).\nWhile much of the methodology is similar, there are certainly some differences which should be noted in the methodology.\nFor example, how the input data was made and set up for the model. What are the parameters used in the model? Any\ndifferences from the default parameter values? This could be in the form of a table so we know the parameters varied, as\nfrom reading the former paper before this, there is also no clear indication of parameter values. Also, if there was any\nadditional parameters or modules used to enable the model to be ran through to completion. At the moment, I do not think\nthere is enough for replication purposes.\nAttached is also a commented version of the manuscript for minor points.\nOverall comment and review verdict:\nI personally think the amount of information that could be extracted from such a study is not inline with the reduced amount of\ninformation a nature paper allows the author. While the study and findings are some of the first to have been conducted over\nPatagonia, the outputs need to be described in-depth, with further information also to be presented in the methodology. The\nstudy itself is highly commendable and a pleasure to read, but I do not see Nature Comms being the best fit for such a paper.\nIf the authors wish to resubmit to another journal, I would be more then happy to be put down as a reviewer to provide further\nfeedback.", "reviewer_3": "This study produces a new reconstruction of the Patagonian Ice Sheet (PIS) across the last glacial cycle, extending a\nprevious reconstruction which largely covered MIS3. The authors relate changes in their modeled PIS to integrated summer\ninsolation and millennial-scale climate variability. The ice sheet reconstruction is forced by a continuous SST record over\nthe past 140kya from the eastern Pacific offshore of Chile, which records variability in the western extent of the PIS. This\nprovides insight into regional climate dynamics over a significantly longer period than previous nearby records. This SST\nrecord is used to construct a glacial index which guides the interpolation of glacial and pre-industrial climate states to be\nused as model forcing.\nThe results of this study provide important constraints on the history of the PIS beyond what was previously available and\nreveals new millennial-scale dynamics within previously studied intervals. I congratulate the authors on the work that went\ninto this study, and it will undoubtedly provide a great resource to gain further insight into the PIS.\nMy main points of discussion focus on the capability to meaningfully differentiate between the two forcings across MIS5,\nwhich in my view is weakened by the lack of independent verification and uncertainty in the forcing. I will discuss this below,\nalong with the authors three primary findings which I take directly from the abstract:\n(1) “Our analysis suggests PIS advances during the Marine Isotope Stage (MIS) 4 and late MIS3...” This finding is dependent\non the absolute temperature of the Uk37 record, which is the primary difference, rather than high-frequency variability,\nbetween the two forcings across MIS5.\n(2) “We show that millennial-scale climate variability played a role in the PIS dynamics…” This is an important observation\nfor this study, as it shows that the model can sustain a dynamic PIS in-phase with millennial-scale climate variability.\n(3) “…while the long-term PIS evolution can be attributed to changes in the integrated summer insolation” This is also an\nimportant observation. However, it is unclear if this observation is limited to this reconstruction or is a consistent feature\nacross all of the Southern Hemisphere.\nOverall, this study provides valuable new information surrounding the history of the PIS. A strength of this study is the\nidentification of millennial-scale PIS variability across MIS3, which warrants further description. In my view, this paper should\nprovide more conservative claims surrounding the MIS5-4 history and additional insight into the dynamics of newly observed\nmillennial-scale variations.\nMain Discussion Points:\n1) My main discussion point is the impact of the of glacial index in the PIS reconstruction. This is particularly evident in the\ndifference in extent during MIS5, where the reconstruction based on Uk’37 rapidly expands at ~115 kya. I presume this is\ndue to the systematically lower temperatures in Uk’37 as compared to Uk37.\nHow well understood is the temperature uncertainty using these SST reconstruction methods? Neither this study, nor the\ndata source Hagemann et al., 2024, discuss this. The difference between records is ~1.25-1.5oC, which is in the range of\nabsolute uncertainty in comparable reconstruction methods. Based on the uncertainty in the SST’s, can the different PIS\nreconstruction scenarios be differentiated?\nWhile I agree with the author’s interpretation that other lines of evidence (primarily the MAR/biomarker data of Hagemann et\nal., 2024 and the glacial advancement markers cited within Hagemann et al., 2024) suggest that the PIS reconstruction\nbased on Uk’37 can be excluded, I feel that discussion of uncertainty within the SST forcing and the sensitivity of the model\nto that forcing is a necessary discussion in order to justify an advancement across MIS4 rather than MIS5. As it stands, the\nresults only suggest that PIS initiation occurred by ~115ka and reached its approximate LGM extent by at least the onset of\nMIS4.\n2) For the observed millennial-scale variations in extent and volume, where do these changes manifest? I feel this is key\nsource of information which is not discussed. Are these changes in extent isolated to specific regions, or are they well\ndistributed across the PIS? This could be included as a supplemental spatial anomaly map.\nHow do the changes in this new reconstruction compare to other reconstructions like PATICE or that of Darvill et al., 2016?\nDoes the difference in extent in Figure 3 in the 35ka panels imply this reconstruction underestimates the geological\nevidence?\n3) In a recent study (Castillo-Llarena et al., 2024), Antarctic ice core data is used as a forcing over the last 70 kya. Has this\nbeen done for comparison over the last glacial cycle here? This could provide a necessary source of independent\nverification for the orbital-scale variability of the PIS.\nThe structure and magnitude across MIS3 is also consistent between the EDC forced reconstruction of Castillo-Llarena et\nal., 2024 and Uk37 shown here. Including a comparison between the two could provide insight into their similarities and\ndifferences across both millennial and orbital timescales. Furthermore, how do the other reconstructions of\nCastillo-Llarena et al., 2024 which utilize other SST records across MIS3/2 compare to the ones presented here? If the PIS\nadvancement during late MIS3 is a regional phenomenon, is it expressed across the wider SST forcings and in their respective\nice sheet reconstructions?\n4) In a similar vein, is statement in lines 136-138: “The local minimum in integrated summer energy during the beginning of the\nMIS4 aligns well with the timing of the PIS growth, while the positive anomaly in early MIS3 helps explain the gradual\ndisintegration of the ice sheet (Fig. 2b).” further supported by previous reconstructions of the PIS?\n5) Does this ice sheet reconstruction support glacial sediment delivery to the site? The western margin is reduced as\ncompared to previous reconstructions, with more significant advances across the eastern margin. Does this reduced western\nextent still result in marine terminating conditions, or ice extent across the continental margin as suggested by\nHagemann et al? In other words, is the PIS reconstruction presented here consistent with the conditions necessary to reproduce the\nsediment forcing record?\nMinor Comments:\nLine 87-88: References Figure 2e here.\nLine 99-101: Does the Be10/Be9 contradict the results of PIS advancement at the start of MIS4? The Be10/Be9 is indicative\nof terrigenous sedimentation beginning at ~80kya, with oscillations preceding it. The record across MIS5 is of a similar\nmagnitude across MIS3, despite the differing ice extents in the Uk37 reconstruction.\nLine 104-105: Which records of Hagemann et al., (2024) suggest that Uk37 is the preferred record? The Uk37 SST record is\nused as the forcing, so that is circular. Is it the MAR record? If so, I recommend being specific about which records provide\nsupport.\nLine 115: Is the PATICE reconstruction meant to be shown in Figure 2b?"}, "author_response": "Final Response to Reviewer #1:\nReview of “Orbital and millennial-scale forcing of the Patagonian Ice Sheet throughout the\nLast Glacial Cycle”\nAndrés et al present a modeling study of the Patagonian Ice Sheet (PIS) through the last\nglacial cycle that explores why the PIS reached it maximum likely before the global LGM,\nand to a lesser extent, why this is out of phase with Northern Hemisphere ice sheets. This\nstudy adds to an expansive body of literature over the last several decades that have all\nfocused on 1) when were glaciers/ice sheets in the Southern Hemisphere expanded, 2)\nwhen did they retreat, 3) what forced these changes, and 4) how/why is this signature\ndifferent from the Northern Hemisphere.\nI am not an ice-sheet modeler, but I am a Quaternary paleoclimatologist/geochronologist\nwho thinks about these questions above. As such, I found this paper quite interesting and it\nhad me thinking about my own research and sifting through papers and furthermore, I think\nthis manuscript is generally well written. I honestly cannot comment if this manuscript is\n“correct”, but the methods seem sound, the discussion reasonable, and I do think it a worth\naddition to the broader discussion surrounding Southern Hemisphere glaciation. I only have\na few suggestions and questions:\n1) At first, I struggled a bit with “What is the main selling point here?”, which frankly is\nprobably not a good sign for a manuscript like this. I think in the end the main selling points\nare a) a scheme that uses integrated insolation as a driver of long-term change provides the\nbest data-model fit, and 2) more traditional treatments of insolation are unable to explain the\npattern of PIS glaciation. I encourage the authors to really emphasize integrated insolation.\nFor example, its mentioned in the abstract, but the first thing I read (and thought was the\nmain takeaway) was that the PIS advanced during MIS/3…to which my response was “we’ve\nknown that for awhile.”\nResponse: We are pleased with the positive response and the valuable feedback provided\nby Reviewer 1. The integrated summer energy enables to explain not only the earlier\nmaximum extension of the PIS to the detriment of more traditional treatment of insolation\nsuch as the peak summer insolation, peak winter insolation or summer duration, but also\nexplaining the glacial advance over the MIS4. Additional info to support the integrated\nsummer energy as the best model-data fit has been included in the supplementary\ninformation material. Changes to stress the main findings have been implemented in the\nabstract as suggested.\n2) Selection of sediment core MR16-09 PC03. This seems a bit arbitrary, and I do wonder\nabout how reproducible these results are using other temp records. I strongly suspect this\nrecord was used because it’s the longest record from the region and spans the last glacial\ncycle (and probably the most recent record?). But it may be worth exploring how other\nrecords affect the results. ODP sites 1233 and 1234 slightly north of the record used here\nalso have alkenone SST records. I believe at least one of those also spans the entirety of\nthe last glacial. How easily can these be incorporated. At the bare minimum, I do not think\nthese records can be completely ignored….are they different/similar to MR16-09 PC03?\nShould other records be plotted on a figure as well?\nResponse: The selection of the sediment core MR16-09 PC03 was due to being the longest\nand most recent record from the region with a frequency sampling that allows for the study of\nthe PIS variability over the Last Glacial Cycle with thorough accuracy. Additionally, it is\nlocated in the center of the former latitudinal extension of the PIS. However, we\nacknowledge that the discussion of other records is beneficial to enhance the robustness of\nthe findings. The orbital scale variability has imprints not only in MR16-09 PC03 but also in\ndifferent cores over the South Pacific. The marine records located north of the study zone\nODP-1233 and ODP-1234 have imprints of the orbital scale variability, suggesting two main\nperiods of relative minimum temperature: during the MIS4 and late MIS3 and early MIS2,\nboth periods agree in the timings of the orbital-scale variability induced by the ISE. South of\nMR16-09 PC03, the core MD07-3128 does not cover the MIS4, therefore, limited information\ncan be extracted. Far off-shore records covering the time period here studied reveal that\nminimum sea-surface temperatures earlier than LGM are captured, however they have a\nlimitation in the sampling frequency. Additionally, EDML also has imprints of such variability.\nFigures 1 and 2 of the main text were modified, including the location of the offshore records\n18\nmentioned above and by including the δ 𝑂 record of EDML, respectively. Additionally, a\nfigure with the temperature reconstruction is now shown in the supplemental material\n(Supplementary Figure 2). An earlier study (Castillo-Llarena et al., 2024) has shown have\nconcluded that the timings of maximum glacial advance and the consequent asynchronicity\nwith the global mean behaviour are not captured by Antarctic cores, and local proxy data are\nneeded to capture the local variability (Castillo-Llarena et al., 2024). However, imprints of the\norbital-scale variability are indeed captured. A discussion on the imprints in other offshore\nrecords in the southeastern Pacific Ocean has been included.\n3) Like #2, I find it hard to believe I just read an entire Southern Hemisphere glaciation paper\nwithout seeing a single ice-core record plotted, let alone used. The authors routinely mention\nthat “millennial-scale Antarctic cold stadials played a critical role….” I think the authors nee\nto at the bare minimum show these on a figure somewhere. If I’m playing devil’s advocate,\nwhat would happen if you simulated the last glacial cycle using a temp record from\nAntarctica, or at least WAIS (realizing that WAIS does not go back 140ka)? Has this been\ndone before? Maybe what I’m also getting at here is how do the authors simulations\ncompare to previous efforts. Since this is Nature Communications I think you have a little\nmore leeway with article length.\nResponse: As commented in the response of #2, insights of the Orbital Scale variability are\nbeing visible not only in the core initially considered in this study but also in further offshore\nrecords analysed and the Antarctic ice core record EDML. We thank the reviewer for this\n18\nsuggestions that help the interpretation of the results. The δ 𝑂 record from EDML and the\nAntarctic Isotope Maxima are now shown in Figure 2. The influence of the millennial-scale\nvariability on the core record MR1609 PC03 has been largely discussed by Hagemann et al.\n(2024), and it has been highlighted now in Figure 2 to facilitate the discussion. Previous\nmodelling efforts focusing on the PIS have shown that the Antarctic cores failed to reproduce\nthe timing of the maximum extent of Patagonia during the late MIS3, highlighting the\nimportance of considering local records to capture the variability of the PIS with enough\naccuracy (Castillo-Llarena et al., 2024).\n4) I like how the authors made a decent attempt to bring in the geologic record. I think many\nmodeling groups often don’t even bother, so this is a nice touch. My only comment on this\naspect would be how the authors evaluate model-data success or mismatches. There are a\nlot of references to how the author’s efforts *broadly* do a good job at reproducing the\ngeological record of PIS change, but it is a bit tough to gauge from the text and figures. I will\nsay that having the 35-ka outline in Figure 3 is very helpful although I was wondering a bit\nabout the slight mismatch between model output and the geologic record at 35 ka. Might be\nnice for the authors to comment on that. Also, why does the model output appear shifted\neast at 35 ka? Or maybe another way to put it is that it undershoots ice on the west side but\novershoots on the east side. I think this paper would be stronger if the comparison between\nmodel output and the geological record were a bit more robust, maybe there is some way to\nquantify this. The most obvious targets/benchmarks are MIS 4 (PIS at max near max as\nconstrained by geologic record) and MIS 2 where there are numerous dated moraines.\nResponse: We are pleased with the positive feedback provided by Reviewer 1 with regard to\nthe geologic records. On the one hand, our simulations are able to reproduce the timing of\nthe two main periods of glacial advance proposed by the available literature, during the MIS4\nand during the late MIS3. On the other hand, our simulations are also capable of simulating\na maximum ice extent earlier than the LGM. Despite MPI-ESM standing as the\nbest-performing climate model (Yan et al., 2022; Castillo-Llarena et al., 2024), mismatches\nbetween the reconstructed and modeled ice extent have been found not only in this research\nwork but also in previous works, being largely attributed to the coarse resolution of the\natmospheric forcing and the consequent prescribed topography (Castillo-Llarena et al.,\n2024). Coarse and inaccurate orographic forcing in Patagonia induces wetter and colder\ntemperatures over the leeside of the Andes (Damseaux et al., 2019) promoting the ice extent\nbeyond its geologically reconstructed margin. We have compared the total modelled area\nwith the existing PIS geochronological reconstruction PATICE (Davies et al., 2020) showing\na good agreement despite the mismatch with the geological shape due to the atmospheric\nforcing of the paleoclimate model used. The total area of PATICE for each time slide has\nbeen added to figure 2 for comparison.\nReferences:\nCastillo-Llarena, A., Retamal-Ramírez, F., Bernales, J., Jacques-Coper, M., Prange, M., &\nRogozhina, I. (2024). Climate and ice sheet dynamics in Patagonia throughout marine\nisotope stages 2 and 3. Climate of the Past, 20(7), 1559-1577.\nDamseaux, A., Fettweis, X., Lambert, M., & Cornet, Y. (2019). Representation of the rain\nshadow effect in Patagonia using an orographic‐derived regional climate model.\nInternational Journal of Climatology.\nDavies, B. J., Darvill, C. M., Lovell, H., Bendle, J. M., Dowdeswell, J. A., Fabel, D., ... &\nThorndycraft, V. R. (2020). The evolution of the Patagonian Ice Sheet from 35 ka to the\npresent day (PATICE). Earth-Science Reviews, 204, 103152.\nHagemann, J. R., Lamy, F., Arz, H. W., Lembke-Jene, L., Auderset, A., Harada, N., ... &\nTiedemann, R. (2024). A marine record of Patagonian ice sheet changes over the past\n140,000 years. Proceedings of the National Academy of Sciences, 121(12), e2302983121.\nYan, Q., Wei, T., & Zhang, Z. (2022). Modeling the climate sensitivity of Patagonian glaciers\nand their responses to climatic change during the global last glacial maximum. Quaternary\nScience Reviews, 288, 107582.\nFinal Response to Reviewer #2:\nThis manuscript is one of the first to attempt to model the Patagonian Ice Field / Ice Sheet\nthrough the Last Glacial Period using a new ocean core that provides a alkenone-based sea\nsurface temperature record (MR16-09 PC03) from Hagemann et al. (2024). The ability to\nuse such a record is really interesting and fascinating, as such attempts have previously not\nbeen able to do so. I do commend the manuscript for attempting to do so and being at the\nforefront of modelling the Patagonian Ice Field using such recent and new data.\nNoteworthy results\n- Modelling the Patagonian Ice Field through the Last Glacial Period itself is noteworthy\n- Multiple advances of the Ice Field both within MIS 4 and 3 that are varied due to\ninter-millennial variability\n- That, using the records provided, the simulations agree with PATICE (Davies et al., 2020),\nbut disagree with a smooth deglaciation but rather varied deglaciated with interannual\nfluctuations\n- Summer and winter isolation do not provide a clear explanation for ice variability, but more\nclosely follow summer energy.\nDoes the work supports the conclusions?\nThe evidence provided in Figure 2, the model output, and the comparison with other records\nprovides clear evidence for the conclusions presented within the paper. Would have been\nnice to have seen more of a comparison with more outlines of the Patagonian Ice Sheet,\nto see it through different timesteps, and if the model with the new SST record could match\nmore then just the 35ka transition.\nResponse: We are pleased with the positive feedback provided by Reviewer 2. As\nsuggested, PATICE total area has been incorporated into Figure 2 for comparison. Additional\nfigures showing our model results and its comparison with PATICE at 35, 30, 25, 20 and 15\nka outlines with both proxies (supplementary figures 4 and 5). Moreover, figures showing the\nice sheet velocities for 80, 65, 57, 35 and 21 ka can now be found in supplementary figure 3.\nFlaws in data analysis:\nWhile the output itself of the ice extent and thickness are valuable in and of themselves, I\nthink more could bet said from your model which you have negated to do so. It would be\nreally interesting and useful to understand what is happening wit the velocity fields, how the\nice is evacuating itself through the ice sheet. Does this change as the ice sheet involves\nthrough the Last Glacial Period? Is the thickness realistic, and does it match what other\nmodels have done over Patagonia for the LGM for example? While it may capture 35ka well,\ndoes it actually match other time periods well?\nResponse: We thank the reviewer for this thoughtful and constructive comment. Numerical\nreconstructions of the Ice Thickness of the Patagonian Ice Sheet diverge considerably, being\na matter of debate between the existing literature. Ice sheet thickness estimations are\nranging between 1500 and 3200 m (Sudgen et al., 2002; Yan et al., 2022; Wolf et al., 2023;\nCastillo-Llarena et al., 2024). Those discrepancies could be attributed to model resolution,\ninput data and parametrizations. Our results regarding the ice sheet thickness are within the\nexisting range of values proposed by previous studies. Additional figures showing the\ncomparison between our model simulations and the PATICE reconstruction at 35, 30, 25, 20\nand 15 ka have been included (Supplementary Figure 4). Figures with the ice velocity for the\nsame time slices are included in the Supplemental Material (Supplementary Figure 5). The\nmodeled PIS is mainly drained by fast-flowing outlet glaciers to the east, while even larger\nvelocities can be found in the western margin (Supplementary Figure 3 and 5). The\nvelocities of ice streams that terminate in the Pacific Ocean range over 1000 m/year,\nwhereas the velocities of land-terminating fast-flowing outlet glaciers at the eastern margins\nof the PIS are generally within 100 m/year. The ice divide mainly extends north-south along\nthe Patagonian Andes between. Variations in the ice velocities can be found throughout the\ndifferent time steps, particularly focused in the margins of the ice sheet (Supplementary\nFigure 3 and 5).\nMethodology:\nThis modelling is a continuation of a modelling studies from the same authors over\nPatagonia (Castillo-Llarena et al., 2024). While much of the methodology is similar, there are\ncertainly some differences which should be noted in the methodology. For example, how the\ninput data was made and set up for the model. What are the parameters used in the model?\nAny differences from the default parameter values? This could be in the form of a table so\nwe know the parameters varied, as from reading the former paper before this, there is also\nno clear indication of parameter values. Also, if there was any additional parameters or\nmodules used to enable the model to be ran through to completion. At the moment, I do not\nthink there is enough for replication purposes.\nResponse: The model setup is the same as in Castillo-Llarena et al. (2024). The parameters\nused to perform the simulations shown in this paper are now included in the supplemental\nmaterial of the manuscript (Supplementary table 1). More detailed description of the model\nconfiguration has been included into the method section of the revised manuscript.\nAttached is also a commented version of the manuscript for minor points.\nResponse: All comments have been addressed.\nOverall comment and review verdict:\nI personally think the amount of information that could be extracted from such a study is not\ninline with the reduced amount of information a nature paper allows the author. While the\nstudy and findings are some of the first to have been conducted over Patagonia, the outputs\nneed to be described in-depth, with further information also to be presented in the\nmethodology. The study itself is highly commendable and a pleasure to read, but I do not\nsee Nature Comms being the best fit for such a paper. If the authors wish to resubmit to\nanother journal, I would be more then happy to be put down as a reviewer to provide further\nfeedback.\nWe sincerely thank the reviewer for their thoughtful and encouraging feedback. We are\nparticularly grateful for the kind remarks regarding the novelty of our work over Patagonia, as\nwell as the acknowledgment that the study was a pleasure to read and scientifically\ncommendable. We remain confident that Nature Communications is a suitable venue for this\nwork due to its relevance to global climate processes, the originality of the dataset, and the\nimplications for understanding ice sheet behavior outside polar regions. We believe that the\nmain message of the paper can be effectively conveyed within the constraints of Nature\nCommunications format.\nReferences\nCastillo-Llarena, A., Retamal-Ramírez, F., Bernales, J., Jacques-Coper, M., Prange, M., &\nRogozhina, I. (2024). Climate and ice sheet dynamics in Patagonia throughout marine\nisotope stages 2 and 3. Climate of the Past, 20(7), 1559-1577.\nWolff, I. W., Glasser, N. F., Harrison, S., Wood, J. L., & Hubbard, A. (2023). A steady-state\nmodel reconstruction of the patagonian ice sheet during the last glacial maximum.\nQuaternary Science Advances, 12, 100103.\nSugden, D. E., Hulton, N. R., & Purves, R. S. (2002). Modelling the inception of the\nPatagonian icesheet. Quaternary International, 95, 55-64.\nYan, Q., Wei, T., & Zhang, Z. (2022). Modeling the climate sensitivity of Patagonian glaciers\nand their responses to climatic change during the global last glacial maximum. Quaternary\nScience Reviews, 288, 107582.\nFinal Response to Reviewer #3:\nThis study produces a new reconstruction of the Patagonian Ice Sheet (PIS) across the last\nglacial cycle, extending a previous reconstruction which largely covered MIS3. The authors\nrelate changes in their modeled PIS to integrated summer insolation and millennial-scale\nclimate variability. The ice sheet reconstruction is forced by a continuous SST record over\nthe past 140kya from the eastern Pacific offshore of Chile, which records variability in the\nwestern extent of the PIS. This provides insight into regional climate dynamics over a\nsignificantly longer period than previous nearby records. This SST record is used to\nconstruct a glacial index which guides the interpolation of glacial and pre-industrial climate\nstates to be used as model forcing.\nThe results of this study provide important constraints on the history of the PIS beyond what\nwas previously available and reveals new millennial-scale dynamics within previously\nstudied intervals. I congratulate the authors on the work that went into this study, and it will\nundoubtedly provide a great resource to gain further insight into the PIS.\nMy main points of discussion focus on the capability to meaningfully differentiate between\nthe two forcings across MIS5, which in my view is weakened by the lack of independent\nverification and uncertainty in the forcing. I will discuss this below, along with the authors\nthree primary findings which I take directly from the abstract:\n(1) “Our analysis suggests PIS advances during the Marine Isotope Stage (MIS) 4 and late\nMIS3...” This finding is dependent on the absolute temperature of the Uk37 record, which is\nthe primary difference, rather than high-frequency variability, between the two forcings\nacross MIS5.\n(2) “We show that millennial-scale climate variability played a role in the PIS dynamics…”\nThis is an important observation for this study, as it shows that the model can sustain a\ndynamic PIS in-phase with millennial-scale climate variability.\n(3) “…while the long-term PIS evolution can be attributed to changes in the integrated\nsummer insolation” This is also an important observation. However, it is unclear if this\nobservation is limited to this reconstruction or is a consistent feature across all of the\nSouthern Hemisphere.\nOverall, this study provides valuable new information surrounding the history of the PIS. A\nstrength of this study is the identification of millennial-scale PIS variability across MIS3,\nwhich warrants further description. In my view, this paper should provide more conservative\nclaims surrounding the MIS5-4 history and additional insight into the dynamics of newly\nobserved millennial-scale variations.\nResponse: We sincerely thank the reviewer for their thoughtful and constructive comments,\nas well as for their positive assessment of our study. We are particularly grateful for the\nrecognition of the effort invested in developing this new reconstruction of the Patagonian Ice\nSheet and for highlighting the value of our work in extending the temporal framework and\nidentifying orbital and millennial-scale variability. We appreciate the reviewer’s suggestions\nregarding a more cautious interpretation of the MIS5-4 interval and the need for further\nelaboration on the dynamics observed during MIS3. In response, we have revised the\nrelevant sections to better reflect the uncertainties associated with MIS5-4 and have\nexpanded our discussion of the millennial-scale signals during MIS3, as detailed below.\nMain Discussion Points:\n1) My main discussion point is the impact of the of glacial index in the PIS reconstruction.\nThis is particularly evident in the difference in extent during MIS5, where the reconstruction\nbased on Uk’37 rapidly expands at ~115 kya. I presume this is due to the systematically\nlower temperatures in Uk’37 as compared to Uk37.\nHow well understood is the temperature uncertainty using these SST reconstruction\nmethods? Neither this study, nor the data source Hagemann et al., 2024, discuss this. The\ndifference between records is ~1.25-1.5oC, which is in the range of absolute uncertainty in\ncomparable reconstruction methods. Based on the uncertainty in the SST’s, can the different\nPIS reconstruction scenarios be differentiated?\nWhile I agree with the author’s interpretation that other lines of evidence (primarily the\nMAR/biomarker data of Hagemann et al., 2024 and the glacial advancement markers cited\nwithin Hagemann et al., 2024) suggest that the PIS reconstruction based on Uk’37 can be\nexcluded, I feel that discussion of uncertainty within the SST forcing and the sensitivity of the\nmodel to that forcing is a necessary discussion in order to justify an advancement across\nMIS4 rather than MIS5. As it stands, the results only suggest that PIS initiation occurred by\n~115ka and reached its approximate LGM extent by at least the onset of MIS4.\nResponse: The simulations have been performed by using a transient forcing derived from a\nglacial index method. The method takes as reference values the pre-industrial conditions\nand the Last Glacial Maximum (LGM, 23,000 to 19,000 years before present). During the\nLGM, the average temperature of Uk’37 is around 1.5ºC warmer than for Uk37. Such a\ndifference is a consequence of the millennial-scale variability captured by Uk37, inducing\nmuch more extreme conditions over the LGM and the late MIS3. During MIS5, the\ntemperatures derived from Uk’37 are close to those derived by the core during the LGM,\ntherefore, the forcing during those periods mimics the LGM forcing. The glacial index method\nand therefore the modelled ice sheet behaviour are independent of the uncertainty in the\ncomparison of the temperature records and the offsets of both proxies.\n2) For the observed millennial-scale variations in extent and volume, where do these\nchanges manifest? I feel this is key source of information which is not discussed. Are these\nchanges in extent isolated to specific regions, or are they well distributed across the PIS?\nThis could be included as a supplemental spatial anomaly map.\nHow do the changes in this new reconstruction compare to other reconstructions like\nPATICE or that of Darvill et al., 2016? Does the difference in extent in Figure 3 in the 35ka\npanels imply that this reconstruction underestimates the geological evidence?\nResponse: Our simulations are performed using the glacial index method, which enables\ngenerating climate states at any given time through the weights interpolation of the LGM\nand PI snapshots of the best-fitted paleoclimate model in Patagonia based on previous modelling\nstudies (Yan et al., 2022; Castillo-Llarena et al., 2024). The imprints of both the orbital and\nthe millennial-scale variability are present in the core record used to create the glacial index\ncurve. While we reckon the limitations of the method which assumes that the perturbations\nare applied homogeneously over the full study zone, the method has been largely used for\npaleoclimate studies in Patagonia and elsewhere (eg., Charbit et al., 2007; Niu et al., 2019,\nNiu et al., 2021; Mas e Braga et al., 2021; Castillo-Llarena et al., 2024). The marks of the\nmillennial-scale variability are present in the vicinity of the offshore records, in the western\nand eastern margins of the PIS, particularly north of 44ºS and south of 50ºS (Supplementary\nFigure 6).\nThe obtained geometry of the maximum extension of PIS is then closely linked with the near\nsurface air temperature of the cold member of the glacial index, in this case, the full LGM\nconditions of the MPI-ESM1-2-LR, which due to the coarse atmospheric model resolution\nand the topographic misrepresentation leads to mismatches with the geochronologically\nreconstructed geometry of the PIS (Castillo-Llarena et al., 2024). Matching the simulated ice\nsheet with the geochronological reconstruction under different timesteps is not the main\nfocus of this study but to understand the overall evolution of the ice mass over the last\nGlacial Cycle and why its asynchronous behaviour with the northern hemisphere ice sheets.\n3) In a recent study (Castillo-Llarena et al., 2024), Antarctic ice core data is used as a forcing\nover the last 70 kya. Has this been done for comparison over the last glacial cycle here?\nThis could provide a necessary source of independent verification for the orbital-scale\nvariability of the PIS.\nThe structure and magnitude across MIS3 is also consistent between the EDC forced\nreconstruction of Castillo-Llarena et al., 2024 and Uk37 shown here. Including a comparison\nbetween the two could provide insight into their similarities and differences across both\nmillennial and orbital timescales. Furthermore, how do the other reconstructions of\nCastillo-Llarena et al., 2024 which utilize other SST records across MIS3/2 compare to the\nones presented here? If the PIS advancement during late MIS3 is a regional phenomenon, is\nit expressed across the wider SST forcings and in their respective ice sheet reconstructions?\nResponse: The core MR16-09 PC03 was selected due to being the longest and most recent\nrecord from the regions with a frequency sampling that allows the study of the variability over\nthe Last Glacial Cycle accurately. However, we acknowledge that the discussion of other\nrecords and moreover those used in previous modelling efforts is beneficial in enhancing the\ndiscussion and robustness of the findings. Existing modelling efforts of the PIS have\nconcluded that the timings of maximum glacial advance and the consequent asynchronicity\nwith the global mean behaviour is not captured by Antarctic cores, and local proxy data is\nneeded to capture the local variability (Castillo-Llarena et al., 2024). However, imprints of the\norbital-scale variability are indeed captured. An analysis comparing different SST\nreconstructions based on other off-shore records over the Last Glacial Cycle was included.\nAlso, the results of Castillo-Llarena et al. (2024) have been discussed. The orbital-scale\nvariability has imprints not only in MR16-09 PC03 but also in different cores over the south\nPacific (Supplementary Figure 2). The near offshore records located north of the study zone\nODP-1233 and ODP-1234 have imprints of the orbital scale variability, suggesting two main\nperiods of relative minimum temperature: during the MIS4 and during late MIS3 and early\nMIS2, both periods agree in the timings of the orbital-scale variability induced by the ISE.\nSouth of MR16-09 PC03, the core MD07-3128 does not cover the MIS4, therefore, limited\ninformation can be extracted. The far offshore records GeoB3327-5 and PS75-034/2,\ndespite their limitation in the sampling frequency, reveal that minimum sea-surface\ntemperatures earlier than LGM are captured with a variability in line with the ISE.\nAdditionally, EDC also has imprints of such variability.\n4) In a similar vein, is statement in lines 136-138: “The local minimum in integrated summer\nenergy during the beginning of the MIS4 aligns well with the timing of the PIS growth, while\nthe positive anomaly in early MIS3 helps explain the gradual disintegration of the ice sheet\n(Fig. 2b).” further supported by previous reconstructions of the PIS?\nResponse: Periods of glacial advance are broadly in line with glacial advance suggested by\nprevious research efforts (Hagemann et al., 2024; Sproson et al., 2024). Additional\nindependent verification of the timings can be found based on geological records for\nPatagonia (Moreno et al., 2015; Garcia et al., 2018; Garcia et al., 2021; Glasser et al., 2011;\nHein et al., 2010; Kaplan et al., 2008; Peltier et al., 2021; Mendelova et al., 2020) and the\nsouthern hemisphere mid latitudes outside Patagonia (Schaefer et al., 2015; Hall et al.,\n2020; Rudolph et al., 2020; Lesic et al., 2022). Additional offshore records have been\nincluded in the supplemental material.\n5) Does this ice sheet reconstruction support glacial sediment delivery to the site? The\nwestern margin is reduced as compared to previous reconstructions, with more significant\nadvances across the eastern margin. Does this reduced western extent still result in marine\nterminating conditions, or ice extent across the continental margin as suggested by\nHagemann et al? In other words, is the PIS reconstruction presented here consistent with the\nconditions necessary to reproduce the sediment forcing record?\nResponse: To isolate the atmospheric dynamics and the ocean dynamics, the model is set in\na way that as soon as it reaches the coast, it calves. Our model outputs are therefore\nconsidering a simple calving scheme. However, we interpret the timings of advance and\nretreat of the ice in agreement with Hagemann et al. (2024), but also with the geologically\nconstrained margin of the PIS. Moreover, our model simulations are forced by climate model\ninput that due to its coarse resolution and the poorly constrained topography are unable to\ncapture the climate variability with enough proficiency to enable us to reproduce in a fully\nrealistically way the geochronologically reconstructed margins and the sediments forcing\nrecord.\nMinor Comments:\nLine 87-88: References Figure 2e here.\nDone\nLine 99-101: Does the Be10/Be9 contradict the results of PIS advancement at the start of\nMIS4? The Be10/Be9 is indicative of terrigenous sedimentation beginning at ~80kya, with\noscillations preceding it. The record across MIS5 is of a similar magnitude across MIS3,\ndespite the differing ice extents in the Uk37 reconstruction.\nResponse: Sproson et al. (2024) have used Be10/Be9 as an indicator of the glacial variation\nof the PIS, identifying two main periods of glacial advance. The periods broadly align with\nour timings of the glacial advance of both integrated summer insolation and MR19-09 PC03.\nThe age model uncertainties increase particularly during the proposed glacial advance\naround the MIS4. This has been included in the new manuscript.\nLine 104-105: Which records of Hagemann et al., (2024) suggest that Uk37 is the preferred\nrecord? The Uk37 SST record is used as the forcing, so that is circular. Is it the MAR record?\nIf so, I recommend being specific about which records provide support.\nWe referred to the MAR record. This was corrected in the manuscript.\nLine 115: Is the PATICE reconstruction meant to be shown in Figure 2b?\nPATICE reconstruction has been shown for comparison\nReferences\nCastillo-Llarena, A., Retamal-Ramírez, F., Bernales, J., Jacques-Coper, M., Prange, M., &\nRogozhina, I. (2024). Climate and ice sheet dynamics in Patagonia throughout marine\nisotope stages 2 and 3. Climate of the Past, 20(7), 1559-1577.\nCharbit, S., Ritz, C., Philippon, G., Peyaud, V., & Kageyama, M. (2007). Numerical\nreconstructions of the Northern Hemisphere ice sheets through the last glacial-interglacial\ncycle. Climate of the Past, 3(1), 15-37.\nGarcía, J. L., Lüthgens, C., Vega, R. M., Rodés, Á., Hein, A. S., & Binnie, S. A. (2021). A\ncomposite 10 Be, IR-50 and 14 C chronology of the pre-Last Glacial Maximum (LGM) full ice\nextent of the western Patagonian Ice Sheet on the Isla de Chiloé, south Chile (42° S). E&G\nQuaternary Science Journal, 70(1), 105-128.\nGlasser, N. F., Jansson, K. N., Goodfellow, B. W., De Angelis, H., Rodnight, H., & Rood, D.\nH. (2011). Cosmogenic nuclide exposure ages for moraines in the Lago San Martin Valley,\nArgentina. Quaternary Research, 75(3), 636-646.\nHagemann, J. R., Lamy, F., Arz, H. W., Lembke-Jene, L., Auderset, A., Harada, N., ... &\nTiedemann, R. (2024). A marine record of Patagonian ice sheet changes over the past\n140,000 years. Proceedings of the National Academy of Sciences, 121(12), e2302983121.\nHall, B. L., Lowell, T. V., & Brickle, P. (2020). Multiple glacial maxima of similar extent at∼\n20–45 ka on Mt. Usborne, East Falkland, South Atlantic region. Quaternary Science\nReviews, 250, 106677.\nHein, A. S., Hulton, N. R., Dunai, T. J., Sugden, D. E., Kaplan, M. R., & Xu, S. (2010). The\nchronology of the Last Glacial Maximum and deglacial events in central Argentine\nPatagonia. Quaternary Science Reviews, 29(9-10), 1212-1227.\nKaplan, M. R., Fogwill, C. J., Sugden, D. E., Hulton, N. R. J., Kubik, P. W., & Freeman, S. P.\nH. T. (2008). Southern Patagonian glacial chronology for the Last Glacial period and\nimplications for Southern Ocean climate. Quaternary Science Reviews, 27(3-4), 284-294.\nLešić, N. M., Streuff, K. T., Bohrmann, G., & Kuhn, G. (2022). Glacimarine sediments from\nouter Drygalski Trough, sub-Antarctic South Georgia–evidence for extensive glaciation\nduring the Last Glacial Maximum. Quaternary Science Reviews, 292, 107657.\nMas e Braga, M., Bernales, J., Prange, M., Stroeven, A. P., & Rogozhina, I. (2021).\nSensitivity of the Antarctic ice sheets to the warming of marine isotope substage 11c. The\nCryosphere, 15(1), 459-478.\nMendelová, M., Hein, A. S., Rodés, Á., & Xu, S. (2020). Extensive mountain glaciation in\ncentral Patagonia during Marine Isotope Stage 5. Quaternary Science Reviews, 227,\n105996.\nP. I. Moreno et al., Radiocarbon chronology of the last glacial maximum and its termination in\nnorthwestern Patagonia. Quaternary Science Reviews 122, 233-249 (2015).\nNiu, L. U., Lohmann, G., Hinck, S., Gowan, E. J., & Krebs-Kanzow, U. T. A. (2019). The\nsensitivity of Northern Hemisphere ice sheets to atmospheric forcing during the last glacial\ncycle using PMIP3 models. Journal of Glaciology, 65(252), 645-661.\nNiu, L., Lohmann, G., Gierz, P., Gowan, E. J., & Knorr, G. (2021). Coupled climate-ice sheet\nmodelling of MIS-13 reveals a sensitive Cordilleran Ice Sheet. Global and Planetary Change,\n200, 103474.\nPeltier, C., Kaplan, M. R., Birkel, S. D., Soteres, R. L., Sagredo, E. A., Aravena, J. C., ... &\nSchaefer, J. M. (2021). The large MIS 4 and long MIS 2 glacier maxima on the southern tip\nof South America. Quaternary Science Reviews, 262, 106858.\nRudolph, E. M., Hedding, D. W., Fabel, D., Hodgson, D. A., Gheorghiu, D. M., Shanks, R., &\nNel, W. (2020). Early glacial maximum and deglaciation at sub-Antarctic Marion Island from\ncosmogenic 36Cl exposure dating. Quaternary Science Reviews, 231, 106208.\nSproson, A. D., Yokoyama, Y., Miyairi, Y., Aze, T., Clementi, V. J., Riechelson, H., ... &\nChildress, L. B. (2024). Near-synchronous Northern Hemisphere and Patagonian Ice Sheet\nvariation over the last glacial cycle. Nature Geoscience, 17(5), 450-457.\nSugden, D. E., Hulton, N. R., & Purves, R. S. (2002). Modelling the inception of the\nPatagonian icesheet. Quaternary International, 95, 55-64.\nWolff, I. W., Glasser, N. F., Harrison, S., Wood, J. L., & Hubbard, A. (2023). A steady-state\nmodel reconstruction of the patagonian ice sheet during the last glacial maximum.\nQuaternary Science Advances, 12, 100103.\nYan, Q., Wei, T., & Zhang, Z. (2022). Modeling the climate sensitivity of Patagonian glaciers\nand their responses to climatic change during the global last glacial maximum. Quaternary\nScience Reviews, 288, 107582."}, {"round": 2, "reviewer_comments": {"reviewer_1": "I've gone through my initial review, the author's response, and the kupdated manuscript. Everything looks good to me. I think the authros did a great job addressing everything. Nothing further from me.\n\n(Remarks on code availability)\nI looked at the code, but did not \"review\" it. I'm not a modeller or work with code, so I cant rally comment on it. Seemed like the other reviewers didnt find any faults with the code.", "reviewer_2": "Thank you to the authors for taking on my concerns and comments. I have taken a look through all the edits which have been provided, along with the additions and believe all my concerns have been addresses and I am happy to accept it as it is.\n\n(Remarks on code availability)\nGood looks good and is accessible for the build and code needed for the model.", "reviewer_3": "The manuscript has been revised based on comments from 3 reviewers. My first review primarily focused on how the uncertainty of the two Uk37 temperature reconstructions impacts the glacial index forcing, and subsequently the two PIS reconstructions. I have a few remaining minor comments surrounding this, but otherwise the authors have adequately addressed all of my other points. Thanks to the authors for the detailed responses and discussion. Well done on an excellent manuscript!\n\nMain point:\nI agree that the difference between the two PI-LGM anomalies is increased by the presence of millennial-scale variability in the Uk37 record, which then determines whether MIS5 is more LGM-like or more PI-like. This is the key aspect to note, and I’m glad that the authors have included this. However, I disagree that the glacial index is completely independent from the uncertainty in the SST reconstructions. For example: If the 1.5oC difference in PI-LGM anomaly is enough to drive more LGM-like conditions during MIS 5, and the uncertainty within one of the SST reconstructions is on the order of 1.5oC, then it’s challenging to distinguish the scenarios based at the face-value of the SST reconstructions. Additionally, the PI interval (average of 0-3ka) is determined by a very small amount of data from Hagemann et al. 2024, which differs by ~1.5oC between the two reconstruction methods. If the PI data is systematically biased by either variability in the calibration, during measurement, or from other environmental factors, this could reproduce very different ice sheet histories.\n\nThis all being said, I absolutely agree with the authors that the Uk37 reconstruction is supported based on MAR at the site, other landform markers, and Be10 at a nearby site. I think that adequately addressing this is outside the scope of this paper, and that this difference does not impact the author’s main points surrounding insolation and millennial-scale variability. My only recommendation is that lines 152-153:\n“The glacial index method and therefore the modelled ice sheet behavior are independent of the uncertainty in the comparison of the temperature records and offsets of both proxies”\nBe either removed, or revised to something along the lines of:\n“The glacial index method and therefore the modelled ice sheet behavior are influenced by the presence of millennial-scale climate variability present in the MR16-09 PC03-based climate reconstructions.”\nSuch that it notes the influence of the data and supports the importance of millennial-scale variability on the reconstruction.\n\nMinor points:\n1) Figure 2: I would plot EDML d18O on normal y-direction rather than inverted, such that it qualitatively aligns with SST.\n2) Line 96: Should “distinction” be “distinct”?\n3) Line 170: Should “Leeside” be “leeside”?"}, "author_response": "Final Response to Reviewer #1\nI've gone through my initial review, the author's response, and the kupdated manuscript.\nEverything looks good to me. I think the authros did a great job addressing everything.\nNothing further from me.\nReviewer #1 (Remarks on code availability):\nI looked at the code, but did not \"review\" it. I'm not a modeller or work with code, so I cant\nrally comment on it. Seemed like the other reviewers didnt find any faults with the code.\nResponse: We sincerely thank Reviewer #1 for the positive feedback.\nFinal Response to Reviewer #2\nThank you to the authors for taking on my concerns and comments. I have taken a look\nthrough all the edits which have been provided, along with the additions and believe all my\nconcerns have been addresses and I am happy to accept it as it is.\nReviewer #2 (Remarks on code availability):\nGood looks good and is accessible for the build and code needed for the model.\nResponse: We sincerely thank Reviewer #2 for the positive feedback. We are pleased that\nthe revisions and additions have satisfactorily addressed all concerns.\nFinal Response to Reviewer #3\nThe manuscript has been revised based on comments from 3 reviewers. My first review\nprimarily focused on how the uncertainty of the two Uk37 temperature reconstructions\nimpacts the glacial index forcing, and subsequently the two PIS reconstructions. I have a few\nremaining minor comments surrounding this, but otherwise the authors have adequately\naddressed all of my other points. Thanks to the authors for the detailed responses and\ndiscussion. Well done on an excellent manuscript!\nResponse: We sincerely thank Reviewer #3 for the feedback. We have implemented all\nchanges suggested by the reviewer.\nMain point:\nI agree that the difference between the two PI-LGM anomalies is increased by the presence\nof millennial-scale variability in the Uk37 record, which then determines whether MIS5 is\nmore LGM-like or more PI-like. This is the key aspect to note, and I’m glad that the authors\nhave included this. However, I disagree that the glacial index is completely independent from\nthe uncertainty in the SST reconstructions. For example: If the 1.5oC difference in PI-LGM\nanomaly is enough to drive more LGM-like conditions during MIS 5, and the uncertainty\nwithin one of the SST reconstructions is on the order of 1.5oC, then it’s challenging to\ndistinguish the scenarios based at the face-value of the SST reconstructions. Additionally,\nthe PI interval (average of 0-3ka) is determined by a very small amount of data from\nHagemann et al. 2024, which differs by ~1.5oC between the two reconstruction methods. If\nthe PI data is systematically biased by either variability in the calibration, during\nmeasurement, or from other environmental factors, this could reproduce very different ice\nsheet histories.\nThis all being said, I absolutely agree with the authors that the Uk37 reconstruction is\nsupported based on MAR at the site, other landform markers, and Be10 at a nearby site. I\nthink that adequately addressing this is outside the scope of this paper, and that this\ndifference does not impact the author’s main points surrounding insolation and\nmillennial-scale variability.\nResponse: We thank the reviewer for this thoughtful and constructive remark. We\nacknowledge that the limited PI data, and potential systematic biases in reconstructions, add\nto this uncertainty. As the reviewer notes, however, a full quantification of these uncertainties\nand their implications for ice-sheet histories lies beyond the scope of the present study. We\nare reassured that the reviewer concurs that our main conclusions, particularly regarding the\nrole of insolation and millennial-scale variability, remain robust and well supported by the\nindependent lines of evidence considered.\nMy only recommendation is that lines 152-153:\n“The glacial index method and therefore the modelled ice sheet behavior are independent of\nthe uncertainty in the comparison of the temperature records and offsets of both proxies”\nBe either removed, or revised to something along the lines of:\n“The glacial index method and therefore the modelled ice sheet behavior are influenced by\nthe presence of millennial-scale climate variability present in the MR16-09 PC03-based\nclimate reconstructions.”\nSuch that it notes the influence of the data and supports the importance of millennial-scale\nvariability on the reconstruction.\nResponse: Thanks for those thoughtful comments. The sentence, as suggested, was\nremoved.\nMinor points:\n1) Figure 2: I would plot EDML d18O on normal y-direction rather than inverted, such that it\nqualitatively aligns with SST.\nResponse: Corrected. Figure has been modified accordingly.\n2) Line 96: Should “distinction” be “distinct”?\nResponse: Corrected.\n3) Line 170: Should “Leeside” be “leeside”?\nResponse: Corrected."}]

Decision letter