Articles | Volume 14, issue 4
Methods for assessment of models
09 Apr 2021
Methods for assessment of models |  | 09 Apr 2021

Analytical solutions for mantle flow in cylindrical and spherical shells

Stephan C. Kramer, D. Rhodri Davies, and Cian R. Wilson

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.

1 Introduction

Mantle convection transports Earth's internal heat to its surface: it is the “engine” driving our dynamic Earth (e.g. Davies1999). 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. Morgan1972; Mitrovica et al.1989; Gurnis1993; 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 Davies1986; Davies and Stevenson1992; Moresi and Solomatov1995; 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. Baumgardner1985; Tackley et al.1993; Bunge et al.1996, 1997; Zhong et al.2000; Oldham and Davies2004; McNamara and Zhong2005; Choblet et al.2007; Zhong et al.2008; Tackley2008; Davies and Davies2009; 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 Gassmoller2018; 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 Peltier1994; van Keken and Yuen1995), the 2-D cylinder (e.g. Jarvis1993; van Keken and Ballentine1998, 1999; van Keken2001; Nakagawa and Tackley2005), and the spherical annulus (e.g. Hernlund and Tackley2008). 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; Roache2002). 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 Richards1989). 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 (Backus1986). 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; Thieulot2017; 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; Kramer2020).

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.

2 Analytical Solutions

2.1 Equations

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 -r^, and a dimensionless density deviation ρ: ρ=ρ0(1+ρ), 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 R-rR+, where r is the radial distance to the origin.

2.2 Cylindrical

In 2-D, any incompressible velocity field u can be written as the skew gradient of a streamfunction ψ:

(4) u r = - 1 r ψ φ , u φ = ψ r ,

where φ is the angle from the x axis in polar coordinates, x=rcos(φ),y=rsin(φ), 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

(9) 2 ψ = 1 r r r ψ r + 1 r 2 2 ψ φ 2 .

For a derivation of Eqs. (5)–(8), see Appendix A.

The curl of the momentum equation is obtained by summation of the operators -1rφ and 1rrr applied to the radial (Eq. 7) and transverse (Eq. 8) components, respectively, which leads to

(10) - 4 ψ = g ν 1 r ρ φ .

The general, real-valued solution to the biharmonic equation 4ψ=0 is given by

(11) ψ ( r , φ ) = n > 1 A n r n + B n r - n + C n r n + 2 + D n r - n + 2 e n sin ( n φ ) + f n cos ( n φ ) + A 1 r + B 1 r - 1 + C 1 r 3 + D 1 r ln r ( e 1 sin ( φ ) + f 1 cos ( φ ) ) + A 0 + B 0 ln r + C 0 r 2 + D 0 r 2 ln r ,

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=∇22. In the following, we will, for simplicity, focus on sin (nφ) solutions for a single n>1 and set en=1,fn=0.

An equation for pressure can be derived by taking the divergence of the momentum equation:

(12) 2 p = - g r r ρ r .

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. ρ=0, are thus the standard harmonic solutions (again neglecting the n=0,1 solutions)

(13) p ( r , φ ) = n > 1 G n r n + H n r - n cos ( n φ ) ,

where Gn and Hn are constant coefficients, n>1.

After substitution of p and ψ in Eqs. (7)–(8), the following relations

(14) G n = - 4 ν C n ( n + 1 ) , H n = - 4 ν D n ( 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:

(15) ρ = r k R + k cos ( n φ ) ,

with k>0. It is easily verified that an inhomogeneous solution to Eq. (10) exists of the form

(16) ψ = E r k + 3 sin ( n φ ) , E = g R + - k n ν ( k + 3 ) 2 - n 2 ( k + 1 ) 2 - n 2 ,

assuming kn-3 and kn-1, so that a general solution can be written as

(17) ψ ( r , φ ) = A r n + B r - n + C r n + 2 + D r - n + 2 + E r k + 3 sin ( n φ ) .

Note that for the remainder of this derivation we drop the subscript n in the coefficients for A,B,C, and D.

An inhomogeneous solution for pressure of Eq. (12) is given by

(18) p ( r , φ ) = F r k + 1 cos ( n φ ) , F = - g R + - k ( k + 1 ) ( k + 1 ) 2 - n 2 .

The general solution for pressure is thus

(19) p ( r , φ ) = G r n + H r - n + F r k + 1 cos ( n φ ) ,

with G and H given by Eq. (14).

The four remaining coefficients (A,B,C, and D) are fixed by a choice of boundary conditions at the inner and outer surfaces of the cylindrical domain at r=R+ and r=R-, respectively. At both, no-normal-flow conditions are imposed via ψφ=0. Two further equations are found by imposing either τrφ=0 (free slip) or ψr=0 (zero slip) at both boundaries.

The solution coefficients for free-slip, no normal flow at both boundaries are given by


where we use α=R-/R+.

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,

(20) ρ = δ ( r - r ) cos ( n φ ) ,

representing an infinitely thin density anomaly at r=r with R-<r<R+. Since ρ=0 for rr, the solution is described by combining two homogeneous solutions

(21) ψ ( r , φ ) = ψ - ( r , φ ) = A - r n + B - r - n + C - r n + 2 + D - r - n + 2 sin ( n φ )  for  R - r < r , ψ + ( r , φ ) = A + r n + B + r - n + C + r n + 2 + D + r - n + 2 sin ( n φ )  for  r < r R + .

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:

(22) ψ - ( r , φ ) = ψ + ( r , φ ) ,  and  ψ - r ( r , φ ) = ψ + r ( r , φ ) .

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

(23) 2 ψ - r 2 ( r , φ ) = 2 ψ + r 2 ( r , φ ) .

Finally, the eighth equation is obtained by integrating Eq. (10) over a small strip r-ϵrr+ϵ, between two arbitrary angles φ1φφ2:


Taking the limit ϵ→0, the second term on the left-hand side disappears, whereas the first term becomes


Again using Eqs. (22) and (23), only the jump term for the third radial derivative of ψ± remains:


Thus, for this to hold for arbitrary φ1 and φ2, we need

(26) 3 ψ + r 3 ( r , φ ) - 3 ψ - r 3 ( r , φ ) = g n sin ( n φ ) ν r .

The solution coefficients for free-slip boundary conditions at r=R- and r=R+ are


where we use α±=R±/r.

Zero-slip conditions at both boundaries lead to


where, in addition, we use γ=R-/R+.

Figure 1Density perturbation (δρ) field for smooth cylindrical cases across a range of n and k. As the wavenumber n increases, the wavelength of the perturbation decreases. As k increases, the perturbation becomes more concentrated towards the domain's outer boundary.


2.3 Spherical

In this section, we derive the equivalent of the four cylindrical cases in a 3-D spherical domain, R-rR+. 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

(27) u = × r × P + r × T

using poloidal and toroidal scalar fields 𝒫 and 𝒯 (Backus1986), where r=rr^. In spherical coordinates,

(28) u r = 1 r Λ 2 P , u θ = - 1 r 2 r P r θ - 1 sin ( θ ) T φ , u φ = - 1 r sin ( θ ) 2 r P r φ + T θ ,

where φ is the longitude and θ is the co-latitude, and

(29) Λ 2 P = 1 sin ( θ ) θ sin ( θ ) P θ + 1 sin ( θ ) 2 2 P φ 2 ,  so that  2 = 1 r 2 r r 2 r + Λ 2 .

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,

(36) - θ ^ ν sin ( θ ) 4 P φ + φ ^ ν 4 P θ = - θ ^ 1 r sin ( θ ) g ρ φ + φ ^ 1 r g ρ θ .

In a spherical domain, this implies that ν4P-gρ/r varies in the radial direction only. Without loss of generality, we may therefore seek solutions 𝒫 that satisfy

(37) 4 P = g ρ ν r ,

since for any other solution P+P 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

(38) P ( r , θ , φ ) = l 0 m = - l l A l m r l + B l m r - l - 1 + C l m r l + 2 + D l m r - l + 1 Y l m ( θ , φ ) ,

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

(39) Λ 2 Y l m = - l ( l + 1 ) Y l m .

As in 2-D, an equation for the pressure is obtained by taking the divergence of the momentum equation, giving

(40) 2 p = - g r 2 r r 2 ρ .

Homogeneous solutions for pressure are written as

(41) p ( r , θ , φ ) = l 0 m = - l l G l m r l + H l m r - l - 1 Y l m ( θ , φ ) ,

where the substitution of Eqs. (38) and (41) in Eq. (33) gives

(42) G l m = - 2 ν ( l + 1 ) ( 2 l + 3 ) C l m , H l m = - 2 ν l ( 2 l - 1 ) D l m .

Figure 2Illustrations of the density perturbation (δρ) field for smooth spherical cases across a range of l, m, and k.


2.3.1 Smooth density profile – spherical

We consider a density perturbation of the following form:

(43) ρ = r k R + k Y l m ( θ , φ ) .

An inhomogeneous solution of Eq. (37) is given by

(44) P = E r k + 3 Y l m , E = g R + - k / ν ( k + 1 ) ( k + 2 ) - l ( l + 1 ) ( k + 3 ) ( k + 4 ) - l ( l + 1 ) ,

with a more generic solution written as

(45) P = A r l + B r - l - 1 + C r l + 2 + D r - l + 1 + E r k + 3 Y l m ,

where we have dropped the lm subscripts of coefficients A, B, C, and D. The pressure solution for Eq. (40) is

(46) p = G r l + H r - l - 1 + F r k + 1 Y l m , F = - g ( k + 2 ) R + - k ( k + 1 ) ( k + 2 ) - l ( l + 1 ) ,

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:

no normal flow: (47)ur=1rΛ2P=0P=0,at r=R±,

where we use Λ2P=-l(l+1)P. The two no-normal-flow conditions at r=R- and r=R+ are combined with two further conditions:

 free slip: τrθ=τrφ=01rΛ2P-2Pr2(48)+2r2P=02Pr2=0,at r=R±, or  zero slip: uθ=uφ=0(rP)r=(49)0Pr=0,at r=R±.

For free-slip conditions at both r=R- and r=R+, 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

(50) ρ = δ ( r - r ) Y l m ( θ , φ ) ,

by combining two homogeneous solutions:

(51) P ( r , θ , φ ) = P - ( r , θ , φ ) = A - r l + B - r - l - 1 + C - r l + 2 + D - r - l + 1 Y l m ( θ , φ )  for  R - r < r , P + ( r , θ , φ ) = A + r l + B + r - l - 1 + C + r l + 2 + D + r - l + 1 Y l m ( θ , φ )  for  r < r R + .

The eight coefficients are found by imposing the same four constraints derived from the boundary conditions at r=R- and r=R+ 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:

continuity ofur:(52)P-(r,θ,φ)=P+(r,θ,φ), continuity of uθ and uφ:(rP-)r|r=r=(rP+)r|r=rP-r(r,θ,φ)(53)=P+r(r,θ,φ),zero-shear condition: (54)2P-r2(r,θ,φ)=2P+r2(r,θ,φ),normal-shear condition: (55)3P+r3(r,θ,φ)-3P-r3(r,θ,φ)=gYlm(θ,φ)νr,

where Eq. (55) is derived from Eq. (37) in the same way as Eq. (26), and Eqs. (53)–(55) assume Eq. (52).

The free-slip solution coefficients are given by


The zero-slip solution coefficients are given by

3 Fluidity

The test cases in the previous section have been examined using Fluidity, a finite element, control-volume computational modelling framework (Davies et al.2011; Kramer et al.2012).

3.1 Discretisation

The numerical solutions for velocity and pressure, u and p, are written as a linear combination of basis functions Nj and Ml:

(56) u = j u j N j , p = l p l M l .

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 Guermond2004; Boffi et al.2013).

The Stokes equations are written in the weak form, using the same Nj and Ml basis as test functions. After integrating by parts, Eqs. (1) and (3) then become (omitting boundary terms)

ΩhνNijNjTuj+ujTNj(57)+NilplMl=-ΩhNigρr^for all Ni,(58)ΩhMkjujNj=0for all Mk.

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

Figure 3Convergence for 2-D cylindrical cases with free-slip and zero-slip boundary conditions at a series of different wavenumbers, n, as indicated in the legend. Note that cases with a smooth forcing are run at k=2 and k=8, as indicated. Convergence rate is indicated by dashed lines, with the order of convergence provided in the legend.


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 Xlin:ξXlin(ξ) 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:

(59) X quad ( ξ ) quad. interp. r lin ( ξ ) X lin ( ξ ) X lin ( ξ ) ,

at the standard Lagrange nodes of the quadratic function space. Here, rlin(ξ) is the linear interpolation of the radius, i.e. rlin=Xlin 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.

Figure 4Convergence of velocity and pressure for 3-D spherical cases with free-slip and zero-slip boundary conditions, at a range of degrees l and orders m. Note that all cases with a smooth forcing are run at k=l+1.


Figure 5Convergence of velocity (a, b) and pressure (c, d), respectively, for smooth (k=2) 2-D cylindrical cases with free-slip and zero-slip boundary conditions. Note that these cases do not incorporate an isoparametric approximation of the domain; hence, the reduced convergence is relative to comparable cases in Fig. 3.


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 r=r oriented such that its normal n=r^ 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:

(62) p p - Ω h p Ω h 1 .

Similarly, for free-slip cases, in 2-D, we may add an arbitrary rotation of the form (-y,x)=rθ^ to the velocity solution. We therefore apply the following projection to the numerical solution:

(63) u u - Ω h r θ ^ u Ω h r 2 r θ ^ ,

which ensures that the angular momentum rθ^u 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

(64) x i x i - λ , x i λ .

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.

4 Convergence Results

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 128×16×2=4096 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 R-=1.22 and R+=2.22, and the delta-function cases used r=(R-+R+)/2. 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:

(65) u - u 2 u 2 = Ω h | u - u | 2 Ω h | u | 2 , p - p 2 p 2 = Ω h ( p - p ) 2 Ω h ( p ) 2 .

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 k=l+1 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 Discussion

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 King2019). 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:

(66) δ r : v Ω δ ( r - r ) v ( r , ϕ ) r d r d ϕ = Ω v ( r ) r d ϕ ,

which maps functions vH1(Ω), 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 m-1, 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 m=-1 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

(67) sup v V Ω v q v 1 β q 2  for all  q W ,

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 ||1 is the semi-norm

(68) | u | 1 = i = 1 dim Ω i u 2 .

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

(69) inf q W p - q 2 C 2 h 1 2 p 2

(this bound can be derived from the order of convergence results in Bernardi1989). This therefore limits the convergence of the method.

Figure 6Convergence of velocity (a, b) and pressure (c, d), respectively, for delta-function 2-D cylindrical cases with free-slip and zero-slip boundary conditions with the P2bubble–P1DG element pair.


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 Guermond2004; 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.

Figure 7Convergence of velocity and radial stresses at top (a, c) top and (b, d) bottom surfaces (r=R±), respectively, for delta-function 3-D spherical cases with free-slip boundary conditions.


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 p-p2 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 p¯, second-order convergence can sometimes be observed in p-p¯2. 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.

6 Conclusions

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; Kramer2020).

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.

Appendix A: Equations in polar coordinates

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

(A1) u r = - 1 r ψ φ , u φ = ψ r .

We make use of the following expressions for the derivatives of the unit vectors r^ and φ^ with respect to r and φ:


Using these, we can work out the different components of stress:


Note that, as expected, τrr+τφφ=0. In the same way, we derive the following expression for the vorticity, or curl of the velocity:

(A7) curl u = r ^ u φ ^ - φ ^ u r ^ = 2 ψ r 2 + 1 r φ 1 r ψ φ + 1 r ψ r = 2 ψ .

The viscosity term in the Stokes equations can be written as


In addition to Eq. (A3), we use the following identities:

(A10) r ^ = 1 r , φ ^ = 0 ,

and the fact that τφφ=-τrr. After reordering to group the radial and transverse components, this leads to


In combination with the pressure gradient term in polar coordinates

(A14) p = p r r ^ + 1 r p φ φ ^ ,

we obtain the radial and transverse components of the Stokes momentum equation in Eqs. (7) and (8).

Code availability

The Python package Assess, which implements the analytical solutions and evaluates them at arbitrary locations in the domain, is available from (last access: 12 June 2020) (see (last access: 12 June 2020). An archived version is available from (Kramer2020). 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 Developers2019), 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 (last access: 17 August 2020). All cases in this paper have been run with tag version 4.1.17, which is archived at (Kramer et al.2020).

Author contributions

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.

Competing interests

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.

Financial support

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

Review statement

This paper was edited by Thomas Poulet and reviewed by Marcus Mohr and one anonymous referee.


Backus, G.: Poloidal and toroidal fields in geomagnetic field modeling, Rev. Geophys., 24, 75–109, 1986. a, b

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

Baumgardner, J. R.: Three-dimensional treatment of convective flow in the Earth's mantle, J. Stat. Phys., 39, 501–511,, 1985. 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,, 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.: The effect of depth–dependent viscosity on the planform of mantle convection, Nature, 279, 436–438,, 1996. 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,, 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,, 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,, 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,, 2007. a, b

Dannberg, J. and Gassmoller, R.: Chemical trends in ocean islands explained by plume-slab interaction, P. Natl. Acad. Sci. USA, 115, 4351–4356,, 2018. a

Davies, D. R. and Davies, J. H.: Thermally–driven mantle plumes reconcile multiple hotspot observations, Earth Planet. Sc. Lett., 278, 50–54,, 2009. a

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,, 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,, 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,, 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,, 2019. a

Davies, G. F.: Dynamic Earth: plates, plumes and mantle convection, Cambridge University Press, Cambridge, 1999. a

Davies, J. H. and Stevenson, D. J.: Physical Model of Source Region of Subduction Zone Volcanics, J. Geophys. Res., 97, 2037–2070,, 1992. a

Ern, A. and Guermond, J.-L.: Theory and Practice of Finite Elements, Springer, Princeton, 2004. a, b

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,, 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

Gurnis, M.: Phanerozoic marine inundation of continents driven by dynamic topography above subducting slabs, Nature, 364, 589–593,, 1993. a

Gurnis, M. and Davies, G. F.: Mixing in numerical-models of mantle convection incorporating plate kinematics, J. Geophys. Res., 91, 6375–6395,, 1986. 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,, 2017. a, b

Hernlund, J. W. and Tackley, P. J.: Modeling mantle convection inside the spherical annulus, Phys. Earth Planet. Int., 171, 48–54,, 2008. a

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

Horbach, A., Mohr, M., and Bunge, H.-P.: A semi-analytic accuracy benchmark for Stokes flow in 3-D spherical mantle convection codes, GEM-International J. Geomath., 11, 1–35, 2020. a, b, c, d

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,, 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,, 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,, 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,, 2019. a

Kramer, S. C.: Assess (v1.0), Analytical Solutions for the Stokes Equations in Spherical Shells in python, Zenodo,, 2020. a, b, c

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,, 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,, 2020. a

Kronbichler, M., Heister, T., and Bangerth, W.: High accuracy mantle convection simulation through modern numerical methods, Geophys. J. Int., 191, 12–29,, 2012. a, b

Liu, S. and King, S. D.: A benchmark study of incompressible Stokes flow in a 3-D spherical shell using ASPECT, Geophys. J. Int., 217, 650–667, 2019. a, b, c

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,, 2012. a

McKenzie, D. P., Roberts, J. M., and Weiss, N. O.: Convection in the Earth's mantle: towards a numerical simulation, J. Fluid Mech., 62, 465–538,, 1974. a

McNamara, A. K. and Zhong, S.: Thermo–chemical structures beneath Africa and the Pacific Ocean, Nature, 437, 1136–1139,, 2005. 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,, 2017. a

Mitrovica, J. X., Beaumont, C., and Jarvis, G. T.: Tilting of continental interiors by the dynamical effects of subduction, Tectonics, 8, 1079–1094,, 1989. a

Moresi, L. N. and Solomatov, V. S.: Numerical investigations of 2D convection with extremely large viscosity variations, Phys. Fluid, 7, 2154–2162,, 1995. 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,, 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,, 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,, 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,, 2016. a

Ribe, N. M.: Analytical Approaches to Mantle Dynamics, in: Mantle Dynamics, edited by: Bercovici, D. and Schubert, G., vol. 7 of Treatise on Geophysics, Elsevier, New York, 167–226,2009. a, b

Roache, P. J.: Code Verification by the Method of Manufactured Solutions, J. Fluid. Eng.-T. ASME, 124, 4–10,, 2002. 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,, 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,, 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,, 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,, 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,, 2011. a

The Sage Developers: SageMath, the Sage Mathematics Software System (Version 8.6), available at: (last access: 1 June 2020), 2019. a

Thieulot, C.: Analytical solution for viscous incompressible Stokes flow in a spherical shell, Solid Earth, 8, 1181–1191,, 2017. a, b

Tosi, N. and Martinec, Z.: Semi-analytical solution for viscous Stokes flow in two eccentrically nested spheres, Geophys. J. Int., 170, 1015–1030, 2007. a, b

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,, 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,, 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,, 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,, 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,, 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,, 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,, 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,, 2008. a, b, c, d, e

Short summary
Computational models of Earth's mantle require rigorous verification and validation. Analytical solutions of the underlying Stokes equations provide a method to verify that these equations are accurately solved for. However, their derivation in spherical and cylindrical shell domains with physically relevant boundary conditions is involved. This paper provides a number of solutions. They are provided in a Python package (Assess) and their use is demonstrated in a convergence study with Fluidity.