Summary:  This manuscript describes, implements, and tests an adaptive mesh refinement (AMR) scheme for the 2D time-dependent shallow shelf approximation (SSA) models of marine ice sheets.  The scheme is based on metric-based mesh adaptation, using fields derived a posteriori from ice thickness, ice flow speed, the magnitude of basal shear stress, or intersections thereof, as the local refinement criterion.  The essential idea is that the adapted mesh comes from solving an optimization problem for a Riemannian metric field, constrained by a fixed number of mesh nodes, using Hessians (2nd derivatives) of the current numerical solution.  The scheme has an open-source Firedrake-based Python implementation, with mesh construction from the computed metric ("mesh adaptation") via the library Mmg2d.  The AMR scheme is tested on the MISMIP+ marine ice sheet idealized experiment.  With the caveat that no MISMIP+ exact solution is known, clear reduction in mesh complexity is observed, relative to uniform mesh refinement, for a given numerical error level.  Modest speedup, in a time-to-error sense, occurs in the right circumstances, but the large cost of computing solution Hessians is noted.

Recommendation:  The ideas here, and I believe the implementation (though I have not tested it), represent valuable work.  Concretely, an open source, Firedrake-based toolchain for adding principled (metric and anisotropic) AMR to marine ice sheet simulations is a *good thing*.  The reported performance results suggest a meaningful speed-up over a uniform mesh, though apparently less than an order of magnitude.  (In this regard, speedup must at least compensate for the increased complexity.)  This work should be published in a well-written and self-contained form.  However, while the authors intend that the current manuscript is "self-contained", in fact it is fundamentally lacking in this regard.  Almost all readers will have difficulty in following the technical details; only those with prior metric-based adaptation experience could be comfortable.  Serious readers will split into those who willing use the authors' implementation, trusting that things work in some reasonable manner, and a (smaller?) group who must reconstruct for themselves, from the literature and the source code, many aspects of how things actually work, either in high-level outline or in detail.  I care about the latter group; please see the suggestions following.  I would suggest that the paper be returned for major revisions.

Suggested presentation, so as to be self-contained:  The content in sections 2.1--2.4 is reasonably self-contained, but it would be more so if it explicitly described how things work for the familiar Poisson equation, presumably with the solution (absolute value of the solution?) as the "sensor field".  A picture of how a mesh is actually generated from a metric is currently missing as far as I can tell, and it could be easily added to this Poisson case; it is the main point.  The fact that actual mesh adaptation occurs inside Mmg2d should be made clearer to the casual reader.  A clear visual image of metric and mesh anisotropy is also missing, and this could be added.  Now a new section 3 could be created containing the marine ice sheet equations (which should be stated as a coherent, essentially complete block), the new time-dependent ideas, and a much clearer discussion of the nontrivial sensor field choice for marine ice sheets.  These new things are the novel ones, relative to prior anisotropic metric-based AMR, as far as I can tell.  Such important new concepts should not be confusingly scattered in the way they currently are.  It must be clearly indicated to the glaciologically-inclined reader, presumably in the new section 3, and in any case well *before* the end of the results section, that vector-valued SSA-based marine ice sheet simulation requires a nontrivial physics choice regarding AMR, and that this choice is not a settled matter.  That is, one needs to *choose* between surface speed, ice thickness, basal stress, intersected speed/thickness, and etc as the sensor field.  (The current manuscript assumes glaciological readers will happily read all the way to page 20 and 3.4.2 before having any clear statement of the physics basis on which the mesh is refined!  This is crazy.)  The same choice should be clearly sketched in the Introduction, and perhaps even the abstract.  Then, in section 4 on results, the particular MISMIP+ experimental configuration could be described, and finally the actual results given.  Finally, many of the Figures are elaborate in the modern Nature/Science style, and thus no reader understands them; Figures 5,9,10 are this way and I suggest they be broken into discrete digestable form.  The spaghetti needs to be unwound.

Specific Comments on Manuscript

line 2,3 (in abstract):  Avoiding the boring sales pitch about "comprehensive overview of state-of-the-art" stuff.  Perhaps "In this paper we propose feature-based anisotropic mesh adaptation ..."

line 6:  Remove "more optimal"; it does not mean anything.  In the same line, "50% fewer" than what?  Say briefly.

line 11:  Almost no readers (me included) know what a "sensor field" is.  This needs a rewrite in the abstract.  I suggest avoiding "sensor field" here (but paraphrase it), and then carefully define it later, e.g. at line 181, which seems to be the next use of the phrase.

line 11:  What does "Due primarily to the iterative nature" mean here?  The reader is unlikely to believe that just because it is iterative it is therefore faster or better.  Perhaps this last sentence should be something like: "In practice the adapted mesh technique translates into a ..."

lines 14-81:  I think the Introduction is quite well-written, and it accords with my understanding of the development of models suitable for marine ice sheets.

line 22:  For the reader to understand the meaning here I suggest: "... the majority have focussed on reducing computational costs by numerical and meshing techniques, while maintaining the model physics."

line 27:  Recent significant progress has been made on this contact problem.  See:  de Diego, G. G., Farrell, P. E., & Hewitt, I. J. (2022). Numerical approximation of viscous contact problems applied to glacial sliding. Journal of Fluid Mechanics, 938, A21.

line 29:  It is not clear what "Alternative" means here.  This lead has the same meaning?:  "Mesh adaptation methods select finer ..."

line 65:  Suggest "allows to equidistribute" --> "equidistributes".

line 72:  "velocity norm" --> "speed"

lines 73--75:  I had to reread the last two sentence of this paragraph several times.  I suggest this 3 sentence replacement:  "These methods generate only a single mesh that is used throughout the simulation.  However, they do not consider the nonlinear and time-dependent coupling between the solution and the underlying mesh.  This coupling suggests an iterative mesh adaptation procedure, which we will implement."

line 76:  As noted earlier, my main concern about the current manuscript is that it is *not* a sufficiently self-contained description, given the reader audience, including myself.

line 79:  "but implement" --> "and we implement"

line 83:  Presumably \Omega has a polygonal boundary.

line 88-89:  It is true that there is not a "monotonic" relationship, but that is hardly the point.  AMR does not guarantee monotonicity here either, as far as I know.  The point is that uniform refinement can lead to only small, thus expensive, reductions in the error.

equation (1), page 4:  The displayed equation should be taken more seriously.  It should say something like "find \H_opt solving min E(\H) over \H with N_v vertices".  Then the sentence after can be shortened.

lines 96-97:  Suggest: "In anisotropic mesh adaptation the optimal mesh is only found approximately, and by iteration."

line 102:  "a way" --> "a local way"

lines 103-104:  Actually put the definition after the sentence which promises the definition: "A Riemannian metric is a smooth function yielding a positive-definite 2x2 matrix M(x) at each x in the domain."

somewhere around equation (2):  State that \lambda_1 \ge \lambda_2 somewhere, as I am pretty sure you are using that.

line 115:  "a metric tensor M is a continuous analogue of a single mesh element" makes no literal sense.  Probably: "a metric tensor field M(x) is the continous analog of the changing shape [dimensions?] of elements as one moves through the mesh."  Yes?  Otherwise "continuous" is not modifying anything?

lines 119-120:  "fully" makes no sense here except as vacuous advertising.  How about "... in R), thus the metric-based approach allows anisotropy in the resulting mesh."

caption Figure 1:  "is 16 times" --> "is in fact 16 times"

caption Figure 1:  ", which is not reflected in the figure for readability" --> "; the domain is laterally shortened for legibility."

line 121:  "metric space" actually means something, which is not here.  How about "... the Riemannian space \M = (\Omega, M(x)) models ..."  Surely one can call M(x) a metric tensor field, or just a metric field, and not conflate the ideas of "metric space" and "Riemannian manifold"?

line 126:  "is of" --> "has"?

line 127:  Since I actually know what a Riemannian manifold is, and what a metric space is, I am again irritated that you could say "Riemanninan space" or "Riemannian manifold" here without harm, but instead you have a jumble.  The metric space structure of a Riemannian space/manifold comes from minimizing the integral of the Riemannian metric over curves, right?  But you don't need that metric space structure, I am quite sure.  In conclusion, avoid "metric space" unless that is really what you mean.

following up:  I see that the fundamental confusion is already present in Alauzet 2010.  That author simply decides to measure distances between points in a "Riemannian metric space" in a way which is different from the metric space structure of a Riemannian space.  He simply decides that you integrate over Euclidean lines.  Indeed he does not need the correct metric space, and just overwrites it with his stuff using the same words.  It is sad and destructive that that author has now generated a branch of mathematics where words mean different things from the rest of math.  (You will never be able to do metric-based adaptation on a Riemannian manifold without undoing his language choice.)  Why did he call his thing "Riemannian" when he announces that it is not?  I would be cool with "Alauzet metric space".

line 135:  Just strike "which appear to be ...".  The shortening should be clearly stated in the Figure caption, e.g. as suggested above.  Then do not refight that battle; it is confusing to do so.

line 148:  ||u - u_h|| is *not* the "local interpolation error".  Correct this; I think you are merely talking about ||u-\Pi_h u||?

end of line 153:  "interpolant", right?  (I see Alauzet already makes this English mistake.)

equation (6):  How about

"find \M_{L^p} which minimizes \E_{L^p}(\M) = \int_\Omega ..."

That is, the input variable to the objective is (caligraphic) \M.  Putting the L^p subscript on "\M_{L^p}" then labels this as the minimizer of \E_{L^p}.

line 164:  What does "smooth" mean?  Drop that word?  This is an L^2 projection into a continuous FE space, right?

around lines 155-164:  I think it should be said that you don't actually have access to u and you are applying the idea to u_h a posteriori.

line 166:  At this point the reader has no idea what you will later mean by "one partition of the simulation time interval".  Just refer to a "time interval" here, and don't burden the reader with your later optimizations.

line 176: "may necessitate" --> "needs"

line 181:  At some point you need to define "sensor field".  As done here it is either too late or too awkwardly said.  Maybe connect to equation (6) at the beginning of this paragraph; say that in the scalar Poisson context one uses u for adaptation, and one calls this use the sensor field.  Now say that there are choices in the time-dependent coupled context.

line 184:  "mean" --> "method"

equation (10):  Here I think this is only well defined if lambda_1 >= lambda_2, right?  See comment about equation (2).

line 190:  "eigenvalues" --> "eigenvalue estimates"?  I don't see how they can be eigenvalues unless the two metrics are assumed to be simultaneously diagonalizable by the {e_i}?

lines 191, 192:  strike "popular" and "commonly"

Figure 2:  Note this is a detailed figure about a fiddly detail.  It should go later, as part of putting novel ideas later.  A couple of early figures emphasizing the main ideas of metric-based anisotropic mesh adaptation should be substituted.

section 2.4:  This immediate descent into details of postprocessing misses the main message: from the a posteriori computed metric estimate M(x), or M^i(x) in the time-dependent case, you are asking Mmg2d to generate mesh(es).  *Then* talk post processing details, though probably put later.

line 238:  Speaking as a reader, I would want the "several shortcomings" to be explained here.  Evidently I am too stupid to remember the brilliant exposition from the Introduction.

line 245:  Interleaving the explanation of the algorithm and advertising its benefits is not helpful.  Do the former first.

line 247:  Probably clearer without "before they are adapted".

line 248:  Again, advertising is mixed into algorithm description.  Not helpful.

line 252:  On first pass I definitely did not understand the meaning of "The fixed-point algorithm iterates ..." *here*.  (I wrote: """Is this recomputing solutions u along the way?  Are you returning to the start of the time partitions to revise the metric, or is the metric getting updated only in a feed-forward manner.  Iteration is an ambiguous word, but your use of "fixed-point" suggests things are not feed-forward only.  Does the "x" symbol in Figure 2 mean "compute a new metric but don't change the mesh"?""")

section 2.5 generally:  The first 4 paragraphs are quite good.  Starting in the 5th paragraph the job is to describe unambiguously an algorithm.  A pseudocode may be better than the paragraph-with-interleaved-commentary-and-advertising current form.  Figures showing the action on coarse meshes might help.

generally:  These time-dependent mesh adaptation schemes seem to be awkward, and early-development.  Do you really recommend them for production simulations?

section 2.6:  This part makes sense to me, but how understandable is it if you are not already comfortable with supermesh ideas, interpolation, and Galerkin projection?

line 290:  The avoidance of negative values via P^1 makes sense, but I suppose this begs the question of whether one should, instead of h-refinement and the costs of AMR, instead p-refine in a way which avoids the out-of-bounds problem.  See recent work:  Kirby & Shapero 2023.  High-order bounds-satisfying approximation of partial differential equations via finite element variational inequalities arXiv:2311.05880

line 305:  Give some citation for "PETSc's Riemannian metric functionality"?

line 308:  As a nonexpert I don't know what "mesh optimization" means here.  Which of the many ideas relating to building better meshes does it refer to?

sections 2.6, 2.7:  Easy to read.

line 340:  "the equations of glacier flow" should actually be stated.  (The "the" is ridiculous here; it is a field with lots of models.)

line 342:  Indeed this alternation is "as with most current ... models".  The scheme is called the explicit Euler approximation to the coupled model, that is, the scheme is explicit and first-order!  As with current ice sheet models, attaching the "prognostic" and "diagnostic" language, a leftover of pre-numeric weather forecasting, is not helpful.  The current authors are not to blame here, but they don't have to climb into the dumpster too.

line 353:  Arg.  Is "P" the pressure or the effective pressure?  Effective pressure is the *difference* between the ice pressure (overburden pressure) and the pressure of subglacial water.  At the bottom of floating ice this value is zero; assuming h is ice thickness  and d_W is the depth below the water line of the bottom of the ice, then rho_i g h = rho_W g d_W because of flotation.  Please don't use "P" for the effective pressure.  Fix all this.

line 356:  Please do not assert that the scheme in IcePack is unconditionally stable.  It is not.  Importantly, Shapero (2021) does not assert that it is.  They use a scheme which is unconditionally stable *for the linear advection equation*.  Then they say: "We note that, while the stability properties of different schemes for the linear advection equation guided our choice of method, the coupled system for both thickness and velocity is not linear and not hyperbolic."  The scheme in IcePack is implicit for the linear advection equation but *not* implicit for the coupled system.  In this context it completely suffices to say that "Finally, IcePack implements a second-order Lax-Wendroff scheme for transport (Shapero, 2021)."

Figure 4:  If this result is worth showing for uniform refinement, then also show it for a preferred AMR strategy?  Presumably, with fewer dofs than the 0.25km uniform?

line 363-364:  Me too.  You don't need to include this detail particularly, but the cp linesearch in PETSc never seems to work for me, and I fall back to backtracking.

line 390:  "emulate a more general problem".  Huh?  This sentence probably means "In the mesh adaptation experiments to follow, neither exact nor high-fidelity solutions are assumed to be available."  On the other hand, a couple of sentences later you treat the 0.25km uniform result as high-fidelity.  Rewrite this paragraph?

line 396:  I don't know what "effectively ... unbounded" means, just after you give bounds.  Frankly, I have a sense that something is numerically unstable etc., and you are special pleading.

section 3.3 generally:  I think you can write up your experiment choices more declaratively and less defensively.  To first approximation, just say what you did.

Figure 5:  A busy figure trying to show several different things, and with different horizontal axes.  (Likewise Fig 9, Fig 10.)  Decide to show one thing and show it?

section 3.4:  Use "H" for ice thickness.  "h" already has its usual meshing meaning of mesh size.

line 431:  This casual mention of using \tau_b for mesh adaptation seems to be the first contact with glaciological interests.  On page 17.

subsection 3.4.1:  Speaking as a reader who is making an effort, I am a bit lost at this point on what the difference between global and hybrid is.  I see hybrid is better (= closer to uniform ref).  Is hybrid more expensive?  Generally much of this subsection reads as chain of thoughts of the numerical experimenter, not a purposeful statement of results.

line 442:  "equal observations" means what?  Re-word.

subsection 3.4.2:  I have already made the point that this physics choice, of which field to adapt upon, should be put much earlier.  Needless to say I am not trying to make this point again.

equations (13), (14):  Also an equation for \tau_b, right?  (If not I am definitely missing something ... and it is your fault.)

line 500-501:  Since basal stress nulls as you cross the grounding line, perhaps intersecting tau_b and the strain rate magnitude as sensor fields?  The latter is telling you about the parts of ice shelves which are interesting.  Maybe the Hessian of speed is all you need, not strain rate magnitude (i.e. as sensor field); so intersect tau_b and |u|?  The thickness will follow?  In any case, the equations of a glacier are (at core) balancing stresses and strain rates to produce smooth velocity and thickness fields.

Figure 8:  This dark dense figure tells me little and makes me want to look elsewhere.  Maybe just reverse the shading.

line 503:  Again, "not reflected in Fig 7 to not hinder readability" ... is hindering readability.  Since this can be clarified already for Figure 1, I don't think you need to say anything here.

Figure 9 caption: State here that errors e_h, e_u are relative to 0.25km uniform ref.  The reader jumping to this Figure needs that knowledge first.

lines 534-537:  I appreciate this honest assessment of the technique.

line 561:  "utilized libraries" --> "libraries used here"; redundant statement of Firedrake usage

lines 566-568:  Mention this fact, that basal shear stress is the preferred sensor field for efficiency (error vs time) in abstract?  It is actually a result, not just a "we did some simulations" statement.

line 570:  "permille" is not widely used in English literature?  But perfectly clear.  Interesting.

lines 573-574:  suggest replacing "primarily due to the iterature procedure of mesh adaptation algorithms" with "reflecting the cost of our iterative mesh adaptation procedure"

line 575:  "utilising" --> "using"

line 576:  "utilised" --> "used"; suggest replacing "features of metric-based adaptation" with "available metric-base adaptation techniques"

line 582:  "utilised" --> "used"  ... and etc ... you get the point

line 589:  "necessitated" --> "required"

lines 597-599:  My conclusion has been the opposite, and it is driven by experience switching from idealized experiements with smooth bed topography to modeling whole ice sheets in the BedMachine era.  In the later case one has bumpy bed essentially everywhere, and bed topography always locally influences basal shear stress through thermodynamic and liquid water means.  Thus everywhere needs about comparable resolution.  For example, as we get good bed elevation observations in East Antarctica, then (at least in a 2D sense) everywhere will need good res., with only a slightly enhanced need for resolution near e.g. grounding lines and shear margins.  Said a different way, mesh adaptation is effective when part of the domain is easy, and I don't think that is meaningfully true, except for the small area of ice shelves, for real marine ice sheets with real beds.  The contrast with mesh adaptation for fluid problems within engineered geometries is strong.

line 600:  My favorite line of the paper.

Dundovic and co-authors present a new ice flow model that represents significant progress in adaptive mesh methods applied to ice flow problems. The key novelty is in an iterative adaption technique that does not simply produce a time-sequence of meshes, each based on the ice sheet state at a given time, revises the mesh at each time step to explicitly obtain optimal solutions at later times. The paper also presents some generally applicable ideas in that type of mesh refinement in plain language (with some illustrative math) and I enjoyed reading those parts in particular. Overall a good paper that should be published and fits GMD well. My comments are minor.

General Comments

The paper claims at a couple of points that adaptive mesh (AMR) methods are underused in ice flow modelling and should be used more on the basis of this paper. These are not justified (or technical) points and should be removed (e.g last paragraph of the conclusion). There are several adaptive mesh ice flow models: Ua (often with AMR), BISICLES (always), ISSM (sometimes), UFEMISM, in use on large scale, realistic problems today (the paper references these codes). The model here is tested only on a toy problem with a tiny domain (MISMP+), and even then, halves the domain. There are also well-regarded uniform resolution models (PISM, CISM and others). Clearly, many modellers find AMR not utterly compelling even when it is available.  This paper will not change that, because the performance improvements noted relative to fixed-in-time meshes are modest.

There is no mention of the expected or observed asymptotic rates of error estimate convergence with 1/N. Figure 9 shows error estimates, which look to me to be O(N^-1) in panel (a) and perhaps tending to O(N^-2) in panel (b). Figure 9 could be modified to show those reference rates, and some discussion made – what is leading order error do you expect from your discretization? The discontinuity Tb (or rather, in dTb/dx in this case, I think) likely reduces you to first order in space. It also sounds from your description (compute the velocity, then advance the thickness) as though you are first order in time.

Specific Comments

L84 H of \Omega -> H(\Omega) ?

L119 \rho and r_1 are adequate to describe the mesh parameters, but cannot ‘fully support anisotropic mesh adaption’ – for that you need the other things you describe (sensor fields etc). Removing the sentence would retain the meaning of the paragraph.

Figure 1. I would like to see a zoom of the GL regions – that way it is easier for the reader to see what is meant by the signs of r_1 & r_2

L126: given that the unit size is impossible, why not say `approximately unit size’?

L134-135. I see what you are saying there, but I think you can say more -in some regions the elements are stretched along x, in other along y.

L147: both h- and p-adaption are ‘active and impactful’. Agreed, but citation needed. Also, p-adaption is in many problems the more efficient approach, you might comment on that.

L170 : say a bit more about Cst. To me it seems to be space complexity, controlling the mesh density in a ‘uniform’ fashion in space (crudely, you can obtain a mesh with 2X Cst by splitting each element in two, I know that is not what you actually do). But how is it a time-complexity?

Fig 2 and section 2.5: I did not find Figure 2 helpful – only partly understanding it by reading the (quite wordy) text of section 2.5, rather than the figure illustrating the text. I would like to see these algorithms presented as pseudo-code, which would make section 2.5 more accessible.

Section 2.6 Interpolation between meshes. Indeed, this is important, at least in principle. It is for the sake of conservation and avoidance of (a type of) numerical diffusion when transferring thickness/temperature between meshes that we use the more restrictive block-structured meshes in BISICLES (which relies on a long-standing AMR library, Chombo, originally developed for shock problems in hyperbolic PDEs, where numerical diffusion is anathema). You might note this.

L318-319 Not quite. Viscous (longitudinal and lateral) stresses matter in at least some parts of ice streams, close to the GL, or where basal traction is low, or in shear margins.

L344 – this is what makes me think your methods is first order in time. Is that the case?

L349 / eqn 12. Is this really done just to smooth the GL? In which case, you are following the approach of Leguy 2014 as well as AMR. Or is it in fact because eq 12 thought to better represent the physics than the discontinuous laws? The pure ‘Weertman’ problem (P -> infinity) is harder (Tb is discontinuous across the GL), than the modified form (Tb is continuous across the GL, at least when P is given by the formula in L354* but dTb/dx is not), and the best justification for not addressing it is physics. Also, should \beta^2 be \beta^{-2} in eqn 12. So that Tb = \beta^2 |u|^{1/3}. * You say that L354is the formula for water pressure under floating ice, but you only need such a formula in the basal traction under grounded ice. Or do you compute P there differently?

L356 ‘second order’ – but only if you supply it with velocity at t + dt/2, rather than t? Seems not to agree with L344.

Fig. 4 / sec 3.2 I am not convinced that you are seeing convergence here. As you note earlier, you need to show a pattern of successively closer solutions, which might apply to the 1km-500m-250 sequence, but you need one more refinement (125m) to look convincing. The gap between 1km and 500 m looks much larger (10x?) than the gap between 500 and 250 m, but once the leading order term in truncation error dominates, you hope to see the gap successively halving (for O(N^-1)) / quartering (O(N^-2)) etc.

L368 ‘The Ice1 experiment starts from *an* initial state’. Just a typo

L382 – why this time-step in particular? CFL on your finest mesh?

L389 – all simulations were run in serial – comment on this, was that just to measure CPU time, or is your method not yet suitable for parallel runs.

Fig 5 captions – define the symbols used on the axes e.g deviation in V.A.F ($||\Delta V_f ||_\infty$)

Fig 6 (and 7) Zoom in on the GL – as it is I can see there is some difference in your meshes, but not in detail

L 487: ‘Ice bends etc’ – in SSA? There is no adjusting to buoyant forces in SSA, which assumes hydrostatic balance. The horizontal stretching occurs because all the driving stress from the front / local slopes is resisted by viscous stresses and hence results in horizontal longintudinal strain-rates.

L496 ‘basal stress is discontinuous across the GL’ – not according to eqn 12. P -> 0 as you approach the GL from the upstream side (or is P computed differently from the formula in L354?).  Also L500: as u->0, P -> \infty eqn 12 becomes| Tb| = | \beta^2 u^1/3| - the discontinuity you are seeing at the inflow boundary seems like an arithmetic artefact, with u/|u| being undefined (but also irrelevant since |Tb| = 0)

L575: This paragraph seems spurious.

L597: I disagree. You have demonstrated progress, but it is up to the reader to decide whether they should adopt the methods.