| 13 Feb 2026
Volume of Fluid method applied to free surface boundaries in numerical geodynamic models
Abstract. Tracking the evolution of interfaces in numerical geodynamic models, particularly the rock-air/water boundary when imposing free surface boundary conditions, presents significant computational challenges. Traditional marker-in-cell methods, while widely used in mantle convection codes like StagYY, suffer from inherent noise due to random tracer distributions and require high tracer densities for accurate topography resolution, leading to substantial computational costs. This study presents the implementation of the Volume of Fluid (VOF) method into StagYY using the open-source gVOF package, offering a volume-conservative alternative for interface tracking in geodynamic simulations.
Two implementations are developed and tested: a basic method with one VOF cell per computational cell, and an improved method utilizing subcells with full coupling to the Stokes solver through density and viscosity fields. Both methods employs the CLCIR (Conservative Level Contour-based Interface Reconstruction) reconstruction scheme combined with FMFPA (Face-Matched Flux Polyhedron Advection) for optimal accuracy. A consolidation approach for translating VOF fields to Eulerian surface locations ensures volume conservation and numerical stability.
Benchmarking in 2D and 3D across multiple geometries demonstrates that while the basic VOF implementation provides adequate results for simple 2D cases, the improved method is essential for accurate tracking in complex flows and three-dimensional applications. The improved VOF method successfully eliminates topographic noise inherent in tracer-based approaches and decouples surface quality from tracer density, enabling high-quality results with reduced computational overhead from tracer advection.
Despite increased computational costs and memory consumption associated with the method, particularly with the improved implementation, the VOF method offers distinct advantages including explicit volume conservation, sharp interface representation, and versatility for tracking various geophysical boundaries beyond free surfaces. This work establishes a foundation for future applications in sea level modelling, continental margin tracking, and coupled planetary system simulations, advancing toward global-scale biogeodynamic modelling.
General comments:
The study “Volume of Fluid method applied to free surface boundaries in numerical geodynamic models” by Gray et al. describes the implementation of the “Volume of Fluid” method in the geodynamic code StagYY. The goal of the method is to improve the accuracy of the free surface by replacing the use of the marker-in-cell method for this purpose. As such, the topic of the manuscript fits well into the scope of GMD.
Unfortunately, the presentation of the content is in very poor shape. It looks like a first paper draft with little thought and effort put into effective communication and proof reading. This starts with a number of harmless mistakes and a lack of supporting references. More importantly, many explanations are missing or incomplete. Figure captions are often uninformative, Figure legends poor, and there are too many individual Figures that should be combined into one. In the Results sections, captions are mostly interpretation instead of description. At the same time, there is hardly any text in the Results section. There are also too many subsections which are extremely short. Some concepts are introduced and then never mentioned again (sampling method, tracer bouncing). Details are given below.
The content might be worthy of a publication but the presentation speaks volumes about the lack of care and effort put into this submission. I encourage the authors to add missing explanations, combine Figures, make useful legends, write meaningful captions, and entirely rewrite section 4. Then, the quality of the manuscript can be reevaluated.
Overarching issues:
The main method and its implementation are not well explained. To me it is unclear how Figure 2 relates to the actual surface tracked in the code? Several horizontal coordinates have multiple vertical coordinates like between the third and forth cell from the left. How is the colour function (eq. 1) actually implemented and discretized? On a tracer level? Then advected and summed up per cell? How is parameter n reconstructed? Instead of explaining these things, sections 2.1, 2.2, 2.3 mostly just list names of options in gVOF and which ones where picked.
One of the main messages is the coupling of the surface to the stokes equations in the “improved VOF method”. Yet, it is unclear what this coupling actually means. The only explanation is section 3.5 and equation 5. Both of which sound like density and viscosity are weighted averages of the present phases in the cell or control volume. If this is unique to the improved method, then how is this dealt with in the other methods? Simply copy the closest tracer?
The paper is missing a description of how the marker-in-cell approach handles the free surface. This makes statements about “floating” cells like line 114 difficult to understand. The introduction should explain how the free surface is handled traditionally, and what issues this method has. Then, it is easier to understand how the VOF method improves upon it. Furthermore, the authors appear to use the terms “tracer” and “marker” synonymously which should be clarified once.
The sampling method is mentioned in section 3.3, plotted in Figure 4, but never explained or mentioned again. Same for tracer bouncing.
One of the biggest issues for geodynamic models with free surface is the “drunken sailor problem” (e.g., Rose et al., 2017, Kaus et al., 2010). This is never mentioned. How does the VOF method affect this? Some Figures give hints but none of them show velocities over time and the text does not discuss it at all.
Almost all Figure captions are problematic. Early on, they are too short and don’t describe everything that is relevant. In the results section, they contain too much discussion and interpretation. These things belong in the main text. Captions should only be descriptive. A side effect of this is that there is hardly any text in the results section. What text is left is constantly broken by new headings.
Figures 4, 15, 17, 20 show surface topography for presumably the same model setup. It is unclear at which time in the model these are.
When the VOF method is evaluated in Figures 17 and 20, the results are called consistent which is simply not the case. If anything, Figure 20 is concerning as a change in tracer density appears to result in completely different topographic patterns while the text uses independence from tracer density as one of the main benefits of the VOF method. Biogeodynamics is often mentioned as a field that would benefit from this improved free surface implementation, but looking at Figure 20 suggests that a numerical parameter can alter the position of mountain chains and therefore climate zones. Also, how does the marker-in-cell approach perform in these cases in comparison?
The marker-in-cell approach is often interpreted to suffer from tracer noise but the manuscript never explains what this term means.
There are too many Figures. A lot of them could easily be combined to make thematic Figures with subfigure panels. E.g., 1+2+3, 6+7+8, 13+14 and so on
References to equations and Figures are not in line with the journal guidelines (https://www.geoscientific-model-development.net/submission.html)
Specific comments:
Affiliation: “Zürich” and “Zurich”
L. 8: remove “s” at the end of “employs”
L. 9: “reconstruction” is repeated
L.15 and L.17: One mentions “reduced computational overhead” and the other “increased computational cost”, but both appear to describe the VOF method. Please clarify.
L.26: provide references that support this statement.
L.27: If this is supposed to be a general statement, remove “the” in the middle of the line and mention other codes like ASPECT, LaMEM, Underworld, pTatin. If the statement is specific to StagYY and i3elvis, restructure the citations.
L.31: Seems like an unnecessary statement. Is there any planet where this would not be the case?
L.33: Seems like a good spot to cite Crameri, 2012.
L.43: References or links for OpenFOAM and ANSYS?
L.45: same for ASPECT.
L.46: where are the applications for other interfaces?
L.49: missing comma after “models”
Eq. 2: define “t”
Fig.2: I assume this is approximating the the interface from Fig. 1. If that’s the case, please plot it here as well.
L.67: Interface normal reconstruction.
L.69: So alpha is not constant in time? Also not in space?
L.75: Explain how n and alpha are reconstructed.
L.78: Should be V_i as we are talking about a single cell?
L.81: Which method are you using?
L.93: Please clarify why subcells allow coupling with the Stokes solver. I assume density and viscosity are also available on a cell-size grid and not sub-cell-sized grid. This is eventually cleared up by Figure 8. Explain it here.
L.95: wrong citation style
Fig. 3: Only shows how it works in 2D despite caption and referencing text stating 3D as well.
L.102: Right cell in Fig. 3 has 3 planes, not 2.
L.106: Sampling method is introduced and then never explained.
L.108: Without 3.3.2, there should be no 3.3.1.
Fig. 4: This should be in the results section, not in the methods section. Furthermore, are both methods used for the same problem here? The results look completely different. Consolidation method also appears to produce an undershoot in the lower left part. Colorbar seems to be missing minus signs. Could use a diverging colormap.
L.118: good place to cite Duretz et al., 2016. Also how did StagYY deal with this problem before? How about other codes?
Fig. 5: Figure is quite ambiguous. Left: Why are arrows going both ways? Which tracer of what phase is going where? Right: What meaning do inner and outer color have?
L.120: What does reflect mean in this context? If a rock tracer was 1 km above the surface, it is moved 1 km below it? Why not just a minimal distance below?
L.122: Why not use a rock tracer within the cell, but below the surface?
L.123: As far as I can see, tracer bouncing is never mentioned again. Why even introduce it?
Fig.6: The point of this Figure is to show where viscosity, density and volume fractions are defined in space. None of them are listed in the legend….
L.129: what is z-momentum equation? Vertical momentum? Where is the equation in the manuscript?
L.129-135: This begs the question how the code generally handles cells at any interface where tracers with different phases are present within the same cell. I assume the values of the phases are averaged which is what Eq. 5 does assuming rho_air = 0. So it would not be unique to the presented method and not even to StagYY.
L.136: If there is no 3.5.2, there should not be a 3.5.1. And it is pointless to have a sub-section that is two lines long.
Fig.7: Same problem as Fig. 6 + random label “v_z momentum equation” + Caption does not describe the Figure. Finally, why is this not combined with Figure 8 and drawn in the same style?
L.142: As mentioned before, the “coupling” seems to not be more than averaging of quantities onto locations in the grid. Is this not automatically part of the code for all equations and cells?
L.144: The 3 benchmarks should have already been mentioned in the introduction or methodology. Tab.1: Table title appears again underneath table. Third benchmark is mantle convection without crust? No definition of Ra.
L.153: viscosity = 10^18 Pa s, also: This is just the air (what viscosity are the other units then?) or all of it?
L.153: what are symmetric boundary conditions?
L.154: gravitational acceleration
Fig.9: Explain dashed line in caption. Red and green should be avoided together. This also affects many of the following Figures.
L.155: The wording implies that the result would be different with regularly spaced tracers, is this the case?
L.159-163: Same issues as previous benchmark + it should be Myr, there are many of these in the rest of the manuscript.
Fig.10: Caption should describe what the Figure shows.
L.164-165: To me, this sounds like the problem is not marker-in-cell in general, but more so marker-in-cell in StagYY. Or do none of the models in Crameri et al., 2012 use marker-in-cell? Please clarify.
L.166: Unclear.
Fig. 11: Why is this not part of Figure 10? Again, caption is interpretation and discussion instead of description of the Figure. Nothing in this Figure suggests the existence of noise. Missing citation.
L.170: Paragraphs should have more than one sentence. Are the results consistent between the approaches? Was marker-in-cell used here as well?
Fig.12: What does reconstruction refer to?
L.181: Material parameters.
Tab.2: Table title appears again underneath table. Figure 13 suggests, that the benchmarks were not only done for VOF.
Fig.14: Should be part of Fig. 13. Caption: There are no small-scale variations for the red line to capture.
L.183: Delete “high-quality”
Fig.15: See comment for Fig. 4.
Fig.16: Should be part of Fig.15. Density is not the correct label for the y-axis. Caption lacks essential information and has interpretation again. If steep peak indicates superior accuracy, then Figure 18 shows that low resolution produces superior accuracy. Explain.
L.185-187: If 2.5 lines is all there is to say here, then you might as well delete this section.
L.189-197: Text is closer to a bullet point list. Does not even describe model behavior. Rewrite.
Fig.17: Same problem as Fig.4. Topography is not consistent across resolutions. This needs to be discussed.
Fig.18: Same issues as before.
Fig.20: Caption is simply not true. Same for discussion in the text.
Sec.5.1: I am not generally against bullet-style lists for summaries, but a numbered list would look better. Same for 5.3.
L.206: Abstract mentions that the presented content is also applicable to other interfaces beyond free surface but consolidation method does not work for those. Discuss.
L.208: Figure 10 suggests that basic VOF is already insufficient in 2D.
L.213: None of the Figures show a performance or memory comparison with marker-in-cell. Also, lines 224-226 seem to contradict this statement.
L.227: Would that not remove the advantages of the improved method again?
L.228: So currently the VOF method is applied throughout the entire model domain? Why? This seems straight forward to avoid.
L.229: Elaborate. What information is currently stored redundantly?
L.230: Unclear, please elaborate.
L.231: Please explain what yin-yang geometry means
Sec.5.4: Using “sub-headers” again. Either make subsections or integrate the applications into the text.
L.236: Would be problematic with the consolidation method, no?
L.244: References for biogeodynamics?
L.245: This seems to make the same point as line 236.
References:
Duretz, T., May, D. A., & Yamato, P. (2016). A free surface capturing discretization for the staggered grid finite difference scheme. Geophysical Journal International, 204(3), 1518-1530.
Kaus, B. J., Mühlhaus, H., & May, D. A. (2010). A stabilization algorithm for geodynamic numerical simulations with a free surface. Physics of the Earth and Planetary Interiors, 181(1-2), 12-20.
Rose, I., Buffett, B., & Heister, T. (2017). Stability and accuracy of free surface time integration in viscous flows. Physics of the Earth and Planetary Interiors, 262, 90-100.