I find the manuscript entitled “Modelling chemical advection during magma ascent” by Dominguez et al. quite interesting and helpful for the geoscience community. In my view, the problem of numerical advection is often underappreciated, and the presented study clearly shows the pitfalls of using too simplistic advection schemes. In their work, the authors introduce and compare the performance of four selected numerical advection schemes. In conclusion, they advocate the usage of the WENO type of schemes for petrological applications (chemical advection). Below I give my general and detailed suggestions that the authors may consider to follow before the final publication of the manuscript.

The numerical schemes selected for the study are introduced as simplified 1D variants and later tested using 2D setups. Although the 2D generalizations are rather straight-forward, in my view, it would be worth to give some clues to the readers on how they are developed. In particular, the 2D generalization of the WENO scheme that builds on a 5-point stencil in 1D is perhaps least obvious in this respect (not to mention the treatment of the boundary conditions). In fact, in the final parts of the manuscript, the authors themselves mention an alternative formulation for the 2D WENO scheme that builds on a 3x3 stencil.

I was a bit confused about the MIC scheme that was used in the study. The authors mention the initial step of interpolating the field values from grid nodes to markers. Considerable focus was given to integrating the trajectories of the markers. The step of interpolating the field values back from markers to nodes was sketched out. Unless I got it wrong, the advected fields in this study are essentially phase fields, so no calculations and updates need to be performed on the Eulerian grid. Thus, there is no need to interpolate back these updates to markers in such cases. However, interpolating the base field values from grid nodes to markers would be an undesired and suboptimal step leading to unnecessary numerical diffusion. I think that it would be helpful to clarify the design of the MIC scheme that was used in the tests. What is the reason for the observed mass loss in the second test using the MIC scheme? Could it be related to the non-conservative reseeding and removal of the markers? I would suspect that the mass conservation of the standard MIC should be almost ideal.

The author presented the rigid-body rotation test as their base test. In my view, a bit more complex synthetic 2D test such as for example the shear cell test would be a great addition to the study. On a different note, I actually wonder how to design a balanced comparison of the accuracy of grid-based advection schemes against marker-based methods, when their numerical resolutions largely differ (for a fixed grid resolution, a much larger number of markers is used than nodes in this study). In addition, the authors took advantage of the extended stability of the MIC and SL schemes and used a larger value of the time step than in the case of the upwind and WENO schemes, however, this may have an effect on the numerical performance of the studied schemes.

I was confused by the fact that it was not possible to obtain MIC results using a 128 GB machine for grid resolutions such as 500x1000, when the number of markers is on the order of 10 million. What was the exact reason for such limit in the current MIC implementation? The final conclusions build on the statement that the MIC scheme lacks mass conservation and it is expensive in terms of computation and memory. I already expressed my concerns about the observed lack of mass conservation in the case of the MIC scheme. I think that it would be also useful if the authors could elaborate more on the computational and memory performance of their MIC implementation. 

Marcin Dabrowski

Detailed comments:

38 What is exactly the notion of triviality in this case?

43 “This brings limitation to the resolution of the model …” - It is indeed correct that for any given computational resources 1D numerical models can be studied with a higher resolution than 2D models. On the other hand, looking from a 3D physical perspective using approximate 2D numerical models could be considered as relaxing computational limitations to the resolution of the model.

eq.5 Is it necessary to use the brackets around C_e^f?

84 Perhaps early rather than earlier?

87 “It consists at tracking individual particles on a Lagrangian frame and to reinterpolate them when needed on an Eulerian stationary mesh grid”. consist at or consist of? to reinterpolate or reinterpolating?

94 “..to the position of a fixed Eulerian grid..” It would be perhaps worth to mention grid nodes here.

106 Please change spacial to spatial.

106 It might be a bit confusing when the term “element” is mentioned here.

eq. 6 It would be perhaps more transparent if the f symbol, which stands for fluid, was used in the superscript rather than subscript. It could be mentioned that v denotes the x component of the velocity.

109 Please consider using “grid spacing” rather than “increment in space”.

eq. 9 Please consider using the f symbol in the superscript. The j index that denotes the spatial component of the velocity could be then used in the subscript.

129 Here and in fig. 1 caption, “Upwind” is capitalized, while it is not in l. 111.

129 Please consider removing the comma.

133 Perhaps “for a single grid node” rather than “1 element” would suit better here. This issue reappears in other parts of the manuscript.

153-54: What is meant by “complex problems” here? I suspect that for any problem boundary conditions need a special or careful treatment when a 5-point stencil is used. I don’t think that it is necessary to present in details how BC are treated in this case, but a short description would be useful.

158 It is not that important, but I’m wondering about the notion of “fully capture” in this case.

164 Please replace spacial by spatial.

171 What is “therefore” referring to?

182 “ … is not easy to determine xd. A common approach to overcome this limit…” In my vewi, the term limit may be misleading in this context. I would suggest something like: A common approach to accurately determine …

193 “a stencil point on the grid” Perhaps it could be described as “a grid node”.

196 It could be several (neighboring) cells depending on the interpolation order (here cubic, so four nodes are needed to span the interpolation space). On the other hand, it depends on the adopted definition of the cell. Anyway, I think that it would be worth to clarify this issue.

214 “…interpolate their values…”. Perhaps “…interpolate field values associated with them…”

220 “The initial value of each marker …“ Perhaps “The initial field value in each marker …”

222-23 The problem of deteriorating marker resolution is not only limited to highly divergent flows. It may also affect incompressible, but strongly stretching flows. In such cases, the density of markers can actually depend on the direction (markers may be locally jammed in a certain direction and rarefied in the direction that is perpendicular). Thus, simple scalar measures of marker density may fail in such cases.

246-54 In my view, based on the presented description, it is not exactly clear how the field values are interpolated from markers to grid nodes. Is it the least-square type of linear interpolation? How many surrounding cells are used?

256-261 The domain size is given as dimensionless, while the time is in seconds.

268 Isn’t S_k constant for all the internal cells?

352-354 Is this the description of the ODE solver implemented in DifferentialEquations.jl (l. 348)? Please clarify.

363 The standard deviation is given as a dimensionless quantity, while the size of the model domain is given in meters.

366-67 I might be missing something here, but I don’t exactly understand how the initial distribution of the melt chemistry is governed. In fig. 8 it appears just bi-modal at t=0 Ma.

Table 3 Please fix the value of the shear modulus (an incorrect exponent is used).

374-383 This part might be a bit confusing. Although the link to eq. 5 is mentioned in l. 320, I think that it would be useful to clarify which physical fields are advected in this model. Please explain why the oxide content may not be conserved, despite the melt fraction is claimed to be conserved due to the chosen numerical scheme. If the total mass of the melt composition (summed weight percentages?) is renormalized in each time step, then I’m not sure if it could be treated as a conserved quantity.

389 What is meant by mixture here?

397 I am a bit confused by the observed lack of mass conservation in the upwind scheme. Is it due to the fact that the renormalization procedure is used? I would also expect MIC to be essentially mass conserving and non-diffusive. Are the observed mass loss effects due to marker reseeding and removal?

407 Using mid-point for computing trajectories in MIC and SL results in some iterations.

423 What is exactly meant by complex boundary conditions?

Fig. 9 & 10 The axes are logarithmic and there is no need to add log(.) in the labels.

“Modelling chemical advection during magma ascent” by Dominguez et al. is an interesting manuscript that compares four different numerical schemes for phase field advection, emphasizing their suitability within the context of two-phase flow models. ​I believe the manuscript is valuable for the earth sciences modeling community, as it assesses the behaviour of these advection schemes in terms of accuracy -here mass conservation- and computational cost. 

The upwind and SL schemes are quickly disregarded as they introduce excessive numerical diffusion and have a low accuracy. The authors find WENO as the most suitable advection scheme as it conserves mass at a decent computational cost, arguing that MIC does not conserve mass efficiently well in one of their two numerical experiments and is far too expensive in terms of performance and memory consumption.

Below are some general comments regarding the manuscript, followed by detailed line by line comments.

Albert de Montserrat

General comments:

(1) The authors should provide more details about the interpolants used in both the Semi-Lagrangian and MIC schemes. They also mention that complex boundary conditions may require some extra care in the WENO scheme. It is not clear to me what these complex boundary conditions are. Do commonly used boundary conditions such as free-slip or no-slip require some special treatment? Since the periodic boundary conditions used in the examples given in the manuscript are not the only common ones, a more detailed discussion would be helpful.

(2) The rotating circle benchmark shows that WENO and MIC produce similar results in terms of mass conservation. However, WENO performs better in conserving mass in the two-phase flow  example. It's important to note that in the first benchmark, WENO SL, and MIC are run with the same constant time step. In the other example, the Courant number for MIC was more than twice that of WENO's. To ensure a fair comparison of how well these schemes conserve mass, I suggest running the two-phase flow example with the same time step for WENO, MIC and SL. This will make the comparison of mass conservation fairer. If the mass conservation of MIC and SL improves by using the same time step as WENO, it would be helpful to run both MIC and SL with increasing time steps to determine when their accuracy starts to deteriorate.

(3) The manuscript concludes that, besides better mass conservation, the WENO scheme is the best candidate because of the poor performance of MIC and its excessive memory consumption. However, the comment stating that the MIC does not run at high resolution because it does not fit in 128 GB of RAM is surprising and unclear. Let’s assume that one allows a maximum of 30 particles per cell, then in the higher resolution case you have an upper bound of 15e6 particles (500 x 1000 x 30). Further assuming that you are using FP64, this means that only 120MBs (in Julia: sizeof(Float64)*15e6 / 1e6 ) are needed to store a particle field. With this in mind, there is plenty of RAM available to store multiple fields in the particles, and even for a much larger number of particles. Last, all the necessary arrays can be pre-allocated so that all the MIC kernels (advection, interpolation, injection….) can be performed in-place without requiring additional memory

Inspecting your MIC algorithm with packages such as Cthulhu.jl and JET.jl will reveal type instabilities in the code, which result in dynamic dispatch and potentially generate a large number of unnecessary allocations. Available profiling tools can be used to track down where the remaining allocations come from.

The main point of this paper is not addressing the optimization of the implementation of any of the advection schemes here present. However, since performance is one of the criterions chosen by the authors to select the best candidate scheme, I would suggest that the authors consider resolving the main performance issues in their codes. By doing so, the MIC‘s performance is also likely to improve, enabling it to run the two-phase flow example at higher resolutions. This would allow for a better comparison between MIC and WENO, and potentially turn MIC into a more viable candidate.

(4) Attached is a PDF with a series of typo, grammar and style corrections and suggestions.

Line by line comments

L5: “weighted essentially non-oscillatory”. Throughout the manuscript WENO is spelled out with either lower or upper case initials (e.g. line 125). For consistency, I recommend to use either one or the the other throughout the entire manuscript

L87 “It consists at” -> “It consists of”

L88 “to interpolate” -> “interpolating”

L88 “mesh grid”. I would call it either mesh or grid.

L89 “to be” -> “being”

L133 “element” what is the element here?

L147 You could consider using the machine precision epsilon (function eps() in Julia) instead of fixing epsilon=1e-6

L153 “Careful consideration must then be given to boundary conditions for complex problems.”  What kind of considerations have to be made for the boundary conditions? It would be helpful to elaborate on this. I don’t think it is obvious.

L166, Isn’t the time integration third order in space and first order in time? If I’m reading correctly eqs. 11-13, C is evaluated at the same time step in all the integration stages.

L176 “rectilinear grid” perhaps “uniform” or “regular” are more appropriate

L176 ”this reduces the complexity of the implementation and the numerical cost of the interpolation function.” Could you briefly describe why this is so (finding the parent cell, cheap mapping to the reference cell, etc…) ? It may not be obvious for a reader who never tried to implement any of these advection schemes.

L177 “From a particle point of view, the goal is to find the starting points at the previous timestep for each grid point.” Perhaps this could benefit from some rephrasing, what you want to do is to find where a particle that is at a grid point (i,j) at time t^n+1 was at the time t^n

L183 What about higher order methods such as Runge-Kutta 4 ? would they work better to back track the particles?

L190 How is vf interpolated to t^n+½? is it just the average between t^n and t^n+1? I guess this requires storing the velocity at two time-steps, this should be mentioned somewhere in the manuscript.

L197 Are bi/tri-linear interpolations not good enough?

L205 “clipping”, I think “clamping” is a better word for this.

L205-210 The manuscript lacks the definition of C^H. It would be helpful to include this information.

L214 “interpolate their values on a fixed grid.” Interpolation often goes in both directions (e.g. temperature), not only from the markers onto the grid.

L220 Are the initial positions of the markers randomized? If they are, it should be mentioned later on the benchmarks descriptions.

L220 “The initial value of each marker can then be obtained by linear interpolation from the initial conditions of the Eulerian grid.” Most of the time, a field can be initialized directly at the particles, which provides a more accurate starting point

L223 “a non-conservative strategy”

L246-147 “This is more complex than for SL schemes because the markers are not uniformly distributed for a non-trivial velocity field.” Another disadvantage is that the  parallelizing this interpolation is susceptible to race conditions in shared memory systems, requiring the use of atomic operations

L252-253 Could you add eq. of the linear interpolant?

L263 How is the computational time of one time step measured? Is it the time of the first step? If so, does it include compilation time? Or is it the average step time throughout the whole model?

As the number of particles varies throughout the model run time, it may be better to average the time of each time step (excluding the first one to avoid measuring compilation time, or doing a warm up run) and add a confidence interval.

Would be helpful to clarify what is actually done and measured in the MIC case. I believe that you are interpolating from particles to grid (plus reseeding?) in each time step. This is technically not needed in this case, since you are merely advecting passive markers. It may be helpful to include in the manuscript a breakdown of the percentage of time spent on all the MIC stages done in the benchmark.

Table 2. 

• Here you call it “SL QM”, however, it’s written as “QM SL” in some other places. Please choose one naming convention for the whole manuscript.
• I get very different timings for SL and MIC (Julia 1.10, Windows 11, CPU: 24 × AMD Ryzen 9 5900X 12-Core Processor). Different Julia versions are known to have different performances, it would be helpful to add to the caption what Julia version was used. Later on you mention what CPU you used, but I think it would be nice to have this information also in this caption.

Section 4.2 Perhaps it is better to define the scaling lengths in a small table instead of spelling them out line by line.

Section 4.3 The time step of QMSL and MIC is roughly ~x2 with respect to WENO and Upwind. It may be beneficial to run QMSL and MIC with the same Courant number as in the other two cases and include the mass conservation results to Fig 9.. This could potentially reduce the reseeding in the MIC case and improve the mass conservation.

L349 “This method” -> “This package”

L354 As a curiosity, what kind of automatic differentiation is being used?

L360 “assert” means to confirm something. Perhaps “examine”, “test”, or similar, work better here.

L365 It is stated later that all the runs are running serially in a single core, perhaps it’s useful to mention that fact here instead/as well.

L399 “that is linked on how often the two-phase flow solver is called” perhaps you could rephrase it so that it is more clear that being able to increase the time step results in fewer time steps being required.

L400 “...QMSL shows the best scaling…” What do you mean by scaling here? Scaling is usually referred to as how well one algorithm scales with the number of cores/CPUs. However, all the models were run serially in a single core.

L409 “...has no stability condition” However, it is common to limit time steps so that a particle does not move more than one cell  (or a fraction of it). This reduces significantly the amount of effort to find the new parent cell of a given particle, in particular in the case of non-regular grids.

L423 “If complex boundary conditions are required” As in my previous comment, I think it is not obvious what these complex boundary conditions are.