Patterson et al. present coupled ice sheet – climate simulations of the last two glacial maxima (last and penultimate glacial maxima, LGM and PGM). They run an ensemble of simulations varying model parameters for the LGM and the PGM. The methodology follows what their group has done in several previous papers for the LGM. From the large ensemble, Patterson et al. keep only the members that match the estimates of LGM ice sheet volume and extent. They then discuss the simulated ice sheets at the PGM given the LGM constraints.

The fact of using coupled ice sheet – climate model with a fast GCM is certainly an improvement with respect to most of the existing literature. However I am not sure that the methodology chosen is the best approach to study PGM and LGM differences. It is unclear to me so far what the main message of the paper should be since I do not currently really see what is the major new finding. I provide my comments below.

Major comments

1) It seems that the major conclusion of this study is that there is a very strong sensitivity of the simulated ice sheets to the initial ice sheet configuration. This is also what stands out from the abstract. This is an unsurprising result given the fact that the feedbacks linked to albedo (surface mass balance and temperature) are one of the most influential for climate and ice sheet dynamics. Although unsurprising I agree that it is worth studying and documenting. However this point first appears in Sec. 3.4, so relatively late in the manuscript, and as a sensitivity analysis. I would have expected more analysis and discussion on these simulations. For example: spatial pattern of ice gain/loss under PGM/LGM conditions, using the same ice sheet boundary conditions? Where these spatial pattern differences come from?… I have the feeling that it would have been much more informative to reshape the manuscript in order to present the experiments with identical ice sheet boundary conditions as the main results. In this way the respective impact of climate and ice sheet changes, the most interesting thing of this paper, could have been presented more clearly and thoroughly.

2) There is one thing that is not very clear to me in the experimental setup. The climate model uses PMIP4 ice sheet boundary conditions representative of LGM and PGM, as explained from L138. For Glimmer, the initial ice sheets (North America and Greenland) are the one at 18.2 ka BP from a previous experiments. Since all the simulations are bi-directionnally coupled as shown in Fig. 1 it means that the climate model in fact does not use the PMIP4 boundary conditions for North America and Greenland but use the 18.2 ka BP ice sheets. Is this correct? If yes, where does the difference in albedo discussed in Sec. 3.3 come from?

2) Sec. 3.1 is mostly methodological, which is fine, but I don’t see what can we learn from it. The best ensemble members for the LGM produce also realistic PGM ice sheets, but the fact the PGM ice sheets are smaller are linked to the chosen ice sheet boundary conditions as shown later…

In Sec. 3.2 we have not enough ensemble members to draw any conclusion I think. The only parameters that show some clear impacts are Rho and AV_GR. For the other parameters, there are good ensemble members that span the whole range. So here again I am not sure what conclusions can be drawn for the reader.

Sec 3.3 does not explain anything. Basically it is said that the LGM as a higher albedo in the saddle region which explains the more positive SMB. But why is that (see my point 2)? If this is a result of ice sheet boundary conditions (as discussed in Sec. 3.4) I have the feeling that it should be presented later.

Given this, It  would have made sense I think to group together Sec. 3.1 and 3.2 and to present Sec. 3.3 after what is currently Sec. 3.4.

3) Sensitivity to the oceanic forcing is hardly discussed. I think it is a quite strong limitation of this study. First, there is no real justification on the fact of using HadCM3. Then, how the ice sheet evolution would be impacted by different SST?

4) There is insufficient background to be able to link the parameters values to actual SMB changes. Since the equations are not shown and very little description is given we don’t know how the different parameters play in the model.

5) In the introduction there is no review of published works using ice sheet – climate models. There are now a relatively large literature for glacial inceptions, glacial terminations or the whole cycle with various climate models of intermediate complexity (CLIMBER-2, CLIMBER-X, LOVECLIM, iLOVECLIM, BERN3D, UVic). They have generally documented the impact of initial ice sheet configurations and the importance of albedo for ice sheet evolution. As such, I think they deserve at least a dedicated paragraph to understand how this paper is participating to knowledge increment with respect to this literature.

Minor comments / questions

L13-14. It is a subjective but strong statement. The ice sheet response to different insolation/GHG pathways through the last two deglaciations might also be “crucial”...

L48-49. I don’t think this is a strong constraint given the uncertainties in term of timing of the maximal extent of the Eurasian ice sheet.

L61-70. I would suggest to remove this part as it has very little link to the general purpose of the study. In addition, the differences listed here might be linked to differing ice sheet and climate configurations at the glacial maxima but they are also most likely linked to different insolation evolutions.

L112. Unclear. 50 decades of climate years, meaning 500 years simulated per day? Seems very quick for a low-res GCM.

Fig. 1. Explain better what is represented. Horizontal line in top left map? Solid line in bottom left graph?

L122-123. Strictly speaking you do not follow the PMIP4 protocols since you use interactive ice sheets that overwrite the ice sheets (as shown in Fig. 1). Also Menviel et al. (2019) present a protocol for deglaciation with prescribed ice sheet.

Table 1. You should add a column with the reference for the ice sheets.  

L138-157. Please show the ice sheet boundary conditions used in the climate model (including over Eurasia) for the PGM and LGM.

L145. “constant” but with a seasonal cycle right? Daily forcing?

L146. Reference for these simulations? Why these forcings and not FAMOUS computed SST and sea ice for consistency?

L148. Show summer SST as well since it seems important.

L172-173. Show difference between HadCM3 SST and reconstructions?

Table 2. Rho seems to be Fsnow in Gandy et al. (2023). Be consistent (at least in the paper).

Table 2. Description is generally too vague. For most parameters we cannot guess in which direction the parameters can influence the simulated climate, SMB or ice sheets. 

Table 2. Please include the range tested for each parameters.

Fig. 5. Draw the 1:1 line in b and c.

L258. 4 parameters are listed here, including basal sliding. While in L458 the flow factor is mentioned and not basal sliding. Why?

L258. From the plot I clearly see a tendency for Rho and AV_GR but for the two others it is much harder. For example for Daice we see good ensemble members on both side of the tested range. For basal sliding there might be a tendency but given the fact that we don’t have a lot of ensemble members here I do not think that we draw any strong statement.

L267. There is something unclear, I guess in the representation. On Fig. D3 we see ice volume difference of -1 to 1 e7 km3 which seems not minor with respect to the ice volume of about 3 e7 km3.

L267-269. I don’t see this result in the plot. Please clarify this.

Fig. 7. The good ensemble members are always on the lower hand of the reconstructions. What about the modern bias of FAMOUS-ice?

L295. “passing all reductions”, what does that mean? Aging for instance? Please clarify.

L296-297. There is something I don’t understand in the set-up. L152-155 it is said that the same initial ice sheet is used for the LGM and PGM in GLIMMER. Since you use a coupling as in Fig. 1 the climate model also sees the same initial ice sheet in the saddle region. Please clarify this. 

Table 4. Add a column with V, Vc, Vi, Vci.

Table 4. FAMOUS initial ice sheet is not GLIMMER initial ice sheet? You are talking about the ice sheets outside the GLIMMER region (Eurasia)? Unclear.

L340-348. Remove this part. I don’t understand why there is this discussion here while you don’t account for vegetation changes.

L352. This should be shown!

L359-375. I enjoy this section but it should be in a separate discussion section.

L369-370. Why comparing the insolation peak of 172 and 148 ka BP to the ones of MIS4? In terms of relative timing they should be compared with the ones of MIS3 (55 and 30 ka BP).

Fig. 11. You should use the same x-axis scale. Here there is a distortion (longer period preceding the LGM than PGM) that makes the comparison difficult to do. You could group the two cycles in one graph only.

Fig. 11. In terms of insolation in the Northern Hemisphere 21 ka BP is more comparable to 137 ka BP than 140 ka BP. You should perhaps comment on this as your results could have been slightly different if using the 137 kaBP orbital and GHG configuration.

L389. Why vegetation-albedo feedback is mentioned here since it is not tackled here?

L390. I think it is not necessarily true. It is just that the initial ice sheet configuration is more important.

Fig. A1. Relatively minor impact but with a large trend.

Fig. B1. Show difference LGM-PGM as well (and summer SST).

Fig. B1. The colour scale is not appropriate (SST of -20 degreeC are relatively rare).

L447. Average SMB over the ice sheet?

L454. Since you start from smaller ice sheets it was expected that a positive SMB was required…

L457-458. Be consistent with the parameters names in Tab. 2.

L486. Not observations.

Fig. D3. Poor quality figure.

Technical corrections

L102. Typo, “this allows to model”

L150. “The HadCM3 LGM SST”

L161. Add reference of SST here.

L170. Appendix C is mentioned before B.

L176,L177,L178. Set-up x 3

L230 Fig. 5 has not yet been mentioned. 

L263. Typo, two dots.

L327. Keep one notation: 10**6 but not 10**7

L367. Why reference to Bonelli et al. (2009) here?

L420. Opening parenthesis missing.

L420. Define nu.

In this article, Patterson et al, perform coupled ice sheet and climate simulations. They run a wide ensemble of simulations to assess the parametric uncertainty. The subsequent analyses allows then to gain some conclusions about the effects of the different orbital configurations and initial states on the final ice sheet configurations.

Coupling a GCM to an ice sheet model is in itself of great value and informative for the community. The manuscript is very well written. The introduction adequately deals with the existing knowledge of the subject. The analysis of the results is very exhaustive and clear. And the conclusions appear generally justified with respect to what is shown in the rest of the article. Therefore, I find this work is well suited for Climate of the Past, and I recommend publication subjected to some clarifications of the experimental set up and their potential implications for the main conclusions of the study.

More specifically, I found the strategy concerning initialization a bit strange and not clearly described. Thus, my main concern is about the experimental set up and is the following:

The paragraph starting at lines 123 reads: “Our FAMOUS-ice simulations are set up following the Paleoclimate Modelling Intercomparison Project Phase 4 (PMIP4) protocols for the LGM (Kageyama et al., 2017) and PGM (Menviel et al., 2019).”

Around line 139:"In the climate model, the global orography (including ice sheets) and land-sea mask for the LGM are calculated from the GLAC1D 21 ka BP reconstruction (Tarasov et al., 2012) which is one of the two options in the PMIP4 protocol (Kageyama et al., 2017). For the PGM simulations we used the 140 ka BP combined ice sheet reconstruction (Tarasov et al., 2012; Abe-Ouchi et al., 2013; Briggs et al., 2014)

And line 153 reads: “In the ice sheet model, we use the same ice sheet domain and initial condition for the LGM and PGM, [...] and the initial ice sheet extent, thickness and bedrock elevation is from a previous Last Deglaciation ensemble of the NAIS, at 18.2 ka BP”

So the reader can easily wonder why using different initial conditions for the ice–sheet and the climate models. It is not clear whether this is the best way to address the influence of the initial conditions on the final ice sheet configurations, as stated in the abstract and conclusions.

Reciprocally, concluding that the climate boundary conditions, if considered in isolation, imply a larger PGM might be dependent on the way the ice sheet initial conditions are managed under the current experimental set up. 

In other words, if the ice sheet model was initialized with an ice sheet configuration close to the PGM reconstruction (which, as far as I understood, has been used by the climate model as a boundary condition) it is conceivable that the climate does not react in the same manner than using a 18.2 kyr reconstruction, so that at the end, both the climate and the final ice sheet configurations widely differ with respect to what has been concluded here.

As a modeler, I am aware that there is not a perfect strategy for initializing the ice sheet model when the focus is on two single time snapshots. It is understandable then that using a previous deglaciation run at 18.2 kyrs has the advantage that the temperature profiles and thus viscosity have at least some internal consistency. 

However, someone could also wonder why not initializing with the ice sheet configurations that have been used as boundary conditions for the climate model (particularly so if SSTs and sea ice are fixed). You could then let the ice sheet model run to achieve internal equilibrium with the initial climate for several thousand years and subsequently “liberate” the coupled system and see where it goes.

If you have done something in these lines, I recommend incorporating it into the manuscript. If not, and you consider this suggestion unfeasible or out of the scope, please state why (there might be some subtle technical arguments I am not considering). I would still encourage the authors to include a discussion on how the choices of the initialization of the experimental set up could alter the main findings of the current study.   

Minor/technical comments:

Why using Tarasov’s reconstruction for the LGM and the combined one for the PGM?

Lines 14 and 15 of the abstract read: “Therefore, a better understanding of how and why these two glacial maxima differed is crucial for developing the full picture on why the Last Interglacial sea level was up to 9 meters higher than today, and thus may help constrain future sea level rise.”

This makes sense but is not addressed at all in the rest of the manuscript. Therefore, I suggest removing it or elaborate something in the discussion on the potential implications of your findings on this matter.