Analytical solution for viscous incompressible Stokes flow in a spherical shell

I present a new family of analytical flow solutions to the incompressible Stokes equation in a spherical shell. The velocity is tangential to both inner and outer boundaries, the viscosity is radial and of the power-law type, and the solution has been designed so that the expressions for velocity, pressure, and body force are simple polynomials and therefore simple to implement in (geodynamics) codes. Various flow average values, e.g., the root mean square velocity, are analytically computed. This forms the basis of a numerical benchmark for convection codes and I have implemented it in two finite-element codes: ASPECT and ELEFANT. I report error convergence rates for velocity and pressure.


Introduction
The mantle is a very complex region characterized by large variations in temperature, viscosity, composition (Schubert et al., 2001;van Keken et al., 2002), phase changes, melting, and anisotropic structures as revealed by seismic imaging (Lay and Garnero, 2011).As a consequence, the numerical modeling of mantle convection has grown in complexity (e.g., Tackley, 2012;van Heck et al., 2016;Dannberg and Heister, 2016) with the advent of ever more refined modeling techniques and powerful computers.It has also been recognized that the mantle exerts a primary control on the evolution of tectonic plates and that both should be simulated together if one is to build a fully dynamic Earth model (e.g., van Hinsbergen et al., 2011;Bull et al., 2014;Bower et al., 2015;Crameri and Tackley, 2016).
While inter-code comparisons are useful for problems without an analytical solution (e.g., Arrial et al., 2014;Tosi et al., 2015), such benchmark studies rely on the comparisons between a handful of scalar values (e.g., root mean square velocity, Nusselt number), which account for the global character of the solution at steady state but do not lend themselves to error convergence measurements.Fully analytical solutions have also been recently proposed, attempting to represent a mid-ocean ridge (Burstedde et al., 2013) or being more abstract in nature (Zhong, 1996;Popov et al., 2014;Blinova et al., 2016).Actually, any analytical solution to the Stokes equation in three dimensions could be used in a spherical shell provided that the velocity is appropriately applied on the inner and outer boundaries (e.g., Burstedde et al., 2013), but such solutions usually do not satisfy the condition v • n = 0 on the inner and outer boundaries; i.e., there is flow through the boundaries.Not only does the presence of a material flow through the boundaries make the flow not Earthlike, but it also precludes its use for particle-in-cell advection benchmarking (particles would leave and enter the domain).
I propose a new analytical solution for viscous incompressible Stokes flow in a spherical shell.It has been designed with three constraints in mind: (1) the boundary conditions, buoyancy forces, and viscosity fields should be simple to implement and exact; (2) the solution should also be smooth and simple enough to be usable; (3) and it should satisfy tangential slip boundary conditions on both surfaces.I present in Sect. 2 the simplified Stokes equation in spherical coordinates under these flow assumptions and outline the procedure to arrive at the analytical solution for pressure and velocity in Sect.3, for both constant and depth-dependent viscosity profiles.I compute in Sect. 4 the exact analytical values for the root mean square velocity and various other averages.In Sect. 5 the two numerical codes are introduced and the results obtained with these are shown.Finally, these results are discussed in Sect.6.

The Stokes equations in spherical coordinates
The domain is a spherical shell parameterized by its inner radius R 1 and outer radius R 2 .For an incompressible fluid, the Stokes flow equations are given by where v is the velocity vector, ρ is the mass density, g is the gravitational acceleration vector, and σ is the full stress tensor, which can be split: where p is the pressure, 1 is the unit tensor, and s is the deviatoric stress tensor.Equation (2) then becomes In spherical coordinates, Eqs. ( 1) and ( 4) become In this work, the following spherical coordinate conventions are used: r is the radial distance, θ ∈ [0, π ] is the polar angle, and φ ∈ [0, 2π ] is the azimuthal angle.In the case of an incompressible fluid, the deviatoric stress is simply where µ is the dynamic viscosity, which can depend on space coordinates, and ˙ is the (deviatoric) strain rate tensor: In spherical coordinates, the components of the deviatoric stress tensor are given by Equations ( 5)-( 8) supplemented by Eqs. ( 11)-( 16) form a closed set of PDEs that can be solved given an appropriate set of boundary conditions.

Assumptions on the flow
In order to derive an analytical solution for the flow velocity and pressure in the domain, the following assumptions are made.
-All quantities v r , v θ , v φ , p, ρ, and µ are independent of the azimuthal angle φ.As a consequence, all the terms containing partial derivatives with respect to φ can be discarded.This is one of the most stringent limitations in this work since it implies rotational symmetry with respect to the vertical axis.
-The polar and azimuthal components of the velocity are equal and of the form -The radial component of the velocity is null on the inside r = R 1 and outside r = R 2 of the domain, thereby ensuring a tangential flow on the boundaries, i.e., -The viscosity is a function of the radial distance only and takes the form where m is an integer (positive or negative).Note that m = −1 yields a constant viscosity.
-The gravity vector is set to g = −e r and is therefore of unit norm, i.e., |g| = 1.

The set of simplified partial differential equations
Under the above assumptions, it is easy to show that the Stokes equation in spherical coordinates is given by where is the Laplacian operator: 3 Derivation of the analytical solution of the Stokes equations My goal is to determine the f (r) function in the context of all the assumptions made about the flow nature.The following four steps will then be taken: 1. use the continuity Eq. ( 20) to arrive at v r (r, θ ) = g(r) cos θ ; 2. use the θ component of Stokes Eq. ( 22) to arrive at p(r, θ ) = h(r) cos θ ; 3. use the φ component of Stokes Eq. ( 23) to arrive at f (r) and g(r) using boundary conditions on v r ; 4. and use the r component of Stokes Eq. ( 21) to arrive at the density ρ(r, θ ).
3.1 Using the continuity equation to arrive at v r Inserting Eq. ( 17) into Eq.( 20) yields where the integration constant has been set to zero for simplicity.I then proceed to compute the Laplacian of each velocity component using Eq. ( 24), and once again all partial derivatives with respect to φ are neglected: 3.2 Using the θ component of Stokes equations to arrive at the pressure p Inserting Eqs. ( 17) and (25) into Eq.( 22) yields where the integration constant has once again been omitted for simplicity.
3.3 Using the φ component of Stokes equations to arrive at f (r) and g(r) Inserting Eqs. ( 17) and (25) into Eq.( 23) yields I now make use of Eq. ( 19) so that Eq. ( 33) becomes This equation has to hold for all r values so that f (r) is the solution of the second-order ODE: I postulate that the solution is of the form f (r) = r a , which yields The only two acceptable values for a are the roots of this second-order polynomial, i.e., a = 1 or a = −(m + 3).The general solution of the ODE is then where α and β are two constants yet to be determined.Having obtained f (r), one can now compute g(r).However, as will become obvious, one must make a distinction between m = −1 and m = −1.
Note that the case m = −3 should be considered separately since Eq. ( 35) then becomes f (r) = 0 and yields a solution f (r) = βr.

Case
and from Eq. ( 26) one obtains where γ = 0 is a constant.From Eqs. ( 18) and ( 25) it follows that g(R 1 ) = g(R 2 ) = 0, which yields In this case the function g(r) takes the form and the boundary conditions impose Note that this imposes m = −4.

Using the r component of Stokes equations to arrive at density ρ
Equation ( 21) contains the term ∂p/∂r = ∂h/∂r cos θ, which needs to be addressed beforehand: In this case, µ = 0 and µ = 0 so that ∂h/∂r cos θ simplifies to and Eq. ( 21) becomes and by using Eq. ( 27), the radial function F is given by Inserting the radial functions f (r) and g(r) given in Eqs. ( 38) and (39) into Eq.( 49) yields Then Eq. ( 21) becomes I postulate here that and arrive at where I have used Inserting f (r) and g(r) expressions into the above equation yields F(r) so that in the end

The pressure field
The pressure is defined in Eq. ( 32) and h(r) can now be computed.

Case
4 Additional measurements

Root mean square velocity
Many benchmark studies (e.g., Blankenbach et al., 1989;Tosi et al., 2015) report the root mean square velocity quantity, defined as follows: This is a convenient quantity as it captures (in an average sense) the nature of the velocity field in a single scalar value, which allows for an easy comparison either with a known analytical value or across multiple codes.Since the velocity is known in all of the domain, it is a simple although tedious exercise to compute the RMS velocity.I find with The values of A and B depend on the f and g function, so the distinction must once again be made between m = −1 and m = −1.

Radial averages
The radial average of a quantity χ (r, θ, φ) is defined as follows: and due to symmetry it is trivial to show that < p> R = 0 and < v r > R = 0. Likewise, one easily arrives at

Volume averages
The volume average of a quantity χ (r, θ, φ) is defined as follows: Here again, < p> V = 0 and < v r > V = 0, while One can also look at the volume averages of the Cartesian coordinate components of the velocity.
Note that < v x > and < v y > are zero because of symmetry ( 2π 0 cos φdφ = 2π 0 sin φdφ = 0), while < v z >=0 (for all values of m) after tedious calculations using the definitions of α and β.

Surface averages
The average of a quantity χ (r, θ, φ) on a surface of radius R is simply the radial average function evaluated at a given radial distance R. Rather importantly, we have i.e., the average pressure at the outside surface is zero.

Moment of inertia
Because of the expression of the density field, i.e., ρ(r, θ, φ) = F(r) cos θ , it is trivial to show that the moment of inertia of the system with respect to the x, y, and z axis are identically zero.

Stress field
Since the velocity and pressure fields are known, I can also compute an analytical expression for the stress field, and the stress tensor is given by

Implementation and results
As mentioned earlier, this flow solution was designed with a geodynamics application in mind.It has therefore been implemented in the state-of-the-art open-source code ASPECT1 (Kronbichler et al., 2012;Heister et al., 2017) and in the ELEFANT2 code (Thieulot, 2014;Tosi et al., 2015;Lavecchia et al., 2017).Both codes solve the incompressible flow Stokes equations in spherical shell domains but use Cartesian coordinates.

ASPECT
ASPECT is a finite-element code intended to solve the equations that describe thermally driven convection with a focus on doing so in the context of convection in Earth's mantle.The default element type Q 2 Q 1 (quadratic velocity, linear pressure) has been used in this work, but since ASPECT is based on the deal.iilibrary (Bangerth et al., 2007(Bangerth et al., , 2016)), one can easily change the element type from the ASCII input file (".prm" file), and the Q 1 P 0 (linear velocity, constant pressure) was also used (the same element is used in the CIT-COM code; Zhong et al., 2008).I make use of the plugin architecture of the code which allows users and developers to easily add or switch between already-implemented features.The flow velocity and pressure solutions, as well as the viscosity and body force expressions, are all encapsulated in a single piece of code alongside an ASCII input file in which resolution, element type, boundary conditions, and other parameters are set.This benchmark is now part of the mainline since version 2.0.0-pre(see AS-PECT manual).

ELEFANT
ELEFANT borrows largely from the FANTOM code (Thieulot, 2011), but it also brings a number of critical improvements compared to its predecessor, such as spherical shell geometry and the use of a preconditioned conjugate gradient scheme for both inner and outer iterations (Braess, 2007;Elman, 1996).It is a finite-element code based on Q 1 P 0 elements, which is developed and maintained by the author.It has successfully been benchmarked against a wide range of analytical problems and also against other codes (including ASPECT) in the case of visco-plastic thermal convection (Tosi et al., 2015).

Setup
In what follows, I set the inner and outer radii to R 1 = 0.5 and R 2 = 1, respectively.The functions f (r), g(r), h(r), and the viscosity function are shown in Fig. 1.
The ASPECT computational grid consists of 12 blocks making up a hollow sphere (Zhong et al., 2000), with each block subdivided into (2 n ) 3 elements, where n = 2, 3, . ... Since elements vary in size in the radial direction, I have chosen to report the average element size in convergence plots, and it is computed as follows: where V is the volume of the domain and N el is the total number of elements.The ELEFANT mesh is also based on the same 12 blocks, but each block can be subdivided into nel 3 elements, where nel is a positive integer that is not bound to be a power of 2. The analytical velocity solution is prescribed on both internal (r = R 1 ) and external (r = R 2 ) boundaries.Given these boundary conditions, the pressure is determined up to a con-   stant.Both codes allow for a surface normalization (the average pressure along the surface has a prescribed value, in this case zero) and a volume normalization (the average pressure over the whole volume has a prescribed zero value).Since I have shown above that both are identically null for the pressure field, the choice of pressure normalization does not matter.Also, the boundary conditions preclude the presence of a pure rotational mode of numerical origin (Zhong et al., 2008).Both codes were run for values of m = −1 (constant viscosity) and m = 3 (viscosity varies by a factor 16 from the inside to the outside).The density, velocity, and pressure fields for both cases are shown in Fig. 2.
Figure 3 shows the relative root mean square velocity error as a function of the (average) element size for both codes and both m values.The error is found to quadratically decrease with resolution for both codes.Likewise, the L 2 -norm of the   velocity error is found to decrease with resolution, linearly for Q 1 P 0 elements and quadratically for Q 2 Q 1 , as shown in Fig. 4. Looking at the pressure error convergence, it is found to decrease quadratically with the resolution for both types of elements, as shown in Fig. 5. ELEFANT routinely outputs all three average quantities < u >, < v >, and < w >.All three values were found to be zero within machine precision (oscillating around 10 −15 ).

Conclusions
I have derived in this paper a family of analytical solutions to incompressible Stokes flow in a spherical shell under a few assumptions, such as tangential velocity on the boundaries and a radial viscous profile.The velocity, pressure, density, and viscous fields that satisfy the flow equations at every point in space are then used to benchmark two multipurpose geodynamics codes.The L 2 -norms of the velocity and pressure errors were reported and shown to decrease when the resolution is increased.Furthermore, various analytical expressions for flow averages were derived, and it was shown that the computed solutions converged to these expected values.
A number of previous studies (Popov et al., 2014) use an exponential viscosity of the form µ(r) = µ 0 exp(α•r), where α is a parameter controlling the amplitude of the viscosity www.solid-earth.net/8/1181/2017/Solid Earth, 8, 1181-1191, 2017 variations in the system.This approach was tried during the preparation of the paper, but Eq. ( 35) then becomes r 2 f + (2 + mr)rf − (2 + mr)f = 0. (96) Although this equation can be solved, the form of the solution f (r) involves the exponential integral function Ei(r) = − ∞ −r e −t /t dt, which would (a) render the derivations of g(r) and all subsequent quantities very cumbersome and (b) make the solution only semi-analytical.This approach was then abandoned.
Looking at the pressure equation, or rather at Fig. 1c, we see that the pressure is zero for r = R 1 and r = R 2 .One simple interpretation is that the pressure in this work should be interpreted as an over pressure with regards to a background lithostatic pressure.Likewise, the (complex) density profiles have to be interpreted as density variations with regards to a background density profile corresponding to the above-mentioned lithostatic pressure.
Finally, it must be mentioned that the driving density field and flow solution both consist azimuthally of spherical harmonic degree 1 and order 0 combined with an azimuthally constant viscosity.Presumably, the presented solutions are part of a family of solutions that could also be derived for other spherical harmonics.

Figure 3 .
Figure 3. Root mean square velocity as a function of the average element size for both codes.

Figure 4 .
Figure 4. Velocity L 2 -norm error vs. average resolution for ELEFANT and ASPECT.

Figure 5 .
Figure 5. Pressure L 2 -norm error vs. average resolution for ELEFANT and ASPECT.