the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 28 Dec 2025
Variational Stokes method applied to free surface boundaries in numerical geodynamic models
Abstract. Accurately and efficiently modelling topographic evolution is a key challenge in geodynamic modelling, which requires the solution of the Stokes equations with free surface boundary conditions. While finite difference methods on staggered grids, as used in geodynamic modelling codes such as StagYY, I3ELVIS and LaMEM, offer strong computational performance and compatibility with multigrid solvers, the use of fixed Eulerian grids complicates the implementation of realistic, deformable free surfaces. Two existing methods are available to model free surface boundary conditions in StagYY: the commonly used sticky-air method, which suffers from limitations relating to high viscosity contrasts, and the "staircase" method, which improves upon the sticky air method by imposing free surface boundary conditions at cell boundaries.
To address the limitations of existing methods of implementing free surface boundary conditions, this study investigates an alternative variational discretisation of the Stokes equations that uses volume fractions to represent a smooth surface within a fixed Eulerian grid, allowing the imposition of accurate free surface boundary conditions while allowing it to bypass the limitations of existing free surface discretisation methods.
The variational Stokes method is demonstrated to be an accurate and computationally efficient alternative to existing methods. It reproduces results comparable to existing methods while reducing computational cost and enabling broader applications, including non-zero surface tractions, complex surface loading, and compatibility with 3D spherical geometries.
- Preprint
(3525 KB) - Metadata XML
- BibTeX
- EndNote
Status: open (until 22 Feb 2026)
-
CEC1: 'Comment on egusphere-2025-6354 - No compliance with the policy of the journal', Juan Antonio Añel, 07 Jan 2026
reply
-
AC1: 'Reply on CEC1', Timothy Gray, 10 Jan 2026
reply
Dear Juan,
A version of StagYY was created to accompany this manuscript, which was not yet ready at the time of submission. The code can be found here: https://zenodo.org/records/18096250
Kind regards,
Timothy Gray
Citation: https://doi.org/10.5194/egusphere-2025-6354-AC1 -
CEC2: 'Reply on AC1', Juan Antonio Añel, 11 Jan 2026
reply
Dear authors,
Many thanks for your reply. However, it would be good to clarify if there is a confusion in what you have published in the above mentioned repository. There, you claim to have published the StagFS code, but in the manuscript you refer to the StagYY code. It looks like if both are different things, and we where requesting the StagYY code that you use to develop your manuscript. Therefore, please, clarify if the StagFS code that you have published is the StagYY code, and if yes, the causes for the different names.
Also, we understand that you did not have the code ready at the submission time. However, in such case, you should have not submitted your manuscript to the journal until you had the code released. Publication of code is not an additional requirement in GMD, but the first part that must be accomplished before submitting a manuscript, and your manuscript should have never been accepted for Discussions or review given such flaw (unfortunately, it was overlooked in this case). Please, be aware of it for potential future submissions.
Juan A. Añel
Geosci. Model Dev. Executive Editor
Citation: https://doi.org/10.5194/egusphere-2025-6354-CEC2 -
AC2: 'Reply on CEC2', Timothy Gray, 20 Jan 2026
reply
Dear Juan,
Thank you for your input. We agreed that there may be some confusion regarding the supplied code StagFS (Stag Free Surface), the name of which was intended to convey that it was a limited feature set version of StagYY intended for benchmarking of free surface methodologies. To avoid further confusion, we have renamed the supplied code with the somewhat more descriptive name StagYYFreeSurface (StagYYFS).
Kind regards,
Timothy Gray
Citation: https://doi.org/10.5194/egusphere-2025-6354-AC2
-
AC2: 'Reply on CEC2', Timothy Gray, 20 Jan 2026
reply
-
CEC2: 'Reply on AC1', Juan Antonio Añel, 11 Jan 2026
reply
-
AC1: 'Reply on CEC1', Timothy Gray, 10 Jan 2026
reply
-
RC1: 'Comment on egusphere-2025-6354', Wolfgang Bangerth, 08 Feb 2026
reply
Review of
Variational Stokes method applied to free surface boundaries in
numerical geodynamic models
by
T. Gray, P. Tackley, and T. Gerya
This manuscript investigates the use of a variational statement of the
inclusion of free surfaces with the Stokes equations, in the context
of staggered finite difference methods. While I think that in general
such a paper would be of interest, *this* paper is not ready and
the authors need to fundamentally rethink how they want to present
this material. Specifically, my opinion is based on the following
thoughts that I will elaborate upon below:1/ The paper is far too long and poorly structured.
2/ The paper does not actually describe the method in question in a
way that would be understandable by readers.
3/ It is unclear how the method relates to others in the field because
most of the many other methods are not mentioned or referenced.
4/ It is unclear how the method shown in the paper differs from the
one in Larianov et al. that is mentioned in various places as the basis
for the work herein. As a consequence, I cannot determine what is
actually new.Let me address these points in more detail:
1/ The paper is 41 pages long. It reads like a chapter of a thesis,
which I originally thought might just be because the author who wrote
most of it is inexperienced. But based on the first sentence of
section 6, it may really just literally be a chapter of a
thesis. Regardless, the key issue is that a thesis is generally a
complete record of the work a student has done as part of their
research, and so it will typically show *every* numerical experiment
that was performed. A paper, on the other hand, is a way for
researchers to *teach* other researchers about what they figured
out. For this, it is not necessary -- and often counterproductive --
to show everything that was done. Rather, one synthesizes the
knowledge gained and picks the experiments that most clearly show what
one has learned. One also does not have to explain in detail
alternative methods, and it should be possible to introduce the method
one wants to use earlier than just on page 8.For the paper here, given that the method is pretty straightforward, I
see no reason why the paper should be much more than 20 pages long --
about half of what it has now. I can see how the material the authors
show here could be made into a paper that is useful for the community,
but I think it would be useful to first create an outline (i.e.,
basically just the section and subsection headings) of a paper into
which a certain amount of thought is put to determine what really is
necessary and what is not. An outline that is a copy of a thesis
chapter is not usually a good starting point.
2/ I cannot figure out what the method actually is. The method is
introduced in Section 3, starting on p. 8, and is based on the
variational formulation provided in eq. (11) as an integral over
the fluid domain (which may change from time step to time
step). Following eq. (12), the authors then simply state
"The variational formulation (Equation 11) may be discretised."
by which they mean that they replace integrals by products between
vectors and matrices. The key question everyone will face at this
point, but that the paper despite its length does not actually answer
is: How exactly are they to be discretized? From earlier sections, one
can infer that the authors use a finite difference discretization,
which defines the solution only at points on a staggered grid. But how
does that extend to the definition of integrals then? In the finite
element context, one defines the solution as *functions* over the
domain, and so it is clear how integrals are computed. If one
re-interpreted the stagged mesh finite difference method as a mixed
finite element method, then that would apply as well. But in the current
context, it is not clear to me how one does that. At some point, the
authors introduce a volume fraction factor phi for each point, and
then simply multiply certain terms by it, but it is not clear whether
phi is a scalar, a matrix, or anything else, and in any case, the
transition from integrals to vectors/matrices remains unclear.Given that a complete description of the method to be investigates
should be a central aspect of the paper, the omission of such a
description appears like a key oversight.
3/ There are many other methods, predominantly in the finite element
field, that deal with boundaries that intersect cells. To name a few,
one could look at the CutFEM and immersed boundary
methods. Fundamentally, these start from the same point, namely by
defining the variational formulation on only a subset of the
computational domain, and as a consequence I think it would have been
useful to contrast and compare them to the method presented
herein. But beyond just the comparison, methods such as the CutFEM or
immersed boundary methods have shown that a key issue is when only
small slivers of cells are intersected by the boundary; in those
cases, the condition number deteriorates because it is proportional to
one over the smallest cell fraction cut off by the boundary. One would
suspect that the same happens for the method here, and that perhaps
that is a reason for the otherwise unexplained failures of numerical
solvers? Given that that is a well-known and well-understood issue
(also with remedies) in the finite element context, it seems like a
missed point to not make that connection in the current paper.
4/ Finally, the manuscript suggests that it implements the method by
Larionov et al. (2017). But how exactly does the current paper differ
from Larionov? In other words, what is it that is *new* herein?
Beyond this, here are number of other, minor issues:* I think it would be useful to make clear in the title that the method
is specifically intended for the *finite difference* method.* Typographically, formulas are just a regular part of sentences, even
if they are offset. As a consequence, if the sentence ends with the
formula, the formula should end with a period. If the sentence
continues after that, a comma may or may not be appropriate, but in
any case the text that follows should start with a lower-case letter
and not be indented as if a new paragraph starts. (In LaTeX, that
means that there should not be an empty line between the formula and
the following text since empty lines indicate paragraph breaks.) I
have marked this up in the attached PDF file in many places, but
eventually gave up -- there are likely many more places.* Section 3.1 introduces the method for the case without body forces,
and then spends another half page in Section 3.3 on the case of body
forces. Since you're always short on space when writing papers, just
include the general case in the exposition of 3.1. The same could be
said for the use of traction boundaries.* Section 3.5 seems unnecessary. It takes up 2/3 of a page, but
nothing that precedes this section was specific about Cartesian
meshes and it seems obvious to a reader that the method should also
work on non-Cartesian meshes; in fact, I suspect that that is also
the case for general unstructured meshes.* The issue of volume fraction factors is key to the method. Section
4.1 starts by saying that there are multiple methods to compute phi,
but then only shows one. What are the others, and why are they not
discussed here?* Section 4.2 about air tracers seems unconnected to the rest of the
paper. It discusses some practical concerns, but I took the paper as
discussing and evaluating a numerical method, and the questions
discussed in 4.2 do not seem germane to this goal. (This falls into
the description of how to write papers above: Create an outline
first in which you wonder what it is that needs to be shown and
discussed, and then write about that. Everything that is not
*necessary* is by definition *unnecessary* and should be
omitted.)* Section 5: This section shows far too many numerical experiments. I
did not read through all of them because it was not clear to me how
they all differed and why they all needed to be shown. I understand
the impulse to show everything that was worked on during a project,
but that does not make for good papers. Ask which experiment *needs*
to be shown to demonstrate *specific* aspects of the method, to
illustrate when it works and when it doesn't, and how well. I
suspect that you can get away with just 2-3 examples between the 2d
and 3d sections.* In Figures such as Fig. 10, show errors and mesh sizes on the axes,
not their logarithms. All visualization tools have modes where they
can show values on logarithmic scales, rather than logarithms on
linear scales.* In a similar vein, show Fig. 24 on a log scale.
* If a method does not converge, that is a major issue. In Section 5.3
(and perhaps others), you encounter such a case and just brush it
aside by saying "in this case noise associated with the tracer
distribution prevents the surface from fully relaxing". I don't
understand this, and I think that as a reader this is quite
unsatisfying: What exactly is this noise, and how would I prevent
that from happening? In essence, what this example shows is that
when you use a small number of particles per cell, the method
converges unreliably, and when you use a large number of particles,
it does not converge at all. This warrants a discussion of more than
a half sentence.* Section 6 (Summary and conclusions) consists of bullet points
(either directly, or in the form of very short sections). This does
not make for exciting reading. I also don't think that the section
should 3.5 pages.
Data sets
Software accompanying manuscript "Variational Stokes method applied to free surface boundaries in numerical geodynamic models" Timothy Gray https://doi.org/10.5281/zenodo.17956246
Model code and software
Software accompanying manuscript "Variational Stokes method applied to free surface boundaries in numerical geodynamic models" Timothy Gray https://doi.org/10.5281/zenodo.17956246
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 170 | 105 | 23 | 298 | 24 | 42 |
- HTML: 170
- PDF: 105
- XML: 23
- Total: 298
- BibTeX: 24
- EndNote: 42
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Dear authors,
Unfortunately, after checking your manuscript, it has come to our attention that it does not comply with our "Code and Data Policy".
https://www.geoscientific-model-development.net/policies/code_and_data_policy.html
Specifically, you do not provide a repository for the StagYY code, which you use in your work. We can not accept this. You must publish openly all the code and data used in the development of your manuscript.
The GMD review process depends on reviewers and community commentators being able to access, during the discussion phase, the code and data on which a manuscript depends. Please, therefore, publish the StagYY code in one of the appropriate repositories and reply to this comment with the relevant information (link and a permanent identifier for it (e.g. DOI)) as soon as possible. We cannot have manuscripts under discussion that do not comply with our policy.
The 'Code and Data Availability’ section must also be modified to cite the new repository locations, and corresponding references added to the bibliography.
I must note that if you do not fix this problem, we cannot continue with the peer-review process or accept your manuscript for publication in GMD.
Juan A. Añel
Geosci. Model Dev. Executive Editor