The authors present a particle-based Lagrangian scheme for advecting damage along glaciers. While the method is interesting, and has certain advantages over existing schemes, the description of how this novel scheme is coupled to existing Eulerian finite element methods could be improved. Additionally, the authors compare results to existing methods for simple cases, but this comparison is missing for application-relevant cases, making it unclear under which circumstances the presented scheme is an improvement. More specific comments, in order of occurrence:

24, "can be solved using different approaches": It might be more appropriate to refer to this as "using different reference frames", as Eulerian and Lagrangian does not refer to a specific method/approach, but merely as to whether the mesh deforms with the material (or in the case of particle-based Lagrangian methods, follows the particles) or is static as material moves through it.

29-30, "In this case, the transport equations are expressed in terms of partial differential equations (PDEs) that describe the spatial and temporal variations of the transported quantities": This is also the case for Lagrangian descriptions, and has nothing to do with the difference between Lagrangian vs Eulerian reference frames. This is the difference between continuum models (either Lagrangian or Eulerian) and particle-based methods (e.g. discrete element methods)

34-35, "Finite Element Models": The acronym FEM stands for Finite Element Methods

37-41: SUPG does not introduce a diffusion term, it modifies the test functions to prioritize information from upstream. One of the terms resulting from this corresponds to a diffusion-like term, but as it is consistently applied to all terms within the PDE, this added diffusion goes towards zero as the residual of the solution approaches zero (e.g. upon mesh and time step refinement). As such, it does not add "significant numerical diffusion", but instead adds the bare minimum amount of diffusion to prevent oscillations as a result of advective terms, with the magnitude of this diffusion-like-term scaling with the velocity and used mesh size.

53-55, "The Lagrangian approach takes a different perspective. Rather than solving equations on a fixed grid or mesh, it directly tracks the motion of individual particles throughout their trajectories defined by a flow field": This is only tangentially related to the use of a Lagrangian reference frame, it might be worth clarifying that the authors mean particle/discrete element based approaches, not the use of a Lagrangian reference frame.

56, "However, the Lagrangian approach may face challenges when dealing with large-scale simulations due to the large number of particles that need to be tracked.": Please provide a reference on this, as this does not seem considerably different from the challenges faced with Eulerian formulations where they are restricted due to the size of meshes used. Particle-based approaches are typically extremely efficient, scaling well to high-performance computing systems.

66, "https://csc.fi/elmer/": This link leads to an error 404 page. The authors might want to replace this to use a proper doi, e.g. using https://doi.org/10.5281/zenodo.7892181 (or adding a proper citation). 

Section 2.1: Could the authors clarify if the particles are initialized once at the nodes, and then tracked throughout the simulation, or if they are initialized at every time increment. Additionally, as the scheme locates the particles at the nodes, not at integration points, please clarify how the fields advected via the particles are used within the velocity updates (which requires integration-point values). Are the nodal values interpolated using DG shape functions? (and if so, wouldn't this add a considerable amount of diffusion every time increment by interpolating between the nodal values each time step?)

Figure 1: Is this figure correct in showing that the velocity and time increment are sufficiently high to advect particles over multiple elements? Wouldn't this cause instability issues related to CFL conditions within the solver for the fluid?

Figure 2 and 4: These results indicate that the presented scheme is not suitable for advecting quantities over any reasonable amount of time increments. While the DG scheme allows for time refinement to obtain more accurate solutions (which is typically expected within FEM, and follows a similar trend to fluid mechanics solvers improving in accuracy with smaller increments), it seems SL offers no possibility to get more accurate results. Could the authors add a brief discussion on which cases the behaviour of the SL scheme would be preferable over the DG scheme?

Section 4: Would the authors be able to also present results for this case using DG? As there is no "correct" solution to compare the predictions to, it would be nice to compare the two to see how much the results differ for realistic cases, especially for cases where the time increment of the flow and advection solvers match in time increment (e.g. Fig 6c).

Line 264, "small time steps (e.g. Δt ∼ 10−2 km)": I presume the authors mean 10^-2 a?

Line 267, "In this study, we do not account for feedbacks between damage and viscosity (Eq. (B4)), so the reduced update frequency of damage does not impact the flow solution.": This seems a major shortcoming of the presented scheme. Was the decision to have these fields remain uncoupled decided a priori, or did the authors try this coupling and run into issues?

Lines 331-336: Could the authors expand on this more? Typically Lagrangian formulations should be perfectly conservative. Is this related to how the Lagrangian particles are coupled and updated (in which case, please expand on how this is done, and why it makes the scheme non-conservative)?

Line 340, "and ensuring the stability of the advection with respect to Eulerian FEM implementations of SUPG methods": Eulerian methods with SUPG are perfectly stable in advective cases. Do the authors mean limiting the diffusion compared to SUPG instead?

Line 346, "increasing time steps and increasing spatial resolution": Please clarify if the authors mean "increased time step size" or "increased amount of time steps" (presumably, the former)

Line 362 onwards: Why would the SL scheme benefit from increased computational power? It seems DG would benefit from this (smaller and more time increments) whereas SL is more appropriate for when only a reduced number of time steps are taken, thus having less computational cost.

Eqs. B7 and B9: Please use \cdot to indicate dot products (as is done in B1)

The authors might want to check their reference list more carefully, e.g. on line 545 "slope limiter for p<math><mi is="true">p</mi></math>-adaptive discontinuous", line 547 where a full title is given in capitals, or the many entries where journal articles are formatted as references to book chapters 

General Comments




The manuscript discusses the Discontinuous Galerkin (DG) and Semi-Lagrangian (SL) advection methods. However, the Particles-in-Cell (PiC) approach (also known as Markers-in-Cell) is another particle-based technique that is widely applied in fields such as fluid and plasma dynamics as well as computational geodynamics. This method is known to produce substantially less numerical diffusion and can capture certain subgrid-scale features more accurately. Its implementation would likely be straightforward, as many components of the current SL scheme could be reused. I therefore recommend that the authors include a short discussion explaining why the PiC method was not considered, or outlining why it may not be suitable for their particular application.

It is not entirely clear whether the Stokes solver also employs DG or another FEM formulation. Providing a brief and explicit description of the Stokes solver in Elmer would improve the manuscript’s clarity. If DG is also used for the flow solution, the velocity field would be discontinuous across element boundaries. In this case, when combining DG with the SL method, it would be useful to clarify how this discontinuity is handled. Do the authors use the discontinuous field directly, or is the velocity averaged (e.g., at element nodes) before advection?

The manuscript refers to “linear square elements” and “linear edge elements” in the numerical experiments. Are these the same for both velocity and pressure fields? Additionally, line 329 mentions that wedge elements have 44 integration points—are these the same “linear” elements referred to earlier? It would be beneficial to provide a concise description of these elements along with the relevant Elmer implementation details.

In Section 3.2.2, the authors state that the SL scheme seldom randomly leads to either an artificial increase or decrease of the advected field accumulated over the whole domain, and this potentially being caused by the number of internal time steps or loss of particles. Figure 5 shows ± oscillations over time for all time steps used, so  assuming the number of internal time steps remains constant throughout the simulation, it is unclear why such oscillations would occur as a result of the number of internal time steps. Could the authors elaborate on this behavior? Similarly, the observed particle loss when moving across partitions requires clarification. If particles are correctly reassigned to their respective MPI ranks after each internal time step, such losses should not occur.

In Section 4, only the SL method is applied. It would strengthen the study to include a comparison using the DG method for the same experiment, thereby demonstrating the relative performance of both methods in a realistic application, rather than only in synthetic tests. This would better illustrate the practical trade-offs between the two schemes.

Specific Comments

Abbreviation: Throughout the manuscript, the unit “a” is used to denote “year.” The standard abbreviation in scientific literature is yr. I suggest the authors to consider using the latter.

Line 35: The acronym FEM conventionally refers to the Finite Element Method. Please correct this usage.

Line 96: Did the authors test a Runge–Kutta 4 scheme? Would you expect any significant improvement compared to the current approach?

Lines 113–114: Not sure I know what “traditional cell search” is. There are multiple algorithms for identifying the host element of a particle, depending on element shape and structure. Please clarify.

Line 159: The root mean square error is defined here as RMS, but referred to elsewhere as RMSE. Please use consistent terminology throughout.

Line 165: It is customary to include a brief introductory sentence before a subsection rather than leaving it blank.

Line 200: Please add a concise description of what is meant by “3D slab flow” in Section 3.2.

Section 3.2.1: Is this setup composed of a single layer of 3D square elements? Please clarify.

Line 209: Is the MPI communication of the particles done after every internal step of the time integration scheme?

Lines 210–213: This applies specifically to distributed parallelization. Are the models executed using only MPI, or a hybrid OpenMP/MPI approach? Please provide a brief comment on the parallelization strategy.

Table 1: Is the time per N approximately constant? It appears that time/N increases slightly with N. 

Section 4.1: Include a concise setup description—specifying the grid resolution and element types used would be helpful.Section 4.1 (Damage computation): Where is the damage variable computed? If it is evaluated at the integration points and then interpolated to the nodes, please clarify this in the text.

Attached is a annotated version of the manuscript with other minor comments, typos, and style suggestions.

This manuscript evaluates a semi-Lagrangian (SL) advection scheme that seems to be already developed and implemented in Elmer/Elmer-Ice. The results from the SL scheme are then compared with those from the discontinuous Galerkin (DG) method, which is also already implemented in Elmer/Elmer-Ice. The authors evaluate both approaches using two benchmark test (bodies in rotation and donut transport) and apply the SL method to simulate ice-damage transport in Amundsen Sea Sector of Antarctica. The study addresses a persistent numerical challenge in ice sheet modeling – to maintain accuracy and reduce spurious diffusion while lowering computational cost and ensuring stability – particularly when advecting sharp crevasse features.

My main concern/comment is that the amount of new contribution made in this article is not significant in terms of model development as the SL and DG methods were previously developed and implemented in Elmer/Elmer-Ice. Also, the results from the two benchmark tests visually seem to show excessive diffusion although the RMSE is less than 5%, which I found a bit confusing. Moreover, the Antarctic ice damage study is too simple as it lacks the coupling between ice flow and damage processes. Perhaps, the most significant finding of this article is that SL scheme suffers from less spurious diffusion than the DG method with larger time steps, which is useful information for the ice sheet modelers. However, it is not clear how one would estimate the optimal time step size with the SL scheme in a general continental-scale ice sheet simulation. Is there a theoretical way of estimating the optimal time step size or does one have to figure by trial and error.

Looking at the other two reviews, it is fair that the authors get a chance to respond to reviewer comments and make substantial revisions. I have a few detailed comments for the authors to consider below. Based on the authors' response regarding the paper's new contribution, I could be convinced that the article warrants publication.

– Ravindra Duddu

Detailed comments:

1. Page 2, line 45 – I would assume the SL scheme and SUPG have a long history of application in geosciences and glaciology. The limitations of SUPG can be overcome with newer methods based on variational multiscale (VMS) stabilization, which reduces cross-wind instability and robustly handles non-uniform meshes. The authors could at the least improve the literature review here and cite the following works.

Masud, A., & Khurram, R. A. (2004). A multiscale/stabilized finite element method for the advection–diffusion equation. Computer Methods in Applied Mechanics and Engineering, 193(21-22), 1997-2018.

Cheng, G., Morlighem, M., & Gudmundsson, G. H. (2024). Numerical stabilization methods for level-set-based ice front migration. Geoscientific Model Development, 17(16), 6227-6247.

1. Page 2, line 56 – Lagrangian approaches can also involve moving mesh methods (Jimenez et al., 2017), it does not have to be particle-based approaches. Semi-Lagrangian or hybrid methods also include Arbitrary Lagrangian-Eulerian (ALE) and coupled Eulerian–Lagrangian (CEL) formulations, including the material point method (MPM) that we utilized. In Huth et al. (2021a), we proposed an advection benchmark for the MISMIP+ configuration that is more relevant to ice sheet modeling than the Zalesack disc test. In Huth et al. (2021b) we show that MPM along with creep damage model can handle the feedback between ice flow and fracture, damage production and advection and ice-ocean boundary evolution. In Huth et al. (2023) we show that the creep damage model can reproduce the last two years of rift propagation in Larsen C ice shelf that led to the calving of A-68 iceberg, which is perhaps the only such prognostic modeling study in the literature. Bassis and Berg (2022) also use tracers using particles in cell method to track crevasse advection. Citing these works and putting them in context with the current work would improve the discussion of this paper and inform the reader about the pros and cons of different SL schemes.

Jiménez, S., Duddu, R., & Bassis, J. (2017). An updated-Lagrangian damage mechanics formulation for modeling the creeping flow and fracture of ice sheets. Computer Methods in Applied Mechanics and Engineering, 313, 406-432.

Huth, A., Duddu, R., & Smith, B. (2021). A generalized interpolation material point method for shallow ice shelves. 1: Shallow shelf approximation and ice thickness evolution. Journal of Advances in Modeling Earth Systems, 13(8), e2020MS002277.

Huth, A., Duddu, R., & Smith, B. (2021). A generalized interpolation material point method for shallow ice shelves. 2: Anisotropic nonlocal damage mechanics and rift propagation. Journal of Advances in Modeling Earth Systems, 13(8), e2020MS002292.

Huth, A., Duddu, R., Smith, B., & Sergienko, O. (2023). Simulating the processes controlling ice-shelf rift paths using damage mechanics. Journal of Glaciology, 69(278), 1915-1928.

Berg, B., & Bassis, J. (2022). Crevasse advection increases glacier calving. Journal of Glaciology, 68(271), 977-986.

1. Figures 1–4 can be improved along with the discussion.

In Figure 1, there is not much difference in the two subfigures except for showing two different quantities – position vector r and integral I_t. Perhaps, they can be combined. Also, the SL scheme is explained with a triangular mesh in Figure 1, whereas the mesh used for Figure 2 is noted as a linear square finite element. Was there any triangulation done to apply the SL scheme.

In Figure 2, the subfigures (a) – (g) are not labelled but are referred to in the caption. The result for the SL scheme with reinitialization looks quite diffused compared to the reference solution, however, the error is less than 5% in Figure 3 which is puzzling. Please explain this better.

Figure 3, please increase the font size of the text and axis labels and tick values.

In Section 3.2.2, the authors refer to Fig. 4a and Fig. 4b, but there is such labeling in Figure 4, only left and right.

1. Page 5 line 129 and Page 6 line 131 – I did not quite follow the argument of linear and logarithmic scaling. Isn’t linear scaling better than logarithmic. Please explain.
2. Figure 3 – Typically, a numerical method for a time-dependent advection problem is designed to yield more accurate results as the time step is reduced. Depending on the rate of decay of error with the time step size one would call it first or second order and so on. With the DG method, the error decreases with decreasing time step size, which makes sense, whereas it does not make sense that with the SL method, the error increases with decreasing time step size.  Moreover, the results in Figure 4 look basically juxtaposed. The result for DG with dt = 1 looks the same as SL with dt = 0.01. This requires more mathematical analysis and better explanation, perhaps it is explained in the literature and would be useful to discuss here.
3. Page 16, line 329 – How does it help to use 44 points as the integration of quadratic elements just needs 9 Gauss points. Over-integration sometimes causes issues including numerical instability or increasing computational cost. Perhaps, it is ok to do this with SL scheme but not with DG. Please add a comment on this.
4. In Appendix A, the authors mention parallelization and the benefits of the SL scheme implementation. Table 1 states simulations were conducted with 8 partitions. Adding strong/weak scaling studies and discussing communication overhead would improve the discussion in the paper. One issue maybe that for problems with small number of degrees of freedom, parallelization may not show optimal scaling.
5. I found more than a few typos in the paper. Also, some words have been unnecessarily capitalized, for example you can use lower case for “finite element method” “semi-Lagrangian” “discontinuous Galerkin” “streamline upwind Petrov-Galerkin” Only if it is based on the name of a person like Lagrange or Galerkin it warrants capitalization.