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.

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”