This paper describes a new feature developed for the Ice-sheet and Sea-level System Model (ISSM), which adds stochastic parameterizations for particular forcings and model parameters. The new model framework is called The Stochastic Ice-Sheet and Sea-Level System Model v1.0 (StISSM v1.0). Overall, the paper is very clear in its structure and explanation of the purpose of the new functionality, an overview of how to use the new functionality, as well as interpretation of some initial experiments that used the new functionality. The authors should be commended for explaining a complicated concept (adding stochasticity to an already existing model framework) in a very clear and straightforward way. Additionally, the paper makes use of existing "standard" ice sheet model setups (e.g., MISMIP+ and IQIS) and makes good connections with prior related literature (e.g., the comparison with Roe and Baker, 2016).

I have some general comments as well as editorial comments further below, amounting to minor revisions needed before resubmission. I will emphasize again that the paper is very well written and presented in a clear manner. The StISSM model will be very impactful for future studies and I am pleased to see this significant step forward in ice sheet modeling.

General comments:

1. In the introduction, I suggest some additional text to state that, although StISSM will alleviate the need to run large ensembles of climate models if it is possible to correctly parameterize the structure of internal climate variability from other sources of information, such as observations. In this current manuscript draft, the wording in the introduction too strongly claims that StISSM will eliminate the need to run large GCM ensembles. I think there still may be a need, although that need could very well be reduced by StISSM. This will also connect nicely with the discussion text on lines 553-554.
2. In the future, is it planned to add functionality to specify different stochastic time steps for different input variables (if they are uncorrelated)? I suggest adding some text to either the methods or the discussion sections to touch on this.
3. I suggest making the connection between "y" and "eta" more clear. I think it's good to have two separate symbols to make it clear that one is coming from an autoregressive process. But it would be good to state very explicitly: "At each simulation time step, a value for either y or \eta is calculated and used for that particular time step and subsequent simulation time steps until the next stochastic update." That wording is a little clunky and can definitely be improved. My suggestion is just to make it clear that "y" and "eta" represent a similar thing: the realization of a random variable that is used as the value for a particular forcing or parameter by ISSM.
4. I suggest changing the presentation of the changes in ice mass in Section 4 from showing the initial mass and final mean masses in each ensemble to showing the mean changes in mass. In other words, show just the differences between the initial mass and the final mean masses in Gt, as well as the percent changes (as you have already shown). I don't see the need to show the initial and final masses themselves; the differences will illustrate the results more clearly. This would also make it easier to compare the mean mass change against the deterministic drift in Table 6. Additionally, please change how this is shown in the figures (e.g., Fig 4c, Fig 6c, Fig 8b).
5. The "Code and data availability" section states that "the simulation results, and the scripts to reproduce all the figures are available" on Zenodo. Are the scripts user to initialize, configure, and run the ISSM simulations also available there. The policy states that "preprocessing, run control and postprocessing scripts" and, I do see postprocessing scripts in the Zenodo archive but I don't see preprocessing and run control scripts. If these are there, please ignore this comment. If they are missing, please provide these in the same Zenodo repository.

Editorial comments:

1. line 9: Change "of" to "for"
2. line 29: It'll be a mouthful but I would spell out "CMIP6" here.
3. line 31: Change "AR6" to "Assessment Report 6 (AR6)"
4. line 32: Change "inclusion" to "selection"
5. lines 82-84: This statement is brought up as motivation for the paper but it's not really addressed: "Finally, climate model simulations are generally not coupled to ISMs, which neglects possible impacts from ice-sheet changes on the climate system, such as surface elevation changes and modified ice discharge into the ocean." I suggest removing this from the intro because it doesn't directly motivate the need for a stochastic ice sheet model. Alternatively, if you'd like to keep it in, I suggest adding discussion text on how StISSM can help address the coupling issue.
6. line 96: Possibly change "underline" to "emphasize" for clarity
7. line 101: Add "The new ..." before the start of this first sentence to make it clear that this is referring to the new stochastic functionality within ISSM and not something that had already existed.
8. line 106: I might suggest changing "ocean forcing" to "frontal ablation" or something like that. I think of "ocean forcing" as a climatic forcing (i.e., ocean thermal forcing), whereas the way that you convert ocean forcing to frontal melt (via parameterization) is a glaciological process.
9. line 151: This refers to \epsilon_t but there's no such variable in the preceding section.
10. lines 151-152: It should be made clearer here and in the preceding section which user-specified variables can vary in (1) space and (2) time. In the regular version of ISSM, I believe that all of the mean fields (melt rate, water pressure, and TF) can be specified as varying in both space and time. My understanding of Section 2.2 indicates that, in StISSM, the user can specify \sigma_y to vary in both space and time, as well. I suggest adding this statement to the paragraph on lines 115-122 and reference Section 2.2 from there. That would make it clear up front to the reader.
11. lines 347-348: The extrapolation of C_B is a bit unclear to me. Is this needed because the ice sheet will grow in extent during the transient simulation to get to steady-state? In other words, this is an extrapolation of C_B beyond the extent of the present-day ice sheet, where there are no ice velocities available to invert for C_B, correct? Please clarify.
12. line 360: "Free-flow boundary condition" is the same as "Neumann boundary condition", correct? If so, please state this.
13. line 383: I suggest changing "to quantify the minimal amount of deterministic model drift" to "to quantify the amount of deterministic model drift, which is minimal", if that is indeed what is meant here.
14. Figure 3: Please add ticks and axis labels and also make sure that the axes are equal so that Greenland doesn't appear stretched in one direction or the other.
15. line 408: Please add a very brief explanation for what the Shapiro-Wilk test signifies and how to interpret the p-value.
16. line 410: Please provide support for the statement: "and the PDFs of final glacier state have not yet converged to statistical steady-states." Is this determined by looking at the changes in the PDF statistics (mean, std dev, skew) over the last X years? Please specify.
17. line 413: Change "combined to" to "due to"
18. line 570: Would it be fair to make "experiments" more specific by changing to "laboratory experiments"? If so, please make that change.
19. line 593: I suggest elaborating on the statement that stochastic forcing causes "asymmetry in the response." From the experiments presented, it seems to me that there is asymmetry during the transient but that the asymmetry decreases towards the end of each simulation and, as demonstrated by the Shapiro-Wilk p-values, ends up being fairly close to a symmetric normal distribution. If I am misinterpreting, it is because of my lack of familiarity with the Shapiro-Wilk test and that should be addressed in the paper (I have a comment about this above). But if what I wrote is correct, I suggest expanding this conclusion to state something similar to what I have suggested.

The goal of the paper is to present the recently developed Stochastic Ice-Sheet and Sea-Level System Model (StISSM), which adds stochastic parametrizations and tools to handle ensemble runs to the ISSM model. Available parametrizations are based on Gaussian noise, added at each (stochastic) time-step, and autoregressive processes, which can be used to represent surface mass balance and ocean forcing. The model capabilities are showcased targeting idealized experiments (including MISMIP+) and the Greenland ice sheet (GrIS).

Accounting for uncertainty in ice-sheet modeling is of paramount importance and tools like the one presented here are important and worthy of publication. However, I find the exposition hard to follow. Details about the stochastic model are at times buried in the presentation of the numerical experiments and some important details are missing. In my understanding the StISSM provides two stochastic processes in time (potentially different in each sub-domain): an autoregressive process and an one based on Gaussian noise. It also allows to run ensemble members in parallel, although it is not clear how the computational resources (cores, memory) are used. Statistical quantities (e.g. moments, p-values) are computed and reported in the Results section, but it is not clear whether these are computed directly by StISSM or how the data from the ensembles are collected. Finally, how is the stochastic layer implemented? How is it coupled to ISSM? Is it a driver to ISSM? Is it implemented in C++, Python or other languages?

A good part of the paper is devoted to using the StISSM for applying different stochastic parametrizations to two synthetic ice problems and a real ice sheet (Greenland). The numerical experiments are well thought out but I don't think it is particularly useful to target three different applications. I think it would have been better to target only one application (maybe a glacier) and show the effect of different choices of the parametrizations on the glacier evolution and mass balance. Targeting different (and complex) problems makes it harder to understand the impact of different parametrizations, without adding much in terms of explaining or demonstrating the stochastic model. Moreover, given that the main novelty introduced by StISSM is that it provides parametrizations, I would have expected more emphasis on 1) why to choose a specific stochastic (e.g autoregressive) processes for modeling, e.g., the surface mass balance (SMB), 2) how the stochastic process compares, in a statistical sense, with available time-series of SMB, 3) what is the impact of using a first-order versus an higher-order autoregressive process, and so on.

I'm not asking the authors to completely change the numerical experiments, but I would encourage them to significantly improve the presentation, clearly exposing the additional capabilities introduced by the StISSM, adding some important software and computational details and expanding on how to choose (calibrate?) the stochastic parametrizations for different forcing/parameters. Some additional details of the numerical experiments, that are important for completeness and reproducibility but that are not essential for explaining the new capabilities of StISSM ,could be moved to the appendix.

Here are some additional comments:

• Eq. (1). Is the Gaussian noise uniform in space even if the mean value is not? Please specify this in the text and discuss this choice. Please specify that the "mean" is intended in time, not in space (if I understand correctly)
• line 153: would it make sense to have different stochastic time steps for different parameters?
• line 164: OK, so the spatial stochasticity is introduced at the subdomain level. This should have been explained before, in the introduction and, in more detail in section 2.1 where it should be explained that eq. (1) is at the subdomain level.
• eq. (7): I do not fully understand the purpose of the intercept and trend terms. Also, what is the choice for beta_0 and beta_1 in the numerical experiments in sections 3?
• line 235: can you detail how you manage resources (nodes, cores, memory) when you run in parallel multiple members of the ensemble (each of the members might need to be distributed on several ranks). Do you use any strategy to reduce I/O and storage when running large ensembles? Any strategy to monitor the runs (e.g. what happens if a few of the 500 simulations in the ensemble fail?)
• eq (8) and (11): this is very minor, but I think that the use of squared terms "C_W²" and "C_B²" is poor notation. I know it is somewhat common, but it is misleading because it sort of implies that  C_W and C_B have some physical meaning. Using the square to denote positive quantities (if that's the reason for the square) is hardly defensible because there are a lot of other physical variables (e.g. thickness) or coefficients (flow factor) that are positive (or nonnegative) and they are not denoted with a square of some other quantity. I would suggest dropping the square and using directly the coefficient C_W and C_B.
• Sections 3 and 4: The rigid separation of the "Model experiments" and "Results" sections makes it harder to follow the exposition. I think that the Results part should follow the Model Experiments part for each of the three examples.