Analytical solutions for mantle flow in cylindrical and spherical shells

Computational models of mantle convection must accurately represent curved boundaries and the associated boundary conditions within 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 5 source terms that adhere to natural 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 10 codes across the wider community, we provide a python package, Assess, that evaluates the analytical solutions at arbitrary points of the domain.


Introduction
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 is reflected in near-surface phenomena such as plate tectonics, 15 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.
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 20 be hard to choose solutions: (i) that satisfy natural boundary conditions, in particular, in less trivial domains; and (ii) with a velocity that is divergence free, 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 25 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. (2000Zhong et al. ( , 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 30 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. In addition we derive analytical solutions for both the delta-function and smooth forcing cases with either free-slip or zero-slip boundary conditions in a 2-D cylindrical domain. 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 are 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 Section 2 we derive the analytical solutions in cylindrical (section 5 2.2) and spherical (section 2.3) geometries. Smooth solutions are first provided, for free-slip and zero-slip cases, followed by the delta-function solutions. In Section 3, we briefly describe the Fluidity computational modelling framework, which is used in Section 4 to obtain convergence results of the P2-P1 finite element discretisation to the analytical solutions. In Section 5 we discuss these results, in particular an issue with the delta-function cases for discretisations that use continuous pressure. To demonstrate, we also show results with a P2 bubble -P1 DG discretisation that overcomes this issue.

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 (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 − ≤ r ≤ R + , where r is the radial distance to 20 the origin.

Cylindrical
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, x = r cos(ϕ), y = r sin(ϕ), and u r and u ϕ are the radial and trans-25 verse components of velocity respectively. The normal deviatoric stress and shear stress components are given by The momentum equation (1) can be decomposed into radial and transverse components ν r −ν ∂ ∂r where the Laplacian, in polar coordinates, is given by For a derivation of Equations 5-8, see appendix A.
The curl of the momentum equation is obtained by summation of the operators − 1 r ∂ q ∂ϕ and 1 r ∂ ∂r (r q ) applied to the radial (7) and transverse (8) components respectively, which leads to The general, real-valued, solution to the biharmonic equation ∇ 4 ψ = 0 is given by 10 ψ(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(ϕ)) where A n , B n , C n , D n , e n , and f n are constant coefficients, n ≥ 0. In the following we will, for simplicity, focus on sin(nϕ)solutions for a single n > 1 and set e n = 1, f n = 0, ∀n.
An equation for pressure can be derived by taking the divergence of the momentum equation The homogeneous solutions, i.e. ρ = 0, are thus the standard harmonic solutions where G n and H n are constant coefficients, n > 1. From (7) it can be seen that sin(nϕ) solutions for ψ are associated with cos(nϕ) solutions of p and ρ . After substitution of p and ψ in (7)-(8), the following relation 20 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.

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 (10) exists of the form 5 assuming k = n − 3 and k = n − 1, so that a general solution can be written as ψ(r, ϕ) = Ar n + Br −n + Cr n+2 + Dr −n+2 + Er 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 (12) is given by 10 The general solution for pressure is thus p(r, ϕ) = Gr n + Hr −n + F r k+1 cos(nϕ) with G and H given by (14).
The four remaining coefficients A, B, C, and D, are fixed by a choice of boundary conditions at the top and bottom of the cylindrical domain at r = R + and r = R − , respectively. At both, no normal flow conditions are imposed via ∂ψ ∂ϕ = 0. Two 15 further equations are found by imposing either τ rϕ = 0 (free slip), or ∂ψ ∂r = 0 (no slip) at both boundaries. The solution for free slip, no normal flow at both boundaries is given by

Green's function solution -cylindrical
Another set of useful solutions is found considering the following perturbation density representing an infinitely thin density anomaly at r = r with R − ≤ r ≤ R + . Since ρ = 0 for r = r, the solution is described by combining two homogeneous solutions We thus have 8 coefficients that are fixed by boundary conditions and conditions at the interface. Boundary conditions at the top and bottom 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 a continuity of the shear stress (6) which in combination 15 with the above implies Finally the eighth equation is obtained by integrating (10) over a small strip r − ≤ r ≤ r + https://doi.org /10.5194/gmd-2020-194 Preprint. Discussion started: 24 August 2020 c Author(s) 2020. CC BY 4.0 License.
Taking the limit → 0, and again using (22) and (23), we obtain The solution with free slip boundary conditions at r = R − and r = R + is where we use α ± = R ± /r . No slip conditions at these boundaries leads to where in addition we use γ = R − /R + .

Spherical
In this section we derive the equivalent of the four cylindrical cases in a 3-D spherical domain, R − ≤ r ≤ R + . A more detailed 15 derivation of the equations can be found in Ribe (2009). In 3-D, an incompressible velocity field u can be decomposed as using poloidal and toroidal scalar fields P and T where r = rr. In spherical coordinates where ϕ is the longitude and θ is the co-latitude, and Figure 1. Density 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 top boundary.
We further derive the components of stress normal to a spherical surface Using these, we can work out spherical components of the momentum equation (1) to be As can be seen from these equations, the toroidal part of the solution is independent of the density distribution. A nonzero toroidal component is only introduced through a particular choice of boundary conditions. For the test cases in this paper we 10 only consider free slip and no-slip boundary conditions with radial forcing so we may assume T = 0.
Taking the curl of the momentum equation In a connected spherical domain this implies that ν∇ 4 P − gρ /r varies in the radial direction only. Without loss of generality we may therefore seek solutions P that satisfy 5 since for any other solution P + P of (36) we know that ∇ 4 P is purely radial, and therefore P can be written as a sum of biharmonic solutions and a purely radial function. The biharmonic solutions are included in (37) and the purely radial function is discarded by (27).
Solutions to the biharmonic equation ∇ 4 P = 0 in 3-D are given by 10 where A lm , B lm , C lm , and D lm are constant coefficients, and Y lm are 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 15 Homogeneous solutions for pressure are written where substitution of (38) and (41) in (33) gives 2.3.1 Smooth density profile -spherical 20 We consider a density perturbation of the following form An inhomogeneous solution of (37) is given by with a more generic solution written as where we have dropped the lm subscripts of the coefficients A, B, C, D. The pressure solution for (40) is and G and H are given by (42).

5
As before, the four coefficients A, B, C, and D are fixed by the boundary conditions no normal flow: where we use Λ 2 P = −l(l + 1)P. The two no-normal flow conditions at r = R − and r = R + are combined with two further conditions free slip: 10 no-slip: https://doi.org /10.5194/gmd-2020-194 Preprint. Discussion started: 24 August 2020 c Author(s) 2020. CC BY 4.0 License.
For free slip conditions at both r = R − and r = R + , the solution is given by No slip conditions at both boundaries lead to the following solution

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 4 constraints derived from the boundary conditions at r = R − and r = R + as in the previous section, and by imposing a further 4 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 of u r : P − (r , θ, ϕ) = P + (r , θ, ϕ) (52) 20 continuity of u θ and u ϕ : zero shear condition: ∂ 2 P − ∂r 2 (r , θ, ϕ) = ∂ 2 P + ∂r 2 (r , θ, ϕ), normal shear condition: where (55) is derived from (37) in the same way as (26), and (53)-(55) assume (52).
The free slip solution is given by The no-slip solution is given by

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

15
The numerical solutions for velocity and pressure, u and p, are written as a linear combination of basis functions N j and M l We use either the P2-P1 (Taylor Hood) or P2 bubble -P1 DG element pairs. In both cases the velocity and pressure basis functions The Stokes equations are written in the weak form, using the same N j and M l basis as test-functions. After integrating by parts (1) and (3) then become (omitting boundary terms) Note that we apply strong Dirichlet boundary conditions, so that the boundary integrals can indeed be neglected. In the free 5 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).

Isoparametric representation of the domain
For an accurate representation of the quadratic approximation of velocity, we need to also approximate the curved cylin-10 drical/spherical domain quadratically. This means that rather than each cell in the mesh being described by a linear map X lin : ξ → X lin (ξ) from local coordinates ξ in a reference element to physical coordinates X, we use a quadratic map X quad (ξ), which maps to a curved triangle/tetrahedron that better represents the domain. This map can be obtained from a linear mesh with coordinate mappings, X lin , through quadratic interpolation 15 at the standard Lagrange nodes of the quadratic function space. Here r lin (ξ) is the linear interpolation of the radius, i.e. r lin = X lin 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.

Forcing term
The density perturbation ρ on the RHS of (57) is a prescribed analytical expression in each of the test cases. For the Green's 20 function solutions in 2-D and 3-D, using (20) and (50) respectively, we get where Γ is an internal boundary at r = r oriented such that its normal n =r points outwards.

Solving the linear system and dealing with nullspaces
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). 5 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 (−y, x) = rθ to the velocity solution. We 10 therefore apply the following projection to the numerical solution 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 modes lead to zero modes (eigenvectors) for the linear system based on (57) and (58), rendering the resulting matrix singular. In preconditioned Krylov 15 methods we typically need to subtract the zero modes from the approximate solution at every iteration. With iterative approximation x i and zero eigenvector λ we get This l 2 -projection, based on the l 2 inner product ·, · , is analogous but not equivalent to the projections in (62) and (63) (which are L 2 -projections). The L 2 -projections should, therefore, be applied as a separate step after the iterative solvers have 20 completed.

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 25 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 3D mesh consisting of 61440 tetrahedra. Again, resolution is doubled in all directions for subsequent refinement levels. 15 https://doi.org /10.5194/gmd-2020-194 Preprint. Discussion started: 24 August 2020 c Author(s) 2020. CC BY 4.0 License. Figure 4. Convergence 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.
Free Slip Zero Slip Figure 5. Convergence 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 relative to comparable cases in Fig. 3.
In all figures, errors are given as relative errors, comparing the numerical solution, u and p, with the analytical solutions, u * and p * in the L 2 -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 5 zero-slip boundary conditions. Cases with lower wave-number 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 10 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
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 10 dual of the Sobolev space H 1 (Ω). This means that the delta function can be thought of as a continuous function which maps functions v ∈ H 1 (Ω), the space of square integrable functions with square integrable weak derivatives, to R. Girault and Raviart (2012) demonstrate that even with the very loose regularity condition that the right-hand side f is in For velocity-pressure finite element pairs that satisfy the standard inf-sup, or 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 20 method converges and in fact where C 1 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 25 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 a discontinuous function, such as our analytical pressure solution p * , the best approximation by continuous, piecewise-linear polynomials, i.e. W = P1, is bounded by This therefore limits the convergence of the method.
A solution to this problem is found by allowing for discontinuities in the discrete pressure space. We demonstrate this here by considering W = P1 DG the space of piecewise linear but discontinuous functions. To satisfy the inf-sup condition this 5 requires enriching the quadratic function space P 2 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) and Burstedde et al. (2013), who ran the same spherical cases with a delta function forcing, and found second order convergence for velocity and second order convergence for 10 pressure related diagnostics. It should be noted, however, that the aforementioned 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) were obtained with CitcomS based on a Q1-P0 (continuous piecewise trilinear velocity and piecewise constant discontinuous pressure). Although our analysis above limits the convergence of p − p * 2 for a P0 pressure p to first order, second order super-convergence can be obtained in some cases by evaluating the analytical solution 5 only in the cell centre. In other words, by comparing to a filtered piecewise constant analytical approximationp * , second order convergence can sometimes be observed in p −p * 2 . For continuous pressure approximations, such as the P2-P1 results in this paper, and the Q1-Q1 discretisation used in Burstedde et al. (2013), reduced convergence in the interior of the domain is to be expected.

10
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 available from https://doi.org/10.5281/zenodo.3891545. To ensure correctness, both in the manuscript as in the python package, the coefficients for the various solution are extracted from the L A T E X-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/. 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.

15
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.
Acknowledgements. DRD and SCK acknowledge support from the Australian Research Council, under grant numbers FT140101262 and 20 DP170100058. Numerical simulations were undertaken on the NCI National Facility in Canberra, Australia, which is supported by the Australian Commonwealth Government.