This is a well-written technical paper that analyses the numerical errors of different advection methods commonly used in geodynamic modelling and compares them to higher-order particle methods. The paper also presents two very nice new advection benchmarks for transient flow, which highlight the importance of not only benchmarking but also developing new benchmarks. A pdf with some minor comments and suggestions is attached.

General comments

I have read the manuscript of Gassmöller et al. with interest. The study investigates the accuracy of the particle in celle methods, widely used in geodynamic modelling. 

New analytical solutions are provided (and used) and I'm sure they will be very useful for the community.  

The main conclusion is that geodynamic simulators should include second-order in time particle advection solvers. 

I think the manuscript is worth publication after addressing the moderate comments listed below.

Addressing these comments should help with:

- providing a better overview of existing literature

- clarifying whether their conclusions are applicable in the general case (Main comment #1)

- demonstrating their points in a more convincing manner (Main comment #2)

section 1

Reading the first part of the introduction gives the feeling that geodynamic modelling started in the lated 2010s. I would suggest to also cite older work. Numerical geodynamic models can be traced back to the late 1970s.

l 36. you may want to oppose Eulerian to Lagrangian schemes (instead of field vs particle-based schemes). 

l 37-46 the intro would better flow if this paragraph would be better located elsewhere (e.g. section 3). As such it reads weird. One would expect the intro to focus on the topic of advection (which is introduce in the line above this paragraph and discussed right after)

l 58 yes, this is right, there is no systematic study of this type of error. Some examples are shown in Gerya and Yuen 2003 (e.g. Fig 12) and a potential solution is proposed (subgrid diffusion), the scheme is further extended to visco-elasticity in Gerya and Yuen (2007). This could be mentioned.

I would mention Weinberg & Schmeling (1992) that uses particles with RK4 and a predictor-corrector scheme

https://www.sciencedirect.com/science/article/abs/pii/0191814192901034

section 2

eq 3 I’m not sure there is a point to include a diffusion term which is not used later on. Otherwise inactive terms could also appear in other governing equations (inertia, compressibility…).

l 98 shear heating is not restricted to compressible flow. You probably want to refer to adiabatic heating instead.

section 3

I would definitely put the paragraph about Aspect history here (it’s currently in the middle of the introduction).

In the description of the each scheme (FE, RK2, RK2FOT) I would clearly indicated how many full (non-linear) mechanical solves are needed to integrate one time step.

By reading the last paragraph, it may be worth mentioning somewhere the work of Furuichi and May (2015) on implicit advection:

https://www.sciencedirect.com/science/article/pii/S0010465515000569

There is also intersecting work on time integration of non-linear Stokes problems in Adamuszek et al (2016):

https://www.sciencedirect.com/science/article/abs/pii/S0191814116300013

section 4.1

I think this both the analytic solution and its description are very useful.

section 4.2

Here it reads a bit weird in the sense that the authors describe a modification of a solution used in a previous study and the difficulties encountered. As such it reads a bit like personal notes. Maybe you could state the complete solution in Cartesian coordinates directly and simply mention that it’s different than the one used in the previous study? This would decrease the dependance of the current manuscript on this preceding study.

eq 19 is already defined in the text of section 4.1.

section 5.1

l. 308 here it’s not clear why RK2 delivers third order convergence for velocity.

l 333 not sure what the last sentence actually means.

section 5.1

MAIN COMMENT #1:

l. 365 Here you argue that the spatial accuracy of the mechanical solution and the temporal accuracy of the particle advection scheme should be on the same order of accuracy, right?

If yes, then one should opt for RK2 since the temporal resolution matches the spatial accuracy of a Q2 element.

My worry is that, in practice, model configurations generally include material boundaries with coefficient jumps (e.g. viscosity). In these cases, the order of accuracy of the mechanical solutions (L1, L2) will degrade (down to 1st order?) unless the jump is actually meshed. This last option is most often not selected in geodynamic models that rely on particles. Therefore, despite using second order stencils in space, the solutions are effectively of lower order. These cases are the most generally relevant but are not accounted by the analytical solution nor in the next section. 

So, in practice, is it still beneficial to use a second order time integration order if spatial accuracy is effectively of lower order?   

section 6

MAIN COMMENT #2:

The comparison here is interesting but this section is bit weak to my opinion.

- Is there any random noise on the initial marker position? Is it exactly the same initial condition for both models?

- The RF2FOT is first order but it should also converge to a physical solution if time step is refined. Ideally, the same pattern as for RK2, right?

- What is missing here is a temporal resolution test for the 2 schemes. It would be nice to see wether RF2FOT converges to similar patter as RK2 if time step is refined. 

In fact it’s not even clear if the RK2 is close to any physical solution here. The addition of a time convergence study would greatly improve the content of this study. 

Moreover, if the 2 schemes converge to the same physical solution, you could clearly evaluate the computing effort needed to reach the “converged” solution pattern with both methods. This could clearly show the benefit of RK2 on a practical case.

Ideally, an application that includes material jump would have been ideal (e.g. D’’, magma body, lithosphere stratification) but I realise it might be out of the scope of the current study. 

p 17 Physical units need not to be italicised.