Analytical solutions for mantle flow in cylindrical and spherical shells
- 1Department of Earth Science and Engineering, Imperial College London, London, UK
- 2Research School of Earth Sciences, The Australian National University, Canberra, Australia
- 3Earth and Planets Laboratory, Carnegie Institution for Science, Washington, DC, USA
Correspondence: Stephan C. Kramer (email@example.com)
Computational models of mantle convection must accurately represent curved boundaries and the associated boundary conditions of a 3-D spherical shell, bounded by Earth's surface and the core–mantle boundary. This is also true for comparable models in a simplified 2-D cylindrical geometry. It is of fundamental importance that the codes underlying these models are carefully verified prior to their application in a geodynamical context, for which comparisons against analytical solutions are an indispensable tool. However, analytical solutions for the Stokes equations in these geometries, based upon simple source terms that adhere to physically realistic boundary conditions, are often complex and difficult to derive. In this paper, we present the analytical solutions for a smooth polynomial source and a delta-function forcing, in combination with free-slip and zero-slip boundary conditions, for both 2-D cylindrical- and 3-D spherical-shell domains. We study the convergence of the Taylor–Hood (P2–P1) discretisation with respect to these solutions, within the finite element computational modelling framework Fluidity, and discuss an issue of suboptimal convergence in the presence of discontinuities. To facilitate the verification of numerical codes across the wider community, we provide a Python package, Assess, that evaluates the analytical solutions at arbitrary points of the domain.
Mantle convection transports Earth's internal heat to its surface: it is the “engine” driving our dynamic Earth (e.g. Davies, 1999). The structure, composition, and flow regime within the mantle are reflected in near-surface phenomena such as plate tectonics, mountain building, dynamic topography, sea-level change, volcanism, and the activity of Earth's magnetic field (e.g. Morgan, 1972; Mitrovica et al., 1989; Gurnis, 1993; Olson et al., 2013; Kloecking et al., 2019; Davies et al., 2019). The grand challenge is to understand the operation of this giant heat engine over geologic time and its relationship to the surface geological record.
Computational modelling is one of the primary tools available for tackling this challenge. Whilst 2- and 3-D numerical models of mantle convection processes in Cartesian domains have provided important insights into a range of mantle processes (e.g. McKenzie et al., 1974; Gurnis and Davies, 1986; Davies and Stevenson, 1992; Moresi and Solomatov, 1995; van Keken et al., 2002; Hunt et al., 2012; Garel et al., 2014; Davies et al., 2016; Jones et al., 2016, 2019), 3-D spherical geometry is required to simulate global mantle dynamics. Global 3-D spherical mantle convection models, and studies focussing on their application, are now in common use (e.g. Baumgardner, 1985; Tackley et al., 1993; Bunge et al., 1996, 1997; Zhong et al., 2000; Oldham and Davies, 2004; McNamara and Zhong, 2005; Choblet et al., 2007; Zhong et al., 2008; Tackley, 2008; Davies and Davies, 2009; Wolstencroft et al., 2009; Stadler et al., 2010; Tan et al., 2011; Kronbichler et al., 2012; Davies et al., 2013; Burstedde et al., 2013; Heister et al., 2017; Dannberg and Gassmoller, 2018; Stotz et al., 2018). However, the use of this geometry for calculations at a realistic convective vigour remains expensive. As a consequence, simplifying geometries are often used, including the axisymmetric spherical shell (e.g. Solheim and Peltier, 1994; van Keken and Yuen, 1995), the 2-D cylinder (e.g. Jarvis, 1993; van Keken and Ballentine, 1998, 1999; van Keken, 2001; Nakagawa and Tackley, 2005), and the spherical annulus (e.g. Hernlund and Tackley, 2008). Such 2-D models can also provide a rapid and broad parameter space appraisal, allowing one to focus on a targeted and more sparse study in a full 3-D spherical geometry.
Recent decades have seen extensive validation, verification, and benchmarking of Cartesian mantle dynamics codes, in both 2- and 3-D. Verification is typically achieved via comparisons of numerical predictions against analytical solutions (e.g. Zhong et al., 1993; Kramer et al., 2012), with benchmarking typically undertaken against solutions from other comparable codes, for incompressible (Blankenbach et al., 1989; Travis et al., 1990; Busse et al., 1994; van Keken et al., 1997, 2008; Tosi et al., 2015) and compressible (King et al., 2009) convection. See also Popov et al. (2014) for an overview of geodynamical benchmarks. The number of comparable studies, within a 2-D cylindrical or 3-D spherical geometry, however, is more limited. Given a recent surge in the state-of-the-art tools available for simulating mantle dynamics in these geometries (e.g. Kronbichler et al., 2012; Logg et al., 2012; Burstedde et al., 2013; Rathgeber et al., 2016; Heister et al., 2017; Wilson et al., 2017), it is important that these tools are verified and validated against a range of analytical and benchmark solutions.
A popular method to obtain analytical solutions is the method of manufactured solutions (MMS; Roache, 2002). With this approach, an arbitrary analytical solution is chosen beforehand and the necessary forcing terms on the right-hand side are derived by substitution of the solution into the left-hand side of the flow equations. A drawback of this approach is that it can be hard to choose solutions (i) that satisfy physically realistic boundary conditions, in particular, in less trivial domains; and (ii) with a velocity that is divergence free, in order to avoid unnatural source terms. Alternatively, one can choose a simple analytical expression for the forcing but the derivation of the corresponding analytical solutions is often laborious.
For the Stokes equations, a well-known set of solutions in the latter category is based on a forcing term in the form of a delta function, corresponding to an infinitely thin density anomaly at a certain depth. The solutions to these are also used in the propagator matrix method, where a convolution of delta-function solutions at different depths is used to obtain the response to arbitrary density anomalies (e.g. Hager and Richards, 1989). Solutions for the Cartesian case have been published (e.g. Zhong et al., 1993; Kramer et al., 2012). Spherical solutions can be found in, for example, Ribe (2009). They have previously been used for validation by, for example, Zhong et al. (2000, 2008), Choblet et al. (2007), Davies et al. (2013), and Burstedde et al. (2013). The derivation of these solutions is non-trivial, and often only cases with simpler free-slip boundary conditions are explicitly presented. Here, we derive solutions for both free- and zero-slip boundary conditions. In our convergence analyses, where we compare the analytical solutions with numerical solutions obtained using the Fluidity computational modelling framework (Davies et al., 2011), we discuss a limitation of this set of solutions for the benchmarking of geodynamic codes.
For this reason, we also present solutions based on a smooth forcing term, with radial dependence formed by a polynomial of arbitrary order, again for free- and zero-slip cases. This set of solutions provides a flexible way to test mantle dynamics codes with physically relevant solutions, where, for instance, a high-order polynomial forcing can be used to obtain solutions with a strong gradient near the surface. The radial dependence can be combined with spherical harmonics of arbitrary degree and order.
A key step in the derivation of these benchmarks is a decomposition of the solution into poloidal and toroidal components in the Mie representation (Backus, 1986). This results in a biharmonic equation for the poloidal scalar function. Combined with a set of conditions for the radial dependence of this poloidal function, analytical Stokes solutions can be obtained that satisfy desired free-slip or zero-slip conditions. Similar techniques have been used in Tosi and Martinec (2007) and Horbach et al. (2020). In Sect. 5.1, we will discuss previously published analytical benchmark cases in shell domains and how they differ from those presented here (e.g. Blinova et al., 2016; Thieulot, 2017; Horbach et al., 2020).
Finally, in addition to the delta-function and smooth cases with either free-slip or zero-slip boundary conditions in a spherical-shell domain, we present the solutions for the corresponding four cases in a 2-D cylindrical-shell domain (annulus). Although, ultimately, mantle convection is a 3-D phenomenon, a number of processes can be modelled adequately in two dimensions, and accordingly access to benchmark cases for 2-D numerical models is equally important. The number of published analytical Stokes solutions in 2-D cylindrical-shell domains, which are suitable as geodynamical benchmarks and include a complete derivation, is limited. By presenting this extensive set of explicit analytical solutions, we provide a suite of verification cases for use by the wider community of mantle dynamics code developers. An implementation of the solutions is provided through the Python package Assess (Analytical Solutions for the Stokes Equations in Spherical Shells; Kramer, 2020).
The remainder of this paper is structured as follows. In Sect. 2, we derive the analytical solutions in cylindrical (Sect. 2.2) and spherical (Sect. 2.3) geometries. Smooth solutions are first provided, for free-slip and zero-slip cases, followed by the delta-function solutions. In Sect. 3, we briefly describe the Fluidity computational modelling framework, which is used in Sect. 4 to obtain convergence results of the P2–P1 finite element discretisation to the analytical solutions. In Sect. 5, we discuss these results and relate the analytical solutions presented here to those that have previously been published. In particular, we discuss an issue with the delta-function cases for discretisations that use continuous pressure. To demonstrate, we also show results with a P2bubble–P1DG discretisation that overcomes this issue.
The following derivation is applicable to the incompressible Stokes equations:
where the unknowns are velocity u and pressure p, with an assumed constant kinematic viscosity ν. The buoyancy force on the right-hand side (RHS) of Eq. (1) is based on a gravity of constant magnitude g, directed in the inward radial direction , and a dimensionless density deviation ρ′: , where ρ is the density and the reference density ρ0 is constant. The equations are solved in a 2-D cylindrical or 3-D spherical domain, bounded by , where r is the radial distance to the origin.
In 2-D, any incompressible velocity field u can be written as the skew gradient of a streamfunction ψ:
where φ is the angle from the x axis in polar coordinates, , and ur and uφ are the radial and transverse components of velocity, respectively. The normal deviatoric stress and shear stress components are given by
The momentum equation (Eq. 1) can be decomposed into radial and transverse components:
where the Laplacian, in polar coordinates, is given by
For a derivation of Eqs. (5)–(8), see Appendix A.
The general, real-valued solution to the biharmonic equation ∇4ψ=0 is given by
where An, Bn, Cn, Dn, en, and fn are constant coefficients, n≥0. The An and Bn terms are the standard solutions to the homogeneous harmonic equation, and the Cn and Dn terms are obtained as inhomogeneous solutions of the harmonic equation with the homogeneous harmonic solutions (An and Bn terms) as the right-hand side. The fact that these (the Cn and Dn terms) are homogeneous biharmonic solutions then follows from ∇4=∇2∇2. In the following, we will, for simplicity, focus on sin (nφ) solutions for a single n>1 and set .
An equation for pressure can be derived by taking the divergence of the momentum equation:
From Eq. (7), it can be seen that sin (nφ) solutions for ψ are associated with cos (nφ) solutions of p and ρ′. The homogeneous solutions, i.e. , are thus the standard harmonic solutions (again neglecting the solutions)
where Gn and Hn are constant coefficients, n>1.
between the coefficients of the homogeneous solutions for ψ and p can be derived.
2.2.1 Smooth density profile – cylindrical
In the first test case, we consider a density perturbation of the following form:
with k>0. It is easily verified that an inhomogeneous solution to Eq. (10) exists of the form
assuming and , so that a general solution can be written as
Note that for the remainder of this derivation we drop the subscript n in the coefficients for , and D.
An inhomogeneous solution for pressure of Eq. (12) is given by
The general solution for pressure is thus
with G and H given by Eq. (14).
The four remaining coefficients ( and D) are fixed by a choice of boundary conditions at the inner and outer surfaces of the cylindrical domain at and , respectively. At both, no-normal-flow conditions are imposed via . Two further equations are found by imposing either τrφ=0 (free slip) or (zero slip) at both boundaries.
The solution coefficients for free-slip, no normal flow at both boundaries are given by
where we use .
The zero-slip solution coefficients are given by
2.2.2 Green's function solution – cylindrical
Another set of useful solutions is found considering the following perturbation density,
representing an infinitely thin density anomaly at with . Since for , the solution is described by combining two homogeneous solutions
We thus have eight coefficients that are fixed by boundary conditions and conditions at the interface. Boundary conditions at the inner and outer boundaries again provide four equations. Furthermore, we impose continuity of both velocity components at the interface:
Since no lateral force is being applied at the interface, we also expect continuity of the shear stress (6) which, in combination with the above, implies
Finally, the eighth equation is obtained by integrating Eq. (10) over a small strip , between two arbitrary angles :
Taking the limit ϵ→0, the second term on the left-hand side disappears, whereas the first term becomes
Thus, for this to hold for arbitrary φ1 and φ2, we need
The solution coefficients for free-slip boundary conditions at and are
where we use .
Zero-slip conditions at both boundaries lead to
where, in addition, we use .
In this section, we derive the equivalent of the four cylindrical cases in a 3-D spherical domain, . A more detailed derivation of the equations can be found in Ribe (2009). In 3-D, a solenoidal velocity field u can be decomposed as
where φ is the longitude and θ is the co-latitude, and
We further derive the components of stress normal to a spherical surface:
Using these, we can work out spherical components of the momentum equation (Eq. 1) to be
As can be seen from these equations, the toroidal part of the solution is independent of the density distribution. A non-zero toroidal component is only introduced through a particular choice of boundary conditions. For the test cases in this paper, we only consider free-slip and zero-slip boundary conditions with radial forcing so we may assume 𝒯=0. The assumption of the velocity field being solenoidal, which underlies the Mie representation Eq. (27), follows directly from the incompressibility condition and the no-normal-flow condition on the boundary.
Taking the curl of the momentum equation,
In a spherical domain, this implies that varies in the radial direction only. Without loss of generality, we may therefore seek solutions 𝒫 that satisfy
since for any other solution of Eq. (36) we know that ∇4𝒫′ is purely radial, and therefore 𝒫′ can be written as a sum of biharmonic solutions and a purely radial function. The biharmonic solutions are included in Eq. (37) and the purely radial function is discarded by Eq. (27).
Solutions to the biharmonic equation ∇4𝒫=0 in 3-D are given by
where Alm, Blm, Clm, and Dlm are constant coefficients, and Ylm indicates Laplace's spherical harmonic functions of degree l and order m, which have the property that
As in 2-D, an equation for the pressure is obtained by taking the divergence of the momentum equation, giving
Homogeneous solutions for pressure are written as
2.3.1 Smooth density profile – spherical
We consider a density perturbation of the following form:
An inhomogeneous solution of Eq. (37) is given by
with a more generic solution written as
where we have dropped the lm subscripts of coefficients A, B, C, and D. The pressure solution for Eq. (40) is
and G and H are given by Eq. (42).
As before, the four coefficients (A, B, C, and D) are fixed by the boundary conditions:
where we use . The two no-normal-flow conditions at and are combined with two further conditions:
For free-slip conditions at both and , the solution coefficients are given by
Zero-slip conditions at both boundaries lead to the following solution coefficients:
2.3.2 Green's function solution – spherical
As in two dimensions, we find solutions for the case where
by combining two homogeneous solutions:
The eight coefficients are found by imposing the same four constraints derived from the boundary conditions at and as in the previous section and by imposing a further four conditions: continuity of all components of u, no shear force between the two halves of the domain, and a normal force that is proportional to the density anomaly:
The free-slip solution coefficients are given by
The zero-slip solution coefficients are given by
The numerical solutions for velocity and pressure, u and p, are written as a linear combination of basis functions Nj and Ml:
We use either the P2–P1 (Taylor–Hood) or P2bubble–P1DG element pairs. In both cases, the velocity and pressure basis functions Nj and Ml are piecewise quadratic and linear, respectively, on a triangular (2-D) or tetrahedral (3-D) mesh of the domain Ω. Because the curved boundaries of the cylindrical- and spherical-shell domains can only be approximated by the mesh, the numerical domain is denoted by Ωh. When using the P2–P1 element pair, the basis functions are continuous between cells. For the P2bubble–P1DG element pair, the piecewise-linear pressure is treated as discontinuous between cells and the continuous quadratic velocity basis functions are enriched by an extra cubic “bubble” function with a corresponding cell-centred degree of freedom (Ern and Guermond, 2004; Boffi et al., 2013).
Note that we apply strong Dirichlet boundary conditions, so that the boundary integrals can indeed be neglected. In the free-slip case, a local rotation is applied to the velocity vectors, so that the degrees of freedom correspond to velocity components in either the normal or tangential directions. This allows us to enforce a zero-normal component, while leaving the tangential components free. For additional details about Fluidity and its implementation, see Davies et al. (2011).
3.2 Isoparametric representation of the domain
For an accurate representation of the quadratic approximation of velocity, we need to also approximate the curved cylindrical/spherical domain quadratically. This means that rather than each cell in the mesh being described by a linear map from local coordinates ξ in a reference element to physical coordinates X, we use a quadratic map Xquad(ξ), which maps to a curved triangle/tetrahedron that better represents the domain. This map can be obtained from a linear mesh with coordinate mappings, Xlin, through quadratic interpolation:
at the standard Lagrange nodes of the quadratic function space. Here, rlin(ξ) is the linear interpolation of the radius, i.e. at the vertices of the linear mesh. This particular choice ensures that for an equal-radius boundary, with the boundary vertices of the linear mesh exactly on the boundary, the quadratic Lagrange nodes also lie exactly on this boundary.
3.3 Forcing term
The density perturbation ρ′ on the RHS of Eq. (57) is a prescribed analytical expression in each of the test cases. For Green's function solutions in 2-D and 3-D, using Eqs. (20) and (50), respectively, we get
where Γ′ is an internal boundary at oriented such that its normal points outwards.
3.4 Solving the linear system and dealing with null spaces
Equations (57) and (58) form a saddle
point linear system which is solved by applying a Schur decomposition technique
where the outer iteration, which solves for the pressure degrees of freedom, is
solved with a flexible Krylov subspace method,
FGMRES. The inner solve associated with
the velocity degrees of freedom, is solved with the conjugate gradient method
preconditioned with an algebraic multigrid method
GAMG available through PETSc; Balay et al., 1997).
In all cases, the pressure solution is only defined up to an arbitrary constant. The analytical pressure solution has the property that its mean is zero. For comparative purposes, we therefore subtract the volume-averaged pressure from the obtained numerical pressure solution:
Similarly, for free-slip cases, in 2-D, we may add an arbitrary rotation of the form to the velocity solution. We therefore apply the following projection to the numerical solution:
which ensures that the angular momentum is zero, as it is for the analytical solutions. In the same way, in three dimensions, we subtract the three rotational (rigid body) modes.
It should be noted that the same velocity and pressure modes lead to zero modes (eigenvectors) for the linear system based on Eqs. (57) and (58), rendering the resulting matrix singular. In preconditioned Krylov methods, we typically need to subtract the zero modes from the approximate solution at every iteration. With iterative approximation xi and zero eigenvector λ, we get
This l2 projection, based on the l2 inner product , is analogous but not equivalent to the projections in Eqs. (62) and (63) (which are L2 projections). Therefore, despite the l2 projections fixing the null modes during the iterative solve, the L2 projections should be applied as an additional step after the iterative solvers have completed to ensure convergence to the analytical solution.
In this section, we show the convergence of the numerical solutions obtained with Fluidity, using the P2–P1 element pair, towards the analytical solutions. For 2-D cylindrical cases, the series of meshes start at refinement level 1, where the mesh consists of 128 divisions in the horizontal, and 16 layers, giving triangles. At each subsequent level, the mesh is refined, doubling the resolution in both directions. For the spherical cases, the mesh at refinement level 1 is obtained from an icosahedron refined three times, starting with 1280 triangles in the horizontal, which is extruded radially to 16 layers, giving a 3-D mesh consisting of 61 440 tetrahedra. Again, resolution is doubled in all directions for subsequent refinement levels. In all cases, non-dimensionalised coordinates were used with and , and the delta-function cases used . This choice of r′ ensures the density anomaly coincides with a grid layer at all mesh resolutions considered herein.
In all figures, errors are given as relative errors, comparing the numerical solution, u and p, with the analytical solutions, u∗ and p∗ (interpolated into the P2 and P1 function spaces, respectively) in the L2 norm:
Convergence plots for the 2-D cylindrical cases are presented in Fig. 3. With a smooth density profile, we see optimal convergence for the P2–P1 element pair at third and second order for velocity and pressure, respectively, with both free-slip and zero-slip boundary conditions. Cases with lower wavenumber n show smaller relative error than those at higher n, as expected. The same observation holds for lower and higher polynomial order, k=2 and k=8, for the radial density profile. For the free-slip and zero-slip delta-function cases, however, convergence drops to 1.5 and 0.5 for velocity and pressure, respectively. Furthermore, cases with lower n do not consistently show smaller relative error than those at higher n.
We see similar results for the spherical results illustrated in Fig. 4: third and second order for velocity and pressure for the cases with a smooth density profile, with smaller relative errors for lower wave numbers l and m. Note that here, the smooth vertical profile for density uses in all cases. Again, for cases with a delta-function density anomaly, we observe a reduced order of convergence of 1.5 and 0.5 for velocity and pressure, respectively.
To examine the importance of an isoparametric approximation of the domain by a quadratic mesh, we ran the same cases with a linear mesh. The results are shown in Fig. 5, which demonstrates that the order of convergence of velocity in the smooth cylindrical cases is indeed limited to second order. The convergence of pressure remains at second order.
5.1 Existing analytical benchmarks in shell domains
As indicated in the introduction, spherical delta-function cases, like those presented herein, have previously been used to validate global mantle convection codes (e.g. Zhong et al., 2008; Burstedde et al., 2013; Davies et al., 2013; Liu and King, 2019). These will be discussed in more detail in the following section.
The derivation for all 3-D cases in this paper relies on the Mie representation that decomposes the velocity solution into poloidal and toroidal components, through which, under the assumption of purely poloidal flow, the Stokes equations can be reduced to a biharmonic equation (10). Any solution to the inhomogeneous solution can then be combined with four linearly independent homogeneous solutions to this equation, the coefficients of which can be derived through the imposition of boundary conditions. For the smooth case with a generic monomial forcing term, the same decomposition (i.e. Eqs. 45 and 44) was used in Tosi and Martinec (2007) to derive the analytical solution for Stokes flow in two eccentrically nested spheres.
In Horbach et al. (2020), similar techniques were employed to derive benchmarks in spherical-shell domains that satisfy zero-slip and free-slip conditions and, in addition, a mixed zero-slip/free-slip case. Here, the derivation starts by simply selecting four, in principle arbitrary, linearly independent solutions for the radially dependent part of the poloidal scalar function. Again the imposition of boundary conditions fixes the coefficients of this linear combination. The corresponding right-hand side forcing term is then obtained by substitution.
The number of published benchmarks for 2-D cylindrical-shell domains is more limited (e.g. Buffa et al., 2011; Blinova et al., 2016; Hoang et al., 2017). The derivation of the equivalent cases in 2-D cylindrical-shell domains is somewhat simpler, but also ultimately relies on combining four independent homogeneous solutions and one inhomogeneous solution to a biharmonic equation.
In Blinova et al. (2016), analytical solutions in both cylindrical and spherical domains are presented for the Stokes equations with a radially dependent viscosity. Because these are solutions in cylindrical or spherical coordinates without reference to any specific domain, they do not satisfy no-normal-flow conditions on the boundary of a shell domain. They can be used in such domains as a numerical benchmark by specifying all velocity components of the analytical solution as a Dirichlet condition for the model. Analytical solutions for radially dependent viscosity were also presented in Thieulot (2017). Their solutions (in 3-D only) do satisfy no-normal-flow conditions in a spherical-shell domain, but the tangential components are non-zero at the boundary and thus still require inhomogeneous Dirichlet conditions. Spatially varying viscosity is of course an import aspect of mantle convection models for which these are effective benchmarks. The isoviscous solutions presented here, and those in Horbach et al. (2020), however, allow for the testing of zero-slip and free-slip conditions, where in particular free-slip conditions may pose various numerical challenges such as rotational modes and, depending on the discretisation used, the (non-)alignment of velocity components with normal and tangential directions at the boundary.
5.2 Reduced order of convergence with discontinuous pressures
At first sight, the reduced order of convergence for the delta-function cases seems at odds with those expected for the P2–P1 element pair. However, the mathematical proofs for the ideal order of convergence to solutions of the Stokes equations rely on certain regularity assumptions of the right-hand side forcing term and, related to that, on the regularity of the velocity and pressure solutions. The regularity of the delta function can be classified as being a member of the Sobolev space H−1(Ω) the dual of the Sobolev space H1(Ω), where for the sake of simplicity we assume Ω=Ωh in this section. This means that the delta function can be thought of as a continuous function:
which maps functions v∈H1(Ω), the space of square integrable functions with square integrable weak derivatives, to ℝ. Girault and Raviart (2012) demonstrate that even with the very loose regularity condition that the right-hand side f is in Hm(Ω) with , the Stokes equations in the weak form have a unique solution (given sufficient integral constraints) with velocity in Hm+2(Ω) and pressure in Hm+1(Ω). The analytical solutions derived here for the delta case with indeed have a discontinuous pressure in H0(Ω)=L2(Ω) and a velocity with discontinuous normal derivative in H1(Ω).
For velocity–pressure finite element pairs that satisfy the standard inf-sup, or Ladyzhenskaya–Babuška–Brezzi (LBB) condition
where V and W are the discrete vector and scalar function spaces, respectively, and β is a constant, it can be shown that the method converges and in fact
where C1 is a constant independent of h, u∗ and p∗ are the exact solutions, u and p the numerical solutions in the discrete function spaces V and W based on a mesh with mesh distance h, and is the semi-norm
In other words, the convergence of u and p to u∗ and p∗ is bounded by the best possible approximation of u∗ and p∗ in the discrete spaces V and W. For bounded functions with a discontinuity along a smooth interface, such as our analytical pressure solution p∗, the best approximation by continuous, piecewise-linear polynomials, i.e. W=P1, is bounded by
A solution to this problem is found by allowing for discontinuities in the discrete pressure space. We demonstrate this here by considering W=P1DG the space of piecewise-linear but discontinuous functions. To satisfy the inf-sup condition, this requires enriching the quadratic function space P2 for velocity with a cubic bubble (Ern and Guermond, 2004; Boffi et al., 2013). Convergence results for this element pair are shown in Fig. 6. For the 2-D cylindrical delta-function cases, we observe the expected orders of convergence: third order for velocity and second order for pressure.
Finally, we compare our results with those presented in Zhong et al. (2008), Davies et al. (2013), Burstedde et al. (2013), and Liu and King (2019), who ran the same spherical cases with a delta-function forcing and found second-order convergence for velocity and second-order convergence for pressure related diagnostics. It should be noted, however, that these studies only examined surface diagnostics, such as the velocity divergence and the normal stress. When we examine comparable diagnostics, specifically, the relative error in velocity at the boundary and the boundary normal stress (illustrated in Fig. 7), we find velocity convergence at third order and normal stress at second order. It therefore appears that the reduced order of accuracy in the interior of the domain does not affect the surface response which still converges at the same order as for the smooth case.
The results of Zhong et al. (2008) (using CitcomS) and Davies et al. (2013) (using TERRA) were based on a Q1–P0 discretisation with a continuous piecewise-trilinear velocity and piecewise-constant, discontinuous pressure. Although our analysis above limits the convergence of for a P0 pressure p to first order, second-order super-convergence can be obtained in some cases by evaluating the analytical solution only in the cell centre. In other words, by comparing to a filtered piecewise-constant analytical approximation , second-order convergence can sometimes be observed in . For continuous pressure approximations, such as the P2–P1 results in this paper, the Q1–Q1 discretisation of the Rhea model in Burstedde et al. (2013), and the Q2–Q1 discretisation of ASPECT in Liu and King (2019), reduced convergence in the interior of the domain is to be expected.
We have presented a series of 2-D cylindrical and 3-D spherical analytical solutions for the purpose of verifying mantle dynamics codes. These solutions are based upon either a delta-function density perturbation or a smooth forcing term, and we provide solutions for both free-slip and zero-slip boundary conditions. The combinations of dimension, forcing, and boundary conditions provide a series of eight analytical solutions that can be used as a basis for validating existing and future numerical codes, in cylindrical and spherical geometries. To facilitate this, we provide solutions in the form of a Python package (Assess: Analytical Solutions for the Stokes Equations in Spherical Shells; Kramer, 2020).
We verify the convergence of the P2–P1 (Taylor–Hood) finite element discretisation using Fluidity (Davies et al., 2011; Kramer et al., 2012). The continuous approximation of pressure can lead to a reduced order of convergence in the presence of discontinuities, which can be overcome using a discontinuous numerical approximation of pressure. It is important to note that this reduced order of convergence was only observed by comparing the numerical solution with the entire analytical solution in the interior of the domain. A comparison based on surface response only failed to highlight this issue.
In this appendix, we work out the incompressible Stokes equations in polar coordinates in terms of a streamfunction ψ, where the components of velocity are given by
We make use of the following expressions for the derivatives of the unit vectors and with respect to r and φ:
Using these, we can work out the different components of stress:
Note that, as expected, . In the same way, we derive the following expression for the vorticity, or curl of the velocity:
The viscosity term in the Stokes equations can be written as
In addition to Eq. (A3), we use the following identities:
and the fact that . After reordering to group the radial and transverse components, this leads to
In combination with the pressure gradient term in polar coordinates
The Python package Assess, which implements the analytical solutions and evaluates them at arbitrary locations in the domain, is available from http://github.com/stephankramer/assess (last access: 12 June 2020) (see https://assess.readthedocs.io (last access: 12 June 2020). An archived version is available from https://doi.org/10.5281/zenodo.3891545 (Kramer, 2020). To ensure correctness, in the paper as well as in the Python package, the coefficients for the various solutions are extracted from the LaTeX source automatically using SymPy (Meurer et al., 2017), verified to adhere to the equations using SageMath (The Sage Developers, 2019), and substituted in the Python package. Assess has also been used to compute the errors in the Fluidity results in this paper. The Fluidity model, including source code and documentation, is available from https://fluidityproject.github.io/ (last access: 17 August 2020). All cases in this paper have been run with tag version 4.1.17, which is archived at https://doi.org/10.5281/zenodo.3988620 (Kramer et al., 2020).
SCK derived the analytical solutions presented herein and updated Fluidity to run these cases. DRD ran the cases and prepared the corresponding figures. SCK wrote the manuscript with input from DRD and CRW. CRW implemented, set up, and ran the P2bubble–P1DG cases. All authors had significant input on the design and development of this research.
The authors declare that they have no conflict of interest.
D. Rhodri Davies and Stephan C. Kramer acknowledge support from the Australian Research Council, under grant nos. FT140101262 and DP170100058. Stephan C. Kramer also acknowledges support from the UK Engineering and Physical Sciences Research Council (EPSRC) under grant no. EP/R029423/1. Numerical simulations were undertaken on the NCI National Facility in Canberra, Australia, which is supported by the Australian Commonwealth Government. We would also to thank Marcus Mohr for his very thorough review of the paper and thank him and a second anonymous reviewer for their constructive feedback.
This research has been supported by the Australian Research Council (grant nos. DP170100058 and FT140101262) and the UK Engineering and Physical Sciences Research Council (EPSRC) (grant no. EP/R029423/1).
This paper was edited by Thomas Poulet and reviewed by Marcus Mohr and one anonymous referee.
Balay, S., Gropp, W. D., McInnes, L. C., and Smith, B. F.: Efficient management of parallelism in object-oriented numerical software libraries, in: Modern software tools for scientific computing, Birkhauser, Boston Inc., 163–202, 1997. a
Bernardi, C.: Optimal finite-element interpolation on curved domains, SIAM J. Numer. Anal., 26, 1212–1240, 1989. a
Blankenbach, B., Busse, F., Christensen, U., Cserepes, L., Gunkel, D., Hansen, U., Harder, H., Jarvis, G., Koch, M., Marquart, G., Moore, D., Olson, P., Schmeling, H., and Schnaubelt, T.: A benchmark comparison for mantle convection codes, Geophys. J. Int., 98, 23–38, https://doi.org/10.1111/j.1365-246X.1989.tb05511.x, 1989. a
Blinova, I., Makeev, I., and Popov, I.: Benchmark solutions for stokes flows in cylindrical and spherical geometry, Bulletin of the Transilvania University of Brasov., Mathematics, Informatics, Physics. Series III, 9, 11–16, 2016. a, b, c
Boffi, D., Brezzi, F., and Fortin, M.: Incompressible Materials and Flow Problems, in: Mixed Finite Element Methods and Applications, Springer Series in Computational Mathematics, Springer, Berlin, Heidelberg, chap. 8, 459–538, 2013. a, b
Buffa, A., De Falco, C., and Sangalli, G.: Isogeometric analysis: stable elements for the 2D Stokes equation, Int. J. Numer. Meth. Fl., 65, 1407–1422, 2011. a
Bunge, H.-P., Richards, M. A., and Baumgardner, J. R.: A sensitivity study of 3-D-spherical mantle convection at 108 Rayleigh number: effects of depth-dependent viscosity, heating mode and an endothermic phase change, J. Geophys. Res., 102, 11991–12007, https://doi.org/10.1029/96JB03806, 1997. a
Burstedde, C., Stadler, G., Alisic, L., Wilcox, L. C., Tan, E., Gurnis, M., and Ghattas, O.: Large-scale adaptive mantle convection simulation, Geophys. J. Int., 192, 889–906, https://doi.org/10.1093/gji/ggs070, 2013. a, b, c, d, e, f
Busse, F. H., Christensen, U., Clever, R., Cserepes, L., Gable, C., Giannandrea, E., Guillou, L., Houseman, G., Nataf, H. C., Ogawa, M., Parmentier, M., Sotin, C., and Travis, B.: 3D convection at infinite Prandtl number in Cartesian geometry – a benchmark comparison, Geophys. Astrophys. Fluid Dyn., 75, 39–59, https://doi.org/10.1080/03091929408203646, 1994. a
Choblet, G., Cadek, O., Couturier, F., and Dumoulin, C.: OEDIPUS: a new tool to study the dynamics of planetary interiors, Geophys. J. Int., 170, 9–30, https://doi.org/10.1111/j.1365-246X.2007.03419.x, 2007. a, b
Davies, D. R., Wilson, C. R., and Kramer, S. C.: Fluidity: a fully unstructured anisotropic adaptive mesh computational modeling framework for geodynamics, Geochem. Geophy. Geosy., 120, Q06001, https://doi.org/10.1029/2011GC003551, 2011. a, b, c, d
Davies, D. R., Davies, J. H., Bollada, P. C., Hassan, O., Morgan, K., and Nithiarasu, P.: A hierarchical mesh refinement technique for global 3-D spherical mantle convection modelling, Geosci. Model Dev., 6, 1095–1107, https://doi.org/10.5194/gmd-6-1095-2013, 2013. a, b, c, d, e
Davies, D. R., Le Voci, G., Goes, S., Kramer, S. C., and Wilson, C. R.: The mantle wedge's transient 3-D flow regime and thermal structure, Geochem. Geophy. Geosy., 17, 78–100, https://doi.org/10.1002/2015GC006125, 2016. a
Davies, D. R., Valentine, A. P., Kramer, S. C., Rawlinson, N., Hoggard, M. J., Eakin, C., and Wilson, C.: Earth's multi-scale topographic response to global mantle flow, Nat. Geosci., 12, 845–850, https://doi.org/10.1038/s41561-019-0441-4, 2019. a
Davies, G. F.: Dynamic Earth: plates, plumes and mantle convection, Cambridge University Press, Cambridge, 1999. a
Garel, F., Goes, S., Davies, D. R., Davies, J. H., Kramer, S. C., and Wilson, C. R.: Interaction of subducted slabs with the mantle transition-zone: A regime diagram from 2-D thermo-mechanical models with a mobile trench and an overriding plate, Geochem. Geophy. Geosy., 15, 1739–1765, https://doi.org/10.1002/2014GC005257, 2014. a
Girault, V. and Raviart, P.-A.: Finite element methods for Navier-Stokes equations: theory and algorithms, vol. 5, Springer Science & Business Media, Berlin, Heidelberg, 2012. a
Hager, B. H. and Richards, A. M.: Long–wavelength variations in Earth's geoid: Physical models and dynamical implications, Philos. T. R. Soc. London, 328, 309–327, 1989. a
Heister, T., Dannberg, J., Gassmoller, R., and Bangerth, W.: High accuracy mantle convection simulation through modern numerical methods: II: realistic models and problems, Geophys. J. Int., 210, 833–851, https://doi.org/10.1093/gji/ggx195, 2017. a, b
Hoang, T., Verhoosel, C. V., Auricchio, F., van Brummelen, E. H., and Reali, A.: Mixed isogeometric finite cell methods for the stokes problem, Comput. Meth. Appl. M., 316, 400–423, 2017. a
Hunt, S. A., Davies, D. R., Walker, A. M., McCormack, R. J., Wills, A. S., Dobson, D. P., and Li, L.: On the increase in thermal diffusivity caused by the perovskite to post-perovskite phase transition and its implications for mantle dynamics, Earth Planet. Sc. Lett., 319, 96–103, https://doi.org/10.1016/j.epsl.2011.12.009, 2012. a
Jarvis, G. T.: Effects of curvature on two–dimensional models of mantle convection: cylindrical polar coordinates, J. Geophys. Res., 98, 4477–4485, 1993. a
Jones, T. D., Davies, D. R., Campbell, I. H., Wilson, C. R., and Kramer, S. C.: Do mantle plumes preserve the heterogeneous structure of their deep-mantle source?, Earth Planet. Sc. Lett., 434, 10–17, https://doi.org/10.1016/j.epsl.2015.11.016, 2016. a
Jones, T. D., Davies, D. R., and Sossi, P. A.: Tungsten isotopes in mantle plumes: Heads it's positive, tails it's negative, Earth Planet. Sc. Lett., 506, 255–267, https://doi.org/10.1016/j.epsl.2018.11.008, 2019. a
King, S. D., Lee, C., van Keken, P. E., Leng, W., Zhong, S., Tan, E., Tosi, N., and Kameyama, M. C.: A community benchmark for 2-D Cartesian compressible convection in Earth's mantle, Geophys. J. Int., 179, 1–11, 2009. a
Kloecking, M., White, N. J., Maclennan, J., McKenzie, D., and Fitton, J. G.: Quantitative Relationships Between Basalt Geochemistry, Shear Wave Velocity, and Asthenospheric Temperature Beneath Western North America, Geochem. Geophy. Geosy., 19, 3376–3404, https://doi.org/10.1029/2018GC007559, 2019. a
Kramer, S. C., Wilson, C. R., and Davies, D. R.: An implicit free-surface algorithm for geodynamical simulations, Phys. Earth Planet. Int., 194, 25–37, https://doi.org/10.1016/j.pepi.2012.01.001, 2012. a, b, c, d
Kramer, S. C., Wilson, C. R., Davies, D. R., et al.: FluidityProject/fluidity: New test cases Analytical solutions for mantle flow in cylindrical and spherical shells (Version 4.1.17), Zenodo, https://doi.org/10.5281/zenodo.3988620, 2020. a
Kronbichler, M., Heister, T., and Bangerth, W.: High accuracy mantle convection simulation through modern numerical methods, Geophys. J. Int., 191, 12–29, https://doi.org/10.1111/j.1365-246X.2012.05609.x, 2012. a, b
Logg, A., Mardal, K.-A., and Wells, G.: Automated Solution of Differential Equations by the Finite Element Method: The FEniCS Book, Lecture Notes in Computational Science and Engineering, Springer, Berlin, 84, https://doi.org/10.1007/978-3-642-23099-8, 2012. a
Meurer, A., Smith, C. P., Paprocki, M., Čertík, O., Kirpichev, S. B., Rocklin, M., Kumar, A., Ivanov, S., Moore, J. K., Singh, S., Rathnayake, T., Vig, S., Granger, B. E., Muller, R. P., Bonazzi, F., Gupta, H., Vats, S., Johansson, F., Pedregosa, F., Curry, M. J., Terrel, A. R., Roučka, v., Saboo, A., Fernando, I., Kulal, S., Cimrman, R., and Scopatz, A.: SymPy: symbolic computing in Python, PeerJ, 3, e103, https://doi.org/10.7717/peerj-cs.103, 2017. a
Morgan, W. J.: Deep mantle convection plumes and plate motions, Am. Assoc. Petr. Geol. B., 56, 203–213, 1972. a
Nakagawa, T. and Tackley, P. J.: The interaction between the post-perovskite phase change and a thermo-chemical boundary layer near the core-mantle boundary, Earth Planet. Sc. Lett., 238, 204–216, 2005. a
Oldham, D. N. and Davies, J. H.: Numerical investigation of layered convection in a three-dimensional shell with application to planetary mantles, Geochem. Geophy. Geosy., 5, Q12C04, https://doi.org/10.1029/2003GC000603, 2004. a
Olson, P., Deguen, R., Hinnov, L. A., and Zhong, S. J.: Controls on geomagnetic reversals and core evolution by mantle convection in the Phanerozoic, Phys. Earth Planet. Int., 214, 87–103, https://doi.org/10.1016/j.pepi.2012.10.003, 2013. a
Popov, I. Yu., Lobanov, I. S., Popov, S. I., Popov, A. I., and Gerya, T. V.: Practical analytical solutions for benchmarking of 2-D and 3-D geodynamic Stokes problems with variable viscosity, Solid Earth, 5, 461–476, https://doi.org/10.5194/se-5-461-2014, 2014. a
Rathgeber, F., Ham, D. A., Mitchell, L., Lange, M., Luporini, F., Mcrae, A. T. T., Bercea, G.-T., Markall, G. R., and Kelly, P. H. J.: Firedrake: Automating the Finite Element Method by Composing Abstractions, ACT T. Math. Softw., 43, 1–24, https://doi.org/10.1145/2998441, 2016. a
Solheim, L. P. and Peltier, W. R.: Avalanche effects in phase transition modulated convection – a model of Earth's mantle, J. Geophys. Res., 99, 6997–7018, 1994. a
Stadler, G., Gurnis, M., Burstedde, C., Wilcox, L. C., Alisic, L., and Ghattas, O.: The dynamics of plate tectonics and mantle flow: from local to global scales, Science, 329, 1033–1038, https://doi.org/10.1126/science.1191223, 2010. a
Stotz, I. L., Iaffaldano, G., and Davies, D. R.: Pressure-Driven Poiseuille Flow: A Major Component of the Torque-Balance Governing Pacific Plate Motion, Geophys. Res. Lett., 45, 117–125, https://doi.org/10.1002/2017GL075697, 2018. a
Tackley, P. J.: Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the Yin-Yang grid, Phys. Earth Planet. Int., 171, 7–18, https://doi.org/10.1016/j.pepi.2008.08.005, 2008. a
Tackley, P. J., Stevenson, D. J., Glatzmaier, G. A., and Schubert, G.: Effects of an endothermic phase transition at 670 km depth in a spherical model of convection in the Earth's mantle, Nature, 361, 699–704, https://doi.org/10.1038/361699a0, 1993. a
Tan, E., Leng, W., Zhong, S., and Gurins, M.: On the location of plumes and mobility of thermo–chemical structures with high bulk modulus in the 3-D compressible mantle, Geochem. Geophy. Geosy., 12, Q07005, https://doi.org/10.1029/2011GC003665, 2011. a
Tosi, N., Stein, C., Noack, L., Hüttig, C., Maierová, P., Samuel, H., Davies, D. R., Wilson, C. R., Kramer, S. C., Thieulot, C., Glerum, A., Fraters, M., Spakman, W., Rozel, A., and Tackley, P. J.: A community benchmark for viscoplastic thermal convection in a 2-D square box, Geochem. Geophy. Geosy., 16, 2175–2196, https://doi.org/10.1002/2015GC005807, 2015. a
Travis, B. J., Anderson, C., Baumgardner, J. R., Gable, C. W., Hager, B. H., O'Connell, R. J., Olson, P., Raefsky, A., and Schubert, G.: A benchmark comparison of numerical methods for infinite Prandtl number thermal convection in two–dimensional Cartesian geometry, Geophys. Astrophys. Fluid Dyn., 55, 137–160, 1990. a
van Keken, P. E.: Cylindrical scaling for dynamical cooling models of the Earth, Phys. Earth. Planet. Int., 124, 119–130, 2001. a
van Keken, P. E. and Ballentine, C. J.: Whole–mantle versus layered mantle convection and the role of a high–viscosity lower mantle in terrestrial volatile evolution, Earth Planet. Sc. Lett., 156, 19–32, 1998. a
van Keken, P. E. and Ballentine, C. J.: Dynamical models of mantle volatile evolution and the role of phase transitions and temperature–dependent rheology, J. Geophys. Res., 104, 7137–7151, 1999. a
van Keken, P. E. and Yuen, D. A.: Dynamical influences of high viscosity in the lower mantle induced by the steep melting curve of perovskite: effects of curvature and time–dependence, J. Geophys. Res., 100, 15233–15248, 1995. a
van Keken, P. E., King, S. D., Schmeling, H., Christensen, U. R., Neumeister, D. and Doin, M. P.: A comparison of methods for the modeling of thermo–chemical convection, J. Geophys. Res., 102, 22477–22495, 1997. a
van Keken, P. E., Kiefer, B., and Peacock, S.: High resolution models of subduction zones: Implications for mineral dehydration reactions and the transport of water into the deep mantle, Geochem. Geophy. Geosy., 3, 1056, https://doi.org/10.1029/2001GC000256, 2002. a
van Keken, P. E., Currie, C., King, S. D., Behn, M. D., Cagnioncle, A., He, J., Katz, R. F., Lin, S., Parmentier, E. M., Spiegelman, M., and Wang, K.: A community benchmark for subduction zone modeling, Phys. Earth Planet. Int., 171, 187–197, https://doi.org/10.1016/j.pepi.2008.04.015, 2008. a
Wilson, C. R., Spiegelman, M., and van Keken, P. E.: TerraFERMA: The Transparent Finite Element Rapid Model Assembler for multiphysics problems in Earth sciences, Geochem. Geophy. Geosy., 18, 769–810, https://doi.org/10.1002/2016GC006702, 2017. a
Wolstencroft, M., Davies, J. H., and Davies, D. R.: Nusselt-Rayleigh number scaling for spherical shell Earth mantle simulation up to a Rayleigh number of 109, Phys. Earth Planet. Int., 176, 132–141, https://doi.org/10.1016/j.pepi.2009.05.002, 2009. a
Zhong, S., Gurnis, M., and Hulbert, G.: Accurate determination of surface normal stress in viscous flow from a consistent boundary flux method, Phys. Earth Planet. Int., 78, 1–8, https://doi.org/10.1016/0031-9201(93)90078-N, 1993. a, b
Zhong, S., Zuber, M. T., Moresi, L., and Gurnis, M.: Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection, J. Geophys. Res., 105, 11063–11082, https://doi.org/10.1029/2000JB900003, 2000. a, b
Zhong, S., McNamara, A., Tan, E., Moresi, L., and Gurnis, M.: A benchmark study on mantle convection in a 3-D spherical shell using CitcomS, Geochem. Geophy. Geosy., 9, Q10017, https://doi.org/10.1029/2008GC002048, 2008. a, b, c, d, e