Summary 

-------

The paper presents a method to synthetically generate surface mass balance for the Greenland ice sheet based on existing process model output. It allows to statistically produce SMB with the same overall temporal and spatial characteristics, but with different inter-annual variability. The parameterisation may be used to generate a large number of forcing to study the effect of internal variability.       

I didn't review the provided example code in detail, but note that it is easy to access, well documented and facilitates reproduction of the presented material and the figures. This is a very nice example of open science. Thank you to the authors.    

General comments

----------------

My main general comment is that the paper suggests (explicitly and implicitly, e.g. l3, l38, l58 and by the referenced papars) that the provided parameterisation of SMB could be used for future projections. However, the presented methods focus on stationary signals over the historical period. While we can imagine similar methods could be applied for projections, no evidence is given that this is actually the case. Additional considerations would certainly be needed to get there. This is mentioned in the discussion, but until then the reader may be mislead. I therefore suggest to either 1) make it clear from the beginning (abstract, introduction) that the given approach is only valid for SMB close to the historical, or 2) provide examples of an application for the multi-decadal to centennial times scale.     

I am concerned if the set of regions (Fig 1) is a meaningful and best choice for that type of analysis. The catchment delineation was motivated by the ice flow, not by SMB at all. The sizes of the catchments are all different, they encompass regions of largely different SMB, the relative share of accumulation vs ablation areas is different from region to region. 

It is not immediately clear how an optimal delineation would look for SMB like instead, but it seems that something more adapted in at least some of the aspects above could be constructed. I am thinking e.g. of one of the higher order mascon delineations in https://doi.org/10.1093/gji/ggy242 or similar.  

My main concern is in how much your results depend on that specific delineation. You mentioned that problem in in the discussion (l300), but did not provide further analysis about it. I think it would be useful to compare to at least one other definition of regions to show robustness of the results.

The model fitting process is not very transparent, because it relies on existing software packages and is not well described beyond that. I belive my confusing about the fitting process (comments l126 and l155) is related to that lack of description. I'd suggest to add more detail about how the fitting process works, what characteristics of the original time series are supposedly matched with the selected models and to illustrate these relationships with a figure if needed.  

And maybe a bit speculative, does any of the models significantly outperform a random (white noise) process scaled to the right amplitude per basin and model? 

Also, while I can see how the noise characteristics may slightly differ and match better or worse, how do projected mass changes produced by an ice sheet model under these different forcing attempts compare. Since this method is already integrated in ISSM, could you show some results? How does your parameterisation with the best statistical model compare to sampling individual years randomly, which is what people have tried in the absence of an appropriate statistical tool?   

Specific comments

-----------------

Title: Given that the presented parameterisation is specifically developed for the Greenland ice sheet, I would suggest to add 'Greenland' in the title. Clarifying this may also be important in light of contrasting results found from Antarctica as discussed.  

l1 "Many scientific and societal questions that draw on ice sheet modelling could be best addressed by sampling a wide range of potential climatic changes and realizations of internal climate variability"

I can think of many questions that were better addressed by focussing on a specific scenario or where internal climate variability is not that relevant. Suggest to remove the superlative "best" from this sentence. 

l5 "The wide sampling required to address such questions is computationally infeasible with sophisticated numerical climate models at the resolution required to accurately force ice sheet models. Stochastic generation of climate forcing of ice sheets offers a complementary alternative."

The presented parameterisation is limited by the available high-resolution climate model simulations for a specific time-period and scenario. The added value is mainly in generating internal variability, which is an important but not sufficient ingredient to produce a wide range of forcing. The limitations of the method should be recognised and made clear in the text.    

l11 "while including feedbacks to ice sheet surface elevation"

More than one feedback? It is basically the SMB-height feedback. Reformulate.

l38 "providing too little information to estimate the probability distribution of output variables such as future sea level"

This is the reason process-based ice sheet sea-level projections had to be emulated for the AR6 assessment. Could be added with a reference to Edwards et al. (2021).  

l43 "ice sheet models are not typically forced directly by climate model output."

Suggest to add "global" before "climate model output" to distinguish between GCMs/ESMs and RCMs. The latter also being climate models that are used to force ice sheet models directly. 

l44 Maybe mention "downscaling" here and be more specific. "Rather, global climate model output must be downscaled to construct an SMB field of high enough spatial resolution and quality, often through use of a specialized mass and energy balance model that accounts for processes at the snow/ice surface and in the snowpack".

l48 "Stochastic methods provide a low-cost alternative to sophisticated process models".

I find this sentence somewhat in contrast to l292 " ... this framework can be applied to new SMB process model outputs as they become available". While l48 can be read to mean that stochastic methods can replace sophisticated process models altogether, l292 recognises that your method is based on available process model results. Suggest to reformulate l48.

l54 Suggest to distinguish between the generator and the SMB product it generates. The latter is what the 'which' in the sentence can refer to: "to construct a statistical generator that can produce a SMB to force an ice sheet model, which should include ..."

l56 Wouldn't the sentence make more sense when 'a' and 'one' are exchanged? "We approximate the output of one process-based SMB model as a realization of a stochastic process"

l56 "The statistical model that produces a realization best fit to the process-based model output can then be used to generate hundreds of other realizations"

To get the concept clear already here: Do you produce one statistical model for each process model, or one statistical model that generalises over an ensemble of process models? Clarify. Maybe "produces a realization best fit to a process-based model"  

l58 "representing a range of possibilities for future SMB"

Most of your examples focus on the historical period and it is not clear how that extends e.g. to a future under strong warming. From looking at available RCM SMB until 2100, it seems that the amplitude of inter-annual variability increases already by 2050. How will this be incorporated into the statistical model? 

l60 "we base our method entirely on open-source software packages and provide our own open-source code where necessary." 

That is very well appreciated. Thanks again.

l62 add "surface": "stochastic surface mass balance generator"

l68 "First, we will examine ice core reconstructions of SMB in Greenland"

I am missing more detail and I am also concerned about the usefulness of ice core data for this study in the first place. Where are these ice cores located? What kind of information do they provide (before processing autocorrelations)? If they are from regions in the accumulation area, the constraining potential is probably limited to large scale snowfall. Can they be assumed representative for the SMB you are interested in, which may instead be dominated by variations in melt and coastal precipitation? I believe a critical evaluation and better documentation of this data source is important.

l75 Not sure what makes them "benchmark SMB models". Maybe remove "benchmark".

Consider to also list the model names here.

l80 "catchment outlines (Figure 1) provided by Mouginot and Rignot (2019)"

I am concerned about the choice of region delineation. See general comment.

l80 "Delaunay triangulation to produce a covering of each catchment area"

I don't understand the need and purpose to triangulate the domain. Would it not be enough to select and sum/average over the 1km grid cells that fall into a specific catchment area? 

l81 "We average SMB over the grid point closest to (or within) each triangle of a catchment".

Does each catchment have one triangle or several? Sorry I don't understand the process. Can you explain this better? 

l82 "We then average over annual time scales"

Is this a moving average or just summing monthly values to get annual SMB?

l92 Maybe you base this on formal mathematical theory, but I found calling E(t) here an "error term" confusing. Since you describe the general form of the equation, maybe something similar to l106 "spatial correlations between catchments are captured in ε(t)" would be better here. 

l96 How does m=30 map to the 33 year data set 1980 to 2012?

l97 This sentence is confusing to me. Is f(t) or b1(t) the forcing variable? 

Maybe, if that is correct: "The temporal trend μ(t) includes historical mean SMB for each catchment β0, and the forcing variable f(t) with a linear coefficient β1. The forcing variable, f(t) ..."

l105 "The error term ε(t)"

See comment l92, maybe call it something else.

l117 ARFIMA is both used to refer to the family of models (AR, ARIMA, ARFIMA) and to a specific choice. This should be resolved to avoid confusion. 

l125 I don't see what choice for f(t) is used here. Is this one of the parameters that are optimized?

l126 I am confused about the model fitting process and request some clarification. The SMB time series have a certain inter-annual variability pattern in time. Is the purpose of model fitting to fit that specific time series or to find a model that can produce a time-series with the same characteristics (maybe amplitude, frequency)? I would assume the latter. If so, how do you characterise the produced time series? Did you look at the range, amplitude distribution and spectrum? The example fits shown in Figure 3 have similar time evolution, but that may not be the criterium used for a model selection. 

l155 This seems like a surprising result and conclusion to me. How can the AR(1) process be a good choice, if the amplitude of the variability is (guessing from the figure) a factor 20 lower than expected? Wasn't the aim of finding a good model to be able to reproduce the inter-annual variability of the process model? If "residual variability may be captured in the spatial noise", why bother finding a model for inter-annual variability at all? Then the AR(1) process could maybe be replaced by a constant and all the variability captured in the spatial noise?

157 Isn't it a problem for the covariance analysis that the catchments all have different sizes and shapes and the statistical SMB is fully correlated on the catchment level? This is related to my general comment about optimal delineation.  

l165 How does the method get information about where the basins are spatially located and if they are close to each other or not? Are neighbouring basins now more likely to see the same SMB in one year? Can basins that touch along the divide be considered neighbours in this approach. 

l166 "We calculate the empirical correlation matrix Cˆ , which is an approximation of C, from the residuals of per-catchment best-fit models described in Section 3.2"

The sentence seems to suggest that each catchment has one time vector of residuals. This goes back to my question above (l126). Ultimately each statistical model can produce a number of realisations with different inter-annual variability. Wouldn't each of these produce another residual. Is the best-fit model one realisation or the model that can produce the best realisation among the infinite number of possible cases? 

 

l169 Not sure showing the code here is useful, given that everything else until now is described in mathematical notation. I would prefer to continue in that way and maybe move the code to an appendix if it is deemed important. 

l187 What is missing in the description for me here is a visualisation of the effect of the spatial correlation in the model. Could you add a figure. 

l192 I agree with this interpretation, except for the part "for ensemble simulations of ice sheet change given uncertain SMB". Uncertainty in SMB (in particular in the future) is due to many factors aside from internal variability. The method is useful to produce/introduce inter-annual variability, but needs to rely on existing SMB for a given scenario, GCM and downscaling. 

l202 A similar catchment-based elevation-SMB relationship has also been used in Goelzer et al., 2020b:  https://doi.org/10.5194/tc-14-1747-2020. Could add a reference here. 

l205 "generally quite good for Greenland SMB"

It could be good to mention that the SMB-elevation relationship is very strong in the ablation area and lower accumulation area, but mostly breaks down in the upper accumulation area. This is because melt is fundamentally controlled by temperature, while snow accumulation is also controlled by wind. This is the reason why the SMB-height relationship in ISMIP6 was based on runoff gradients instead of gradients in SMB itself.   

l208 "to capture the known feedback between ice sheet surface elevation change and SMB change"

Equation 7 (and 5) do not capture that feedback. I think in its present form it produces information only at the same elevation as the process model. It should be explained how the equation needs to be adapted to map the SMB to the evolving ice sheet geometry.   

l214 The potential problem I see with this approach in a non-stationary situation is that the number of months/years used to construct the gradients could be rather small. At least there is a conflict between producing robust gradients and capturing the transient SMB change with time. This issue may be raised here. 

l211 'SMB anomaly' is an often used term to describe the SMB difference relative to a reference climate. Here it is a difference of the local SMB to the catchment mean SMB. Suggest to make that explicit as: "the local SMB anomaly relative to the catchment mean". Same in l223

l224 'The same principle could be adapted for training data provided at even finer temporal resolution (i.e., weekly or daily).'

I think this statement is too easily made, it needs some additional work. It is not obvious in that case how to build the gradient lookups. If we go for daily time resolution, selecting all January 1st and June 23rd would probably produce very inconsistent gradient maps, because of the weather in the system. Finding the right averaging period and enough data is not easy in that case. See also comment l214.

l256 Why did you not test white noise as an option? 

L266 In addition, Antarctic SMB variability is dominated by snowfall (hardly any runoff) and therefore quite different from Greenland. If you were to fit a large central accumulation area region in Greenland instead of the elongated drainage basin regions, maybe results would be more comparable to findings for Antarctica? That is something you could easily test to maybe add some substance to this more speculative paragraph.

269 "seven SMB process models"

I am missing discussion of the differences between those seven models. Did emulating them with the statistical model produce similar results?

277 "the decadal timescale on which we focus" 

This does not become clear from reading the manuscript until here and is also in contrast to l330. If it is true, I think this sentence should be in abstract and introduction to clarify the focus of the paper. 

287 "Recent analyses also suggest that parts of the Greenland Ice Sheet may be approaching a destabilization"

I am not aware of any evidence for tipping behaviour on that time scale aside from this one paper. So I would be hesitant to propagate it. 

l315 "there is too little data available from high-fidelity models of SMB to pursue spatial methods beyond what we have shown here"

I think there is quite a range of process model output available by now and suggest to make a statement more humble to potential other research in this area. 

l325 What do you mean by "coupled simulations"? Maybe "Regional Climate Model simulations"?

Summary

 Future projections of the sea-level contributions of the Greenland and Antarctic ice sheets depend on accurate projections of the surface mass balance of these ice sheets. The atmospheric processes responsible for the surface mass balance are affected not just by long-term climate change, but also by natural variability. As both observations and model simulations of the atmosphere contain signals from both long-term trends and stochastic variability, separating the two is problematic.

The authors present a novel method for efficiently generating many realisations of the surface mass balance, using a stochastic parameterisation. The parameters of this model are fit to realisations of different process-based SMB models, so that the statistical properties of the results from the process models are reflected in the results of the stochastic model. This allows ice-sheet models to run large ensembles of simulations, exploring the phase-space of SMB variability without having to run the computationally expensive process-based SMB models.

I think this an interesting and valuable development. In general, I think the authors present their methods and results well. Being an ice-sheet modeller myself, my knowledge of statistical mathematics is insufficient to have an informed opinion about the details of the stochastic model presented here. However, I have some more general questions that I think could be answered in the manuscript.

General comments

 1) Relative contribution of uncertainty in SMB.

In ISMIP6 and related projects focusing on future projections of ice-sheet mass loss, there is a general tendency to focus either on exploring the uncertainty in the ice-sheet models themselves (poorly constrained physical quantities and processes of the ice-sheet, model numerics, resolution, etc.), or on exploring the uncertainty from the climate/SMB forcing, which is usually done by using different CMIP models as forcing. It would be interesting to see how large the contribution from variability in SMB projections to the uncertainty in sea-level projections is, compared to the uncertainties already explored in ISMIP6. I think there is already some work on this in the Verjans et al. (2022) StISSM paper, and I think it would be informative to devote a bit more text to this in the introduction.

2) Relevance for ice-sheet modellers

The stochastic models presented here are trained on different members of the GrSMBMIP model ensemble. The authors mention that an advantage of their stochastic model is its low computational cost, allowing for large ensembles of simulations. However, about half of the models in the GrSMBMIP ensemble are also relatively cheap in terms of computation time (specifically the energy balance models and the positive degree day models), to the point where they would not noticeably slow down an ice-sheet model. What exactly is the advantage of using the stochastic model over any of those “cheap” SMB models? Does it allow for a more comprehensive sampling of the SMB phase space? Does it produce a better agreement with the RCMs like RACMO or MAR (which are generally agreed to provide much more accurate output than the “cheap” SMB models)? I think the relevance of the stochastic model could be explained better.

3) Downscaling

The authors present a relatively simple downscaling scheme to convert the catchment basin-averaged SMB time-series produce by the stochastic model, to 2-D data fields that can be used in an ice-sheet model. The problem of down-scaling low-resolution (climate) model output to a higher-resolution ice-sheet model grid is not new, and also not entirely trivial. Right now, the only evaluation of the downscaling method presented here, is a visual inspection of output for a single catchment basin. The timeseries in the top panel of Fig. 6B suggest that the stochastic model greatly overestimates the variability in the accumulation at this location, at least when compared to the process-model output (is the process model here ANICE-ITM again? If so, mention this – and why you chose to compare to this one, instead of the much more accurate RACMO or MAR data). I’d like to see a bit more analysis of the performance of this downscaling method. How does the spatial variability in the resulting 2-D field compare to e.g. RACMO? How does this compare to existing downscaling methods, e.g. Noël et al. 2018?

4) Long-term trends

It is not clear to me from the manuscript, how the stochastic model deals with long-term trends in projected SMB. In Fig. 4, ten different realisations of the stochastic model are compared to the ANICE-ITM output (again, why not compare to RACMO or MAR?). That runs to 2014, whereas the stochastic models continue to 2050. Does this mean the trend in the “observed” SMB (which ANICE-ITM, and the other models in GrSMBMIP, attempt to reproduce) is extrapolated? How would that work when I want to simulate ice-sheet mass loss under different RCP scenarios? As this is ultimately the main goal of ice-sheet models, I think this deserves some attention here.

Specific comments

Line 80: “We overlay … and use Delaunay triangulation to produce a covering of each catchment area.” I’m unclear what you mean here. As far as I know, the Delaunay triangulation is a unique triangulation for a given set of points in the plane. How does this relate to the catchment basins shown here?

Figure 5: it is not immediately obvious what the small inset panels mean, until you notice that only the first one has tick values (but no axis label). Maybe consider only showing one or two months (I think that would be enough to get the point across), and promote the inset panels with the downscaling lapse-rate fit to full panels?

Lines 303 – 306: “We suggest … near-terminus SMB variability.” This is an oversimplification. Overestimating the accumulation in the interior will not immediately affect the flow of the outlet glaciers, but it will strongly reduce the projected sea-level contribution of the ice sheet (this is especially apparent in the ISMIP6 Antarctic projections, where many ensemble members show a net mass gain due to increased snowfall in the East Antarctic interior).

References

Noël et al., 2018: “A daily, 1 km resolution data set of downscaled Greenland ice sheet surface mass balance (1958–2015)”.