Please, find my review attached as supplement

Review of The Utrecht Finite Volume Ice-Sheet Model (UFEMISM version 2.0) – part 1: description and idealised experiments

by Constantijn J. Berends, Victor Azizi, Jorge A. Bernales and Roderik S. W. van de Wal

Summary:

The manuscript by Berends et al. provides a description of the updated and new model components of the open source “Utrecht Finite Volume Ice-Sheet Model”, UFEMISM v2.0. The model is designed for long glacial-scale applications, as well es for ensembles of ice sheet projections and has been first described in Berends et al., 2021 (v1.0), in GMD. It is freely available on GitHub and provides automated unit tests and benchmark experiments using GitHub Action to test for reproducibility.

The key feature of UFEMISM, the dynamic adaptive grid, has been improved, and allows user-constrained refinement at indicated points, along line segments or within polygons. For an adaptive mesh, changing over time, second-order conservative remapping is used internally to limit information loss during frequent remapping. The same remapping tools are used automatically for input data available on different regular or irregular grids.

The authors have implemented a higher-order approximations of the Stokes stress balance (BPA) and a depth-integrated solver (DIVA), next to the previously availbale (hybrid) SIA/SSA solver. DIVA is the new default stress balance solver. In order to increase numerical stability Berends and colleagues make use of a least squares-based discretization scheme and a predictor/corrector time-stepping scheme for the ice thickness evolution. By running simplified experiments with idealized geometries according (or similar) to model intercomparison studies (ISMIP-HOM, a circular-symmetric version of MISMIP and MISMIP+) Berends et al. test for differences in accuracy compared to the published ensembles of simulations using Stokes or different approximations.

A few basal sliding schemes were added to the model code and were described in detail, but a comparison of the induced effects on ice sheets was not the focus of this study. There is a follow up study (part 2) in preparation by Bernales et al., covering the “application to the Greenland and Antarctic ice sheets”, which was not online available at the time of the review. The authors claim that the use of distributed memory architecture would enable for higher-resolution applications, but the provided scaling tests suggest the need for more improvements.

General assessment:

The paper provides a detailed overview about implemented approximations of the Stokes momentum balance, sliding laws, grounding line interpolation, time discretization of the ice thickness evolution and a generalized adaptive time stepping scheme. The theory behind has all been described elsewhere, but the authors put a lot of effort in writing out equations and discretizations. This is truly valuable, as the implementation on an unstructured triangular mesh adds quite some complexity.

This manuscript focuses on the DIVA solver (the new default stress balance solver), as the authors claim the higher-order BPA solver might be too computational expensive for the real-world long-time scale applications for given resources. I think that ISMIP-HOM is a valuable experiment, and DIVA seems to perform better than the hybrid SIA/SSA. I would encourage the authors to better motivate the acceptance level with respect to the Stokes or higher-order ensemble results (why is 20km still accepted, but 10km is not?). However, DIVA can be thought of a modified SSA approximation (considering membrane stresses in the plane), in which the effective viscosity is treated in a different way (also accounting for vertical velocity gradients). I suspect that the non-linearly diffusive component of the ice sheet flow (e.g. Bueler et al., 2007), which is relevant in many (purely shear-stress-driven) parts of the ice sheets, may be underrepresented in DIVA. Maybe, benefits of DIVA over SIA/SSA will become clearer in real world (or regional) application, as planned for the mentioned follow-up study. I am also missing in the description how (lateral/marginal) boundary conditions have been treated.

Apart from the DIVA solver, some model components (e.g. thermodynamics validated with EISMINT) have already been described in the v1.0 paper by Berends et al., 2021. It would be good to mention more clearly what has not changed since v1.0. Also in the v1.0 paper, the same MISMIP-inspired experiment has been performed, but with a SSA/SIA hybrid stress balance, flux correction, an explicit first-order finite volume upwind scheme for the ice thickness evolution and much coarser resolution. Readers may be interested in a more quantitative comparison of the new model version compared to the older version, or to other similar models (e.g. flow-line MISMIP or MISMIP2d). This hold also true for a comparison to the results of the v1.0 paper from the 10 000-year Antarctic retreat simulation with a 4 km grounding line resolution, which may be addressed in the follow-up study.

The authors also mention that they switched for the discretization from neighbor functions to a least squares-based scheme. The v1.0 paper already uses an averaged-gradient (numerical stencils) approach similar to an unweighted least-squares approach (Syrakos et al., 2017). Is this the same as in v1.0?

Berends and colleagues put a lot of effort in the updated model version and provide a detailed description of the implementation on an irregular mesh (13 pages of Appendix), and many features, which should simplify the access for new users (yet, I tried to install UFEMISM (without nix), but failed for the PETSc dependency, while searching for help in an installation manual or README). As a model description this manuscript fits well into the GMD Journal category. The six main figures and captions would benefit from some minor adjustments (see detailed comments below). Sliding laws are mentioned but not further discussed in the manuscript. The paper focuses on the DIVA solver and necessary adjustments to the time stepping, tested in some idealized benchmark tests. But from the provided idealized experiments it is hard to judge if DIVA is not missing an essential quality of ice sheet flow. The higher-order BP approximation is introduced but only tested in one of the benchmark tests. Numerical efficiency seems to be a motivation for the changed parallelization architecture and the the time discretization, but the model does not scale as well as expected. Also the authors just provide some relative speed-up factors, but not really explicit metrics. These are only minor issues which I encourage the authors to address in a revised manucript (see detailed comments, questions and recommendations below)

Detailed comments:

L10: “... irreducible uncertainty in many of these processes...”

Not sure what this means.

L43: “...have recently directed their efforts at creating new, more powerful ice-sheet models (e.g. Pattyn, 2017; Hoffman et al., 2018; Quiquet et al., 2018; Lipscomb et al., 2019; Robinson et al., 2020; Berends et al., 2022)”

I am not sure if all models on this list had the intention to become “more powerful”, as some are used for improving process understanding etc.

L47: “Version 1.0 (Berends et al., 2021) was the second ice-sheet model to use a dynamic adaptive mesh...”

There should be previous experience with AMR, e.g. dos Santos, 2019 for ISSM or Gladstone et al., 2010a, just to give two examples. Better rephrase or provide a complete list.

L62: “Part 2, which is submitted for review and publication separately (Bernales et al, in prep.)...”

It would have been great to had both as companion.

L98: “...with no significant loss of accuracy”

Hard to find quantitative numbers here or in Berends et al. (2021). If this remapping step is performed regularly also a small information loss (diffusion) can accumulate.

L135ff: What about the lateral boundaries (margins, calving fronts)?

L145: “DIVA produces velocities that agree well with the Stokes solution down to horizontal scales for basal topographical features of about 20 km”

How can “agree well” be quantified? In the mentioned references, the solution for 20km is not within the higher order or Stokes ensemble.

L164: Maybe it is worth mentioning the differences and similarity of SSA to DIVA.

L186: “but starts to deviate significantly from the Stokes solution earlier than the DIVA as the length scale decreases (Berends et al., 2022; this study).

“In the ISMIP-HOM experiment”. Generally SIA describes a (non-linear) diffusive process that is characteristic for large parts of the ice sheet flow and not represented in SSA, and likely not (or only limited) in DIVA. Please refer to the respective section in this study.

L228: “... flux condition has been replaced by a sub-grid friction scaling scheme, following the approach used in PISM (Feldmann et al., 2014),”.

Is this really a “replacement”, as both methods may have different effects in controlling grounding line flux? I think the linear GL interpolation goes back to Gladstone et al., 2010b, and then there are 2D extension (bi-linear), e.g. Feldmann et al., 2014. Would it be possible to align the mesh with the grounding line, such that meshes are either grounded or floating/icefree?

L247: “ an implicit scheme”

The is an interesting paper which may be cited as well (Bueler, 2023). Is the (semi-) implicit scheme still mass conservative (as it should be as a natural property of finite volume models)?

L284: Eq. 33: I understand the notation is consistent with Robinson et al. 2022 and Cheng et al. 2017, but τ is often associated with stresses or time scales. Better use another epsilon consistent with Eq. 34. From Eq. 33 it seems that τ has dimensions of m/s. How does this fit with the default value for ε =3 m in L287? I would have expected it to be dimensionless.

L293ff: “memory chips” and “processors”

For HPC architectures terms like “multi-core CPU nodes” are used, and the standard is rather 2x64 cores per node.

L299: “32 processors”

Do you mean CPU cores here, as mentioned in the Fig. 2 caption?

L299: “...often accounting for more than 80 % of the total computation time of a simulation.”

This very much depends on the application and the size of the computational domain. For instance, I/O can be a bottleneck for high-resolution application.

L307: “…, which could be improved by paying more attention to the way the model domain is partitioned over the processes, and the way PETSc decides which data should be communicated.”

In deed, domain decomposition as well as matrix factorization and preconditioning have quite some potential for performance speed-up. Berends et al., 2021 describes in v1.0 already a load-balanced processor domain decomposition, which would be worth the refer to.

L322: “appropriate remapping function”

You could already refer to the kind of remapping (bilinear, conservative) in Sect. ?. What input files not covering all of the computational domain, are there missing values or extrapolation applied?

L329: “full list of the 100+ fields..”

Does fields“ imply that only 2D variables are available as diagnostic output? I guess for the current size of the applications, parallel I/O is not yet considered?

L338: “The UFEMISM Github repository also features integration with the nix package manager...”

Does also imply also other options? I found an EasyBuild example in the “templates” folder in the Git repo. Is there a general installation manual or a website with instructions? What open source license is used, I found in the “CITATIONS.cff” the entry for “Apache-2.0”?

L365: “In experiment C (Fig. 4), which concerns sliding over a bed with spatially varying roughness, all three approximations result in velocities that agree well with the ensemble.”

Please define “agree well”, what is the acceptance range with respect to the FS ensemble? SIA/SSA seems to perform better than DIVA and BPA , why?

L368: “The DIVA remains accurate to spatial scales of about 20 km,...”

Gain, please define “accurate”. What deviation from the FS ensemble is accepted? Why is the 20km still accepted but 10km not?

L372: “...performed an experiment along the lines of the Marine Ice-Sheet Intercomparison Project (MISMIP; Pattyn et al., 2012). The experiment describes a circular, cone-shaped island, ...”

Please, better motivate how the experiment deviated from the original MISMIP (and from Berends et al., 2021, Sect. 3.3) I assume that it also covers two-dimensional aspects of the stress balance and geometry evolution. However, this prohibits a comparison to the semi-analytical solution provided for the flowline SSA case. The authors should indicate that they only considered the experiment with downward-sloping bed (EXP1+2, without overdeepening EXP3).

L378: “grounding-line hysteresis”

In order to not confuse the reader here, the authors should name this effect “numerical path-dependency” or similar, as this has nothing to do with the (intrinsic) hysteresis associated with EXP3 in MISMIP.

L405: “using the Schoof sliding law”

How is this choice motivated?

L425: “...in a fraction of the required computation time, or running a simulation in the same amount of time as before, but at a much higher resolution.”

I encourage the authors to provide some rough numbers, in comparison to v1.0 and maybe to other similar models? PISM with 4km resolution in Antarctica (1521 x 1521 x 221) would need for one model year about 10 wall clock minutes on one 2x64 CPU node, with 16km and 64 CPU cores (one socket) it would be about 30 wall clock hours (from Albrecht et al., in review).

L465: “The polygon-based routine can be used to increase the mesh resolution over a certain ice-sheet section, e.g. the Pine Island Glacier drainage basin. This is illustrated in Fig. A2.”

Where is the polygon-based refinement in PIG illustrated in Fig. A2?

L674: “The number of vertical layers is configurable, and is by default set to 12.”

But this would mean for 3 km thick ice more than 600m at the top and about 80m at the base? In order to resolve temperate ice at the base, this is quite coarse (Kleiner et al., 2015).

L714: “matrix whose coefficients depend on the mesh geometry and the ice velocities”

Is L^ij in Eq. F7 related to mesh geometry?

Typos and recommendations:

L8: “...man-made climate-change-caused mass loss...”, also L26

better say “human-made” or “anthropogenic”

L18: “The version control system..”

you mean git? Not mentioned before.

L20: “The i/o...”, also L55, L334

Define and better write capital I/O.

L23: “earth”

Earth

L29: comma before “range”

L56: “It also includes a version control system that includes...”

uses

L68: “ It solves an approximation of the Stokes equations...”

“different approximations”, or “for approximations”

L94: “...thus defeating the purpose of the adaptive mesh.”

Sounds a bit harsh, maybe use “offsetting the benefits of the adaptive mesh”?

L100: “… adapted into computer code...”

“the way it has been implemented”

L104: “The most complete is the Blatter-Pattyn approximation...”

Maybe mention L109 earlier, ”can all be derived by neglecting increasingly more terms in the Stokes equation”

L180: Eq. 13 and Eq. 20 use the horizontal nabla operator, should be introduced at some point, maybe in Table 1? It is defined later in Eq. E2 in the appendix.

L216: “... (which is not the square of the basal friction coefficient β, but a confusingly named separate entity, which we maintain for the sake of consistency with earlier literature) for the Weertman-type part.”

Then better use β* or some different index (e.g. T).

L264: Better not start a sentence with an abbreviation.

L319: “square grids”

Do you mean a regular “Cartesian grids”? Can the domain be also rectangular, or is a square with L_x = L_y required?

L321: “...UFEMISM will automatically detect the type of grid from the dimensions of the NetCDF file”

Does this mean, the model checks if lon/lat or x/y are used as dimensions? Or do you also use projection parameters from the metadata in the HDF5 headers, proj string?

L324: “The sparse matrices representing the remapping operators...”

This is often called remapping “weights” (e.g. YAC based on CDO).

L332: “Github Actions”

Better use “GitHub” with capital “H”, and refer to the website for the GitHub Actions software development framework (https://docs.github.com/en/actions).

L417: “ verified its performance”

What does this mean? I would associate “performance” with numerical efficiency.

L436: “define regions where the ice thickness should not change.”

For instance, along flow divides. What about defined boundary velocities?

L488: “ nV-by-nV matrix”

with V the number of vertices?

Figures

Fig. 1: Maybe provide a scale as measure, or mention the length of the domain for reference. What are the used distance measures here, mentioned later as config variables?

Fig. 2: I would assume you used the default DIVA? What kind of architecture did you use, CPU nodes with 2x32 cores? Hence, 128 cores would be associated with 2 CPU nodes using inter communication?

Figs. 3+4: In the labels, I would expect the FS/HO mean to be the line, and in transparent the ensemble range?

Fig. 5 How does this experiment differ from Berends et al., 2021, Fig10?

Fig. 6: The dashed lines in panel a are hard to distinguish, e.g. where is the red 5km solution? The solutions of MISMIP+ ice1r seem to converge for increasing resolution against a solution above the ensemble mean from Cornford et al., 2020. I am assuming the authors used the DIVA stress balance approximation here? This is interesting, as the 750m solution seems to be an exceptions of this convergence? What could be the reason for this resolution dependence. Fig. 11b in Cornford et al., 2020 suggests that higher resolution may provide solutions below the mean, same for the HO contributions in Fig. 9b in Cornford et al., 2020.

Fig. A2: I personally like the humor in this figure but 3 rows would be sufficient to show how line-based, point-based and polygon-based routines work. Also, the refinement seems to be more dense in the half-circle case compared to the circle? I guess this refers to the different config parameter named in L476? If yes, this information would be helpful in the figure caption.

References (not already mentioned in the manuscript)

Albrecht, T., Bagge, M. and Klemann, V., 2023. Feedback mechanisms controlling Antarctic glacial cycle dynamics simulated with a coupled ice sheet–solid Earth model. EGUsphere, 2023, pp.1-31. https://doi.org/10.5194/egusphere-2023-2990

Bueler, E., 2023. Performance analysis of high-resolution ice-sheet simulations. Journal of Glaciology, 69(276), pp.930-935. https://doi.org/10.1017/jog.2022.113

Bueler, E., Brown, J. and Lingle, C., 2007. Exact solutions to the thermomechanically coupled shallow-ice approximation: effective tools for verification. Journal of Glaciology, 53(182), pp.499-516. https://doi.org/10.3189/002214307783258396

Cheng, G., Lötstedt, P., and von Sydow, L.: Accurate and stable time stepping in ice sheet modeling, J. Comput. Phys., 329, 2947, https://doi.org/10.1016/j.jcp.2016.10.060, 2017.

Dos Santos, T.D., Morlighem, M., Seroussi, H., Devloo, P.R.B. and Simões, J.C., 2019. Implementation and performance of adaptive mesh refinement in the Ice Sheet System Model (ISSM v4. 14). Geoscientific Model Development, 12(1), pp.215-232. https://doi.org/10.5194/gmd-12-215-2019

Gladstone, R.M., Lee, V., Vieli, A. and Payne, A.J., 2010a. Grounding line migration in an adaptive mesh ice sheet model. Journal of Geophysical Research: Earth Surface, 115(F4). https://doi.org/10.1029/2009JF001615

Gladstone, R.M., Payne, A.J. and Cornford, S.L., 2010b. Parameterising the grounding line in flow-line ice sheet models. The Cryosphere, 4(4), pp.605-619. https://doi.org/10.5194/tc-4-605-2010

Kleiner, T., Rückamp, M., Bondzio, J.H. and Humbert, A., 2015. Enthalpy benchmark experiments for numerical ice sheet models. The Cryosphere, 9(1), pp.217-228. https://doi.org/10.5194/tc-9-217-2015