The tracer nudging method for correcting and preventing uneven tracer distributions in geodynamical models

Tackley, Paul James

doi:10.5194/egusphere-2025-1354

Preprints

https://doi.org/10.5194/egusphere-2025-1354

Preprints

31 Mar 2025

| 31 Mar 2025

The tracer nudging method for correcting and preventing uneven tracer distributions in geodynamical models

Paul James Tackley

Abstract. Tracers/markers/particles are commonly used in geodynamical models to track composition and sometimes temperature throughout the domain. A common problem is that over time, gaps in the tracer distribution can develop, often resulting in cells with no tracers as well as bunching of tracers. Here a correction method that perturbs or “nudges” the positions of tracers in such a way as to close gaps and eliminate bunching is presented. Test results show that this tracer nudging method is highly effective. Starting from an extremely heterogeneous tracer distribution with large regions of the domain devoid of tracers, it can produce an even distribution in only a few nudge iterations. In a time-stepping situation with a nudge every time-step, the amplitudes of the nudges are small yet sufficient to prevent gaps and bunches, allowing a low-order tracer advection method to be used while maintaining a tracer distribution that is more even than that obtained using higher-order advection methods alone. The nudge essentially corrects any non-conservation error inherent in an advection method. The computational cost is small because the method simply requires solving a Poisson equation.

Received: 21 Mar 2025 – Discussion started: 31 Mar 2025

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.

Download & links

Preprint (PDF, 1447 KB)

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 (1447 KB)

Download & links

The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.

Journal article(s) based on this preprint

17 Nov 2025

The tracer nudging method for correcting and preventing uneven tracer distributions in geodynamical models

Paul James Tackley

Geosci. Model Dev., 18, 8651–8662, https://doi.org/10.5194/gmd-18-8651-2025,https://doi.org/10.5194/gmd-18-8651-2025, 2025

Short summary

Paul James Tackley

Interactive discussion

Status: closed

RC1:
'Comment on egusphere-2025-1354', Anonymous Referee #1, 05 May 2025
This manuscript tackles an important numerical difficulty in geodynamic modelling: the development of uneven tracer distributions during Lagrangian particle advection. The authors propose a tracer-nudging algorithm, derived from the requirement that material density remain constant, which iteratively redistributes tracers until a uniform spatial density is achieved. The idea is elegant and, if widely adopted, could mitigate one of the longest-standing practical problems in high-resolution mantle-convection and lithosphere-deformation studies. While I acknowledge the novelty and potential impact of the work, several aspects of the presentation and validation need to be strengthened before the paper is suitable for publication.

Major comments
Demonstrate the method on realistic geodynamic problems

The manuscript shows simple circulation tests only. Please include at least one geologically meaningful application—e.g. a high-viscosity-contrast convection benchmark or a 2-D subduction experiment—to illustrate how tracer nudging behaves in complex, time-dependent flow fields and when and how often the nudging is needed.

2. Compare with established schemes

The paper states that the computational cost of tracer nudging is small, yet no comparison is provided. As far as I can see from the manuscript, the additional computational cost is probably higher than the existing remedies that only add a correction items to the velocity interpolation. It would be nice if the author make a comparison with other method and emphasize the advantages (and limitations) of this method
3. Clarify applicability to compressible versus incompressible flow

Lines 34–35 imply the method corrects non-divergence-free advection errors. However, many geodynamic models employ compressible Stokes flow, where ∇·v ≠ 0 by definition.

State explicitly that the current formulation targets incompressible Stokes problems.

Discuss whether, and how, the nudging algorithm could be adapted to compressible flow, and what errors might arise if it were used without modification. 4. Streamline the figures

Figures 2 and 3 convey nearly the same information. I suggest keeping the one that best illustrates the tracer-density evolution and moving the other to the supplement or removing it entirely.
Citation: https://doi.org/10.5194/egusphere-2025-1354-RC1
RC2:
'Comment on egusphere-2025-1354', Anonymous Referee #2, 19 May 2025
The study by P. Tackley presents a new algorithm for tracer advection for the marker-in-cell method. The algorithm corrects the displacement of tracers by solving for mass conservation of tracers. It is a neat, physics-based approach; however, I find that the choice of the potential function requires more discussion, and, if this algorithm is to be useful to the geodynamics community—or more broadly, to anyone using the marker-in-cell method, the writing can be improved and more tests would be useful. I detail these points below. Considering that the revision work I propose is considerable, I recommend a major revision of the manuscript.

1. Choice of the potential function for displacement vector in Eq. 3.
It is not clear why the author chose Eq. 4 (instead of Eq. 7) besides convenience to solve a linear Poisson equation. The justification is not there. I agree that it is more convenient to define dx=grad(phi)/rho_t, but what is the physical meaning of this equation? What would be the physical units of such a displacement? Maybe the author has thought about this, but the displacement should not depend on tracer density, and only Eq. 7 is a physical choice. Neither Eq. 4 nor Eq. 9 is balanced in terms of units. Has the author tried solving the non-linear Poisson equation in Eq. 8?

L77, Eqs. 7-8: Why is the non-linear equation problematic in areas where rho_t is zero? Not sure this is an issue because in the marker-in-cell method, zero tracer density in control volume is a violation, so this is avoided (which makes the models in Fig. 2 somewhat artificial).

Minor comments regarding mathematical theory:
L46, Eq. 1: Need definition (equation) for tracer density in the main text. It was obtained from the code but was not obvious. Density = #nlocal_mark * nx*ny*nz / #ntracers

Equally, for the mean density, which is taken equal to 1 here.
Mean density = nmarkx_cell*nmarkz_cell * nx*ny*nz / #ntracers = 1

L47: tracer density is "defined as the mass (or number) of tracers per unit volume". This statement is correct only if it is assumed that all tracers have the same mass, and the density relates to the number of tracers. This needs to be clarified.

Eq. 2 to Eq. 3: the discrete version is missing for v=dx/dt.

L55-57, Eq. 1: rho_c is mean density of tracers in this study. To generalise, better to use 'initial' instead of 'correct' (i.e., rho_0). Also for the purpose of this study, rho (without the t subscript) can be used for the tracer density.

2. Generalisation and stopping criterion for nudging
L111: would the equations change for variable grid spacing? Also, is the scheme generalisable for other boundary conditions, other than impermeable (L72)? Can the boundary conditions for Poisson be generalised from the velocity boundary conditions?

The number of 'nudges' required seems arbitrary and it is not clear what is the computational cost of each nudge. I wonder if a stopping criterion can be derived? Example, norm(density error)<=tolerance, where tolerance = 1e-1 or 1e-2 from Fig. 4. Is that an overkill for using the algorithm?

3. Geodynamic and performance tests
The tests presented are highly simplified, and maybe artificial (Fig. 2). To demonstrate better the impact of the algorithm, I recommend adding another test showing a geodynamic problem—with sharp velocity interfaces or rotational flows near corners. It would be interesting to see, for example, the problem cited in L27 on eruption/intrusion, which was stated, but not explained, and how this algorithm deals with it.

Then, there is the issue of the computational cost of this algorithm (with nudges). Simple Matlab scripts run fast, but the marker-in-cell method becomes expensive for 2-D and 3-D problems. What is the cost of this algorithm? How does it vary with grid size (resolution)?
L201-202: says the computational cost of solving Poisson is small, but need to demonstrate. For example, using a geodynamics test case, show time to solve for Stokes and time for tracer correction (Poisson*number_nudges).

4. Focus and style of writing
The algorithm corrects the displacement of markers by solving for conservation of tracer mass. This does not transpire from the nudging and bunching described. The writing is too focused on the algorithm (which comes out almost as a trick to deal with marker dispersion/clustering), and not so much on the physics or why is it better than other methods. I think it is not doing a great job advertising why this could be a useful technique for geodynamics, and the informal style of writing does not help (e.g. usage of regular+, Euler+nudge+random without clear definition).

Other examples,
L40: the goal shouldn't be to nudge the tracers, but to correct their displacement in order to keep uniform density.
L46: the goal is to preserve the initial tracer distribution.

The Abstract should also be revised because I had the following questions:
L9: is the correction method physics-based? A sentence describing the method is lacking. Only results are described.
L14: what does the author refer to as non-conservation errors? How do they occur? They are introduced somewhat late, even in the abstract.

The details on the content of the Matlab routines in Section 3 (L112-118) and Fig 1. should be in an appendix. Instead, the author can outline a pseudo-algorithm such as
Initialise tracer distribution, rho = rho_0

Time loop:

- advect tracers with velocity field v,
- calculate rho and density error,
- correct tracer location by solving for displacement dx (Poisson equation).

Minor comments
L6-7, 17: tracers may track other fields besides composition and temperature throughout the domain. Their advantage is that they may perform better at advecting these fields compared to grid-based methods.
L7: why does the problem occur? i.e., rotating fields in a box domain or at sharp interfaces.
L8: Suggest to use 'clustering' or 'accumulation' instead of 'bunching'.
L20: (e.g., \cite not \citep)

L21: why all the major codes use marker-in-cell? It is better to justify advantages rather than state a common practice.
L24: (e.g. 5-50) not necessary. The accuracy of the marker in cell technique relies on having a high density of tracers in each cell.
L29: why the errors occur in the first place?
L30: physics-based remedy is required
L31: one previous solution is to create new tracers, but why is this not a good remedy?

L44, L96: just because it is already implemented in StagYY, it is not a statement of validation or verification. StagYY should cite the work in the current manuscript to demonstrate the method, not the other way around.
L109: first time density bounds are given; this should occur earlier when density is defined.
L149: unclear statement. Tracers at the boundaries were nudged only a fraction of the calculated displacement to avoid crossing the boundaries.

Eqs. 11-12: could plot either the stream function or plot the velocity field to understand why the center particles might get dispersed.
Section 4.2 notation: need partial differentials (\partial).
Eq. 12: mistake Vx = 1/Lz sin() cos()
L172: define Euler+nudge+random
L174: calculated as in Eq. 12 (not above)

Figures 2 and 3 do not have colorbars or state that blue points represent tracers.
Fig 5: what is the computational cost or runtime between the cases? Add colorbars.
Fig 6: all cases should start with the same initial condition, and then additionally a case with completely random distribution.

Matlab: Add in readme how to run the script. I.e. change ntest = 1; choice between 1-11.
Citation: https://doi.org/10.5194/egusphere-2025-1354-RC2
AC1: 'Comment on egusphere-2025-1354', Paul Tackley, 23 Jun 2025

I'd like to thank the anonymous referees for their time taken to prepare RC1 and RC2.
My response is in the attached pdf. (Due to having equations and special characters it's easier to prepare it locally than in this tool).

Citation: https://doi.org/10.5194/egusphere-2025-1354-AC1

Interactive discussion

Status: closed

RC1:
'Comment on egusphere-2025-1354', Anonymous Referee #1, 05 May 2025
This manuscript tackles an important numerical difficulty in geodynamic modelling: the development of uneven tracer distributions during Lagrangian particle advection. The authors propose a tracer-nudging algorithm, derived from the requirement that material density remain constant, which iteratively redistributes tracers until a uniform spatial density is achieved. The idea is elegant and, if widely adopted, could mitigate one of the longest-standing practical problems in high-resolution mantle-convection and lithosphere-deformation studies. While I acknowledge the novelty and potential impact of the work, several aspects of the presentation and validation need to be strengthened before the paper is suitable for publication.

Major comments
Demonstrate the method on realistic geodynamic problems

The manuscript shows simple circulation tests only. Please include at least one geologically meaningful application—e.g. a high-viscosity-contrast convection benchmark or a 2-D subduction experiment—to illustrate how tracer nudging behaves in complex, time-dependent flow fields and when and how often the nudging is needed.

2. Compare with established schemes

The paper states that the computational cost of tracer nudging is small, yet no comparison is provided. As far as I can see from the manuscript, the additional computational cost is probably higher than the existing remedies that only add a correction items to the velocity interpolation. It would be nice if the author make a comparison with other method and emphasize the advantages (and limitations) of this method
3. Clarify applicability to compressible versus incompressible flow

Lines 34–35 imply the method corrects non-divergence-free advection errors. However, many geodynamic models employ compressible Stokes flow, where ∇·v ≠ 0 by definition.

State explicitly that the current formulation targets incompressible Stokes problems.

Discuss whether, and how, the nudging algorithm could be adapted to compressible flow, and what errors might arise if it were used without modification. 4. Streamline the figures

Figures 2 and 3 convey nearly the same information. I suggest keeping the one that best illustrates the tracer-density evolution and moving the other to the supplement or removing it entirely.
Citation: https://doi.org/10.5194/egusphere-2025-1354-RC1
RC2:
'Comment on egusphere-2025-1354', Anonymous Referee #2, 19 May 2025
The study by P. Tackley presents a new algorithm for tracer advection for the marker-in-cell method. The algorithm corrects the displacement of tracers by solving for mass conservation of tracers. It is a neat, physics-based approach; however, I find that the choice of the potential function requires more discussion, and, if this algorithm is to be useful to the geodynamics community—or more broadly, to anyone using the marker-in-cell method, the writing can be improved and more tests would be useful. I detail these points below. Considering that the revision work I propose is considerable, I recommend a major revision of the manuscript.

1. Choice of the potential function for displacement vector in Eq. 3.
It is not clear why the author chose Eq. 4 (instead of Eq. 7) besides convenience to solve a linear Poisson equation. The justification is not there. I agree that it is more convenient to define dx=grad(phi)/rho_t, but what is the physical meaning of this equation? What would be the physical units of such a displacement? Maybe the author has thought about this, but the displacement should not depend on tracer density, and only Eq. 7 is a physical choice. Neither Eq. 4 nor Eq. 9 is balanced in terms of units. Has the author tried solving the non-linear Poisson equation in Eq. 8?

L77, Eqs. 7-8: Why is the non-linear equation problematic in areas where rho_t is zero? Not sure this is an issue because in the marker-in-cell method, zero tracer density in control volume is a violation, so this is avoided (which makes the models in Fig. 2 somewhat artificial).

Minor comments regarding mathematical theory:
L46, Eq. 1: Need definition (equation) for tracer density in the main text. It was obtained from the code but was not obvious. Density = #nlocal_mark * nx*ny*nz / #ntracers

Equally, for the mean density, which is taken equal to 1 here.
Mean density = nmarkx_cell*nmarkz_cell * nx*ny*nz / #ntracers = 1

L47: tracer density is "defined as the mass (or number) of tracers per unit volume". This statement is correct only if it is assumed that all tracers have the same mass, and the density relates to the number of tracers. This needs to be clarified.

Eq. 2 to Eq. 3: the discrete version is missing for v=dx/dt.

L55-57, Eq. 1: rho_c is mean density of tracers in this study. To generalise, better to use 'initial' instead of 'correct' (i.e., rho_0). Also for the purpose of this study, rho (without the t subscript) can be used for the tracer density.

2. Generalisation and stopping criterion for nudging
L111: would the equations change for variable grid spacing? Also, is the scheme generalisable for other boundary conditions, other than impermeable (L72)? Can the boundary conditions for Poisson be generalised from the velocity boundary conditions?

The number of 'nudges' required seems arbitrary and it is not clear what is the computational cost of each nudge. I wonder if a stopping criterion can be derived? Example, norm(density error)<=tolerance, where tolerance = 1e-1 or 1e-2 from Fig. 4. Is that an overkill for using the algorithm?

3. Geodynamic and performance tests
The tests presented are highly simplified, and maybe artificial (Fig. 2). To demonstrate better the impact of the algorithm, I recommend adding another test showing a geodynamic problem—with sharp velocity interfaces or rotational flows near corners. It would be interesting to see, for example, the problem cited in L27 on eruption/intrusion, which was stated, but not explained, and how this algorithm deals with it.

Then, there is the issue of the computational cost of this algorithm (with nudges). Simple Matlab scripts run fast, but the marker-in-cell method becomes expensive for 2-D and 3-D problems. What is the cost of this algorithm? How does it vary with grid size (resolution)?
L201-202: says the computational cost of solving Poisson is small, but need to demonstrate. For example, using a geodynamics test case, show time to solve for Stokes and time for tracer correction (Poisson*number_nudges).

4. Focus and style of writing
The algorithm corrects the displacement of markers by solving for conservation of tracer mass. This does not transpire from the nudging and bunching described. The writing is too focused on the algorithm (which comes out almost as a trick to deal with marker dispersion/clustering), and not so much on the physics or why is it better than other methods. I think it is not doing a great job advertising why this could be a useful technique for geodynamics, and the informal style of writing does not help (e.g. usage of regular+, Euler+nudge+random without clear definition).

Other examples,
L40: the goal shouldn't be to nudge the tracers, but to correct their displacement in order to keep uniform density.
L46: the goal is to preserve the initial tracer distribution.

The Abstract should also be revised because I had the following questions:
L9: is the correction method physics-based? A sentence describing the method is lacking. Only results are described.
L14: what does the author refer to as non-conservation errors? How do they occur? They are introduced somewhat late, even in the abstract.

The details on the content of the Matlab routines in Section 3 (L112-118) and Fig 1. should be in an appendix. Instead, the author can outline a pseudo-algorithm such as
Initialise tracer distribution, rho = rho_0

Time loop:

- advect tracers with velocity field v,
- calculate rho and density error,
- correct tracer location by solving for displacement dx (Poisson equation).

Minor comments
L6-7, 17: tracers may track other fields besides composition and temperature throughout the domain. Their advantage is that they may perform better at advecting these fields compared to grid-based methods.
L7: why does the problem occur? i.e., rotating fields in a box domain or at sharp interfaces.
L8: Suggest to use 'clustering' or 'accumulation' instead of 'bunching'.
L20: (e.g., \cite not \citep)

L21: why all the major codes use marker-in-cell? It is better to justify advantages rather than state a common practice.
L24: (e.g. 5-50) not necessary. The accuracy of the marker in cell technique relies on having a high density of tracers in each cell.
L29: why the errors occur in the first place?
L30: physics-based remedy is required
L31: one previous solution is to create new tracers, but why is this not a good remedy?

L44, L96: just because it is already implemented in StagYY, it is not a statement of validation or verification. StagYY should cite the work in the current manuscript to demonstrate the method, not the other way around.
L109: first time density bounds are given; this should occur earlier when density is defined.
L149: unclear statement. Tracers at the boundaries were nudged only a fraction of the calculated displacement to avoid crossing the boundaries.

Eqs. 11-12: could plot either the stream function or plot the velocity field to understand why the center particles might get dispersed.
Section 4.2 notation: need partial differentials (\partial).
Eq. 12: mistake Vx = 1/Lz sin() cos()
L172: define Euler+nudge+random
L174: calculated as in Eq. 12 (not above)

Figures 2 and 3 do not have colorbars or state that blue points represent tracers.
Fig 5: what is the computational cost or runtime between the cases? Add colorbars.
Fig 6: all cases should start with the same initial condition, and then additionally a case with completely random distribution.

Matlab: Add in readme how to run the script. I.e. change ntest = 1; choice between 1-11.
Citation: https://doi.org/10.5194/egusphere-2025-1354-RC2
AC1: 'Comment on egusphere-2025-1354', Paul Tackley, 23 Jun 2025

I'd like to thank the anonymous referees for their time taken to prepare RC1 and RC2.
My response is in the attached pdf. (Due to having equations and special characters it's easier to prepare it locally than in this tool).

Citation: https://doi.org/10.5194/egusphere-2025-1354-AC1

Peer review completion

AR – Author's response | RR – Referee report | ED – Editor decision | EF – Editorial file upload

AR by Paul Tackley on behalf of the Authors (05 Sep 2025) Author's response Author's tracked changes Manuscript

ED: Referee Nomination & Report Request started (09 Sep 2025) by Boris Kaus

RR by Anonymous Referee #2 (24 Sep 2025)

ED: Publish subject to technical corrections (06 Oct 2025) by Boris Kaus

AR by Paul Tackley on behalf of the Authors (15 Oct 2025) Author's response Manuscript

Journal article(s) based on this preprint

17 Nov 2025

The tracer nudging method for correcting and preventing uneven tracer distributions in geodynamical models

Paul James Tackley

Geosci. Model Dev., 18, 8651–8662, https://doi.org/10.5194/gmd-18-8651-2025,https://doi.org/10.5194/gmd-18-8651-2025, 2025

Short summary

Paul James Tackley

Model code and software

Software accompanying manuscript "The tracer nudging method for correcting and preventing uneven tracer distributions in geodynamical models" ETH Zurich and Paul J. Tackley https://doi.org/10.5281/zenodo.15065273

Paul James Tackley

Viewed

Total article views: 1,690 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	BibTeX	EndNote
1,624	46	20	1,690	28	35

HTML: 1,624
PDF: 46
XML: 20
Total: 1,690
BibTeX: 28
EndNote: 35

Views and downloads (calculated since 31 Mar 2025)

Month	HTML	PDF	XML	Total
Mar 2025	21	2	0	23
Apr 2025	42	16	4	62
May 2025	34	7	5	46
Jun 2025	39	2	4	45
Jul 2025	30	3	1	34
Aug 2025	336	5	1	342
Sep 2025	1,083	2	3	1,088
Oct 2025	31	8	1	40
Nov 2025	8	1	1	10

Cumulative views and downloads (calculated since 31 Mar 2025)

Month	HTML	PDF	XML	Total
Mar 2025	21	2	0	23
Apr 2025	42	16	4	62
May 2025	34	7	5	46
Jun 2025	39	2	4	45
Jul 2025	30	3	1	34
Aug 2025	336	5	1	342
Sep 2025	1,083	2	3	1,088
Oct 2025	31	8	1	40
Nov 2025	8	1	1	10

Viewed (geographical distribution)

Total article views: 1,688 (including HTML, PDF, and XML) Thereof 1,688 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 17 Nov 2025

Download

The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.

Preprint (1447 KB)
Metadata XML

Short summary

Tracers are commonly used in geodynamical models to track composition, but a common problem is that over time, gaps in the tracer distribution can develop, as well as bunches. Here a method to correct such problems is presented and tested. The method perturbs or “nudges” the positions of tracers in such a way as to close gaps and eliminate bunching. Test results show that this tracer nudging method is highly effective. The computational cost is small.


Total:	0
HTML:	0
PDF:	0
XML:	0

The tracer nudging method for correcting and preventing uneven tracer distributions in geodynamical models

Journal article(s) based on this preprint

Interactive discussion

Interactive discussion

Peer review completion

Suggestions for revision or reasons for rejection

Journal article(s) based on this preprint

Model code and software

Viewed

Viewed (geographical distribution)