A report on the manuscript “FastIsostasy v1.0 - An accelerated ... by  Jan Swierczek-Jereczek, Marisa Montoya and 4 additional authors.

General Comments.

This submission to Geoscientific Model Development proposes a method and code to tackle one of the most difficult of all geodynamical problems in solid Earth geophysics: how to efficiently solve the fundamental equations of glacial isostatic adjustment with an Earth model domain that has both lithospheric thickness and an underlying Newtonian viscosity that are generally 3-D.  The problem has attracted considerable interest nowadays and has been solved in various mathematical-numerical approaches since at least the mid 1990’s.  A fundamental problem for the codes is the computational resources that are required, and the shear wall-clock time involved in obtaining and analyzing a result, and then proceeding with additional runs are formidable.

A method is proposed by the authors that has various levels of approximations.  The claim is that these yield a code that runs much faster, hence attacking, potentially, the problem of the additional parametric unknowns in the 3-D earth model setup.  This is a very worthwhile goal and it is admirable that authors set out to expand the work of Ed Bueler et al. (2007, Ann. Glaciology) which employed ‘shallow layer’ approximations to derive a code that efficiently interacts with a changing load (ice sheet) to predict the vertical deflection of the surface of the crust/lithosphere such that the latter can be used in ice sheet history/evolution modeling studies.  The latter code is for a radially stratified (1-d) Earth model.  Should the proposed (3-d) code and formulation be shown to give results that can reliably be incorporated into ice sheet models, as in the (1-d) case (Bueler), this would be a great accomplishment and would gain great traction in the ice sheet modeling community.  This sets a very high bar for the evaluation of the submitted manuscript.

Unfortunately, I have concluded that the high bar has not been achieved, and hence recommend either major revision or outright rejection of the manuscript in its current form.  Should rejection be recommended by GMD, authors are given great encouragement by this reviewer to resubmit, possibly as a theoretical treatment in Ann. Glaciology, where it might be thought of as a companion paper to the 2007 effort of Bueler.  I give my impressions of how such a manuscript could be restructured and provide a series of tests that would be more convincing than those presented.

I give a list of detailed comments, section by section.  First, however, it is important to summarize my assessment in more general terms. Here are three major areas to consider.

1. The introduction to the manuscript is poorly written. Perhaps, it’s symptomatic of a larger issue.  While Figure 1 might certainly be relevant to a follow-on paper, these frames are not relevant to this research paper.  The heart of this GMD submission is to solve a generally time-dependent 3-D viscoelastic flow problem with complex 3-D distribution of material properties, and then show that the prediction of vertical time dependent deflection is reliable, albeit with some error.  The pertinence of the problem to ice sheet stability and past and future sea-level change can be addressed in a single paragraph with appropriate referencing.  What is needed in this manuscript are more meaningful descriptions of the approximations and the tests against the Seakon code. I elaborate more in item 2 below.  The paper should give the details at the same level as might be seen in journals such as the Journal of Computational Physics (https://www.sciencedirect.com/journal/journal-of-computational-physics). In fact, such a journal might also be a logical option.
2. Fluid behavior, and elasticity for that matter, are poorly described, both from the standpoint of physics and basic mathematics. The approximation for the lithosphere is that of planar elasticity, an approximation in which a reduced stress tensorial balance is assumed.  This is never mentioned, and the severity of the approximation also goes without explanation.  Equally, for the fluid part of the physical space the assumption of the layering is that it can be ‘lumped’ and furthermore each layer acts in a manner like the shallow water approximation in physical oceanography.  This puts severe restrictions on the type of interactions a laterally heterogeneous sub-lithosphere mantle may possess (i.e., lateral flow may occur but not like in a 3-D Stokes flow that we are familiar with in ice sheet models and/or the nominal 1-D structure / 3-D GIA simulation of the poloidal flow that employs the full equation system governing the gravitational-viscoelastic response to surface forcing.  This could be accomplished by writing out the full equations, then dropping out the terms that yield the theory that is then developed in Section 2.1.
3. The referencing is poor. This is usually a problem that can easily be sorted out: referee asks for a new reference, or to delete one, and the author response is simply to agree or argue the contrary.  This is not what I mean.  The referencing as it stands in the submitted manuscript is quizzical, to the point at which it is clear, in my opinion, that the authors might not have read the papers that they so haphazardly refer to.  For example, Section 2.4 discusses the sea-level computations.  There is reference to the nonlinearity in the system that is required to be solved but has nothing to do with the traditional problem of nonlinearity treated, for example, by Spada and Stocchi (Computers & Geosciences 33, 2007, 538–562) and that is treated by Mitrovica and Milne (2003) and Kendall et al. (2005).  The latter authors also treat, in an elegant way, the problem of migrating boundaries, which, to my reckoning, is NOT a nonlinearity but non-perturbative part of the global GIA problem, i.e., the ocean boundaries need to be updated, and this is simply a trial-and-error process.   Mixing the reference to the method of Goezler (2020) and the references in the last three sentences on Page 12 is both confusing and wrong.  Another glaring example of confusion, is the fact that Table 2 shows values for Young’s modulus and Poisson’s ratio, implying that we are dealing with a compressible model.  But then the Appendix A discusses necessary corrections for the assumption of incompressiblity and claims that the Maxwell time needs to be adjusted to account for this.  The Maxwell time (incorrectly written with Young’s modulus) is not influenced by compressibility, since it is a material parameter of the constitutive equation.  It is difficult to understand what the authors are up to, but as I interpret what is written in Appendix A, there is confusion between Maxwell time and relaxation time(s), the latter of which can only be derived by solving the eigenvalue problem in the full inverse Laplace transform (though there are alternative approaches to obtaining the decay spectra).  I am afraid Appendix A is simply wrong. Some justification is needed as to why this is implemented in the author-touted Julia code.

Some additional remarks:  Some parts of the paper are well-organized, and I encourage the authors to strengthen these in a revision. The section “FastIsostasy in the model hierarchy” (which seems not to have a number ??) is well written and informative.  While the benchmark cases and comparisons in Section 4 are also good, but they are also incomplete.  I think there is an underlying contradiction, and I allude to this in the opening sentence of item 1 above.  The revised manuscript needs a table to define all acronyms (BSL, etc.).

The authors, it appears, are striving to have an efficient way to compute the coupled ice-sheet solid Earth problems.  It appears they may have made progress toward this goal.  I mention a ‘contradiction’.  It has to do with the time scales involved.  The opening discussion gives a description of an ice sheet retreating on a retrograde slope (Figure 1, which I recommend be deleted).  The time scale relevant to such analysis is comparable to the time scale of melt-water pulse 1A, or about 50-500 years, yet the benchmarks of the analysis involve time scales that are more than 1-2 orders of magnitude larger than this.  So, this needs to be clarified.  In my detailed remarks and conclusion to this review I suggest a way to better address this.  As it stands the new code is tested against an FE model on time scales of a full glacial cycle.  It’s not clear that any of the faster time-scale phenomenon has been addressed.  Again, my detailed remarks are designed to help the authors better address this.

Detailed Comments.

Abstract.

I find the phrase “solve its underlying partial differential equation” odd because Section 2.1 just jumps immediately to the approximations.

Introduction.

It appears from the 1st paragraph that the authors are interested in ice-solid Earth coupling.  In the next three paragraphs the authors are interested in Antarctica, vertical deflection of the surface, lateral heterogeneity and run time per simulation.  It’s not clear that the title fits the goals.  Secondly, it is odd that the authors choose not to reference Sasgen (2017), Kachuck (2020) or Weerdsteijn (2022), since they have goals that align tightly with goals of this paper, especially the former since he attacks lateral heterogeneity in a very simple way that is analogous to what is set up in FI3D.

In any revision, there needs to be an explanation of the disparity of modeled time scales (instability vs glacial cycle).

Figure 1 needs to be removed since this geometry and time scale does not exist in the remainder of the paper.  The physics is an important background to the motivation for the study.  This manuscript is about developing an efficient code strategy and finding how much errors exists in the simplifications needed to speed up the run times.

Sasgen, I., et al. (2017) Geophys. J. Int., 211, (3), 1534- 1553, doi:10.93/gji.ggx368.

Weerdesteijin, M.F. M., et al. (2022) Solid earth uplift due to contemporary ice melt above low-viscosity regions of the upper mantle. Geophys. Res. Lett., 49.

Kachuck, S.B., et al (2020) Geophys. Res. Lett., 47.

FastIsostasy .. (unnumbered)

This section is well organized, but here the authors need to provide a table of acronyms employed. In describing ELRA it is an oversight not to mention that using the relaxation time to estimate topography change is a method assumed by ice sheet models for a long time now (George Denton, for example).

Line 82. ‘as e.g.’ -> ‘for example as in the case of’

Line 88.  I am confused by the ‘neighboring points’.  It’s a continuous medium.

Line 93. Cite a few ice sheet modelers that use ELRA.

Line 106.  How is tau (x, y) determined? As a relaxation time it will generally always longer than the local Maxwell time, and at any x , y position there are multiple relaxation times, each with their own wave number dependencies.

Line 110 ‘vertical ..’ -> ‘gravity, vertical …’

Line 110 The phrasing ‘by means of spherical harmonics’ is comical. Better to say: ‘by solving the partial differential equations of the time dependent boundary value problem after spherical harmonic expansion of the dependent variables.’

Lines 119-121.  There is no mention of van der Wal’s 3-d Finite Volume method.  Also, perturbation methods are not really used in practice, though they can provide insight.

Line 125.  No mention is made of the ease with which a 1-D model can be structured in the context of a radially stratified earth model (just as in a seismic model) versus the great difficulty of parameterizing a Newtonian viscosity from a 3-D seismic model.

Line 141. ‘ … against analytical, 1D …’.  I don’t know of any analytical 3-D solutions.  Is this just a typo?

Section 2.1

Equations (2) and (3) need to be derived from the full governing equations in their vector and tensor forms.

Section 2.2

Line 170.  If eta sub l (x,y) is confined to a layer, then the model really is not 3-D.  In the real earth we should imagine that flow from a weak zone centered at x1, y1, z1 can flow into/out of another neighboring (weaker/stronger) zone cantered at x2, y2, z2.  This is kind of fundamental to a laterally heterogenous mantle under conditions of gravitational disequilibrium.

Line 180.  It is unclear how this equation applies to layers at different depths sine the wave number dependencies should vary as a function of depth.

Line 196.  As stated above I do not understand this scaling.  I don’t think this is rationalized correctly.

Line 215.  The reference to Farrell is incomplete.  What Green’s function?  A spherical model?  Agreement with Test 3?  That hasn’t been introduced to the reader at this point in the manuscript.

Line 221.  A Stokes flow will always have a dynamic pressure term.  The logic presented is flawed.

Section 2.3

Line 225.  We don’t normally see a partial u / partial t term in a momentum balance equation of GIA.  Please derive this equation and Eq. (14).  This is my first time to read a reference to a ‘place holder matrix’.  This perhaps this illuminates the problem with this paper for me as a reviewer.  Perhaps it needs to be submitted to a journal wherein this choice of vernacular is familiar.  J. Comp. Phys, or Physical Review, or others might be options.

Section 2.4

Comments on Figure 4 and discussion of Sea Level.  See my remarks above.

Section 3.

Rationalize the use of Julia beyond flowery language like ‘vast ecosystem’ and ‘good performance’.  There must be something at the heart of this language that has advantages over C or C++.

Section 4.1

Equation (23) has some assumption of a time history.  Please state this assumption.

Should remind the reader that wave number kappa is a cylindrical wave number that only at high value corresponds to the wave number of a Cartesian system.

Section 4.2

Figure 6.  It seems ‘sea surface’ is the depth of water column which changes as a function of time.  This needs to be explained much earlier in the paper, maybe even in the Introduction.

Section 4.3

Figure 7.  The results are really mixed here.  Sometimes FI can receive a ‘pass’ and other times a ‘fail’.

Line 457.  This is a good point about the final quasi-equilibrium snapshot being misleading.

Line 462.  Concerning shorter time steps: this is common knowledge among those of us who do these computations and formulate models.

Section 4.4

Line 486.  I am confused about how to interpret ‘mean’. “Spread around unity” leaves me in the cold.

Figure 8 is generally good.  The differences in the 3D approaches might be interpreted by some as being unacceptably large. But this depends on the specific geological or glaciological question being asked of the geoscientific modeler.

Again, there is no reference to the important time scale for instability: 50 – 500 years.

Appendix A

See above comments.

Appendix B

If material in Appendix A has any logical rationale (and I doubt it), then it needs folding into the main text.  Appendix B can be added as a Figure companion to Figure 8 in the main text.  I recommend removing the Appendices altogether.

 Suggestion for an additional test

To capture the physics of solid earth LV and stability of an ice sheet margin, I suggest the following test that samples shorter scales in both time and space.

Allow the x, y system at its various layers to have a set heterogeneity, as for example in one direction across the space.  Then consider linear (in time) surface load as a half-torus (ring) that has negative changes in mass over 500 model years, testing the model prediction of vertical deflection at 50-year intervals at the surface across the entire space.  This would excite both high and low values of kappa as the surface load would have both short and long spatial scales of loading.

In this study, the authors present FastIsostasy, a regional GIA model capturing lateral variations in lithosphere thickness and mantle viscosity, as well as gravitationally consistent sea-surface changes. The key advantage of FastIsostasy is its ability to bypass the computational expense of 3-D full self-gravitating earth models. This makes it a very promising tool, particularly for ice-sheet modellers, as it facilitates the consideration of lateral variations in solid Earth properties often overlooked in sea-level projections due to computational constraints. FastIsostasy hence holds promise for significantly enhancing the accuracy of sea-level predictions, particularly for the Antarctic ice sheet. Spatial variations, especially beneath the West Antarctic ice sheet with its low mantle viscosity, have been identified as crucial in Antarctica. The potential for triggering negative feedbacks that limit and delay mass loss adds to the interest of this novel model.

Overall, I find the manuscript to be well-written and effectively presented. FastIsostasy appears to be a compelling tool, advancing beyond previous models that addressed a similar exercise (such as the Elementary GIA model). The benchmarking in section 4 against 1D and 3G GIA solutions is very interesting and provides valuable insights into the model’s behaviour. The open-source and well-documented nature of the model adds immense value. However, the manuscript could benefit from a few clarifications:

• If my understanding is correct, FastIsostasy is essentially a 2D model due to the lumping of the depth dimension. I think this needs to be emphasised more clearly in the manuscript. This aspect is sometimes presented in a misleading way, and a more explicit acknowledgement of this approximation and its implications is needed.
• FastIsostasy allows to include gravitationally consistent sea-surface changes in a regional domain. Again, this allows for a significant improvement in ice-sheet projections, given that most ice-sheet models consider the sea surface to be uniform, missing important feedback influencing the stability of grounding lines. Unfortunately, the section introducing the sea-level model and its improvements compared to Coulon et al. (to account for time-varying ocean area) lacks clarity. A restructuring and clarification of this section, with a focus on explaining how the improvements offer a key advantage, would enhance the manuscript (see my suggestions in the specific comments below).
• Table 1 provides valuable insights for comparing FastIsostasy with other existing approaches. However, the current manuscript could benefit from a more explicit discussion in the main text regarding the similarities and differences between FastIsostasy and the Elementary GIA model. This is particularly important since both models aim to address similar challenges, namely, offering a computationally-efficient approach to incorporate LV solid Earth properties and sea-surface changes. While it is evident that FastIsostasy is more complex, there are noteworthy similarities between the two approaches. One notable advantage of FastIsostasy, briefly highlighted in the conclusion section, is its ability to avoid translating viscosities into relaxation times. Expanding on this point in the main text, in addition to the information presented in Table 1, would contribute to a more comprehensive understanding of the strengths and distinctions of FastIsostasy compared to the Elementary GIA model.

Once these clarifications are addressed, along with the specific comments provided below, I believe this paper is suitable for publication and would make a very valuable contribution to the field.

Specific comments :

• Introduction, l.24-26: I would suggest reformulating this, as it implies that enhanced melting is the driver of bedrock uplift and sea-level drop, not the grounding-line retreat itself.
• Introduction, l.25: Maybe it is worth clarifying what is meant by ‘sea-level’ here, as a panel of definitions exist. I suspect you mean the relative sea level, i.e., the difference between the bedrock and the sea surface.
• Introduction, l.31: Maybe clarify here that it is the creation of pinning points that explains the influence of GIA on the buttressing ice shelves?
• Introduction, l.34: Replace ‘parameters’ with ‘properties’? To check throughout the manuscript.
• Introduction, l.56: I would remove ‘of the ice-sheet and GIA models’ from the sentence as you do not discuss uncertainties in ice-sheet models themselves. For example, ‘Such uncertainty thus needs to be propagated…’
• Introduction, l.59-60: Coulon et al. (2021) addressed this parametric uncertainty by running large ensembles of simulations with a computationally efficient GIA model. This is worth mentioning here.
• Introduction, l.62-64: This is not entirely true, as Coulon et al. (2021) proposed a way to account for the lateral variability in Antarctic rheological properties at a low computational cost. I know that it is presented later, but I think that, even though their model is of lower complexity than FastIsostasy, it is worth mentioning their work here, and more generally in this section, especially as their motivation was very similar to the one presented here.
• Introduction, l. 65-69: Some of these concepts have not yet been introduced or only very briefly, e.g., the depth dependence of the mantle viscosity (only shown in Fig.2 but not mentioned in the text), or the dependence of the response time scale on the spatial scale of the load. Either introduce these before or simplify the description here and give these details later.
• Introduction, l.78-79: I am not sure what is meant by ‘sufficiently constrained in literature’.
• Introduction, ELRA description: I think this section lacks a few references, for example, to illustrate studies that use such a model.
• Table 1: Why does the Elementary GIA model have a ‘~’ symbol in LV? This would be worth explaining in the main text or caption. Same for the ‘sea-level’ for LV-ELRA and FastIsostasy, clarification in the caption might be useful.
• Section 2.2, 198: While I find the idea of lumping the depth-varying viscosity profile at a location into an effective viscosity interesting and understand its value for the models’ computational efficiency, I think that it remains important to discuss the limitations in more detail. In particular, it is important to acknowledge that (similar to previous LV-ELRA attempts) because you end up with a unique effective viscosity value, it does not allow to capture the full multi-normal mode response of the Earth to surface loading which is accounted for in 1-D and 3-D GIA models due to their viscoelastic layering. In fact, the larger the load, the deeper the deformation reaches into the mantle. The ease with which the mantle relaxes thus depends on the radial viscosity profile, with, e.g., the shallower layers being the more relevant at the local spatial scale.
• Section 2.3, l.209-210: Is this similar to what is performed in Bueler et al. (2007)? If yes, this should be acknowledged. If not, maybe explain how it differs.
• Section 2.3, l.215: refer to section 4.3 here for clarity.
• Section 2.4: I find section 2.4 along with Figure 4 a little confusing. In particular, I find the motivation for the extension to Goelzer et al. (2020) unclear. Overall, your regional sea-level model is largely based on Coulon et al. (2021), except that you propose an improvement to account for the time-varying ocean surface. The motivation and significance of this improvement are not clearly expressed. I would suggest that you start this section by providing more detail on what influences sea-surface changes, i.e., explaining that you calculate the regional sea-level field using equation (20). This will give the reader more context. I would then introduce the gravitationally consistent sea-surface changes, i.e., the sea-surface height perturbation, which is what has been emphasised so far in the manuscript and especially in the introduction, and which will dominate the sea-level signal. Unless my understanding is wrong, equations (17-18) are required to improve the estimation of s(t). The section would also benefit from a better introduction to s(t) and how it is defined.
• Section 2.4, l. 276, and Figure 4: I believe that SLC has not been introduced so far.
• Section 4, l.356-362: Is it really necessary to define ‘acceptable’ error bounds a priori? Unless you can infer them from actual studies using the same 3D GIA model with different viscosity fields (if so, please provide a reference), the values proposed here seem rather arbitrary. I think it is sufficient to say that larger errors comparable to parametric uncertainties are acceptable.
• Section 4.2, l.412-414: What do you think is the influence of the regional versus global domain? Could it influence the larger differences in N towards the edge of the load? It might be worth mentioning it.
• Section 4.2, l.415: I think I would use ‘comparable’ instead of ‘excellent’.
• Section 4.3, l.429: I find the name FI3D misleading. If my understanding is correct, the parameter fields in FastIsostasy are 2D and not 3D, given that the depth dimension is lumped. Please clarify.
• Section 4.3, l.440: To what do you attribute this underestimation of the forebulge?
• Section 4.3, l.453: I find the reference to ELRA misleading given that it also has a computationally-efficient LV version (Coulon et al., 2021).
• Section 4.3, l.455-459: This is an interesting point, but I do not understand the reference to Le Meur and Huybrechts here since (i) they do not only look at the final uplift map but also at the transient evolution (for the mean bedrock evolution) and (ii) they do not consider heterogeneous solid Earth configurations.
• Section 4.4, l.473: observed where? Please provide context or a reference.
• Figure A3: the colorbar is missing.
• Section 4.4, l.477-479 (Figure 8): I find this choice of the masking corresponding to the LGM extent questionable, especially as the area that matters for marine ice sheets is the area around the grounding line. I would suggest taking the whole domain, or an area larger than the ice sheet extend to include the vicinity of the grounding line. Typically, you have shown in the previous tests that FI underestimates the peripheral forebulge. The masking applied here may therefore ignore this signal.
• Figure 8: The panel subscripts are arranged in a confusing way.
• Section 4.4, l.501: I believe that t=-14kyr is not shown on the figure.
• Section 4.4, l.502: I am not sure where to look at. From Figures 8e and 8h, it seems to me that FI3D underestimates around West Antarctica and overestimates around East Antarctica.
• Section 4.4, l.507: But does FastIsostasy take into account the loading influence of sea-surface changes? I don’t think this is mentioned in the manuscript?
• Conclusion, l.525-532: I am not sure whether these limitations have their place in the conclusion, as I do not think they (or at least for some of them) have been discussed before. Perhaps instead in a ‘discussion’ or ‘limitations’ section, or in section 2, which presents the model?

Minor comments/typos:

• Introduction, l.44-48: This long sentence is a bit hard to follow, maybe try to split it for clarity?
• Introduction, l.56: ‘ensemble simulations’ -> ‘ensemble of simulations’.
• Figure 1, caption: ‘from (Whitehouse et al., 2019)’ -> from Whitehouse et al. (2019). I spotted this issue at other locations in the manuscript.
• Figure 1, caption: ‘enhanced melting at the grounding line, leading to larger thickness…’ -> ‘enhanced sub-shelf melting, leading to grounding-line retreat, and therefore larger thickness and increased outflow at the grounding line’.
• Section 2.2, l.179: ‘the following scaling factor’?
• Section 2.3, l.232: ‘a an ad-hoc’
• Section 4, l.346: ‘ice loading history’, or even instead ‘by an ice loading history over a full glacial cycle’?
• Section 4.1, l.393: ‘present the following differences’
• Section 4.3, l.421: ‘in (Gomez et al., 2018)’

Review of:

FastIsostasy v1.0 – An accelerated regional GIA model accounting for the lateral variability of the solid Earth

Jan Swierczek-Jereczek, Marisa Montoya, Konstantin Latychev, Alexander Robinson, Jorge Alvarez-Solas, and Jerry Mitrovica

Development and technical paper egusphere-2023-2869, submitted on 30 Nov 2023 to GMD

Summary:

Swierczek-Jereczek and colleagues present a regional Earth model that improves previous regional approaches by adding capabilities to account for the heterogeneous structure of the Earth as well as for barystatic and gravitational feedbacks on relative sea level. FastIsostacy is mainly build on the fast Fourier transform (FFT) formulation (Bueler et al., 2007) of the Lingle-Clark bed deformation model (Lingle and Clark, 1985), where two linear models for purely elastic and viscous displacement were superimposed. This results in a single time-dependent differential equation for the combined vertical displacement (Cathles 1975). Numerically, this equation is discretized in time using a finite difference scheme, while spatial derivatives are computed in Fourier space (Fourier spectral collocation method), allowing for using explicit integration methods instead of solving for large systems of linear equations.

The proposed Earth model also applies a regional sea level model introduced by Coulon et al., (2021), which here additionally accounts for barystatic sea-level dependent ocean area (shoreline migration), which is required in the conversion of ice volume changes into sea-level changes. The authors show that on glacial time scales this correction can explain about 4% of the sea level change.

The presented Earth model is of intermediate complexity between the simplified models used in current ice-sheet simulations and the full spherical GIA models, nicely categorized in Table 1 by the authors. FastIsostacy (or LV-ELVA in their notation) has not been coupled to an ice sheet model (at least not shown here) but from an ice modeling perspective it potentially fills a gap. The authors raise the problem, that GIA models are rarely open source and computational expensive due to the numerical integration against the global load history. Parameter ensemble simulations, as often used to estimate uncertainties in sea-level projections, require computational costs of the GIA model comparably to the ice sheet model or lower. As standalone model, FastIsostacy makes use of Julia programing language allowing for GPU parallelization for the FFT-related operations and additional functionalities. The main speed-up of FastIsostacy compared to full spherical GIA models, however, is due to the Fourier spectral collocation method, as already used in ice sheet model coupling (e.g. PISM, https://research-software-directory.org/software/pism).

General assessment:

FastIsostacy combines the previous ELVA approach (Bueler et al., 2007) with a regional sea level approach (Coulon et al., 2021) and extends those by adding a correction for ocean area and for projection (distortion factor K) and by including hydrostatic force induced by elastic displacement. The authors suggest to account for ocean and sediment loads, which has not been treated in the previous ELVA, as this may influence the boundary conditions (load perturbations should vanish at the boundaries assuming spatial periodicity) and therefore the FFT solution, even for an extended computational domain. For ice load changes in the domain center with distance to the lateral boundaries this assumption is more likely fulfilled.

The authors claim to account for a higher vertical resolution in Earth structure using a “lumped” weighted average defining an effective viscosity, a concept commonly used in terms of 2D Shallow Approximations in ice sheet modeling. I am not sure if the channelized (lateral) viscous flow of mantle material (e.g. in a low viscosity channel) can be really represented this way, as the systematic underestimation of peripheral forebulges compared to the 3D GIA model suggests. Also self-gravitational effects are ignored in the approach. The viscosity (Maxwell time) tuning, which is motovated by representing compressibility and shear modulus variations, may compensate for the effects of reduced complexity. The comparison to the 3D GIA effects also shows the effect of the missing rotational sea level feedback in FastIsostacy, which can be of the order of 10s of meters for glacial cycles simulations. This can be significant in terms of marine ice sheet instability and the onset of deglaciation in coupled ice sheet – GIA models.

Nevertheless, the modeling work done by the authors is remarkable, and add another piece to bridging between ice sheet model and GIA model community. The performed standalone test simulations cover different aspects of the implementation (verification, 1D benchmark, 3D and glacial cycle application for Antarctica), which makes sense to me. The presented misfits are below 20%, which can still be of the order of 100m. The authors argue that this is comparable to parametric uncertainties, but those misfits can be amplified in coupled simulations, especially for a more realistic deglacial ice sheet grounding line retreat (with faster time scales).

This model approach may still fill a gap, but ice sheet – GIA coupling is just improving significantly with reduced computational costs and comparably small coupling time steps (e.g. up to 10yr with a GIA model based on spherical harmonics as in Albrecht et al., preprint), while providing a more comprehensive GIA response covering all implied GRD (gravitational, rotational, deformational) feedbacks, with a globally conserved water budget, considering different sea level contributions and their interaction (e.g. Greenland) and global relative sea level fingerprints (geoid), also accounting for horizontal displacements (for validation against GNSS data). Hence, sea level projections and uncertainties could be directly linked to potential impacts along the global coastlines. This is what decision makers want to know in the end, and maybe it is time for ice sheet modelers to do this step forward, to think their models more globally.

Detailed comments:

l. 58 I think that SELEN (Spada & Melini, 2019) come with a BSD3 open source software license, which allows commercial usage. This is comparable to the permissive MIT license used for FastIsostacy, both without copyleft.

l. 162 What would be the induced effect of neglecting the distortion factor, just a rough estimate? In terms of estimating sea-level contributions, the same weighting terms should be used (Goelzer et al., 2020) in the regional sea-level model (Eq. 16). I encourage the authors to double-check for double-accounting?

l. 251 Does u_i,j in Eq. 15 also depend on time t?

l. 252 To my knowledge this requirement in the current PISM implementation has been improved by using an extended computational domain for the Earth model, at least 4 times the size of the ice domain. The authors also mention an extended domain in l. 292, but without any details.

l. 258 Maybe refer to Gregory et al., (2019) for definition of ‘barystatic sea level’ change. I think that sea level change induced by vertical land movement (V_pov) is treated separately.

l. 262 For gradual sea level changes the slightly different approach by Adhikari et al., (2020) may be a valid alternative.

l. 263 the sea level depends on the ocean area as function of the sea level, which is implicit, but not necessarily nonlinear.

l. 289 Eq. 20 is an approximation of gravitationally consistent geoid changes, allowing to approximate near-field relative sea-level changes, but the complete sea-level equation (Farrell & Clark, 1976) is not solved here, as the deformation of the whole Earth surface is not considered.

l. 295 This predefined mask should be shown somewhere. Does this also apply to the ocean mass changes in Eq. 1? Peripheral forebulge effects at the edge of the LGM ice sheet extent can provide important feedbacks and should be represented (see Albrecht et al., 2023). Does “extended domain” then mean, without the mask?

l. 306 What are FFT plans?

Fig. 5: It would be helpful to indicate in the figure where the load is located (vertical dashed line, or illustrated at top of panel a) at 1000km. Similar in Fig. 6 for the margin (10°) or shape and for Fig. 7 for the anomaly extent.

l. 440 Albrecht et al., (2023) highlight a possible forebulge feedback for ice sheet growth.

Fig. 7: It would be nice if the experiment names (150→50km, 150→250km, 21→20, 21→22) were put into a subtitle, legend or figure caption. Also the forebulge region could be highlighted (arrow?). I like the detailed error statistics.

l. 462 Doesn’t ‘stiffer’ mean a larger Maxwell or relaxation time?

l. 466 How is the forcing applied, the model uses the ICE6G ice history ΔH_ice in Antarctica and the barystatic sea level change from the northern hemisphere (or outside of the computational domain) from ICE6G as s(t)?

l. 472 This overestimation suggest that the weighted mean (Eq. 4) may not be valid for n layers.

Fig. 8: panel labels may be twisted. Panel e) and f) clearly show the order of magnitude of the rotational sea level feedback of the order of 10s of meters (as discussed in l.506ff).

l. 500 80m maximum error sounds quite a lot, even though it is below the 20% acceptance limit.

l. 501 mentions peak maximal error of FI3D at -14kyr, while the figure subtitle and panel b) says -16kyr.

l. 515 The authors may confuse model time steps with coupling time steps or coupling intervals. The more comprehensive 3D GIA models likely use time steps lower than 10 years as well.

l. 518 The factor 70000 for Seakon vs. FastIsostacy might be a bit exaggerated, as it is only valid if you compare CPUh with GPUh, while GPU also use internal parallelization. If only wall clock time or energy consumption was compared (GPU: 175W, 128 CPU core node: 280W), it would be still a factor of around 1000. So more importantly here is the fact, that computational costs for the GIA model should be comparable to the ice sheet model (or lower). And this is already the case (see Albrecht et al., 2023, Table 1, which can be associated with a factor of still around 100 for PISM-VILMA compared to FI3D).

Fig. A3 colorbar is missing

The affiliation of Konstantin Latychev states ‘Seakon’, or should it rather be the Physics Department of the University of Toronto?

Typos are covered by the other reviewers already.

References not mentioned in manuscript:

Adhikari, S., Ivins, E. R., Larour, E., Caron, L., and Seroussi, H.: A kinematic formalism for tracking ice–ocean mass exchange on the Earth's surface and estimating sea-level change, The Cryosphere, 14, 2819–2833, https://doi.org/10.5194/tc-14-2819-2020, 2020.

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

Gregory, J. M., Griffies, S. M., Hughes, C. W., Lowe, J. A., Church, J. A., Fukimori, I., Gomez, N., Kopp, R. E., Landerer, F., Cozannet, 610 G. L., Ponte, R. M., Stammer, D., Tamisiea, M. E., and van de Wal, R. S. W.: Concepts and Terminology for Sea Level: Mean, Variability and Change, Both Local and Global, Surveys in Geophysics, 40, 1251–1289, https://doi.org/10.1007/s10712-019-09525-z, 2019.

Spada, G. and Melini, D.: SELEN4 (SELEN version 4.0): a Fortran program for solving the gravitationally and topographically self-consistent sea-level equation in glacial isostatic adjustment modeling, Geosci. Model Dev., 12, 5055–5075, https://doi.org/10.5194/gmd-12-5055-2019, 2019.