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 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. Davies, 1999). The structure, composition, and flow regime within the mantle are reflected in near-surface phenomena such as plate tectonics, mountain building, dynamic topography, sea-level change, volcanism, and the activity of Earth's magnetic field (e.g. Morgan, 1972; Mitrovica et al., 1989; Gurnis, 1993; Olson et al., 2013; Kloecking et al., 2019; Davies et al., 2019). The grand challenge is to understand the operation of this giant heat engine over geologic time and its relationship to the surface geological record.

Computational modelling is one of the primary tools available for tackling this challenge. Whilst 2- and 3-D numerical models of mantle convection processes in Cartesian domains have provided important insights into a range of mantle processes (e.g. McKenzie et al., 1974; Gurnis and Davies, 1986; Davies and Stevenson, 1992; Moresi and Solomatov, 1995; van Keken et al., 2002; Hunt et al., 2012; Garel et al., 2014; Davies et al., 2016; Jones et al., 2016, 2019), 3-D spherical geometry is required to simulate global mantle dynamics. Global 3-D spherical mantle convection models, and studies focussing on their application, are now in common use (e.g. Baumgardner, 1985; Tackley et al., 1993; Bunge et al., 1996, 1997; Zhong et al., 2000; Oldham and Davies, 2004; McNamara and Zhong, 2005; Choblet et al., 2007; Zhong et al., 2008; Tackley, 2008; Davies and Davies, 2009; Wolstencroft et al., 2009; Stadler et al., 2010; Tan et al., 2011; Kronbichler et al., 2012; Davies et al., 2013; Burstedde et al., 2013; Heister et al., 2017; Dannberg and Gassmoller, 2018; Stotz et al., 2018). However, the use of this geometry for calculations at a realistic convective vigour remains expensive. As a consequence, simplifying geometries are often used, including the axisymmetric spherical shell (e.g. Solheim and Peltier, 1994; van Keken and Yuen, 1995), the 2-D cylinder (e.g. Jarvis, 1993; van Keken and Ballentine, 1998, 1999; van Keken, 2001; Nakagawa and Tackley, 2005), and the spherical annulus (e.g. Hernlund and Tackley, 2008). Such 2-D models can also provide a rapid and broad parameter space appraisal, allowing one to focus on a targeted and more sparse study in a full 3-D spherical geometry.

Recent decades have seen extensive validation, verification, and benchmarking of Cartesian mantle dynamics codes, in both 2- and 3-D. Verification is typically achieved via comparisons of numerical predictions against analytical solutions (e.g. Zhong et al., 1993; Kramer et al., 2012), with benchmarking typically undertaken against solutions from other comparable codes, for incompressible (Blankenbach et al., 1989; Travis et al., 1990; Busse et al., 1994; van Keken et al., 1997, 2008; Tosi et al., 2015) and compressible (King et al., 2009) convection. See also Popov et al. (2014) for an overview of geodynamical benchmarks. The number of comparable studies, within a 2-D cylindrical or 3-D spherical geometry, however, is more limited. Given a recent surge in the state-of-the-art tools available for simulating mantle dynamics in these geometries (e.g. Kronbichler et al., 2012; Logg et al., 2012; Burstedde et al., 2013; Rathgeber et al., 2016; Heister et al., 2017; Wilson et al., 2017), it is important that these tools are verified and validated against a range of analytical and benchmark solutions.

A popular method to obtain analytical solutions is the method of manufactured solutions (MMS; Roache, 2002). With this approach, an arbitrary analytical solution is chosen beforehand and the necessary forcing terms on the right-hand side are derived by substitution of the solution into the left-hand side of the flow equations. A drawback of this approach is that it can be hard to choose solutions (i) that satisfy physically realistic boundary conditions, in particular, in less trivial domains; and (ii) with a velocity that is divergence free, in order to avoid unnatural source terms. Alternatively, one can choose a simple analytical expression for the forcing but the derivation of the corresponding analytical solutions is often laborious.

For the Stokes equations, a well-known set of solutions in the latter category is based on a forcing term in the form of a delta function, corresponding to an infinitely thin density anomaly at a certain depth. The solutions to these are also used in the propagator matrix method, where a convolution of delta-function solutions at different depths is used to obtain the response to arbitrary density anomalies (e.g. Hager and Richards, 1989). Solutions for the Cartesian case have been published (e.g. Zhong et al., 1993; Kramer et al., 2012). Spherical solutions can be found in, for example, Ribe (2009). They have previously been used for validation by, for example, Zhong et al. (2000, 2008), Choblet et al. (2007), Davies et al. (2013), and Burstedde et al. (2013). The derivation of these solutions is non-trivial, and often only cases with simpler free-slip boundary conditions are explicitly presented. Here, we derive solutions for both free- and zero-slip boundary conditions. In our convergence analyses, where we compare the analytical solutions with numerical solutions obtained using the Fluidity computational modelling framework (Davies et al., 2011), we discuss a limitation of this set of solutions for the benchmarking of geodynamic codes.

For this reason, we also present solutions based on a smooth forcing term, with radial dependence formed by a polynomial of arbitrary order, again for free- and zero-slip cases. This set of solutions provides a flexible way to test mantle dynamics codes with physically relevant solutions, where, for instance, a high-order polynomial forcing can be used to obtain solutions with a strong gradient near the surface. The radial dependence can be combined with spherical harmonics of arbitrary degree and order.

A key step in the derivation of these benchmarks is a decomposition of the solution into poloidal and toroidal components in the Mie representation (Backus, 1986). This results in a biharmonic equation for the poloidal scalar function. Combined with a set of conditions for the radial dependence of this poloidal function, analytical Stokes solutions can be obtained that satisfy desired free-slip or zero-slip conditions. Similar techniques have been used in Tosi and Martinec (2007) and Horbach et al. (2020). In Sect. 5.1, we will discuss previously published analytical benchmark cases in shell domains and how they differ from those presented here (e.g. Blinova et al., 2016; Thieulot, 2017; Horbach et al., 2020).

Finally, in addition to the delta-function and smooth cases with either free-slip or zero-slip boundary conditions in a spherical-shell domain, we present the solutions for the corresponding four cases in a 2-D cylindrical-shell domain (annulus). Although, ultimately, mantle convection is a 3-D phenomenon, a number of processes can be modelled adequately in two dimensions, and accordingly access to benchmark cases for 2-D numerical models is equally important. The number of published analytical Stokes solutions in 2-D cylindrical-shell domains, which are suitable as geodynamical benchmarks and include a complete derivation, is limited. By presenting this extensive set of explicit analytical solutions, we provide a suite of verification cases for use by the wider community of mantle dynamics code developers. An implementation of the solutions is provided through the Python package Assess (Analytical Solutions for the Stokes Equations in Spherical Shells; Kramer, 2020).

The remainder of this paper is structured as follows. In Sect. 2, we derive the analytical solutions in cylindrical (Sect. 2.2) and spherical (Sect. 2.3) geometries. Smooth solutions are first provided, for free-slip and zero-slip cases, followed by the delta-function solutions. In Sect. 3, we briefly describe the Fluidity computational modelling framework, which is used in Sect. 4 to obtain convergence results of the P2–P1 finite element discretisation to the analytical solutions. In Sect. 5, we discuss these results and relate the analytical solutions presented here to those that have previously been published. In particular, we discuss an issue with the delta-function cases for discretisations that use continuous pressure. To demonstrate, we also show results with a P2_bubble–P1_DG discretisation that overcomes this issue.

2 Analytical Solutions

2.1 Equations

The following derivation is applicable to the incompressible Stokes equations:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }-\mathrm{\nabla }\cdot \mathit{\tau }+\mathrm{\nabla }p=-g{\mathit{\rho }}^{\prime }\widehat{\mathit{r}},\text{(2)}& {\displaystyle }& {\displaystyle }\mathit{\tau }=\mathit{\nu }\left[\mathrm{\nabla }\mathit{u}+\mathrm{\nabla }{\mathit{u}}^T\right],\text{(3)}& {\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0},\end{array}$$

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 $-\widehat{\mathit{r}}$, and a dimensionless density deviation ρ^′: $\mathit{\rho }={\mathit{\rho }}_{\mathrm{0}}(\mathrm{1}+{\mathit{\rho }}^{\prime })$, where ρ is the density and the reference density ρ₀ is constant. The equations are solved in a 2-D cylindrical or 3-D spherical domain, bounded by $R_{-}\le r\le R_{+}$, 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 ψ:

$$\begin{array}{}\text{(4)}& u_r=-{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\phi }}},u_{\mathit{\phi }}={\displaystyle \frac{\partial \mathit{\psi }}{\partial r}},\end{array}$$

where φ is the angle from the x axis in polar coordinates, $x=r\cos \left(\mathit{\phi }\right),y=r\sin \left(\mathit{\phi }\right)$, and u_r and u_φ are the radial and transverse components of velocity, respectively. The normal deviatoric stress and shear stress components are given by

$$\begin{array}{}\text{(5)}& {\displaystyle }{\mathit{\tau }}_{rr}=\mathrm{2}\mathit{\nu }\widehat{\mathit{r}}\cdot \left[\mathrm{\nabla }\mathit{u}\right]\cdot \widehat{\mathit{r}}=-\mathrm{2}\mathit{\nu }{\displaystyle \frac{\partial }{\partial r}}\left({\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\phi }}}\right),\text{(6)}& {\displaystyle }\begin{array}{rl}{\mathit{\tau }}_{r\mathit{\phi }}& =\mathit{\nu }\widehat{\mathit{r}}\cdot \left[\mathrm{\nabla }\mathit{u}\right]\cdot \widehat{\mathit{\phi }}+\mathit{\nu }\widehat{\mathit{\phi }}\cdot \left[\mathrm{\nabla }\mathit{u}\right]\cdot \widehat{\mathit{r}}\\ & =\mathit{\nu }\left({\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial r^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial r}}-{\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial {\mathit{\phi }}^{\mathrm{2}}}}\right).\end{array}\end{array}$$

The momentum equation (Eq. 1) can be decomposed into radial and transverse components:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathit{\nu }}{r}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\psi }+{\displaystyle \frac{\partial p}{\partial r}}& {\displaystyle }=-g{\mathit{\rho }}^{\prime },\text{(8)}& {\displaystyle }-\mathit{\nu }{\displaystyle \frac{\partial }{\partial r}}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\psi }+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial p}{\partial \mathit{\phi }}}& {\displaystyle }=\mathrm{0},\end{array}$$

where the Laplacian, in polar coordinates, is given by

$$\begin{array}{}\text{(9)}& {\mathrm{\nabla }}^{\mathrm{2}}\mathit{\psi }={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\displaystyle \frac{\partial \mathit{\psi }}{\partial r}}\right)+{\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial {\mathit{\phi }}^{\mathrm{2}}}}.\end{array}$$

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

The curl of the momentum equation is obtained by summation of the operators $-\frac{\mathrm{1}}{r}\frac{\partial \cdot }{\partial \mathit{\phi }}$ and $\frac{\mathrm{1}}{r}\frac{\partial }{\partial r}\left(r\cdot \right)$ applied to the radial (Eq. 7) and transverse (Eq. 8) components, respectively, which leads to

$$\begin{array}{}\text{(10)}& -{\mathrm{\nabla }}^{\mathrm{4}}\mathit{\psi }={\displaystyle \frac{g}{\mathit{\nu }}}{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial {\mathit{\rho }}^{\prime }}{\partial \mathit{\phi }}}.\end{array}$$

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

$$\begin{array}{}\text{(11)}& \begin{array}{rl}\mathit{\psi }(r,\mathit{\phi })=\sum _{n>\mathrm{1}}& \left(A_n{r}^n+B_n{r}^{-n}+C_n{r}^{n+\mathrm{2}}+D_n{r}^{-n+\mathrm{2}}\right)\\ & \left(e_n\sin \left(n\mathit{\phi }\right)+f_n\cos \left(n\mathit{\phi }\right)\right)\\ +& \left(A_{\mathrm{1}}r+B_{\mathrm{1}}{r}^{-\mathrm{1}}+C_{\mathrm{1}}{r}^{\mathrm{3}}+D_{\mathrm{1}}r\ln r\right)\\ & \left(e_{\mathrm{1}}\sin \right(\mathit{\phi })+f_{\mathrm{1}}\cos (\mathit{\phi }\left)\right)\\ +& A_{\mathrm{0}}+B_{\mathrm{0}}\ln r+C_{\mathrm{0}}{r}^{\mathrm{2}}+D_{\mathrm{0}}{r}^{\mathrm{2}}\ln r,\end{array}\end{array}$$

where A_n, B_n, C_n, D_n, e_n, and f_n are constant coefficients, n≥0. The A_n and B_n terms are the standard solutions to the homogeneous harmonic equation, and the C_n and D_n terms are obtained as inhomogeneous solutions of the harmonic equation with the homogeneous harmonic solutions (A_n and B_n terms) as the right-hand side. The fact that these (the C_n and D_n terms) are homogeneous biharmonic solutions then follows from ∇⁴=∇²∇². In the following, we will, for simplicity, focus on sin (nφ) solutions for a single n>1 and set $e_n=\mathrm{1},f_n=\mathrm{0}$.

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

$$\begin{array}{}\text{(12)}& {\mathrm{\nabla }}^{\mathrm{2}}p=-{\displaystyle \frac{g}{r}}{\displaystyle \frac{\partial \left(r{\mathit{\rho }}^{\prime }\right)}{\partial r}}.\end{array}$$

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. ${\mathit{\rho }}^{\prime }=\mathrm{0}$, are thus the standard harmonic solutions (again neglecting the $n=\mathrm{0},\mathrm{1}$ solutions)

$$\begin{array}{}\text{(13)}& p(r,\mathit{\phi })=\sum _{n>\mathrm{1}}\left(G_n{r}^n+H_n{r}^{-n}\right)\cos \left(n\mathit{\phi }\right),\end{array}$$

where G_n and H_n are constant coefficients, n>1.

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

$$\begin{array}{}\text{(14)}& G_n=-\mathrm{4}\mathit{\nu }{C}_n(n+\mathrm{1}),H_n=-\mathrm{4}\mathit{\nu }{D}_n(n-\mathrm{1}),\end{array}$$

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:

$$\begin{array}{}\text{(15)}& {\mathit{\rho }}^{\prime }={\displaystyle \frac{r^k}{R_{+}^k}}\cos \left(n\mathit{\phi }\right),\end{array}$$

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

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& \mathit{\psi }=E{r}^{k+\mathrm{3}}\sin \left(n\mathit{\phi }\right),\\ & E={\displaystyle \frac{g{R}_{+}^{-k}n}{\mathit{\nu }\left((k+\mathrm{3}{)}^{\mathrm{2}}-n^{\mathrm{2}}\right)\left((k+\mathrm{1}{)}^{\mathrm{2}}-n^{\mathrm{2}}\right)}},\end{array}\end{array}$$

assuming $k\ne n-\mathrm{3}$ and $k\ne n-\mathrm{1}$, so that a general solution can be written as

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& \mathit{\psi }(r,\mathit{\phi })=\left(A{r}^n+B{r}^{-n}+C{r}^{n+\mathrm{2}}\right.\\ & \left.+D{r}^{-n+\mathrm{2}}+E{r}^{k+\mathrm{3}}\right)\sin \left(n\mathit{\phi }\right).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(18)}& p(r,\mathit{\phi })=F{r}^{k+\mathrm{1}}\cos \left(n\mathit{\phi }\right),F=-{\displaystyle \frac{g{R}_{+}^{-k}(k+\mathrm{1})}{(k+\mathrm{1}{)}^{\mathrm{2}}-n^{\mathrm{2}}}}.\end{array}$$

The general solution for pressure is thus

$$\begin{array}{}\text{(19)}& p(r,\mathit{\phi })=\left(G{r}^n+H{r}^{-n}+F{r}^{k+\mathrm{1}}\right)\cos \left(n\mathit{\phi }\right),\end{array}$$

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 $\frac{\partial \mathit{\psi }}{\partial \mathit{\phi }}=\mathrm{0}$. Two further equations are found by imposing either τ_rφ=0 (free slip) or $\frac{\partial \mathit{\psi }}{\partial }r=\mathrm{0}$ (zero slip) at both boundaries.

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

$$\begin{array}{ll}{\displaystyle }A& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{-n+\mathrm{3}}}{\mathit{\nu }}}{\displaystyle \frac{{\mathit{\alpha }}^{k+n+\mathrm{3}}-{\mathit{\alpha }}^{\mathrm{2}}}{\mathrm{4}\left(\mathit{\alpha }+{\mathit{\alpha }}^n\right)\left({\mathit{\alpha }}^n-\mathit{\alpha }\right)\left(k+n+\mathrm{1}\right)\left(k-n+\mathrm{3}\right)}}\\ {\displaystyle }B& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{n+\mathrm{3}}}{\mathit{\nu }}}{\displaystyle \frac{{\mathit{\alpha }}^{k+n+\mathrm{3}}-{\mathit{\alpha }}^{\mathrm{2}n+\mathrm{2}}}{\mathrm{4}\left({\mathit{\alpha }}^{n+\mathrm{1}}+\mathrm{1}\right)\left({\mathit{\alpha }}^{n+\mathrm{1}}-\mathrm{1}\right)\left(k+n+\mathrm{3}\right)\left(k-n+\mathrm{1}\right)}}\\ {\displaystyle }C& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{-n+\mathrm{1}}}{\mathit{\nu }}}{\displaystyle \frac{-{\mathit{\alpha }}^{k+n+\mathrm{3}}+\mathrm{1}}{\mathrm{4}\left({\mathit{\alpha }}^{n+\mathrm{1}}+\mathrm{1}\right)\left({\mathit{\alpha }}^{n+\mathrm{1}}-\mathrm{1}\right)\left(k+n+\mathrm{3}\right)\left(k-n+\mathrm{1}\right)}}\\ {\displaystyle }D& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{n+\mathrm{1}}}{\mathit{\nu }}}{\displaystyle \frac{-{\mathit{\alpha }}^{k+n+\mathrm{3}}+{\mathit{\alpha }}^{\mathrm{2}n}}{\mathrm{4}\left(\mathit{\alpha }+{\mathit{\alpha }}^n\right)\left({\mathit{\alpha }}^n-\mathit{\alpha }\right)\left(k+n+\mathrm{1}\right)\left(k-n+\mathrm{3}\right)}},\end{array}$$

where we use $\mathit{\alpha }=R_{-}/R_{+}$.

The zero-slip solution coefficients are given by

$$\begin{array}{ll}{\displaystyle }A& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{-n+\mathrm{3}}n}{\mathit{\nu }}}[\left({\mathit{\alpha }}^{k+n+\mathrm{3}}+{\mathit{\alpha }}^{\mathrm{2}n}\right)\left(k+n+\mathrm{1}\right)\left(n+\mathrm{1}\right)\\ {\displaystyle }& {\displaystyle }-\left({\mathit{\alpha }}^{k+n+\mathrm{1}}+{\mathit{\alpha }}^{\mathrm{2}n+\mathrm{2}}\right)\left(k+n+\mathrm{3}\right)n-\left({\mathit{\alpha }}^{k+\mathrm{3}n+\mathrm{3}}+\mathrm{1}\right)\\ {\displaystyle }& {\displaystyle }\left(k-n+\mathrm{1}\right)]/[\mathrm{2}\left({\left({\mathit{\alpha }}^{n+\mathrm{1}}-{\mathit{\alpha }}^{n-\mathrm{1}}\right)}^{\mathrm{2}}{n}^{\mathrm{2}}-{\left({\mathit{\alpha }}^{\mathrm{2}n}-\mathrm{1}\right)}^{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }\left({\left(k+\mathrm{3}\right)}^{\mathrm{2}}-n^{\mathrm{2}}\right)\left({\left(k+\mathrm{1}\right)}^{\mathrm{2}}-n^{\mathrm{2}}\right)]\\ {\displaystyle }B& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{n+\mathrm{3}}n}{\mathit{\nu }}}[-\left({\mathit{\alpha }}^{k+\mathrm{3}n+\mathrm{3}}+{\mathit{\alpha }}^{\mathrm{2}n}\right)\left(k-n+\mathrm{1}\right)\left(n-\mathrm{1}\right)\\ {\displaystyle }& {\displaystyle }+\left({\mathit{\alpha }}^{k+\mathrm{3}n+\mathrm{1}}+{\mathit{\alpha }}^{\mathrm{2}n+\mathrm{2}}\right)\left(k-n+\mathrm{3}\right)n-\left({\mathit{\alpha }}^{k+n+\mathrm{3}}+{\mathit{\alpha }}^{\mathrm{4}n}\right)\\ {\displaystyle }& {\displaystyle }\left(k+n+\mathrm{1}\right)]/[\mathrm{2}\left({\left({\mathit{\alpha }}^{n+\mathrm{1}}-{\mathit{\alpha }}^{n-\mathrm{1}}\right)}^{\mathrm{2}}{n}^{\mathrm{2}}-{\left({\mathit{\alpha }}^{\mathrm{2}n}-\mathrm{1}\right)}^{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }\left({\left(k+\mathrm{3}\right)}^{\mathrm{2}}-n^{\mathrm{2}}\right)\left({\left(k+\mathrm{1}\right)}^{\mathrm{2}}-n^{\mathrm{2}}\right)]\\ {\displaystyle }C& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{-n+\mathrm{1}}n}{\mathit{\nu }}}[\left({\mathit{\alpha }}^{k+n+\mathrm{1}}+{\mathit{\alpha }}^{\mathrm{2}n}\right)\left(k+n+\mathrm{3}\right)\left(n-\mathrm{1}\right)\\ {\displaystyle }& {\displaystyle }-\left({\mathit{\alpha }}^{k+n+\mathrm{3}}+{\mathit{\alpha }}^{\mathrm{2}n-\mathrm{2}}\right)\left(k+n+\mathrm{1}\right)n+\left({\mathit{\alpha }}^{k+\mathrm{3}n+\mathrm{1}}+\mathrm{1}\right)\\ {\displaystyle }& {\displaystyle }\left(k-n+\mathrm{3}\right)]/[\mathrm{2}\left({\left({\mathit{\alpha }}^{n+\mathrm{1}}-{\mathit{\alpha }}^{n-\mathrm{1}}\right)}^{\mathrm{2}}{n}^{\mathrm{2}}-{\left({\mathit{\alpha }}^{\mathrm{2}n}-\mathrm{1}\right)}^{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }\left({\left(k+\mathrm{3}\right)}^{\mathrm{2}}-n^{\mathrm{2}}\right)\left({\left(k+\mathrm{1}\right)}^{\mathrm{2}}-n^{\mathrm{2}}\right)]\\ {\displaystyle }D& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{n+\mathrm{1}}n}{\mathit{\nu }}}[-\left({\mathit{\alpha }}^{k+\mathrm{3}n+\mathrm{1}}+{\mathit{\alpha }}^{\mathrm{2}n}\right)\left(k-n+\mathrm{3}\right)\left(n+\mathrm{1}\right)\\ {\displaystyle }& {\displaystyle }+\left({\mathit{\alpha }}^{k+\mathrm{3}n+\mathrm{3}}+{\mathit{\alpha }}^{\mathrm{2}n-\mathrm{2}}\right)\left(k-n+\mathrm{1}\right)n+\left({\mathit{\alpha }}^{k+n+\mathrm{1}}+{\mathit{\alpha }}^{\mathrm{4}n}\right)\\ {\displaystyle }& {\displaystyle }\left(k+n+\mathrm{3}\right)]/[\mathrm{2}\left({\left({\mathit{\alpha }}^{n+\mathrm{1}}-{\mathit{\alpha }}^{n-\mathrm{1}}\right)}^{\mathrm{2}}{n}^{\mathrm{2}}-{\left({\mathit{\alpha }}^{\mathrm{2}n}-\mathrm{1}\right)}^{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }\left({\left(k+\mathrm{3}\right)}^{\mathrm{2}}-n^{\mathrm{2}}\right)\left({\left(k+\mathrm{1}\right)}^{\mathrm{2}}-n^{\mathrm{2}}\right)].\end{array}$$

2.2.2 Green's function solution – cylindrical

Another set of useful solutions is found considering the following perturbation density,

$$\begin{array}{}\text{(20)}& {\mathit{\rho }}^{\prime }=\mathit{\delta }(r-r^{\prime })\cos \left(n\mathit{\phi }\right),\end{array}$$

representing an infinitely thin density anomaly at $r=r^{\prime }$ with $R_{-}<r^{\prime }<R_{+}$. Since ${\mathit{\rho }}^{\prime }=\mathrm{0}$ for $r^{\prime }\ne r$, the solution is described by combining two homogeneous solutions

$$\begin{array}{}\text{(21)}& \begin{array}{rl}& \mathit{\psi }(r,\mathit{\phi })=\\ & \left\{\begin{array}{l}{\mathit{\psi }}_{-}(r,\mathit{\phi })=\left(A_{-}{r}^n+B_{-}{r}^{-n}+C_{-}{r}^{n+\mathrm{2}}\right.\\ \left.+D_{-}{r}^{-n+\mathrm{2}}\right)\sin \left(n\mathit{\phi }\right)\text{ for }{R}_{-}\le r<r^{\prime },\\ {\mathit{\psi }}_{+}(r,\mathit{\phi })=\left(A_{+}{r}^n+B_{+}{r}^{-n}+C_{+}{r}^{n+\mathrm{2}}\right.\\ \left.+D_{+}{r}^{-n+\mathrm{2}}\right)\sin \left(n\mathit{\phi }\right)\text{ for }{r}^{\prime }<r\le R_{+}.\end{array}\right.\end{array}\end{array}$$

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:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& {\mathit{\psi }}_{-}(r^{\prime },\mathit{\phi })={\mathit{\psi }}_{+}(r^{\prime },\mathit{\phi }),\text{ and }\\ & {\displaystyle \frac{\partial {\mathit{\psi }}_{-}}{\partial r}}(r^{\prime },\mathit{\phi })={\displaystyle \frac{\partial {\mathit{\psi }}_{+}}{\partial r}}(r^{\prime },\mathit{\phi }).\end{array}\end{array}$$

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

$$\begin{array}{}\text{(23)}& {\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\psi }}_{-}}{\partial r^{\mathrm{2}}}}(r^{\prime },\mathit{\phi })={\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\psi }}_{+}}{\partial r^{\mathrm{2}}}}(r^{\prime },\mathit{\phi }).\end{array}$$

Finally, the eighth equation is obtained by integrating Eq. (10) over a small strip $r^{\prime }-\mathit{ϵ}\le r\le r^{\prime }+\mathit{ϵ}$, between two arbitrary angles ${\mathit{\phi }}_{\mathrm{1}}\le \mathit{\phi }\le {\mathit{\phi }}_{\mathrm{2}}$:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{r=r^{\prime }-\mathit{ϵ}}{\overset{r^{\prime }+\mathit{ϵ}}{\int }}\underset{\mathit{\phi }={\mathit{\phi }}_{\mathrm{1}}}{\overset{{\mathit{\phi }}_{\mathrm{2}}}{\int }}\left({\mathrm{\nabla }}^{\mathrm{4}}\mathit{\psi }\right)r\mathrm{d}r\mathrm{d}\mathit{\phi }\\ \text{(24)}& {\displaystyle }& {\displaystyle }=\underset{r=r^{\prime }-\mathit{ϵ}}{\overset{r^{\prime }+\mathit{ϵ}}{\int }}\underset{\mathit{\phi }={\mathit{\phi }}_{\mathrm{1}}}{\overset{{\mathit{\phi }}_{\mathrm{2}}}{\int }}{\displaystyle \frac{gn\sin \left(n\mathit{\phi }\right)}{\mathit{\nu }r}}\mathit{\delta }(r-r^{\prime })r\mathrm{d}r\mathrm{d}\mathit{\phi },{\displaystyle }& {\displaystyle }{\left[\underset{\mathit{\phi }={\mathit{\phi }}_{\mathrm{1}}}{\overset{{\mathit{\phi }}_{\mathrm{2}}}{\int }}\widehat{\mathit{r}}\cdot \mathrm{\nabla }\left({\mathrm{\nabla }}^{\mathrm{2}}\mathit{\psi }\right)rd\mathit{\phi }\right]}_{r=r^{\prime }-\mathit{ϵ}}^{r=r^{\prime }+\mathit{ϵ}}\\ {\displaystyle }& {\displaystyle }+{\left[\underset{r=r^{\prime }-\mathit{ϵ}}{\overset{r^{\prime }+\mathit{ϵ}}{\int }}\widehat{\mathit{\phi }}\cdot \mathrm{\nabla }\left({\mathrm{\nabla }}^{\mathrm{2}}\mathit{\psi }\right)r\mathrm{d}r\right]}_{\mathit{\phi }={\mathit{\phi }}_{\mathrm{1}}}^{{\mathit{\phi }}_{\mathrm{2}}}\\ \text{(25)}& {\displaystyle }& {\displaystyle }=\underset{\mathit{\phi }={\mathit{\phi }}_{\mathrm{1}}}{\overset{{\mathit{\phi }}_{\mathrm{2}}}{\int }}{\displaystyle \frac{gn\sin \left(n\mathit{\phi }\right)}{\mathit{\nu }}}\mathrm{d}\mathit{\phi }.\end{array}$$

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

$$\begin{array}{rl}& \underset{\mathit{\phi }={\mathit{\phi }}_{\mathrm{1}}}{\overset{{\mathit{\phi }}_{\mathrm{2}}}{\int }}{\displaystyle \frac{\partial }{\partial r}}{\left.\left({\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\displaystyle \frac{\partial {\mathit{\psi }}_{+}}{\partial r}}\right)+{\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\psi }}_{+}}{\partial {\mathit{\phi }}^{\mathrm{2}}}}\right)\right|}_{r=r^{\prime }}{r}^{\prime }\mathrm{d}\mathit{\phi }\\ & -\underset{\mathit{\phi }={\mathit{\phi }}_{\mathrm{1}}}{\overset{{\mathit{\phi }}_{\mathrm{2}}}{\int }}{\displaystyle \frac{\partial }{\partial r}}{\left.\left({\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial r}}\left(r{\displaystyle \frac{\partial {\mathit{\psi }}_{-}}{\partial r}}\right)+{\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\psi }}_{-}}{\partial {\mathit{\phi }}^{\mathrm{2}}}}\right)\right|}_{r=r^{\prime }}{r}^{\prime }\mathrm{d}\mathit{\phi }.\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{\mathit{\phi }={\mathit{\phi }}_{\mathrm{1}}}{\overset{{\mathit{\phi }}_{\mathrm{2}}}{\int }}\left({\displaystyle \frac{{\partial }^{\mathrm{3}}{\mathit{\psi }}_{+}}{\partial r^{\mathrm{3}}}}(r^{\prime },\mathit{\phi })-{\displaystyle \frac{{\partial }^{\mathrm{3}}{\mathit{\psi }}_{-}}{\partial r^{\mathrm{3}}}}(r^{\prime },\mathit{\phi })\right)r^{\prime }\mathrm{d}\mathit{\phi }=\\ {\displaystyle }& {\displaystyle }\underset{\mathit{\phi }={\mathit{\phi }}_{\mathrm{1}}}{\overset{{\mathit{\phi }}_{\mathrm{2}}}{\int }}{\displaystyle \frac{gn\sin \left(n\mathit{\phi }\right)}{\mathit{\nu }}}\mathrm{d}\mathit{\phi }.\end{array}$$

Thus, for this to hold for arbitrary φ₁ and φ₂, we need

$$\begin{array}{}\text{(26)}& {\displaystyle \frac{{\partial }^{\mathrm{3}}{\mathit{\psi }}_{+}}{\partial r^{\mathrm{3}}}}(r^{\prime },\mathit{\phi })-{\displaystyle \frac{{\partial }^{\mathrm{3}}{\mathit{\psi }}_{-}}{\partial r^{\mathrm{3}}}}(r^{\prime },\mathit{\phi })={\displaystyle \frac{gn\sin \left(n\mathit{\phi }\right)}{\mathit{\nu }{r}^{\prime }}}.\end{array}$$

The solution coefficients for free-slip boundary conditions at $r=R_{-}$ and $r=R_{+}$ are

$$\begin{array}{ll}{\displaystyle }{A}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^{-n+\mathrm{2}}}{\mathit{\nu }}}{\displaystyle \frac{\pm \left({\mathit{\alpha }}_{\mp }^{\mathrm{2}n-\mathrm{2}}-\mathrm{1}\right)}{\mathrm{8}\left({\mathit{\alpha }}_{\pm }^{\mathrm{2}n-\mathrm{2}}-{\mathit{\alpha }}_{\mp }^{\mathrm{2}n-\mathrm{2}}\right)\left(n-\mathrm{1}\right)}}\\ {\displaystyle }{B}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^{n+\mathrm{2}}}{\mathit{\nu }}}{\displaystyle \frac{\pm \left({\mathit{\alpha }}_{\mp }^{\mathrm{2}n+\mathrm{2}}-\mathrm{1}\right){\mathit{\alpha }}_{\pm }^{\mathrm{2}n+\mathrm{2}}}{\mathrm{8}\left({\mathit{\alpha }}_{\pm }^{\mathrm{2}n+\mathrm{2}}-{\mathit{\alpha }}_{\mp }^{\mathrm{2}n+\mathrm{2}}\right)\left(n+\mathrm{1}\right)}}\\ {\displaystyle }{C}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^{-n}}{\mathit{\nu }}}{\displaystyle \frac{\pm \left({\mathit{\alpha }}_{\mp }^{\mathrm{2}n+\mathrm{2}}-\mathrm{1}\right)}{\mathrm{8}\left({\mathit{\alpha }}_{\mp }^{\mathrm{2}n+\mathrm{2}}-{\mathit{\alpha }}_{\pm }^{\mathrm{2}n+\mathrm{2}}\right)\left(n+\mathrm{1}\right)}}\\ {\displaystyle }{D}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^n}{\mathit{\nu }}}{\displaystyle \frac{\pm \left({\mathit{\alpha }}_{\mp }^{\mathrm{2}n-\mathrm{2}}-\mathrm{1}\right){\mathit{\alpha }}_{\pm }^{\mathrm{2}n-\mathrm{2}}}{\mathrm{8}\left({\mathit{\alpha }}_{\mp }^{\mathrm{2}n-\mathrm{2}}-{\mathit{\alpha }}_{\pm }^{\mathrm{2}n-\mathrm{2}}\right)\left(n-\mathrm{1}\right)}},\end{array}$$

where we use ${\mathit{\alpha }}_{\pm }=R_{\pm }/r^{\prime }$.

Zero-slip conditions at both boundaries lead to

$$\begin{array}{ll}{\displaystyle }{A}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^{-n+\mathrm{2}}}{\mathit{\nu }}}\left[\left(\left({\mathit{\alpha }}_{+}^{\mathrm{2}}-{\mathit{\alpha }}_{-}^{\mathrm{2}}\right)n+{\mathit{\alpha }}_{+}^{-\mathrm{2}n}-{\mathit{\alpha }}_{-}^{-\mathrm{2}n}\pm (n+\mathrm{1})\right)\right.\\ {\displaystyle }& {\displaystyle }\left.\left(n-\mathrm{1}\right)-\left({\mathit{\alpha }}_{+}^{-\mathrm{2}n}{\mathit{\alpha }}_{-}^{\mathrm{2}}-{\mathit{\alpha }}_{+}^{\mathrm{2}}{\mathit{\alpha }}_{-}^{-\mathrm{2}n}\right)n\pm \left({\mathit{\gamma }}^{\pm \mathrm{2}n}-{\mathit{\gamma }}^{\mp \mathrm{2}}{n}^{\mathrm{2}}\right)\right]\\ {\displaystyle }& {\displaystyle }/\left[\mathrm{8}\left({\left(\mathit{\gamma }-{\mathit{\gamma }}^{-\mathrm{1}}\right)}^{\mathrm{2}}{n}^{\mathrm{2}}-{\left({\mathit{\gamma }}^{-n}-{\mathit{\gamma }}^n\right)}^{\mathrm{2}}\right)\left(n-\mathrm{1}\right)\right]\\ {\displaystyle }{B}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^{n+\mathrm{2}}}{\mathit{\nu }}}\left[\left(\left({\mathit{\alpha }}_{+}^{\mathrm{2}}-{\mathit{\alpha }}_{-}^{\mathrm{2}}\right)n-{\mathit{\alpha }}_{+}^{\mathrm{2}n}+{\mathit{\alpha }}_{-}^{\mathrm{2}n}\pm (n-\mathrm{1})\right)\left(n+\mathrm{1}\right)\right.\\ {\displaystyle }& {\displaystyle }\left.+\left({\mathit{\alpha }}_{+}^{\mathrm{2}n}{\mathit{\alpha }}_{-}^{\mathrm{2}}-{\mathit{\alpha }}_{+}^{\mathrm{2}}{\mathit{\alpha }}_{-}^{\mathrm{2}n}\right)n\pm \left({\mathit{\gamma }}^{\mp \mathrm{2}n}-{\mathit{\gamma }}^{\mp \mathrm{2}}{n}^{\mathrm{2}}\right)\right]\\ {\displaystyle }& {\displaystyle }/\left[\mathrm{8}\left({\left(\mathit{\gamma }-{\mathit{\gamma }}^{-\mathrm{1}}\right)}^{\mathrm{2}}{n}^{\mathrm{2}}-{\left({\mathit{\gamma }}^{-n}-{\mathit{\gamma }}^n\right)}^{\mathrm{2}}\right)\left(n+\mathrm{1}\right)\right]\\ {\displaystyle }{C}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^{-n}}{\mathit{\nu }}}\left[\left(\left({\mathit{\alpha }}_{-}^{-\mathrm{2}}-{\mathit{\alpha }}_{+}^{-\mathrm{2}}\right)n+{\mathit{\alpha }}_{+}^{-\mathrm{2}n}-{\mathit{\alpha }}_{-}^{-\mathrm{2}n}\mp (n-\mathrm{1})\right)\right.\\ {\displaystyle }& {\displaystyle }\left.\left(n+\mathrm{1}\right)-\left({\mathit{\alpha }}_{+}^{-\mathrm{2}n}{\mathit{\alpha }}_{-}^{-\mathrm{2}}-{\mathit{\alpha }}_{+}^{-\mathrm{2}}{\mathit{\alpha }}_{-}^{-\mathrm{2}n}\right)n\mp \left({\mathit{\gamma }}^{\pm \mathrm{2}n}-{\mathit{\gamma }}^{\pm \mathrm{2}}{n}^{\mathrm{2}}\right)\right]\\ {\displaystyle }& {\displaystyle }/\left[\mathrm{8}\left({\left(\mathit{\gamma }-{\mathit{\gamma }}^{-\mathrm{1}}\right)}^{\mathrm{2}}{n}^{\mathrm{2}}-{\left({\mathit{\gamma }}^{-n}-{\mathit{\gamma }}^n\right)}^{\mathrm{2}}\right)\left(n+\mathrm{1}\right)\right]\\ {\displaystyle }{D}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^n}{\mathit{\nu }}}\left[\left(\left({\mathit{\alpha }}_{-}^{-\mathrm{2}}-{\mathit{\alpha }}_{+}^{-\mathrm{2}}\right)n-{\mathit{\alpha }}_{+}^{\mathrm{2}n}+{\mathit{\alpha }}_{-}^{\mathrm{2}n}\mp (n+\mathrm{1})\right)\left(n-\mathrm{1}\right)\right.\\ {\displaystyle }& {\displaystyle }\left.+\left({\mathit{\alpha }}_{+}^{\mathrm{2}n}{\mathit{\alpha }}_{-}^{-\mathrm{2}}-{\mathit{\alpha }}_{+}^{-\mathrm{2}}{\mathit{\alpha }}_{-}^{\mathrm{2}n}\right)n\mp \left({\mathit{\gamma }}^{\mp \mathrm{2}n}-{\mathit{\gamma }}^{\pm \mathrm{2}}{n}^{\mathrm{2}}\right)\right]\\ {\displaystyle }& {\displaystyle }/\left[\mathrm{8}\left({\left(\mathit{\gamma }-{\mathit{\gamma }}^{-\mathrm{1}}\right)}^{\mathrm{2}}{n}^{\mathrm{2}}-{\left({\mathit{\gamma }}^{-n}-{\mathit{\gamma }}^n\right)}^{\mathrm{2}}\right)\left(n-\mathrm{1}\right)\right],\end{array}$$

where, in addition, we use $\mathit{\gamma }=R_{-}/R_{+}$.

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 outer boundary.

2.3 Spherical

In this section, we derive the equivalent of the four cylindrical cases in a 3-D spherical domain, $R_{-}\le r\le R_{+}$. 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

$$\begin{array}{}\text{(27)}& \mathit{u}=\mathrm{\nabla }\times \left(\mathit{r}\times \mathrm{\nabla }\mathcal{P}\right)+\mathit{r}\times \mathrm{\nabla }\mathcal{T}\end{array}$$

using poloidal and toroidal scalar fields 𝒫 and 𝒯 (Backus, 1986), where $\mathit{r}=r\widehat{\mathit{r}}$. In spherical coordinates,

$$\begin{array}{}\text{(28)}& \begin{array}{rl}& u_r={\displaystyle \frac{\mathrm{1}}{r}}{\mathrm{\Lambda }}^{\mathrm{2}}\mathcal{P},u_{\mathit{\theta }}=-{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\left(r\mathcal{P}\right)}{\partial r\partial \mathit{\theta }}}-{\displaystyle \frac{\mathrm{1}}{\sin \left(\mathit{\theta }\right)}}{\displaystyle \frac{\partial \mathcal{T}}{\partial \mathit{\phi }}},\\ & u_{\mathit{\phi }}=-{\displaystyle \frac{\mathrm{1}}{r\sin \left(\mathit{\theta }\right)}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\left(r\mathcal{P}\right)}{\partial r\partial \mathit{\phi }}}+{\displaystyle \frac{\partial \mathcal{T}}{\partial \mathit{\theta }}},\end{array}\end{array}$$

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

$$\begin{array}{}\text{(29)}& \begin{array}{rl}& {\mathrm{\Lambda }}^{\mathrm{2}}\mathcal{P}={\displaystyle \frac{\mathrm{1}}{\sin \left(\mathit{\theta }\right)}}{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left(\sin \left(\mathit{\theta }\right){\displaystyle \frac{\partial \mathcal{P}}{\partial \mathit{\theta }}}\right)+{\displaystyle \frac{\mathrm{1}}{\sin (\mathit{\theta }{)}^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathcal{P}}{\partial {\mathit{\phi }}^{\mathrm{2}}}},\\ & \text{ so that }{\mathrm{\nabla }}^{\mathrm{2}}={\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}\left({\displaystyle \frac{\partial }{\partial r}}{r}^{\mathrm{2}}{\displaystyle \frac{\partial }{\partial r}}+{\mathrm{\Lambda }}^{\mathrm{2}}\right).\end{array}\end{array}$$

We further derive the components of stress normal to a spherical surface:

$$\begin{array}{ll}\text{(30)}& {\displaystyle }{\mathit{\tau }}_{rr}& {\displaystyle }=\mathrm{2}\mathit{\nu }{\displaystyle \frac{\partial u_r}{\partial r}}=\mathrm{2}\mathit{\nu }{\displaystyle \frac{\partial }{\partial r}}\left({\displaystyle \frac{\mathrm{1}}{r}}{\mathrm{\Lambda }}^{\mathrm{2}}\mathcal{P}\right),{\displaystyle }{\mathit{\tau }}_{r\mathit{\theta }}& {\displaystyle }=\mathit{\nu }r{\displaystyle \frac{\partial }{\partial r}}\left({\displaystyle \frac{u_{\mathit{\theta }}}{r}}\right)+{\displaystyle \frac{\mathit{\nu }}{r}}{\displaystyle \frac{\partial u_r}{\partial \mathit{\theta }}}=\\ {\displaystyle }& {\displaystyle }\mathit{\nu }{\displaystyle \frac{\partial }{\partial \mathit{\theta }}}\left({\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{\mathrm{\Lambda }}^{\mathrm{2}}\mathcal{P}-{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathcal{P}}{\partial r^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{2}}{r^{\mathrm{2}}}}\mathcal{P}\right)-\\ \text{(31)}& {\displaystyle }& {\displaystyle }\mathit{\nu }r{\displaystyle \frac{\partial }{\partial r}}\left({\displaystyle \frac{\mathrm{1}}{r\sin \left(\mathit{\theta }\right)}}{\displaystyle \frac{\partial \mathcal{T}}{\partial \mathit{\phi }}}\right),{\displaystyle }{\mathit{\tau }}_{r\mathit{\phi }}& {\displaystyle }=\mathit{\nu }r{\displaystyle \frac{\partial }{\partial r}}\left({\displaystyle \frac{u_{\mathit{\phi }}}{r}}\right)+{\displaystyle \frac{\mathit{\nu }}{r\sin \left(\mathit{\theta }\right)}}{\displaystyle \frac{\partial u_r}{\partial \mathit{\phi }}}=\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{\mathit{\nu }}{\sin \left(\mathit{\theta }\right)}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}\left({\displaystyle \frac{\mathrm{1}}{r^{\mathrm{2}}}}{\mathrm{\Lambda }}^{\mathrm{2}}\mathcal{P}-{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathcal{P}}{\partial r^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{2}}{r^{\mathrm{2}}}}\mathcal{P}\right)\\ \text{(32)}& {\displaystyle }& {\displaystyle }+\mathit{\nu }r{\displaystyle \frac{\partial }{\partial r}}\left({\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \mathcal{T}}{\partial \mathit{\theta }}}\right).\end{array}$$

Using these, we can work out spherical components of the momentum equation (Eq. 1) to be

$$\begin{array}{}\text{(33)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\mathit{\nu }}{r}}{\mathrm{\Lambda }}^{\mathrm{2}}{\mathrm{\nabla }}^{\mathrm{2}}\mathcal{P}+{\displaystyle \frac{\partial p}{\partial r}}=-g{\mathit{\rho }}^{\prime },\text{(34)}& {\displaystyle }& {\displaystyle \frac{\mathit{\nu }}{r}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\left(r{\mathrm{\nabla }}^{\mathrm{2}}\mathcal{P}\right)}{\partial r\partial \mathit{\theta }}}+{\displaystyle \frac{\mathit{\nu }}{\sin \left(\mathit{\theta }\right)}}{\displaystyle \frac{\partial {\mathrm{\nabla }}^{\mathrm{2}}\mathcal{T}}{\partial \mathit{\phi }}}+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial p}{\partial \mathit{\theta }}}=\mathrm{0},\text{(35)}& {\displaystyle }& {\displaystyle \frac{\mathit{\nu }}{r\sin \left(\mathit{\theta }\right)}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\left(r{\mathrm{\nabla }}^{\mathrm{2}}\mathcal{P}\right)}{\partial r\partial \mathit{\phi }}}+\mathit{\nu }{\displaystyle \frac{\partial {\mathrm{\nabla }}^{\mathrm{2}}\mathcal{T}}{\partial \mathit{\theta }}}+{\displaystyle \frac{\mathrm{1}}{r\sin \left(\mathit{\theta }\right)}}{\displaystyle \frac{\partial p}{\partial \mathit{\phi }}}=\mathrm{0}.\end{array}$$

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,

$$\begin{array}{}\text{(36)}& \begin{array}{rl}& -\widehat{\mathit{\theta }}{\displaystyle \frac{\mathit{\nu }}{\sin \left(\mathit{\theta }\right)}}{\displaystyle \frac{\partial {\mathrm{\nabla }}^{\mathrm{4}}\mathcal{P}}{\partial \mathit{\phi }}}+\widehat{\mathit{\phi }}\mathit{\nu }{\displaystyle \frac{\partial {\mathrm{\nabla }}^{\mathrm{4}}\mathcal{P}}{\partial \mathit{\theta }}}=-\widehat{\mathit{\theta }}{\displaystyle \frac{\mathrm{1}}{r\sin \left(\mathit{\theta }\right)}}{\displaystyle \frac{\partial g{\mathit{\rho }}^{\prime }}{\partial \mathit{\phi }}}\\ & +\widehat{\mathit{\phi }}{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial g{\mathit{\rho }}^{\prime }}{\partial \mathit{\theta }}}.\end{array}\end{array}$$

In a spherical domain, this implies that $\mathit{\nu }{\mathrm{\nabla }}^{\mathrm{4}}\mathcal{P}-g{\mathit{\rho }}^{\prime }/r$ varies in the radial direction only. Without loss of generality, we may therefore seek solutions 𝒫 that satisfy

$$\begin{array}{}\text{(37)}& {\mathrm{\nabla }}^{\mathrm{4}}\mathcal{P}={\displaystyle \frac{g{\mathit{\rho }}^{\prime }}{\mathit{\nu }r}},\end{array}$$

since for any other solution $\mathcal{P}+{\mathcal{P}}^{\prime }$ of Eq. (36) we know that ∇⁴𝒫^′ 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 ∇⁴𝒫=0 in 3-D are given by

$$\begin{array}{}\text{(38)}& \begin{array}{rl}& \mathcal{P}(r,\mathit{\theta },\mathit{\phi })=\sum _{l\ge \mathrm{0}}\sum _{m=-l}^l\left(A_{lm}{r}^l+B_{lm}{r}^{-l-\mathrm{1}}\right.\\ & \left.+C_{lm}{r}^{l+\mathrm{2}}+D_{lm}{r}^{-l+\mathrm{1}}\right)Y_{lm}(\mathit{\theta },\mathit{\phi }),\end{array}\end{array}$$

where A_lm, B_lm, C_lm, and D_lm are constant coefficients, and Y_lm indicates Laplace's spherical harmonic functions of degree l and order m, which have the property that

$$\begin{array}{}\text{(39)}& {\mathrm{\Lambda }}^{\mathrm{2}}{Y}_{lm}=-l(l+\mathrm{1})Y_{lm}.\end{array}$$

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

$$\begin{array}{}\text{(40)}& {\mathrm{\nabla }}^{\mathrm{2}}p=-{\displaystyle \frac{g}{r^{\mathrm{2}}}}{\displaystyle \frac{\partial }{\partial r}}\left(r^{\mathrm{2}}{\mathit{\rho }}^{\prime }\right).\end{array}$$

Homogeneous solutions for pressure are written as

$$\begin{array}{}\text{(41)}& p(r,\mathit{\theta },\mathit{\phi })=\sum _{l\ge \mathrm{0}}\sum _{m=-l}^l\left(G_{lm}{r}^l+H_{lm}{r}^{-l-\mathrm{1}}\right)Y_{lm}(\mathit{\theta },\mathit{\phi }),\end{array}$$

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

$$\begin{array}{}\text{(42)}& \begin{array}{rl}{G}_{lm}=& -\mathrm{2}\mathit{\nu }(l+\mathrm{1})(\mathrm{2}l+\mathrm{3})C_{lm},\\ H_{lm}=& -\mathrm{2}\mathit{\nu }l(\mathrm{2}l-\mathrm{1})D_{lm}.\end{array}\end{array}$$

Figure 2 Illustrations 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:

$$\begin{array}{}\text{(43)}& {\mathit{\rho }}^{\prime }={\displaystyle \frac{r^k}{R_{+}^k}}{Y}_{lm}(\mathit{\theta },\mathit{\phi }).\end{array}$$

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

$$\begin{array}{}\text{(44)}& \begin{array}{rl}& \mathcal{P}=E{r}^{k+\mathrm{3}}{Y}_{lm},\\ & E=\left[g{R}_{+}^{-k}\right]/\left[\mathit{\nu }\left((k+\mathrm{1})(k+\mathrm{2})-l(l+\mathrm{1})\right)\right.\\ & \left.\left((k+\mathrm{3})(k+\mathrm{4})-l(l+\mathrm{1})\right)\right],\end{array}\end{array}$$

with a more generic solution written as

$$\begin{array}{}\text{(45)}& \mathcal{P}=\left(A{r}^l+B{r}^{-l-\mathrm{1}}+C{r}^{l+\mathrm{2}}+D{r}^{-l+\mathrm{1}}+E{r}^{k+\mathrm{3}}\right)Y_{lm},\end{array}$$

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

$$\begin{array}{}\text{(46)}& \begin{array}{rl}p=& \left(G{r}^l+H{r}^{-l-\mathrm{1}}+F{r}^{k+\mathrm{1}}\right)Y_{lm},\\ F=& -{\displaystyle \frac{g(k+\mathrm{2})R_{+}^{-k}}{(k+\mathrm{1})(k+\mathrm{2})-l(l+\mathrm{1})}},\end{array}\end{array}$$

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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{no normal flow: }\\ \text{(47)}& {\displaystyle }& {\displaystyle }{u}_r={\displaystyle \frac{\mathrm{1}}{r}}{\mathrm{\Lambda }}^{\mathrm{2}}\mathcal{P}=\mathrm{0}⟹\mathcal{P}=\mathrm{0},& {\displaystyle }\text{at }r=R_{\pm },\end{array}$$

where we use ${\mathrm{\Lambda }}^{\mathrm{2}}\mathcal{P}=-l(l+\mathrm{1})\mathcal{P}$. The two no-normal-flow conditions at $r=R_{-}$ and $r=R_{+}$ are combined with two further conditions:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{ free slip: }{\mathit{\tau }}_{r\mathit{\theta }}={\mathit{\tau }}_{r\mathit{\phi }}=\mathrm{0}⟹{\displaystyle \frac{\mathrm{1}}{r}}{\mathrm{\Lambda }}^{\mathrm{2}}\mathcal{P}-{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathcal{P}}{\partial r^{\mathrm{2}}}}\\ \text{(48)}& {\displaystyle }& {\displaystyle }+{\displaystyle \frac{\mathrm{2}}{r^{\mathrm{2}}}}\mathcal{P}=\mathrm{0}⟹{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathcal{P}}{\partial r^{\mathrm{2}}}}=\mathrm{0},\text{at }r=R_{\pm },\text{ or }{\displaystyle }& {\displaystyle }\text{ zero slip: }{u}_{\mathit{\theta }}=u_{\mathit{\phi }}=\mathrm{0}⟹{\displaystyle \frac{\partial \left(r\mathcal{P}\right)}{\partial r}}=\\ \text{(49)}& {\displaystyle }& {\displaystyle }\mathrm{0}⟹{\displaystyle \frac{\partial \mathcal{P}}{\partial r}}=\mathrm{0},\text{at }r=R_{\pm }.\end{array}$$

For free-slip conditions at both $r=R_{-}$ and $r=R_{+}$, the solution coefficients are given by

$$\begin{array}{ll}{\displaystyle }A& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{-l+\mathrm{3}}\left({\mathit{\alpha }}^{k+\mathrm{3}}-{\mathit{\alpha }}^{-l+\mathrm{1}}\right)}{\mathrm{2}\mathit{\nu }\left({\mathit{\alpha }}^l-{\mathit{\alpha }}^{-l+\mathrm{1}}\right)\left(k+l+\mathrm{2}\right)\left(k-l+\mathrm{3}\right)\left(\mathrm{2}l+\mathrm{1}\right)}}\\ {\displaystyle }B& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{l+\mathrm{4}}-{\mathit{\alpha }}^{k+\mathrm{4}}+{\mathit{\alpha }}^{l+\mathrm{3}}}{\mathrm{2}\mathit{\nu }\left({\mathit{\alpha }}^{-l}-{\mathit{\alpha }}^{l+\mathrm{3}}\right)\left(k+l+\mathrm{4}\right)\left(k-l+\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{1}\right)}}\\ {\displaystyle }C& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{-l+\mathrm{1}}{\mathit{\alpha }}^{k+\mathrm{4}}-{\mathit{\alpha }}^{-l}}{\mathrm{2}\mathit{\nu }\left({\mathit{\alpha }}^{-l}-{\mathit{\alpha }}^{l+\mathrm{3}}\right)\left(k+l+\mathrm{4}\right)\left(k-l+\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{1}\right)}}\\ {\displaystyle }D& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{l+\mathrm{2}}-{\mathit{\alpha }}^{k+\mathrm{3}}+{\mathit{\alpha }}^l}{\mathrm{2}\mathit{\nu }\left({\mathit{\alpha }}^l-{\mathit{\alpha }}^{-l+\mathrm{1}}\right)\left(k+l+\mathrm{2}\right)\left(k-l+\mathrm{3}\right)\left(\mathrm{2}l+\mathrm{1}\right)}}.\end{array}$$

Zero-slip conditions at both boundaries lead to the following solution coefficients:

$$\begin{array}{ll}{\displaystyle }A& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{-l+\mathrm{3}}}{\mathit{\nu }}}\left[\left({\mathit{\alpha }}^{k+\mathrm{2}}+{\mathit{\alpha }}^{l-\mathrm{1}}\right)\left(k+l+\mathrm{2}\right)\left(\mathrm{2}l+\mathrm{3}\right)\right.\\ {\displaystyle }& {\displaystyle }-\left({\mathit{\alpha }}^k+{\mathit{\alpha }}^{l+\mathrm{1}}\right)\left(k+l+\mathrm{4}\right)\left(\mathrm{2}l+\mathrm{1}\right)\\ {\displaystyle }& {\displaystyle }\left.-\mathrm{2}\left({\mathit{\alpha }}^{k+\mathrm{2}l+\mathrm{3}}+{\mathit{\alpha }}^{-l-\mathrm{2}}\right)\left(k-l+\mathrm{1}\right)\right]/\mathrm{\Gamma }\\ {\displaystyle }B& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{l+\mathrm{4}}}{\mathit{\nu }}}\left[\left({\mathit{\alpha }}^{k+\mathrm{2}l+\mathrm{1}}+{\mathit{\alpha }}^{l+\mathrm{1}}\right)\left(k-l+\mathrm{3}\right)\left(\mathrm{2}l+\mathrm{1}\right)\right.\\ {\displaystyle }& {\displaystyle }-\left({\mathit{\alpha }}^{k+\mathrm{2}l+\mathrm{3}}+{\mathit{\alpha }}^{l-\mathrm{1}}\right)\left(k-l+\mathrm{1}\right)\left(\mathrm{2}l-\mathrm{1}\right)\\ {\displaystyle }& {\displaystyle }\left.-\mathrm{2}\left({\mathit{\alpha }}^{k+\mathrm{2}}+{\mathit{\alpha }}^{\mathrm{3}l}\right)\left(k+l+\mathrm{2}\right)\right]/\mathrm{\Gamma }\\ {\displaystyle }C& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{-l+\mathrm{1}}}{\mathit{\nu }}}\left[-\left({\mathit{\alpha }}^{k+\mathrm{2}}+{\mathit{\alpha }}^{l-\mathrm{3}}\right)\left(k+l+\mathrm{2}\right)\left(\mathrm{2}l+\mathrm{1}\right)\right.\\ {\displaystyle }& {\displaystyle }+\left({\mathit{\alpha }}^k+{\mathit{\alpha }}^{l-\mathrm{1}}\right)\left(k+l+\mathrm{4}\right)\left(\mathrm{2}l-\mathrm{1}\right)\\ {\displaystyle }& {\displaystyle }\left.+\mathrm{2}\left({\mathit{\alpha }}^{k+\mathrm{2}l+\mathrm{1}}+{\mathit{\alpha }}^{-l-\mathrm{2}}\right)\left(k-l+\mathrm{3}\right)\right]/\mathrm{\Gamma }\\ {\displaystyle }D& {\displaystyle }={\displaystyle \frac{g{R}_{+}^{l+\mathrm{2}}}{\mathit{\nu }}}\left[-\left({\mathit{\alpha }}^{k+\mathrm{2}l+\mathrm{1}}+{\mathit{\alpha }}^{l-\mathrm{1}}\right)\left(k-l+\mathrm{3}\right)\left(\mathrm{2}l+\mathrm{3}\right)\right.\\ {\displaystyle }& {\displaystyle }+\left({\mathit{\alpha }}^{k+\mathrm{2}l+\mathrm{3}}+{\mathit{\alpha }}^{l-\mathrm{3}}\right)\left(k-l+\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{1}\right)\\ {\displaystyle }& {\displaystyle }\left.+\mathrm{2}\left({\mathit{\alpha }}^k+{\mathit{\alpha }}^{\mathrm{3}l}\right)\left(k+l+\mathrm{4}\right)\right]/\mathrm{\Gamma }\\ {\displaystyle }\mathrm{\Gamma }& {\displaystyle }=\left(\left({\mathit{\alpha }}^{l+\mathrm{1}}+{\mathit{\alpha }}^{l-\mathrm{3}}\right){\left(\mathrm{2}l+\mathrm{1}\right)}^{\mathrm{2}}-\mathrm{2}{\mathit{\alpha }}^{l-\mathrm{1}}\left(\mathrm{2}l+\mathrm{3}\right)\left(\mathrm{2}l-\mathrm{1}\right)-\mathrm{4}{\mathit{\alpha }}^{\mathrm{3}l}\right.\\ {\displaystyle }& {\displaystyle }\left.-\mathrm{4}{\mathit{\alpha }}^{-l-\mathrm{2}}\right)\left(k+l+\mathrm{4}\right)\left(k+l+\mathrm{2}\right)\left(k-l+\mathrm{3}\right)\left(k-l+\mathrm{1}\right).\end{array}$$

2.3.2 Green's function solution – spherical

As in two dimensions, we find solutions for the case where

$$\begin{array}{}\text{(50)}& {\mathit{\rho }}^{\prime }=\mathit{\delta }(r-r^{\prime })Y_{lm}(\mathit{\theta },\mathit{\phi }),\end{array}$$

by combining two homogeneous solutions:

$$\begin{array}{}\text{(51)}& \begin{array}{rl}& \mathcal{P}(r,\mathit{\theta },\mathit{\phi })=\\ & \left\{\begin{array}{l}{\mathcal{P}}_{-}(r,\mathit{\theta },\mathit{\phi })=\left(A_{-}{r}^l+B_{-}{r}^{-l-\mathrm{1}}+C_{-}{r}^{l+\mathrm{2}}\right.\\ \left.+D_{-}{r}^{-l+\mathrm{1}}\right)Y_{lm}(\mathit{\theta },\mathit{\phi })\text{ for }{R}_{-}\le r<r^{\prime },\\ {\mathcal{P}}_{+}(r,\mathit{\theta },\mathit{\phi })=\left(A_{+}{r}^l+B_{+}{r}^{-l-\mathrm{1}}+C_{+}{r}^{l+\mathrm{2}}\right.\\ \left.+D_{+}{r}^{-l+\mathrm{1}}\right)Y_{lm}(\mathit{\theta },\mathit{\phi })\text{ for }{r}^{\prime }<r\le R_{+}.\end{array}\right.\end{array}\end{array}$$

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:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\text{continuity of}{u}_r:\\ \text{(52)}& {\displaystyle }& {\displaystyle }{\mathcal{P}}_{-}(r^{\prime },\mathit{\theta },\mathit{\phi })={\mathcal{P}}_{+}(r^{\prime },\mathit{\theta },\mathit{\phi }),{\displaystyle }& {\displaystyle }\text{ continuity of }{u}_{\mathit{\theta }}\text{ and }{u}_{\mathit{\phi }}:\\ {\displaystyle }& {\displaystyle }{\displaystyle \frac{\partial \left(r{\mathcal{P}}_{-}\right)}{\partial r}}{|}_{r=r^{\prime }}={\displaystyle \frac{\partial \left(r{\mathcal{P}}_{+}\right)}{\partial r}}{|}_{r=r^{\prime }}⟹{\displaystyle \frac{\partial {\mathcal{P}}_{-}}{\partial r}}(r^{\prime },\mathit{\theta },\mathit{\phi })\\ \text{(53)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{\partial {\mathcal{P}}_{+}}{\partial r}}(r^{\prime },\mathit{\theta },\mathit{\phi }),{\displaystyle }& {\displaystyle }\text{zero-shear condition: }\\ \text{(54)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathcal{P}}_{-}}{\partial r^{\mathrm{2}}}}(r^{\prime },\mathit{\theta },\mathit{\phi })={\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathcal{P}}_{+}}{\partial r^{\mathrm{2}}}}(r^{\prime },\mathit{\theta },\mathit{\phi }),{\displaystyle }& {\displaystyle }\text{normal-shear condition: }\\ \text{(55)}& {\displaystyle }& {\displaystyle }{\displaystyle \frac{{\partial }^{\mathrm{3}}{\mathcal{P}}_{+}}{\partial r^{\mathrm{3}}}}(r^{\prime },\mathit{\theta },\mathit{\phi })-{\displaystyle \frac{{\partial }^{\mathrm{3}}{\mathcal{P}}_{-}}{\partial r^{\mathrm{3}}}}(r^{\prime },\mathit{\theta },\mathit{\phi })={\displaystyle \frac{g{Y}_{lm}(\mathit{\theta },\mathit{\phi })}{\mathit{\nu }{r}^{\prime }}},\end{array}$$

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

$$\begin{array}{ll}{\displaystyle }{A}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^{-l+\mathrm{2}}}{\mathit{\nu }}}{\displaystyle \frac{\pm \left({\mathit{\alpha }}_{\mp }^{\mathrm{2}l-\mathrm{1}}-\mathrm{1}\right)}{\mathrm{2}\left({\mathit{\alpha }}_{\pm }^{\mathrm{2}l-\mathrm{1}}-{\mathit{\alpha }}_{\mp }^{\mathrm{2}l-\mathrm{1}}\right)\left(\mathrm{2}l-\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{1}\right)}}\\ {\displaystyle }{B}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^{l+\mathrm{3}}}{\mathit{\nu }}}{\displaystyle \frac{\pm \left({\mathit{\alpha }}_{\mp }^{-\mathrm{2}l-\mathrm{3}}-\mathrm{1}\right)}{\mathrm{2}\left({\mathit{\alpha }}_{\pm }^{-\mathrm{2}l-\mathrm{3}}-{\mathit{\alpha }}_{\mp }^{-\mathrm{2}l-\mathrm{3}}\right)\left(\mathrm{2}l+\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{3}\right)}}\\ {\displaystyle }{C}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^{-l}}{\mathit{\nu }}}{\displaystyle \frac{\pm \left(-{\mathit{\alpha }}_{\mp }^{\mathrm{2}l+\mathrm{3}}+\mathrm{1}\right)}{\mathrm{2}\left({\mathit{\alpha }}_{\pm }^{\mathrm{2}l+\mathrm{3}}-{\mathit{\alpha }}_{\mp }^{\mathrm{2}l+\mathrm{3}}\right)\left(\mathrm{2}l+\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{3}\right)}}\\ {\displaystyle }{D}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^{l+\mathrm{1}}}{\mathit{\nu }}}{\displaystyle \frac{\pm \left(-{\mathit{\alpha }}_{\mp }^{-\mathrm{2}l+\mathrm{1}}+\mathrm{1}\right)}{\mathrm{2}\left({\mathit{\alpha }}_{\pm }^{-\mathrm{2}l+\mathrm{1}}-{\mathit{\alpha }}_{\mp }^{-\mathrm{2}l+\mathrm{1}}\right)\left(\mathrm{2}l-\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{1}\right)}}.\end{array}$$

The zero-slip solution coefficients are given by

$$\begin{array}{ll}{\displaystyle }{A}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^{-l+\mathrm{2}}}{\mathit{\nu }}}[{\mathit{\alpha }}_{+}^{\mathrm{2}}-{\mathit{\alpha }}_{-}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{2}}{\mathrm{2}l+\mathrm{1}}}\left({\mathit{\alpha }}_{+}^{-\mathrm{2}l-\mathrm{1}}-{\mathit{\alpha }}_{-}^{-\mathrm{2}l-\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }\pm {\displaystyle \frac{\mathrm{2}l+\mathrm{3}}{\mathrm{2}l+\mathrm{1}}}+{\displaystyle \frac{\mathrm{2}}{\mathrm{2}l-\mathrm{1}}}\left({\mathit{\alpha }}_{+}^{\mathrm{2}}{\mathit{\alpha }}_{-}^{-\mathrm{2}l-\mathrm{1}}-{\mathit{\alpha }}_{+}^{-\mathrm{2}l-\mathrm{1}}{\mathit{\alpha }}_{-}^{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }\pm {\displaystyle \frac{\mathrm{4}{\mathit{\gamma }}^{\pm (\mathrm{2}l+\mathrm{1})}}{\left(\mathrm{2}l-\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{1}\right)}}\mp {\displaystyle \frac{{\mathit{\gamma }}^{\mp \mathrm{2}}\left(\mathrm{2}l+\mathrm{1}\right)}{\mathrm{2}l-\mathrm{1}}}]/\\ {\displaystyle }& {\displaystyle }[-\mathrm{8}{\mathit{\gamma }}^{-\mathrm{2}l-\mathrm{1}}-\mathrm{8}{\mathit{\gamma }}^{\mathrm{2}l+\mathrm{1}}+\mathrm{2}\left({\mathit{\gamma }}^{\mathrm{2}}+{\mathit{\gamma }}^{-\mathrm{2}}\right){\left(\mathrm{2}l+\mathrm{1}\right)}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }-\mathrm{4}\left(\mathrm{2}l-\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{3}\right)]\\ {\displaystyle }{B}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^{l+\mathrm{3}}}{\mathit{\nu }}}[{\mathit{\alpha }}_{+}^{\mathrm{2}}-{\mathit{\alpha }}_{-}^{\mathrm{2}}-{\displaystyle \frac{\mathrm{2}}{\mathrm{2}l+\mathrm{1}}}\left({\mathit{\alpha }}_{+}^{\mathrm{2}l+\mathrm{1}}-{\mathit{\alpha }}_{-}^{\mathrm{2}l+\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }\pm {\displaystyle \frac{\mathrm{2}l-\mathrm{1}}{\mathrm{2}l+\mathrm{1}}}-{\displaystyle \frac{\mathrm{2}}{\mathrm{2}l+\mathrm{3}}}\left({\mathit{\alpha }}_{+}^{\mathrm{2}}{\mathit{\alpha }}_{-}^{\mathrm{2}l+\mathrm{1}}-{\mathit{\alpha }}_{+}^{\mathrm{2}l+\mathrm{1}}{\mathit{\alpha }}_{-}^{\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }\pm {\displaystyle \frac{\mathrm{4}{\mathit{\gamma }}^{\mp (\mathrm{2}l+\mathrm{1})}}{\left(\mathrm{2}l+\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{3}\right)}}\mp {\displaystyle \frac{{\mathit{\gamma }}^{\mp \mathrm{2}}\left(\mathrm{2}l+\mathrm{1}\right)}{\mathrm{2}l+\mathrm{3}}}]/\\ {\displaystyle }& {\displaystyle }[-\mathrm{8}{\mathit{\gamma }}^{-\mathrm{2}l-\mathrm{1}}-\mathrm{8}{\mathit{\gamma }}^{\mathrm{2}l+\mathrm{1}}+\mathrm{2}\left({\mathit{\gamma }}^{\mathrm{2}}+{\mathit{\gamma }}^{-\mathrm{2}}\right){\left(\mathrm{2}l+\mathrm{1}\right)}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }-\mathrm{4}\left(\mathrm{2}l-\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{3}\right)]\\ {\displaystyle }{C}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^{-l}}{\mathit{\nu }}}[{\mathit{\alpha }}_{-}^{-\mathrm{2}}-{\mathit{\alpha }}_{+}^{-\mathrm{2}}+{\displaystyle \frac{\mathrm{2}}{\mathrm{2}l+\mathrm{1}}}\left({\mathit{\alpha }}_{+}^{-\mathrm{2}l-\mathrm{1}}-{\mathit{\alpha }}_{-}^{-\mathrm{2}l-\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }\mp {\displaystyle \frac{\mathrm{2}l-\mathrm{1}}{\mathrm{2}l+\mathrm{1}}}+{\displaystyle \frac{\mathrm{2}}{\mathrm{2}l+\mathrm{3}}}\left({\mathit{\alpha }}_{+}^{-\mathrm{2}}{\mathit{\alpha }}_{-}^{-\mathrm{2}l-\mathrm{1}}-{\mathit{\alpha }}_{+}^{-\mathrm{2}l-\mathrm{1}}{\mathit{\alpha }}_{-}^{-\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }\mp {\displaystyle \frac{\mathrm{4}{\mathit{\gamma }}^{\pm (\mathrm{2}l+\mathrm{1})}}{\left(\mathrm{2}l+\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{3}\right)}}\pm {\displaystyle \frac{{\mathit{\gamma }}^{\pm \mathrm{2}}\left(\mathrm{2}l+\mathrm{1}\right)}{\mathrm{2}l+\mathrm{3}}}]/\\ {\displaystyle }& {\displaystyle }[-\mathrm{8}{\mathit{\gamma }}^{-\mathrm{2}l-\mathrm{1}}-\mathrm{8}{\mathit{\gamma }}^{\mathrm{2}l+\mathrm{1}}+\mathrm{2}\left({\mathit{\gamma }}^{\mathrm{2}}+{\mathit{\gamma }}^{-\mathrm{2}}\right){\left(\mathrm{2}l+\mathrm{1}\right)}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }-\mathrm{4}\left(\mathrm{2}l-\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{3}\right)]\\ {\displaystyle }{D}_{\pm }& {\displaystyle }={\displaystyle \frac{g{r^{\prime }}^{l+\mathrm{1}}}{\mathit{\nu }}}[{\mathit{\alpha }}_{-}^{-\mathrm{2}}-{\mathit{\alpha }}_{+}^{-\mathrm{2}}-{\displaystyle \frac{\mathrm{2}}{\mathrm{2}l+\mathrm{1}}}\left({\mathit{\alpha }}_{+}^{\mathrm{2}l+\mathrm{1}}-{\mathit{\alpha }}_{-}^{\mathrm{2}l+\mathrm{1}}\right)\\ {\displaystyle }& {\displaystyle }\mp {\displaystyle \frac{\mathrm{2}l+\mathrm{3}}{\mathrm{2}l+\mathrm{1}}}-{\displaystyle \frac{\mathrm{2}}{\mathrm{2}l-\mathrm{1}}}\left({\mathit{\alpha }}_{+}^{-\mathrm{2}}{\mathit{\alpha }}_{-}^{\mathrm{2}l+\mathrm{1}}-{\mathit{\alpha }}_{+}^{\mathrm{2}l+\mathrm{1}}{\mathit{\alpha }}_{-}^{-\mathrm{2}}\right)\\ {\displaystyle }& {\displaystyle }\mp {\displaystyle \frac{\mathrm{4}{\mathit{\gamma }}^{\mp (\mathrm{2}l+\mathrm{1})}}{\left(\mathrm{2}l-\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{1}\right)}}\pm {\displaystyle \frac{{\mathit{\gamma }}^{\pm \mathrm{2}}\left(\mathrm{2}l+\mathrm{1}\right)}{\mathrm{2}l-\mathrm{1}}}]/\\ {\displaystyle }& {\displaystyle }[-\mathrm{8}{\mathit{\gamma }}^{-\mathrm{2}l-\mathrm{1}}-\mathrm{8}{\mathit{\gamma }}^{\mathrm{2}l+\mathrm{1}}+\mathrm{2}\left({\mathit{\gamma }}^{\mathrm{2}}+{\mathit{\gamma }}^{-\mathrm{2}}\right){\left(\mathrm{2}l+\mathrm{1}\right)}^{\mathrm{2}}\\ {\displaystyle }& {\displaystyle }-\mathrm{4}\left(\mathrm{2}l-\mathrm{1}\right)\left(\mathrm{2}l+\mathrm{3}\right)].\end{array}$$

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 N_j and M_l:

$$\begin{array}{}\text{(56)}& \mathit{u}=\sum _j{\mathit{u}}_j{N}_j,p=\sum _l{p}_l{M}_l.\end{array}$$

We use either the P2–P1 (Taylor–Hood) or P2_bubble–P1_DG element pairs. In both cases, the velocity and pressure basis functions N_j and M_l 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 P2_bubble–P1_DG element pair, the piecewise-linear pressure is treated as discontinuous between cells and the continuous quadratic velocity basis functions are enriched by an extra cubic “bubble” function with a corresponding cell-centred degree of freedom (Ern and Guermond, 2004; Boffi et al., 2013).

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, Eqs. (1) and (3) then become (omitting boundary terms)

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\underset{{\mathrm{\Omega }}_h}{\int }\mathit{\nu }\left(\mathrm{\nabla }{N}_i\right)\cdot \sum _j\left[\mathrm{\nabla }{N}_j^T{\mathit{u}}_j+{\mathit{u}}_j^T\mathrm{\nabla }{N}_j\right]\\ \text{(57)}& {\displaystyle }& {\displaystyle }+N_i\sum _l{p}_l\mathrm{\nabla }{M}_l=-\underset{{\mathrm{\Omega }}_h}{\int }{N}_{i}g{\mathit{\rho }}^{\prime }\widehat{\mathit{r}}\text{for all }{N}_i,\text{(58)}& {\displaystyle }& {\displaystyle }\underset{{\mathrm{\Omega }}_h}{\int }\left(\mathrm{\nabla }{M}_k\right)\cdot \sum _j{\mathit{u}}_j{N}_j=\mathrm{0}\text{for all }{M}_k.\end{array}$$

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 3 Convergence 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 $X_{\text{lin}}:\mathit{\xi }\mapsto X_{\text{lin}}\left(\mathit{\xi }\right)$ 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:

$$\begin{array}{}\text{(59)}& X_{\text{quad}}\left(\mathit{\xi }\right)\stackrel{\text{quad. interp.}}{\approx }{\displaystyle \frac{r_{\text{lin}}\left(\mathit{\xi }\right)}{\Vert X_{\text{lin}}\left(\mathit{\xi }\right)\Vert }}{X}_{\text{lin}}\left(\mathit{\xi }\right),\end{array}$$

at the standard Lagrange nodes of the quadratic function space. Here, r_lin(ξ) is the linear interpolation of the radius, i.e. $r_{\text{lin}}=\Vert X_{\text{lin}}\Vert$ 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 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+\mathrm{1}$.

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

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }-\underset{{\mathrm{\Omega }}_h}{\int }{N}_{i}g{\mathit{\rho }}^{\prime }\widehat{\mathit{r}}=-\underset{{\mathrm{\Omega }}_h}{\int }{N}_i(r,\mathit{\theta })g\cos \left(n\mathit{\theta }\right)\mathit{\delta }(r-r^{\prime })\widehat{\mathit{r}}=\\ \text{(60)}& {\displaystyle }& {\displaystyle }-\underset{{\mathrm{\Gamma }}^{\prime }}{\int }{N}_i(r,\mathit{\theta })g\cos \left(n\mathit{\theta }\right)\mathit{n},{\displaystyle }& {\displaystyle }-\underset{{\mathrm{\Omega }}_h}{\int }{N}_i(r,\mathit{\theta })g{\mathit{\rho }}^{\prime }\widehat{\mathit{r}}=-\underset{{\mathrm{\Omega }}_h}{\int }{N}_i(r,\mathit{\theta },\mathit{\phi })g{Y}_{lm}(\mathit{\theta },\mathit{\phi })\mathit{\delta }(r-r^{\prime })\widehat{\mathit{r}}\\ \text{(61)}& {\displaystyle }& {\displaystyle }=-\underset{{\mathrm{\Gamma }}^{\prime }}{\int }{N}_i(r,\mathit{\theta },\mathit{\phi })g{Y}_{lm}(\mathit{\theta },\mathit{\phi })\mathit{n},\end{array}$$

where Γ^′ is an internal boundary at $r=r^{\prime }$ oriented such that its normal $\mathit{n}=\widehat{\mathit{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:

$$\begin{array}{}\text{(62)}& p\to p-{\displaystyle \frac{{\int }_{{\mathrm{\Omega }}_h}p}{{\int }_{{\mathrm{\Omega }}_h}\mathrm{1}}}.\end{array}$$

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

$$\begin{array}{}\text{(63)}& \mathit{u}\to \mathit{u}-{\displaystyle \frac{{\int }_{{\mathrm{\Omega }}_h}r\widehat{\mathit{\theta }}\cdot \mathit{u}}{{\int }_{{\mathrm{\Omega }}_h}{r}^{\mathrm{2}}}}r\widehat{\mathit{\theta }},\end{array}$$

which ensures that the angular momentum $\int r\widehat{\mathit{\theta }}\cdot \mathit{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 xⁱ and zero eigenvector λ, we get

$$\begin{array}{}\text{(64)}& {\mathit{x}}^{\mathrm{i}}\to {\mathit{x}}^{\mathrm{i}}-\langle \mathit{\lambda },{\mathit{x}}^{\mathrm{i}}\rangle \mathit{\lambda }.\end{array}$$

This l₂ projection, based on the l₂ inner product $\langle \cdot ,\cdot \rangle$, is analogous but not equivalent to the projections in Eqs. (62) and (63) (which are L₂ projections). Therefore, despite the l₂ projections fixing the null modes during the iterative solve, the L₂ 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 $\mathrm{128}\times \mathrm{16}\times \mathrm{2}=\mathrm{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_{-}=\mathrm{1.22}$ and $R_{+}=\mathrm{2.22}$, and the delta-function cases used $r^{\prime }=(R_{-}+R_{+})/\mathrm{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 L₂ norm:

$$\begin{array}{}\text{(65)}& \begin{array}{rl}{\displaystyle \frac{\Vert \mathit{u}-{\mathit{u}}^{\ast }{\Vert }_{\mathrm{2}}}{\Vert {\mathit{u}}^{\ast }{\Vert }_{\mathrm{2}}}}=& {\displaystyle \frac{\sqrt{{\int }_{{\mathrm{\Omega }}_h}|\mathit{u}-{\mathit{u}}^{\ast }{|}^{\mathrm{2}}}}{\sqrt{{\int }_{{\mathrm{\Omega }}_h}|{\mathit{u}}^{\ast }{|}^{\mathrm{2}}}}},\\ {\displaystyle \frac{\Vert p-p^{\ast }{\Vert }_{\mathrm{2}}}{\Vert p^{\ast }{\Vert }_{\mathrm{2}}}}=& {\displaystyle \frac{\sqrt{{\int }_{{\mathrm{\Omega }}_h}(p-p^{\ast }{)}^{\mathrm{2}}}}{\sqrt{{\int }_{{\mathrm{\Omega }}_h}(p^{\ast }{)}^{\mathrm{2}}}}}.\end{array}\end{array}$$

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+\mathrm{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 King, 2019). These will be discussed in more detail in the following section.

The derivation for all 3-D cases in this paper relies on the Mie representation that decomposes the velocity solution into poloidal and toroidal components, through which, under the assumption of purely poloidal flow, the Stokes equations can be reduced to a biharmonic equation (10). Any solution to the inhomogeneous solution can then be combined with four linearly independent homogeneous solutions to this equation, the coefficients of which can be derived through the imposition of boundary conditions. For the smooth case with a generic monomial forcing term, the same decomposition (i.e. Eqs. 45 and 44) was used in Tosi and Martinec (2007) to derive the analytical solution for Stokes flow in two eccentrically nested spheres.

In Horbach et al. (2020), similar techniques were employed to derive benchmarks in spherical-shell domains that satisfy zero-slip and free-slip conditions and, in addition, a mixed zero-slip/free-slip case. Here, the derivation starts by simply selecting four, in principle arbitrary, linearly independent solutions for the radially dependent part of the poloidal scalar function. Again the imposition of boundary conditions fixes the coefficients of this linear combination. The corresponding right-hand side forcing term is then obtained by substitution.

The number of published benchmarks for 2-D cylindrical-shell domains is more limited (e.g. Buffa et al., 2011; Blinova et al., 2016; Hoang et al., 2017). The derivation of the equivalent cases in 2-D cylindrical-shell domains is somewhat simpler, but also ultimately relies on combining four independent homogeneous solutions and one inhomogeneous solution to a biharmonic equation.

In Blinova et al. (2016), analytical solutions in both cylindrical and spherical domains are presented for the Stokes equations with a radially dependent viscosity. Because these are solutions in cylindrical or spherical coordinates without reference to any specific domain, they do not satisfy no-normal-flow conditions on the boundary of a shell domain. They can be used in such domains as a numerical benchmark by specifying all velocity components of the analytical solution as a Dirichlet condition for the model. Analytical solutions for radially dependent viscosity were also presented in Thieulot (2017). Their solutions (in 3-D only) do satisfy no-normal-flow conditions in a spherical-shell domain, but the tangential components are non-zero at the boundary and thus still require inhomogeneous Dirichlet conditions. Spatially varying viscosity is of course an import aspect of mantle convection models for which these are effective benchmarks. The isoviscous solutions presented here, and those in Horbach et al. (2020), however, allow for the testing of zero-slip and free-slip conditions, where in particular free-slip conditions may pose various numerical challenges such as rotational modes and, depending on the discretisation used, the (non-)alignment of velocity components with normal and tangential directions at the boundary.

5.2 Reduced order of convergence with discontinuous pressures

At first sight, the reduced order of convergence for the delta-function cases seems at odds with those expected for the P2–P1 element pair. However, the mathematical proofs for the ideal order of convergence to solutions of the Stokes equations rely on certain regularity assumptions of the right-hand side forcing term and, related to that, on the regularity of the velocity and pressure solutions. The regularity of the delta function can be classified as being a member of the Sobolev space H⁻¹(Ω) the dual of the Sobolev space H¹(Ω), 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:

$$\begin{array}{}\text{(66)}& {\mathit{\delta }}_{r^{\prime }}:v\mapsto \underset{\mathrm{\Omega }}{\int }\mathit{\delta }(r-r^{\prime })v(r,\mathit{\varphi })r\mathrm{d}r\mathrm{d}\mathit{\varphi }=\underset{\mathrm{\Omega }}{\int }v\left(r^{\prime }\right)r^{\prime }\mathrm{d}\mathit{\varphi },\end{array}$$

which maps functions v∈H¹(Ω), 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 H^m(Ω) with $m\ge -\mathrm{1}$, the Stokes equations in the weak form have a unique solution (given sufficient integral constraints) with velocity in H^m+2(Ω) and pressure in H^m+1(Ω). The analytical solutions derived here for the delta case with $m=-\mathrm{1}$ indeed have a discontinuous pressure in H⁰(Ω)=L₂(Ω) and a velocity with discontinuous normal derivative in H¹(Ω).

For velocity–pressure finite element pairs that satisfy the standard inf-sup, or Ladyzhenskaya–Babuška–Brezzi (LBB) condition

$$\begin{array}{}\text{(67)}& \underset{\mathit{v}\in V}{sup}{\displaystyle \frac{{\int }_{\mathrm{\Omega }}\mathit{v}\cdot \mathrm{\nabla }q}{\Vert \mathit{v}{\Vert }_{\mathrm{1}}}}\ge \mathit{\beta }\Vert q{\Vert }_{\mathrm{2}}\text{ for all }q\in W,\end{array}$$

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

$$|{\mathit{u}}^{\ast }-\mathit{u}{|}_{\mathrm{1}}+\Vert p^{\ast }-p{\Vert }_{\mathrm{2}}\le C_{\mathrm{1}}\left\{\underset{\mathit{v}\in V}{inf}|{\mathit{u}}^{\ast }-\mathit{v}{|}_{\mathrm{1}}+\underset{q\in W}{inf}\Vert p^{\ast }-q{\Vert }_{\mathrm{2}}\right\},$$

where C₁ 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 $|\cdot {|}_{\mathrm{1}}$ is the semi-norm

$$\begin{array}{}\text{(68)}& |\mathit{u}{|}_{\mathrm{1}}=\sqrt{\sum _{i=\mathrm{1}}^{\text{dim}}\underset{\mathrm{\Omega }}{\int }\Vert {\partial }_i\mathit{u}{\Vert }^{\mathrm{2}}}.\end{array}$$

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

$$\begin{array}{}\text{(69)}& \underset{q\in W}{inf}\Vert p^{\ast }-q{\Vert }_{\mathrm{2}}\le C_{\mathrm{2}}{h}^{\frac{\mathrm{1}}{\mathrm{2}}}\Vert p^{\ast }{\Vert }_{\mathrm{2}}\end{array}$$

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

Figure 6 Convergence 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 P2_bubble–P1_DG 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=P1_DG the space of piecewise-linear but discontinuous functions. To satisfy the inf-sup condition, this requires enriching the quadratic function space P₂ 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.

Figure 7 Convergence of velocity and radial stresses at top (a, c) top and (b, d) bottom surfaces ($r=R_{\pm }$), 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 $\Vert p-p^{\ast }{\Vert }_{\mathrm{2}}$ 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 ${\overline{p}}^{\ast }$, second-order convergence can sometimes be observed in $\Vert p-{\overline{p}}^{\ast }{\Vert }_{\mathrm{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; Kramer, 2020).

We verify the convergence of the P2–P1 (Taylor–Hood) finite element discretisation using Fluidity (Davies et al., 2011; Kramer et al., 2012). The continuous approximation of pressure can lead to a reduced order of convergence in the presence of discontinuities, which can be overcome using a discontinuous numerical approximation of pressure. It is important to note that this reduced order of convergence was only observed by comparing the numerical solution with the entire analytical solution in the interior of the domain. A comparison based on surface response only failed to highlight this issue.

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

$$\begin{array}{}\text{(A1)}& u_r=-{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\phi }}},u_{\mathit{\phi }}={\displaystyle \frac{\partial \mathit{\psi }}{\partial }}r.\end{array}$$

We make use of the following expressions for the derivatives of the unit vectors $\widehat{\mathit{r}}$ and $\widehat{\mathit{\phi }}$ with respect to r and φ:

$$\begin{array}{}\text{(A2)}& {\displaystyle }\widehat{\mathit{r}}\cdot \mathrm{\nabla }\widehat{\mathit{r}}& {\displaystyle }={\displaystyle \frac{\partial \widehat{\mathit{r}}}{\partial r}}=\mathrm{0},& {\displaystyle }\widehat{\mathit{\phi }}\cdot \mathrm{\nabla }\widehat{\mathit{r}}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \widehat{\mathit{r}}}{\partial \mathit{\phi }}}={\displaystyle \frac{\mathrm{1}}{r}}\widehat{\mathit{\phi }},\text{(A3)}& {\displaystyle }\widehat{\mathit{r}}\cdot \mathrm{\nabla }\widehat{\mathit{\phi }}& {\displaystyle }={\displaystyle \frac{\partial \widehat{\mathit{\phi }}}{\partial r}}=\mathrm{0},& {\displaystyle }\widehat{\mathit{\phi }}\cdot \mathrm{\nabla }\widehat{\mathit{\phi }}& {\displaystyle }={\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \widehat{\mathit{\phi }}}{\partial \mathit{\phi }}}=-{\displaystyle \frac{\mathrm{1}}{r}}\widehat{\mathit{r}}.\end{array}$$

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

$$\begin{array}{}\text{(A4)}& {\displaystyle }\begin{array}{rl}& {\mathit{\tau }}_{rr}=\mathrm{2}\mathit{\nu }\widehat{\mathit{r}}\cdot \left[\mathrm{\nabla }\mathit{u}\right]\cdot \widehat{\mathit{r}}=\mathrm{2}\mathit{\nu }\widehat{\mathit{r}}\cdot \mathrm{\nabla }\left(\widehat{\mathit{r}}\cdot \mathit{u}\right)-\mathrm{2}\mathit{\nu }\left(\widehat{\mathit{r}}\cdot \mathrm{\nabla }\widehat{\mathit{r}}\right)\\ & \cdot \mathit{u}=\mathrm{2}\mathit{\nu }{\displaystyle \frac{\partial u_r}{\partial r}}-\mathrm{0}=-\mathrm{2}\mathit{\nu }{\displaystyle \frac{\partial }{\partial r}}\left({\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\phi }}}\right),\end{array}\text{(A5)}& {\displaystyle }\begin{array}{rl}& {\mathit{\tau }}_{\mathit{\phi }\mathit{\phi }}=\mathrm{2}\mathit{\nu }\widehat{\mathit{\phi }}\cdot \left[\mathrm{\nabla }\mathit{u}\right]\cdot \widehat{\mathit{\phi }}=\mathrm{2}\mathit{\nu }\widehat{\mathit{\phi }}\cdot \mathrm{\nabla }\left(\widehat{\mathit{\phi }}\cdot \mathit{u}\right)\\ & -\mathrm{2}\mathit{\nu }\left(\widehat{\mathit{\phi }}\cdot \mathrm{\nabla }\widehat{\mathit{\phi }}\right)\mathit{u}={\displaystyle \frac{\mathrm{2}\mathit{\nu }}{r}}{\displaystyle \frac{\partial u_{\mathit{\phi }}}{\partial \mathit{\phi }}}+{\displaystyle \frac{\mathrm{2}\mathit{\nu }}{r}}\widehat{\mathit{r}}\cdot \mathit{u}\\ & ={\displaystyle \frac{\mathrm{2}\mathit{\nu }}{r}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial r\partial \mathit{\phi }}}-{\displaystyle \frac{\mathrm{2}\mathit{\nu }}{r^{\mathrm{2}}}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\phi }}}=\mathrm{2}\mathit{\nu }{\displaystyle \frac{\partial }{\partial r}}\left({\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\phi }}}\right),\end{array}\text{(A6)}& {\displaystyle }\begin{array}{rl}& {\mathit{\tau }}_{r\mathit{\phi }}=\mathit{\nu }\widehat{\mathit{r}}\cdot \left[\mathrm{\nabla }\mathit{u}\right]\cdot \widehat{\mathit{\phi }}+\mathit{\nu }\widehat{\mathit{\phi }}\cdot \left[\mathrm{\nabla }\mathit{u}\right]\cdot \widehat{\mathit{r}}=\mathit{\nu }\widehat{\mathit{r}}\cdot \mathrm{\nabla }\left(\mathit{u}\cdot \widehat{\mathit{\phi }}\right)\\ & -\mathit{\nu }\left(\widehat{r}\cdot \mathrm{\nabla }\widehat{\mathit{\phi }}\right)\cdot \mathit{u}+\mathit{\nu }\widehat{\mathit{\phi }}\cdot \mathrm{\nabla }\left(\mathit{u}\cdot \widehat{\mathit{r}}\right)-\mathit{\nu }\left(\widehat{\mathit{\phi }}\cdot \mathrm{\nabla }\widehat{\mathit{r}}\right)\cdot \mathit{u}\\ & =\mathit{\nu }{\displaystyle \frac{\partial }{\partial r}}\left({\displaystyle \frac{\partial \mathit{\psi }}{\partial r}}\right)-\mathrm{0}+{\displaystyle \frac{\mathit{\nu }}{r}}{\displaystyle \frac{\partial u_r}{\partial \mathit{\phi }}}-{\displaystyle \frac{\mathit{\nu }}{r}}\widehat{\mathit{\phi }}\cdot \mathit{u}\\ & =\mathit{\nu }\left({\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial r^{\mathrm{2}}}}-{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}\left({\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\phi }}}\right)-{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial r}}\right).\end{array}\end{array}$$

Note that, as expected, ${\mathit{\tau }}_{rr}+{\mathit{\tau }}_{\mathit{\phi }\mathit{\phi }}=\mathrm{0}$. In the same way, we derive the following expression for the vorticity, or curl of the velocity:

$$\begin{array}{}\text{(A7)}& \begin{array}{rl}& \mathrm{curl}\mathit{u}=\widehat{\mathit{r}}\cdot \left[\mathrm{\nabla }\mathit{u}\right]\cdot \widehat{\mathit{\phi }}-\widehat{\mathit{\phi }}\cdot \left[\mathrm{\nabla }\mathit{u}\right]\cdot \widehat{\mathit{r}}=\\ & {\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial r^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}\left({\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial \mathit{\phi }}}\right)+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial r}}={\mathrm{\nabla }}^{\mathrm{2}}\mathit{\psi }.\end{array}\end{array}$$

The viscosity term in the Stokes equations can be written as

$$\begin{array}{ll}\text{(A8)}& {\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot \mathit{\tau }=\mathrm{\nabla }\cdot \left(\widehat{\mathit{r}}{\mathit{\tau }}_{rr}\widehat{\mathit{r}}+\widehat{\mathit{r}}{\mathit{\tau }}_{r\mathit{\phi }}\widehat{\mathit{\phi }}+\widehat{\mathit{\phi }}{\mathit{\tau }}_{\mathit{\phi }r}\widehat{\mathit{r}}+\widehat{\mathit{\phi }}{\mathit{\tau }}_{\mathit{\phi }\mathit{\phi }}\widehat{\mathit{\phi }}\right){\displaystyle }& {\displaystyle }=\left(\mathrm{\nabla }\cdot \widehat{\mathit{r}}\right)\left({\mathit{\tau }}_{rr}\widehat{\mathit{r}}+{\mathit{\tau }}_{r\mathit{\phi }}\widehat{\mathit{\phi }}\right)+\widehat{\mathit{r}}\cdot \mathrm{\nabla }\left({\mathit{\tau }}_{rr}\widehat{\mathit{r}}+{\mathit{\tau }}_{r\mathit{\phi }}\widehat{\mathit{\phi }}\right)\\ \text{(A9)}& {\displaystyle }& {\displaystyle }+\left(\mathrm{\nabla }\cdot \widehat{\mathit{\phi }}\right)\left({\mathit{\tau }}_{\mathit{\phi }r}\widehat{\mathit{r}}+{\mathit{\tau }}_{\mathit{\phi }\mathit{\phi }}\widehat{\mathit{\phi }}\right)+\widehat{\mathit{\phi }}\cdot \mathrm{\nabla }\left({\mathit{\tau }}_{\mathit{\phi }r}\widehat{\mathit{r}}+{\mathit{\tau }}_{\mathit{\phi }\mathit{\phi }}\widehat{\mathit{\phi }}\right).\end{array}$$

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

$$\begin{array}{}\text{(A10)}& \mathrm{\nabla }\cdot \widehat{\mathit{r}}={\displaystyle \frac{\mathrm{1}}{r}},\mathrm{\nabla }\cdot \widehat{\mathit{\phi }}=\mathrm{0},\end{array}$$

and the fact that ${\mathit{\tau }}_{\mathit{\phi }\mathit{\phi }}=-{\mathit{\tau }}_{rr}$. After reordering to group the radial and transverse components, this leads to

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot \mathit{\tau }=\left[{\displaystyle \frac{\mathrm{2}}{r}}{\mathit{\tau }}_{rr}+{\displaystyle \frac{\partial {\mathit{\tau }}_{rr}}{\partial r}}+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}{\mathit{\tau }}_{\mathit{\phi }r}\right]\widehat{\mathit{r}}\\ \text{(A11)}& {\displaystyle }& {\displaystyle }+\left[{\displaystyle \frac{\mathrm{2}}{r}}{\mathit{\tau }}_{r\mathit{\phi }}+{\displaystyle \frac{\partial {\mathit{\tau }}_{r\mathit{\phi }}}{\partial r}}-{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial {\mathit{\tau }}_{rr}}{\partial \mathit{\phi }}}\right]\widehat{\mathit{\phi }}{\displaystyle }& {\displaystyle }=\left[-{\displaystyle \frac{\mathrm{2}\mathit{\nu }}{r}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial r^{\mathrm{2}}}}+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}{\mathit{\tau }}_{\mathit{\phi }r}\right]\widehat{\mathit{r}}+\left[\mathit{\nu }{\displaystyle \frac{{\partial }^{\mathrm{3}}\mathit{\psi }}{\partial r^{\mathrm{3}}}}+{\displaystyle \frac{\mathit{\nu }}{r}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial r^{\mathrm{2}}}}-{\displaystyle \frac{\mathit{\nu }}{r^{\mathrm{2}}}}{\displaystyle \frac{\partial \mathit{\psi }}{\partial r}}\right.\\ \text{(A12)}& {\displaystyle }& {\displaystyle }\left.+{\displaystyle \frac{\mathit{\nu }}{r^{\mathrm{2}}}}{\displaystyle \frac{{\partial }^{\mathrm{3}}\mathit{\psi }}{\partial {\mathit{\phi }}^{\mathrm{2}}\partial r}}-{\displaystyle \frac{\mathrm{2}\mathit{\nu }}{r^{\mathrm{3}}}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\mathit{\psi }}{\partial {\mathit{\phi }}^{\mathrm{2}}}}\right]\widehat{\mathit{\phi }}\text{(A13)}& {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\mathit{\nu }}{r}}\left[{\displaystyle \frac{\partial }{\partial \mathit{\phi }}}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\psi }\right]\widehat{\mathit{r}}+\mathit{\nu }\left[{\displaystyle \frac{\partial }{\partial r}}{\mathrm{\nabla }}^{\mathrm{2}}\mathit{\psi }\right]\widehat{\mathit{\phi }}.\end{array}$$

In combination with the pressure gradient term in polar coordinates

$$\begin{array}{}\text{(A14)}& \mathrm{\nabla }p={\displaystyle \frac{\partial p}{\partial r}}\widehat{\mathit{r}}+{\displaystyle \frac{\mathrm{1}}{r}}{\displaystyle \frac{\partial p}{\partial \mathit{\phi }}}\widehat{\mathit{\phi }},\end{array}$$

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