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.
-
Notice on discussion status
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
-
Preprint
(3525 KB)
-
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(3525 KB) - Metadata XML
- BibTeX
- EndNote
- Final revised paper
Journal article(s) based on this preprint
Interactive discussion
Status: closed
-
CEC1: 'Comment on egusphere-2025-6354 - No compliance with the policy of the journal', Juan Antonio Añel, 07 Jan 2026
-
AC1: 'Reply on CEC1', Timothy Gray, 10 Jan 2026
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
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
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
-
CEC2: 'Reply on AC1', Juan Antonio Añel, 11 Jan 2026
-
AC1: 'Reply on CEC1', Timothy Gray, 10 Jan 2026
-
RC1: 'Comment on egusphere-2025-6354', Wolfgang Bangerth, 08 Feb 2026
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.-
AC3: 'Reply on RC1', Timothy Gray, 22 Mar 2026
Dear Wolfgang,
Thank you for your thorough and constructive review of our manuscript. We have carefully considered all comments and have made substantial revisions in response. Below we outline how each point has been addressed.
1. Paper length and structure
We have substantially restructured and shortened the manuscript from 41 to 26 pages. The thesis-like presentation has been replaced with a more focused paper structure. Background material on existing methods has been condensed, and the number of numerical experiments has been reduced to those that most clearly demonstrate the key properties of the method.2. Description of the method
Section 3 has been rewritten to more clearly explain how the continuous variational formulation is discretised on the staggered grid. We now explicitly describe how the volume integrals are approximated cellwise, how the volume fraction phi enters each term, and how the discrete operators G, D, and M relate to the standard staggered-grid discretisation. The stress elimination procedure leading to the final pressure-velocity system is also set out explicitly.3. Relationship to other methods
We acknowledge that a fuller comparison to related finite element approaches would strengthen the paper, and we have added a brief discussion of these methods and how they relate to the approach presented here.4. Relationship to Larionov et al. (2017)
We have clarified in the introduction the specific contributions of this paper relative to Larionov et al. (2017), which are the implementation in a geodynamic FD code with multigrid, and systematic comparison against sticky air and staircase methods.Minor points:
- The title has been updated to make clear that the method is specifically for finite difference discretisations.
- Equation punctuation has been corrected throughout the manuscript.
- Sections 3.1 and 3.3 have been merged so that the general case including body forces and traction boundary conditions is presented from the outset.
- Section 3.5 on non-Cartesian meshes has been removed.
- Section 4.1 has been revised; the claim that multiple methods exist for computing volume fractions is retained but the other methods are now briefly named, with a note that only Lagrangian tracking is used in this study.
- Section 4.2 on air tracers has been removed.
- The results section has been substantially reduced. We now present key benchmarks in 2D and 3D rather than the full set of experiments previously shown.
- Error and mesh size axes in convergence plots now use logarithmic scales rather than showing logarithms on linear scales.
- The non-convergence issue in the tracer distribution test has been discussed more fully.
- The summary and conclusions section has been rewritten in prose rather than bullet points and shortened to a more appropriate length.We hope that these revisions adequately address your concerns. We would be happy to make any further changes if needed.
Yours sincerely,
Timothy Gray, Paul Tackley, and Taras Gerya
Institute of Geophysics, ETH ZürichCitation: https://doi.org/10.5194/egusphere-2025-6354-AC3
-
AC3: 'Reply on RC1', Timothy Gray, 22 Mar 2026
-
RC2: 'Comment on egusphere-2025-6354', Thibault Duretz, 12 Feb 2026
Review of Gray et al: “Variational Stokes method applied to free surface boundaries in numerical geodynamic model”
I have read the manuscript of Gray et al. with interest. Accurately capturing planetary surfaces is a particularly relevant question in geodynamic modelling, with far-reaching implications.
This study discusses the application of a new method to treat free-surface flow in geodynamics: the variational finite-difference approach. The method is shown to provide several improvements over existing approaches, which is a strong and positive outcome of the paper.
Nevertheless, the current form of the manuscript is not ideal. It is too long, and while the focus should be on the novelty of the approach, the main message becomes somewhat diluted. At several points, the text reads more like a review article, which detracts from the emphasis on the new contributions. I have made a number of comments below that I believe could help improve the structure and clarity of the manuscript and shorten it where appropriate.
Thibault Duretz
——————————————————————
Variational form of Stokes equations can be used in various contexts, the paper title should be more specific and clearly state the variational finite differences approach.
l 32 “on on the latter”
l 40 also: fast flow of air can lead to reduction of time steps.
l 49-51 I am not aware of this problems and did not experience this in my GMG implementation (which was only 2D). This may be specific to the StagYY implementation.
l 53 kinds of confirm what I state just above.
Section 2.1 the authors maybe curious about reading this: https://api.unil.ch/iris/server/api/core/bitstreams/54c253b1-d357-455c-922a-3f3e9cd14aae/content
Section 2.1 is a review of the sticky air method. I find some text to be redundant with the intro (disadvantages of sticky air). This section could be shortened, there a lot of papers on sticky air and space should be left to set the focus on the variational FD scheme.
Section 2.2 is a review of the “staircase” approach. It coud maybe be shortened or rather being used to introduce the variational scheme. As far as I see, the variational formulation ends up with a similar discretisation but features the additional weights. These weights increases the accuracy, thereby improving the “staircase” method. Wouldn’t it be a way to frame it?
2.2.1 and 2.2.2 take a lot of space and delay the moment to which the reader reaches the core of the manuscript: the variational FD approach
2.2.2 repeats some information that is already given earlier in the test (limitation of staircase). Instead of making it sound like a review of methods for free surface, one could actually parts of the text to actually motivate the development of the variational FD scheme, where the focus is supposed to be set.
around l 200. again I am not convinced with this statements about “numerical instabilities”. I ran successful GMG models in my 2016 paper. Also, we have plenty of studies that relied on this approach (basically all papers published using MDoodz), which definitely do not feature instability cause by free surface discretisation. I think you may just have some implementation problems and I would tone down this statement.
Section 3 Variational FD Stokes method arrives too late (it’s already page 8!). Everything which is above can be condensed and rephrased in the way that motivates the development and testing of the variational FD Stokes method.
Section 3.1 Stokes equation are already stated as Eq 1. Note that there is no body force anymore in Section 3.1 for some reason.
Reading this, it feels that this section should be in fact the start of the study (right after intro) and that the above text about other (established methods) could be used for discussion purposes
Figure 3 I would gray out all the dots that are inactive for this material configuration - or are they all active
l 289 you’re discussing benefits over other approach, I would do that in the discussion.
Section 3.3 is about right hand side but the body force disappears form equation 9 onwards.
l 299 this is confusing: there is an air layer although there is no sticky air?
Eq 24 and 25 feature no equality. Maybe introduce symbols for these quantities.
section 4. I’m not sure we need all the details of the implementation within StagYY. Better focus on the results
Section 5. Results section only land on page 16. Definitely, what is above should be shortened.
around l 412 this is a lot of discussion about the sticky air method, that I barely can follow. It can be interesting in a supporting material but here the reader expects to focus on the new scheme and its benefits.
Fig 7-8-9 numerical resolution should be mentioned.
Which interpolation is being used for the sticky air method?
Fig 11-12 can be merged
Sec 5.3.2 feels over the top because the volume tracking is yet another problem (interesting too)
l 473 this conclusion is the same as the one from Tackley and King mentioned earlier? Maybe move this section a supp mat if it reproduces previous/known results?
Fig 16 feels like it does not belong to the study or should be moved to an appendix
Fig 18-19 when solid lines overlap, consider using markers/symbols instead.
For the subduction section, I think you only need Fig 21 merge with 23 (21 is too much)
l 543 you only identify here a problem that you mentioned twice in the introduction (convergence of MG solver with staircase). Therefore in the introduction, you cannot assume the problem is known to the reader.
Figure 25 is a weird figure, that appears to show 3 times the same thing. If results look all the same, one should heat for differences, not show three times the same thing. I would simply remove this figure
maybe you could remove this 3D plume benchmark section. It’s not clear why it’s sitting in between 2 scaling tests, and it’s not useful to identify differences between approaches.
Figure 27, caption: the sentence “the staircase method shows variable performance results.” is imprecise as all methods show variable performance.
To my opinion all scaling tests could be combine in one section and 1 figure.
Section 5.10 is really cool and useful application of the approach of Larionov. Somehow it ends up a bit lost between the scaling test and the discussion. I am not sure you need both 3D and 2D tests though (put the 2D in a sump mat?)
Section I would avoid using bullet points here.
6.5.1 “Applicability to incompressible flows” I think it should be “Applicability to compressible flows”
Citation: https://doi.org/10.5194/egusphere-2025-6354-RC2 -
AC4: 'Reply on RC2', Timothy Gray, 22 Mar 2026
Dear Thibault,
Thank you for your thorough and constructive review of our manuscript. We have carefully considered all comments and have made substantial revisions in response.
We have substantially restructured and shortened the manuscript from ~41 to ~26 pages. The thesis-like presentation has been replaced with a more focused paper structure. Background material on existing methods has been condensed, and the number of numerical experiments has been reduced to those that most clearly demonstrate the key properties of the method.
More specific comments have been addressed as such:
- The title has been made more specific to clearly reflect the variational finite difference approach.
- The typographic error "on on the latter" has been corrected.
- The discussion of fast air flow and its effect on timesteps has been acknowledged and the relevant text adjusted.
- Statements about numerical instabilities have been toned down to avoid overstating their generality beyond the specific StagYY implementation context.
- Sections 2.1 and 2.2 on existing methods have been shortened and reframed so that they motivate the development of the variational scheme rather than reading as a standalone review. The suggested reference has been incorporated.
- Section 3 now begins earlier in the manuscript, immediately following a condensed introduction to existing methods.
- The body force term is now included consistently from Section 3.1 onwards.
- Implementation details specific to StagYY have been reduced in Section 4, with the focus shifted to results.
- The results section has been streamlined; several figures have been merged or removed, and figures that previously showed overlapping solid lines now use markers or symbols for clarity.
- Scaling tests have been consolidated into fewer figures.
- The subduction benchmark section has been revised; the discussion of near-surface viscosity errors has been expanded and connected more explicitly to the known sensitivity documented by Schmeling et al. (2008).
- The 3D plume benchmark section has been removed.
- Bullet points in the conclusions have been replaced with prose.
- The typo in the conclusions ("incompressible" corrected to "compressible") has been fixed.We hope that these revisions address your concerns satisfactorily. We would be happy to make any further changes if needed.
Yours sincerely,
Timothy Gray, Paul Tackley, and Taras Gerya
Institute of Geophysics, ETH ZürichCitation: https://doi.org/10.5194/egusphere-2025-6354-AC4
-
AC4: 'Reply on RC2', Timothy Gray, 22 Mar 2026
Peer review completion
Interactive discussion
Status: closed
-
CEC1: 'Comment on egusphere-2025-6354 - No compliance with the policy of the journal', Juan Antonio Añel, 07 Jan 2026
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 EditorCitation: https://doi.org/10.5194/egusphere-2025-6354-CEC1 -
AC1: 'Reply on CEC1', Timothy Gray, 10 Jan 2026
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
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
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
-
CEC2: 'Reply on AC1', Juan Antonio Añel, 11 Jan 2026
-
AC1: 'Reply on CEC1', Timothy Gray, 10 Jan 2026
-
RC1: 'Comment on egusphere-2025-6354', Wolfgang Bangerth, 08 Feb 2026
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.-
AC3: 'Reply on RC1', Timothy Gray, 22 Mar 2026
Dear Wolfgang,
Thank you for your thorough and constructive review of our manuscript. We have carefully considered all comments and have made substantial revisions in response. Below we outline how each point has been addressed.
1. Paper length and structure
We have substantially restructured and shortened the manuscript from 41 to 26 pages. The thesis-like presentation has been replaced with a more focused paper structure. Background material on existing methods has been condensed, and the number of numerical experiments has been reduced to those that most clearly demonstrate the key properties of the method.2. Description of the method
Section 3 has been rewritten to more clearly explain how the continuous variational formulation is discretised on the staggered grid. We now explicitly describe how the volume integrals are approximated cellwise, how the volume fraction phi enters each term, and how the discrete operators G, D, and M relate to the standard staggered-grid discretisation. The stress elimination procedure leading to the final pressure-velocity system is also set out explicitly.3. Relationship to other methods
We acknowledge that a fuller comparison to related finite element approaches would strengthen the paper, and we have added a brief discussion of these methods and how they relate to the approach presented here.4. Relationship to Larionov et al. (2017)
We have clarified in the introduction the specific contributions of this paper relative to Larionov et al. (2017), which are the implementation in a geodynamic FD code with multigrid, and systematic comparison against sticky air and staircase methods.Minor points:
- The title has been updated to make clear that the method is specifically for finite difference discretisations.
- Equation punctuation has been corrected throughout the manuscript.
- Sections 3.1 and 3.3 have been merged so that the general case including body forces and traction boundary conditions is presented from the outset.
- Section 3.5 on non-Cartesian meshes has been removed.
- Section 4.1 has been revised; the claim that multiple methods exist for computing volume fractions is retained but the other methods are now briefly named, with a note that only Lagrangian tracking is used in this study.
- Section 4.2 on air tracers has been removed.
- The results section has been substantially reduced. We now present key benchmarks in 2D and 3D rather than the full set of experiments previously shown.
- Error and mesh size axes in convergence plots now use logarithmic scales rather than showing logarithms on linear scales.
- The non-convergence issue in the tracer distribution test has been discussed more fully.
- The summary and conclusions section has been rewritten in prose rather than bullet points and shortened to a more appropriate length.We hope that these revisions adequately address your concerns. We would be happy to make any further changes if needed.
Yours sincerely,
Timothy Gray, Paul Tackley, and Taras Gerya
Institute of Geophysics, ETH ZürichCitation: https://doi.org/10.5194/egusphere-2025-6354-AC3
-
AC3: 'Reply on RC1', Timothy Gray, 22 Mar 2026
-
RC2: 'Comment on egusphere-2025-6354', Thibault Duretz, 12 Feb 2026
Review of Gray et al: “Variational Stokes method applied to free surface boundaries in numerical geodynamic model”
I have read the manuscript of Gray et al. with interest. Accurately capturing planetary surfaces is a particularly relevant question in geodynamic modelling, with far-reaching implications.
This study discusses the application of a new method to treat free-surface flow in geodynamics: the variational finite-difference approach. The method is shown to provide several improvements over existing approaches, which is a strong and positive outcome of the paper.
Nevertheless, the current form of the manuscript is not ideal. It is too long, and while the focus should be on the novelty of the approach, the main message becomes somewhat diluted. At several points, the text reads more like a review article, which detracts from the emphasis on the new contributions. I have made a number of comments below that I believe could help improve the structure and clarity of the manuscript and shorten it where appropriate.
Thibault Duretz
——————————————————————
Variational form of Stokes equations can be used in various contexts, the paper title should be more specific and clearly state the variational finite differences approach.
l 32 “on on the latter”
l 40 also: fast flow of air can lead to reduction of time steps.
l 49-51 I am not aware of this problems and did not experience this in my GMG implementation (which was only 2D). This may be specific to the StagYY implementation.
l 53 kinds of confirm what I state just above.
Section 2.1 the authors maybe curious about reading this: https://api.unil.ch/iris/server/api/core/bitstreams/54c253b1-d357-455c-922a-3f3e9cd14aae/content
Section 2.1 is a review of the sticky air method. I find some text to be redundant with the intro (disadvantages of sticky air). This section could be shortened, there a lot of papers on sticky air and space should be left to set the focus on the variational FD scheme.
Section 2.2 is a review of the “staircase” approach. It coud maybe be shortened or rather being used to introduce the variational scheme. As far as I see, the variational formulation ends up with a similar discretisation but features the additional weights. These weights increases the accuracy, thereby improving the “staircase” method. Wouldn’t it be a way to frame it?
2.2.1 and 2.2.2 take a lot of space and delay the moment to which the reader reaches the core of the manuscript: the variational FD approach
2.2.2 repeats some information that is already given earlier in the test (limitation of staircase). Instead of making it sound like a review of methods for free surface, one could actually parts of the text to actually motivate the development of the variational FD scheme, where the focus is supposed to be set.
around l 200. again I am not convinced with this statements about “numerical instabilities”. I ran successful GMG models in my 2016 paper. Also, we have plenty of studies that relied on this approach (basically all papers published using MDoodz), which definitely do not feature instability cause by free surface discretisation. I think you may just have some implementation problems and I would tone down this statement.
Section 3 Variational FD Stokes method arrives too late (it’s already page 8!). Everything which is above can be condensed and rephrased in the way that motivates the development and testing of the variational FD Stokes method.
Section 3.1 Stokes equation are already stated as Eq 1. Note that there is no body force anymore in Section 3.1 for some reason.
Reading this, it feels that this section should be in fact the start of the study (right after intro) and that the above text about other (established methods) could be used for discussion purposes
Figure 3 I would gray out all the dots that are inactive for this material configuration - or are they all active
l 289 you’re discussing benefits over other approach, I would do that in the discussion.
Section 3.3 is about right hand side but the body force disappears form equation 9 onwards.
l 299 this is confusing: there is an air layer although there is no sticky air?
Eq 24 and 25 feature no equality. Maybe introduce symbols for these quantities.
section 4. I’m not sure we need all the details of the implementation within StagYY. Better focus on the results
Section 5. Results section only land on page 16. Definitely, what is above should be shortened.
around l 412 this is a lot of discussion about the sticky air method, that I barely can follow. It can be interesting in a supporting material but here the reader expects to focus on the new scheme and its benefits.
Fig 7-8-9 numerical resolution should be mentioned.
Which interpolation is being used for the sticky air method?
Fig 11-12 can be merged
Sec 5.3.2 feels over the top because the volume tracking is yet another problem (interesting too)
l 473 this conclusion is the same as the one from Tackley and King mentioned earlier? Maybe move this section a supp mat if it reproduces previous/known results?
Fig 16 feels like it does not belong to the study or should be moved to an appendix
Fig 18-19 when solid lines overlap, consider using markers/symbols instead.
For the subduction section, I think you only need Fig 21 merge with 23 (21 is too much)
l 543 you only identify here a problem that you mentioned twice in the introduction (convergence of MG solver with staircase). Therefore in the introduction, you cannot assume the problem is known to the reader.
Figure 25 is a weird figure, that appears to show 3 times the same thing. If results look all the same, one should heat for differences, not show three times the same thing. I would simply remove this figure
maybe you could remove this 3D plume benchmark section. It’s not clear why it’s sitting in between 2 scaling tests, and it’s not useful to identify differences between approaches.
Figure 27, caption: the sentence “the staircase method shows variable performance results.” is imprecise as all methods show variable performance.
To my opinion all scaling tests could be combine in one section and 1 figure.
Section 5.10 is really cool and useful application of the approach of Larionov. Somehow it ends up a bit lost between the scaling test and the discussion. I am not sure you need both 3D and 2D tests though (put the 2D in a sump mat?)
Section I would avoid using bullet points here.
6.5.1 “Applicability to incompressible flows” I think it should be “Applicability to compressible flows”
Citation: https://doi.org/10.5194/egusphere-2025-6354-RC2 -
AC4: 'Reply on RC2', Timothy Gray, 22 Mar 2026
Dear Thibault,
Thank you for your thorough and constructive review of our manuscript. We have carefully considered all comments and have made substantial revisions in response.
We have substantially restructured and shortened the manuscript from ~41 to ~26 pages. The thesis-like presentation has been replaced with a more focused paper structure. Background material on existing methods has been condensed, and the number of numerical experiments has been reduced to those that most clearly demonstrate the key properties of the method.
More specific comments have been addressed as such:
- The title has been made more specific to clearly reflect the variational finite difference approach.
- The typographic error "on on the latter" has been corrected.
- The discussion of fast air flow and its effect on timesteps has been acknowledged and the relevant text adjusted.
- Statements about numerical instabilities have been toned down to avoid overstating their generality beyond the specific StagYY implementation context.
- Sections 2.1 and 2.2 on existing methods have been shortened and reframed so that they motivate the development of the variational scheme rather than reading as a standalone review. The suggested reference has been incorporated.
- Section 3 now begins earlier in the manuscript, immediately following a condensed introduction to existing methods.
- The body force term is now included consistently from Section 3.1 onwards.
- Implementation details specific to StagYY have been reduced in Section 4, with the focus shifted to results.
- The results section has been streamlined; several figures have been merged or removed, and figures that previously showed overlapping solid lines now use markers or symbols for clarity.
- Scaling tests have been consolidated into fewer figures.
- The subduction benchmark section has been revised; the discussion of near-surface viscosity errors has been expanded and connected more explicitly to the known sensitivity documented by Schmeling et al. (2008).
- The 3D plume benchmark section has been removed.
- Bullet points in the conclusions have been replaced with prose.
- The typo in the conclusions ("incompressible" corrected to "compressible") has been fixed.We hope that these revisions address your concerns satisfactorily. We would be happy to make any further changes if needed.
Yours sincerely,
Timothy Gray, Paul Tackley, and Taras Gerya
Institute of Geophysics, ETH ZürichCitation: https://doi.org/10.5194/egusphere-2025-6354-AC4
-
AC4: 'Reply on RC2', Timothy Gray, 22 Mar 2026
Peer review completion
Journal article(s) based on this preprint
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 | |
|---|---|---|---|---|---|
| 1,817 | 1,050 | 223 | 3,090 | 191 | 338 |
- HTML: 1,817
- PDF: 1,050
- XML: 223
- Total: 3,090
- BibTeX: 191
- EndNote: 338
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Cited
Timothy Stephen Gray
Paul James Tackley
Taras Gerya
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(3525 KB) - Metadata XML
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