On the sensitivity of 3-D thermal convection codes to numerical discretization: a model intercomparison

. Fully 3-D numerical simulations of thermal convection in a spherical shell have become a standard for studying the dynamics of pattern formation and its stability under perturbations to various parameter values. The question arises as to how the discretization of the governing equations affects the outcome and thus any physical interpretation. This work demonstrates the impact of numerical discretization on the observed patterns, the value at which symmetry is broken, and how stability and stationary behavior is dependent upon it. Motivated by numerical simulations of convection in the Earth’s mantle, we consider iso-viscous Rayleigh–Bénard convection at inﬁnite Prandtl number, where the aspect ratio between the inner and outer shell is 0.55. We show that the subtleties involved in developing mantle convection models are considerably more delicate than has been previously appreciated, due to the rich dynamical behavior of the system. Two codes with different numerical discretization schemes – an established, community-developed, and benchmarked ﬁnite-element code (CitcomS) and a novel spectral method that combines Chebyshev with radial basis functions (RBFs) – are compared.


Abstract.
Fully 3-D numerical simulations of thermal convection in a spherical shell have become a standard for studying the dynamics of pattern formation and its stability under perturbations to various parameter values. The question arises as to how the discretization of the governing equations affects the outcome and thus any physical interpretation. This work demonstrates the impact of numerical discretization on the observed patterns, the value at which symmetry is broken, and how stability and stationary behavior is dependent upon it. Motivated by numerical simulations of convection in the Earth's mantle, we consider isoviscous Rayleigh-Bénard convection at infinite Prandtl number, where the aspect ratio between the inner and outer shell is 0.55. We show that the subtleties involved in developing mantle convection models are considerably more delicate than has been previously appreciated, due to the rich dynamical behavior of the system. Two codes with different numerical discretization schemes -an established, communitydeveloped, and benchmarked finite-element code (CitcomS) and a novel spectral method that combines Chebyshev polynomials with radial basis functions (RBFs) -are compared. A full numerical study is investigated for the following three cases. The first case is based on the cubic (or octahedral) initial condition (spherical harmonics of degree = 4). How this pattern varies to perturbations in the initial condition and Rayleigh number is studied. The second case investigates the stability of the dodecahedral (or icosahedral) initial condition (spherical harmonics of degree = 6). Although both methods first converge to the same pattern, this structure is ultimately unstable and systematically degenerates to cubic or tetrahedral symmetries, depending on the code used. Lastly, a new steady-state pattern is presented as a combination of third-and fourth-order spherical harmonics leading to a fivecell or hexahedral pattern and stable up to 70 times the critical Rayleigh number. This pattern can provide the basis for a new accuracy benchmark for 3-D spherical mantle convection codes.

Introduction
For 3-D Rayleigh-Bénard convection in a spherical shell at infinite Prandtl number, analytical studies by Busse (1975) and Riahi (1988, 1982), using weakly nonlinear perturbation theory, predicted a set of solutions that exhibited steady-state polyhedral pattern formations that would also persist into stronger nonlinear regimes. Later, these solutions were numerically verified by Bercovici et al. (1989Bercovici et al. ( , 1991, Ratcliff and Schubert (1996), and Machetel et al. (1986) for up to 100 times the critical Rayleigh number (Ra = 712) (such as the cubic symmetry test case -which forms the corner structure of an octahedron). Some studies, such as Bercovici et al. (1991), have questioned the properties of these steady-state solutions by considering the influence of the nondominant spherical harmonic modes on modifying boundary layer thickness as the Rayleigh number increased. However, the stability of these polyhedral patterns to perturbations in the initial conditions, i.e the dominant spherical harmonic modes that actually define them remains unclear. Furthermore, the dynamical behavior of steady-state solutions with higher orders of polyhedral symmetry predicted by Busse (1975)   examined. From a computational standpoint, each numerical scheme will handle unstable steady states, nonuniqueness in the solution, and bifurcations differently, depending on how the continuous eigenvalue spectrum has been discretely represented when linearized around the steady state.
In this light, the goal of this paper is to illustrate that the subtleties involved in the development of numerical mantle convection models are considerably more delicate than has been previously appreciated, due to the rich dynamical behavior of the system. For fully nonlinear large-scale systems with millions of unknowns, as considered in this paper, using classical eigenvalue stability analysis to understand the influence of numerical discretization is not an option as (1) the analytical solution and thus the continuous eigenvalue spectrum is not available and (2) calculating the eigenvalues for such systems is computationally not feasible. Although recent advancements have been made in developing iterative schemes to detect Hopf bifurcations in large-scale systems (Meerbergen and Spence, 2010;Elman et al., 2012), the following study exhibits a much richer pattern of dynamical instability and transitional behavior, leading to a variety of end states. Therefore, we will perform an intensive computational investigation of the stationary behavior and stability of three different types of symmetries with regard to perturbations on the initial condition and as a function of Ra, observing how both transitional and end states are strongly dependent on numerical discretization.
The numerical studies are done using two state-ofthe-art models, CitcomS-3.1.1 (http://www.geodynamics. org/cig/software/citcoms) and a pseudospectral radial basis function-Chebyshev model (Wright et al., 2010) (RBF-PS), with the former funded on a ongoing basis by the USA National Science Foundation. Section 2 provides an overview of the system of partial differential equations (PDEs) to be solved and the computational methods used. Section 3 numerically studies the sensitivity of the steady-state solution to perturbations in the cubic initial condition for both low and higher Ra number. Section 4 explores the stability regimes of a higher order initial spherical symmetry, studying the transition between steady states as a function of Ra. Section 5 introduces a new initial condition mode, leading to the observation of a novel steady-state pattern and future benchmark for assessing model performance.

Governing equations and computational models
The governing equations describe a Boussinesq fluid at infinite Prandtl number in a 3-D spherical shell that is heated from below and cooled from above: where u = (u r , u θ , u λ ) is the velocity field in spherical coordinates (θ: latitude; λ: longitude), p is pressure, T is temperature,r is the unit vector in the radial direction, η is the viscosity, and Ra is the Rayleigh number defined below. The boundary conditions on the fluid velocity at the inner and outer surfaces of the spherical shell are and r ∂ ∂r where R i = 11/9 is the radius of the inner surface of the 3-D spherical shell and R o = 20/9 is the radius of the outer surface as measured from r = 0. The boundary conditions on the temperature are Equations (1)-(3) are nondimensionalized with the length scale chosen as the approximate thickness of the mantle, R = R o −R i = 1; the timescale chosen as the thermal diffusion time of mantle minerals, t = ( R) 2 /κ (noting a nondimensional time t = 1 corresponds to 265 billion years, i.e., 58 times the age of the Earth); and the temperature scale chosen as the difference between the temperature at the inner and outer boundaries, T = 1. The fluid is treated as isoviscous: η = constant. Thus, the dynamics of the fluid are governed entirely by the Ra, which can be interpreted as a ratio of the destabilizing force due to the buoyancy of the heated fluid to the stabilizing force due to the viscosity of the fluid and heat transfer by conduction. It takes the specific form of where ρ 0 is a reference density of the fluid, g is the acceleration due to gravity, and α is the coefficient of thermal expansion. The initial condition for the temperature is specified as with The first term in Eq. (6) represents a purely conductive temperature profile, while the second term T P is a perturbation to this profile, determining the final patterns of polyhedral symmetry. Y m denotes the normalized spherical harmonic of degree and order m (Eq. 8) and the nonaxisymmetric perturbation ε will play an important role in studying transitional pattern formations in the cubic case.
where P m are the (unnormalized) associated Legendre functions and δ m0 is the Kronecker delta. It should be noted that the stability of preferred patterns in purely axisymmetric convective flows has been studied by Zebib et al. (1980Zebib et al. ( , 1983.

CitComS
CitcomS is a second-order finite-element code written in C. Its purpose is to explore mantle convection problems in 3-D spherical geometry (Zhong et al., 2000;Tan et al., 2006). Developed from the software Citcom (Moresi and Solomatov, 1995;Moresi et al., 1996), a code structured for 3-D Cartesian geometry, CitcomS employs an Uzawa algorithm to solve the momentum equation coupled with the incompressibility constraints (Ramage and Wathen, 1994). The energy equation is solved with a streamline upwind Petrov-Galerkin method (Brooks, 1981). We used version 3.1.1, available from the Computational Infrastructure for Geodynamics (http://www.geodynamics.org/cig/software/citcoms). The global mesh is obtained by first dividing the spherical shell into 12 caps of approximatively equal size. Then each cap is divided into N × N elements in the angular directions and M elements in the vertical direction, forming a layered brick-like structure. For each 3-D element, eight velocity nodes with trilinear interpolation functions and one constant pressure node are used. Per cap, we will be using 48 elements in each dimension, resulting in 12 × 48 × 48 × 48 total elements.

RBF-PS
Here, an overview of the spectral RBF-PS model is given; for a detailed description of the numerical method, see Wright et al. (2010). To spatially discretize the 3-D spherical shell, a "2(θ, λ) + 1(r)" layered approach is used. In the radial direction, M + 2 Chebyshev nodes (corresponding to M interior points and two boundary points) and N "scattered" nodes (for example, see Womersley andSloan, 2003/2007) are placed on each of the resulting M spherical surfaces. This gives a tensor product structure between the radial and lateral directions, which allows the spatial operators to be computed While all radial derivatives are discretized using Chebyshev polynomials, differential operators in the latitudinal direction (θ) and longitudinal direction (λ) are approximated discretely on each spherical surface using radial basis functions (RBFs). Within a given limit, RBFs reproduce spherical harmonics Piret, 2007, 2008). However, they generally give higher accuracy than spherical harmonics for nonlinear systems of PDEs (Wright et al., 2010;Wright, 2007, 2009;Flyer and Fornberg, 2011;Flyer et al., 2012) (for examples of how to implement RBFs on spherical surfaces, see Wright, 2007, 2009). For all cases in the paper, N = 4096 nodes are used on each sphere, with M = 43 Chebyshev nodes used in the radial direction. The time discretization of the energy equation uses a semi-implicit method. All terms that involve radial derivatives are time-stepped with a Crank-Nicolson method, while terms involving latitudinal and longitudinal derivatives are time-stepped with a third-order Adams-Bashforth method.

Stability of cubic steady state with regard to perturbations in the initial condition at varying Ra
The cubic initial condition temperature profile used in many 3-D spherical convection studies -such as Ratcliff and Schubert (1996), Kameyama et al. (2008), Zhong et al. (2000), Yoshida and Kageyama (2004), Stemmer et al. (2006), Choblet et al. (2007), Zhong et al. (2008), and Kameyama et al. (2008) -is specified by letting T P in Eq. (7) be equal to where δ = 0. A perturbation parameter δ has been introduced to allow us to slowly perturb the amplitude of the nonaxisymmetric mode. The θ −λ temperature dependence of Eq. (9) on a spherical shell surface can be seen in Fig. 1a for δ = 0. As δ increases the initial condition slowly tends to a pure Y 0 4 initial condition, with the amplitude of the four plumes along the equatorial region decreasing and progressively merging together as seen in Fig. 1b for δ = 0.30. It should be noted here that δ = 0 does not correspond to perfect cubic symmetry but has, however, become the standard in modern geophysical and astrophysical simulations such as those cited above. Indeed, the maximum amplitude of the plumes in Fig. 1a varies slightly between the poles and the equator. Perfect cubic symmetry, as predicted by Busse (1975), numerically discovered by Young (1974), and with early simulations by Machetel et al. (1986) and Bercovici et al. (1989), is obtained with ε = 5 7 instead of 5 7 . In the next two subsections, we examine how transitions from the cubic steady state to axisymmetric patterns of lower order occur as a function of perturbing the nonaxisymmetric mode of the initial condition, and more interestingly how these transitions differ depending on the numerical discretization of the governing equations.

Sensitivity to amplitude perturbations in the initial condition at low Ra
At low Rayleigh number, 5000 ≤ Ra ≤ 10 000, the cubic steady-state pattern is stable for both models up to δ = 0.30. In other words, the cubic pattern is maintained if the ratio of the spectral coefficient of Y 4 4 to Y 0 4 is at least 1/2. This can be seen in Fig. 2, where the isosurfaces of residual temperature (see caption for further details, as this is how 3-D convection will be illustrated in the paper) are plotted as a function of δ.
Incrementing δ by 0.01, the RBF-PS model displays a clear transition between the cubic steady state and an order = 4 axisymmetric pattern. In contrast, CitcomS converges to a transitional steady-state pattern for 0.31 ≤ δ ≤ 0.32, in which the four plumes along the equator grow and merge together two by two, but the process is not completed. This is never observed with the RBF-PS discretization (see Fig. 2). At higher values of δ, CitcomS and the RBF-PS method converge to the same pattern. Thus, at the parameter value of destabilization (δ = 0.30), the numerical discretization plays an important role in what pattern emerges. Also, the transition point at which the Y 0 4 spherical harmonic mode completely dominates and the Y 4 4 part of the initial condition no longer influences the final pattern of convection differs between the two models. Figure 3 shows the evolution of the volume-averaged temperature (< T >) for both models at Ra = 7000. As just discussed, the figure illustrates that CitcomS converges to three different steady states, depending on the value of δ. In contrast, for δ > 0.30, the figure shows that the RBF-PS solution is attracted to the = 4 axisymmetric mode. In either case, the solution, once destabilized, transitions to patterns characterized by a higher < T >.

Sensitivity to amplitude perturbations in the initial condition at high Ra
As would be expected, at higher Ra, the cubic steady state is much more sensitive to small perturbations in the initial condition. For Ra = 70 000, the Y 4 4 mode of the initial condition was very slowly perturbed in increments of δ = 5 × 10 −3 , as shown in Fig. 4. The cubic steady state is destabilized at δ ≥ 0.065 for CitcomS and δ ≥ 0.070 for the RBF-PS method with different transitional patterns.
With CitcomS, the destabilization shows a transitional pattern between a cubic steady state to an unsteady axisymmetric pattern at δ = 0.065 and δ = 0.007, characterized by two diametrically opposed upwelling plumes in the equatorial region with a great circle of downwelling that encompasses the polar regions. It develops via a two-by-two merging of upwelling plumes on the equator; initial upwelling plumes at the poles are destabilized and migrate to the equatorial region. The end state for perturbations of δ ≥ 0.75 is also an unsteady axisymmetric pattern. However, the pattern of convection has been completely rearranged, with upwelling now occurring in the polar regions and downwelling in the equatorial region, yielding a strong dominance of an oscillating = 2 mode. The quasi-uniform oscillation of this end state can be seen in the time traces of the outer N u and volumeaveraged root-mean-square (rms) velocity in Fig. 5, where the region for t ≥ 0.2 has been enlarged for better viewing.
With the RBF-PS model, the cubic steady state also eventually evolves into an unsteady axisymmetric pattern for δ ≥ 0.085, similar to that of the CitcomS as shown in Fig. 4. However, the transition between these two states is very different than what was observed with the CitcomS model. For δ = 0.07, the cubic pattern is only partially destabilized. Two plumes on one side merge and begin to pulsate. Although this structure is unsteady, it stays stable with no other changes in the general pattern of convection observed. At δ > 0.075, the cubic geometry is fully destabilized and the model begins to converge to the unsteady axisymmetric pattern, seen for δ >= 0.09.
For the two methods, the stability of the cubic symmetry pattern as a function of the Rayleigh number and the amount of perturbation δ to the initial condition is summarized in Fig. 6. The amount of perturbation needed to destabilize the steady-state cubic symmetry pattern begins to decrease rapidly after Ra ≈ 20 000. The shaded blue and pink regions depict where transition states are observed for the CitcomS model and RBF-PS model, respectively. Generally, the evolution of the transition is well defined using both methods. Ra > 30 000. In all cases, using RBF-PS, the transition is not characterized by a single pattern, as in CitcomS, but by a progressive transition as a function of the perturbation (δ). Surprisingly, this transitional regime broadens for large Rayleigh numbers (see red shaded area with Ra ≥ 50 000), implying larger perturbations are required to fully diminish the influence of the = 4 modes. These results clearly demonstrate how numerical discretization impacts pattern formation and its interpretation in simulations of 3-D convective flow.

Stability at higher orders of symmetry: a dodecahedral initial condition
In Busse (1975), a steady-state higher-order convection pattern corresponding to dodecahedral symmetry is predicted.
Here, for the first time (to the authors' knowledge), the stability of this pattern for low Ra is studied, with surprising results on how the numerical discretization scheme severely affects the interpretation of steady-state stability ranges. The initial condition is given by Eq. (6) with The θ − λ temperature dependence on a sphere is shown in Fig. 7. It has 12 initial plumes of upwelling, forming the faces of a dodecahedron, where the strongest downwelling (in dark blue) occurs at the vertices of the pentagons.
The evolution of convection with a dodecahedral initial condition at a Ra = 7000 is presented in Fig. 8. Both methods first converge to a steady-state dodecahedral pattern; however, this convection pattern is unstable. The symmetry is broken at different times for RBF-PS and CitcomS models. Plumes begin to merge after t = 0.7 with CitcomS, while for the RBF-PS model, plumes do not merge until t = 2.7. Surprisingly, the final stable stationary state differs between the two numerical discretizations: RBF-PS converges to a tetrahedral pattern, dominated by a = 3 mode, while Cit-comS reaches the cubic pattern studied in the previous section. In order to reduce the possible effect of spatial discretization error, mesh resolution in CitcomS was increased by a factor of 8 to 12 × 96 3 and doubled in the RBF-PS to 51(r) × 6561(θ, λ). The results are displayed in Fig. 8. The same final patterns are observed, with the only difference being that the dodecahedral pattern is maintained for a longer period. These results imply that there are at least two stable branches of solutions that correspond to these patterns; however, which branch manifests itself in simulations is dependent on the numerical discretization. We will see more evidence of this later in the discussion.  Fig. 2 for 5000 ≤ Ra ≤ 10 000 and Fig. 4 for Ra = 70 000. For each model, the bottom curve is the maximum δ value for which it converges to the cubic steady state, while the top curve is the minimum δ value which converges to an axisymmetric pattern (steady or unsteady). The black dotted line is the value of the Rayleigh number that marks the transition between the = 4 steady state and the unsteady state. For CitcomS, the mesh discretization shows a symmetrical effect. The shell is initially divided in 12 caps. Each cap is diametrically opposite to another one (Zhong et al., 2000). Thus during the transition, we can observe that destabilization occurs in symmetrical pairs with respect to the caps. As a result it is reasonable to presume that mesh discretization and the cap divisions influence the distribution of numerical errors and favor even modes. Under these conditions, CitcomS will not reproduce the tetrahedral or the five-cell pattern observed with the RBF method without adding an additional initial perturbation representing these odd modes.
The stability of the dodecahedral steady-state solution for 2000 ≤ Ra ≤ 10 000 can also be seen in the time evolution of the volume-averaged rms velocity and the inner and outer Nusselt numbers as given in Fig. 9. In all cases of the Ra, the dodecahedral convection pattern is initially observed and stationary. This pattern is identical in both methods, whether one considers its geometry, the convergence of rms velocity, average temperature, or Nusselt numbers before the transition (Table 1). However, weakly unstable modes of lower spherical harmonic degree become excited and cause the solution to transition to a second steady state. When this transition occurs in the time evolution is clearly dependent on the model. For instance at t = 2, CitcomS has already reached a steady-state cubic pattern while RBF-PS is still in the weakly unstable steady-state dodecahedral pattern.
As the Rayleigh number increases from 2000, the final stationary pattern observed varies greatly between the two models, also showing how preferred patterns of convection in numerical simulations are dependent on the spatial discretization scheme. Figure 10 illustrates these end states for both numerical methods, starting from the dodecahedral initial condition for 2000 ≤ Ra ≤ 70000 for each of the models. The RBF-PS model shows a clear transition from the dodecahedral pattern to a variety of steady states, depending on the Rayleigh number. For 3000 ≤ Ra ≤ 5750, end-state convection is dominated by the cubic steady-state pattern discussed in the previous section. In CitcomS, this pattern is not seen until Ra = 5000, and persists to Ra = 7000. In fact, the regime within 5000 ≤ Ra ≤ 5750 is the only range where both the CitcomS and RBF-PS transition to the same final steady-state convection pattern. At 5775 ≤ Ra ≤ 6025,  a newly observed five-cell pattern emerges as the end stationary state in the RBF-PS model. It results from a mixed-mode interaction between the = 3 and = 4 modes, as will be discussed in the next section. For 6000≤ Ra ≤ 10 000, the final pattern of convection for RBF-PS is the tetrahedral pattern observed in Fig. 8. In contrast, CitcomS transitions to a stable steady-state axisymmetric = 2 pattern. For Ra > 10 000, the final patterns become unsteady, yet maintain a resemblance to the axisymmetric and tetrahedral patterns seen in CitcomS and RBF-PS, respectively.  Figure 10. Transition of dodecahedral plume formation to different steady states for 2000 ≤ Ra ≤ 10 000 and to an unsteady state for 10 000 < Ra ≤ 70 000 with the RBF-PS method (right column) and CitcomS (left column).

A new convection mode: five cells
At 5750 ≤ Ra ≤ 6050 with RBF-PS method, the weakly unstable dodecahedral pattern relaxed into a steady-state fivecell convection pattern. This structure is characterized by five upwelling plumes: two at the poles, each surrounded by a triangular region of downwelling, and three along the equator, each surrounded by a square region of downwelling. The pattern appeared for a narrow range of Rayleigh numbers, between the cubic pattern at lower Rayleigh number and the tetrahedral pattern at higher Rayleigh number. This observation along with the fact that the convective regions of descending motion are defined by both the vertices of a triangle in the polar regions (the case for the tetrahedral pattern) and those of a square in the equatorial regions (the case for the cubic pattern) leads us consider a mixed-mode interaction between the = 3 and an = 4 modes for an initial condition. Previous studies of mixed-mode patterns bifurcating from spherically symmetric ones have been predicted in Busse and Riahi (1988) and numerically observed in Feudel et al. (2011). However, these studied reported a seven cell pattern resulting from an interaction of a = 4 and = 5 modes. In Chossat and Beltrame (2014), the authors investigated = 3, 4 mode interactions in a context compatible with Rayleigh-Bénard convection without having highlighted the occurrence of a five-cell structure. Here, we focus on the formation of a steady-state five-cell pattern that is stable at large Rayleigh number, Ra = 50 000, approximately 70 times the critical Rayleigh number (Ra c = 712, the onset of convection). Through numerical experimentation, we discovered that a combination of Y 3 3 and Y 0 4 spherical harmonics will yield a five-cell pattern. However, in order to determine the volumeaveraged spectral energies or variances between the two modes that yield the fastest stabilization in a five-cell steadystate pattern, a parameter γ on the Y 0 4 mode is introduced. It will be varied slowly from γ equals 0 to 1, in increments of 10 −2 . The initial condition is now given by Eq. (6) with Figure 11 displays the initial conditions for two different values of γ that will lead to two different steady states.
We begin at a low Rayleigh number, such as Ra = 7000. However, the results hold for even lower Rayleigh numbers, down to Ra = 1000, just above the onset of convection. Figure 12 shows the evolution of the volume-averaged temperature and the final convection patterns (isosurfaces of δT = ±0.15; yellow: ascending motion; blue: descending motion) as γ is varied. As can be seen, depending on the value of γ , the model converges to three distinct steady states. For γ ≤ 0.20, the = 4 mode has no influence and the models converge to a steady state defined by the Y 3 3 spherical harmonic mode. This pattern is similar to that found in Busse and Riahi (1988), except that there is a merging of the ascending motion in the polar regions. The steady-state five-cell pattern shown in Fig. 12d manifests itself in both models for 0.2 ≤ γ < 0.3, with the fastest stabilization to this state for γ = 0.5. As a result, this is what will be used when observing the stability of the five-cell pattern as a function of Rayleigh number. With the RBF-PS model, once the volumeaveraged spectral energies between the two modes are equal (i.e., γ = 1), the flow reverts to an axisymmetric steady state, dominated by the = 4 mode. With the CitcomS model, the ratio of the modes only have to be within 10 % of one another (i.e., γ = 0.9) for this to occur. Lastly, Fig. 13 shows that this convection pattern is not only steady but stable with respect to perturbing the Rayleigh number for values at least up to Ra = 50 000, 70 times the critical Rayleigh number. Both models obtained this result. Also, as the Rayleigh number increases, the boundary layer thickness decreases, as would be expected with increased convection.

Conclusions
In time-dependent fully nonlinear systems, when numerical simulations are performed, a great variety of complex spatiotemporal regimes can be observed depending on parameter values. However, what this paper has illustrated is that what patterns are actually observed and at which parameter values they manifest themselves is definitely impacted by the numerical discretization used. Since computation has become a third arm of physical understanding, along with experimentation and analysis, it is important to highlight this fact so that a discretization scheme is not blindly applied just because it is commonly used, as in the case of spherical harmonics. Here, we have compared an RBF-Chebyshev discretization (RBF-PS) to a finite-element discretization. The latter is a community-based model called CitcomS, specially designed for studying thermal convection in a 3-D spherical shell. For simpler spherical symmetries as the cubic pattern (sometimes referred to as the octahedral pattern), the results at low Rayleigh number were more similar between the models, both destabilizing when the contribution of the nonaxisymmetric = 4 spherical harmonic mode in the initial condition fell below 50 %. However, CitcomS showed a transition to three steady states as this mode was perturbed, while RBF-PS went directly to the = 4 axisymmetric mode. At higher Rayleigh number, the difference in the transitional states manifested between the two models was more drastic.
The effect of the numerical discretization on pattern formation at higher orders of symmetry, such as dodecahedral symmetry where the initial condition is defined by a combination of = 6 spherical harmonic modes, was of even greater interest. Although deemed a stable state by Busse (1975) for Rayleigh numbers near the onset of convection (Ra c = 712), it was shown to be unstable (after a long computational period -equivalent to 25 times the age of the Earth) for a Rayleigh number of just 2.5 times Ra c at extremely high resolutions for both models. However, regardless of the Rayleigh number, the convection evolved completely differently for each model, with the end steady state also being very different. For example, at Ra = 7000, the RBF-PS model evolved into a tetrahedral symmetry, and Cit-comS into a cubic symmetry.
Another outcome of differences in numerical discretization can be the discovery of a stable convection pattern (with regard to perturbations in the Rayleigh number) that does not seem to have been highlighted in the literature. In studying the dodecahedral convection pattern, in a narrow range of the Rayleigh number, the RBF-PS model stabilized to a five-cell steady-state pattern that was never seen in the Cit-comS model regardless of the Rayleigh number. This led the authors to investigate its formation, discovering that it is a strongly stable steady-state pattern of convection up to Ra = 50 000. Both models agreed that it forms by the interaction of the Y 3 3 and Y 0 4 modes. As a general observation, both methods show a good match in the cubic and five-cell steady-state patterns, and even for the stationary dodecahedral pattern before the transition. However, the above in-depth computational study strongly illustrates how numerical discretization can impact both the resulting patterns of convection and the transitional states that occur. This is particularly true when scientists have to rely on such simulations in cases of strongly nonlinear systems with over a million unknowns. In such cases, eigenvalue stability analysis is simply not an option. Furthermore, we hope to have shed some light on cases of higherorder symmetry (such as the dodecahedral case), as well as nonsymmetric cases such as the five-cell pattern discussed. Although these patterns of convection are not expected to be found in the Earth, they can further aid the verification, validation and comparison of new numerical methods, algorithms, and codes, as applied to mantle convection in the Earth and other terrestrial planets. We also hope that this paper will stimulate further investigation into how the type and order of numerical discretization affects pattern formation in the context of benchmarking community codes.