the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Benchmarking the accuracy of higherorder particle methods in geodynamic models of transient flow
Juliane Dannberg
Wolfgang Bangerth
Elbridge Gerry Puckett
Cedric Thieulot
Numerical models are a powerful tool for investigating the dynamic processes in the interior of the Earth and other planets, but the reliability and predictive power of these discretized models depends on the numerical method as well as an accurate representation of material properties in space and time. In the specific context of geodynamic models, particle methods have been applied extensively because of their suitability for advectiondominated processes and have been used in applications such as tracking the composition of solid rock and melt in the Earth's mantle, fluids in lithospheric and crustalscale models, light elements in the liquid core, and deformation properties like accumulated finite strain or mineral grain size, along with many applications outside the Earth sciences.
There have been significant benchmarking efforts to measure the accuracy and convergence behavior of particle methods, but these efforts have largely been limited to instantaneous solutions, or timedependent models without analytical solutions. As a consequence, there is little understanding about the interplay of particle advection errors and errors introduced in the solution of the underlying transient, nonlinear flow equations. To address these limitations, we present two new dynamic benchmarks for transient Stokes flow with analytical solutions that allow us to quantify the accuracy of various advection methods in nonlinear flow. We use these benchmarks to measure the accuracy of our particle algorithm as implemented in the ASPECT geodynamic modeling software against commonly employed field methods and analytical solutions. In particular, we quantify if an algorithm that is higherorder accurate in time will allow for better overall model accuracy and verify that our algorithm reaches its intended optimal convergence rate. We then document that the observed increased accuracy of higherorder algorithms matters for geodynamic applications with an example of modeling smallscale convection underneath an oceanic plate and show that the predicted place and time of onset of smallscale convection depends significantly on the chosen particle advection method.
Descriptions and implementations of our benchmarks are openly available and can be used to verify other advection algorithms. The availability of accurate, scalable, and efficient particle methods as part of the widely used opensource code ASPECT will allow geodynamicists to investigate complex timedependent geodynamic processes such as elastic deformation, anisotropic fabric development, melt generation and migration, and grain damage.
 Article
(4491 KB)  Fulltext XML
 BibTeX
 EndNote
Numerical models have been a key tool for geoscientists investigating the processes governing plate tectonics and mantle convection. Among the many one could cite, a crosssection of publications include studies of the evolution of mantle heterogeneities over time (e.g., Kellogg and Turcotte, 1990; McNamara and Zhong, 2005; Brandenburg et al., 2008; Gülcher et al., 2021; Jones et al., 2021), the initiation and evolution of plate boundaries (e.g., Tackley, 1998; Bercovici and Ricard, 2014; Baes et al., 2020; Schierjott et al., 2020), the fate of subducted slabs (e.g., Gurnis and Hager, 1988; Billen, 2008; Faccenna et al., 2017; Grima et al., 2020), plume dynamics (e.g., Farnetani and Richards, 1994; Lin and van Keken, 2006; Dannberg and Gassmöller, 2018; Arnould et al., 2020), the dynamics of microplates (e.g., Glerum et al., 2020; Neuharth et al., 2021), and the seismic cycle (e.g., Van Dinther et al., 2013; Van Zelst et al., 2019). Obviously, the usefulness of such dynamic models relies on the accurate approximation of solutions of the equations that describe the processes under consideration. For geodynamic models of the solid Earth, this usually requires solving the Stokes equations governing the flow, and advection(diffusion) equations governing the transport of thermodynamic properties like temperature or entropy, chemical composition, and trace elements, as well as deformation properties like damage, or mineralogical properties like grain size. Established methods for solving the Stokes equations typically treat the fluid as a continuum and are based on the finiteelement (e.g., Moresi et al., 2022), finitedifference (e.g., Gerya and Yuen, 2003; Kaus et al., 2016), and finitevolume methods (e.g., Tackley, 2008). In contrast, there is a wide variety of methods for solving the advection equations (van Keken et al., 1997; Puckett et al., 2018) such as particle methods, continuous or discontinuous field methods, markerchain methods, or volumeoffluid methods.
Due to their inherent suitability for modeling advectiondominated problems, different variants of particle methods have become popular in the geodynamic modeling community (Weinberg and Schmeling, 1992; van Keken et al., 1997; Tackley and King, 2003; Gerya and Yuen, 2003; McNamara and Zhong, 2004; Gassmöller et al., 2018; Samuel, 2018; Sime et al., 2021). The main advantage of particles in geodynamic applications is that particles advected with the material flow keep their associated material properties; that is, these properties do not diffuse in space as is the case for many fieldbased methods. It also means that the differential equations for each particle's location are simply ordinary differential equations for which many good solution approaches are available. On the other hand, while errors in particle methods are less apparent than for field methods, they still exist (Tackley and King, 2003; Gassmöller et al., 2019). In particular, previous work has discussed the influence of errors due to interpolating properties from particles to finiteelement functions representing Stokes discretizations (Thielmann et al., 2014), the influence of the divergence of the computed velocity field on particle distributions (Wang et al., 2015; Pusok et al., 2017; Sime et al., 2021), and the advection of particles over time in spatially variable flow (Gassmöller et al., 2019). However, a source of error in particle advection methods that has, to the best of our knowledge, not been systematically discussed is the error in advecting the particle position in transient, rapidly changing flow. (Some examples of this can be seen in Gerya and Yuen, 2003.) This type of flow is common in geodynamic models of the upper mantle or lithosphere, because employing a viscoplastic or stressdependent rheology can cause strong nonlinear feedbacks between the current solution and material properties, and therefore can cause fast changes over time. While the presence of these errors is known, only few studies systematically investigate its influence on geodynamic applications (Trim et al., 2023a, b). This is largely due to the difficulty in quantifying their influence, as one needs a timedependent model solution against which to compare numerical results, and most currently available benchmarks either rely on instantaneous solutions (Duretz et al., 2011; Zhong, 1996; Zhong et al., 2008; Schmid and Podladchikov, 2003; Kramer et al., 2021), a steadystate solution (Zhong et al., 2008; Gassmöller et al., 2019), or a comparison between several numerical methods without known exact solutions (van Keken et al., 1997; Tackley and King, 2003).
In this work, we measure the particle advection error in transient flow using the particle architecture we have developed as part of our work on the Advanced Solver for Planetary Evolution, Convection, and Tectonics (ASPECT; Kronbichler et al., 2012; Heister et al., 2017). We start with a description of the mathematical problem we would like to solve in Sect. 2 and then present an analysis of the numerical errors that result from the advection of particles in transient flow (Sect. 3). We develop new benchmarks for transient flow in a box and spherical shell that have known analytical solutions (Sect. 4) and use these benchmarks to measure the accuracy of the discussed particle advection methods and quantify their influence on the results of the Stokes equations (Sect. 5). Finally, we illustrate why focusing on the accuracy of particle methods matters for practical geodynamic applications with a model example of smallscale convection developing underneath oceanic lithosphere (Sect. 6). We conclude in Sect. 7. Appendix A contains the derivation of the analytical solution for the spherical shell benchmark and Appendix B contains a more detailed discussion of some benchmark results.
For the models in this work, we will consider the incompressible Stokes equations using the Boussinesq approximation. They consist of a force balance and a mass continuity equation:
where bold script represents vector quantities, $\mathit{u}=\mathit{u}(\mathit{x},t)$ is the velocity, $p=p(\mathit{x},t)$ the pressure, ρ the density, η the viscosity, and g the gravity. $\dot{\mathit{\epsilon}}\left(\mathit{u}\right)=\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{\nabla}\mathit{u}+\mathrm{\nabla}{\mathit{u}}^{\mathrm{T}})$ denotes the strain rate tensor. In other words, because inertia is absent, the equations above describe an instantaneous equilibrium.
The models evolve in time because the density and viscosity may depend on time through additional variables ${\mathit{\varphi}}_{c}={\mathit{\varphi}}_{c}(\mathit{x},t)$, $c=\mathrm{1},\mathrm{\dots},{N}_{c}$, where N_{c} is the number of additional quantities. Examples include the temperature or other thermodynamic quantities, or chemical compositions. These quantities typically satisfy advection–diffusion equations and may be solved through either field or particlebased approaches. For the purposes of this paper, let us specifically focus on those properties with negligible diffusion (such as chemical compositions or grain sizes); particle methods for applications with nonnegligible diffusivity and reactions have been described elsewhere (Gerya and Yuen, 2003; Sime et al., 2022). In other words, we consider equations of the form
where the H_{c} are source terms. In Eqs. (1)–(3), material properties η and ρ, as well as source terms H_{c}, are then assumed to depend (perhaps nonlinearly) on the solution variables u, p, and ϕ_{c}, resulting in a coupled and timedependent system of equations.
We end this section by noting that while for simplicity we use the incompressible Stokes equations, the usefulness of the benchmark models we present below do not rely on this assumption and will be transferable to compressible models. In fact, an accurate solution to the advection equation may matter more in compressible models, because they often contain more coupled terms such as adiabatic heating (depending on the pressure gradient), the pressure dependence of the density, and additional processes like phase transitions caused by pressure or temperature changes.
Over the past years we have developed a flexible, scalable, and efficient particle architecture (Gassmöller et al., 2018). This work is open source, and performs well in modern highperformance computing environments. In particular, it supports advanced computational methods such as an adaptively refined, unstructured, and dynamically changing background mesh, parallelization beyond tens of thousands of (CPU) compute cores, storing of arbitrary particle properties, and complex nonlinear solvers. The underlying particle infrastructure is application agnostic and independent of the used discretization for the fieldbased quantities. Its methods are integrated into the general purpose opensource finiteelement software library deal.II (Arndt et al., 2023) and have been used to model a wide range of geoscientific applications as well as Navier–Stokes flow, mixing of granular materials, solid–fluid interaction, and laser metal deposition of metallic particles (Popov and Marchevsky, 2022; Arndt et al., 2020; Golshan et al., 2022; El Geitani et al., 2023; Golshan and Blais, 2022; Murer et al., 2022).
We have discussed the numerical methods for most steps of our particle algorithm (Gassmöller et al., 2018, 2019) and Stokes solver (Kronbichler et al., 2012; Heister et al., 2017) in earlier work and refer the reader there for details on the finiteelement method, timestepping algorithm, particle generation, advection, and interpolation from particles to grid. Here we will extend this earlier work by developing transient solutions and focus on how the temporal accuracy of advection methods controls the overall accuracy of a coupled geodynamic model.
3.1 Particle advection
In particle methods, the values of fields ϕ_{c} are approximated by advecting particles that carry these field values as “properties”. Particles move with the velocity u(x,t) that results from solving the Stokes Eqs. (1)and (2), and the properties carried by a particle evolve based on H_{c} in Eq. (3). In other words, the solution of the partial differential Eq. (3) is approximated by solving an ordinary differential equation (ODE) tracking the position x_{i}=x_{i}(t) for each particle i and a separate ODE tracking the evolution of the properties carried by the particle:
In practice, the exact velocity u(x,t) is not available but only a numerical approximation in space u_{h}(x,t) to u(x,t). Furthermore, this approximation is only available at discrete time steps, ${\mathit{u}}_{h}^{n}\left(\mathit{x}\right)={\mathit{u}}_{h}(\mathit{x},{t}^{n})$ and it needs to be interpolated between time steps if the advection algorithm for integrating Eq. (4) requires one or more evaluations at intermediate times between t^{n} and t^{n+1}. If we denote this interpolation in time by ${\stackrel{\mathrm{\u0303}}{\mathit{u}}}_{h}(\mathit{x},t)$ where ${\stackrel{\mathrm{\u0303}}{\mathit{u}}}_{h}(\mathit{x},{t}^{n})={\mathit{u}}_{h}^{n}\left(\mathit{x}\right)$, then the equation being solved is actually
where ${\stackrel{\mathrm{\u0303}}{\mathit{x}}}_{i}\left(t\right)$ is the exact solution of this equation using the “wrong” velocity field. If ${\stackrel{\mathrm{\u0303}}{\mathit{u}}}_{h}$ is a good approximation to u, then we hope that $\stackrel{\mathrm{\u0303}}{\mathit{x}}\left(t\right)$ is a good approximation of x(t). In practice, however, we can not even compute $\stackrel{\mathrm{\u0303}}{\mathit{x}}\left(t\right)$ but need to further approximate it via time stepping.
3.2 Convergence of particle advection methods
The particle positions contain error contributions from the inexactly known velocity field discussed in the previous subsection, as well as the error introduced by time stepping the ODEs describing particle position and properties. If we denote by ${\stackrel{\mathrm{\u0303}}{\mathit{x}}}_{i,h}\left(t\right)$ the numerical approximation to the solution of Eq. (6), then the error at some time t will typically satisfy a relationship like
where Δt_{p} is the time step used by the ODE solver, which is often an integer fraction of the time step Δt_{u} used to advance the velocity field u. In our application we will choose Δt_{p}=Δt_{u}; q is the convergence order of the method and C(t) is a (generally unknown) constant that depends on the end time t at which one compares the solutions as well as on $\stackrel{\mathrm{\u0303}}{\mathit{u}}$. We want to compare this computed solution against exact trajectories using the exact velocity as in Eq. (4) and then assess the error as $\Vert {\mathit{x}}_{i}\left(t\right){\stackrel{\mathrm{\u0303}}{\mathit{x}}}_{i,h}\left(t\right)\Vert $. This quantity will, in the best case, only satisfy an estimate of the form
with appropriately chosen norms for the second and third terms, which represent how accurately the flow field is discretized in space and time. All of these terms can converge to zero at different rates with the mesh size h and the timestep size Δt. As a consequence, each of these terms may be the limiting factor for the overall accuracy of the ODE integrator.
3.3 Common particle integrators
Given these considerations and that ODE integrators require the expensive step of evaluating the velocity field ${\stackrel{\mathrm{\u0303}}{\mathit{u}}}_{h}$ at arbitrary points in time and space, choosing a simpler, less accurate scheme can significantly reduce the computation time. In our work, we have implemented the Forward Euler, Runge–Kutta 2 (RK2) and Runge–Kutta 4 (RK4) schemes (Hairer and Wanner, 1991; Gassmöller et al., 2018), although other methods are available and have been used in geodynamic applications (e.g., Heun's method in Zhong and Hager (2003) and Sime et al. (2021), Runge–Kutta schemes with additional predictor–corrector steps in Weinberg and Schmeling (1992), implicit Euler and BDF2 methods in Furuichi and May (2015), or Adams Bashforth methods in Adamuszek et al. (2016)). We will briefly discuss our selected methods below and will limit ourselves to a discussion of two different variants of the RK2 integrator, which is sufficient to support our conclusions.
For simplicity, we will omit the particle index i from formulas in the remainder of this section and will assume that the ODE and partial differential equation (PDE) time steps Δt_{p} and Δt_{u} are equal. We will therefore simply denote them as Δt. This is often the case in practice because the velocity field is typically computed with a method that requires a Courant–Friedrichs–Lewy (CFL) number around or smaller than one, implying that particles move no more than by one cell diameter per (PDE) time step. In such cases, even explicit time integrators for particle trajectories can be used without leading to instabilities, and all of the methods below fall in this category. The formulas in the remainder of this section are, however, obvious to generalize to cases where Δt_{p}<Δt_{u}. We will also assume that we have already solved the velocity field up to time t^{n+1} and are now updating particle locations from x^{n} to x^{n+1}. In cases where one wants to solve for particle locations before updating the velocity field, ${\stackrel{\mathrm{\u0303}}{\mathit{u}}}_{h}$ can be extrapolated beyond t^{n} from previous time steps, or particle advection and velocity computation could be iterated in a nonlinear solver scheme. Because of this choice, the number of Stokes solutions which have to be computed is independent of the choice of particle advection scheme.
In the following, we briefly describe some of the common time stepping algorithms, including those we use in this work:

Forward Euler (FE) – the simplest method often used is the forward Euler scheme,
$$}{\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n+\mathrm{1}}={\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n}+\mathrm{\Delta}t\phantom{\rule{0.25em}{0ex}}{\stackrel{\mathrm{\u0303}}{\mathit{u}}}_{h}({t}^{n},{\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n}).$$It is only of first order (that is, the exponent in Eq. (7) is q=1) but cheap to evaluate because it only requires evaluating the velocity solution at an alreadycomputed time point.

Runge–Kutta second order (RK2) – accuracy and stability can be improved by using a secondorder Runge–Kutta scheme (that is, q=2 in Eq. (7)). Among the many secondorder Runge–Kutta methods, we specifically choose what is commonly referred to as the (explicit) midpoint method in which the new particle position is computed as
$$\begin{array}{ll}{\displaystyle}{\mathit{k}}_{\mathrm{1}}& {\displaystyle}={\displaystyle \frac{\mathrm{\Delta}t}{\mathrm{2}}}{\stackrel{\mathrm{\u0303}}{\mathit{u}}}_{h}({t}^{n},{\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n}),\\ {\displaystyle}{\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n+\mathrm{1}}& {\displaystyle}={\mathit{x}}^{n}+\mathrm{\Delta}t\phantom{\rule{0.25em}{0ex}}{\stackrel{\mathrm{\u0303}}{\mathit{u}}}_{h}\left({t}^{n}+{\displaystyle \frac{\mathrm{\Delta}t}{\mathrm{2}}},{\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n}+{\displaystyle \frac{{\mathit{k}}_{\mathrm{1}}}{\mathrm{2}}}\right).\end{array}$$This method requires evaluating the computed velocity at two different locations and two different points in time, including a time point intermediate between (velocity) time steps. Note that there are other RK2 methods, such as Ralston's method, which reduce the theoretical truncation error of the method while maintaining the order of convergence. However, in our benchmarks the difference in error is small, and the midpoint method allows us to reduce the memory requirement of the algorithm.

Runge–Kutta secondorder space, firstorder time (RK2FOT) – in practice, many geodynamic modeling packages only store a single velocity solution at a time, which prevents the interpolation of the velocity field at ${t}^{n}+\frac{\mathrm{\Delta}t}{\mathrm{2}}$ used in RK2 from adjacent solutions at t^{n} and t^{n+1}. However, reasonable accuracy can still often be achieved when ignoring the time dependence of the velocity (Gerya and Yuen, 2003; McNamara and Zhong, 2004). Here we implement such an advection scheme as a modification to RK2, in which the new particle position is computed as
$$\begin{array}{ll}{\displaystyle}{\mathit{k}}_{\mathrm{1}}& {\displaystyle}={\displaystyle \frac{\mathrm{\Delta}t}{\mathrm{2}}}{\stackrel{\mathrm{\u0303}}{\mathit{u}}}_{h}({t}^{n},{\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n}),\\ {\displaystyle}{\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n+\mathrm{1}}& {\displaystyle}={\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n}+\mathrm{\Delta}t\phantom{\rule{0.25em}{0ex}}{\stackrel{\mathrm{\u0303}}{\mathit{u}}}_{h}\left({t}^{n},{\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n}+{\displaystyle \frac{{\mathit{k}}_{\mathrm{1}}}{\mathrm{2}}}\right).\end{array}$$Note how, compared with the RK2 scheme, the velocity is evaluated at the (wrong) time t^{n} instead of ${t}^{n}+\frac{\mathrm{\Delta}t}{\mathrm{2}}$ but at the correct location ${\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n}+\frac{{\mathit{k}}_{\mathrm{1}}}{\mathrm{2}}$. RK2FOT still requires the evaluation of the velocity at two different points in space but only a single point in time.

Runge–Kutta fourth order (RK4) – a further improvement in particle advection can be achieved by a fourthorder Runge–Kutta scheme. We choose the most commonly used scheme that computes the new position as
$$\begin{array}{ll}{\displaystyle}{\mathit{k}}_{\mathrm{1}}& {\displaystyle}=\mathrm{\Delta}t\phantom{\rule{0.25em}{0ex}}{\stackrel{\mathrm{\u0303}}{\mathit{u}}}_{h}\left({t}^{n},{\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n}\right),\\ {\displaystyle}{\mathit{k}}_{\mathrm{2}}& {\displaystyle}={\displaystyle \frac{\mathrm{\Delta}t}{\mathrm{2}}}\phantom{\rule{0.25em}{0ex}}{\stackrel{\mathrm{\u0303}}{\mathit{u}}}_{h}\left({t}^{n}+{\displaystyle \frac{\mathrm{\Delta}t}{\mathrm{2}}},{\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n}+{\displaystyle \frac{{\mathit{k}}_{\mathrm{1}}}{\mathrm{2}}}\right),\\ {\displaystyle}{\mathit{k}}_{\mathrm{3}}& {\displaystyle}={\displaystyle \frac{\mathrm{\Delta}t}{\mathrm{2}}}\phantom{\rule{0.25em}{0ex}}{\stackrel{\mathrm{\u0303}}{\mathit{u}}}_{h}\left({t}^{n}+{\displaystyle \frac{\mathrm{\Delta}t}{\mathrm{2}}},{\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n}+{\displaystyle \frac{{\mathit{k}}_{\mathrm{2}}}{\mathrm{2}}}\right),\\ {\displaystyle}{\mathit{k}}_{\mathrm{4}}& {\displaystyle}=\mathrm{\Delta}t\phantom{\rule{0.25em}{0ex}}{\stackrel{\mathrm{\u0303}}{\mathit{u}}}_{h}\left({t}^{n+\mathrm{1}},{\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n}+{\mathit{k}}_{\mathrm{3}}\right),\\ {\displaystyle}{\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n+\mathrm{1}}& {\displaystyle}={\stackrel{\mathrm{\u0303}}{\mathit{x}}}^{n}+{\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}{\mathit{k}}_{\mathrm{1}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\mathit{k}}_{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\mathit{k}}_{\mathrm{3}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{6}}}{\mathit{k}}_{\mathrm{4}}.\end{array}$$RK4 requires the evaluation of the velocity at four points in space and three points in time.
The primary expense in all of the methods above is the evaluation of the velocity field ${\mathit{u}}_{h}^{n}$ and ${\mathit{u}}_{h}^{n+\mathrm{1}}$ at arbitrary positions x. Given that the velocity fields u_{h} that we consider here are often finiteelement fields defined with shape functions whose values are determined by mapping a reference cell $\widehat{K}$ to each cell K using a transformation $\mathit{x}={\mathrm{\Phi}}_{K}\left(\widehat{\mathit{x}}\right)$, the evaluation at arbitrary points requires the inversion of Φ_{K}, which is an expensive operation for nonlinear mappings such as those used in deformed or curved geometries.
3.4 Particle integrators used in the benchmarks
Based on our earlier work measuring the convergence properties of the integrators described above in analytically known flow (see the Supplement in Gassmöller et al., 2018), we expect FE and RK2FOT to converge with first order (in Δt) in time variable flow, while RK2 and RK4 are expected to converge with second order in time. RK2FOT is limited from reaching the potential of RK2 by the use of only a single velocity solution in time, and RK4 in our specific implementation (though not in general) is limited by only storing two velocity solutions, which only allows for a linear extrapolation from t^{n} to ${t}^{n}+\frac{\mathrm{\Delta}t}{\mathrm{2}}$ and t^{n+1}. Therefore, while there are valid reasons to choose either FE or RK4, we will limit our benchmark results to RK2 and RK2FOT, because we expect them to illustrate the significant difference between algorithms that are first or secondorder accurate in time.
For our benchmarks we want to reduce the coupling between Eqs. (1) and (2) with (3) to a minimum in order to precisely measure the influence of exactly one coupled property. This step also simplifies the construction of the benchmarks. Therefore, we focus on problems with constant viscosity (η=1) and no source terms (H_{c}=0). The advected material property ϕ_{c} we consider here is the density ρ. The final set of equations for our benchmarks will therefore be
The set of benchmarks we will consider is an extension of previously published benchmarks (Gassmöller et al., 2019). In order to explain the extensions to transient flow, here we will briefly revisit our approach to derive steady flow fields. In our earlier work we have considered incompressible flow fields u that were derived based on a known and timeindependent stream function Ψ. Under the assumption that viscosity is known and constant, and that boundary conditions are chosen to match the desired solution, we were able to compute the righthandside terms to Eq. (8) that satisfied the set of equations and therefore created an analytical benchmark for the whole system of equations. However, this only guarantees a consistent solution for the distribution of the density ρ at the current instant in time. It is therefore only an instantaneous benchmark solution. In order to create a steadystate flow field – defined as a velocity field u that does not change over time – the righthandside driving force needs to stay constant over time. In other words, the advected property ρ needs to be chosen in such a way that when it is advected with the flow field u, the righthandside ρg does not change over time. In order to find such a density distribution, we can make use of the definition of the streamline, which comprises lines of constant Ψ. If Ψ is independent of time, any property advected with the flow will be advected along the streamlines. Thus, if ρ is constant along the streamlines, the righthandside ρg will not change even if ρ is advected with the flow. Choosing ρ=Ψ is therefore an easy approach to guarantee a steadystate flow field.
The benchmarks below extend these steadystate models with a nonlinear time dependence, which will test how much error the chosen advection scheme accumulates over time when the velocity changes. In order to derive such solutions, we make use of the fact that we can superimpose two independent flow fields. In addition to a steady flow based on a stream function Ψ, we add a timedependent velocity that has two special properties. First, we ensure that this second flow field is purely forced by the boundary conditions instead of internal density forces. This choice opens up control over the exact velocity over time. It also implies that we do not have to modify the density distribution to add this second flow field, i.e., the density is still a function of the steadystate stream function Ψ. Second, we choose the timevariable flow field to be in the nullspace of the Stokes operator, e.g., solid body rotational flow in a spherical shell and translational flow in a box geometry. This ensures that the resulting modification only affects the velocity solution (but not the pressure) and can be interpreted as a (timedependent) coordinate transformation of a steady flow. We will consider one case in a 2D spherical shell and one case in a 2D box geometry, and we will discuss the specific flow fields in the following subsections.
4.1 A benchmark for a 2D spherical shell
As described, we start from an instantaneous solution for Stokes flow in a spherical shell and add a timedependent rotational flow that is enforced using the boundary conditions. A detailed derivation of the benchmark solution is given in Appendix A. It is important to note that while the benchmark is derived in polar coordinates, it is implemented in Cartesian coordinates in our code. Our final benchmark solution consists of a number of convection cells that rotate around the center of a spherical shell and is described by
The constants A and B along with the functions g(r), f(r), h(r), and m(r) are listed in Appendix A. The parameter k controls the number of upwellings and downwellings in the model and is chosen as k=4 for this study. The parameter ω(t) represents the timedependent solid body rotation and is chosen as ω(t)=e^{t}. τ(t) is a phase shift caused by the solid body rotation and is computed as $\mathit{\tau}\left(t\right)={\int}_{\mathrm{0}}^{t}\mathit{\omega}\left(s\right)\phantom{\rule{0.125em}{0ex}}\mathrm{d}s$. The spherical shell has an inner radius of R_{1}=1 and an outer radius of R_{2}=2. The setup of the benchmark and a snapshot of the solution is shown in Fig. 1.
We note that this solution can be interpreted as consistent with a stream function that is variable in time, with a flow field that conveniently advects the density in such a way as to satisfy our Stokes solution at the current point in time. We also note that this solution effectively consists of two parts: a densitydriven internal convection in small convection cells and a forced and analytically known rotational flow of the whole model.
4.2 A benchmark for a 2D box geometry
Our solution for the box benchmark is analogous to the spherical shell case, but we can build directly on our earlier model setup of Gassmöller et al. (2018). We add a solid body translation, and with periodic boundary conditions at both side boundaries this allows us to define a known, transient solution to the incompressible Stokes equations, which is described by
As for the spherical case described above, we will use a nonlinear choice for ω(t), namely ω(t)=e^{t}, and the phase shift τ is computed as before. We choose a box with a width of w=2 and height of h=1. The solution of the benchmark is shown in Fig. 2. The translation of the solution as well as the periodic boundary conditions also represent the main difference between our benchmark solution and the one presented in Trim et al. (2023a), which uses a steadystate flow with a timedependent velocity amplitude.
4.3 How we use these benchmarks in our particle advection algorithms
Adding time dependence to the benchmarks modifies the numerical solution and the accumulated error in distinct ways, depending on which advection method we choose. Here we will consider five cases:

We obtain a computed solution by using the exact density $\mathit{\rho}(x,y,t)$ defined in Eqs. (14) and (20). This solution will act as our baseline benchmark, illustrating the optimal convergence rate for the Stokes solver we used.

We use the (interpolated) exact density as an initial condition for the density advection Eq. (10), whose solution we then approximate using discontinuous, piecewise quadratic (DGQ_{2}) finite elements with a penalty method as described in He et al. (2016).

Same as case 2, but we use continuous, piecewise quadratic (Q_{2}) finite elements and an entropy viscosity stabilization technique (Guermond et al., 2011; Kronbichler et al., 2012). This is the default choice in ASPECT. In both cases 2 and 3 we use a backward differentiation formula (BDF2) timestepping algorithm that is secondorder accurate in time to solve the advection Eq. (10).

We use the exact density as the initial condition for particles whose position we advect using a secondorder accurate Runge–Kutta (RK2) algorithm. Where we need the density for the solution of the Stokes equations, we interpolate properties from particles onto a DGQ_{2} discontinuous finiteelement field and evaluate that at quadrature points as necessary.

Same as case 4, but we use RK2FOT as described in Sect. 3.3.
In order to limit ourselves to examining the accuracy in time of these five benchmark series, we will only consider a single combination of the Stokes finiteelement and particle interpolation algorithm in this paper. We will use a Q_{2}×Q_{1} Stokes element (continuous, piecewise quadratic velocity, and continuous, piecewise linear pressure) and a linear leastsquares particle interpolation algorithm with initially 64 particles per cell. We have described the influence of these choices in earlier work (Gassmöller et al., 2019).
Because the number of particles in a cell can change during the model run, we enforce a minimum of 12 particles per cell, which guarantees that the linear leastsquares interpolation algorithm is always sufficiently constrained. We do not limit the maximum number of particles per cell in these models. In practice, the presented benchmarks never require the addition of particles, and therefore the number of particles stays constant (for the box) or decreases by less than 0.01 % (for the annulus, caused by integration error close to the boundary, and then leading to the loss of particles across the boundary). The two tested integration schemes do not show a significant difference in particle loss in the annulus geometry, even though we have observed such differences between Runge–Kutta algorithms and Forward Euler integrators in our earlier work (Gassmöller et al., 2018).
In the following subsections, let us use the benchmarks derived above for the numerical evaluation of particle schemes.
5.1 Spherical shell benchmark
Figure 3 presents the L_{2} error norm of velocity, pressure, and density for the spherical annulus benchmark at a fixed time as a function of mesh resolution (left column) and at a fixed resolution as a function of time (right column).
The left column illustrates that all advection methods but RK2FOT reach secondorder convergence for the density with increasing resolution (bottom left panel). As expected, RK2FOT is limited by the available time information and only reaches firstorder convergence. An additional detail is that the field methods (Q_{2} and DGQ_{2}) have a larger error constant than the particle method (RK2), even for the same convergence rate. We will revisit the source of this error constant when discussing the error accumulation over time. Starting at moderate resolutions (around $h=\frac{\mathrm{1}}{\mathrm{16}}$) the RK2FOT model only reaches a firstorder convergence rate in velocity, while Q_{2} and DGQ_{2} reach second order. This result is important, because it illustrates that particles do not uniformly generate smaller errors than field methods but can indeed generate larger errors if their advection method is too simple and therefore inaccurate. The RK2 method maintains a thirdorder convergence rate in this metric up to very fine resolutions, which is surprising as the expected convergence order for RK2 is second order. We refer to Sect. 5.2 for a discussion on how RK2 may reach superconvergence for the resolutions shown in this benchmark.
The analysis of error evolution over time (the right column of Fig. 3) illustrates further differences between field and particle methods. Velocity and pressure errors reveal that RK2FOT accumulates the largest errors over time, as expected, followed by Q_{2} and DGQ_{2}. RK2 accumulates the smallest errors. However, the density error norm shows distinct differences between the methods. Both RK2 and RK2FOT start with the same error value, but while the RK2 error remains near constant over the evolution of the model (increasing by less than 2 %), the error of RK2FOT increases by almost two orders of magnitude. The rate of increase in the RK2FOT scheme changes towards the end of the model run. We show in Appendix B that this slowdown is related to the periodicity of our benchmark solution. The field methods Q_{2} and DGQ_{2} behave distinctly different. DGQ_{2} starts at a much smaller error value than all other methods, but it accumulates significant errors towards the end of the model run. Q_{2} already starts at a significantly larger error value than all other methods. This is likely related to the fact that the used entropy viscosity method falls back to a firstorder stabilization scheme for the first time step, which introduces a large amount of numerical diffusion at the model start (and only then). The overall shape of these curves is due to properties of the exact solution, not the method used, but is not of interest to us here.
Summarizing these findings, loworder particle methods show larger errors than the tested field methods, while higherorder particle methods outperform the field methods in our benchmark both with increasing resolution and with increasing model time. Therefore, whenever the other error sources of the solution are sufficiently small (i.e., if the Stokes element and timestepping scheme allow for higherorder accuracy), a higherorder particle scheme can significantly improve accuracy.
5.2 Box benchmark
The box benchmark results follow a similar pattern for the dependence of errors on the methods used (see Fig. 4). First, the solution using the analytical density solution produces a thirdorder convergence in velocity and secondorder convergence in pressure, which proves that the Stokes elements reach their optimal convergence order when given accurate density distributions. Second, and confirming theoretical predictions, the RK2 firstorder time (RK2FOT) advection method creates a firstorder accurate approximation for the density, which also generates a firstorder accurate pressure and velocity solution, thereby significantly limiting the potential accuracy of the Stokes elements. All other advection methods reach secondorder convergence as predicted by their derivations, however with significant differences in the absolute error norm. For all solution variables, particles advected using a full RK2 scheme reach about a one order of magnitude lower error norm at the end time than the Q_{2} and DGQ_{2} finiteelement methods, a value that depends on the chosen end time (compare right column of Fig. 4 and next paragraph). One feature to note is that the velocity error of the RK2 particle advection method starts with a thirdorder convergence rate at low resolutions and transitions to a secondorder convergence rate for higher resolutions (red line with red triangles in Fig. 4a). We assume that this transition is caused by a shift in the dominant error source: at large h the spatial error in the solution dominates, which is consistent with the observation that the particle solution (red line with red triangles in Fig. 4a) is very close to the solution computed with an analytical density (blue line with upsidedown triangles in Fig. 4a) and converges with third order. At smaller h the spatial error is reduced significantly, leaving the time error, which converges at second order as the remaining dominant source of velocity error. We observed such a transition in the dominant error source already in our earlier work (Gassmöller et al., 2019).
When evaluating the error norms of the solution as a function of time for a fixed resolution (right column of Fig. 4), we can gain additional insight into the properties of the advection methods. While it is obvious that the RK2FOT method remains the most inaccurate method at a sufficiently large time, this comparison also makes clear that it shares the same error value as RK2 at the start of the model, a value which is lower than the Q_{2} and DGQ_{2} methods. This is because the error in the first time step is dominated by the accuracy of the spatial approximation of the density. This also means that in benchmarks that are instantaneous or very short, the RK2FOT method will perform at nearly the same accuracy as the RK2 method, leading to misleading conclusions about its accuracy and suitability for timedependent geodynamic models. Both the Q_{2} and DGQ_{2} methods start at significantly larger errors in velocity and pressure, but they accumulate less error over time than RK2FOT, although more than RK2.
For our conclusion it is important to note that even though both particle methods start with a significantly smaller error than finiteelement advection methods, the firstorder accuracy of the RK2FOT scheme produces significantly larger errors, and that this effect becomes more pronounced with increasing resolution and increasing model runtime.
5.3 Accuracy and performance discussion
Summarizing the benchmark results, firstorder particle methods yield larger errors than the tested fieldbased methods, while higherorder particle methods outperform the investigated fieldbased methods both with increasing resolution and increasing model time. Therefore, whenever other error sources of the solution are sufficiently small (i.e., if the Stokes element and timestepping scheme allow for higherorder accuracy), a higherorder particle scheme can significantly enhance the accuracy of the solution. Even though we cannot prove it here, this conclusion is likely also true for the common case of a solution that is not smooth enough to allow for the optimal convergence rate of RK2. For discontinuous solutions the convergence rate of higherorder finite elements can break down to the same rate as for firstorder elements (e.g., Thielmann et al., 2014). However, solutions are rarely fully discontinuous and instead contain a mix of smooth and nonsmooth regions. Additionally, despite showing the same convergence rate, higherorder elements have still delivered higher accuracy in absolute terms than lowerorder methods in many benchmarks (Kronbichler et al., 2012; Thieulot and Bangerth, 2022). We speculate that the same results would be seen for higherorder advection methods in time, although the construction of appropriate benchmarks would be challenging.
Finally, the improved accuracy of higherorder particle methods has to be discussed in the context of their larger memory requirement and computational cost. RK2 requires the storage of two velocity solutions instead of a single solution like RK2FOT; thus, very coarsely (neglecting the memory cost of the particles) one could consider RK2 to be twice as expensive in terms of memory. However, this additional cost has to be compared with the total memory requirement of a modern geodynamic model and is only relevant if models are typically limited by the available memory. Many modern Stokes solvers in geodynamics either rely on matrixbased algorithms or Krylov subspace solvers with a long recurrence relation (e.g., GMRES), both of which can easily require the memory of tens to hundreds of solution vectors. For all of these models, storing an additional velocity solution for the particle advection represents a negligible memory cost. Even if models are computed with modern matrixfree solvers with shortrecurrence relations (e.g., Clevenger and Heister, 2021), in our experience the size of models is typically limited by the available computational power or memory bandwidth, but rarely by the available total memory. We therefore assume that storing an additional velocity solution is not a prohibitive cost and focus for the rest of this section on investigating the performance of RK2 over RK2FOT.
In theory, and at the most granular level, RK2 could be expected to require 50 % more memory bandwidth than RK2FOT, because simplistically speaking its second stage computes ${\mathit{x}}^{n+\mathrm{1}}={\mathit{x}}^{n}+\frac{\mathrm{\Delta}t}{\mathrm{2}}\phantom{\rule{0.25em}{0ex}}\left({\mathit{u}}^{n}+{\mathit{u}}^{n+\mathrm{1}}\right)$ (reading three vector entries from memory) instead of ${\mathit{x}}^{n+\mathrm{1}}={\mathit{x}}^{n}+\mathrm{\Delta}t\phantom{\rule{0.25em}{0ex}}{\mathit{u}}^{n}$ (reading two). If the algorithm is bandwidth limited, it could therefore incur a 33 % performance penalty. However, this is of course only a small part of the total particle advection cost, which also includes algorithms that are independent of the integration scheme (like the first Runge–Kutta stage, or determining the interpolation functions from grid to particle location). Finally, it is important to consider what fraction of the total model cost is used to advect the particles, as even if RK2 is significantly more expensive than RK2FOT, it may deliver higher accuracy for a small fraction of the total model cost.
All performance results were computed on one core of an AMD EPYC 7453 processor with the software listed in the Data Availability statement.
Table 1 illustrates these metrics for a version of the box benchmark that was shortened to a model runtime of t=0.1. As can be seen, particle advection only requires an almost insignificant percentage of the total model runtime (always ≤3.7 %). Additionally, the RK2 advection algorithm only incurs an additional cost compared with RK2FOT on the order of 0.4 % of the total model runtime or a relative increase of 10 %–21 % compared with the cost of RK2FOT. Moreover, it is clearly visible that this additional cost is approximately constant across resolutions, which means that for a constant and relatively small additional cost these models deliver significantly more accurate solutions as shown in the previous sections. In our opinion, these results illustrate that as long as the memory cost of storing an additional velocity solution is acceptable, using RK2 with an extrapolated or interpolated half time step is generally superior to a firstorder method like RK2FOT in applications that depend on accurate solutions.
Above we have illustrated the influence of algorithmic choices on the accuracy and performance of benchmark results. However, this does not by itself justify the increased cost of such an algorithm in practical models: perhaps, in typical geodynamic applications, the error due to a loworder time approximation is negligible compared with other error sources, and therefore a simple advection method may be sufficient. In the following, we use an application model to show that the higher accuracy is indeed important and can influence firstorder outcomes and the interpretation of a geodynamic study.
In order to illustrate this point, we use an example where the property carried on the particles (the grain size d) nonlinearly influences the material properties (the viscosity η) and the corresponding solution of the equations. Our model setup (Fig. 5a) consists of the oceanic lithosphere and asthenosphere, down to a depth of 410 km, moving away from a spreading center at the top left corner of the model and horizontally extending to a distance of 4000 km from the ridge. The grid does not change over time and is spatially uniform with a cell size of 12.8 km (for a total of 320 cells in the x and 32 cells in y direction). The plate speed is prescribed in the horizontal direction to 5 cm yr^{−1} at the top, and right (outflow) boundaries at depths smaller than 100 km. The outflow velocity then linearly decreases with depth starting at 100 km depth towards 0 cm yr^{−1} at the bottom of the model. The left (ridge axis) boundary of the model is closed and stressfree (free slip). The vertical velocity component is not prescribed at the bottom and right boundaries; instead, a depthdependent hydrostatic pressure profile, which is computed at the model start and is constant in time, is prescribed. Therefore, material can flow in beneath the ridge axis and leaves the model either through the bottom or the right boundary. The initial temperature follows an adiabatic profile with a potential temperature of 1623 K and a halfspace cooling profile close to the top boundary with an age consistent with the plate velocity. The initial temperature also includes a small (r=10 km) thermal perturbation at the ridge axis to support the onset of spreading. The initial grain size is set to d=5 mm and also includes a small (r=30 km) Gaussian anomaly, which reduces the grain size to d=2.1 mm close to the ridge axis. Since the temperature is prescribed at the top boundary, the plate cools conductively over time until smallscale convection sets in at the base of the plate.
In this model, we use particles to carry information about the mineral grain size d, which influences the viscosity nonlinearly as
resulting in an effective viscosity of
Here, η_{diff}, η_{dis}, and η_{eff} are the diffusion, dislocation, and effective viscosity, respectively. ${A}_{\mathrm{diff}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{15}}$ m^{3} Pa^{−1} s^{−1} and ${A}_{\mathrm{dis}}={\mathrm{10}}^{\mathrm{15}}$ Pa^{−3.5} s^{−1} are diffusion and dislocation creep prefactors, E_{diff}=375 kJ mol^{−1} and E_{dis}=530 kJ mol^{−1} the activation energy for diffusion and dislocation creep, ${V}_{\mathrm{diff}}=\mathrm{4}\times {\mathrm{10}}^{\mathrm{6}}$ m^{3} mol^{−1} and ${V}_{\mathrm{dis}}=\mathrm{1.4}\times {\mathrm{10}}^{\mathrm{5}}$ m^{3} mol^{−1} the respective activation volumes, R is the universal gas constant, P the pressure, T the temperature, m=3 the grain size exponent of diffusion creep, ε_{dis,II} the square root of the second moment invariant of the dislocation strain rate, and n=3.5 the dislocation creep strain rate exponent. We limit the viscosity computed in Eq. (25) to be between 10^{16} and 10^{23} Pa s.
In addition, particles are not just advected, but both the temperature and the strain rate in the model influence the evolution of the grain size. For a single particle moving along the flow field, we describe this evolution via the equation
which implies that the grain size field d(x,t) satisfies the equation
Here, ${k}_{\mathrm{g}}=\mathrm{1.92}\times {\mathrm{10}}^{\mathrm{10}}$ m^{3} s^{−1} is the grain size growth prefactor, and E_{g}=400 kJ mol^{−1} and V_{g}=0 m^{3} mol^{−1} the activation energy and volume for grain size growth, respectively. p_{g}=3 is the grain size growth exponent, c=3 a geometric constant, λ=0.1 the fraction of deformation work that goes into changing the grain boundary area, and γ=1 the average specific grain boundary energy. The terms on the righthand side of these equations describe how the dynamic grain size increases over time (with a nonlinear dependence on temperature and grain size itself) and how it is decreased by dynamic recrystallization due to strain accommodated by dislocation creep (which depends nonlinearly on stress and temperature). (For a detailed discussion of these terms and all parameter values, we refer the reader to Dannberg et al., 2017.) We solve this model using the particlebased RK2 and RK2FOT advection schemes for different time step lengths expressed as a fraction of the CFL number. We make sure that except for the advection scheme, all other parts of the model and algorithms are identical. In particular, we generate particles in identical and deterministic random locations and we enforce the same minimum (40) and maximum (250) number of particles per cell, and we make sure that all algorithms for particle addition and deletion are deterministic.
As can be seen in Fig. 5b and c, the two advection methods produce noticeably different locations of onset of convection. While the model with a full RK2 advection scheme develops smallscale convection beneath the oceanic plate at a distance of approximately 1900 km from the ridge (from 1893.2 to 1969.2 km), the model with RK2FOT develops the onset at varying distances from 1656.0 to 1916.0 km. These numbers correspond to plate ages of 37.9–39.4 Myr (RK2) and 34.3–38.3 Myr (RK2FOT), respectively. Because of the strong nonlinearity of these models, we do not observe a simple convergence to one solution as in the benchmark results for either of the models. However, it is apparent that the RK2 method produces a much more stable onset location of smallscale convection and a greater similarity of the other downwellings that develop beyond the initial onset. In contrast, the onset of convection varies significantly in the RK2FOT method depending on the time step size. In addition, the downwellings show very different convection structures. One could speculate that the RK2FOT method starts to converge towards the RK2 results for a CFL number of 0.075; however, this is not certain given the strong variations in the RK2FOT results. Considering our benchmark results in the previous sections, we would expect the RK2FOT method to converge to the same solution as RK2, though requiring substantially smaller time steps. At least in this application it would clearly require time steps that make the model prohibitively expensive.
More importantly, one could assume that the shown variations just illustrate temporal variability in the convection pattern, and that the RK2FOT results are only a temporary state at the end of the model (200 Myr). To investigate this question we track the onset of convection over time for both RK2 and RK2FOT, and show the results in Table 2. We find that the different outcomes between the advection methods is not just a temporary state at this precise time, but the models show systematic differences over long ranges of time. The onset of convection is relatively stable in the model with RK2, while it varies significantly in the model with RK2FOT and occurs systematically closer to the ridge as time passes. Even though this may be a temporary development that will eventually reverse, we have not observed similar behavior for RK2.
The exact timing of the onset of convection beneath an oceanic plate is relevant for the argument that smallscale convection causes a flattening of topography in seafloor subsidence datasets, and therefore ultimately relevant also for supporting the plate model of oceanic lithosphere cooling (Stein and Stein, 1992; Huang and Zhong, 2005; Richards et al., 2018). In addition to the difference in the onset of convection, the characteristic length scale at which instabilities develop below the lithosphere is significantly smaller for the RK2FOT method, visible in the larger number and smaller distance between convection cells in Fig. 5b and c. This is especially relevant as the distance between seismic anomalies associated with smallscale convection is a constraint from seismic studies and can be used to validate geodynamic models (Eilon et al., 2022). We want to emphasize here that our model is conceptual and not intended to produce realistic timings or length scales, and that a misprediction of these quantities due to inaccurate particle algorithms in models has concrete consequences for the interpretation of geodynamic features observed on Earth.
Because both the onset of smallscale convection and the length scale of convection cells is governed by the growth of small instabilities in a boundary layer, it is reasonable to assume that the lower accuracy of RK2FOT supports this growth of instabilities and explains the earlier onset of convection. The growth of instabilities in a boundary layer (or internally layered systems) is one of the most common processes for developing flow features in convecting systems like the Earth's mantle and lithosphere. Examples are the generation of plumes at the core–mantle boundary, the stagnation of subducted slabs or plumes at phase transitions, or the initiation of plate boundaries in models of lithosphere dynamics. We therefore infer from our results that models of all of these processes can benefit from incorporating more accurate particle advection methods, and that predictions of models using lowerorder advection schemes may need to be adjusted or reproduced in higherresolution studies.
We have shown in our benchmarks and applications that implementing accurate particle algorithms, in particular higher order in time, can significantly improve the numerical accuracy of geodynamic models. One of the conclusions of our benchmarks is that commonly used particle advection methods that are higher order in space but first order in time acquire significant amounts of numerical error in timevariable flow, which becomes more pronounced the higher the resolution and the longer the model runs. The reason this error is not often discussed in the geodynamic literature is that traditional benchmarks that either rely on instantaneous analytical solutions or on steadystate solutions cannot show this error by their design. Only model comparison studies or benchmarks with analytical solutions in transient flow can point out this error source. Given that many geodynamic finiteelement models already use Stokes elements that allow for higherorder accuracy to ensure stability (e.g., Taylor–Hood Q_{2}×Q_{1} or ${Q}_{\mathrm{2}}\times {P}_{\mathrm{1}}$), it would be straightforward to extend their particle advection algorithms to a secondorderintime method. While this can increase the cost for evaluating velocities at the particle locations, our results show that the increased convergence order and improved accuracy of the model results is well worth the additional cost. Of course, in order to increase the overall model accuracy, all other employed algorithms need to support the same accuracy.
We believe that sharper focus on quantifying the numerical accuracy of geodynamic models will generate more trust in geodynamic model solutions and increase the impact of the discipline of geodynamic modeling as a whole. We provide the reference implementation of our algorithms and benchmarks in the opensource community software ASPECT and hope that they are useful to the community at large.
Extending our previously published spherical benchmark (Gassmöller et al., 2019) seems to be straightforward by adding a timedependent solid body rotation to the existing solution. However, because our earlier solution is already a purely rotational flow, an additional timedependent rotation does not create a transient solution for the density, not allowing to intuitively measure the accuracy of the combined particle–finiteelement algorithm. In other words, an error in the particle position along the streamline because of the time variability in the flow would not change the density distribution and therefore would not translate into an error in the Stokes solution.
We therefore begin by deriving a new exact solution to the stationary, incompressible Stokes equations for an isoviscous, isothermal fluid in a 2D annulus. Given the geometry of the problem, we work in polar coordinates. We denote the orthonormal basis vectors by e_{r} and e_{θ}, the inner radius of the annulus by R_{1}, and the outer radius by R_{2}. Furthermore, we assume that the viscosity is a constant η=1 and set the gravity vector to an inwardpointing vector $\mathit{g}={g}_{r}{\mathit{e}}_{r}$, with g_{r}=1.
Given these assumptions, the incompressible Stokes equations in the annulus are (Schubert et al., 2001)
Equations (A1) and (A2) are the momentum equations in polar coordinates, while Eq. (A3) is the incompressibility constraint.
We then seek solutions whose circumferential velocity can be written as
where k is an integer and where f(r) will be specified later. From Eq. (A3) we then obtain
leading to
where
Since we want to fix the velocity to be tangential at both boundaries we have
for all $\mathit{\theta}\in [\mathrm{0},\mathrm{2}\mathit{\pi}]$. We choose
in analogy to the solution of the Laplace equation in Chap. 6 of Strauss (2007), and thus
where C is a nonzero constant of integration. Given the boundary conditions in Eq. (A8) we find that
In this work we choose $C=\mathrm{1}$. Our earlier choice of f means that
so that Eq. (A2) simplifies to
which together with Eq. (A4) leads to
where l(r) comes from integration with respect to θ and $h\left(r\right)=\left(\mathrm{2}g\right(r)f(r\left)\right)/r$. We now insert Eq. (A15) into Eq. (A1) to obtain
with
since it is easy to verify using Eq. (A7) that ${g}^{\prime}\left(r\right)=(fg)/r$.
Taking k=0 yields $\mathit{\rho}(r,\mathit{\theta})={l}^{\prime}\left(r\right)$, so we choose ${l}^{\prime}\left(r\right)={\mathit{\rho}}_{\mathrm{0}}$. In this case,
which represents a familiar linear pressure increase with depth for constant density and gravity, and where we have imposed $p(r,\mathit{\theta})=\mathrm{0}$ at the outer radius r=R_{2}.
In summary, Eqs. (A4), (A6), and (A15) form a solution of the incompressible Stokes equations, which fully stated reads
with
We can use the velocity solution for v_{r} and v_{θ} to determine a stream function for this flow field, which will be used to derive the stationary benchmark below:
The solution above is time independent and only valid for instantaneous models where the density is not advected. To make it time dependent, we first modify the density and gravity to create a steadystate variant of the benchmark and add a known timedependent component to the velocity as described in Sect. 4. We start by choosing a density field consistent with the streamline $\mathit{\rho}(r,\mathit{\theta})=\mathrm{\Psi}(r,\mathit{\theta})$. In our concrete benchmark this solution no longer satisfies the derived Stokes solution. However, we can recover an analytic solution by exploiting the fact that for the incompressible Stokes equations the density ρ only enters the computation as a product with the gravity g. Therefore, if as described in the example case above m(r)ksin (kθ) is the righthandside force term that satisfies the Stokes equation, we can still choose the density arbitrarily (e.g., $\mathit{\rho}(r,\mathit{\theta})=\mathrm{\Psi}(r,\mathit{\theta})$), as long as we define the gravity to be $\mathit{g}(r,\mathit{\theta})=m\left(r\right)k\mathrm{sin}\left(k\mathit{\theta}\right)/\mathit{\rho}(r,\mathit{\theta})$. This keeps the original forcing term constant, and thus makes the solution independent of time. The steadystate solution therefore is the same as above, except
with all other constants chosen as before.
In order to transform this steadystate benchmark into a known transient solution, we then add a solid body rotation with a nonlinear timedependent rotational velocity to the flow field. Since solid body rotations lie in the nullspace of the incompressible Stokes equations on an annular domain, the resulting flow field will still be a solution of the incompressible Stokes equations. This approach will work as long as we perform an appropriate rotation of all components of the solution, and it is equivalent to defining the solution in a rotating reference frame. We therefore modify the velocity components in θ direction to
Here τ(t) is a phase shift and ω(t) is an angular velocity. The phase shift τ(t) can be computed as the time integral of the angular velocity from the beginning of the model up to the present time t:
In order to not make the problem too simple, we forgo the case of a constant angular velocity and instead choose ω(t)=e^{t}, resulting in $\mathit{\tau}\left(t\right)={e}^{t}\mathrm{1}$.
Since the modification of the velocity in Eq. (A35) by the solid body rotation r ω(t) lies in the nullspace of the Stokes equations, it is straightforward to compute the modifications of the remaining solution variables, which only involves adding the phase shift to the θ coordinate.
The final consideration is how to achieve this prescribed rotation in the model. Since in the incompressible Stokes equations stresses are transmitted instantaneously throughout the entire domain, we can use the exact, known velocities as boundary conditions and expect the motion to apply equally to the entire model domain.
In order to better understand the accuracy of the RK2 method and investigate the source of the error decrease late in the model, we show in Fig. B1 a detailed comparison of only RK2 against the analytical density method. In order to visualize the difference between RK2 and analytical density over resolution (Fig. B1a and c), we no longer plot the absolute error in the L_{2} norm but instead the relative difference in error between RK2 and analytical density; i.e., if ${\mathit{\u03f5}}_{\mathrm{RK}\mathrm{2}}=\Vert \mathit{u}{\mathit{u}}_{h}^{\mathrm{RK}\mathrm{2}}{\Vert}_{{L}_{\mathrm{2}}}$ and ${\mathit{\u03f5}}_{\mathrm{AD}}=\Vert \mathit{u}{\mathit{u}}_{h}^{\mathrm{AD}}{\Vert}_{{L}_{\mathrm{2}}}$, then we plot $({\mathit{\u03f5}}_{\mathrm{RK}\mathrm{2}}{\mathit{\u03f5}}_{\mathrm{AD}})/{\mathit{\u03f5}}_{\mathrm{AD}}$. This way of plotting the error illustrates if both error values converge at the same rate – leading to a constant relative difference between the two errors – or if the RK2 error indeed converges at a lower rate – leading to a linearly (or higherorder) increasing relative error difference towards smaller h. As it turns out, the relative error indeed increases linearly with resolution, meaning that RK2 converges at only second order; however, the secondorder contribution is so small that it is not yet visible in the corresponding plot of Fig. 3. On the other hand, the pressure (Fig. 3c) converges at the same secondorder rate for both analytical density and RK2, leading to a constant relative difference between the two errors.
Turning to the evolution of error over model time (Fig. B1b, d, and e) reveals that what looked like a constant error value in Fig. 3 indeed follows the same trends as the other methods, only at drastically reduced error values. While the error values for velocity and density increased by 1–2 orders of magnitude from model start to end for the other advection methods (RK2FOT, Q_{2}, and DGQ_{2}), they increased by at most ≈2 % for RK2. Additionally, RK2 features the same error reduction close to t=1 as the other methods. Finally, it becomes apparent that even the model using analytical densities features a small but growing velocity error. Because the density error in this model is zero (analytical density), it seems reasonable to assume this error is a result of the Stokes solver. The accuracy of the Stokes solver depends on the absolute value of velocity, which increases exponentially over time. In other words, the blue line in Fig. B1b represents the best possible accuracy any advection method could reach for the given Stokes solver if it transported material information with perfect accuracy. Considering all the results presented in this section, we consider the RK2 scheme to be very close to achieving this theoretical limit.
To understand the reduction in velocity error and density error at certain model times requires us to take a closer look at the benchmark solution. Particularly relevant is that the benchmark solution is rotation symmetric, with four regions of upwelling and four regions of downwelling. Therefore, rotating the density field by 90° at any given time would lead to exactly the same solution. For some reason the reduction in velocity error in RK2 coincides almost exactly with a quarter rotation of the model solution at $t=\mathrm{ln}(\mathrm{1}+\mathit{\pi}/\mathrm{2})\approx \mathrm{0.944}$, the reduction for Q_{2} coincides with a rotation by three quarters of a rotation ($t=\mathrm{ln}(\mathrm{1}+\mathrm{3}\mathit{\pi}/\mathrm{2})\approx \mathrm{1.742}$), and the reduction for RK2FOT coincides with a full rotation ($t=\mathrm{ln}(\mathrm{1}+\mathrm{2}\mathit{\pi})\approx \mathrm{1.985}$). We speculate that the times of rotation symmetry with the starting solution allows for a resonance between the accumulated error in the numerical solution and the analytical solution at a slightly different time. This interaction would allow for an apparent reduction in error that does not actually exist, which is consistent with the observation that all errors rapidly increase again after the minimum. However, while this theory explains why reductions in error could happen at specific times, we have no explanation as to why the anomaly happens at different multiples of the rotation symmetry for different advection methods. We can only speculate that the occurrence depends on a very specific feature in each model, for example, how close individual discrete time steps end at the analytically determined times of rotation symmetry. Independent of the origin of the anomaly, the results of the convergence studies show that it does not influence the measurement of the convergence orders of different methods.
Computations were done using the ASPECT code (Heister et al., 2017; Kronbichler et al., 2012; Bangerth et al., 2023; Bangerth et al., 2022). ASPECT is published under the GPL2 license, and the necessary software version and model configuration to reproduce the results is published on Zenodo at https://doi.org/10.5281/zenodo.10805269 (Gassmöller, 2024).
RG devised the study, devised and ran the benchmark cases, and wrote the majority of the paper. JD provided and described the grain size application model. WB developed the integrator error analysis. EGP provided the particle interpolation algorithm. EGP and CT developed the instantaneous solution of the annulus benchmark case. All authors jointly interpreted the results and improved the paper.
The contact author has declared that none of the authors has any competing interests.
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.
The authors acknowledge University of Florida Research Computing for providing computational resources and support that contributed to the research results reported in this paper. Computations were also run on the Stampede2 system at the Texas Advanced Computing Center (TACC) as part of award TGEAR080022N.
We thank the Computational Infrastructure for Geodynamics – funded by the National Science Foundation under awards EAR1550901 and EAR2149126 – for supporting the development of ASPECT.
Rene Gassmöller and Wolfgang Bangerth were partially supported by the National Science Foundation under award OCI1148116 as part of the Software Infrastructure for Sustained Innovation (SI2) program; and by the Computational Infrastructure in Geodynamics initiative (CIG), through the National Science Foundation under awards EAR0949446, EAR1550901, and EAR2149126. Wolfgang Bangerth was also supported by the National Science Foundation under awards OAC1835673 and EAR1925595.
Elbridge Gerry Puckett was supported by the National Science Foundation under award ACI1440811 as part of the SI2 Scientific Software Elements (SSE) program.
Juliane Dannberg and Rene Gassmöller were supported by the National Science Foundation under awards EAR1925677 and EAR2054605.
This research has been supported by the National Science Foundation (grant nos. EAR1550901, EAR2149126, OCI1148116, EAR0949446, OAC1835673, EAR1925595, ACI1440811, EAR1925677, and EAR2054605).
This paper was edited by Boris Kaus and reviewed by two anonymous referees.
Adamuszek, M., Dabrowski, M., and Schmid, D. W.: Folder: A numerical tool to simulate the development of structures in layered media, J. Struct. Geol., 84, 85–101, 2016. a
Arndt, D., Fehn, N., Kanschat, G., Kormann, K., Kronbichler, M., Munch, P., Wall, W. A., and Witte, J.: ExaDG: HighOrder Discontinuous Galerkin for the ExaScale, in: Software for Exascale Computing – SPPEXA 20162019, edited by: Bungartz, H.J., Reiz, S., Uekermann, B., Neumann, P., and Nagel, W. E., Springer International Publishing, Cham, 189–224, https://doi.org/10.1007/9783030479565_8, 2020. a
Arndt, D., Bangerth, W., Bergbauer, M., Feder, M., Fehling, M., Heinz, J., Heister, T., Heltai, L., Kronbichler, M., Maier, M., Munch, P., Pelteret, J.P., Turcksin, B., Wells, D., and Zampini, S.: The deal.II
Library, Version 9.5, J. Numer. Math., 31, 231–246, https://doi.org/10.1515/jnma20230089, 2023. a
Arnould, M., Coltice, N., Flament, N., and Mallard, C.: Plate tectonics and mantle controls on plume dynamics, Earth Planet. Sc. Lett., 547, 116439, https://doi.org/10.1016/j.epsl.2020.116439, 2020. a
Baes, M., Sobolev, S., Gerya, T., and Brune, S.: PlumeInduced Subduction Initiation: SingleSlab or MultiSlab Subduction?, Geochem. Geophy. Geosys., 21, e2019GC008663, https://doi.org/10.3389/feart.2021.766604, 2020. a
Bangerth, W., Dannberg, J., Fraters, M., Gassmoeller, R., Glerum, A., Heister, T., Myhill, R., and Naliboff, J.: ASPECT: Advanced Solver for Problems in Earth’s ConvecTion, Figshare [data set], https://doi.org/10.6084/m9.figshare.4865333, 2023. a
Bangerth, W., Dannberg, J., Fraters, M., Gassmoeller, R., Glerum, A., Heister, T., Myhill, R., and Naliboff, J.: ASPECT v2.5.0, Zenodo [code], https://doi.org/10.5281/zenodo.8200213, 2023. a
Bercovici, D. and Ricard, Y.: Plate tectonics, damage and inheritance, Nature, 508, 513–516, 2014. a
Billen, M. I.: Modeling the dynamics of subducting slabs, Annu. Rev. Earth Pl. Sci., 36, 325–356, 2008. a
Brandenburg, J., Hauri, E. H., van Keken, P. E., and Ballentine, C. J.: A multiplesystem study of the geochemical evolution of the mantle with forcebalanced plates and thermochemical effects, Earth Planet. Sc. Lett., 276, 1–13, 2008. a
Clevenger, T. C. and Heister, T.: Comparison between algebraic and matrixfree geometric multigrid for a Stokes problem on adaptive meshes with variable viscosity, Numerical Linear Algebr., 28, e2375, https://doi.org/10.1002/nla.2375, 2021. a
Dannberg, J. and Gassmöller, R.: Chemical trends in ocean islands explained by plume–slab interaction, P. Natl. Acad. Sci. USA, 115, 4351–4356, 2018. a
Dannberg, J., Eilon, Z., Faul, U., Gassmöller, R., Moulik, P., and Myhill, R.: The importance of grain size to mantle dynamics and seismological observations, Geochem. Geophy. Geosy., 18, 3034–3061, 2017. a
Duretz, T., May, D. A., Gerya, T. V., and Tackley, P. J.: Discretization errors and free surface stabilization in the finite difference and markerincell method for applied geodynamics: A numerical study, Geochem. Geophy. Geosy., 12, Q07004, https://doi.org/10.1029/2011GC003567, 2011. a
Eilon, Z. C., Zhang, L., Gaherty, J. B., Forsyth, D. W., and Russell, J. B.: SubLithospheric SmallScale Convection Tomographically Imaged Beneath the Pacific Plate, Geophys. Res. Lett., 49, e2022GL100351, https://doi.org/10.1029/2022GL100351, 2022. a
El Geitani, T., Golshan, S., and Blais, B.: Toward HighOrder CFDDEM: Development and Validation, Industrial & Engineering Chemistry Research, 62, 1141–1159, https://doi.org/10.1021/acs.iecr.2c03546, 2023. a
Faccenna, C., Oncken, O., Holt, A. F., and Becker, T. W.: Initiation of the Andean orogeny by lower mantle subduction, Earth Planet. Sci. Lett., 463, 189–201, 2017. a
Farnetani, C. G. and Richards, M. A.: Numerical investigations of the mantle plume initiation model for flood basalt events, J. Geophys. Res.Sol. Ea., 99, 13813–13833, 1994. a
Furuichi, M. and May, D. A.: Implicit solution of the material transport in Stokes flow simulation: Toward thermal convection simulation surrounded by free surface, Comput. Phys. Commun., 192, 1–11, 2015. a
Gassmöller, R.: Benchmarking the accuracy of higher order particle methods in geodynamic models of transient flow: Data, Zenodo [data set], https://doi.org/10.5281/zenodo.10805269, 2024. a
Gassmöller, R., Lokavarapu, H., Heien, E., Puckett, E. G., and Bangerth, W.: Flexible and scalable particleincell methods with adaptive mesh refinement for geodynamic computations, Geochem. Geophy. Geosy., 19, 3596–3604, 2018. a, b, c, d, e, f, g
Gassmöller, R., Lokavarapu, H., Bangerth, W., and Puckett, E. G.: Evaluating the accuracy of hybrid finite element/particleincell methods for modelling incompressible Stokes flow, Geophys. J. Int., 219, 1915–1938, 2019. a, b, c, d, e, f, g, h
Gerya, T. V. and Yuen, D. A.: Characteristicsbased markerincell method with conservative finitedifferences schemes for modeling geological flows with strongly variable transport properties, Phys. Earth Planet. Int., 140, 293–318, 2003. a, b, c, d, e
Glerum, A., Brune, S., Stamps, D. S., and Strecker, M. R.: Victoria continental microplate dynamics controlled by the lithospheric strength distribution of the East African Rift, Nat. Commun., 11, 1–15, 2020. a
Golshan, S. and Blais, B.: LoadBalancing Strategies in Discrete Element Method Simulations, Processes, 10, 79, https://doi.org/10.3390/pr10010079, 2022. a
Golshan, S., Munch, P., Gassmöller, R., Kronbichler, M., and Blais, B.: LetheDEM: An opensource parallel discrete element solver with load balancing, Computational Particle Mechanics, 1–20, https://doi.org/10.1007/s40571022004786, 2022. a
Grima, A. G., LithgowBertelloni, C., and Crameri, F.: Orphaning regimes: the missing link between flattened and penetrating slab morphologies, Front. Earth Sci., 8, 374, https://doi.org/10.3389/feart.2020.00374, 2020. a
Guermond, J.L., Pasquetti, R., and Popov, B.: Entropy viscosity method for nonlinear conservation laws, J. Comput. Phys., 230, 4248–4267, 2011. a
Gülcher, A. J. P., Ballmer, M. D., and Tackley, P. J.: Coupled dynamics and evolution of primordial and recycled heterogeneity in Earth's lower mantle, Solid Earth, 12, 2087–2107, https://doi.org/10.5194/se1220872021, 2021. a
Gurnis, M. and Hager, B. H.: Controls of the structure of subducted slabs, Nature, 335, 317–321, 1988. a
Hairer, E. and Wanner, G.: Solving Ordinary Differential Equations II. Stiff and DifferentialAlgebraic Problems, SpringerVerlag, Berlin, ISBN10: 3540604529, ISBN13: 9783540604525, 1991. a
He, Y., Puckett, E. G., and Billen, M. I.: A Discontinuous Galerkin Method with a Bound Preserving Limiter for the Advection of nonDiffusive Fields in Solid Earth Geodynamics, Phys. Earth Planet. Int., 263, 23–37, https://doi.org/10.1016/j.pepi.2016.12.001, 2016. a
Heister, T., Dannberg, J., Gassmöller, R., and Bangerth, W.: High accuracy mantle convection simulation through modern numerical methods. II: Realistic models and problems, Geophys. J. Int., 210, 833–851, 2017. a, b, c
Huang, J. and Zhong, S.: Sublithospheric smallscale convection and its implications for the residual topography at old ocean basins and the plate model, J. Geophys. Res.Sol. Ea., 110, B05404, https://doi.org/10.1029/2004JB003153, 2005. a
Jones, T. D., Sime, N., and van Keken, P.: Burying Earth's primitive mantle in the slab graveyard, Geochem. Geophy. Geosy., 22, e2020GC009396, https://doi.org/10.1029/2020GC009396, 2021. a
Kaus, B. J., Popov, A. A., Baumann, T., Pusok, A., Bauville, A., Fernandez, N., and Collignon, M.: Forward and inverse modelling of lithospheric deformation on geological timescales, in: Proceedings of nic symposium, John von Neumann Institute for Computing (NIC), NIC Series, 8, 978–983, https://hdl.handle.net/2128/10411 (last access: 9 May 2024) 2016. a
Kellogg, L. and Turcotte, D.: Mixing and the distribution of heterogeneities in a chaotically convecting mantle, J. Geophys. Res.Sol. Ea., 95, 421–432, 1990. a
Kramer, S. C., Davies, D. R., and Wilson, C. R.: Analytical solutions for mantle flow in cylindrical and spherical shells, Geosci. Model Dev., 14, 1899–1919, https://doi.org/10.5194/gmd1418992021, 2021. a
Kronbichler, M., Heister, T., and Bangerth, W.: High Accuracy Mantle Convection Simulation through Modern Numerical Methods, Geophys. J. Int., 191, 12–29, 2012. a, b, c, d, e
Lin, S.C. and van Keken, P. E.: Dynamics of thermochemical plumes: 1. Plume formation and entrainment of a dense layer, Geochem. Geophy. Geosy., 7, Q02006, https://doi.org/10.1029/2005GC001071, 2006. a
McNamara, A. K. and Zhong, S.: Thermochemical structures within a spherical mantle: Superplumes or piles?, J. Geophys. Res., 109, 1–14, https://doi.org/10.1029/2003JB002847, 2004. a, b
McNamara, A. K. and Zhong, S.: Thermochemical structures beneath Africa and the Pacific Ocean, Nature, 437, 1136–1139, 2005. a
Moresi, L., Zhong, S., Han, L., Conrad, C., Tan, E., Gurnis, M., Choi, E., Thoutireddy, P., Manea, V., McNamara, A., Becker, T., Leng, W., and Armendariz, L.: CitcomS v3.3.1, Zenodo [code], https://doi.org/10.5281/zenodo.7271920, 2022. a
Murer, M., Formica, G., Milicchio, F., Morganti, S., and Auricchio, F.: A coupled multiphase LagrangianEulerian fluiddynamics framework for numerical simulation of Laser Metal Deposition process, The International Journal of Advanced Manufacturing Technology, 120, 3269–3286, https://doi.org/10.1007/s00170022087637, 2022. a
Neuharth, D., Brune, S., Glerum, A., Heine, C., and Welford, J. K.: Formation of continental microplates through rift linkage: Numerical modeling and its application to the Flemish Cap and Sao Paulo Plateau, Geochem. Geophy. Geosy., 22, e2020GC009615, https://doi.org/10.1029/2020GC009615, 2021. a
Popov, A. and Marchevsky, I.: MPIBased PFEM2 Method Solver for ConvectionDominated CFD Problems, in: International Conference on Parallel Computational Technologies, Springer, 261–275, https://doi.org/10.1007/9783031116230_18, 2022. a
Puckett, E. G., Turcotte, D. L., He, Y., Lokavarapu, H., Robey, J. M., and Kellogg, L. H.: New numerical approaches for modeling thermochemical convection in a compositionally stratified fluid, Phys. Earth Planet. Int., 276, 10–35, 2018. a
Pusok, A. E., Kaus, B. J., and Popov, A. A.: On the quality of velocity interpolation schemes for markerincell method and staggered grids, Pure Appl. Geophys., 174, 1071–1089, 2017. a
Richards, F., Hoggard, M., Cowton, L., and White, N.: Reassessing the thermal structure of oceanic lithosphere with revised global inventories of basement depths and heat flow measurements, J. Geophys. Res.Sol. Ea., 123, 9136–9161, 2018. a
Samuel, H.: A deformable particleincell method for advective transport in geodynamic modelling, Geophys. J. Int., 214, 1744–1773, 2018. a
Schierjott, J. C., Thielmann, M., Rozel, A. B., Golabek, G. J., and Gerya, T. V.: Can grain size reduction initiate transform faults? Insights from a 3D numerical study, Tectonics, 39, e2019TC005793, https://doi.org/10.1029/2019TC005793, 2020. a
Schmid, D. W. and Podladchikov, Y. Y.: Analytical solutions for deformable elliptical inclusions in general shear, Geophys. J. Int., 155, 269–288, 2003. a
Schubert, G., Turcotte, D. L., and Olson, P.: Mantle Convection in the Earth and Planets, Part 1, Cambridge, ISBN13: 9780521353670, 2001. a
Sime, N., Maljaars, J. M., Wilson, C. R., and van Keken, P. E.: An exactly mass conserving and pointwise divergence free velocity method: Application to compositional buoyancy driven flow problems in geodynamics, Geochem. Geophy. Geosy., 22, e2020GC009349, https://doi.org/10.1029/2020GC009349, 2021. a, b, c
Sime, N., Wilson, C. R., and van Keken, P. E.: A pointwise conservative method for thermochemical convection under the compressible anelastic liquid approximation, Geochem. Geophy. Geosy., 23, e2021GC009922, https://doi.org/10.1029/2021GC009922, 2022. a
Stein, C. A. and Stein, S.: A model for the global variation in oceanic depth and heat flow with lithospheric age, Nature, 359, 123–129, 1992. a
Strauss, W. A.: Partial differential equations: An introduction, John Wiley & Sons, ISBN13 9780470054567, 2007. a
Tackley, P. J.: Selfconsistent generation of tectonic plates in threedimensional mantle convection, Earth Planet. Sc. Lett., 157, 9–22, 1998. a
Tackley, P. J.: Modelling compressible mantle convection with large viscosity contrasts in a threedimensional spherical shell using the yinyang grid, Phys. Earth Planet. Int., 171, 7–18, 2008. a
Tackley, P. J. and King, S. D.: Testing the tracer ratio method for modeling active compositional fields in mantle convection simulations, Geochem. Geophy. Geosy., 4, 8302, https://doi.org/10.1029/2001GC000214, 2003. a, b, c
Thielmann, M., May, D. A., and Kaus, B. J. P.: Discretization Errors in the Hybrid Finite Element Particleincell Method, Pure Appl. Geophys., 171, 2165–2184, 2014. a, b
Thieulot, C. and Bangerth, W.: On the choice of finite element for applications in geodynamics, Solid Earth, 13, 229–249, https://doi.org/10.5194/se132292022, 2022. a
Trim, S. J., Butler, S. L., McAdam, S. S., and Spiteri, R. J.: Manufacturing an exact solution for 2D thermochemical mantle convection models, Geochem. Geophy. Geosy., 24, e2022GC010807, https://doi.org/10.1029/2022GC010807, 2023a. a, b
Trim, S. J., Butler, S. L., and Spiteri, R. J.: The impact of velocity update frequency on time accuracy for mantle convection particle methods, Authorea Preprints, https://doi.org/10.22541/essoar.169444235.56582698/v1, 2023b. a
Van Dinther, Y., Gerya, T., Dalguer, L., Mai, P. M., Morra, G., and Giardini, D.: The seismic cycle at subduction thrusts: Insights from seismothermomechanical models, J. Geophys. Res.Sol. Ea., 118, 6183–6202, 2013. a
van Keken, P., King, S., Schmeling, H., Christensen, U., Neumeister, D., and Doin, M.P.: A comparison of methods for the modeling of thermochemical convection, J. Geophys. Res.Sol. Ea., 102, 22477–22495, 1997. a, b, c
Van Zelst, I., Wollherr, S., Gabriel, A.A., Madden, E. H., and van Dinther, Y.: Modeling megathrust earthquakes across scales: Oneway coupling from geodynamics and seismic cycles to dynamic rupture, J. Geophys. Res.Sol. Ea., 124, 11414–11446, 2019. a
Wang, H., Agrusta, R., and van Hunen, J.: Advantages of a conservative velocity interpolation (CVI) scheme for particleincell methods with application in geodynamic modeling, Geochem. Geophy. Geosy., 16, 2015–2023, https://doi.org/10.1002/2015GC005824, 2015. a
Weinberg, R. F. and Schmeling, H.: Polydiapirs: multiwavelength gravity structures, J. Struct. Geol., 14, 425–436, 1992. a, b
Zhong, S.: Analytic solutions for Stokes' flow with lateral variations in viscosity, Geophys. J. Int., 124, 18–28, 1996. a
Zhong, S. and Hager, B. H.: Entrainment of a dense layer by thermal plumes, Geophys. J. Int., 154, 666–676, 2003. a
Zhong, S., McNamara, A., Tan, E., Moresi, L., and Gurnis, M.: A benchmark study on mantle convection in a 3D spherical shell using CitcomS, Geochem. Geophy. Geosy., 9, Q10017, https://doi.org/10.1029/2008GC002048, 2008. a, b
 Abstract
 Introduction
 Governing equations
 Numerical methodology
 Deriving transient benchmark solutions
 Numerical evaluation of particle schemes
 Application: evolution of the mineral grain size below oceanic lithosphere
 Conclusions
 Appendix A: Derivation of an exact solution in an annulus
 Appendix B: Detailed error investigation of the spherical shell benchmark
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 Governing equations
 Numerical methodology
 Deriving transient benchmark solutions
 Numerical evaluation of particle schemes
 Application: evolution of the mineral grain size below oceanic lithosphere
 Conclusions
 Appendix A: Derivation of an exact solution in an annulus
 Appendix B: Detailed error investigation of the spherical shell benchmark
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References