# Full sphere hydrodynamic and dynamo benchmarks

P. Marti, N. Schaeffer, R. Hollerbach, D. Cébron, C. Nore, F. Luddens, J. -L. Guermond, J. Aubert, S. Takehiro, Y. Sasaki, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

P. Marti, N. Schaeffer, R. Hollerbach, D. Cébron, C. Nore, et al., Full sphere hydrodynamic and dynamo benchmarks. Geophysical Journal International, 2014, 197, pp.119-134. 10.1093/gji/ggt518. insu-03581136

HAL Id: insu-03581136

https://insu.hal.science/insu-03581136

Submitted on 19 Feb 2022

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>
# Full sphere hydrodynamic and dynamo benchmarks

P. Marti,<fnref n="1" />,<fnref n="2" /> N. Schaeffer,<fnref n="3" />,<fnref n="4" /> R. Hollerbach,<fnref n="1" />,<fnref n="5" /> D. Cébron,<fnref n="1" /> C. Nore,<fnref n="6" />,<fnref n="7" /> F. Luddens,<fnref n="6" /> J.-L. Guermond,<fnref n="8" /> J. Aubert,<fnref n="9" /> S. Takehiro,<fnref n="10" /> Y. Sasaki,<fnref n="11" /> Y.-Y. Hayashi,<fnref n="12" /> R. Simitev,<fnref n="13" />,<fnref n="14" />,<fnref n="15" /> F. Busse,<fnref n="14" />,<fnref n="16" /> S. Vantieghem<fnref n="1" /> and A. Jackson<fnref n="1" />

Accepted 2013 December 19. Received 2013 December 12; in original form 2013 March 8

## SUMMARY

Convection in planetary cores can generate fluid flow and magnetic fields, and a number of sophisticated codes exist to simulate the dynamic behaviour of such systems. We report on the first community activity to compare numerical results of computer codes designed to calculate fluid flow within a whole sphere. The flows are incompressible and rapidly rotating and the forcing of the flow is either due to thermal convection or due to moving boundaries. All problems defined have solutions that allow easy comparison, since they are either steady, slowly drifting or perfectly periodic. The first two benchmarks are defined based on uniform internal heating within the sphere under the Boussinesq approximation with boundary conditions that are uniform in temperature and stress-free for the flow. Benchmark 1 is purely hydrodynamic, and has a drifting solution. Benchmark 2 is a magnetohydrodynamic benchmark that can generate oscillatory, purely periodic, flows and magnetic fields. In contrast, Benchmark 3 is a hydrodynamic rotating bubble benchmark using no slip boundary conditions that has a stationary solution. Results from a variety of types of code are reported, including codes that are fully spectral (based on spherical harmonic expansions in angular coordinates and polynomial expansions in radius), mixed spectral and finite difference, finite volume, finite element and also a mixed Fourier–finite element code. There is good agreement between codes. It is found that in Benchmarks 1 and 2, the approximation of a whole sphere problem by a domain that is a spherical shell (a sphere possessing an inner core) does not represent an adequate approximation to the system, since the results differ from whole sphere results.

**Key words:** Numerical solutions; Non-linear differential equations; Dynamo: theories and simulations; Planetary interiors.

## 1 INTRODUCTION

The predominant theory for the generation mechanism of the Earth's magnetic field is that of magnetic field generation by thermal

and compositional convection, creating the so-called self-excited dynamo mechanism. Beginning with the first 3-D self-consistent Boussinesq models of thermal convection (Glatzmaier & Roberts 1995; Kageyama et al. 1995), there has been burgeoning interest

<footnote n="1">¹Institute of Geophysics, ETH Zurich, Zurich 8092, Switzerland</footnote>

<footnote n="2">²Department of Applied Mathematics, Engineering Center, Room ECOT 521, University of Colorado, Boulder, CO 80309-0526, USA.</footnote>

<footnote>E-mail: philippe.marti@colorado.edu</footnote>

<footnote n="3">³Univ. Grenoble Alpes, ISTERRE, F-38041 Grenoble, France</footnote>

<footnote n="4">⁴CNRS, ISTERRE, F-38041 Grenoble, France</footnote>

<footnote n="5">⁵Department of Applied Mathematics, University of Leeds, Leeds LS2 9JT, UK</footnote>

<footnote n="6">⁶Laboratoire d'Informatique pour la Mécanique et les Sciences de l'Ingénieur, CNRS UPR 3251, BP 133, 91403 Orsay Cedex, France, and Université Paris-Sud 11</footnote>

<footnote n="7">⁷Institut Universitaire de France, 103, bd Saint-Michel, 75005 Paris, France</footnote>

<footnote n="8">⁸Department of Mathematics, Texas A&M University 3368 TAMU, College Station, TX 77843-3368, USA</footnote>

<footnote n="9">⁹Dynamique des Fluides Géologiques, Institut de Physique du Globe de Paris, Sorbonne Paris Cité, Université Paris Diderot, UMR 7154 CNRS, F-75005 Paris, France</footnote>

<footnote n="10">¹⁰Research Institute for Mathematical Sciences, Kyoto University, Japan</footnote>

<footnote n="11">¹¹Department of Mathematics, Kyoto University, Japan</footnote>

<footnote n="12">¹²Center for Planetary Science/Department of Earth and Planetary Sciences, Kobe University, Japan</footnote>

<footnote n="13">¹³Solar Physics, Hansen Experimental Physics Laboratory, Stanford University, CA 94305, USA</footnote>

<footnote n="14">¹⁴Department of Earth and Space Sciences, University of California, Los Angeles, CA 90095, USA</footnote>

<footnote n="15">¹⁵School of Mathematics and Statistics, University of Glasgow, Glasgow, G12 8QW, UK</footnote>

<footnote n="16">¹⁶Institute of Physics, University of Bayreuth, 95440 Bayreuth, Germany</footnote>
in numerical solutions to the underlying equations of momentum conservation, magnetic field generation and heat transfer. Given the complexity and non-linearity of the physics, it has been of importance to verify that the codes used correctly calculate solutions to the underlying equations, and also to provide simple solutions that allow newly developed codes to be checked for accuracy. It is now over 11 yr since the undoubtedly successful numerical dynamo benchmark exercise of Christensen et al. (2001); Christensen et al. (2009), hereinafter B1. This benchmark exercise was set in the geometry of a spherical shell, with convection driven by a temperature difference between an inner core and the outer boundary of the spherical shell. In this respect, the computational domain is similar to that of the Earth, possessing as it does a small inner core. Three different benchmarks were devised, the first being purely thermal convection, and the second and third being dynamos (supporting magnetic fields). The latter two benchmarks differed in the treatment of the inner core: in one case the inner core was taken to be electrically insulating and fixed in the rotating frame, and in the other case the core was taken to be electrically conducting and free to rotate in response to torques that are applied to it, arising from the convection in the outer core. Central to these benchmarks was the fact that all of them possessed simple solutions, in the form of steadily drifting convection. As a result, energies are constant and, together with other diagnostics, these provide very clear solutions that could be reproduced by different numerical techniques. A measure of the success of this exercise is given by the fact that it has been used by numerous groups to check their codes.

The present benchmarking exercise is one of two brethren designed to broaden the scope of the original B1 and to provide further accurate solutions for a new generation of computer codes. The associated exercise by Jackson et al. (2013) is also set in a spherical shell as in B1; it is similar to B1 but has been designed to be particularly amenable to computer codes based on local (rather than spectral or global) descriptions of the temperature, magnetic and velocity fields. Thus, that benchmark allows comprehensive checking of finite element, finite volume and similar computer codes, as a result of the implementation of a local rather than global magnetic boundary condition. This paper treats a similar situation to B1 but differs in the removal of the inner core, and thus treats only a whole sphere. Flows in the first two benchmarks thus defined are driven by thermal convection, again under the Boussinesq approximation, and in the third by a boundary forcing. There are two reasons for defining benchmarks set in the whole sphere rather than the spherical shell. First, the whole sphere represents a canonical problem, surely a simpler geometry than the shell. There is one less degree of freedom, since the aspect ratio of the shell is no longer a defined parameter. Secondly, in the context of rapidly rotating fluid dynamics, there is likely to be a simplification in the flow structures generated because of the absence of an inner core. It is well known that the dynamics of rapidly rotating systems is dominated by the Coriolis force, thus leading to the Proudman–Taylor constraint, the alignment of flow structures with the rotation (z) axis. In a spherical shell when viscosity is reduced, as one moves from outside the so-called tangent cylinder (the cylinder that just encloses the inner core) to inside the tangent cylinder, a jump is present in the length of a column in the z-direction. Hence, there is the possibility of the need to resolve very fine shear layers in this region; for a recent discussion, see Livermore & Hollerbach (2012). The presence of very fine structures that need to be resolved can have very deleterious effects on a numerical method, particularly a spectral method based on spherical harmonics (again, for a discussion, see Livermore 2012).

Thus, the choice of a full sphere is likely to be advantageous in the limit in which the viscosity is dropped to insignificant levels.

We note in passing that the whole sphere geometry is particularly relevant to the generation of magnetic fields in the early Earth, prior to the formation of the inner core. In this time period, the convection in the core was driven by secular cooling (and possibly internal heating), and this is precisely the scenario studied here in Benchmarks 1 and 2. Associated with this geometry is a possible numerical obstacle that has perhaps been responsible for the dearth of full-sphere calculations in the literature. Working in a spherical coordinate system (r, θ, φ), that is presumably convenient from the point of view of boundary conditions, the presence of the origin of the spherical coordinate system (r = 0) in the integration domain leads to additional numerical challenges. The results presented here show that the employed methods are able to correctly handle this singularity in coordinates.

The Benchmarks 1 and 2 set up here differ from those of B1 in their use of stress-free boundary conditions, rather than non-slip conditions. This arose purely as a result of our survey of parameter space while searching for whole-sphere dynamos that possess simple solutions with clear diagnostics suitable for a benchmark. In performing this survey, we did not find a dynamo that had a steady character similar to that in B1; this is not to say that one does not exist. The dynamo solution in Benchmark 2 shows an exact periodic character with energy conversion between kinetic and magnetic forms. It thus allows very precise comparison. The use of stress-free boundaries can cause problems with angular momentum conservation (see the discussion in Jones et al. 2011), but these were handled gracefully in the solutions we report.

We mention in passing the other benchmarks that have recently been provided to the community. A new benchmark for anelastic convection has recently been described by Jones et al. (2011) and already used as a comparison for the newly developed code of Zhan et al. (2012). This benchmark again is set in a spherical shell, but has a background state with a very large change in density across the shell. Three solutions are again compared, the first two (pure thermal convection and dynamo action, respectively) possessing simple drifting solutions with very precise diagnostics. A solar mean field benchmark has also recently been provided by Jouve et al. (2008).

The layout of the paper is as follows: in Section 2, we describe the physical problems to be addressed. Benchmarks 1 and 2 are driven by internal heating and Benchmark 3 by boundary forcing. In Section 3, we give a brief overview of the different numerical methods that have been employed by the different contributing teams. In Section 4, we present and discuss the results from the different codes.

## 2 TEST CASES

Three benchmarks for incompressible flows in a rapidly rotating whole sphere are considered. The first two test problems, Benchmarks 1 and 2, are subject to the thermal forcing of a homogeneous distribution of heat sources in the volume. Benchmark 1 is a purely hydrodynamic problem while Benchmark 2 consists of a self-sustained dynamo problem. Benchmark 3 extends the scope of these test cases by considering the mechanical forcing induced by moving boundaries. In all cases, the system consists of a whole sphere of radius $r_o$, filled with a fluid of density ρ and a kinematic viscosity v. The system rotates at a rotation rate Ω. The fluid motion
is described by the velocity field **u** and, for Benchmarks 1 and 2, the temperature field is denoted by T.

## 2.1 Benchmark 1: thermal convection

Benchmark 1 is a purely hydrodynamic problem with the motion of the fluid described in the reference frame of the mantle. The system is described within the frame of the Boussinesq approximation, neglecting the density fluctuations except for the ones responsible for the buoyancy. Under the action of a gravitational field

$$g = g \frac{r}{r_o} \qquad (1)$$

and in the presence of a homogeneous heat source distribution S, the basic state is given by

$$T_b = \frac{\beta}{2} (r_o^2 - r^2) \qquad (2)$$

with $\beta = S/3\kappa$, where $\kappa$ is the thermal diffusivity. The equations are non-dimensionalized using the radius $r_o$ as length scale, the diffusion time $r_o^2/\nu$ as timescale and $\beta r_o^2$ as temperature scale. The three non-dimensional parameters are chosen to be the Ekman number E

$$E = \frac{\nu}{2\Omega r_o^2}, \qquad (3)$$

the Prandtl number Pr

$$Pr = \frac{\nu}{\kappa} \qquad (4)$$

and the modified Rayleigh number Ra

$$Ra = \frac{g\alpha\beta r_o^3}{2\Omega\kappa}, \qquad (5)$$

with $\alpha$ the thermal expansion coefficient. The motion of the fluid is then described by the non-dimensional Navier-Stokes equation and the incompressibility condition for the velocity field **u**

$$E (\partial_t - \nabla^2) u = E u \wedge (\nabla \wedge u) + RaTr - \hat{z} \wedge u - \nabla \pi, \qquad (6)$$

$$\nabla \cdot u = 0 \qquad (7)$$

with $\hat{z}$ being the rotation axis. The evolution of the temperature T is described by the non-dimensional transport equation

$$(P_r \partial_t - \nabla^2) T = S - P_r u \cdot \nabla T, \qquad (8)$$

and the non-dimensional basic state is given by

$$T_b = \frac{1}{2} (1 - r^2). \qquad (9)$$

The system is subject to stress-free and impenetrability mechanical boundary conditions and a fixed temperature at the outer boundary. Thus, while the radial velocity component has to vanish, a non-zero horizontal velocity component is possible at the boundary.

The benchmark solution is obtained for an Ekman number $E = 3 \times 10^{-4}$, a Prandtl number $Pr = 1$, a Rayleigh number $Ra = 95$ and a source term $S = 3$. This choice of parameters is close to the critical values for the onset of convection. More than one solution exists for this choice of parameters. To select the right solution branch, the following initial condition should be used for the temperature field:

$$T = \frac{1}{2} (1 - r^2) + \frac{10^{-5}}{8} \sqrt{\frac{35}{\pi}} r^3 (1 - r^2) (\cos 3\phi + \sin 3\phi) \sin^3 \theta. \qquad (10)$$

For the sake of completeness, the second solution branch might be selected by replacing the spherical harmonic perturbation of degree and order 3 by a spherical harmonic perturbation of degree and order 4. The velocity field can safely be initialized to zero

$$u = 0. \qquad (11)$$

After the initial transient, the solution to Benchmark 1 settles in a quasi-stationary solution with a threefold symmetry. The alternative branch would lead to a similar solution with fourfold symmetry. To illustrate the solution, a few equatorial and meridional slices of the velocity field are provided in Figs 1 and 2.

Once the stable regime is reached, the solution exhibits a constant kinetic energy

$$E_k = \frac{1}{2} \int_V \mathbf{u}^2 dV \qquad (12)$$

providing an ideal diagnostic for the comparison among the different submissions (Fig. 3).

Furthermore, the whole solution is slowly drifting at a constant drift frequency. Similarly to what was done in B1, the velocity field of the solution can be described as

$$u = u(r, \theta, \varphi - 2\pi f_d t), \qquad (13)$$

where $f_d$ is the drift frequency in units of s⁻¹. The drift frequency $f_d$ is related to the angular velocity or drift rate $\omega_d$ in units of rad s⁻¹ by

$$\omega_d = 2\pi f_d. \qquad (14)$$

With this choice of definition, due to the threefold spatial symmetry (see Fig. 1a), the drift frequency $f_d$ represents the frequency at which a given flow pattern will pass through a fixed point in space.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Equatorial slices of (a) the radial component $u_r$ and (b) the azimuthal component $u_\varphi$ of the velocity field for Benchmark 1. The velocity field is equatorially antisymmetric and thus the latitudinal component $u_\theta$ is zero in the equatorial plane.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>Figure 2. Meridional slices of the velocity field **u** for Benchmark 1: (a) radial component $u_r$, (b) latitudinal component $u_\theta$ and (c) azimuthal component $u_\phi$. The slices are chosen such that they contain the maximal amplitude of the velocity field $|\mathbf{u}|$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3. Typical time evolution of the kinetic energy $E_k$ for Benchmark 1. After the initial transient, the kinetic energy reaches a constant value.</figcaption></figure>

The whole solution pattern completes a full rotation at a frequency $\tilde{f}_d = f_d/3$.

To compare the results of the six participants in Benchmark 1, the constant kinetic energy $E_k$ and the drift frequency $f_d$ of their solution were requested.

## 2.2 Benchmark 2: thermally driven dynamo

Benchmark 2 extends the first benchmark by incorporating the generation and evolution of magnetic fields. While still working within the frame of the Boussinesq approximation, the sphere is now filled with a conducting fluid with magnetic diffusivity $\eta$ and magnetic permeability $\mu$. It is still thermally forced through the presence of a homogeneous distribution of heat sources resulting in the basic state given by eq. (2). The system of equations is extended by the induction equation to describe the evolution of the magnetic field **B**. The equations are non-dimensionalized using the radius $r_o$ as length scale, the magnetic diffusion time $r_o^2/\eta$ as timescale, $\beta r_o^2$ as the temperature scale and the magnetic field is rescaled by $\sqrt{2\Omega\rho\mu\eta}$. The four required parameters are the Ekman number $E$

$$E = \frac{\nu}{2\Omega r_o^2}, \qquad (15)$$

the magnetic Rossby number (also referred to as the magnetic Ekman number) $Ro$

$$Ro = \frac{E}{Pm} = \frac{\eta}{2\Omega r_o^2}, \qquad (16)$$

the Roberts number $q$

$$q = \frac{Pm}{Pr} = \frac{\kappa}{\eta}, \qquad (17)$$

with $\kappa$ the thermal diffusivity and the Rayleigh number $Ra$

$$Ra = \frac{g\alpha\beta r_o^3}{2\Omega\kappa}, \qquad (18)$$

with $\alpha$ the thermal expansion coefficient. To ease the conversion to different choices of non-dimensionalizations, the Prandtl number $Pr$

$$Pr = \frac{\nu}{\kappa} = \frac{Pm}{q} = \frac{E}{qRo}, \qquad (19)$$

and the magnetic Prandtl number $Pm$

$$Pm = \frac{\nu}{\eta} = \frac{E}{Ro}, \qquad (20)$$

are also introduced. The motion of the conducting fluid is described by the non-dimensional Navier-Stokes equation and the incompressibility condition for the velocity field **u**

$$(Ro\partial_t - E\nabla^2)\mathbf{u} = Ro\mathbf{u} \wedge (\nabla \wedge \mathbf{u}) + (\nabla \wedge \mathbf{B}) \wedge \mathbf{B} + qRa\mathbf{Tr} - \hat{\mathbf{z}} \wedge \mathbf{u} - \nabla\pi, \qquad (21)$$

$$\nabla \cdot \mathbf{u} = 0. \qquad (22)$$

The magnetic induction equation and the solenoidal condition for the magnetic field **B** read as

$$(\partial_t - \nabla^2) \mathbf{B} = \nabla \wedge (\mathbf{u} \wedge \mathbf{B}), \qquad (23)$$

$$\nabla \cdot \mathbf{B} = 0, \qquad (24)$$

and finally the transport equation for the temperature $T$ is given by

$$(\partial_t - q\nabla^2) T = S - \mathbf{u} \cdot \nabla T. \qquad (25)$$

As for Benchmark 1, the outer boundary is maintained at fixed temperature and a stress-free mechanical boundary condition is imposed. Furthermore, the outer region is considered to be a perfect insulator.

The thermal dynamo solution for Benchmark 2 is obtained for an Ekman number $E = 5 \times 10^{-4}$, a magnetic Rossby number $Ro = \frac{5}{7} \times 10^{-4}$, a Roberts number $q = 7$, a Rayleigh number $Ra = 200$ and a source term $S = 3q = 21$. In terms of the Prandtl numbers, the Benchmark 2 is obtained for a Prandtl number $Pr = 1$ and a magnetic Prandtl $Pm = 7$. This parameter regime is approximately two times supercritical and a magnetic field is generated and sustained by the system. Although the solution for Benchmark 2 can be obtained by starting from a small initial perturbation, the convergence to the final state is extremely slow and requires prohibitively high computational resources. Furthermore, the presence
of several solution branches can not be excluded even if it was not seen during the computations. To reduce the computational load to a reasonable level, a special initial condition exhibiting a much faster convergence has been worked out.

The temperature field $T$ should be initialized with the background conducting state with a small perturbation as a spherical harmonic of degree and order 3

$$T = \frac{1}{2}(1-r^2) + \frac{\epsilon}{8}\sqrt{\frac{35}{\pi}}r^3(1-r^2)(\cos 3\varphi + \sin 3\varphi)\sin^3\theta \quad (26)$$

with $\epsilon = 10^{-5}$. The magnetic field should be initialized with a purely toroidal magnetic field given by

$$B_r = 0, \quad (27)$$

$$B_{\theta} = -\frac{3}{2}r(-1 + 4r^2 - 6r^4 + 3r^6)(\cos \varphi + \sin \varphi), \quad (28)$$

$$B_{\varphi} = -\frac{3}{4}r(-1 + r^2)\cos\theta [3r(2 - 5r^2 + 4r^4)\sin\theta + 2(1 - 3r^2 + 3r^4)(\cos\varphi - \sin\varphi)]. \quad (29)$$

Finally, the velocity field should be initialized with a purely toroidal velocity given by

$$u_r = 0, \quad (30)$$

$$u_{\theta} = -\frac{10r^2}{7\sqrt{3}}\cos\theta[3(-147 + 343r^2 - 217r^4 + 29r^6)\cos\varphi + 14(-9 - 125r^2 + 39r^4 + 27r^6)\sin\varphi], \quad (31)$$

$$u_{\varphi} = -\frac{5r}{5544}\{7[(43\,700 - 58\,113r^2 - 15\,345r^4 + 1881r^6 + 20\,790r^8)\sin\theta + 1485r^2(-9 + 115r^2 - 167r^4 + 70r^6)\sin 3\theta] + 528\sqrt{3}r\cos 2\theta[14(-9 - 125r^2 + 39r^4 + 27r^6)\cos\varphi + 3(147 - 343r^2 + 217r^4 - 29r^6)\sin\varphi]\}. \quad (32)$$

For simulations using spherical harmonics and a toroidal/poloidal decomposition, the expression of the initial conditions can be written in a simpler form. Assuming that the magnetic field is decomposed as $B = \nabla \wedge Tr + \nabla \wedge \nabla \wedge Pr$, the initial field is given by the two scalars

$$\begin{aligned} T = r \left[ \left( \frac{3}{4} - 3r^2 + \frac{9}{2}r^4 - \frac{9}{4}r^6 \right) + \left( \frac{3}{4} - 3r^2 + \frac{9}{2}r^4 - \frac{9}{4}r^6 \right) t \right] Y_1^1 \\ + r^2 \left( \frac{3}{2} - \frac{21}{4}r^2 + \frac{27}{4}r^4 - 3r^6 \right) Y_2^0 + \text{c.c.}, \end{aligned} \quad (33)$$

$$P = 0, \quad (34)$$

where c.c. stands for the complex conjugate without the $m=0$ modes and $Y_1^m$ are Schmidt normalized complex spherical harmonics. Similarly assuming that the velocity field is decomposed as $u = \nabla \wedge Tr + \nabla \wedge \nabla \wedge Pr$, the initial condition is given by

$$\begin{aligned} T = r^2 \left[ \left( 30 + \frac{1250}{3}r^2 - 130r^4 - 90r^6 \right) \right. \\ \left. + \left( 105 - 245r^2 + 155r^4 - \frac{145}{7}r^6 \right) t \right] Y_2^1 \\ + r \left( -\frac{54625}{198} + 350r^2 + \frac{625}{2}r^4 - 325r^6 \right) Y_1^0 \\ + r^3 (45 - 575r^2 + 835r^4 - 350r^6) Y_3^0 + \text{c.c.}, \end{aligned} \quad (35)$$

$$P = 0. \quad (36)$$

Details of the definition and the normalization of the spherical harmonics $Y_1^m$ are given in Appendix A.

The structure of the solution to Benchmark 2 is much more complicated than for Benchmark 1 and no longer exhibits a simple symmetry. Nevertheless, after the initial transient the system settles into a regime with periodic kinetic and magnetic energies. The amplitude and frequency of these oscillations provide a good diagnostic for Benchmark 2 (Fig. 4).

To illustrate the structure of the velocity field, a few equatorial slices are shown in Fig. 5. The time-dependent features of the solutions can be seen in the Hammer projection snapshots of the flow close to the outer boundary (Fig. 6) and the Hammer projections of the radial component of the magnetic field at the outer boundary (Fig. 7). The two comma shaped flux patches of opposite sign

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Kinetic energy $E_k$ and magnetic energy $E_m$ for Benchmark 2. (a) Typical time evolution of $E_k$ and $E_m$. After the initial transient, both energies settle into a periodic oscillation. (b) Detailed view of the oscillations in $E_k$ and $E_m$. The magnetic energy has been rescaled by $\xi = \frac{c_k}{c_m} \approx 39$ (see eqs 37 and 40) to show the phase shift between $E_k$ and $E_m$.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>Figure 5. Equatorial slices of the velocity field $u$ at $t = 4.43241$ for Benchmark 2: (a) radial component $u_r$, (b) latitudinal component $u_\theta$ and (c) azimuthal component $u_\psi$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 6. Hammer projections of three snapshots of the azimuthal velocity component $u_\psi$ at the outer boundary ($r = r_o$) for Benchmark 2.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>Figure 7. Hammer projections of six snapshots, spanning approximately two periods, of the radial magnetic field $B_r$ at the outer boundary $r = r_o$ for Benchmark 2.</figcaption></figure>

emerge periodically. The lower part moves eastwards while rising to higher latitudes until they vanish. Fig 7 shows snapshots over approximately two periods.

The solution to Benchmark 2 will solely be characterized by the computed kinetic and magnetic energies. Their periodic behaviour allows us to define several diagnostic quantities used to compare the solutions from different simulations. The kinetic energy is decomposed into a constant component $C_k$, the amplitude of the leading oscillating term $A_k$, the frequency of this oscillation $f_k$ and a phase shift $\zeta_k$. Using these definitions, the kinetic energy $E_k$ is written as

$$E_k = \frac{1}{2} \int_V \boldsymbol{u}^2 dV = C_k + A_k \cos(2\pi f_k t + \zeta_k) + \dots \quad (37)$$

Furthermore, decomposing the velocity into its equatorially symmetric [labelled as (s)] and equatorially antisymmetric part [labelled as (a)], the velocity field was found to be purely equatorially

symmetric

$$E_k^{(s)} = E_k, \quad (38)$$

$$E_k^{(a)} = 0. \quad (39)$$

Using the same decomposition, the magnetic energy $E_m$ can be written as

$$E_m = \frac{1}{2Ro} \int_V \boldsymbol{B}^2 dV = C_m + A_m \cos(2\pi f_m t + \zeta_m) + \dots \quad (40)$$

Decomposing the magnetic field into its equatorially symmetric and antisymmetric parts, the magnetic field is found to be purely equatorially antisymmetric

$$E_m^{(s)} = 0, \quad (41)$$

$$E_m^{(a)} = E_m. \quad (42)$$
<figure><img src="image_10.png" /><figcaption>**Figure 8.** (a) Frequency spectrum of the kinetic energy $E_k$ and magnetic energy $E_m$ for Benchmark 2 after a flat top window has been applied on the time-series. (b) Details of the peak of the main amplitude of the kinetic and magnetic energies after a Kaiser window with parameter $\beta = 14$ has been applied to the time-series. Note the use of a logarithmic scale for the ordinate, meaning that the frequency localization is extremely good. Both spectra have been obtained on a time-series over approximately 0.35 diffusion times at a sampling rate of $1.75 \times 10^5$ obtained by (MJ) at $N = 31$ and $L = M = 63$ (see Section 3 for details).</figcaption></figure>

It was found and confirmed by all participants that the oscillations in the magnetic energy and kinetic energy have the same frequency. In addition, as can be seen in Fig. 4(b), there is a phase shift between the magnetic and kinetic energy. This relative phase shift ($\zeta_k - \zeta_m$) is found to have value 1.91rad but is not included in the benchmark. The six constants $C_k$, $A_k$, $f_k$, $C_m$, $A_m$ and $f_m$ defined above are the diagnostic values that were requested from all submissions to Benchmark 2.

The computation of decompositions 37 and 40 requires some further explanations as several different approaches are possible. The first one involves computing the Fourier series. To reduce the spectral leakage, a flat top window (Oliphant 2007) is applied on the time-series. The amplitude of the different modes can then easily be extracted as shown in Fig. 8(a). Conversely, to compute an accurate frequency, a Kaiser window (Kaiser & Kuo 1966) with parameter $\beta = 14$ is applied to the time-series. The peaks in the spectrum are then well approximated with a parabolic fit allowing for interpolation between the available frequencies. The peak for the main oscillation is shown in Fig. 8(b). The frequencies can also be computed by counting the zero crossings. Both approaches provide the frequency within a relative error of $10^{-4}$ per cent. There is a simpler approach to obtain the requested data without using a Fourier transform. Assuming $E_i^{\min}$, $E_i^{\max}$ are the minimum and maximum of the time-series for the energy, the constant component is approximately given by

$$C_i \approx \frac{E_i^{\min} + E_i^{\max}}{2}, \qquad (43)$$

and the amplitude of the main oscillation is approximately given by

$$A_i \approx \frac{E_i^{\max} - E_i^{\min}}{2}. \qquad (44)$$

These relations are not exact because the time-series do also include higher frequencies as shown in Fig. 8(a). Using the Fourier series, it is possible to bound the relative errors generated by the simplified approach. By including the second peak with a frequency of $2f_1$ and accounting for the phase shift, a comparison with the approximations 43 and 44 provides the relative errors $\epsilon_{C_k}$ and $\epsilon_{A_k}$ on $C_k$ and $A_k$ for the kinetic energy

$$\epsilon_{C_k} = 3 \times 10^{-2} \text{ per cent}, \qquad (45)$$

$$\epsilon_{A_k} = 5 \times 10^{-2} \text{ per cent}. \qquad (46)$$

The same analysis on the magnetic energy provides the relative errors $\epsilon_{C_m}$ and $\epsilon_{A_m}$ on $C_m$ and $A_m$

$$\epsilon_{C_m} = 0.4 \text{ per cent}, \qquad (47)$$

$$\epsilon_{A_m} = 0.6 \text{ per cent}. \qquad (48)$$

These relative errors have been obtained by assuming that $E_i^{\min}$ and $E_i^{\max}$ are exact. An accurate value for the minimum and maximum requires a time-series sampled at a high frequency but it does in principle only require one cycle. On the other hand, an accurate Fourier series will require a long time-series. The errors made by using this simplified approach depend strongly on the relative magnitude of the higher frequencies. For the kinetic energy, the second oscillation with a frequency of $2f_1$ has an approximately 63 times smaller amplitude which leads to the rather small errors given in eqs (45) and (46). On the other hand, the second oscillation in the magnetic energy is only approximately 10 times smaller leading to the larger errors given in eqs (47) and (48). In both cases, the third oscillation with a frequency $3f_1$ can be neglected as it is more than 20 times smaller than the second oscillation.

## 2.3 Benchmark 3: boundary forced rotating bubble

Benchmark 3 is again a purely hydrodynamic problem. It provides a simple test problem for boundary driven flows in a whole sphere as they might, for example, arise in precession or libration problems. The proposed system describes the motion of an incompressible fluid inside a spherical bubble rising in a rotating fluid. It is an important addition to the first two cases as it replaces the internal thermal forcing by a mechanical forcing due to an imposed tangential flow over the boundary. The bubble is assumed to be of unit radius $r_o = 1$ and is described by the Navier-Stokes equation and incompressibility condition

$$\partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u} + 2\Omega \hat{z} \wedge \mathbf{u} = -\nabla p + \nu \nabla^2 \mathbf{u}, \qquad (49)$$

$$\nabla \cdot \mathbf{u} = 0 \qquad (50)$$

with $\nu$ the viscosity of the fluid and $\Omega$ the rotation rate of the bubble. The rotation axis is parallel to $\hat{z}$. The tangential flow over the boundary imposes a non-homogeneous boundary condition on the fluid at the surface of the bubble

$$u_\theta = -u_0 \cos \theta \cos \phi, \qquad (51)$$

$$u_\phi = u_0 \sin \phi, \qquad (52)$$
<figure><img src="image_11.png" /><figcaption>**Figure 9.** In the top row, velocity field in the x-z plane ($\phi = 0$) and in the bottom row, velocity field in the y-z plane ($\phi = \pi/2$) for Benchmark 3.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Velocity field in the x-y (equatorial) plane for Benchmark 3.</figcaption></figure>

which is the gradient of a pure $l = 1$, $m = 1$ spherical harmonic $Y_l^m(\theta, \phi)$.

The solution to Benchmark 3 is obtained for a boundary velocity $u_0 = \sqrt{3/2\pi}$, a viscosity $\nu = 10^{-2}$ and a rotation rate $\Omega = 10$. It does not require any particular initial condition and can be started with a zero initial velocity field.

The flow converges quickly to a stationary solution with a dominant $m = 1$ component. Interestingly, the solution contains an im-

portant and non-trivial flow through the centre of the bubble. As such, it provides a good diagnostic for this numerically challenging region. To illustrate the flow of the solution, the velocity field components in the x-z and y-z planes are shown in Fig. 9. The velocity field in the x-y plane, which is orthogonal to the rotation axis (i.e. the equatorial plane) is shown in Fig. 10.

The solution to Benchmark 3 was characterized by five diagnostic values. The first diagnostic is the constant kinetic energy
$$E_k = \frac{1}{2} \int_V \mathbf{u}^2 dV \qquad (53)$$

reached after the initial transient. When available, the energy in the spherical harmonic orders $m = 0, m = 1$ and $m = 2$ have also been collected. Secondly, the $\hat{z}$ component of the angular momentum $L_z$ is reported. The three components of the velocity field $\mathbf{u}$ at the centre of the bubble provide that last diagnostic data.

### 3 CONTRIBUTING NUMERICAL CODES

There was a total of nine contributors to the benchmark but each did not necessarily provide results for all three cases. A short description of the algorithm used by each simulation is provided below. Further references are given for a more detailed description.

Marti and Jackson (MJ): Spectral simulation using spherical harmonics for the angular component and polynomials developed by Worland (2004) and Livermore et al. (2007) in radius (see Marti 2012). The Worland polynomials satisfy exactly the parity and regularity conditions required at the origin of the spherical coordinate system. Specifically, the radial basis used is of the form $r^l P_n^{(-1/2, l-1/2)}(2r^2 - 1)$ for a given harmonic degree $l$, with $P_n^{(\alpha, \beta)}(x)$ the Jacobi polynomials. The incompressibility condition is guaranteed by the use of a toroidal/poloidal decomposition of the vector fields. A second-order predictor-corrector scheme is used for the time integration. No specific treatment is required to conserve angular momentum.

Hollerbach (H): An adaptation of the previous spherical shell code described by Hollerbach (2000), based on spherical harmonics and a toroidal/poloidal decomposition of the vector fields. Instead of expanding in the full set of Chebyshev polynomials in radius, regularity and parity conditions at the origin are now accommodated by expanding as $r \cdot T_{2k-1}(r)$ for odd harmonic degrees and $r^2 \cdot T_{2k-1}(r)$ for even harmonic degrees for the toroidal/poloidal scalars. The temperature is similarly expanded as $T_{2k-2}(r)$ for degree $l = 0, r \cdot T_{2k-2}(r)$ for odd degrees, and $r^2 \cdot T_{2k-2}(r)$ for even degrees. Angular momentum conservation is explicitly imposed by a modified stress-free boundary condition.

Aubert (A): Spectral simulation using the code PARODY-JA, which uses spherical harmonics for the angular component and a second-order finite differences in radius (see Dormy et al. 1998; Aubert et al. 2008). The time marching is done with a semi-implicit Crank–Nicolson/Adams–Bashforth scheme. The radial mesh includes a gridpoint exactly at the centre of the sphere. Spectral toroidal and poloidal components of order $l$ behave like $r^l$ at the centre. Angular momentum conservation is achieved by correcting for solid-body rotation at each time step.

Schaeffer (S): Spectral simulation using spherical harmonics for the angular component and second-order finite differences in radius (Monteux et al. 2012). The numerical instability near the origin is overcome by truncating the spherical harmonic expansion at $\ell_{tr}(r)$ before computing the spatial fields that enter the non-linear terms. Specifically, the truncation is $\ell_{tr}(r) = 1 + (\ell_{max} - 1) \left(\frac{r}{r_0}\right)^\alpha$, where $\alpha = 0.5$ gives good results, and also saves some computation time. The time stepping uses a semi-implicit Crank–Nicolson scheme for the diffusive terms, while the non-linear terms can be handled either by an Adams–Bashforth or a predictor-corrector scheme (both second-order in time). The SHTns library (Schaeffer 2013) is used for efficient spherical harmonic transforms. Angular momentum conservation is achieved by adjusting the solid-body rotation component at each time step.

Takehiro, Sasaki and Hayashi (TSH): Spectral simulation using spherical harmonics for the angular components and the polynomials developed by Matsushima & Marcus (1995) and Boyd (2001) in radius (see Sasaki et al. 2012). The radial basis functions satisfy exactly the parity and regularity conditions at the origin of the spherical coordinate system. Specifically, the used radial basis is of the form $r^l P_n^{(\alpha, \beta)}(2r^2 - 1)$ for a given harmonic degree $l$, with $P_n(x)$ the Jacobi polynomials. The incompressibility condition is guaranteed by the use of a toroidal/poloidal decomposition of the vector fields. The time integration is performed with the Crank–Nicolson scheme for the diffusive terms and a second-order Adams–Bashforth scheme for the other terms. No specific treatment is required to conserve angular momentum.

Simitev and Busse (SB): Pseudospectral numerical code using spherical harmonics expansion in the angular variables and Chebyshev polynomials in radius. Time stepping is implemented by a combination of the implicit Crank–Nicolson scheme for the diffusion terms and the explicit Adams–Bashforth scheme for the Coriolis and the non-linear terms; both schemes are second-order accurate. Early versions of the code are described in Tilgner & Busse (1997) and Tilgner (1999). The code has been extensively modified and used for a number of years (Simitev & Busse 2005, 2009, 2012; Busse & Simitev 2006, 2008). This is a spherical shell code and no special effort was made to convert it to the full sphere geometry. Instead, the full sphere is approximated by placing a tiny inner core with radius ratio $r_i/r_o = 0.01$ at the centre of the shell. Angular momentum conservation is achieved by correcting for rigid-body rotation if required.

Cébron (C): Finite elements method simulation using the standard Lagrange element P2–P3, which is quadratic for the pressure field and cubic for the velocity field, and a Galerkin Least-Squares (GLS) stabilization method (Hauke & Hughes 1994). The (unstructured) mesh is made of prisms in the boundary layer and tetrahedrons in the bulk. The incompressibility is imposed using a penalty method. The time stepping uses the Implicit Differential-Algebraic solver (IDA solver), based on variable-coefficient Backward Differenting Formulae (e.g. Hindmarsh et al. 2005). The integration method in IDA is variable-order, the order ranging between 1 and 5. At each time step, the system is solved with the sparse direct linear solver PARDISO (www.pardiso-project.org) or a multigrid GMRES iterative solver. This is all implemented via the commercial code COMSOL Multiphysics®.

Nore, Luddens and Guermond (NLG): Hybrid Fourier and finite element method using a Fourier decomposition in the azimuthal direction and the standard Lagrange elements P1–P2 in the meridian section (with P1 for the pressure and P2 for the velocity field). The meridian mesh is made of quadratic triangles. The velocity and pressure are decoupled by using the rotational pressure-correction method. The time stepping uses the second-order Backward Difference Formula (BDF2). The non-linear terms are made explicit and approximated using second-order extrapolation in time. The code is parallelized in Fourier space and in meridian sections [domain decomposition with METIS (Karypis & Kumar 2009)] using MPI and PETSC (Portable, Extensible Toolkit for Scientific Computation; Balay et al. 1997, 2012a,b). This is implemented in the code SFEMaNS (for Spectral/Finite Element method for Maxwell and Navier–Stokes equations; Guermond et al. 2007, 2009, 2011).

Vantieghem (V): Unstructured finite-volume simulation (see Vantieghem 2011) using a grid of tetrahedral elements with smaller elements close to the wall. The spatial discretization is based on a centred-difference-like stencil that is second-order accurate for
regular tetrahedra. Time stepping is based on a canonical fractional-step method (Kim & Moin 1985), and the equations are integrated in time with a fourth order Runge-Kutta method. A BiCGstab(2) algorithm is used to solve the pressure Poisson equation. The reported Fourier components are obtained by an *a posteriori* interpolation of the results on a regular grid in terms of spherical coordinates ($N_r = 36, N_\theta = N_\phi = 18$), which is subject to considerable additional numerical (interpolation) errors.

<table><caption>**Table 1.** Contributions to Benchmark 1. The labels used for the different codes are defined in Section 3. The values are shown with the number of significant digits provided by the authors. As all these codes are based on a spherical harmonics expansion for the angular component, the resolution is given as the radial resolution *N*, the highest harmonic degree *L* and the highest harmonic order *M*.</caption><thead><tr><th>Code</th><th><i>E</i><sub>k</sub></th><th><i>f</i><sub>d</sub></th><th><i>N</i></th><th><i>L</i></th><th><i>M</i></th></tr></thead><tbody><tr><td>(MJ)</td><td>29.08502</td><td>12.38841</td><td>8</td><td>15</td><td>15</td></tr><tr><td>(MJ)</td><td>29.07661</td><td>12.38860</td><td>8</td><td>23</td><td>23</td></tr><tr><td>(MJ)</td><td>29.12178</td><td>12.38604</td><td>12</td><td>23</td><td>23</td></tr><tr><td>(MJ)</td><td>29.12064</td><td>12.38619</td><td>16</td><td>23</td><td>23</td></tr><tr><td>(MJ)</td><td>29.12064</td><td>12.38619</td><td>16</td><td>31</td><td>31</td></tr><tr><td>(MJ)</td><td>29.12068</td><td>12.38619</td><td>23</td><td>47</td><td>47</td></tr><tr><td>(MJ)</td><td>29.12068</td><td>12.38619</td><td>31</td><td>63</td><td>63</td></tr><tr><td>(H)</td><td>29.11784</td><td>12.3862</td><td>12</td><td>23</td><td>23</td></tr><tr><td>(H)</td><td>29.12054</td><td>12.3862</td><td>16</td><td>31</td><td>31</td></tr><tr><td>(H)</td><td>29.12053</td><td>12.3862</td><td>23</td><td>31</td><td>31</td></tr><tr><td>(H)</td><td>29.12053</td><td>12.3862</td><td>23</td><td>47</td><td>47</td></tr><tr><td>(H)</td><td>29.12053</td><td>12.3862</td><td>31</td><td>63</td><td>63</td></tr><tr><td>(S)</td><td>29.219</td><td>12.388</td><td>120</td><td>31</td><td>31</td></tr><tr><td>(S)</td><td>29.1446</td><td>12.387</td><td>250</td><td>63</td><td>63</td></tr><tr><td>(S)</td><td>29.13501</td><td>12.38648</td><td>320</td><td>85</td><td>85</td></tr><tr><td>(TSH)</td><td>29.03074</td><td>12.3878</td><td>16</td><td>21</td><td>21</td></tr><tr><td>(TSH)</td><td>29.12878</td><td>12.3863</td><td>32</td><td>42</td><td>42</td></tr><tr><td>(TSH)</td><td>29.12878</td><td>12.3863</td><td>48</td><td>85</td><td>85</td></tr><tr><td>(SB)</td><td>29.00617</td><td>11.89445</td><td>33</td><td>42</td><td>42</td></tr><tr><td>(A)</td><td>29.12062</td><td>12.3931</td><td>1600</td><td>63</td><td>63</td></tr></tbody></table>

# 4 RESULTS

There is quite an important diversity in the type of simulations that took part in these benchmarks. All the diagnostics that have been considered for these benchmarks should be straightforward to obtain whenever the simulation code is based on some spectral expansion or on a local method. On the other hand, a direct comparison of the resolution used is a more subtle problem. The comparison will be done by comparing solutions based on the number of degrees of freedom present at the time stepping level. For local methods, the resolution *R* is computed as $R = N_{\text{grid}}^{1/3}$ where $N_{\text{grid}}$ is the number of gridpoints and for spherical harmonic-based codes $R = \{N_r \cdot [L_{\text{max}}(2M_{\text{max}} + 1) - M_{\text{max}}^2 + M_{\text{max}} + 1]\}^{1/3}$. The same approach was used in B1.

## 4.1 Benchmark 1

There were six participants in Benchmark 1 and all of them used a spherical harmonics-based simulation. They agree qualitatively quite well and no important discrepancies were found. The details of all the solutions obtained by the participants is given in Table 1. At the quantitative level, a few interesting observations can be made. The results for the total kinetic energy are summarized in Fig. 11(a). The five codes (MJ), (H), (TSH), (S), (A) do all eventually converge to the same solution within $5 \times 10^{-2}$ per cent. While using a completely different radial expansion, (MJ) and (H) even converge very rapidly within $5 \times 10^{-4}$ per cent. The last code (SB) working in a spherical shell rather than a sphere comes within 0.4 per cent. Note that the results obtained with a very high radial resolution (1600 gridpoints) by (A) matches very closely to the solutions from the fully spectral codes. The other finite differences-based code (S) shows a clear convergence towards the same solution and would most likely have reached it at a higher resolution.

The picture is very similar for the drift frequency (Fig. 11b). The codes (MJ), (H) and (TSH) agree within $6 \times 10^{-4}$ per cent, the solution by (S) is within $2 \times 10^{-3}$ per cent and the solution by (A) is within $6 \times 10^{-2}$ per cent. Finally, the solution by (SB) is within 4 per cent.

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Summary of the solutions of the participants of Benchmark 1. (a) Constant kinetic energy $E_k$, as defined in eq. (12), reached after the initial transient. (b) Drift frequency $f_d$, defined in Section 2, of the threefold symmetric structure of the solution (see Fig. 1). The black horizontal line is the standard value given in Table 2. The error corridor for the standard values is represented by a greyed out area, but as the errors are very small it is only barely visible. The labels used for the different codes are defined in Section 3.</figcaption></figure>
<table><caption>**Table 2.** Summary table and standard values obtained for Benchmark 1. The Ekman number *E*, the Prandtl number *Pr* and the modified Rayleigh number *Ra* as well as the governing equation for the velocity *u* and the temperature *T* are given in Section 2. The kinetic energy *E*<sub>k</sub> is defined in eq. (12) and the drift frequency by eq. (13).</caption><caption>Benchmark 1: Thermal convection</caption><thead><tr><th colspan="3">Parameters</th></tr><tr><th><i>E</i><br>3 × 10<sup>-4</sup></th><th><i>Pr</i><br>1</th><th><i>Ra</i><br>95</th></tr></thead><tbody><tr><td colspan="3"><strong>Boundary conditions</strong></td></tr><tr><td><b><i>u</i>:</b></td><td colspan="2">stress-free</td></tr><tr><td><b><i>T</i>:</b></td><td colspan="2">fixed temperature<br>homogeneous heat sources</td></tr><tr><td colspan="3"><strong>Requested values</strong></td></tr><tr><td><i>E</i><sub>k</sub></td><td>29.1206</td><td>±1 × 10<sup>-4</sup></td></tr><tr><td><i>f</i><sub>d</sub></td><td>12.3862</td><td>±1 × 10<sup>-4</sup></td></tr></tbody></table>

A standard value for the kinetic energy and drift frequency for Benchmark 1 is derived from the three simulations showing the best convergence to a common result. The values as well as a short summary of the problem definition are given in Table 2. The very high convergence to a common solution allow to provide these with roughly six significant digits.

## 4.2 Benchmark 2

All the groups participating in Benchmark 1 also participated in Benchmark 2, except for (A). There are thus again only solutions from spherical harmonics-based simulations. The overall picture is very similar to Benchmark 1. The fully spectral simulations converge very rapidly to their final solution and which are in most cases very close to each other. The finite differences simulation, while requiring higher radial resolution, eventually reaches a similar solution. As was already observed in Benchmark 1, the simulation by (SB) which approximates the whole sphere as a shell with a

tiny inner core, shows the largest discrepancy. The details of all the solutions to Benchmark 2 are given in Table 3.

The solutions obtained by the different simulations are summarized in Figs 12(a)–(e). For the constant part of the kinetic ($C_k$) and magnetic ($C_m$) energies, excluding the solution by (SB), all results lie within $2 \times 10^{-2}$ and 0.3 per cent, respectively. A good convergence to a common solution was also observed for the amplitude of the oscillations in the kinetic energy. The solutions provided by (MJ), (H) and (S) lie within $7 \times 10^{-3}$ per cent. The solution of (TSH) does not seem to converge to exactly the same value but remains within 0.4 per cent of the three other values. The convergence is not as good for the amplitude in the oscillation of the magnetic energy. All the results lie within 1 per cent. The solution by (MJ), (H) and (TSH) seem to show the clearest convergence trend and their solutions lie within 0.2 per cent. Finally, the frequency of the oscillations of the kinetic and magnetic energies have been compared. All groups reported the same frequency *f* for both energies. The summary of the solutions for the frequency is shown in Fig. 12(e). The results for (MJ), (H), (TSH) and (S) lie within 0.3 per cent while (SB) is a little bit further away with 0.6 per cent.

While the spread in the solutions is clearly more important for Benchmark 2, the general situation is very similar to Benchmark 1. The solutions by (MJ) and (H) do nearly overlap for all five diagnostic values and do exhibit a very fast convergence. As was explained in Section 3, the extraction of the different components of the energies requires some post-processing. The choice of methodology by each author, for example, to extract the oscillation amplitude, might explain a part of the somewhat larger discrepancies compared to Benchmark 1. The standard values given in Table 4 are obtained by taking the average of the highest resolution by (MJ) and (H). The error bars are chosen such that at least one additional solution is included in the error corridor. This choice is based on the fast convergence of both codes to essentially the same value for all the requested data.

While it was not part of the actual benchmark, the phase shift between the kinetic and magnetic energy (see Fig. 4b) is also reported in Table 4 to provide a more complete characterization of the solution. The reported value has been computed by (MJ) from

<table><caption>**Table 3.** Spectral method contributions to Benchmark 2. The labels used for the different codes are defined in Section 3. The values are shown with the number of significant digits provided by the authors. The decomposition of the kinetic energy *E*<sub>k</sub> into *C*<sub>k</sub>, *A*<sub>k</sub> and *f*<sub>k</sub> is defined in eq. (37) and the equivalent decomposition of the magnetic energy *E*<sub>m</sub> into *C*<sub>m</sub>, *A*<sub>m</sub> and *f*<sub>m</sub> is defined in eq. (40). As all these codes are based on a spherical harmonic expansion for the angular component, the resolution is given as the radial resolution *N*, the highest harmonic degree *L* and the highest harmonic order *M*.</caption><thead><tr><th>Code</th><th><i>C</i><sub>k</sub></th><th><i>A</i><sub>k</sub></th><th><i>f</i><sub>k</sub></th><th><i>C</i><sub>m</sub></th><th><i>A</i><sub>m</sub></th><th><i>f</i><sub>m</sub></th><th><i>N</i></th><th><i>L</i></th><th><i>M</i></th></tr></thead><tbody><tr><td>(MJ)</td><td>35.141.84</td><td>1836.287</td><td>302.2623</td><td>1153.695</td><td>51.77003</td><td>302.2623</td><td>12</td><td>23</td><td>23</td></tr><tr><td>(MJ)</td><td>35.548.95</td><td>1881.661</td><td>302.6947</td><td>922.3073</td><td>38.54002</td><td>302.6947</td><td>16</td><td>31</td><td>31</td></tr><tr><td>(MJ)</td><td>35.542.15</td><td>1880.460</td><td>302.6858</td><td>924.5757</td><td>38.48190</td><td>302.6858</td><td>23</td><td>31</td><td>31</td></tr><tr><td>(MJ)</td><td>35.551.33</td><td>1880.055</td><td>302.7018</td><td>908.9870</td><td>37.47705</td><td>302.7018</td><td>23</td><td>47</td><td>47</td></tr><tr><td>(MJ)</td><td>35.550.93</td><td>1879.837</td><td>302.7015</td><td>908.8059</td><td>37.45069</td><td>302.7015</td><td>31</td><td>63</td><td>63</td></tr><tr><td>(H)</td><td>35.378</td><td>1855</td><td>302.48</td><td>1043.77</td><td>46.16</td><td>302.48</td><td>12</td><td>23</td><td>23</td></tr><tr><td>(H)</td><td>35.588</td><td>1885</td><td>302.71</td><td>904.30</td><td>37.61</td><td>302.71</td><td>16</td><td>31</td><td>31</td></tr><tr><td>(H)</td><td>35.540</td><td>1881</td><td>302.66</td><td>925.64</td><td>38.50</td><td>302.66</td><td>23</td><td>31</td><td>31</td></tr><tr><td>(H)</td><td>35.551</td><td>1880</td><td>302.67</td><td>909.67</td><td>37.48</td><td>302.67</td><td>23</td><td>47</td><td>47</td></tr><tr><td>(H)</td><td>35.550</td><td>1880</td><td>302.67</td><td>909.46</td><td>37.47</td><td>302.67</td><td>31</td><td>63</td><td>63</td></tr><tr><td>(S)</td><td>35.544</td><td>1878.1</td><td>302.2</td><td>951.59</td><td>41.33</td><td>302.2</td><td>120</td><td>31</td><td>31</td></tr><tr><td>(S)</td><td>35.568</td><td>1881.0</td><td>302.65</td><td>908.69</td><td>37.951</td><td>302.65</td><td>250</td><td>63</td><td>63</td></tr><tr><td>(S)</td><td>35.552</td><td>1880.1</td><td>302.11</td><td>910.75</td><td>38.064</td><td>302.11</td><td>320</td><td>85</td><td>85</td></tr><tr><td>(TSH)</td><td>35.619.63</td><td>1887.200</td><td>303.0303</td><td>881.6272</td><td>36.97913</td><td>303.0303</td><td>32</td><td>42</td><td>42</td></tr><tr><td>(TSH)</td><td>35.564.30</td><td>1872.702</td><td>303.0303</td><td>905.8444</td><td>37.60794</td><td>303.0303</td><td>48</td><td>85</td><td>85</td></tr><tr><td>(SB)</td><td>35.951.5</td><td>1843.38</td><td>304.308</td><td>1046.12</td><td>38.08</td><td>304.308</td><td>41</td><td>96</td><td>96</td></tr></tbody></table>
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Summary of the solutions of the participants to Benchmark 2. The decomposition of the kinetic energy $E_k$ into a constant part $C_k$, the amplitude of the main oscillation $A_k$ and its frequency $f_k$ is defined in eq. (37). The equivalent decomposition for the magnetic energy $E_m$, defining the three constants $C_m$, $A_m$ and $f_m$ is given in eq. (40). The black horizontal line is the standard value given in Table 4. The error corridor for the standard values, defined such that the three most converged solutions are included, is represented by a greyed out area. The labels used for the different codes are defined in Section 3.</figcaption></figure>

time-series of the kinetic and magnetic energy at the highest reported resolution ($N = 31, L = M = 63$). The phase shift has been extracted from the Fourier series shown in Fig. 8.

## 4.3 Benchmark 3

Benchmark 3 had the highest number of participants with eight codes taking part. It is also the only case where results from local methods are available. The computation for Benchmark 3 does not require a very high horizontal resolution. For example, spherical harmonics-based codes exhibit a very good convergence at reso-

lutions as low as $L_{max} = 30$ and $M_{max} = 10$. However, it is more demanding in the radial direction. The centre requires a sufficiently high resolution to describe flow crossing it properly as well as the moving outer boundary which is forcing the system.

The solutions obtained by the different groups are summarized in Figs 13(a)–(e). Note that some of the solutions obtained by spherical harmonics-based codes have been obtained with a reduced longitudinal resolution while a triangular truncation was used for Benchmarks 1 and 2. The details for each solution are given in Table 5 for the spectral methods and in Table 6 for the local methods. A good convergence of the kinetic energy $E_k$ is observed. Surprisingly, it is the solution by (TSH) which shows the largest
<table><caption>**Table 4.** Summary table and standard values obtained for Benchmark 2. The Ekman number *E*, the magnetic Rossby number *Ro*, the Roberts number *q* and the modified Rayleigh number *Ra*, as well as the governing equations for the velocity *u*, the magnetic field *B* and the temperature *T* are given in Section 2. The decomposition of the kinetic energy *E*<sub>k</sub> into constant and oscillating components is defined in eq. (37) and the decomposition of the magnetic energy is defined in eq. (40).</caption><thead><tr><th colspan="4">Benchmark 2: Thermally driven dynamo</th></tr><tr><th colspan="4">Parameters</th></tr><tr><th><i>E</i><br>5 × 10<sup>-4</sup></th><th><i>Ro</i><br>5/7 × 10<sup>-4</sup></th><th><i>q</i><br>7</th><th><i>Ra</i><br>200</th></tr></thead><tbody><tr><td colspan="4">Boundary conditions</td></tr><tr><td><b><i>u</i></b>:</td><td colspan="3">stress-free</td></tr><tr><td><b><i>B</i></b>:</td><td colspan="3">insulating</td></tr><tr><td><b><i>T</i></b>:</td><td colspan="3">fixed temperature<br>homogeneous heat sources</td></tr><tr><td colspan="4">Requested values</td></tr><tr><td colspan="4">Kinetic energy <i>E</i><sub>k</sub></td></tr><tr><td></td><td><i>C</i><sub>k</sub></td><td>35 550.5</td><td>±1.5</td></tr><tr><td></td><td><i>A</i><sub>k</sub></td><td>1879.84</td><td>±0.26</td></tr><tr><td colspan="4">Magnetic energy <i>E</i><sub>m</sub></td></tr><tr><td></td><td><i>C</i><sub>m</sub></td><td>909.133</td><td>±1.62</td></tr><tr><td></td><td><i>A</i><sub>m</sub></td><td>37.4603</td><td>±0.1476</td></tr><tr><td>Frequency <i>f</i> = <i>f</i><sub>k</sub> = <i>f</i><sub>m</sub></td><td></td><td>302.701</td><td>±0.33</td></tr><tr><td colspan="4">Additional characteristic</td></tr><tr><td>Phase shift |<i>ζ</i><sub>k</sub> − <i>ζ</i><sub>m</sub>|</td><td></td><td>1.91rad</td><td></td></tr></tbody></table>

discrepancy (0.4 per cent) while all the other solutions agree within 0.1 per cent at least. (MJ), (H), (S), (C) and (NLG) show the clearest convergence with solutions within $8 \times 10^{-3}$ per cent. The values obtained for the vertical component of the angular momentum show the largest discrepancies among the diagnostics for Benchmark 3. The solutions by (H), (S), (C), (NLG), (V) seem to converge to the same value within $7 \times 10^{-2}$ per cent. While a very high agreement was achieved for the kinetic energy solutions, simulations by (MJ), (TSH) and (A) seem to converge to a lower value of $L_z$ but still within 0.4 per cent. The last two diagnostic values involve the evaluation of the velocity field at the centre of the spherical domain. While all simulations lie within 0.15 per cent for the velocity along the y-axis, the velocity along the x-axis shows a larger discrepancy with values within 0.3 per cent.

Considering the very fast convergence it showed for all diagnostics, the standard value will be taken as the final solution by (H). As for the other two benchmarks, the error bars are chosen such that at least two additional solutions lie within the given bounds. These standard values and error bars for Benchmark 3 are given in Table 7.

## 5 DISCUSSION

The combination of the results for all three test cases paints a uniform and unambiguous picture of a successful benchmarking exercise. With the wide range of classes of problems covered by

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Summary of the solutions of the participants to Benchmark 3. The vertical velocity component $u_z$ is not shown as it was consistently shown to be zero. The black horizontal line is the standard value given in Table 7. The error corridor for the standard values is represented by a greyed out area. Except for (c) which has larger errors, the error corridor is mostly hidden by the black line. The labels used for the different codes are defined in Section 3.</figcaption></figure>
<table><caption>**Table 5.** Spectral method contributions to Benchmark 3. The labels used for the different codes are defined in Section 3. The values are shown with the number of significant digits provided by the authors. As all these codes are based on a spherical harmonics expansion for the angular component, the resolution is given as the radial resolution *R*, the highest harmonic degree *L* and the highest harmonic order *M*.</caption><thead><tr><th>Code</th><th>E<sub>k</sub></th><th>E<sub>k</sub>, m = 0</th><th>E<sub>k</sub>, m = 1</th><th>E<sub>k</sub>, m = 2</th><th>L<sub>z</sub></th><th>U<sub>x</sub></th><th>U<sub>y</sub></th><th>N</th><th>L</th><th>M</th></tr></thead><tbody><tr><td>(MJ)</td><td>6.18485e-2</td><td>4.35430e-4</td><td>6.12948e-2</td><td>1.17244e-4</td><td>2.78499e-2</td><td>-1.06714e-2</td><td>4.00811e-2</td><td>12</td><td>23</td><td>23</td></tr><tr><td>(MJ)</td><td>6.18338e-2</td><td>4.31265e-4</td><td>6.12841e-2</td><td>1.17319e-4</td><td>2.76877e-2</td><td>-8.91803e-3</td><td>3.86386e-2</td><td>16</td><td>23</td><td>23</td></tr><tr><td>(MJ)</td><td>6.18338e-2</td><td>4.31265e-4</td><td>6.12841e-2</td><td>1.17319e-4</td><td>2.76877e-2</td><td>-8.91803e-3</td><td>3.86386e-2</td><td>16</td><td>31</td><td>31</td></tr><tr><td>(MJ)</td><td>6.18291e-2</td><td>4.32506e-4</td><td>6.12783e-2</td><td>1.17195e-4</td><td>2.77151e-2</td><td>-8.21402e-3</td><td>3.83418e-2</td><td>23</td><td>31</td><td>31</td></tr><tr><td>(MJ)</td><td>6.18291e-2</td><td>4.32506e-4</td><td>6.12783e-2</td><td>1.17195e-4</td><td>2.77151e-2</td><td>-8.21402e-3</td><td>3.83418e-2</td><td>23</td><td>47</td><td>47</td></tr><tr><td>(MJ)</td><td>6.18286e-2</td><td>4.32813e-4</td><td>6.12776e-2</td><td>1.17118e-4</td><td>2.77204e-2</td><td>-8.27886e-3</td><td>3.83230e-2</td><td>31</td><td>47</td><td>47</td></tr><tr><td>(MJ)</td><td>6.18286e-2</td><td>4.32813e-4</td><td>6.12776e-2</td><td>1.17118e-4</td><td>2.77204e-2</td><td>-8.27885e-3</td><td>3.83230e-2</td><td>31</td><td>63</td><td>63</td></tr><tr><td>(H)</td><td>6.1832e-2</td><td>4.3521e-4</td><td>6.1278e-2</td><td>1.1762e-4</td><td>2.7783e-2</td><td>-1.0066e-2</td><td>3.7817e-2</td><td>12</td><td>12</td><td>4</td></tr><tr><td>(H)</td><td>6.1831e-2</td><td>4.3515e-4</td><td>6.1277e-2</td><td>1.1754e-4</td><td>2.7796e-2</td><td>-8.3410e-3</td><td>3.8318e-2</td><td>15</td><td>15</td><td>5</td></tr><tr><td>(H)</td><td>6.1831e-2</td><td>4.3514e-4</td><td>6.1277e-2</td><td>1.1754e-4</td><td>2.7796e-2</td><td>-8.2654e-3</td><td>3.8307e-2</td><td>18</td><td>18</td><td>6</td></tr><tr><td>(H)</td><td>6.1831e-2</td><td>4.3514e-4</td><td>6.1277e-2</td><td>1.1754e-4</td><td>2.7796e-2</td><td>-8.2644e-3</td><td>3.8307e-2</td><td>21</td><td>21</td><td>7</td></tr><tr><td>(H)</td><td>6.1831e-2</td><td>4.3514e-4</td><td>6.1277e-2</td><td>1.1754e-4</td><td>2.7796e-2</td><td>-8.2644e-3</td><td>3.8307e-2</td><td>24</td><td>24</td><td>8</td></tr><tr><td>(S)</td><td>6.15237e-2</td><td>4.25401e-4</td><td>6.09814e-2</td><td>1.15913e-4</td><td>2.748153e-2</td><td>-8.121293e-3</td><td>3.778225e-2</td><td>150</td><td>20</td><td>7</td></tr><tr><td>(S)</td><td>6.17682e-2</td><td>4.33033e-4</td><td>6.12169e-2</td><td>1.17198e-4</td><td>2.772647e-2</td><td>-8.230440e-3</td><td>3.820168e-2</td><td>300</td><td>20</td><td>7</td></tr><tr><td>(S)</td><td>6.17537e-2</td><td>4.32675e-4</td><td>6.12028e-2</td><td>1.17133e-4</td><td>2.771565e-2</td><td>-8.226658e-3</td><td>3.817840e-2</td><td>300</td><td>31</td><td>10</td></tr><tr><td>(S)</td><td>6.18047e-2</td><td>4.34295e-4</td><td>6.1252e-2</td><td>1.17400e-4</td><td>2.776804e-2</td><td>-8.251015e-3</td><td>3.826415e-2</td><td>500</td><td>31</td><td>10</td></tr><tr><td>(S)</td><td>6.18239e-2</td><td>4.34907e-4</td><td>6.12705e-2</td><td>1.17501e-4</td><td>2.778725e-2</td><td>-8.262105e-3</td><td>3.829847e-2</td><td>1000</td><td>31</td><td>10</td></tr><tr><td>(TSH)</td><td>5.98943e-2</td><td>3.63627e-4</td><td>5.94187e-2</td><td>1.10888e-4</td><td>2.53084e-2</td><td>-3.29031e-2</td><td>1.22471e-1</td><td>12</td><td>10</td><td>10</td></tr><tr><td>(TSH)</td><td>6.07755e-2</td><td>4.20777e-4</td><td>6.02401e-2</td><td>1.13600e-4</td><td>2.73295e-2</td><td>-8.27489e-3</td><td>3.81288e-2</td><td>24</td><td>21</td><td>21</td></tr><tr><td>(TSH)</td><td>6.23624e-2</td><td>4.42463e-4</td><td>6.17993e-2</td><td>1.19543e-4</td><td>2.80304e-2</td><td>-8.28375e-3</td><td>3.84550e-2</td><td>32</td><td>42</td><td>42</td></tr><tr><td>(TSH)</td><td>6.15750e-2</td><td>4.31638e-4</td><td>6.10257e-2</td><td>1.16576e-4</td><td>2.76827e-2</td><td>-8.25493e-3</td><td>3.82359e-2</td><td>48</td><td>85</td><td>85</td></tr><tr><td>(A)</td><td>6.1771e-2</td><td>4.3704e-4</td><td>6.1220e-2</td><td>1.1753e-4</td><td>2.7702e-2</td><td>-8.359e-3</td><td>3.8316e-2</td><td>200</td><td>63</td><td>63</td></tr></tbody></table>

<table><caption>**Table 6.** Local method contributions to Benchmark 3. The labels used for the different codes are defined in Section 3. The values are shown with the number of significant digits provided by the authors. The kinetic energy $E_k$ in the $m = 0, m = 1$ and $m = 2$ modes has to be computed in a post-processing step which is likely to introduce additional errors in the codes (C) and (V). For this reason, these values were not mandatory for Benchmark 3. The resolution *R* is given as the third root of the total number of gridpoints $N_{\text{grid}}$.</caption><thead><tr><th>Code</th><th>E<sub>k</sub></th><th>E<sub>k</sub>, m = 0</th><th>E<sub>k</sub>, m = 1</th><th>E<sub>k</sub>, m = 2</th><th>L<sub>z</sub></th><th>U<sub>x</sub></th><th>U<sub>y</sub></th><th>N<sup>1/3</sup><sub>grid</sub></th></tr></thead><tbody><tr><td>(C)</td><td>6.1814e-2</td><td>N/A</td><td>N/A</td><td>N/A</td><td>2.7553e-2</td><td>-8.8469e-3</td><td>4.0492e-2</td><td>26.4</td></tr><tr><td>(C)</td><td>6.1839e-2</td><td>N/A</td><td>N/A</td><td>N/A</td><td>2.7787e-2</td><td>-8.8469e-3</td><td>4.0492e-2</td><td>32.8</td></tr><tr><td>(C)</td><td>6.1831e-2</td><td>N/A</td><td>N/A</td><td>N/A</td><td>2.7808e-2</td><td>-8.2943e-3</td><td>3.8329e-2</td><td>45.0</td></tr><tr><td>(C)</td><td>6.1829e-2</td><td>N/A</td><td>N/A</td><td>N/A</td><td>2.7797e-2</td><td>-8.2943e-3</td><td>3.8329e-2</td><td>54.2</td></tr><tr><td>(C)</td><td>6.1830e-2</td><td>N/A</td><td>N/A</td><td>N/A</td><td>2.7797e-2</td><td>-8.2637e-3</td><td>3.8305e-2</td><td>74.5</td></tr><tr><td>(NLG)</td><td>6.1847e-2</td><td>4.4376e-4</td><td>6.1285e-2</td><td>1.1733e-4</td><td>2.8158e-2</td><td>-8.3649e-3</td><td>3.8332e-2</td><td>56.4</td></tr><tr><td>(NLG)</td><td>6.1831e-2</td><td>4.3534e-4</td><td>6.1277e-2</td><td>1.1754e-4</td><td>2.7803e-2</td><td>-8.2946e-3</td><td>3.8308e-2</td><td>82.0</td></tr><tr><td>(NLG)</td><td>6.1831e-2</td><td>4.3515e-4</td><td>6.1277e-2</td><td>1.1754e-4</td><td>2.7796e-2</td><td>-8.2686e-3</td><td>3.8307e-2</td><td>124.</td></tr><tr><td>(V)</td><td>6.19485e-2</td><td>4.2922e-4</td><td>5.6011e-2</td><td>1.1883e-4</td><td>2.77393e-2</td><td>-8.33047e-3</td><td>3.80877e-2</td><td>139.4</td></tr><tr><td>(V)</td><td>6.19288e-2</td><td>4.3333e-4</td><td>5.7559e-2</td><td>1.1808e-4</td><td>2.77620e-2</td><td>-8.31833e-3</td><td>3.82237e-2</td><td>206.3</td></tr><tr><td>(V)</td><td>6.18951e-2</td><td>4.3379e-4</td><td>5.7707e-2</td><td>1.1707e-4</td><td>2.77724e-2</td><td>-8.28240e-3</td><td>3.82734e-2</td><td>280.8</td></tr></tbody></table>

these three test cases, ranging from purely hydrodynamic problems with thermal or boundary forcing to non-linear dynamo simulations, these results do support the confidence that is put into numerical simulations in a full sphere geometry. The different codes used to compute numerical solutions, while based on wide range of numerical methods, all agreed very well with each other. As was observed in similar benchmarking exercises (e.g. Christensen et al. 2001; Jackson et al. 2013) in a spherical shell geometry, the fully spectral simulations showed the fastest convergence to the final solutions followed by the mixed spherical harmonics and finite difference codes. The simulations using local methods exhibited a very good agreement but required a much higher resolution to converge. However, one should keep in mind that the simple spherical geometry and solutions with a simple structure do favour spectral methods.

With at least five different implementations taking part in each benchmark case, the provided standard values and error bounds can be trusted to be accurate. Benchmark 1 showed the strongest convergence among all the solutions proposed. Maybe somewhat suprisingly, Benchmark 2 showed a somewhat larger discrepancy.

Benchmark 3 showed also quite good convergence from all codes, except for the value of the angular momentum where a larger scatter in the solutions was observed also among the spectral codes that agreed well for Benchmarks 1 and 2. Interestingly, the two codes by (MJ) and (H) did exhibit a remarkably similar behaviour and produced nearly the same results for all values in Benchmarks 1 and 2. While both use a spherical harmonic expansion, the radial discretization is quite different. (H) uses a parity constrained Chebyshev expansion while (MJ) uses a basis set that satisfies the regularity conditions at the origin exactly.

Two physical issues have emerged as part of these calculations. The use of stress-free boundary conditions in Benchmarks 1 and 2 imply that angular momentum must be conserved. As was also discovered in Jones et al. (2011), it was not the case in all the codes. Several groups simply monitored the evolution of the angular momentum and reported no problem with the provided resolutions. On the other hand, some of the codes needed to correct every few time steps to avoid building up unphysical angular momentum. The relatively long integration time required to reach Benchmark 2 did
<table><caption>**Table 7.** Table of the standard values obtained for Benchmark 3. The viscosity $\nu$, the rotation rate $\Omega$ and boundary velocity $u_0$ are used to parametrize Benchmark 3. The governing equation for the velocity $\mathbf{u}$ is given in Section 2. The kinetic energy $E_k$ is defined in eq. (53). $L_z$ is the $\hat{z}$ component of the angular momentum. $u_x$ and $u_y$ are the $\hat{x}$ and $\hat{y}$ components of the velocity through the centre of the bubble.</caption><caption>Benchmark 3: Boundary forced rotating bubble</caption><thead><tr><th colspan="3">Parameters</th></tr><tr><th>ν</th><th>Ω</th><th>u<sub>0</sub></th></tr></thead><tbody><tr><td>10<sup>-2</sup></td><td>10</td><td>√(3/2π)</td></tr><tr><td colspan="3">Boundary conditions</td></tr><tr><td colspan="3">Tangential flow:</td></tr><tr><td></td><td>u<sub>bc</sub> = u<sub>0</sub>∇y<sup>1</sup><sub>1</sub>(θ, φ)</td><td></td></tr><tr><td></td><td>Requested values</td><td></td></tr><tr><td>E<sub>k</sub></td><td>6.1831 × 10<sup>-2</sup></td><td>±1 × 10<sup>-6</sup></td></tr><tr><td>L<sub>z</sub></td><td>2.7796 × 10<sup>-2</sup></td><td>±1 × 10<sup>-6</sup></td></tr><tr><td colspan="3">Velocity through centre</td></tr><tr><td>u<sub>x</sub></td><td>-8.2644 × 10<sup>-3</sup></td><td>±2.3 × 10<sup>-6</sup></td></tr><tr><td>u<sub>y</sub></td><td>3.8307 × 10<sup>-2</sup></td><td>±2 × 10<sup>-6</sup></td></tr></tbody></table>

exacerbate the problem as even small errors do accumulate to a sizeable value over a large number of time steps. (H) did follow a different approach and imposed a modified boundary condition to explicitly impose conservation.

The full sphere dynamo problem is of great geophysical importance, as it accurately represents the Early Earth prior to the formation of the inner core (see Jacobs 1953). The results from Benchmarks 1 and 2 show that even a small inner core may result in solutions that strongly differ from the full sphere solutions. Indeed the use of a small inner core systematically produced less accurate solutions. Problems in a full sphere geometry, like the simulation of Early Earth’s dynamo, should be addressed with specialized codes. It is expected that this issue will become even more important in more complex flows.

## ACKNOWLEDGEMENTS

This work was partly funded by ERC grant 247303 ‘MFECE’ to AJ. AJ and PM acknowledge SNF grant 200021-113466, and are grateful for the provision of computational resources by the Swiss National Supercomputing Centre (CSCS) under project ID s225. PM acknowledges the support of the NSF CSEDI Programme through award EAR-1067944. For this work, DC was supported by the ETH Zurich Post-doctoral fellowship programme as well as by the Marie Curie Actions for People COFUND Programme. JLG acknowledges support from the University Paris-Sud, the National Science Foundation (grants NSF DMS-1015984 and DMS-1217262), the Air Force Office of Scientific Research, (grant FA99550-12-0358) and the King Abdullah University of Science and Technology (Award No. KUS-C1-016-04). The computations using SFEMANS were carried out on IBM SP6 of GENCI-IDRIS (project 0254).

## REFERENCES

Aubert, J., Aurnou, J. & Wicht, J., 2008. The magnetic structure of convection-driven numerical dynamos, *Geophys. J. Int.*, **172**(3), 945–956.

Balay, S., Gropp, W.D., McInnes, L.C. & Smith, B.F., 1997. Efficient management of parallelism in object-oriented numerical software libraries, *Modern Software Tools in Scientific Computing*, pp. 163–202, eds Arge, E., Bruaset, A.M. & Langtangen, H. P., Birkhäuser, Boston.

Balay, S. et al., 2012a. PETSc Users Manual Revision 3.3.

Balay, S. et al., 2012b. PETSc Web page, available at: http://www.mcs.anl.gov/petsc, last accessed 12 January 2014.

Boyd, J.P., 2001. Chebyshev and Fourier Spectral Methods, Courier Dover.

Busse, F. & Simitev, R., 2008. Toroidal flux oscillation as possible cause of geomagnetic excursions and reversals, *Phys. Earth planet. Inter.*, **168**(3), 237–243.

Busse, F.H. & Simitev, R.D., 2006. Parameter dependences of convection-driven dynamos in rotating spherical fluid shells, *Geophys. astrophys. Fluid Dyn.*, **100**(4–5), 341–361.

Christensen, U. et al., 2009. Erratum to a numerical dynamo benchmark [Phys. Earth planet. Inter. **128**(1–4) (2001) 25–34], *Phys. Earth planet. Inter.*, **172**(3–4), 356.

Christensen, U.R. et al., 2001. A numerical dynamo benchmark, *Phys. Earth planet. Inter.*, **128**, 25–34.

Dormy, E., Cardin, P. & Jault, D., 1998. MHD flow in a slightly differentially rotating spherical shell, with conducting inner core, in a dipolar magnetic field, *Earth planet. Sci. Lett.*, **160**(1–2), 15–30.

Glatzmaier, G.A. & Roberts, P.H., 1995. A three-dimensional convective dynamo solution with rotating and finitely conducting inner core and mantle, *Phys. Earth planet. Inter.*, **91**(1–3), 63–75.

Guermond, J.-L., Laguerre, R., Léorat, J. & Nore, C., 2007. An interior penalty Galerkin method for the MHD equations in heterogeneous domains, *J. Comput. Phys.*, **221**(1), 349–369.

Guermond, J.-L., Laguerre, R., Léorat, J. & Nore, C., 2009. Nonlinear magnetohydrodynamics in axisymmetric heterogeneous domains using a Fourier/finite element technique and an interior penalty method, *J. Comput. Phys.*, **228**(8), 2739–2757.

Guermond, J.-L., Léorat, J., Luddens, F., Nore, C. & Ribeiro, A., 2011. Effects of discontinuous magnetic permeability on magnetodynamic problems, *J. Comput. Phys.*, **230**(16), 6299–6319.

Hauke, G. & Hughes, T., 1994. A unified approach to compressible and incompressible flows, *Comput. Methods Appl. Mech. Eng.*, **113**(3–4), 389–395.

Hindmarsh, A.C., Brown, P.N., Grant, K.E., Lee, S.L., Serban, R., Shumaker, D.E. & Woodward, C.S., 2005. Sundials: suite of nonlinear and differential/algebraic equation solvers, *ACM Trans. Math. Softw.*, **31**(3), 363–396.

Hollerbach, R., 2000. A spectral solution of the magneto-convection equations in spherical geometry, *Int. J. Numer. Methods Fluids*, **32**(7), 773–797.

Jackson, A. et al., 2013. A spherical shell numerical dynamo benchmark with pseudo-vacuum magnetic boundary conditions, *Geophys. J. Int.*, **196**(2), 712–723.

Jacobs, J.A., 1953. The Earth’s inner core, *Nature*, **172**(4372), 297–298.

Jones, C.A., Boronski, P., Brun, A.S., Glatzmaier, G.A., Gastine, T., Miesch, M.S. & Wicht, J., 2011. Anelastic convection-driven dynamo benchmarks, *Icarus*, **216**(1), 120–135.

Jouve, L. et al., 2008. A solar mean field dynamo benchmark, *A&A*, **483**(3), 949–960.

Kageyama, A. & Sato, T. Complexity Simulation Group, 1995. Computer simulation of a magnetohydrodynamic dynamo. II, *Phys. Plasmas*, **2**(5), 1421–1431.

Kaiser, J. & Kuo, F., 1966. *Systems Analysis by Digital Computer*, Wiley.

Karypis, G. & Kumar, V., 2009. METIS: Unstructured Graph Partitioning and Sparse Matrix Ordering System, Version 4.0, University of Minnesota.

Kim, J. & Moin, P., 1985. Application of a fractional-step method to incompressible navier-stokes equations, *J. Comput. Phys.*, **59**(2), 308–323.

Livermore, P.W., 2012. The spherical harmonic spectrum of a function with algebraic singularities, *J. Fourier Anal. Appl.*, **18**, 1146–1166.

Livermore, P.W. & Hollerbach, R., 2012. Successive elimination of shear layers by a hierarchy of constraints in inviscid spherical-shell flows, *J. Math. Phys.*, **53**(7), 073104, doi:10.1063/1.4736990.

Livermore, P.W., Jones, C.A. & Worland, S.J., 2007. Spectral radial basis functions for full sphere computations, *J. Comput. Phys.*, **227**(2), 1209–1224.

Marti, P., 2012. Convection and boundary driven flows in a sphere, *PhD thesis*, ETH Zurich.
Matsushima, T. & Marcus, P., 1995. A spectral method for polar coordinates, *J. Comput. Phys.*, **120**(2), 365–374.

Monteux, J., Schaeffer, N., Amit, H. & Cardin, P., 2012. Can a sinking metallic diapir generate a dynamo? *J. geophys. Res. Planets*, **117**, E10005, doi:10.1029/2012JE004075.

Oliphant, T.E., 2007. Python for scientific computing, *Compu. Sci. Eng.*, **9**(3), 10–20.

Sasaki, Y., Takehiro, S. & Hayashi, Y.-Y. SPMODEL Development Group, 2012. Project of MHD Dynamo in rotating spheres and spherical shells. Available at: http://www.gfd-dennou.org/library/dynamo/, last accessed 9 January 2014.

Schaeffer, N., 2013. Efficient spherical harmonic transforms aimed at pseudo-spectral numerical simulations, *Geochem. Geophys. Geosyst.*, **14**(3), 751–758.

Simitev, R. & Busse, F., 2005. Prandtl-number dependence of convection-driven dynamos in rotating spherical fluid shells, *J. Fluid Mech.*, **532**, 365–388.

Simitev, R. & Busse, F., 2012. How far can minimal models explain the solar cycle? *Astrophys. J.*, **749**(1), 9, doi:10.1088/0004-637X/749/1/9.

Simitev, R.D. & Busse, F.H., 2009. Bistability and hysteresis of dipolar dynamos generated by turbulent convection in rotating spherical shells, *Europhys. Lett.*, **85**(1), 19001, doi:10.1209/0295-5075/85/19001.

Tilgner, A., 1999. Spectral methods for the simulation of incompressible flows in spherical shells, *Int. J. Numer. Methods Fluids*, **30**(6), 713–724.

Tilgner, A. & Busse, F., 1997. Finite-amplitude convection in rotating spherical fluid shells, *J. Fluid Mech.*, **332**(1), 359–376.

Vantieghem, S., 2011. Numerical simulations of quasi-static magnetohydrodynamics using an unstructured finite volume solver: development and applications, *PhD thesis*, Université Libre de Bruxelles.

Worland, S., 2004. Magnetocvection in rapidly rotating spheres, *PhD thesis*, University of Exeter.

Zhan, X., Schubert, G. & Zhang, K., 2012. Anelastic convection-driven dynamo benchmarks: a finite element model, *Icarus*, **218**(1), 345–347.

## APPENDIX A: SPHERICAL HARMONICS

The Schmidt quasi-normalized spherical harmonic basis was used to provide a simpler expression for the initial conditions. The spherical harmonic $Y_l^m$ of degree $l$ and order $m$ is given by

$$Y_l^m(\theta, \varphi) = P_l^m(\cos \theta)e^{im\varphi}, \qquad (A1)$$

where the $P_l^m(\cos \theta)$ are the Schmidt quasi-normalized associated Legendre functions. The above definition of the spherical harmonic can also be written as function of the normalized associated Legendre functions $\hat{P}_l^m(\cos \theta)$ leading to the expression

$$Y_l^m(\theta, \varphi) = \sqrt{\frac{(l-m)!}{(l+m)!}} \hat{P}_l^m(\cos \theta) e^{im\varphi}. \qquad (A2)$$

The orthogonality relation for the $Y_l^m$ defined above is given by

$$\int_{0}^{\pi} \int_{0}^{2\pi} Y_l^m Y_{l'}^{m'*} d\Omega = \frac{4\pi(2 - \delta_{m0})}{2l + 1} \delta_{ll'} \delta_{mm'}, \qquad (A3)$$

where the * denotes the complex conjugate.