Benchmarking the accuracy of higher-order particle methods in geodynamic models of transient flow

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 advection-dominated 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 crustal-scale 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 time-dependent 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 higher-order 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 higher-order algorithms matters for geodynamic applications with an example of modeling small-scale convection underneath an oceanic plate and show that the predicted place and time of onset of small-scale 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 open-source code ASPECT will allow geodynamicists to investigate complex time-dependent geodynamic processes such as elastic deformation, anisotropic fabric development, melt generation and migration, and grain damage.

1 Introduction

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 cross-section 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 finite-element (e.g., Moresi et al., 2022), finite-difference (e.g., Gerya and Yuen, 2003; Kaus et al., 2016), and finite-volume 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, marker-chain methods, or volume-of-fluid methods.

Due to their inherent suitability for modeling advection-dominated 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 field-based 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 finite-element 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 visco-plastic or stress-dependent 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 time-dependent 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 steady-state 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 small-scale 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.

2 Governing equations

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:

$$\begin{array}{}\text{(1)}& {\displaystyle }-\mathrm{\nabla }\cdot \left(\mathrm{2}\mathit{\eta }\dot{\mathit{\epsilon }}\left(\mathit{u}\right)\right)+\mathrm{\nabla }p& {\displaystyle }=\mathit{\rho }\mathit{g},\text{(2)}& {\displaystyle }\mathrm{\nabla }\cdot \mathit{u}& {\displaystyle }=\mathrm{0},\end{array}$$

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 particle-based 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 non-negligible diffusivity and reactions have been described elsewhere (Gerya and Yuen, 2003; Sime et al., 2022). In other words, we consider equations of the form

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial {\mathit{\varphi }}_c}{\partial t}}+\mathit{u}\cdot \mathrm{\nabla }{\mathit{\varphi }}_c=H_c,\end{array}$$

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 time-dependent 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.

3 Numerical methodology

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 high-performance 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 field-based quantities. Its methods are integrated into the general purpose open-source finite-element 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 finite-element method, time-stepping 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:

$$\begin{array}{}\text{(4)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathit{x}}_i\left(t\right)=\mathit{u}\left({\mathit{x}}_i\right(t),t),\text{(5)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\mathit{\varphi }}_{c,i}\left(t\right)=H_c\left({\mathit{x}}_i\right(t),t,{\mathit{\varphi }}_{c,i}(t\left)\right).\end{array}$$

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ⁿ and tⁿ⁺¹. If we denote this interpolation in time by ${\stackrel{\mathrm{̃}}{\mathit{u}}}_h(\mathit{x},t)$ where ${\stackrel{\mathrm{̃}}{\mathit{u}}}_h(\mathit{x},t^n)={\mathit{u}}_h^n\left(\mathit{x}\right)$, then the equation being solved is actually

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}{\stackrel{\mathrm{̃}}{\mathit{x}}}_i\left(t\right)={\stackrel{\mathrm{̃}}{\mathit{u}}}_h\left({\stackrel{\mathrm{̃}}{\mathit{x}}}_i\right(t),t),\end{array}$$

where ${\stackrel{\mathrm{̃}}{\mathit{x}}}_i\left(t\right)$ is the exact solution of this equation using the “wrong” velocity field. If ${\stackrel{\mathrm{̃}}{\mathit{u}}}_h$ is a good approximation to u, then we hope that $\stackrel{\mathrm{̃}}{\mathit{x}}\left(t\right)$ is a good approximation of x(t). In practice, however, we can not even compute $\stackrel{\mathrm{̃}}{\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{̃}}{\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

$$\begin{array}{}\text{(7)}& {\displaystyle }\Vert {\stackrel{\mathrm{̃}}{\mathit{x}}}_{i,h}\left(t\right)-{\stackrel{\mathrm{̃}}{\mathit{x}}}_i\left(t\right)\Vert \le C\left(t\right)\mathrm{\Delta }{t}_{\mathrm{p}}^q,\end{array}$$

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{̃}}{\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{̃}}{\mathit{x}}}_{i,h}\left(t\right)\Vert$. This quantity will, in the best case, only satisfy an estimate of the form

$$\begin{array}{r}\Vert {\stackrel{\mathrm{̃}}{\mathit{x}}}_{i,h}\left(t\right)-{\mathit{x}}_i\left(t\right)\Vert =∥\left({\stackrel{\mathrm{̃}}{\mathit{x}}}_{i,h}\left(t\right)-{\stackrel{\mathrm{̃}}{\mathit{x}}}_i\left(t\right)\right)+\left({\stackrel{\mathrm{̃}}{\mathit{x}}}_i\left(t\right)-{\mathit{x}}_i\left(t\right)\right)∥\\ \le ∥{\stackrel{\mathrm{̃}}{\mathit{x}}}_{i,h}\left(t\right)-{\stackrel{\mathrm{̃}}{\mathit{x}}}_i\left(t\right)∥+∥{\stackrel{\mathrm{̃}}{\mathit{x}}}_i\left(t\right)-{\mathit{x}}_i\left(t\right)∥\\ \le C_{\mathrm{1}}\left(t\right)\mathrm{\Delta }{t}_{\mathrm{p}}^q+C_{\mathrm{2}}\left(t\right)\Vert \mathit{u}-{\mathit{u}}_h\Vert +C_{\mathrm{3}}\left(t\right)\Vert {\mathit{u}}_h-{\stackrel{\mathrm{̃}}{\mathit{u}}}_h\Vert ,\end{array}$$

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 time-step 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{̃}}{\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ⁿ⁺¹ and are now updating particle locations from xⁿ to xⁿ⁺¹. In cases where one wants to solve for particle locations before updating the velocity field, ${\stackrel{\mathrm{̃}}{\mathit{u}}}_h$ can be extrapolated beyond tⁿ 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:

1. Forward Euler (FE) – the simplest method often used is the forward Euler scheme,

   $${\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{x}}}^{n+\mathrm{1}}={\stackrel{\mathrm{̃}}{\mathit{x}}}^n+\mathrm{\Delta }t{\stackrel{\mathrm{̃}}{\mathit{u}}}_h(t^n,{\stackrel{\mathrm{̃}}{\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 already-computed time point.

2. Runge–Kutta second order (RK2) – accuracy and stability can be improved by using a second-order Runge–Kutta scheme (that is, q=2 in Eq. (7)). Among the many second-order 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{̃}}{\mathit{u}}}_h(t^n,{\stackrel{\mathrm{̃}}{\mathit{x}}}^n),\\ {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{x}}}^{n+\mathrm{1}}& {\displaystyle }={\mathit{x}}^n+\mathrm{\Delta }t{\stackrel{\mathrm{̃}}{\mathit{u}}}_h\left(t^n+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}},{\stackrel{\mathrm{̃}}{\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.

3. Runge–Kutta second-order space, first-order 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ⁿ and tⁿ⁺¹. 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{̃}}{\mathit{u}}}_h(t^n,{\stackrel{\mathrm{̃}}{\mathit{x}}}^n),\\ {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{x}}}^{n+\mathrm{1}}& {\displaystyle }={\stackrel{\mathrm{̃}}{\mathit{x}}}^n+\mathrm{\Delta }t{\stackrel{\mathrm{̃}}{\mathit{u}}}_h\left(t^n,{\stackrel{\mathrm{̃}}{\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ⁿ instead of $t^n+\frac{\mathrm{\Delta }t}{\mathrm{2}}$ but at the correct location ${\stackrel{\mathrm{̃}}{\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.

4. Runge–Kutta fourth order (RK4) – a further improvement in particle advection can be achieved by a fourth-order 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{\stackrel{\mathrm{̃}}{\mathit{u}}}_h\left(t^n,{\stackrel{\mathrm{̃}}{\mathit{x}}}^n\right),\\ {\displaystyle }{\mathit{k}}_{\mathrm{2}}& {\displaystyle }={\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}}{\stackrel{\mathrm{̃}}{\mathit{u}}}_h\left(t^n+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}},{\stackrel{\mathrm{̃}}{\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}}}{\stackrel{\mathrm{̃}}{\mathit{u}}}_h\left(t^n+{\displaystyle \frac{\mathrm{\Delta }t}{\mathrm{2}}},{\stackrel{\mathrm{̃}}{\mathit{x}}}^n+{\displaystyle \frac{{\mathit{k}}_{\mathrm{2}}}{\mathrm{2}}}\right),\\ {\displaystyle }{\mathit{k}}_{\mathrm{4}}& {\displaystyle }=\mathrm{\Delta }t{\stackrel{\mathrm{̃}}{\mathit{u}}}_h\left(t^{n+\mathrm{1}},{\stackrel{\mathrm{̃}}{\mathit{x}}}^n+{\mathit{k}}_{\mathrm{3}}\right),\\ {\displaystyle }{\stackrel{\mathrm{̃}}{\mathit{x}}}^{n+\mathrm{1}}& {\displaystyle }={\stackrel{\mathrm{̃}}{\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 finite-element 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ⁿ to $t^n+\frac{\mathrm{\Delta }t}{\mathrm{2}}$ and tⁿ⁺¹. 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 second-order accurate in time.

4 Deriving transient benchmark solutions

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

$$\begin{array}{}\text{(8)}& {\displaystyle }& {\displaystyle }-\mathrm{\nabla }\cdot \left(\mathrm{2}\dot{\mathit{\epsilon }}\left(\mathit{u}\right)\right)+\mathrm{\nabla }p=\mathit{\rho }\mathit{g},\text{(9)}& {\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0},\text{(10)}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}+\mathit{u}\cdot \mathrm{\nabla }\mathit{\rho }=\mathrm{0}.\end{array}$$

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 time-independent 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 right-hand-side 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 steady-state flow field – defined as a velocity field u that does not change over time – the right-hand-side 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 right-hand-side ρ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 right-hand-side ρg will not change even if ρ is advected with the flow. Choosing ρ=Ψ is therefore an easy approach to guarantee a steady-state flow field.

The benchmarks below extend these steady-state 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 time-dependent 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 steady-state stream function Ψ. Second, we choose the time-variable 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 (time-dependent) 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 time-dependent 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

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }{v}_r(r,\mathit{\theta },t)=g\left(r\right)k\sin \left(k\right(\mathit{\theta }-\mathit{\tau }\left(t\right)\left)\right),\text{(12)}& {\displaystyle }& {\displaystyle }{v}_{\mathit{\theta }}(r,\mathit{\theta },t)=f\left(r\right)\cos \left(k\right(\mathit{\theta }-\mathit{\tau }\left(t\right)\left)\right)+r\mathit{\omega }\left(t\right),\text{(13)}& {\displaystyle }& {\displaystyle }p(r,\mathit{\theta },t)=kh\left(r\right)\sin \left(k\right(\mathit{\theta }-\mathit{\tau }\left(t\right)\left)\right),\text{(14)}& {\displaystyle }& {\displaystyle }\mathit{\rho }(r,\mathit{\theta },t)=-\left({\displaystyle \frac{A}{\mathrm{2}}}{r}^{\mathrm{2}}+B\ln \left(r\right)-\mathrm{1}\right)\cos \left(k\right(\mathit{\theta }-\mathit{\tau }\left(t\right)\left)\right),\text{(15)}& {\displaystyle }& {\displaystyle }{g}_r(r,\mathit{\theta },t)=m\left(r\right)k{\displaystyle \frac{\sin \left(k\right(\mathit{\theta }-\mathit{\tau }\left(t\right)\left)\right)}{\mathit{\rho }(r,\mathit{\theta },t)}},\text{(16)}& {\displaystyle }& {\displaystyle }{g}_{\mathit{\theta }}=\mathrm{0}.\end{array}$$

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 time-dependent 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)\mathrm{d}s$. The spherical shell has an inner radius of R₁=1 and an outer radius of R₂=2. The setup of the benchmark and a snapshot of the solution is shown in Fig. 1.

Figure 1 Solution of a transient spherical shell benchmark. (a) The density field of the benchmark at t=0. (b) Velocity solution at t=0. (c, d) Initial (t=0) and later (t=2.6) particle distributions after almost two full rotations of the model. Particles are colored based on a unique index given to each particle at the beginning.

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 density-driven 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

$$\begin{array}{}\text{(17)}& {\displaystyle }{v}_x(x,y,t)& {\displaystyle }=\sin \left(\mathit{\pi }\right(x-\mathit{\tau }\left(t\right)\left)\right)\cos \left(\mathit{\pi }y\right)+\mathit{\omega }\left(t\right),\text{(18)}& {\displaystyle }{v}_y(x,y,t)& {\displaystyle }=-\cos \left(\mathit{\pi }\right(x-\mathit{\tau }\left(t\right)\left)\right)\sin \left(\mathit{\pi }y\right),\text{(19)}& {\displaystyle }p(x,y,t)& {\displaystyle }=\mathrm{2}\mathit{\pi }\cos \left(\mathit{\pi }\right(x-\mathit{\tau }\left(t\right)\left)\right)\cos \left(\mathit{\pi }y\right),\text{(20)}& {\displaystyle }\mathit{\rho }(x,y,t)& {\displaystyle }=\sin \left(\mathit{\pi }\right(x-\mathit{\tau }\left(t\right)\left)\right)\sin \left(\mathit{\pi }y\right),\text{(21)}& {\displaystyle }{g}_x(x,y,t)& {\displaystyle }=\mathrm{0},\text{(22)}& {\displaystyle }{g}_y(x,y,t)& {\displaystyle }=-\mathrm{4}{\mathit{\pi }}^{\mathrm{2}}{\displaystyle \frac{\cos \left(\mathit{\pi }\right(x-\mathit{\tau }\left(t\right)\left)\right)\sin \left(\mathit{\pi }y\right)}{\mathit{\rho }(x,y,t)}}.\end{array}$$

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 steady-state flow with a time-dependent velocity amplitude.

Figure 2 Solution of a transient box benchmark with known analytic solution. (a) The density field of the benchmark. (b) Velocity solution. (c, d) Particle distributions at model start (t=0) and at t=0.5. Particles are colored as in Fig. 1.

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:

1. 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.

2. 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₂) finite elements with a penalty method as described in He et al. (2016).

3. Same as case 2, but we use continuous, piecewise quadratic (Q₂) 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) time-stepping algorithm that is second-order accurate in time to solve the advection Eq. (10).

4. We use the exact density as the initial condition for particles whose position we advect using a second-order 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₂ discontinuous finite-element field and evaluate that at quadrature points as necessary.

5. 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 finite-element and particle interpolation algorithm in this paper. We will use a Q₂×Q₁ Stokes element (continuous, piecewise quadratic velocity, and continuous, piecewise linear pressure) and a linear least-squares 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 least-squares 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).

5 Numerical evaluation of particle schemes

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₂ 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).

Figure 3 Transient spherical annulus benchmark. (a, c, e) L₂ error norms of velocity (a, b), pressure (c, d), and density (e, f) for different cell sizes h at time $t=\ln (\mathrm{1}+\mathrm{4}\mathit{\pi })\approx \mathrm{2.6075}$. Different colors and marker styles show different advection methods; gray lines show ideal first-, second-, and third-order convergence. Note that the line for the exact (benchmark) density overlaps with the RK2 line. (b, d, f) L₂ error norms of velocity (a, b), pressure (c, d), and density (e, f) over time for resolution $h=\mathrm{1}/\mathrm{128}$. Colors as in left column, and the exact benchmark density line is hidden behind the RK2 case. (For more details on the distinction between the RK2 case and the benchmark density case see the Appendix.)

The left column illustrates that all advection methods but RK2FOT reach second-order convergence for the density with increasing resolution (bottom left panel). As expected, RK2FOT is limited by the available time information and only reaches first-order convergence. An additional detail is that the field methods (Q₂ and DGQ₂) 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 first-order convergence rate in velocity, while Q₂ and DGQ₂ 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 third-order 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₂ and DGQ₂. 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₂ and DGQ₂ behave distinctly different. DGQ₂ starts at a much smaller error value than all other methods, but it accumulates significant errors towards the end of the model run. Q₂ 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 first-order 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, low-order particle methods show larger errors than the tested field methods, while higher-order 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 time-stepping scheme allow for higher-order accuracy), a higher-order 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 third-order convergence in velocity and second-order 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 first-order time (RK2FOT) advection method creates a first-order accurate approximation for the density, which also generates a first-order accurate pressure and velocity solution, thereby significantly limiting the potential accuracy of the Stokes elements. All other advection methods reach second-order 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₂ and DGQ₂ finite-element 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 third-order convergence rate at low resolutions and transitions to a second-order 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 upside-down 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).

Figure 4 Transient box benchmark. (a, c, e) L₂ error norms of velocity (a, b), pressure (c, d), and density (e, f) for different cell sizes h at time $t=\ln (\mathrm{1}+\mathrm{2})\approx \mathrm{1.0986}$. Different colors and marker styles show different advection methods; gray lines show ideal first-, second-, and third-order convergence. (b, d, f) L₂ error norms of velocity (a, b), pressure (c, d), and density (e, f) as a function of time for resolution $h=\mathrm{1}/\mathrm{128}$.

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₂ and DGQ₂ 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 time-dependent geodynamic models. Both the Q₂ and DGQ₂ 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 finite-element advection methods, the first-order 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, first-order particle methods yield larger errors than the tested field-based methods, while higher-order particle methods outperform the investigated field-based 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 time-stepping scheme allow for higher-order accuracy), a higher-order 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 higher-order finite elements can break down to the same rate as for first-order elements (e.g., Thielmann et al., 2014). However, solutions are rarely fully discontinuous and instead contain a mix of smooth and non-smooth regions. Additionally, despite showing the same convergence rate, higher-order elements have still delivered higher accuracy in absolute terms than lower-order methods in many benchmarks (Kronbichler et al., 2012; Thieulot and Bangerth, 2022). We speculate that the same results would be seen for higher-order advection methods in time, although the construction of appropriate benchmarks would be challenging.

Finally, the improved accuracy of higher-order 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 matrix-based 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 matrix-free solvers with short-recurrence 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}}\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{\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.

Table 1 Performance comparison of RK2 and RK2FOT for the box benchmark with a shortened end time of t=0.1 in dependence of resolution. The table shows the total model runtime (Model) and particle advection time (RK2 and RK2FOT) for the two advection algorithms, as well as the fraction of total model runtime used for particle advection (RK2/Model and RK2FOT/Model) and the relative performance cost of RK2 vs RK2FOT. Note that parts of the algorithm, such as sorting particles into mesh cells or interpolating properties back to the mesh, are not listed because their cost is independent of the advection algorithm.

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 first-order method like RK2FOT in applications that depend on accurate solutions.

6 Application: evolution of the mineral grain size below oceanic lithosphere

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 low-order 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 first-order 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⁻¹ 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⁻¹ at the bottom of the model. The left (ridge axis) boundary of the model is closed and stress-free (free slip). The vertical velocity component is not prescribed at the bottom and right boundaries; instead, a depth-dependent 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 half-space 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 small-scale 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

$$\begin{array}{}\text{(23)}& {\displaystyle }{\mathit{\eta }}_{\mathrm{diff}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}_{\mathrm{diff}}^{-\mathrm{1}}{d}^m\exp \left({\displaystyle \frac{E_{\mathrm{diff}}+P{V}_{\mathrm{diff}}}{RT}}\right),\text{(24)}& {\displaystyle }{\mathit{\eta }}_{\mathrm{dis}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}_{\mathrm{dis}}^{-\frac{\mathrm{1}}{n}}{\dot{\mathit{\epsilon }}}_{\mathrm{dis},\mathrm{II}}^{\frac{\mathrm{1}-n}{n}}\exp \left({\displaystyle \frac{E_{\mathrm{dis}}+P{V}_{\mathrm{dis}}}{nRT}}\right),\end{array}$$

resulting in an effective viscosity of

$$\begin{array}{}\text{(25)}& {\displaystyle }{\mathit{\eta }}_{\mathrm{eff}}{\displaystyle }={\displaystyle \frac{{\mathit{\eta }}_{\mathrm{diff}}{\mathit{\eta }}_{\mathrm{dis}}}{{\mathit{\eta }}_{\mathrm{diff}}+{\mathit{\eta }}_{\mathrm{dis}}}}.\end{array}$$

Here, η_diff, η_dis, and η_eff are the diffusion, dislocation, and effective viscosity, respectively. $A_{\mathrm{diff}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{15}}$ m³ Pa⁻¹ s⁻¹ and $A_{\mathrm{dis}}={\mathrm{10}}^{-\mathrm{15}}$ Pa^−3.5 s⁻¹ are diffusion and dislocation creep prefactors, E_diff=375 kJ mol⁻¹ and E_dis=530 kJ mol⁻¹ the activation energy for diffusion and dislocation creep, $V_{\mathrm{diff}}=\mathrm{4}\times {\mathrm{10}}^{-\mathrm{6}}$ m³ mol⁻¹ and $V_{\mathrm{dis}}=\mathrm{1.4}\times {\mathrm{10}}^{-\mathrm{5}}$ m³ mol⁻¹ 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¹⁶ and 10²³ 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

$$\begin{array}{}\text{(26)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}d\left(t\right)& =p_{\mathrm{g}}^{-\mathrm{1}}{d}^{\mathrm{1}-p_{\mathrm{g}}}{k}_{\mathrm{g}}\exp \left(-{\displaystyle \frac{E_{\mathrm{g}}+P{V}_{\mathrm{g}}}{RT}}\right)\\ & -\mathrm{4}{\dot{\mathit{\epsilon }}}_{\mathrm{II}}{\dot{\mathit{\epsilon }}}_{\mathrm{dis},\mathrm{II}}{\mathit{\eta }}_{\mathrm{eff}}{\displaystyle \frac{\mathit{\lambda }{d}^{\mathrm{2}}}{c\mathit{\gamma }}},\end{array}\end{array}$$

which implies that the grain size field d(x,t) satisfies the equation

$$\begin{array}{}\text{(27)}& \begin{array}{rl}\left({\displaystyle \frac{\partial d}{\partial t}}+\mathit{u}\cdot \mathrm{\nabla }d\right)& =p_{\mathrm{g}}^{-\mathrm{1}}{d}^{\mathrm{1}-p_{\mathrm{g}}}{k}_{\mathrm{g}}\exp \left(-{\displaystyle \frac{E_{\mathrm{g}}+P{V}_{\mathrm{g}}}{RT}}\right)\\ & -\mathrm{4}{\dot{\mathit{\epsilon }}}_{\mathrm{II}}{\dot{\mathit{\epsilon }}}_{\mathrm{dis},\mathrm{II}}{\mathit{\eta }}_{\mathrm{eff}}{\displaystyle \frac{\mathit{\lambda }{d}^{\mathrm{2}}}{c\mathit{\gamma }}}.\end{array}\end{array}$$

Here, $k_{\mathrm{g}}=\mathrm{1.92}\times {\mathrm{10}}^{-\mathrm{10}}$ m³ s⁻¹ is the grain size growth prefactor, and E_g=400 kJ mol⁻¹ and V_g=0 m³ mol⁻¹ 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 right-hand 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 particle-based 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.

Figure 5 (a) Setup of the application model. Background colors illustrate temperature and solid black lines are streamlines. Arrows indicate velocity. (b, c) Grain size at the end of the application model after 200 million years for models using RK2 (b), and RK2FOT (c) for different time step lengths (rows). The onset of convection below the plate is marked for each model and measured as the distance from the ridge to the −50 K non-adiabatic temperature contour at 150 km depth. The corresponding age of the oceanic plate at this distance is also shown.

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 small-scale 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 small-scale 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.

Table 2 Onset of convection for RK2 and RK2FOT for the application model with a CFL number of 0.15 over time. The table shows distance from ridge and plate age of onset of convection listed for different model times for both RK2 and RK2FOT.

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 small-scale 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 small-scale 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 small-scale 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 lower-order advection schemes may need to be adjusted or reproduced in higher-resolution studies.

7 Conclusions

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 time-variable 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 steady-state 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 finite-element models already use Stokes elements that allow for higher-order accuracy to ensure stability (e.g., Taylor–Hood Q₂×Q₁ or $Q_{\mathrm{2}}\times P_{-\mathrm{1}}$), it would be straightforward to extend their particle advection algorithms to a second-order-in-time 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 open-source community software ASPECT and hope that they are useful to the community at large.

Appendix A: Derivation of an exact solution in an annulus

Extending our previously published spherical benchmark (Gassmöller et al., 2019) seems to be straightforward by adding a time-dependent solid body rotation to the existing solution. However, because our earlier solution is already a purely rotational flow, an additional time-dependent rotation does not create a transient solution for the density, not allowing to intuitively measure the accuracy of the combined particle–finite-element 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₁, and the outer radius by R₂. Furthermore, we assume that the viscosity is a constant η=1 and set the gravity vector to an inward-pointing 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)

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{{\partial }^{\mathrm{2}}{v}_r}{\partial r^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial v_r}{\partial r}}+{\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{v}_r}{\partial {\mathit{\theta }}^{\mathrm{2}}}}-{\displaystyle \frac{v_r}{r^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{2}}{r^{\mathrm{2}}}}{\displaystyle \frac{\partial v_{\mathit{\theta }}}{\partial \mathit{\theta }}}\\ & -{\displaystyle \frac{\partial p}{\partial r}}-\mathit{\rho }{g}_r=\mathrm{0}\end{array},\text{(A2)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}{v}_{\mathit{\theta }}}{\partial r^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial v_{\mathit{\theta }}}{\partial r}}+{\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{v}_{\mathit{\theta }}}{\partial {\mathit{\theta }}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{2}}{r^{\mathrm{2}}}}{\displaystyle \frac{\partial v_r}{\partial \mathit{\theta }}}-{\displaystyle \frac{v_{\mathit{\theta }}}{r^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial p}{\partial \mathit{\theta }}}=\mathrm{0},\text{(A3)}& {\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \left(r{v}_r\right)}{\partial r}}+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial v_{\mathit{\theta }}}{\partial \mathit{\theta }}}=\mathrm{0}.\end{array}$$

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

$$\begin{array}{}\text{(A4)}& v_{\mathit{\theta }}(r,\mathit{\theta })=f\left(r\right)\cos \left(k\mathit{\theta }\right),\end{array}$$

where k is an integer and where f(r) will be specified later. From Eq. (A3) we then obtain

$$\begin{array}{}\text{(A5)}& {\displaystyle \frac{\partial \left(r{v}_r\right)}{\partial r}}=-{\displaystyle \frac{\partial v_{\mathit{\theta }}}{\partial \mathit{\theta }}}=kf\left(r\right)\sin \left(k\mathit{\theta }\right),\end{array}$$

leading to

$$\begin{array}{}\text{(A6)}& v_r(r,\mathit{\theta })=g\left(r\right)k\sin \left(k\mathit{\theta }\right),\end{array}$$

where

$$\begin{array}{}\text{(A7)}& g\left(r\right)={\displaystyle \frac{\mathrm{1}}{r}}\stackrel{r}{\int }f\left(r^{\prime }\right)\mathrm{d}{r}^{\prime }.\end{array}$$

Since we want to fix the velocity to be tangential at both boundaries we have

$$\begin{array}{}\text{(A8)}& v_r(R_{\mathrm{1}},\mathit{\theta })=v_r(R_{\mathrm{2}},\mathit{\theta })=\mathrm{0}\end{array}$$

for all $\mathit{\theta }\in [\mathrm{0},\mathrm{2}\mathit{\pi }]$. We choose

$$\begin{array}{}\text{(A9)}& f\left(r\right)=Ar+B/r\end{array}$$

in analogy to the solution of the Laplace equation in Chap. 6 of Strauss (2007), and thus

$$\begin{array}{}\text{(A10)}& g\left(r\right)={\displaystyle \frac{A}{\mathrm{2}}}r+{\displaystyle \frac{B}{r}}\ln r+{\displaystyle \frac{C}{r}},\end{array}$$

where C is a non-zero constant of integration. Given the boundary conditions in Eq. (A8) we find that

$$\begin{array}{}\text{(A11)}& {\displaystyle }A=-C{\displaystyle \frac{\mathrm{2}(\ln R_{\mathrm{1}}-\ln R_{\mathrm{2}})}{R_{\mathrm{2}}^{\mathrm{2}}\ln R_{\mathrm{1}}-R_{\mathrm{1}}^{\mathrm{2}}\ln R_{\mathrm{2}}}},\text{(A12)}& {\displaystyle }B=-C{\displaystyle \frac{R_{\mathrm{2}}^{\mathrm{2}}-R_{\mathrm{1}}^{\mathrm{2}}}{R_{\mathrm{2}}^{\mathrm{2}}\ln R_{\mathrm{1}}-R_{\mathrm{1}}^{\mathrm{2}}\ln R_{\mathrm{2}}}}.\end{array}$$

In this work we choose $C=-\mathrm{1}$. Our earlier choice of f means that

$$\begin{array}{}\text{(A13)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}f}{\partial r^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial f}{\partial r}}-{\displaystyle \frac{f}{r^{\mathrm{2}}}}=\mathrm{0},\end{array}$$

so that Eq. (A2) simplifies to

$$\begin{array}{}\text{(A14)}& {\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{v}_{\mathit{\theta }}}{\partial {\mathit{\theta }}^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{2}}{r^{\mathrm{2}}}}{\displaystyle \frac{\partial v_r}{\partial \mathit{\theta }}}-{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial p}{\partial \mathit{\theta }}}=\mathrm{0},\end{array}$$

which together with Eq. (A4) leads to

$$\begin{array}{}\text{(A15)}& p(r,\mathit{\theta })=kh\left(r\right)\sin \left(k\mathit{\theta }\right)+l\left(r\right),\end{array}$$

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

$$\begin{array}{}\text{(A16)}& \begin{array}{rl}\mathit{\rho }(r,\mathit{\theta })& ={\displaystyle \frac{{\partial }^{\mathrm{2}}{v}_r}{\partial r^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial v_r}{\partial r}}+{\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{v}_r}{\partial {\mathit{\theta }}^{\mathrm{2}}}}-{\displaystyle \frac{v_r}{r^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{2}}{r^{\mathrm{2}}}}{\displaystyle \frac{\partial v_{\mathit{\theta }}}{\partial \mathit{\theta }}}-{\displaystyle \frac{\partial p}{\partial r}}\\ & =k{g}^{\prime \prime }\left(r\right)\sin \left(k\mathit{\theta }\right)+k{\displaystyle \frac{g^{\prime }\left(r\right)}{r}}\sin \left(k\mathit{\theta }\right)-k^{\mathrm{3}}{\displaystyle \frac{g\left(r\right)}{r^{\mathrm{2}}}}\sin \left(k\mathit{\theta }\right)\\ & -k{\displaystyle \frac{g\left(r\right)}{r^{\mathrm{2}}}}\sin \left(k\mathit{\theta }\right)+k{\displaystyle \frac{\mathrm{2}f\left(r\right)}{r^{\mathrm{2}}}}\sin \left(k\mathit{\theta }\right)\\ & -k{h}^{\prime }\left(r\right)\sin \left(k\mathit{\theta }\right)-l^{\prime }\left(r\right)\\ & =m\left(r\right)k\sin \left(k\mathit{\theta }\right)-l^{\prime }\left(r\right)\end{array}\end{array}$$

with

$$\begin{array}{}\text{(A17)}& \begin{array}{rl}m\left(r\right)& =g^{\prime \prime }-{\displaystyle \frac{g^{\prime }}{r}}-{\displaystyle \frac{g}{r^{\mathrm{2}}}}(k^{\mathrm{2}}-\mathrm{1})+{\displaystyle \frac{f}{r^{\mathrm{2}}}}+{\displaystyle \frac{f^{\prime }}{r}}\\ & =-{\displaystyle \frac{f-g}{r^{\mathrm{2}}}}+{\displaystyle \frac{f^{\prime }-g^{\prime }}{r}}-{\displaystyle \frac{f-g}{r^{\mathrm{2}}}}-{\displaystyle \frac{g}{r^{\mathrm{2}}}}(k^{\mathrm{2}}-\mathrm{1})\\ & +{\displaystyle \frac{f}{r^{\mathrm{2}}}}+{\displaystyle \frac{f^{\prime }}{r}}\\ & =-{\displaystyle \frac{f-g}{r^{\mathrm{2}}}}+{\displaystyle \frac{f^{\prime }}{r}}-{\displaystyle \frac{f-g}{r^{\mathrm{2}}}}-{\displaystyle \frac{f-g}{r^{\mathrm{2}}}}\\ & -{\displaystyle \frac{g}{r^{\mathrm{2}}}}(k^{\mathrm{2}}-\mathrm{1})+{\displaystyle \frac{f}{r^{\mathrm{2}}}}+{\displaystyle \frac{f^{\prime }}{r}}\\ & =-\mathrm{3}{\displaystyle \frac{f-g}{r^{\mathrm{2}}}}+{\displaystyle \frac{f^{\prime }}{r}}-{\displaystyle \frac{g}{r^{\mathrm{2}}}}(k^{\mathrm{2}}-\mathrm{1})+{\displaystyle \frac{f}{r^{\mathrm{2}}}}+{\displaystyle \frac{f^{\prime }}{r}}\end{array}\end{array}$$

since it is easy to verify using Eq. (A7) that $g^{\prime }\left(r\right)=(f-g)/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,

$$\begin{array}{}\text{(A18)}& p(r,\mathit{\theta }){|}_{k=\mathrm{0}}=l\left(r\right)={\mathit{\rho }}_{\mathrm{0}}{g}_r(R_{\mathrm{2}}-r),\end{array}$$

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₂.

In summary, Eqs. (A4), (A6), and (A15) form a solution of the incompressible Stokes equations, which fully stated reads

$$\begin{array}{}\text{(A19)}& {\displaystyle }{v}_{\mathit{\theta }}(r,\mathit{\theta })& {\displaystyle }=f\left(r\right)\cos \left(k\mathit{\theta }\right),\text{(A20)}& {\displaystyle }{v}_r(r,\mathit{\theta })& {\displaystyle }=g\left(r\right)k\sin \left(k\mathit{\theta }\right),\text{(A21)}& {\displaystyle }p(r,\mathit{\theta })& {\displaystyle }=kh\left(r\right)\sin \left(k\mathit{\theta }\right)+{\mathit{\rho }}_{\mathrm{0}}{g}_r(R_{\mathrm{2}}-r),\end{array}$$

with

$$\begin{array}{}\text{(A22)}& {\displaystyle }\mathit{\rho }(r,\mathit{\theta })=m\left(r\right)k\sin \left(k\mathit{\theta }\right)+{\mathit{\rho }}_{\mathrm{0}},\text{(A23)}& {\displaystyle }{g}_r=\mathrm{1},\text{(A24)}& {\displaystyle }{g}_{\mathit{\theta }}=\mathrm{0},\text{(A25)}& {\displaystyle }f\left(r\right)=Ar+B/r,\text{(A26)}& {\displaystyle }g\left(r\right)={\displaystyle \frac{A}{\mathrm{2}}}r+{\displaystyle \frac{B}{r}}\ln r+{\displaystyle \frac{C}{r}},\text{(A27)}& {\displaystyle }h\left(r\right)={\displaystyle \frac{\mathrm{2}g\left(r\right)-f\left(r\right)}{r}},\text{(A28)}& {\displaystyle }\begin{array}{rl}m\left(r\right)& =g^{\prime \prime }\left(r\right)-{\displaystyle \frac{g^{\prime }\left(r\right)}{r}}-{\displaystyle \frac{g\left(r\right)}{r^{\mathrm{2}}}}(k^{\mathrm{2}}-\mathrm{1})\\ & +{\displaystyle \frac{f\left(r\right)}{r^{\mathrm{2}}}}+{\displaystyle \frac{f^{\prime }\left(r\right)}{r}},\end{array}\text{(A29)}& {\displaystyle }A=-C{\displaystyle \frac{\mathrm{2}(\ln R_{\mathrm{1}}-\ln R_{\mathrm{2}})}{R_{\mathrm{2}}^{\mathrm{2}}\ln R_{\mathrm{1}}-R_{\mathrm{1}}^{\mathrm{2}}\ln R_{\mathrm{2}}}},\text{(A30)}& {\displaystyle }B=-C{\displaystyle \frac{R_{\mathrm{2}}^{\mathrm{2}}-R_{\mathrm{1}}^{\mathrm{2}}}{R_{\mathrm{2}}^{\mathrm{2}}\ln R_{\mathrm{1}}-R_{\mathrm{1}}^{\mathrm{2}}\ln R_{\mathrm{2}}}},\text{(A31)}& {\displaystyle }C=-\mathrm{1}.\end{array}$$

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:

$$\begin{array}{}\text{(A32)}& \mathrm{\Psi }(r,\mathit{\theta })=-\left({\displaystyle \frac{A}{\mathrm{2}}}{r}^{\mathrm{2}}+B\ln \left(r\right)+C\right)\cos \left(k\mathit{\theta }\right).\end{array}$$

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 steady-state variant of the benchmark and add a known time-dependent 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 right-hand-side 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\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 steady-state solution therefore is the same as above, except

$$\begin{array}{}\text{(A33)}& {\displaystyle }\mathit{\rho }(r,\mathit{\theta })& {\displaystyle }=\mathrm{\Psi }(r,\mathit{\theta }),\text{(A34)}& {\displaystyle }{g}_r(r,\mathit{\theta })& {\displaystyle }=m\left(r\right)k\sin \left(k\mathit{\theta }\right)/\mathit{\rho }(r,\mathit{\theta }),\end{array}$$

with all other constants chosen as before.

In order to transform this steady-state benchmark into a known transient solution, we then add a solid body rotation with a nonlinear time-dependent 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

$$\begin{array}{}\text{(A35)}& v_{\mathit{\theta }}(r,\mathit{\theta },t)=f\left(r\right)\cos \left(k\right(\mathit{\theta }-\mathit{\tau }\left(t\right)\left)\right)+r\mathit{\omega }\left(t\right).\end{array}$$

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:

$$\begin{array}{}\text{(A36)}& \mathit{\tau }\left(t\right)=\underset{\mathrm{0}}{\overset{t}{\int }}\mathit{\omega }\left(s\right)\mathrm{d}s.\end{array}$$

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.

Appendix B: Detailed error investigation of the spherical shell benchmark

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₂ norm but instead the relative difference in error between RK2 and analytical density; i.e., if ${\mathit{ϵ}}_{\mathrm{RK}\mathrm{2}}=\Vert \mathit{u}-{\mathit{u}}_h^{\mathrm{RK}\mathrm{2}}{\Vert }_{L_{\mathrm{2}}}$ and ${\mathit{ϵ}}_{\mathrm{AD}}=\Vert \mathit{u}-{\mathit{u}}_h^{\mathrm{AD}}{\Vert }_{L_{\mathrm{2}}}$, then we plot $({\mathit{ϵ}}_{\mathrm{RK}\mathrm{2}}-{\mathit{ϵ}}_{\mathrm{AD}})/{\mathit{ϵ}}_{\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 higher-order) 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 second-order 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 second-order rate for both analytical density and RK2, leading to a constant relative difference between the two errors.

Figure B1 Transient spherical annulus benchmark. (a, c) L₂ error norms of velocity (a, b) and pressure (c, d) for different cell sizes h at time $t=\ln (\mathrm{1}+\mathrm{4}\mathit{\pi })\approx \mathrm{2.6075}$. The pink line shows the relative difference in the error between RK2 and the analytical density model. Gray lines indicate the same convergence order (dashed line) or one convergence order lower (dash-dotted line) than the analytical density model. (b, d, e) L₂ error norms of velocity (a, b), pressure (c, d), and density (e) as a function of time for resolution $h=\mathrm{1}/\mathrm{128}$ for analytical density and the RK2 model.

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₂, and DGQ₂), 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=\ln (\mathrm{1}+\mathit{\pi }/\mathrm{2})\approx \mathrm{0.944}$, the reduction for Q₂ coincides with a rotation by three quarters of a rotation ($t=\ln (\mathrm{1}+\mathrm{3}\mathit{\pi }/\mathrm{2})\approx \mathrm{1.742}$), and the reduction for RK2FOT coincides with a full rotation ($t=\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.