the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 20 Aug 2025
Development of a Semi-Lagrangian advection scheme in the Finite Element Model Elmer (v9.0): Application to Ice dynamics
Abstract. Transport processes are of great importance in geophysical applications, including atmospheric, oceanic, or ice flow dynamics. An Eulerian view is commonly adopted in models dealing with the mathematical representation of fluid dynamics. In such a framework, transport processes are accounted for by prescribing advection terms within the partial differential equations (PDEs) of the model. Yet, advection terms are prone to cause instabilities in the numerical solution of these equations, notably when using the Finite Element Method with a standard Galerkin approach. Various methods have been developed to overcome these instabilities, but often at the price of spurious artificial diffusion. To avoid such unwanted numerical smoothing, a commonly used technique is the Discontinuous Galerkin method, which allows for discontinuous solutions; hence, a better tracking of fine features with steep gradients without relying on artificial diffusion. In this study, we explore an alternative approach which lies in Semi-Lagrangian methods, combining elements of both Eulerian and Lagrangian approaches by updating particle positions based on the Eulerian velocity field from the previous time step. The method does not rely on any artificial diffusion and can accurately capture advection. Here, we present a computationally efficient Semi-Lagrangian algorithm to track the motion of particles in complex 3D geometries that is suitable for highly parallel computing. We illustrate the accuracy and power of the Semi-Lagrangian (SL) algorithm by comparing it to a Discontinuous Galerkin (DG) method developed within the open-source multi-physics code Elmer. We show that both the DG and SL methods provide accurate transport solutions with a different sensitivity to resolution. We conclude that, for practical use, the choice between the SL and the DG methods will depend on specific simulation requirements and the trade-off between acceptable diffusion and computational efficiency.
- Preprint
(14226 KB) - Metadata XML
- BibTeX
- EndNote
Status: open (until 03 Dec 2025)
- CEC1: 'Comment on egusphere-2025-3039', Astrid Kerkweg, 10 Sep 2025 reply
-
RC1: 'Comment on egusphere-2025-3039', Anonymous Referee #1, 30 Oct 2025
reply
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
Citation: https://doi.org/10.5194/egusphere-2025-3039-RC1 -
RC2: 'Comment on egusphere-2025-3039', Albert de Montserrat Navarro, 31 Oct 2025
reply
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.
Data sets
cmosbeux/SEMI-LAGRANGIAN-PUBLICATION: Semi-Lagrangian Advection Scheme in Elmer/Ice – Benchmark and Damage Tests Cyrille Mosbeux https://doi.org/10.5281/zenodo.15741827
Model code and software
cmosbeux/SEMI-LAGRANGIAN-PUBLICATION: Semi-Lagrangian Advection Scheme in Elmer/Ice – Benchmark and Damage Tests Cyrille Mosbeux https://doi.org/10.5281/zenodo.15741827
Video supplement
Amundsen Damage Evolution for a 50 year simulation – No Feedback Cyrille Mosbeux et al. https://github.com/cmosbeux/SEMI-LAGRANGIAN-PUBLICATION/blob/main/Videos/AMU_d1int_dt_test_2000-2050.mp4
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 1,843 | 94 | 19 | 1,956 | 22 | 22 |
- HTML: 1,843
- PDF: 94
- XML: 19
- Total: 1,956
- BibTeX: 22
- EndNote: 22
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Dear authors,
I would like to clarify the GMD requirements for publication and the content of you code and data availability section.
Note, that GMD asks that exactly that code, which was used to produce the results of the paper needs to be archived permanently (i.e., should have a DOI e.g., from zenodo).
From your code availability section I get that commit 97773a1f2 was used for the production. However, this commit is from November 2024, while the cited zenodo link dates back to 2023.
Could you please clarify this, and, if I am right, it is required, that you also create a zenodo archive for the code belonging to commit hash 97773a1f2.
Best regards, Astrid Kerkweg (GMD Executive Editor)