Reviewer 1: 

Response to reviewers’ comments on manuscript “Surface mass balance and climate of the Last Glacial Maximum Northern hemisphere ice sheets: simulations with CESMv2.1.

The authors would like to thank the reviewer for their positive comments and useful suggestions. Each comment is shown in italic and our response follows. Any changes to the manuscript are given in blue.

Main comments 

1: The authors use a plethora of figures, which negatively impacts the readability of the paper. At the same time, I acknowledge that this is not easy to resolve, but I kindly ask the authors to critically assess whether the volume of text and figures could be reduced and or the structure changes to make it more accessible.

The reviewer suggested that we revise the layout of the manuscript and quality and consistency of the figures. We have reduced the complexity of all figures and updated them to be consistent in style. For example, figure 8-10 are now focused on the Northern Hemisphere ice sheet region. 

2: My main concern is that the paper is a bit of a hodgepodge of topics. The many ice sheets and variables, in combination with two different time slices, make for a huge number of permutations of results to describe. Moreover, the results are mainly described in a qualitative fashion, which reduces the impact of the study

We do not agree that the paper is a hodgepodge of topics, however we have revised the structure of the paper to make reading clearer. SMB is not a quantity for the LGM (21ka or 26ka) that can be quantified as there are no direct observations with which to compare. There are a limited number of studies which have produced SMB results for 21ka (Kapsch et al., 2021) but these results are not open access. If we removed the sections of the paper describing the global climate (SAT, SST, AMOC) this would reduce the paper to be only focused on SMB, and we note that the title is SMB and climate.

3:  24: For this run, did the ice sheets in the simulation feed freshwater into the oceans? If so, how much and how was it spatially distributed?

For this setup of the CESMv2.1 model we do not run with interactive ice sheets - so the ice sheets are not evolving and do not contribute to the freshwater budget in the ocean. However, the land model (CLM5) does contribute to the freshwater budget, via the evolution of the snowpack (i.e melting or refreezing) for example. These surface and subsurface fluxes are routed to the ocean model (POP2) via the river runoff model (MOZART). Once these fluxes ‘reach’ the ocean they are spatially distributed using the estuary balance model recently developed by Sun et al., 2019, which is a physically based parameterisation of estuarine mixing processes. 

4: The comparison with Kapsch et al. (2021) would be more valuable if (spatially integrated) quantities were quantitatively compared.

We agree this would be great, however the results in the Kapsch et al., 2021 paper are not open access nor does the study provide spatially integrated quantities. For example, the detailed analysis we have produced from our study in Figure 6 of the SMB across each of the Northern hemisphere ice sheets is not replicated in the Kapsch et al., 2021 study. We are not sure what quantities which are in Kapsch study we would compare to?.We believe it is beyond the scope of our study to access other authors' work and process.  We have revised our SMB results (Fig.8 and Fig.10) to be downscaled onto the CISM (ice sheet model grid) which allows a direct comparison to the Kapsch figure. 

4: " a colder atmosphere can hold less water " This is wrong terminology. Even in the absence of other gases, the saturated water vapor pressure would be the same.

We have removed this sentence

5: l. 293: " Two ice sheets have a extremely negative SMB: NAISC and BIIS". This may be true for BIIS but not for NAIS, which is only weakly negative.

We agree that the NAISC is only weakly negative (negative ( -200 mm/yr-1) when calculated as an area average but across the ablation region (see Fig8a), it reaches > -10000mm/yr. We have edited this sentence to make this statement clearer.

“Two ice sheets have an extremely negative SMB across the ablation area: NAISC and BIIS”

6: l. 2: Please use 'higher' and 'lower' temperatures rather than 'warmer' and 'colder/cooler' temperatures throughout; I realize it is a rear-guard battle but isn’t that the reviewer's prerogative

We were using warmer/colder to typically refer to regions of the ice sheet that are “warmer” (i.e the interior) relative to the margins which are colder. However, we have added throughout the manuscript clarification when we have used warmer/colder: for example, warmer (higher) temperature, or colder (lower temperature) to ensure clarity. 

7: l. 5: These acronyms will not be familiar to all readers

We have checked acronyms throughout the text and expanded model names when first used. For example, in the abstract:

Here, we present two simulations of the Northern Hemisphere ice sheet climate and surface mass balance (SMB) using the Community Earth System model v2.1 (CESM2.1) using the Community Atmosphere Model v5 (CAM5).

8: l. 7: "snowfall accumulation "; The term 'snowfall' is commonly reserved for the solid

precipitation fraction and 'accumulation' for snowfall minus sublimation. So, to avoid confusion, preferably use 'snow accumulation'. 

We have changed this terminology as suggested and now throughout the manuscript we have used “snow accumulation”.

9: NAISC and SIS are not defined. NAISC is only defined in line 44

Thank you for pointing this out - we have checked all acronyms throughout. 

10a: l. 18: " along the equilibrium line " I suppose you mean "in the lower accumulation zone". 

10b: l. 297: " High refreezing rates are simulated along the equilibrium line altitude " See previous comment.

We were referring to the fact that refreezing is highest at the equilibrium line altitude, which can be seen in Fig10d. The equilibrium line altitude is a known feature, which separates the accumulation zone and ablation zone. But we have revised the term to be equilibrium line altitude 

Refreezing is largest along the equilibrium line altitude for all ice sheets…

High refreezing rates are simulated along the equilibrium line altitude, at the transition between accumulation and ablation zone along the southern margin.

11:l. 114: In order of preference: NOT surface temperature BUT 2 m temperature OR near-surface air temperature OR surface air temperature

Our terminology was not consistent, but we have revised to near-surface air temperature (SAT) throughout. 

12: 121: significance. Consider reducing the number of decimals, i.e., "6.8C" instead of "6.84C"

 We have checked our numbers to be consistent to 1d.p.

13: l. 135: " the former of which coincides with highest standard deviations from the model ensemble, up to 9◦C. " Unclear, please clarify/reformulate.

We agree this sentence was not clear and have revised as follows:

Note that these regions coincide with the highest standard deviations (up to 9°C) in the SAT from the ensemble of models performed by Osman et al., (2021).

14: l. 139: " 4◦C and 2◦C 140 cooling at GRIP and VOSTOK ice core sites, respectively " It is not intuitively clear whether you refer to the observed (ice core-based) temperature changes, or model changes at those ice core sites. If the latter, do they agree with the ice core-derived values?

We have clarified in the table and text that this is the predicted cooling at the location of the ice core sites. The sentence is revised as follows:

“The LG-26ka simulation is 1.5◦C colder than LG-21ka (global average), enhanced at higher latitudes, with a 4◦C and 2◦C cooling at the location of the GRIP and VOSTOK ice core sites”

Figure changes:

Fig.1: These frames are somewhat cluttered; consider removing isotherm contours (dashed/solid) and their labels and only retain shades. Increase font in (d) to match other panels.

This figure has been revised to make it clearer as suggested.

 Fig. 3: The look and feel of these panels are rather different from those in Fig. 1. Please make the format of the figures more uniform and critically assess them for errors or omissions.

            We have edited figure 1 to make the figures within the paper more consistent and simpler. 

Fig. 9: Please increase legend label spacing. changed as suggested:

           Revised

Fig. 5: Please include a legend of red and blue lines in windows a) and c). 

            This has been added.

Fig. A1: Using similar color scales for ice elevation and ocean depth is confusing for this Arctic application. 

            This has been revised. 

Minor changes all addressed:

1. 25: contributor -> contributors. 
2. 38: typo " resolutionclimate "  
3. 42: Suggest: "See Fig. 1 in Ivanovic etc." 
4. 42: solar insolation signal -> insolation  
5. 76 typo: "simulates realistically simulates" 
6. 92: surface air temperature -> 2 m air temperature 
7. 127: correspond ->. Corresponds

1. 127: correspond ->. Corresponds
2. 205: " The interior of the Laurentide " Please continue to use the acronyms previously introduced
3. 265: strength weakens -> weakens
4. 293: a extremely -> an extremely. strength weakens
5. 306: " rapid retreat would occur. " Rapid retreat would likely occur. 

1:363: that -> than can.

Reviewer 2: 

Response to reviewers’ comments on manuscript “Surface mass balance and climate of the Last Glacial Maximum Northern hemisphere ice sheets: simulations with CESMv2.1.

The authors would like to thank the reviewer for their positive comments about our manuscript stating that it was ‘intriguing, especially regarding the analysis of the Surface Mass Balance (SMB) and its decomposition into different components. Additionally, the comparisons of SMBs across six different ice sheets are novel. The model also successfully reproduces global climate and oceanic conditions, including the Atlantic Meridional Overturning Circulation (AMOC)”

Each comment from the reviewer is shown in italic and our response follows. Any changes to the manuscript are given in blue

General comment:

1. Additional analysis on the cloud radiative effect (long wave) and its relation to SMB

A recent study by Gandy et al. (2023) discussed the role of clouds on the SMB over the southern margin of the North American ice sheet. They showed that many of the PMIP models show small cloud cover over the southern margin, which caused an intense downward shortwave radiation at the surface in that region. In this study, the authors have shown that the relatively large cloud cover at the southern margin plays a role in reflecting more sun light. On the other hand, they showed a relatively strong downward long-wave radiation, which could be relative to the existence of clouds or moisture.

Possible additional: I would suggest the authors to conduct further analysis on the role of clouds (especially on the long wave side) and add a discussion on its effect on SMB. Effects of clouds on the SMB is highly uncertain in general, and also thought to be important in the future Greenland (Gregory et al. 2020). Additional analysis on this point would make this study more interesting.

Thank you for highlighting the Gandy et al., 2023 paper, which was published after this manuscript submission. We have used this paper in the later point also in relation to causes of the negative SMB. The earlier figures of the radiation budget were perhaps misleading as we do not have a large cloud cover over the southern margin (see new Fig. C4g), but also have a small (or low) cloud cover as in the PMIP4 models. To discuss the role of clouds in  more detail, we have added an additional figure into the supplementary (Fig. C4) which compares the JJA anomalies (LG-21ka - PI) for a range of variables to investigate cloud forcing in relation to the SMB. Additionally, we have included the results for the comparison study, LGM-Zhu. We have added extra text into section 4.2 to discuss the influence of clouds on the SMB

Compared to the PI simulation, across the ice sheets there is an increase in the cloud fraction (i.e gets cloudier) but the clouds get thinner. These two specific changes in the nature of the clouds can be related to the earlier responses in radiation fluxes (Section 3.3). Thinner clouds act to increase the incoming solar radiation at the surface (Fig.C4a). Conversely, in cloudier areas, the clouds increase in incoming longwave radiation although the clouds are thinner (Fig.C4g.h).

1. 
   
   Implication to the evolution of British Ice sheet during the LGM

Gandy et al. (2018) showed that a warming of climate after 26ka is required to initiate a deglaciation of the British Ice sheet in their ice sheet model simulations. However, the sure about the cause of the warming remained unclear. This study nicely showed that the differences in orbital parameters between 21ka and 26ka can cause a warming and a decrease in SMB at 21ka over the British ice sheet. Adding a discussion on this point citing Gandy et al. (2018) would be beneficial and would also increase the impact of this paper.

Thank you for pointing out this study that we had not included in our paper. We have added a short note to highlight the results from the study as follows:

Gandy et al., 2018 concluded that a warming of the climate after 26 ka, and resultant reduction in SMB was in fact required to initiate the retreat of the BIIS at 21 ka. Therefore, the 1.5°C warming between LG-26ka and LG-21ka due to the change in orbital parameters may be one factor that led to the retreat of the BIIS, due to the increase in melt rate (figure7e).

3: Discussion on the cause of the large negative SMB at the southern margin of North American ice sheet

One of the advantages of this study is the use of a higher horizontal resolution in the AGCM compared to previous studies. First, I thought that this would improve the representation of the so-called stationary wave effect, which cools down the southern margin and increases the SMB. However, even with the higher resolution used in this model, it still simulates a very large negative SMB like other models. Please add a detailed discussion on what could be the potential cause of this negative SMB, e.g. the downscaling method, representations of clouds, biases in the climate model side and so on.

The reviewer has asked us to elaborate on the causes of the large ablation area that we still produce across the NAISC.  What is interesting is that the recent study by Gandy et al., 2023, also found similar results. This this paper, that authors concluded (Section 3: Collapse of the LGM ice sheet with a standard set up”) that if they initiated their LGM (21ka) simulated with a standard large NAISC extent as we have adopted in our study, the ice sheet rapidly deglaciated due to the large ablation area that forms across the southern margin very early on in the spin-up (see Fig.3 Gandy et al., 2023). These results appear to replicate the results from our study. Gandy et al., 2023 concluded that this behaviour was due to the fact ‘FAMOUS-ice’ had been heavily tuned for present-day Greenland and the default set up was not suitable for LGM (or cold) simulations.  This could also be a similar situation for our model set up. The CESMv2.1 model has been shown to have problems when applied to the LGM (Zhu et al., 2021; 2022) and required revised ‘tuning’ to comply with LGM GMST constraints. This implies a similar approach of ‘de-tuning’ as investigated by Gandy et al., 2023 will be required for LGM ice sheet simulations using CESM.  We do simulate a less negative SMB in the 26ka climate, but we have not adjusted the ice sheet extent to represent a smaller ice sheet.  We have added a revised section to highlight this study in Discussion section:

Our simulates SMB for the NAISC and BIIS appears too negative to prevent large marginal retreat if used as a forcing for an ice sheet dynamical model. ...

A recent study by Gandy et al., 2023 also investigated the LGM NAISC in a coupled climate ice sheet model (FAMOUS-ice). This study found that when initiating their simulations from a large NAISC (as in this study), a large ablation area formed across the southern margin of the ice sheet which led to rapid ice sheet retreat (see Fig3 Gandy). This behaviour was attribute to the heavy tuning of their model to present day Greenland ice sheet. The CESMv2.1 model has also been shown to have problems when applied the LGM climate (Zhu2021 ) and required de-tuning to comply with LGM GMST constraints. Future work investigating coupled ice sheet-climate simulation for the LGM with CESMv2.1 may also require de-tuning to correctly simulate the LGM North Hemisphere ice sheets.

4:I would suggest to move subsections 3.2 and 3.3 to the beginning of the section and combine with the analysis of global temperature. In this way, the authors can start the result section from general analysis (larger scale temperature, sea ice & AMOC), followed by detailed analysis on the atmospheric circulation, which then connects well to the SMB section. Hopefully this would increase the readability of the manuscript. –

We have revised the structure of the manuscript as suggested and hope this improves the readability

Minor changes all addressed:

1:L19: It was a bit hard to interpret this sentence. In particular, what do you mean by climate variability? 

 revised as follows:

“ Our SMB results are in qualitative agreement with the variations in temperature across the north hemisphere ice sheets”.

2:L36: Perhaps, you may also cite a recent study by UKESM (Smith et al. 2021) here.

Added

3:L76: Typo, two “stimulates” changed

4: L76: It might be written in a different paper, but could you clarify whether the precipitation field is downscaled or not, please?

We have added a note to clarify as follow on line 94:

“The 4km CISM2.1 grid (Fig. A1) provides high resolution information for CLM5’s elevation class information, as well as downscaled SMB (at 4 km resolution) by horizontal bilinear and vertical interpolation from the elevation classes. (Note at present precipitation is not downscaled)”

5: L80-85: Could you add a sentence explaining the reason for these modifications? Is it just from a practical reason or is there any science behind them.

To try and reduce both the spatial extent and amount of ablation across the ice sheets, we modified parameters within the land model, as listed. We have changed the sentence to add this justification as follows:

“We make two minor modifications to the defaults settings for the elevation class parameterizations with the aim of reducing the magnitude and extent of the ablation zone”

6: L107-109: Please add a sentence on what exactly was modified/adjusted here.

To run the POP2 model, the user must define “regions” or “zones of overflow”. For each overflow zone, the user can define several parameters, including source water thickness, latitude of flow, width of the strait etc, bottom drag coefficient.  For LGM we defined two overflow regions- Denmark strait and the Faroe Bank Channel (which is reduced from the PI configuration). Due to the larger than present day Greenland ice sheet, we were simulating anomalous high temperatures and salinity in the region of the Denmark strait. The overflow setup was modified in this region, the exact changes are very specific to the POP set up but we reduced the size of the overflow region from the PI set up to represent the narrower LGM configuration. We have modified the sentence as follows:

8: L130: “lower surface elevation”. It this a typo? I would expect a warming when there’s a lower surface elevation due to the lapse rate effect.

We agree this sentence was not clear and have revised as follows (see below) and added a new supplementary (Figure C3). which shows the LG-21ka -LGM-Zhu SAT difference and LG-21ka -LGM-Zhu elevation difference as seen by the atmosphere model.

“When comparing the LG-21ka results to LGM-Zhu and Osman, we find some notable differences (Fig.1cd) in the SAT which approximately correspond to the differences in elevation between the two simulations (fig.C2). The lower (colder) temperatures around the margin of the ice sheet are associated with the higher elevations (and vica versa). However, we note that this is not an exact 1-1 relationship”.

9: L135: Do you have a short explanation on why there’s a hemispherically asymmetric temperature differences between LG-21ka and Osman data?

We agree this is an interesting difference. The results across the southern hemisphere are similar as found compared to LGM-Zhu, which we attribute to the differences between the two ice sheet reconstructions (where LGM-Zhu and Osman have the same ice sheet). For this North hemisphere, especially the positive anomaly (over 5°C) across the Arctic we think this can in part be explained by  variation in the Osman ensemble. Our results only using the mean from the 500-model ensemble, but the standard deviation is largest at 21 ka across this region, up to 4°C.

10: L145-: I was confused with this paragraph. My understanding was that the differences in the orbital configuration between LG-21ka and PI do not cause large changes in the radiation budget, but differences in the ice sheet configuration and cloud radiative effect do. In this regard L146 “due to different orbital conditions” doesn’t read right to me. Please restructure the pargraph also following my General comment 1.

> Revising text. 

11: L165: On which height/pressure level is the geopotential height field analyzed?

We assessed the geopotential at an atmospheric pressure of 200hPa. We have added a note in the figure caption to clarify this, 

Note we analyze the geopotential height at an atmospheric pressure of 200 hPa

12:L259: The simulated sea ice field is compared with that of Paul et al. (2021). However, as far as I know, the reconstructions of sea ice area over the Northern Hemisphere during the Last Glacial Maximum (LGM) are not well-constrained. I'll leave it to he authors, but it might be useful to point this out.

We have added a sentence to clarify this as follows”

“However, we note that there is uncertainty in the GLOMAP sea-ice data that has not been quantified fully and we are using it as a guide to assess our simulations.”

Reviewer 2: 

Response to reviewers’ comments on manuscript “Surface mass balance and climate of the Last Glacial Maximum Northern hemisphere ice sheets: simulations with CESMv2.1.

The authors would like to thank the reviewer for their positive comments about our manuscript stating that it was ‘intriguing, especially regarding the analysis of the Surface Mass Balance (SMB) and its decomposition into different components. Additionally, the comparisons of SMBs across six different ice sheets are novel. The model also successfully reproduces global climate and oceanic conditions, including the Atlantic Meridional Overturning Circulation (AMOC)”

Each comment from the reviewer is shown in italic and our response follows. Any changes to the manuscript are given in blue

General comment:

1. Additional analysis on the cloud radiative effect (long wave) and its relation to SMB

A recent study by Gandy et al. (2023) discussed the role of clouds on the SMB over the southern margin of the North American ice sheet. They showed that many of the PMIP models show small cloud cover over the southern margin, which caused an intense downward shortwave radiation at the surface in that region. In this study, the authors have shown that the relatively large cloud cover at the southern margin plays a role in reflecting more sun light. On the other hand, they showed a relatively strong downward long-wave radiation, which could be relative to the existence of clouds or moisture.

Possible additional: I would suggest the authors to conduct further analysis on the role of clouds (especially on the long wave side) and add a discussion on its effect on SMB. Effects of clouds on the SMB is highly uncertain in general, and also thought to be important in the future Greenland (Gregory et al. 2020). Additional analysis on this point would make this study more interesting.

Thank you for highlighting the Gandy et al., 2023 paper, which was published after this manuscript submission. We have used this paper in the later point also in relation to causes of the negative SMB. The earlier figures of the radiation budget were perhaps misleading as we do not have a large cloud cover over the southern margin (see new Fig. C4g), but also have a small (or low) cloud cover as in the PMIP4 models. To discuss the role of clouds in  more detail, we have added an additional figure into the supplementary (Fig. C4) which compares the JJA anomalies (LG-21ka - PI) for a range of variables to investigate cloud forcing in relation to the SMB. Additionally, we have included the results for the comparison study, LGM-Zhu. We have added extra text into section 4.2 to discuss the influence of clouds on the SMB

Compared to the PI simulation, across the ice sheets there is an increase in the cloud fraction (i.e gets cloudier) but the clouds get thinner. These two specific changes in the nature of the clouds can be related to the earlier responses in radiation fluxes (Section 3.3). Thinner clouds act to increase the incoming solar radiation at the surface (Fig.C4a). Conversely, in cloudier areas, the clouds increase in incoming longwave radiation although the clouds are thinner (Fig.C4g.h).

1. 
   
   Implication to the evolution of British Ice sheet during the LGM

Gandy et al. (2018) showed that a warming of climate after 26ka is required to initiate a deglaciation of the British Ice sheet in their ice sheet model simulations. However, the sure about the cause of the warming remained unclear. This study nicely showed that the differences in orbital parameters between 21ka and 26ka can cause a warming and a decrease in SMB at 21ka over the British ice sheet. Adding a discussion on this point citing Gandy et al. (2018) would be beneficial and would also increase the impact of this paper.

Thank you for pointing out this study that we had not included in our paper. We have added a short note to highlight the results from the study as follows:

Gandy et al., 2018 concluded that a warming of the climate after 26 ka, and resultant reduction in SMB was in fact required to initiate the retreat of the BIIS at 21 ka. Therefore, the 1.5°C warming between LG-26ka and LG-21ka due to the change in orbital parameters may be one factor that led to the retreat of the BIIS, due to the increase in melt rate (figure7e).

3: Discussion on the cause of the large negative SMB at the southern margin of North American ice sheet

One of the advantages of this study is the use of a higher horizontal resolution in the AGCM compared to previous studies. First, I thought that this would improve the representation of the so-called stationary wave effect, which cools down the southern margin and increases the SMB. However, even with the higher resolution used in this model, it still simulates a very large negative SMB like other models. Please add a detailed discussion on what could be the potential cause of this negative SMB, e.g. the downscaling method, representations of clouds, biases in the climate model side and so on.

The reviewer has asked us to elaborate on the causes of the large ablation area that we still produce across the NAISC.  What is interesting is that the recent study by Gandy et al., 2023, also found similar results. This this paper, that authors concluded (Section 3: Collapse of the LGM ice sheet with a standard set up”) that if they initiated their LGM (21ka) simulated with a standard large NAISC extent as we have adopted in our study, the ice sheet rapidly deglaciated due to the large ablation area that forms across the southern margin very early on in the spin-up (see Fig.3 Gandy et al., 2023). These results appear to replicate the results from our study. Gandy et al., 2023 concluded that this behaviour was due to the fact ‘FAMOUS-ice’ had been heavily tuned for present-day Greenland and the default set up was not suitable for LGM (or cold) simulations.  This could also be a similar situation for our model set up. The CESMv2.1 model has been shown to have problems when applied to the LGM (Zhu et al., 2021; 2022) and required revised ‘tuning’ to comply with LGM GMST constraints. This implies a similar approach of ‘de-tuning’ as investigated by Gandy et al., 2023 will be required for LGM ice sheet simulations using CESM.  We do simulate a less negative SMB in the 26ka climate, but we have not adjusted the ice sheet extent to represent a smaller ice sheet.  We have added a revised section to highlight this study in Discussion section:

Our simulates SMB for the NAISC and BIIS appears too negative to prevent large marginal retreat if used as a forcing for an ice sheet dynamical model. ...

A recent study by Gandy et al., 2023 also investigated the LGM NAISC in a coupled climate ice sheet model (FAMOUS-ice). This study found that when initiating their simulations from a large NAISC (as in this study), a large ablation area formed across the southern margin of the ice sheet which led to rapid ice sheet retreat (see Fig3 Gandy). This behaviour was attribute to the heavy tuning of their model to present day Greenland ice sheet. The CESMv2.1 model has also been shown to have problems when applied the LGM climate (Zhu2021 ) and required de-tuning to comply with LGM GMST constraints. Future work investigating coupled ice sheet-climate simulation for the LGM with CESMv2.1 may also require de-tuning to correctly simulate the LGM North Hemisphere ice sheets.

4:I would suggest to move subsections 3.2 and 3.3 to the beginning of the section and combine with the analysis of global temperature. In this way, the authors can start the result section from general analysis (larger scale temperature, sea ice & AMOC), followed by detailed analysis on the atmospheric circulation, which then connects well to the SMB section. Hopefully this would increase the readability of the manuscript. –

We have revised the structure of the manuscript as suggested and hope this improves the readability

Minor changes all addressed:

1:L19: It was a bit hard to interpret this sentence. In particular, what do you mean by climate variability? 

 revised as follows:

“ Our SMB results are in qualitative agreement with the variations in temperature across the north hemisphere ice sheets”.

2:L36: Perhaps, you may also cite a recent study by UKESM (Smith et al. 2021) here.

Added

3:L76: Typo, two “stimulates” changed

4: L76: It might be written in a different paper, but could you clarify whether the precipitation field is downscaled or not, please?

We have added a note to clarify as follow on line 94:

“The 4km CISM2.1 grid (Fig. A1) provides high resolution information for CLM5’s elevation class information, as well as downscaled SMB (at 4 km resolution) by horizontal bilinear and vertical interpolation from the elevation classes. (Note at present precipitation is not downscaled)”

5: L80-85: Could you add a sentence explaining the reason for these modifications? Is it just from a practical reason or is there any science behind them.

To try and reduce both the spatial extent and amount of ablation across the ice sheets, we modified parameters within the land model, as listed. We have changed the sentence to add this justification as follows:

“We make two minor modifications to the defaults settings for the elevation class parameterizations with the aim of reducing the magnitude and extent of the ablation zone”

6: L107-109: Please add a sentence on what exactly was modified/adjusted here.

To run the POP2 model, the user must define “regions” or “zones of overflow”. For each overflow zone, the user can define several parameters, including source water thickness, latitude of flow, width of the strait etc, bottom drag coefficient.  For LGM we defined two overflow regions- Denmark strait and the Faroe Bank Channel (which is reduced from the PI configuration). Due to the larger than present day Greenland ice sheet, we were simulating anomalous high temperatures and salinity in the region of the Denmark strait. The overflow setup was modified in this region, the exact changes are very specific to the POP set up but we reduced the size of the overflow region from the PI set up to represent the narrower LGM configuration. We have modified the sentence as follows:

8: L130: “lower surface elevation”. It this a typo? I would expect a warming when there’s a lower surface elevation due to the lapse rate effect.

We agree this sentence was not clear and have revised as follows (see below) and added a new supplementary (Figure C3). which shows the LG-21ka -LGM-Zhu SAT difference and LG-21ka -LGM-Zhu elevation difference as seen by the atmosphere model.

“When comparing the LG-21ka results to LGM-Zhu and Osman, we find some notable differences (Fig.1cd) in the SAT which approximately correspond to the differences in elevation between the two simulations (fig.C2). The lower (colder) temperatures around the margin of the ice sheet are associated with the higher elevations (and vica versa). However, we note that this is not an exact 1-1 relationship”.

9: L135: Do you have a short explanation on why there’s a hemispherically asymmetric temperature differences between LG-21ka and Osman data?

We agree this is an interesting difference. The results across the southern hemisphere are similar as found compared to LGM-Zhu, which we attribute to the differences between the two ice sheet reconstructions (where LGM-Zhu and Osman have the same ice sheet). For this North hemisphere, especially the positive anomaly (over 5°C) across the Arctic we think this can in part be explained by  variation in the Osman ensemble. Our results only using the mean from the 500-model ensemble, but the standard deviation is largest at 21 ka across this region, up to 4°C.

10: L145-: I was confused with this paragraph. My understanding was that the differences in the orbital configuration between LG-21ka and PI do not cause large changes in the radiation budget, but differences in the ice sheet configuration and cloud radiative effect do. In this regard L146 “due to different orbital conditions” doesn’t read right to me. Please restructure the pargraph also following my General comment 1.

> Revising text. 

11: L165: On which height/pressure level is the geopotential height field analyzed?

We assessed the geopotential at an atmospheric pressure of 200hPa. We have added a note in the figure caption to clarify this, 

Note we analyze the geopotential height at an atmospheric pressure of 200 hPa

12:L259: The simulated sea ice field is compared with that of Paul et al. (2021). However, as far as I know, the reconstructions of sea ice area over the Northern Hemisphere during the Last Glacial Maximum (LGM) are not well-constrained. I'll leave it to he authors, but it might be useful to point this out.

We have added a sentence to clarify this as follows”

“However, we note that there is uncertainty in the GLOMAP sea-ice data that has not been quantified fully and we are using it as a guide to assess our simulations.”

Reviewer 3:

Response to reviewers’ comments on manuscript “Surface mass balance and climate of the Last Glacial Maximum Northern hemisphere ice sheets: simulations with CESMv2.1.

The authors would like to thank the reviewer for their positive comments and useful suggestions. Each comment is shown in italic and our response follows. Any changes to the manuscript are given in blue.

Main comments:

1: One of the main conclusions of this study is that "The large, simulated melt for the NAISC suggests potential biases in the climate simulation, ice sheet reconstruction, and/or highly non-equilibrated climate and ice sheet at the LGM time". These are all plausible. But little discussion has been devoted to the simulated melt for the BIIS, which is 3-4 times higher than that for the NAISC, and the simulated BIIS SMB is much more negative.

We agree that the study focused on the NAISC with little discussion of the BIIS. We have added in Section 4.3: LG-26ka versus LG-21ka surface mass balance the following discussion:

Gandy et al., 2018 concluded that a warming of the climate after 26 ka, and resultant reduction in SMB was in fact required to initiate the retreat of the BIIS at 21 ka. Therefore, the 1.5°C warming between LG-26ka and LG-21ka due to the change in orbital parameters may be one factor that led to the retreat of the BIIS, due to the increase in melt rate (figure7e).

2: l.9: It should be mentioned here that the acronym "NAISC" refers to the "North American ice sheet complex". 

Changed and acronym checked throughout

3: l. 37: The previous version CESM2.0 had a very high Equilibrium Climate Sensitivity (ECS). What is the ECS for CESM2.1?

We have revised the sentence as follows:

 “ This excessive cooling is due to a high equilibrium climate sensitivity of 5.3K that has been attributed to updates in the cloud parametrizations introduced in CAM6”

4: l.130: I think "lower surface elevation" (in LG-21ka compared to LGM-Zhu), based on the lapse-rate argument, should result in a warming not cooling over the ice sheets. It was mentioned previously that ICE6G in LGM-Zhu has a smaller GrIS, so should this be "higher surface elevation"?

Thank you for pointing this out. This was not very clear, and we have revised the text as follows and added a new figure into supplementary (Figure C3). which shows the LG-21ka -LGM-Zhu SAT difference and LG-21ka -LGM-Zhu elevation difference as seen by the atmosphere model. Although the relationship is not an exact 1-1, the region of higher elevation (due to the different input ice sheet reconstruction) does approximately correspond to the colder (lower temperatures).

Revised text: page 5: line 131:

“When comparing the LG-21ka results to LGM-Zhu and Osman, we find some notable spatial differences (Fig.1cd) in the SAT which approximately corresponds to the difference in elevation between the two simulations (Fig. C3). The lower (colder) temperatures around the margin of the ice sheet are associated with the higher elevations (and vica versa). However, we note that it is not an exact 1-1 relationship”.

5: l. 145-150: The difference in TOA SW_in between LG-21ka and PI is discussed here, but not shown in Figure 2. I would suggest adding a panel for TOA SW_in in Figure 2, which is useful for understanding the higher surface SW_in in summer months over high latitudes.

This figure (which is now figure 4) has been revised and the TOA SWin anomaly (LG-21ka - PI) added.

6: l.199: "... move, on average, more along the margin ..." doesn't read right.

We have revised this text as follows:

 “Over the SIS the synoptic systems move in relation to the southern margin of the ice sheet, while in the PI they move further north”.

7:l.225: It is not obvious that "In the Southern Hemisphere, differences with GLOMAP are smaller". As seen in Figure 5 (d), differences in sea ice extent are equally substantial for the SH, especially in February.  

We agree that this was not clear and have revised the text in this section as follows:

“Spatially, the differences are more complicated. During the summer in both the Northern and Southern Hemisphere LG-21ka overestimates the spatial extent of the sea ice (Fig.2a and c). For example, there are large areas of sea-ice across the Norwegian and Greenland seas, which are ice free in GLOMAP. In the winter months the simulation underestimates in these regions, but overestimates across the Bering Sea, Baffin Bay and into the Labrador Sea”.

8:l.241: A typo after "shallower". 

Revised

9’ l.298 "snowfall deposition": Should it be frost deposition?  

We have checked this as snowfall deposition is the correct terminology.

10: l.325-330: I think more discussion is needed to understand the features of SHF and LHF in Figure 9 (g, h). Both SHF and LHF are defined so that they are positive downwards, but it is intriguing that SHF is negative in the interior of NAISC -- meaning the ice surface heats the near-surface air in summer. I also have difficulty understanding why the LHF is positive along the southern margins of the ice sheets.

We also found the SHF results intriguing, as it implies that the ice sheet surface has a higher temperature than the air temperature. However, if you calculate the net radiation by adding the SWnet and LWnet from Figure 9, this does produce a positive radiation in those regions. Therefore, this could bring the ice sheet surface up to melting point (0°C) while the air temperature is below 0 (Fig. C2C).

The LHF is positive along the southern margins of the ice sheets due to prolonged bare-ice exposure during the summer. In this situation condensation occurs when the relatively warm and moist air is advected over the bare-ice exposed surface. This was discussed in more detail in Sellevold and Vizcaino, 2020. We have added a sentence in Section 4.2 to highlight this: 

The positive LHF anomaly along the southern margin of the ice sheets is due to prolonged bare-ice exposure; where when relatively warm and moist air is advected over this region condensation occurs (SellevoldandVizcaino2020). 

11:l.329: Please explain what is "ground heat flux (GHF)" and what physical processes are related to this term?

In CLM (the land model) the surface energy balance equation, for example as is given in Sellevold and Vizcaino 2020, is defined as

Melt Energy = SWnet  + LWnet +SHF +LHF + GHF

The GHF is a product of this equation. 

12: l.369-370: I agree that a natural follow-up work is to simulate ice sheet flow in an ice sheet model driven by the simulated SMB of this study. Perhaps the authors could mention here what they plan to do with the simulated ice sheet flow. From my perspective, it could provide an estimation of dynamic ice loss from these ice sheets, and (together with SMB) give a full account of ice sheet mass balance.

We have removed this sentence as it was a little vague.

Figure/table changes: 

Table 2: Please specify the latitudinal range of the "Tropical" region (for calculation of Tropical precipitation an Tropical SST)?

Thank you for pointing this out. We calculated these numbers in a latitude band between -30ºS and  30ºN and have added a note to the table to clarify this.

“Note a latitudinal range of 30°S to 30N was used for the tropical calculations”

Figure 4 (a) subtitle: Missing "water" after "Precipitable". 

changed.

Figure 8, 10: What are the dashed lines in panel (a)?

We have changed this figure so this should be clearer now.

Fig. C1: Based on the downward drifts of Global Surface Temperature, it seems the model was not fully equilibrated after running for 500 or 600 years. I understand that millennia may be needed for the ocean model to fully equilibrate, which requires much more computational resource. The authors may discuss whether this non-equilibrium would result in a substantial bias in the simulated AMOC.

The AMOC (Fig.3) reached a stable (or constant strength) in all three LGM simulations after ~ 300 years so ahead of the global mean surface temperature and we do not believe that AMOC displays trend. It may be the case, as the author states that for the deep ocean to reach a complete equilibrium state would require a simulation time of > 1000 yr (which was beyond the resources for this project). To assess if the atmosphere had reached equilibrium, we used the TOA net radiation (Fig.  C1b), which was a ~ 0 W/m². From Fig. C1 it may appear that the LG-26ka simulation was still cooling, but we assessed the ongoing trend of the GMST as the simulation progressed and the temperature had stabilised. As the discussion in Section 3.2 describes results for the AMOC was very different between different earth system models, even with very similar initial ice sheet extents.

Fig. C4 (a, b, c): There is a polynya-like feature between Greenland and Baffin Island/Labrador. Is that due to the intensive vertical mixing or deep water formation?

We agree that this is an interesting feature, but we were not able to determine the exact origin. We suspected it may be due to winds, however it is present in both LGM-Zhu and LG-21ka/LG-26ka.  We also found that interesting, given the differences between LG-21ka and LGM-Zhu in this region that the feature is present in both simulations.

Reviewer 3:

Response to reviewers’ comments on manuscript “Surface mass balance and climate of the Last Glacial Maximum Northern hemisphere ice sheets: simulations with CESMv2.1.

The authors would like to thank the reviewer for their positive comments and useful suggestions. Each comment is shown in italic and our response follows. Any changes to the manuscript are given in blue.

Main comments:

1: One of the main conclusions of this study is that "The large, simulated melt for the NAISC suggests potential biases in the climate simulation, ice sheet reconstruction, and/or highly non-equilibrated climate and ice sheet at the LGM time". These are all plausible. But little discussion has been devoted to the simulated melt for the BIIS, which is 3-4 times higher than that for the NAISC, and the simulated BIIS SMB is much more negative.

We agree that the study focused on the NAISC with little discussion of the BIIS. We have added in Section 4.3: LG-26ka versus LG-21ka surface mass balance the following discussion:

Gandy et al., 2018 concluded that a warming of the climate after 26 ka, and resultant reduction in SMB was in fact required to initiate the retreat of the BIIS at 21 ka. Therefore, the 1.5°C warming between LG-26ka and LG-21ka due to the change in orbital parameters may be one factor that led to the retreat of the BIIS, due to the increase in melt rate (figure7e).

2: l.9: It should be mentioned here that the acronym "NAISC" refers to the "North American ice sheet complex". 

Changed and acronym checked throughout

3: l. 37: The previous version CESM2.0 had a very high Equilibrium Climate Sensitivity (ECS). What is the ECS for CESM2.1?

We have revised the sentence as follows:

 “ This excessive cooling is due to a high equilibrium climate sensitivity of 5.3K that has been attributed to updates in the cloud parametrizations introduced in CAM6”

4: l.130: I think "lower surface elevation" (in LG-21ka compared to LGM-Zhu), based on the lapse-rate argument, should result in a warming not cooling over the ice sheets. It was mentioned previously that ICE6G in LGM-Zhu has a smaller GrIS, so should this be "higher surface elevation"?

Thank you for pointing this out. This was not very clear, and we have revised the text as follows and added a new figure into supplementary (Figure C3). which shows the LG-21ka -LGM-Zhu SAT difference and LG-21ka -LGM-Zhu elevation difference as seen by the atmosphere model. Although the relationship is not an exact 1-1, the region of higher elevation (due to the different input ice sheet reconstruction) does approximately correspond to the colder (lower temperatures).

Revised text: page 5: line 131:

“When comparing the LG-21ka results to LGM-Zhu and Osman, we find some notable spatial differences (Fig.1cd) in the SAT which approximately corresponds to the difference in elevation between the two simulations (Fig. C3). The lower (colder) temperatures around the margin of the ice sheet are associated with the higher elevations (and vica versa). However, we note that it is not an exact 1-1 relationship”.

5: l. 145-150: The difference in TOA SW_in between LG-21ka and PI is discussed here, but not shown in Figure 2. I would suggest adding a panel for TOA SW_in in Figure 2, which is useful for understanding the higher surface SW_in in summer months over high latitudes.

This figure (which is now figure 4) has been revised and the TOA SWin anomaly (LG-21ka - PI) added.

6: l.199: "... move, on average, more along the margin ..." doesn't read right.

We have revised this text as follows:

 “Over the SIS the synoptic systems move in relation to the southern margin of the ice sheet, while in the PI they move further north”.

7:l.225: It is not obvious that "In the Southern Hemisphere, differences with GLOMAP are smaller". As seen in Figure 5 (d), differences in sea ice extent are equally substantial for the SH, especially in February.  

We agree that this was not clear and have revised the text in this section as follows:

“Spatially, the differences are more complicated. During the summer in both the Northern and Southern Hemisphere LG-21ka overestimates the spatial extent of the sea ice (Fig.2a and c). For example, there are large areas of sea-ice across the Norwegian and Greenland seas, which are ice free in GLOMAP. In the winter months the simulation underestimates in these regions, but overestimates across the Bering Sea, Baffin Bay and into the Labrador Sea”.

8:l.241: A typo after "shallower". 

Revised

9’ l.298 "snowfall deposition": Should it be frost deposition?  

We have checked this as snowfall deposition is the correct terminology.

10: l.325-330: I think more discussion is needed to understand the features of SHF and LHF in Figure 9 (g, h). Both SHF and LHF are defined so that they are positive downwards, but it is intriguing that SHF is negative in the interior of NAISC -- meaning the ice surface heats the near-surface air in summer. I also have difficulty understanding why the LHF is positive along the southern margins of the ice sheets.

We also found the SHF results intriguing, as it implies that the ice sheet surface has a higher temperature than the air temperature. However, if you calculate the net radiation by adding the SWnet and LWnet from Figure 9, this does produce a positive radiation in those regions. Therefore, this could bring the ice sheet surface up to melting point (0°C) while the air temperature is below 0 (Fig. C2C).

The LHF is positive along the southern margins of the ice sheets due to prolonged bare-ice exposure during the summer. In this situation condensation occurs when the relatively warm and moist air is advected over the bare-ice exposed surface. This was discussed in more detail in Sellevold and Vizcaino, 2020. We have added a sentence in Section 4.2 to highlight this: 

The positive LHF anomaly along the southern margin of the ice sheets is due to prolonged bare-ice exposure; where when relatively warm and moist air is advected over this region condensation occurs (SellevoldandVizcaino2020). 

11:l.329: Please explain what is "ground heat flux (GHF)" and what physical processes are related to this term?

In CLM (the land model) the surface energy balance equation, for example as is given in Sellevold and Vizcaino 2020, is defined as

Melt Energy = SWnet  + LWnet +SHF +LHF + GHF

The GHF is a product of this equation. 

12: l.369-370: I agree that a natural follow-up work is to simulate ice sheet flow in an ice sheet model driven by the simulated SMB of this study. Perhaps the authors could mention here what they plan to do with the simulated ice sheet flow. From my perspective, it could provide an estimation of dynamic ice loss from these ice sheets, and (together with SMB) give a full account of ice sheet mass balance.

We have removed this sentence as it was a little vague.

Figure/table changes: 

Table 2: Please specify the latitudinal range of the "Tropical" region (for calculation of Tropical precipitation an Tropical SST)?

Thank you for pointing this out. We calculated these numbers in a latitude band between -30ºS and  30ºN and have added a note to the table to clarify this.

“Note a latitudinal range of 30°S to 30N was used for the tropical calculations”

Figure 4 (a) subtitle: Missing "water" after "Precipitable". 

changed.

Figure 8, 10: What are the dashed lines in panel (a)?

We have changed this figure so this should be clearer now.

Fig. C1: Based on the downward drifts of Global Surface Temperature, it seems the model was not fully equilibrated after running for 500 or 600 years. I understand that millennia may be needed for the ocean model to fully equilibrate, which requires much more computational resource. The authors may discuss whether this non-equilibrium would result in a substantial bias in the simulated AMOC.

The AMOC (Fig.3) reached a stable (or constant strength) in all three LGM simulations after ~ 300 years so ahead of the global mean surface temperature and we do not believe that AMOC displays trend. It may be the case, as the author states that for the deep ocean to reach a complete equilibrium state would require a simulation time of > 1000 yr (which was beyond the resources for this project). To assess if the atmosphere had reached equilibrium, we used the TOA net radiation (Fig.  C1b), which was a ~ 0 W/m². From Fig. C1 it may appear that the LG-26ka simulation was still cooling, but we assessed the ongoing trend of the GMST as the simulation progressed and the temperature had stabilised. As the discussion in Section 3.2 describes results for the AMOC was very different between different earth system models, even with very similar initial ice sheet extents.

Fig. C4 (a, b, c): There is a polynya-like feature between Greenland and Baffin Island/Labrador. Is that due to the intensive vertical mixing or deep water formation?

We agree that this is an interesting feature, but we were not able to determine the exact origin. We suspected it may be due to winds, however it is present in both LGM-Zhu and LG-21ka/LG-26ka.  We also found that interesting, given the differences between LG-21ka and LGM-Zhu in this region that the feature is present in both simulations.