the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 23 Jan 2026
Lagrangian tracking methods applied to free surface boundaries in numerical geodynamic models
Abstract. A desirable characteristic of mantle convection models is the ability to determine surface topography evolution over time on a global scale. This capability is increasingly important due to growing interest in coupling planetary geodynamics with climate, landscape, habitats, and biological evolution systems. A common way of achieving this in numerical geodynamic models with a fixed Eulerian grid is through the implementation of a free-surface boundary condition using the sticky air method in conjunction with an appropriate method of tracking the free surface. Although existing methods for tracking the interface between the air and rock layers are available, they often struggle to provide high-resolution results on a global scale without incurring significant computational costs. We propose a method for representing surfaces directly using Lagrangian markers that can track the location of a free surface in 2D and 3D models and test it using the finite-volume mantle convection code StagYY. This approach offers a direct, high-resolution representation of the surface without the need for a large number of high-cost tracers throughout the model domain. Benchmarks demonstrate the effectiveness of this method compared to the commonly-used marker-in-cell method. The direct representation of the surface enables additional features such as the direct tracking of sea levels over time, potential for coupling with surface process models on the global scale, and enables the implementation of alternative discretisations of the Stokes equations to improve Stokes solver accuracy near the free surface boundary.
- Preprint
(4659 KB) - Metadata XML
- BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2025-6546', Anonymous Referee #1, 08 Mar 2026
-
AC1: 'Reply on RC1', Timothy Gray, 18 May 2026
We thank the reviewer for their positive assessment, interest in our work, and endorsement of our benchmarks. We have extensively revised and restructured the paper in response to the reviewer's comments, and in particular have shortened it considerably from 49 pages to 30 pages by removing unnecessary details and redundant results.
We appreciate the constructive line-by-line feedback, which we address below.
Comment 1: Terminology Confusion ("tracers" vs "markers")
"At line 35 we read: 'For clarity, the Lagrangian points already implemented in StagYY and distributed throughout the model domain are hereafter referred to as "tracers".' I find this choice rather unfortunate: the standard method of tracking materials/fluids inside the entire domain in these codes is called the marker-in-cell technique (e.g. Gerya & Yuen, 2003, 2007) but in this work the authors choose to call such Lagrangian points 'tracers' while referring dozens of times in the text to the marker-in-cell method... As such I suggest a change of names: markers for 'volume' points and tracers for additional surface points."Response: We completely agree with the reviewer that our original choice of terminology caused unnecessary confusion, given the established naming conventions of the marker-in-cell technique. We have updated the text throughout the entire manuscript to follow the reviewer’s suggestion. The standard full-domain Lagrangian points are now consistently referred to as "volumetric markers". The specialized points introduced for the free surface are now referred to as "surface markers".
Comment 2: l83: calculating -> calculated
Response: Section removedComment 3: l88: C -> C's
Response: Section removedComment 4: Section 2.2: Concern about material entrainment
"Although the topic is addressed later with the surface markers techniques in section 4.2, I was wondering when reading this section about entrainment?"
Response: Section modified, material entrainment is addressed laterComment 5: l111: Definition of 'surface mesh refinement' in FDM/FVM codes
"What is meant by 'surface mesh refinement' in the context of FDM codes here?"Response: Section removed, but the refinement in question was referring to decreased vertical grid spacing near the surface.
Comment 6: l130: Reference for the 'divergence-free quadratic spline method'
"I think a reference for the 'divergence-free quadratic spline method' in this context would be welcome."Response: We have added the appropriate references for this method, citing Gerya et al. (2021).
Comment 7: l157: Specify page/chapter for Gerya's textbook
Response: We have updated this citation to specify Chapter 18 of "Introduction to Numerical Geodynamic Modelling" (Gerya, 2019).Comment 8: l169, 257, 265, 316, 321: Equation formatting and sentence structure
"This is not a full sentence, I think a comma should be inserted in Eq 2, and the line should start with 'where ...'. Same at subsequent lines."Response: We have systematically updated the punctuation of all equations highlighted by the reviewer (as well as a general sweep of the rest of the manuscript). Mathematical equations are now treated as integral parts of the text, with correct trailing commas or periods, and subsequent explanatory lines properly starting with a lowercase "where".
Comment 9: l224: Reference for the Yin-Yang grid
Response: We have added a reference to Tackley (2008), which presents the explicit implementation and geometry of the Yin-Yang grid layout within StagYY.Comment 10: l255: 'is' -> in
Response: Section removedComment 11: l376-380: Paragraph on numerical diffusion
"I am not sure of the point of the paragraph. Diffusion is mentioned, pros and cons are mentioned, but it is not clear to me whether this is implemented/used."Response: We have revised this paragraph to make its operational purpose explicitly clear. We now clarify that global reinitialisation inherently introduces a small amount of numerical diffusion, which functions as an implicit smoothing mechanism. We state clearly that this is actively used and implemented automatically within our bilinear method to maintain surface stability, and contrast it with the integral method which lacks this implicit diffusion and therefore requires explicit slope-limiting/smoothing.
Comment 12: l542: Semi-colon replacement
"The semi-colon should probably replaced by a comma + i.e. ?"Response: Corrected as suggested to improve readability.
Comment 13: l670: Word duplication
"The word marker appears twice."Response: We fixed this typographical error.
Comment 14: Section 7.3.1: Citations for geodynamical codes coupled to surface processes
"I think that some key publications of geodynamical codes coupled to surface processes could be cited in the context of this paragraph."Response: We agree that expanding this context strengthens the section. We have incorporated key community references demonstrating the direct coupling of deep mantle/crustal geodynamics with surface process models (e.g., Salles et al., 2018; Yuan et al., 2019; Bahadori et al., 2022).
Comment 15: Reference to the DOUAR code
"Somewhere in the text a reference to DOUAR (Braun et al, PEPI 171, 2008) should be added since the free surface was also tracked by means of a Delaunay-triangulated surface embedded in a Eulerian mesh."Response: We have added an explicit discussion of DOUAR (Braun et al., 2008) in our introduction of alternative techniques, noting that our 3D integral method builds upon a similar conceptual framework of a Delaunay-triangulated surface embedded within a fixed grid environment.
Sincerely,
Timothy Stephen Gray, Paul James Tackley, and Taras GeryaCitation: https://doi.org/10.5194/egusphere-2025-6546-AC1
-
AC1: 'Reply on RC1', Timothy Gray, 18 May 2026
-
RC2: 'Comment on egusphere-2025-6546', Anonymous Referee #2, 13 Mar 2026
General comment:
This paper describes a method for representing free surfaces in mantle convection models and its implementation in the code StagYY. I appreciate that the method is thoroughly tested and the implementation seems carefully done. The new implementation appears to be more computational efficient which is of clear value. I have two major concerns that I think needs to be addressed before this paper is considered for publication:
1) The paper is very long (48 pages, 37 figures) and unstructured. It is therefore hard to get an overview of the paper. My suggestion is to
a) shorten the paper to <30 pages. This can partly be achieved by moving some of the technical details, basic equations and figures to appendix, and trimming the text. b) Work on the structure.The current structure now is:
1. Introduction2. "Limitations of the marker-in-cell method". This still is written like an introduction but there are numerical results mixed in
3 "Surface tracking" - it is not clear from the start of this subsection that this is the novelty. Thes subsections here have no hierachy.
4 "Stabilisation" -seems like this should be a subsection to surface tracking?
5. "Coupling". Write what you mean by coupling already in the title. Coupling between... ?
6. "Results" - Usually results would come after a section explaining what the experiment is.
7. "Summary", but it is actually summary, discussion and outlook. This is way to long and it is hard to get an overview.
Overall the subsection-structure must have a clearer hierachy and more informative headings. Numerical results are scattered throughout the paper which is unorthodox but could perhaps work but only if it is clearer from the beginning what the core numerical experiments will be and what experiments are used to explain the topic.
2) The level of novelty is not clear. How much of the novelty lies in the methodology itself, and how much lies only in its implementation in this particular code and the thorough testing? This must be clarified, and more references to similar earlier work is needed (the reference list is relatively short). What is special with the implemention in this particular code?
Specific comments.Line 24, explain what sticky air is to the reader
Line 74; "the bilinear method and the
integral method, which differ in how they transfer surface information to the underlying grid" are these new methods?The introduction is lacking an overview of similar methods introduced in other codes.
Why is a height function formulation with simple vertical grid deformation not used in this field?
Line 90: Clarify that you mean resolution in the vertical, not 100 m resolution in the horizontal direction: "While this method enables vertical sub-grid level topographic resolution, doing so accurately requires large numbers of
tracers per cell. For example, providing 100 m topography sensitivity on a grid with a vertical resolution of 10 km would require at least 100 tracers per cell."Figure 1: It would help if you rotate the figure so that the air is on top of the rock
Figure 3: It looks like the surface topography it self is very different, not just that it is noisy, and it does not seem to converge?. Can you commment on that?
Figure 4 and 5 are very large can be put together in one figure.
Line 140: is this your method? What similar things has been done in other codes before?
Line 140: Some general introduction text is needed here before all the subsections.
Line 165 - Since this is a method for which you have a reference, you can move some of the details to an appendix.
Equation 5- 6: Is it necesssary to write out these formulas? Can they be moved to appendix? Many of the following equations in section 3 are also standard can can be moved to appendix.
Line 381: Add some reference to slope limiting
Line 392 - The name of this subsection is very general
Line 420 - Explain exactly what you mean by coupling.Figure 15 - Should this be in the result section?
Figure 22: Change "topography with time" to "maximum topography with time"
Line 503: Is section 6.1 necessary? it is very short. Same with section 6.5 If you want to have summaries/overview, I think they should contain a reasoning for why you do all these tests.
Section 6.2. How do you choose the number of tracers to make a fair comparison? In general I would appreciate if you clarify how you ensure your experiments constitute a fair comparison.
Line 533 If the discrepency could be due to the sticky air, what is the value of this benchmark? Clarify
Figure 24: Topography -> maximum topography?
Figure 25: This figure takes too much space
Figure 30,32,34 Also takes too much space
The reference list is quite short and with many papers written be the same authors.
Citation: https://doi.org/10.5194/egusphere-2025-6546-RC2 -
AC2: 'Reply on RC2', Timothy Gray, 18 May 2026
Thank you for taking the time to review this paper. We have substantially restructured and revised the paper in response to the reviews, and have addressed your specific concerns as follows:
Major Concerns:
---------------
Comment 1: Paper Length, Visual Bloat, and Structure
"The paper is very long (48 pages, 37 figures) and unstructured. It is therefore hard to get an overview of the paper. My suggestion is to: a) shorten the paper to <30 pages... by moving technical details to appendix, and trimming text. b) Work on the structure."Response: We have significantly changed the paper's structure to trim it down to a more suitable size. Unnecessary introductory content has been trimmed, the original section 2 has been mostly removed, and standard mathematical derivations have been condensed or removed, leaving only the key content in the body.
Comment 2: Clarification of Novelty vs. Implementation
"The level of novelty is not clear. How much of the novelty lies in the methodology itself, and how much lies only in its implementation in this particular code and the thorough testing? This must be clarified, and more references to similar earlier work is needed..."Response: We have revised the Introduction and Discussion to explicitly demarcate our conceptual contributions from code-specific mechanics:
- Methodological Novelty: The novelty lies in the systematic development and rigorous comparative benchmarking of two distinct surface-to-grid mapping topologies (the node-mapped bilinear method vs. the direct unstructured Delaunay triangulation/integral method) designed specifically for preserving a sharp, sub-grid interface within parallelised Eulerian codes experiencing severe geodynamic deformation (e.g., localized subduction over geological timescales).
- Implementation Focus: We clarify how this is mapped efficiently onto a staggered-grid, parallelised finite-volume framework using MPI ghost/overlap regions.
- Literature Expansion: We have expanded our references to ground this work within a broader community context, referencing alternative formulations in other codes (such as LaMEM, DOUAR, and standard level-set/VOF models), highlighting exactly what distinct advantages our decoupled surface marker framework introduces.Specific Comments:
------------------
Comment 3: Line 24: Explain "sticky air"
Response: We have added a clear definition in the introduction: "sticky air" is a numerical approximation consisting of a low-viscosity, low-density layer placed above the solid planet boundary to simulate a true free surface condition on a rigid Eulerian grid.Comment 4: Line 74: Novelty of the Bilinear and Integral methods
"“the bilinear method and the integral method, which differ in how they transfer surface information to the underlying grid” are these new methods?"Response: Disussion of the novelty of these methods has been revised.
Comment 5: Overview of alternative approaches (e.g., Height Functions, Grid Deformation)
"Why is a height function formulation with simple vertical grid deformation not used in this field?"Response: Such approaches have been used, e.g. by Duretz et al. 2016 and Kaus et al. 2016 - references have been added accordingly.
Comment 6: Line 90: Clarify vertical sub-grid resolution
"Clarify that you mean resolution in the vertical, not 100 m resolution in the horizontal direction..."Response: We have rephrased this text to explicitly specify "vertical sub-grid level topographic resolution" to avoid any horizontal ambiguity.
Comment 7: Figure 1: Visual rotation
Response: This figure has been removed.Comment 8: Figure 3: Discrepancy, noise, and convergence in the marker-in-cell method
"It looks like the surface topography itself is very different, not just that it is noisy, and it does not seem to converge?. Can you comment on that?"Response: The marker-in-cell method is inherently noisy due to the randomness associated with the markers. As can be seen in higher marker-per-cell cases such as 500 markers per cell test, it does indeed slowly converge to the smooth results we see later on. In any case, this figure has been removed and its content moved to the results section.
Comment 9: Figures 4 and 5: Condensation
Response: Figure 4 and 5 have been removed.Comment 10: Line 140: Prior work context
"Is this your method? What similar things has been done in other codes before? Some general introduction text is needed here before all the subsections."Response: We have added introductory context at the head of this section, mapping out the architecture of our tracking routines and clearly cross-referencing how our approach relates to other's work such as LaMEM (Kaus et al., 2016) or Duretz et al. (2016).
Comment 11: Line 165 / Equations 5-6 / Moving equations to Appendix
"Since this is a method for which you have a reference, you can move some of the details to an appendix. Many of the following equations in section 3 are also standard can can be moved to appendix."Response: We have eliminated many of the unnecessary standard equations, preferring to cite them.
Comment 12: Line 381: Reference to slope limiting
Response: This was a novel method - the section has in any case been restructured to hopefully better reflect this.Comment 13: Line 392: Vague subsection title
Response: We have renamed this subsection to be much more descriptive of its technical content (e.g., updating it to clearly reflect the exact stabilisation or reinitialisation routine being detailed).Comment 14: Line 420: Clarify "coupling" definition
Response: We have expanded the text to define exactly what coupling implies in this context: it denotes the algorithmic transfer of the sub-grid free-surface position to the mechanical Stokes solver to continuously modulate cell-averaged density and viscosity fields, ensuring that topography variations dynamically feed back into the model's driving buoyancies.Comment 15: Figure 15: Structural location
Response: We have decided to retain the original locationComment 16: Figures 22 & 24: Label accuracy
"Change 'topography with time' to 'maximum topography with time'"Response: We have updated the figure captions to explicitly specify "maximum topography" rather than generic topography.
Comment 17: Line 503: Consolidating brief sections (6.1 & 6.5)
Response: We have removed these isolated, ultra-short sub-sections and integrated their contents directly into the introductory paragraphs of the Results section, framing our benchmarks with a cohesive rationale for why each numerical test was executed.Comment 18: Section 6.2: Framing a fair comparison
"How do you choose the number of tracers to make a fair comparison? In general I would appreciate if you clarify how you ensure your experiments constitute a fair comparison."Response: We have added an explanatory paragraph in Section 5 (Benchmarks) clarifying our comparative controls. To ensure a fair assessment, all tests are performed under identical spatial Eulerian grid resolutions, velocity fields, and time-stepping criteria. We explain how the volumetric marker densities were chosen to represent common community choices, demonstrating how our surface markers operate independently of full-domain marker counts to prove the scaling efficiency of our new approach.
Comment 19: Line 533: Value of the sticky air benchmark
"If the discrepancy could be due to the sticky air, what is the value of this benchmark? Clarify"Response: We have added text clarifying the value of this benchmark. While the sticky air approximation is known to possess minor, predictable intrinsic limits, this test remains highly valuable because it is the primary benchmark utilized by the broader community (e.g., the Crameri et al., 2012a community comparison). Replicating this test allows us to validate that our tracking methods match established community results perfectly while highlighting the significant noise reduction our method achieves over standard marker-in-cell setups.
Comment 20: Figures 25, 30, 32, 34: Layout space reduction
Response: We have removed several figures and reduced the space taken by others to make for a less cluttered paper.Comment 21: Expansion of the Reference List
"The reference list is quite short and with many papers written be the same authors."Response: We have expanded the reference list where appropriate.
Sincerely,
Timothy Gray, Paul Tackley, and Taras GeryaCitation: https://doi.org/10.5194/egusphere-2025-6546-AC2
-
AC2: 'Reply on RC2', Timothy Gray, 18 May 2026
Data sets
StagFreeSurface (StagFS): A testbed for free-surface methods in geodynamical simulations using the staggered-grid finite volume (difference) discretization. Paul Tackley and Timothy Gray https://doi.org/10.5281/zenodo.18096249
Model code and software
StagFreeSurface (StagFS): A testbed for free-surface methods in geodynamical simulations using the staggered-grid finite volume (difference) discretization. Paul Tackley and Timothy Gray https://doi.org/10.5281/zenodo.18096249
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 999 | 601 | 103 | 1,703 | 101 | 113 |
- HTML: 999
- PDF: 601
- XML: 103
- Total: 1,703
- BibTeX: 101
- EndNote: 113
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
I have read with great interest the manuscript "Lagrangian tracking methods applied to free surface boundaries in numerical geodynamic models" by Gray et al. The authors present two new and/or updated methods to allow for more accurate tracking of the free surface atop geodynamical models.
I found the article well written and structured and the figures clear and helpful. It is a detailed methodological paper which showcases many results obtained on benchmarks which I find convincing. As such I recommend publication with minor modifications.
I hereafter list all my comments:
- At line 35 we read: ``For clarity, the Lagrangian points already implemented in StagYY and distributed throughout the model domain are hereafter referred to as "tracers".''
I find this choice rather unfortunate: the standard method of tracking materials/fluids inside the entire domain in these codes is called the marker-in-cell technique (e.g. Gerya & Yuen, 2003, 2007) but in this work the authors choose to call such Lagrangian points 'tracers' while referring dozens of times in the text to the marker-in-cell method.
More confusion is added in my opinion through the choice of naming the additional set of points introduced to specifically track the surface 'markers' (line 37: ``In contrast, Lagrangian points introduced in this study specifically for tracking the free surface are referred to as “surface markers”.''
The next sentences at line 39 & 47 illustrate my point: ``In StagYY, as in many other codes, the standard technique for interface tracking employs Lagrangian tracers using a marker-in-cell method'' ; ``Using the marker-in-cell method, finer resolution can be achieved trivially by increasing the number of tracers.''
As such I suggest a change of names: markers for 'volume' points and tracers for additional surface points.
- l83: calculating -> calculated
- l88: C -> $C$'s
- Section 2.2: although the topic is addressed later with the surface markers techniques in section 4.2, I was wondering when reading this section about entrainment?
- l111: what is meant by 'surface mesh refinement' in the context of FDM codes here?
- l130: I think a reference for the 'divergence-free quadratic spline method' in this context would be welcome.
- l157: when citing Gerya's book please specify page or chapter
- l169: this is not a full sentence, I think a comma should be inserted in Eq 2, and the line should start with 'where ...'. Same at line 316 after Eq 16, line 321 after Eq 18, line 257 after Eq 11, line 265 after Eq 12.
- l224: the yin-yang grid is here mentioned but it has not been mentioned before. Most readers will know what it pertains to, but a ref could be added.
- l255: 'is' -> in
- l376-380: I am not sure of the point of the paragraph. Diffusion is mentioned, pros and cons are mentioned, but it is not clear to me whether this is implemented/used.
- l542: the semi-colon should probably replaced by a comma + i.e. ?
- l670: the word marker appears twice
- section 7.3.1: I think that some key publications of geodynamical codes coupled to surface processes could be cited in the context of this paragraph.
- somewhere in the text a reference to DOUAR (Braun et al, PEPI 171, 2008) should be added since the free surface was also tracked by means of a Delaunay-triangulated surface embedded in a Eulerian mesh.