General comments

In this paper, the authors develop the first large ensemble calibrated with observations that resolves individual outlet glaciers, using a variable mesh that achieves resolutions of less than 1 km at the margins. They employ the Stochastic Ice-Sheet and Sea-Level System Model (StISSM) to simulate the Greenland Ice Sheet (GrIS). This study quantifies the effect of internal climate variability on the evolution of the GrIS over 185 years under two climate scenarios: RCP8.5 (WARM) and pre-2000 (CTRL) conditions. The authors use the AWI-ESM model for atmospheric and oceanic forcing, the dEBM for downscaling atmospheric forcing to 5 km, and introduce a stochastic component to the surface mass balance (SMB), thermal ocean forcing (TF), and runoff, both individually and in combination with different correlations. The study provides interesting conclusions and certainly relevant to the purpose of TC and is a valuable contribution to understanding the effects of stochastic climate variability on the GrIS evolution. Additionally, its clear and precise writing is appreciated. Therefore, I recommend the publication of this paper in TC after minor revisions, which I outline below.

Regarding the calibration and the initialization, I would like to raise a few points for discussion:

1. The calibration is performed using a deterministic simulation, whereas most of the simulations in this work include a stochastic component in the forcing. Why was this choice made? A brief justification in the text would be greatly appreciated.

2. Figure 2 shows the mass change estimated by the IMBIE (Otosaka et al., 2023) and that simulated by the model. In L317, the authors state that “the total modeled 2007-2017 mass change agrees with the observational record within uncertainty ranges, even though the modeled mass loss rate is over- and under-estimated in the early and later years [..]”. This overestimation in the model's mass change in some years (~2008–2011) appears to exceed the IMBIE uncertainty range by a significant margin (up to ~400 Gt in some instances).

As described by the authors, there is strong agreement between the model and observations in the later years of the simulation (i.e. in the total mass change for this period, which was the goal of the calibration), but the rate of mass loss in the model appears to be steeper than in the observations and this could affect the future evolution of the ice sheet. Indeed, as discussed in the Discussion section (L672), this calibration is responsible for the loss of mass in the early decades of the CTRL-LE simulations, which is the most abrupt change occurring in all the simulations (both in the WARM and in the CTRL scenarios). 

Given these points, I believe it would be beneficial to provide additional justification for the choice of this calibration, clarify the origins of the differences with observations when this difference is maximal (both spatially and in terms of physical mechanisms or model uncertainties.), and explain why this particular calibration approach was adopted despite the noted discrepancies with Otosaka et al. (2023) and the effects in the CTRL-LE simulations in the following years. Additionally, in the discussion (L691–L694), it is mentioned that the results for glaciers with more abrupt retreat should be interpreted with caution. I believe this paragraph could be expanded further, providing a more critical reflection on the results themselves and elaborating more on their limitations.

3. While the glacier retreat shown in Figure 3 is illustrative of the present-day performance, it would be valuable to also assess the model’s pre-2000 performance against observations using a 2D plot of ice thickness and surface velocities, perhaps as an appendix. I am surprised that the western margin experiences almost no ice mass loss until 2050 in any of the simulations, considering that this region is currently one of the areas experiencing the highest ice loss (Mouginot et al., 2019). However, in the Calibration section, the authors note that they are unable to accurately represent the retreat of the SK glacier (Jakobshavn), and in the Discussion section, they further mention that the results for this glacier should be interpreted with caution. It would be helpful if the authors could elaborate on why they believe the ice sheet retreat in this region is not being adequately simulated. 

Many results are presented, with extensive discussion on the spread within different ensembles, the factors that increase dispersion, when internal climate variability is more relevant, and the types of stochastic forcing that have the greatest impact on ice sheet development. These results are important and provide clarity on the problem and quantify the uncertainties associated with internal climate variability. However, one of the main scientific questions of the paper is the actual effect of adding stochastic forcing on the evolution of the ice sheet. For this reason, a more thorough comparison between the ensembles and the deterministic simulations is missing.

Therefore, it would be interesting to assign greater emphasis to the deterministic simulations shown at the beginning of the Results section and to compare the large ensembles with them in greater detail. In this way we could better see the differences between using a purely deterministic forcing and adding a stochastic component. Some recommendations for this include (which could also be incorporated into a separate figure or presented differently):

• Include these simulations in Table 1 where the experiments are summarized, using distinct names (e.g., CTRL-Det and WARM-Det).
• Include their profiles in Figure 6.
• In Figures 7 and 8, include the mass change (or the difference relative to the ensemble mean) for CTRL-Det and WARM-Det in the years 2050, 2100, and 2203, and discuss any spatial differences between the deterministic and stochastic simulations.
• Include in the figure 10 the value of the total ice mass change of CTRL-Det and WARM-Det.

Specific comments

I will now provide a series of more specific comments on certain parts of the article that I believe could be improved.

L20: When I read the abstract, I did not fully understand the mention of the Antarctic Ice Sheet until I read the entire article. Therefore, I recommend either removing it from the abstract or adding a sentence explaining why studying internal variability in Antarctica is important, as done in the Implications section.

Regarding the CTRL simulations, I understand from the description in the Methods section that the forcing applied in these simulations is the same as in the calibration for the years 2007–2017, and that afterward, the mean forcing from 1850–1999 is used. Does this mean that an instantaneous cooling occurs in 2018 (with a drop in TF and runoff and an increase in SMB), returning to pre-2000 values? It would be helpful if authors included the time series of TF, SMB, and runoff applied to the CTRL simulations in Figure F1.

Figure 4c: It is not entirely clear whether the peak around the year 2040 in both curves (WARM and CTRL) coincides exactly. It would be helpful to clarify this in the text or slightly adjust the way it is plotted.

Figure 5: This figure is somewhat unclear, making it difficult to extract information from it. Additionally, compared to Figures 7 and 8, the only extra information it provides is a higher temporal resolution of the results glacier by glacier. Therefore, I would recommend modifying it for greater clarity. Below are some suggestions:

• I am not entirely sure what is shown in the color bars. It represents the terminus retreat, but some glaciers (Figures 5a and 5b) appear in red and then turn gray, which, according to my understanding of the figure caption, would mean they initially retreat and later experience no further advance or retreat. Meanwhile, other glaciers (such as those in the northwest) remain dark red throughout the entire time series until the year 2203. What does it mean when they stay dark red for the entire period? Additionally, it seems that many glaciers retreat by around 50–100 km before 2050, which is a surprisingly large value.
• In L485–L486, it is mentioned that the glaciers in the figure are numbered, but I am unsure what this refers to or whether the numbering is located elsewhere. In any case, it would be helpful to include a clear indication in the figure referencing these glaciers. Additionally, it would be useful to label the most relevant glaciers in the same way they are highlighted in Figure 7.
• The GrIS is divided into the regions indicated in the left margin; including a small map showing this division would help better locate the glaciers and improve the overall understanding of the figure.
• Finally, to maintain consistency with Table 1, I believe the figure titles should be labeled as CTRL-LE and WARM-LE instead of CTRL and WARM.

Regarding the results of the small ensembles shown in Figure 10, it would also be valuable to include the time series of ice mass change from these experiments for better visualization, either in the same figure or in E1. 

In L615–L617, the effects of oceanic variability are discussed. Although they are generally smaller than those of atmospheric variability, they prevail in the western region. It would be interesting to comment on why oceanic variability has a greater influence in these areas. This is expected, as shown in studies such as Slater and Straneo (2022), since although atmospheric forcing currently dominates GrIS mass loss, ocean warming has a greater influence on glaciers in the west and south.

Technical corrections

L755 The creation of an open large ensemble dataset for the community is very much welcomed and useful; however, the link (https://doi.org/10.18739/A2VX0651F) to the repository does not work (last checked on February 22), and this should be fixed before the paper is published.

In some figures I have the impression that the limits in the colorbar have not been well applied and it appears cut in the Petermann and Zachariae Isstrom glaciers (figures 7d, 7e, 7f, 8d, 8e, 8f, 9, 11, 14).

L210: There is a mistake in the sentence “correlation patterns are strong within East- and West-Greenland, but weaker between East- and West-Greenland”. Where exactly are the patterns strong?

L319: “Figure 3 compares the observed and modeled 2007-2017 retreat rates”. However, the figure does not show retreat rates, it shows the total retreat in kilometers.

L177 L196 L316 L319 L374 L405 L413 L418 L426 L430 L482 L492 L526 L563 L777 L802 Instead of “Figure” it should be “Fig.”. 

References

Mouginot, E. Rignot, A.A. Bjørk, M. van den Broeke, R. Millan, M. Morlighem, B. Noël, B. Scheuchl, & M. Wood, Forty-six years of Greenland Ice Sheet mass balance from 1972 to 2018, Proc. Natl. Acad. Sci. U.S.A. 116 (19) 9239-9244, https://doi.org/10.1073/pnas.1904242116 (2019).

Slater, D.A., Straneo, F. Submarine melting of glaciers in Greenland amplified by atmospheric warming. Nat. Geosci. 15, 794–799 (2022). https://doi.org/10.1038/s41561-022-01035-9

In `The Greenland Ice Sheet Large Ensemble: Simulating the future of Greenland under climate variability,' Verjans and co-authors use a stochastic variant of the Ice Sheet System Model to explore the sensitivity of the Greenland Ice Sheet to variability in oceanic and surface mass balance forcing.  In particular, they aim to quantify the relative importance of such so-called `aleatoric' uncertainty relative to other types of uncertainty derived from imperfect or unresolved modeling assumptions and initial conditions.  Through a detailed comparison of ensemble experiments meant to represent both a continuation of contemporary forcing alongside a potential high end warming scenario, they find that the influence of stochastic climate is non-negligible over the coming two or so decades (in terms of total predicted mass change), while these stochastic effects become relatively unimportant over century-scales.  This result is interesting (albeit not particularly surprising) in that it illuminates a principal challenge for short term sea-level prediction, while providing some important guidance as to whether short time-scale variability represents a source of uncertainty that needs to be better quantified for long term projection (thankfully not, it seems!).  

This manuscript represents an impressive and insightful culmination of several methodological threads that seem to have been `in the works' for a few years -- the development of StISSM and its ensemble generation tools, the statistical characterization of climate variability in a generative sense, and the coupling of ice dynamics to downscaled surface mass balance and frontal ablation paramterizations.  The current work is undoubtedly at the vanguard of ensemble methods for ice sheet uncertainty quantification, and a big step forward for understanding Greenland's sensitivity to climate noise.  I have no issues with the paper's general methodology.  I have included below a few comments that I hope can improve the manuscript's clarity and utility.    

L42: `are performed' should be `have been performed' for consistent case.

L46: Perhaps here, perhaps elsewhere, it's maybe worth providing a higher level overview of where climate stochasticity comes from (and where it does not).  In particular, it's worth noting that climate is very likely not actually random, but rather appears that way due to the chaotic dynamics characteristic of the atmosphere (and ocean, to a lesser extent).  Ice sheets do not exhibit such ostensible stochasticity (EISMINT2 and ice streams notwithstanding), so the irreducible uncertainty in the ice sheet context is derived solely from the forcing term.

L84: What question is being referred to here?

L143: A qualitative description of what EN4 is, and why it's helpful for bias correcting the ocean thermal forcing would be very helpful.  

L202: I was expecting a similar interpretation of the moving average component of the fits.  Do these exhibit any interesting patterns?  Does the MA component even matter?  

L203: I spend more time than most glaciologists thinking about covariance, and yet I'm still unclear as to what's going on here.  In particular, after fitting the ARMA model to each time series of TF, SMB, and runoff independently, how are spatio-temporal correlations between them calculated.  Reading the appendix, it seems that there are three layers to this model: Fitting a piecewise linear function, fitting an ARMA model to each basin/variable, and then computing a big covariance matrix between the residuals for all?  Okay, I guess, but I would like a more centralized and coherent justification for why this is a reasonable way to control the spatial relationships.  

Sec 2.2.3: I'm not completely sure that this is the right thing to do, but it might be helpful to lead the section with this (which is essentially the `physics'), so that the reader will have a better idea of what the TF, etc. is going to be used for.  Similarly, you might include here the way that lapse-rates and such enter the SMB calculation.  

L264: I'm sympathetic to the need to use SSA for computational reasons, but it would be worthwhile to briefly describe the implications -- Greenland has a lot of ice that is very much not consistent with the assumptions of that model after all.  

Eq. 3 and lines after: Am I missing a previous point at which $N$ is defined?  How is it computed here?  Constant fraction of overburden?

L271--273: Would it be possible to provide some additional justification with respect to the linear regression step described here?  This isn't something I've seen before, so it would be nice to understand a little bit better how/whether this works.

Sec. 2.3.2: I am confused as to the technical approach for performing this calibration.  Is this done by manually fiddling with $\sigma_{max}$ until the eyeball norm is minimized, or is there an objective (and automated) procedure that is taking place?  

L399: I am surprised that the assertion that the small deviation of the ensemble mean from the deterministic run is a result of noise-induced drift is not backed up by a statistical test.  It would strengthen the argument to include a test of significance here.  

Fig. 4d: This is a challenging metric to use in order to assert the relative importance of uncertainty because the denominator gets so very small close to the start of the simulation period.  I am not sure what the alternative is, but it might be helpful to acknowledge that.  

L465: delete `briefly'.

L544: This is a pretty awkward sentence -- suggest rephrasing.

Discussion: I appreciate the comparison to both Tsai (for the forcing uncertainty comparison) and ISMIP6 (for the model uncertainty comparison), but it might be useful to also compare to some of the previous works that explore parametric uncertainty -- which seems to be of similar size to model uncertainty in some cases.  Would using randomly sampled climate-to-SMB parameters drown out the influence of the stochastic climate?  This would be important to know in making a decision about whether to include stochastic forcing in, say, ISMIP7.