A stochastic parameterization of ice sheet surface mass balance for the Stochastic Ice-Sheet and Sea-Level System Model (StISSM v1.0)

Ultee, Lizz; Robel, Alexander A.; Castruccio, Stefano

doi:https://doi.org/10.5194/egusphere-2023-635

Preprints

https://doi.org/10.5194/egusphere-2023-635

Preprints

28 Apr 2023

| 28 Apr 2023

A stochastic parameterization of ice sheet surface mass balance for the Stochastic Ice-Sheet and Sea-Level System Model (StISSM v1.0)

Lizz Ultee, Alexander A. Robel, and Stefano Castruccio

Abstract. 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. For example, coastal planning literature demonstrates a demand for probabilistic sea-level projections with quantified uncertainty. Further, robust attribution of past and future ice sheet change to specific processes or forcings requires a full understanding of the space of possible ice sheet behaviors. 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. We construct a stochastic generator of Greenland Ice Sheet surface mass balance in time and space. We find that low-order autoregressive models are sufficient to accurately reproduce the interannual variability in process-model simulations of recent Greenland surface mass balance at the glacier-catchment scale. We account for spatial correlations among glacier catchments using sparse covariance techniques, and we apply an elevation-dependent downscaling to recover gridded surface mass balance fields suitable for forcing an ice sheet model while including feedbacks to ice sheet surface elevation. The efficiency gained in the stochastic method supports large ensemble simulations of ice sheet change in a new stochastic ice sheet model. We provide open source Python workflows to support use of our stochastic approach for a broad range of applications.

Received: 31 Mar 2023 – Discussion started: 28 Apr 2023

Download & links

Lizz Ultee et al.

Status: final response (author comments only)

RC1:
'Comment on egusphere-2023-635', Anonymous Referee #1, 06 Jun 2023

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"?

Citation: https://doi.org/10.5194/egusphere-2023-635-RC1
- AC3: 'Response to reviewers', Lizz Ultee, 28 Aug 2023
  
  We thank both reviewers for very constructive comments. As there was some overlap in the comments, we have addressed both point-by-point in the attached PDF.
  
  Citation: https://doi.org/10.5194/egusphere-2023-635-AC3
CEC1:
'Comment on egusphere-2023-635', Juan Antonio Añel, 16 Jun 2023

Dear authors,
Unfortunately, after checking your manuscript, it has come to our attention that it does not comply with our "Code and Data Policy".

https://www.geoscientific-model-development.net/policies/code_and_data_policy.html

You have archived your code on GitHub. However, GitHub is not a suitable repository for scientific publication. GitHub itself instructs authors to use other alternatives for long-term archival and publishing, such as Zenodo. In your manuscript you state that you will do it later; however, this is not acceptable, and our policy is crystal-clear about the fact that all the assets have to be stored in a suitable repository at the very moment of submission.
Therefore, please, publish your code in one of the appropriate repositories, and reply to this comment with the relevant information (link and DOI) as soon as possible, as it should be available for the Discussions stage. Also, please, include the relevant primary input/output data.
In this way, if you do not fix this problem in a prompt manner, we will have to reject your manuscript for publication in our journal. I should note that, actually, your manuscript should not have been accepted in Discussions, given this lack of compliance with our policy. Therefore, the current situation with your manuscript is irregular.
Also, you must include in a potentially reviewed version of your manuscript the modified 'Code and Data Availability' section, the DOI of the code (and another DOI for the dataset if necessary).

Juan A. Añel

Geosci. Model Dev. Exec. Editor

Citation: https://doi.org/10.5194/egusphere-2023-635-CEC1
- AC1:
  'Reply on CEC1', Lizz Ultee, 16 Jun 2023
  
  Dr. Añel:
  As requested, we have generated a tagged release and archived it on Zenodo. The DOI is here: https://www.geoscientific-model-development.net/policies/code_and_data_policy.html
  
  The relevant input and output data have already been stored and were available at the time of submission. Please see https://github.com/ehultee/stoch-SMB/commit/c58190d272f458b41f604e44aeb27cd086b1bb99
  
  Best,
  Lizz Ultee
  
  Citation: https://doi.org/10.5194/egusphere-2023-635-AC1
  - AC2: 'Reply on AC1', Lizz Ultee, 16 Jun 2023
    
    Aha, excuse the error in pasting DOI. The correct DOI: 10.5281/zenodo.8047501
    
    Citation: https://doi.org/10.5194/egusphere-2023-635-AC2
RC2:
'Comment on egusphere-2023-635', Tijn Berends, 31 Jul 2023

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)”.

Citation: https://doi.org/10.5194/egusphere-2023-635-RC2
- AC3: 'Response to reviewers', Lizz Ultee, 28 Aug 2023
  
  We thank both reviewers for very constructive comments. As there was some overlap in the comments, we have addressed both point-by-point in the attached PDF.
  
  Citation: https://doi.org/10.5194/egusphere-2023-635-AC3
AC3: 'Response to reviewers', Lizz Ultee, 28 Aug 2023

We thank both reviewers for very constructive comments. As there was some overlap in the comments, we have addressed both point-by-point in the attached PDF.

Citation: https://doi.org/10.5194/egusphere-2023-635-AC3

Lizz Ultee et al.

Viewed

Total article views: 493 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	BibTeX	EndNote
342	128	23	493	13	10

HTML: 342
PDF: 128
XML: 23
Total: 493
BibTeX: 13
EndNote: 10

Views and downloads (calculated since 28 Apr 2023)

Month	HTML	PDF	XML	Total
Apr 2023	67	14	3	84
May 2023	70	22	1	93
Jun 2023	65	18	8	91
Jul 2023	31	9	1	41
Aug 2023	43	21	4	68
Sep 2023	23	24	3	50
Oct 2023	35	18	3	56
Nov 2023	8	2	0	10

Cumulative views and downloads (calculated since 28 Apr 2023)

Month	HTML	PDF	XML	Total
Apr 2023	67	14	3	84
May 2023	70	22	1	93
Jun 2023	65	18	8	91
Jul 2023	31	9	1	41
Aug 2023	43	21	4	68
Sep 2023	23	24	3	50
Oct 2023	35	18	3	56
Nov 2023	8	2	0	10

Viewed (geographical distribution)

Total article views: 514 (including HTML, PDF, and XML) Thereof 514 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 29 Nov 2023

Short summary

The surface mass balance (SMB) of an ice sheet describes the net gain or loss of mass from ice sheets (such as those in Greenland and Antarctica) through interaction with the atmosphere. We developed a statistical method to generate a wide range of SMB fields that reflect the best understanding of SMB processes. Efficiently sampling the variability of SMB will help us understand sources of uncertainty in ice sheet model projections.


Total:	0
HTML:	0
PDF:	0
XML:	0