Castillo-Llarena and co-authors present numerical experiments to reconstruct the Patagonian ice sheet (PIS) at the last glacial maximum (LGM, 23-19 ka). They use an hybrid shallow ice / shelfy-stream model at 8 km resolution forced by PMIP outputs super-imposed to ERA5. They run 10 kyr under perpetual climate forcing to reach ice sheet equilibrium. Their major finding is that, using this methodology, the PMIP model ensemble does not reproduce the ice sheet extent inferred from geological evidence. 

Numerical simulations of the Patagonian ice sheet are scarce, this is a nice motivation for this work. The methodology itself is sound, although not novel since it has been used multiple times, notably to simulate Northern Hemisphere ice sheets at the LGM. However, I have to admit that I am not convinced that it is appropriate in this case. The authors are trying to reproduce an ice field that is much smaller than the other LGM ice sheets (at max. 1/4th of present-day Greenland ice sheet). The topographic setup is also much more complex than in Greenland (for example) since the whole Patagonian ice sheet (PIS) is sitting over the Cordilleran mountains which display an extremely rough terrain. These lead to a spatially highly heterogeneous climate (altitude dependency, rain shadow effects etc.). The authors do not apply any specific corrections to the PMIP models (only a classical vertical lapse rate correction). Worse, they use ERA5 to get a “better” pre-industrial climate but in doing so they keep the imprint of ERA 5 (pretty much linked to present-day topography) at the LGM. With an ice sheet, the topography becomes necessarily much smoother so the complexities seen in ERA5 do no longer apply. The methodology used here has be followed for other much larger objects (e.g. Greenland), and even in this case it is generally accepted that the Eastern part of the Greenland ice sheet (which has a rough topography) is generally poorly represented with such “anomaly method”. I think that one reason for the small available numerical reconstruction of this ice sheet in the literature is due to the fact that it is a complex thing to model since it requires an adequate climate downscaling.

A second comment on spatial resolution is that the topography is very complex so we expect a LGM ice flow that displays strong horizontal gradients. The authors use a 8 km spatial resolution which seems very coarse given the spatial scales studied here. To my knowledge (even though I have not run the model in its latest revision), SICOPOLIS is not very demanding in terms of computational cost. If really the computing time is an issue I would suggest to reduce the vertical discretisation since 81 points in the ice (while SIA displays a relatively smooth profile) and 41 in the bedrock (where only temperature diffusion is solved!) seem a bit too much. Although I guess the computational cost is mostly linked to horizontal discretisation (since you might need to reduce the timestep to avoid numerical instabilities).

In summary, my major comment is on spatial resolution and, above all, climate downscaling. The authors conclusion (PMIP inability to reproduce LGM PIS) was pretty much expected. This study appears after the one of Yan et al. (2022) but it is much less thorough. Yan et al. (2022) used also the same, potentially controversial, method but they also tested climate sensitivities (temperature and precipitation alone, PMIP anomalies), climate parameters (PDD factors) and ice dynamics (sliding coefficient). I agree that it is better to have different studies on a single scientific question but here it does not bring any new piece of information since the methodology is the same (here with a coarser resolution) but with less discussion on the PMIP outputs.

Major comments

1) Climate downscaling & ice sheet resolution

As explained in the introduction of this review and I am not convinced that the methodology is suited for the LGM PIS. In addition this study arrives later than Yan et al. (2022). For this reason some additional work should have been done to bring some kind of novelty, such as for example, a proper discussion of the limitation of the employed method. It is possible that no regional atmospheric model simulation of Patagonia at the LGM has been run (I am not familiar with the literature) but maybe the authors could have suggested a few possible alternatives? Also, even given the huge biases of the PMIP model, why not running the ice sheet model forced by the PMIP outputs directly with no anomaly (or maybe some domain-wise correction). Or perhaps the authors could have smoothed the ERA 5 climatology to be compatible with the PMIP model resolution before applying the anomaly. Keeping the imprint of ERA5 does not make sense to me. 

Also an ice sheet resolution of at least 2 km (preferably smaller) should have been employed given the rough topography.

2) Climate evaluation.

There is no map of temperature and precipitation under modern condition nor for the LGM. Is ERA5 good for the Patagonian region? Do you simulate the extent of present-day ice fields when using ERA5? If not, perhaps ERA5 needs also some correction. The LGM PMIP anomalies (January temperature and precip) should have been shown (in the supplement if too many maps). I liked the proxy data discussion but I found it too weak: where is located the pollen-based data? Have you selected a model grid point (outside of the LGM ice mask used) close to this site? This could have been a novelty compared to Yan et al. (2022) to include a more detailed proxy-data comparison (but I am not familiar with this data so perhaps it is qualitative than quantitative).

3) At several occasion in the text the authors mention the “chronology” of the PIS maximum. I generally agree with the authors about non-synchronicity of climate change. However I don’t think the present study brings any information about this. Instead of equilibrium simulations, it could have been possible to run transient ice model simulations using an index method (using an index representative of Southern Hemisphere temperature change). This would also have been a real improvement with respect to Yan et al. (2022).

4) Glacial erosion / sediments.

There was large glacial erosion during the last glacial cycle. So it is a simplification to use present-day topography. This is not discussed here. In addition, testing the model response with a different topography map could have been a novelty with respect to Yan et al. (2022). However, the main control here of ice sheet extent is the climate forcing so it is a much minor comment than the other first two.

Minor comments

- I do not agree with your sentence on P1L16: at the LGM the topography is much less complex than for present-day condition since the ice sheet is expected to drastically smooth the topography.

- Sentence P2L48-50 seems a bit out of context here since mostly EMICs use interactive ice sheets and the authors mostly describe GCM models through PMIP experiments (with no interactive ice sheets).

- Sec. 2.1: dragging law used?

- Sec. 2.2: how do you compute the anomalies? You should have put the LGM outputs and the PI outputs to the same topography (since the GCM sees different topographies). You mention the lapse rate on the previous section, is it also used to compute the PMIP anomalies? What has been done to precipitation do you account for the drying linked to temperature decrease?

- Sec. 2.2: show ERA 5 January temperature and annual precipitation. If there is a strong correlation with topography (which I expect) then the anomaly method is not suited (since the topography is much flatter at the LGM given the presence of the ice sheet).

- Sec. 2.2: resolution of ERA5 used?

- Sec. 2.2: should cite previous papers that use the anomaly method such as, e.g., Charbit et al. (2002)

- Sec. 3.x: I am not sure that Fig. 3 is the best illustrative figure. What we see in this figure is that the model strongly disagree, it is hard to see an emerging pattern. Some maps might help to understand the signal better. Maybe maps of surface mass balance for present-day topography or on ICE6G-C topography could have been useful.

- Sec. 3.2: the eastern extent is probably due to a combination of the strong rain shadow effect in ERA 5. Is this rain shadow effect is expected to be so strong at the LGM given the flatter topography and different zonal winds? This is something to discuss.

- P8 L148-150: should we need a model to draw this conclusion?

- Ice mask and topography used in PMIP models. I think this point is really important and it is good that it is brought forward. It should be clearly stated if the PMIP models use an ice mask in agreement with the reconstruction or not. Generally there is a strong albedo effect so the local temperature simulated by the model is strongly impacted by the presence of an ice sheet. If the GCM does not have an ice mask in its boundary conditions it will be hard to reconstruct an ice sheet there. Is this the case here? Do you have a strong correlation of the January temperature anomaly and the ice mask? It seems that the ice mask is shifted to the East in models that do have an ice mask. This is probably the reason of the too large simulated extent there?

Technical comments

- P5L86 Lambeck et al. (2014) have a drop of about 140 m of sea level.

- Fig. 1: it is hard to see the topography with this grey colour gradient (while we see fine the ocean, but we don’t need it). Generally the figures are a bit hard to read.

- Fig. 2: the green line is barely visible when there is a simulated ice sheet.

Reference

Charbit, S., Ritz, C., and Ramstein, G.: Simulations of Northern Hemisphere ice-sheet retreat: sensitivity to physical mechanisms involved during the Last Deglaciation, Quater. Sci. Rev., 21, 243–265, 2002.

This work presents a set of steady-state numerical simulations of the Patagonian Ice Sheet (PIS) during the Last Glacial Maximum (LGM). The applied model is the well-known polythermal ice-sheet model SICOPOLIS, which has been forced by a set of climatic products from the PMIP3 and PMIP4 experiments. The results show that none of the considered forcings are able to build proper ice extent in the northern part of the PIS, in contrast to geological evidence, while there is a tendency to overestimate ice growth in the south. The authors attribute this mismatch mainly to the low resolution of PMIP climate models, which hampers a correct representation of the atmospheric processes and dynamics over the complex topography in Patagonia. 

Paleoclimate modelling studies applied to the PIS are certainly scarce compared to the abundance of comparable studies for polar regions. However they have gained recent attention  with more papers coming out (e.g. see Yan et al., 2022 but also Cuzzone et al., 2023 in review in TC). The motivation behind this new work is well supported by the recent advance in ice-flow and climate modelling as well as by new information on ice sheet area based on geological evidence (Davies et al., 2020), bringing further understanding of the glacial state /deglaciation of the PIS. However,  this work comes after Yan et al., 2022 (LGM steady-state PIS simulations) and mostly at the same time as Cuzzone et al., 2023 (deglaciation simulations of the northern sector of PIS) and it seems it struggles to offer a convincing motivation to be considered as a substantial contribution to the scientific knowledge. 

LGM equilibrium numerical simulations of the PIS are already tackled in Yan et al paper, with a higher spatial resolution (1 km), with a more thorough methodology, 21 PMIP model applied, with additional sensitivity tests both on temperature and precipitation as well as on the choice of some model parameters (PDD factors, enhancement factor for SIA and basal friction law exponent). In comparison to Yan’s paper, the authors here “provide potential reasons for discrepancies” between modelled and reconstructed ice growth at the LGM associated to PMIP forcings, basically based on low model spatial resolution and unresolved topography. However I think the analysis they provide in the submitted version of the manuscript might not be sufficient to consider this work as a substantial novelty to Yan et al., paper. The conclusion “all the ensemble members driven by PMIP products are not able to reproduce the reconstructed ice cover in the northern part of Patagonia” could be already seen in Yan et al., 2022. Still, with a different model setup, they can grow more ice in the north compared to this work (although still underestimating Davies et al., 2020), showing that the conclusion of this paper cannot be drawn that straightforwardly. I still believe that PMIP products provide a bad performance mainly due to the oversimplified terrain as a forcing, which instead is extremely rough in Patagonia, and thus not being able to capture the complex surface atmospheric dynamics shaped by the Andes, but the authors should clearly prove this, tackling other possible sources of mismatch: for example how are their results influenced by the choice of the present climatology and its spatial resolution (ERA5), by the spatial ice-sheet model resolution (8 km), and - less - by the choice of key ice-sheet model parameters (PDD factors, basal friction coefficients)? Following this direction, I strongly suggest the authors to include at least sensitivity tests on the applied present reference climatologies (e.g. ERA5 and CR2MET - with 5 km res), and using a higher spatial resolution for the ice sheet model domain (at least 4 km). By increasing the resolution of both ice sheet model and climate products, despite the PMIP having still a very low resolution, the authors should be able to pinpoint the cause of the bad performance in the north more clearly. Otherwise one could say that the inability of building ice there might be also due to the low spatial resolution of the ice sheet domain, or to the present climatology applied to the anomaly method. Discarding other possible sources of mismatch could help to demonstrate that the unresolved complex topography used to force PMIP climate models is one of the main causes of model-data discrepancy in terms of ice sheet extent.

LGM state and transient deglaciation simulations are tackled by Cuzzone et al. 2023, although only for the northern part of the PIS. They used the climatology from the Trace21-ka experiment, and it performed quite well, although the climate spatial resolution is lower (3.75°), but perhaps based on a better resolved topography (ICE 5G if I am not wrong) (and different experimental design). I would be curious to see how the Trace21-ka product for the LGM performs for this work. If results improved with Trace21-ka (better performance despite the low resolution), the conclusion stating “We largely attribute these discrepancies between the model-based ice geometries and geological evidence to the low resolution of paleoclimate models” would lose significance and causes of the mismatch should be searched elsewhere. 

Another suggestion to make this work more appealing could be developing more aspects of the manuscript that are discussed but not analysed. For example, the authors discuss extensively about the limitations in approximating the LGM in Patagonia from the global LGM (~21 ka) despite  evidence of asynchronous glacial maximum. To overcome this drawback they could think about a a methodology to correct the PMIP climatologies to take into account this spatial and temporal heterogeneity (maybe performing transient simulations between ~35 kyr and ~18 kyr using a climatic index?). This suggestion goes in the direction of the other reviewer’s comment. Also, I think there are now available other PMIP4 LGM climate products from other institutions (IPSL, NCAR). These could be added to the analysis.

Regarding the work presentation quality, I think that generally the manuscript is well written and clear in most of its parts, besides some sentences that need clarifications and figures that might be added and others modified. 

Besides the major concerns I expressed above, here I note down some general and specific comments to the manuscript. 

General comments:

Model set-up

I strongly suggest the authors to include more information about the model setup, such as the basal friction law, the applied oceanic forcing (if any?), how ice shelves and grounding line migration are treated. Results of modelled ice dynamics should be better described (there are almost no comments on this in the manuscript): e.g. are ice streams modelled? Where? How are the velocities distributed? Basal stress plays an important role on the ice-sheet capability to advance, therefore having an impact on the area that is glaciated. Some comment on this is required, maybe presenting sensitivity tests on the basal friction coefficient and/or law exponent. Is the basal hydrology taken into account? If yes, how? I also suggest to add a figure that shows the spatial distribution of the horizontal velocity.  

Ice-sheet model and climate model spatial resolution

As stated before, I think one cannot really prove that the main cause of the LGM data-model mismatch is due to the poor resolution of paleoclimate models as long as the spatial resolution of the ice sheet model itself is too coarse (8 km). This applies also to the applied reference climatology (ERA5). I therefore suggest to increase the spatial resolution of SICOPOLIS to at least 4 km (and possibly to try out a modern climate with a higher resolution such as CR2MET) to better capture the interplay between topography and temperature and precipitation.

Impact of climate forcing on area change

To me the manuscript lacks a clear figure showing how the glaciated area changes with respect to temperature and precipitation thresholds, as inferred from climate outputs for different latitude ranges. Figure 3 provides an interesting perspective on climate specifications at various latitudes, but it is very difficult to estimate from this figure which are the minimum conditions in terms of temperature and precipitation required for ice growth over a certain region. I would like to see a clear figure (maybe a scatterplot) that shows area change versus model precipitation and temperature averaged over a certain latitude interval. This information could also be stored in Figure 3 somehow, so that the models that satisfy the required ice sheet advance are indicated with a different symbol, maybe. Also, how are these thresholds established? I would suggest to think about a more thorough mathematical description to identify such temperature and precipitation thresholds. This could simply be the relative error between the modelled and the reconstructed area being lower than an error bound (e.g. 20%), or the ratio between the number of grids cells where ice is both (or neither) modelled and reconstructed and the total number of grid cells, being higher than a certain value (e.g. the temperature and/or precipitation that ensures that 80% of cells are in agreement with data corresponds to the minimum threshold). 

Clarity of the results

As I wrote before, I think the authors overall do a good job in describing and discussing the results. Still, I see quite confusing the fact that the investigated latitude ranges change over the manuscript: sometimes is the 38-39°S and the 40-42°S, sometimes is the 38-42°S, sometimes is “below 44°S”. I think a more homogeneous description would be beneficial for the paper. 

Specific comments:

• P1 L11-12, but also P20 L409-410 and within the text in several paragraphs. As I wrote in my general comments I think it is difficult to prove that a principal source of mismatch between model and data is due to the low resolution of paleoclimatic models, as long as the low resolution of the ice-sheet domain and the present day climate also contribute to smooth the climate gradients over the complex topography. Please refer to my comments above.  
• P2 L26. "...consequently global sea level dropped to 120-134 m…”, of course the sea level drop from all vanished ice sheet doesn’t sum up to 120-134 m because of the contribution from AIS and GrIS that here is not described. Please clarify this. 
• P2 L29-30. Unclear description, please change to something like “The latter triggered a lowering of the global mean surface air temperature by 3.2°C to 6.7°C with respect to the preindustrial level…”.
• P2 L35-42. Why do you describe in detail the asynchronous occurrence of the glacial maximum when you actually consider the PIS as in a steady state during the global LGM (~21 ka)? Please refer to my comment above about possible transient simulations. 
• P2 L42. Please rephrase “zooming in on the global climate…”
• P3 L56. Which PMIP outputs? Phase?
• P3 L59-63. First, I really don’t see that PMIP4 models perform better than PMIP3 in Yan et al., 2022; please consider a reformulation. Second, if PMIP3 clearly performs worse than PMIP4 from Yan et al., 2022 paper, then why do you also investigate PIMP3 climate outputs? Third, Yan et al., 2022 look in total to 21 model products from PMIP2 to PMIP4, you should write that. 
• P3. Are there other previous model experiments/reconstructions besides what you described here? You could mention for instance the new PIS thickness and volume reconstruction from a perfect plasticity assumption from Wolff et al., 2023. This is actually an interesting paper, that could also be mentioned in the discussion to compare your results against, as it shows a PIS further extended to the north (almost to 36°S, as based on an early reconstruction from Hubbard et al., 2005) with respect to Davies et al., 2020.  
• P5 L85-96 and Table1. How are the PDD factors chosen? Are they calibrated for the Patagonian region at the present or simply taken from the literature?
• P5 L98. Are ice shelves allowed to grow in this model setup? How are the other key model parameterisations considered (basal friction, basal hydrology, enhancement factor, …)? Please refer to my general comment. 
• P6 L103. “During phase 3 and 4 of PMIP” please specify which is the time considered (~21 ka).
• P7 L122. “…relative to the PATICE reconstruction” please add “Figure 1”.  
• P7 L126. Change Fig 2 k, l to Fig 2 k,m.
• P8. As I mentioned in the general comment, sections 3.1, 3.2 and 3.3 should describe the same latitudinal intervals as in Figure 3 c-f for consistency. Also Section 3.2 could be split into 38-42°S, 42-44°S and 44°S-52°S to better describe the performance in the northernmost part. 
• P8 L151. “Performance of these models north of 44°S”, do you mean north of 44°S but south of 42°S?
• P8 L153. “The resulting ice sheet temperature” should be “the resulting ice sheet extent”, I think.
• P8 L153-155. These cited thresholds should be argued with a more precise definition, e.g. considering the relative error of the glaciated area and with a figure showing how the area changes when these minimum conditions are met (see my comments above). 
• P9 L173-176. What do you mean by “northernmost margin”? 38-40°S? What do you mean with “In this part of Patagonia”? Why MPI-ESM1-2-LR “THUS stand out as the only PMIP products providing…”? You should describe better the climatic condition of MPI-ESM1-2-LR as you did with INM-CM4-8.
• P9 L181-184. Again, I struggle to see this from Fig. 3. Also, this is a repetition of P8 L153-155. Finally, is this threshold computed for the region 42-44°S or for north of 42°S or for what? Please clarify. 
• P11 L209. “38-42°S” should be consistent to the sectors described in the results and in Figure 3 c-f.
• P9 L220. Please change 4°C to -4 °C.
• P9 L221. Add a closing bracket to Fig. 5.
• P13 L240-241. “The cooling of ~12 °C observed in INM…during summer months…” I don’t really see this clearly. Maybe Figure 5 should be for 40-42°S or think about producing another figure to be consistent to latitude sections of Fig 3 c-f.
• P14 L244. “Infers a value of around -8°C”, which value? Annual temperature at 40-42°S? I don’t see this from Figure 3d (to me the temperature anomaly is more around ~11 °C). Maybe the 8°C anomaly is calculated at 38-42°S?
• P14 L245-246. Please refer to the figure where we can see this.
• P14 L254-255. “To achieve a good fit with geological evidence, the PDD factors in the SMB were reduced to promote ice sheet growth”. To my understanding this is not true. As I get from Yan et al., 2022 paper, the PDD factors were tuned to reproduce modern glacier geometries. Please rephrase.
• P15 L257. “This evaluation may be biased due to the choice of model parameters”. I would be careful here: they did indeed a sensitivity study to show how the results are affected by the choice of PDD factors, so you cannot really say that the results are biased.
• P15 L259-263. I think this is really difficult to say since their PDD factor for ice (4 mm/d/°C) is very close to yours (3 mm/d/°C). True, the PDD factor for snow is lower than yours to reduce ablation, but their sensitivity test shows that it doesn’t impact significantly the total modelled ice area.
• P15 L273: do you mean “topography” instead of “forcings”? Also could you clarify for which models/PMIP phase the ICE-6G and GLAC 1D reconstructions were used?
• P15 L285. I don’t understand “has undergone visible substantial modifications”. Please rephrase. 
• P17 L309. Please add “Figure 1”.
• P18 L338. Please change “between 700 and 800 (Fig. 7a) m” to “between 700 and 800 m (Fig. 7a)”.
• P18-19 Section 4.3. I see this section interesting to be discussed but it is also somehow a reasoning for the sake of it, as it is put at the end of a whole work based on the assumption of a steady-state condition during the LGM. Perhaps this discussion would gain more interest if you consider investigating the local LGM depending on the latitude through a series of transient LGM runs (see general comments).

Figures

Figure 1: Please change the colour for the present day ice fields, LGM reconstruction and coastal lines to a list of colours that clearly differ from the topography colour palette (e.g. orange, red, magenta…).

Figure 2: Could you show where are the ice shelves located (if there are any?). Also, instead of showing the velocity streamlines I would add a figure showing the 2D velocity fields. Where are the ice streams? 

Figure 3: as pointed out already above, I see figure 3 c-f quite confusing: 1) I would add the values of the x axis to all subplots; 2) I cannot clearly see which are the models from PMIP4 and which are from PMIP3; 3) where is the subplot for latitudes 42-46°S? 4) which are the models that allow ice growth as expected from PATICE? (I would mark them somehow in the plot) 5) I am missing a scatterplot of area change versus temperature and precipitation anomaly (see my general comment) to clearly see which are the climatic thresholds that satisfy ice growth in the north. 

Figure 7: as in figure 3, why latitudes 42-46°S are missing? 

Also, could you add one figure showing the spatial variability of temperature, and another one for  precipitation, for different PMIP climate models (panels like Fig 2). This should add more information to Figure 3, where, for instance, west-east precipitation gradients cannot be seen.

References:

Wolff et al., 2023: https://doi.org/10.1016/j.qsa.2023.100103

Hubbard et al., 2005: A modelling reconstruction of the last glacial maximum ice sheet and its deglaciation in the vicinity of the Northern Patagonian Icefield, South America. Geografiska Annaler: Series A, Physical Geography, 87(2), pp.375-391.

Cuzzone et al., 2023, https://doi.org/10.5194/tc-2023-68.