There have been quite a lot of new developments in 3D spherical mantle convection models in the last few years, particularly in Japan and Europe. It seems to be a good time now to start doing some benchmarks for these codes. Benchmarks for spherical mantle convection were discussed at the Boulder mantle convection workshop in June of 2005. After the workshop, there were some email discussions on defining benchmark problems. This year, Uli Hansen's group published a paper in PEPI in which they included some benchmark calculations from their code and several other codes {Stemmer et al., 2006, PEPI}. The papers by Stemmer et al {2006}, Ratcliff et al. {1996, JGR}, Zhong et al. {2000, JGR}, and Yoshida and Kageyama {2004, GRL} and the previous discussions form a basis for benchmark calculations presented here and for discussions of our future benchmark studies.

Here we present some initial benchmark results for CitcomS (20 cases so far) and comparisons with Stemmer et al {2006} and Ratcliff et al {1996} and Yoshida and Kageyama {2004} when the solutions are available (9 of them) from these codes. We define benchmark cases and present averaged global quantities as well as local properties for benchmark cases. At the moment, all the cases are for entirely basal heating and steady state or weakly time-dependent results, and viscosity variations from temperature-dependence vary from 1 (mobile-lid regime) to 1e6 (stagnant lid regime). And $Ra_{0.5}$ defined by viscosity at $T=0.5$ varies from 7e3 to 1e5 (i.e., to achieve steady state or weakly time-dependent state).

We hope that the community will join our effort and discuss what other benchmark cases should be included in our future benchmarks (e.g., internal heating, time-dependent cases at relatively high Ra, …).

All the cases assume the Boussinesq approximation and constant material properties except for viscosity. There are no phase transitions, no compositional anomalies, and no internal heating (i.e., entirely basal heating) for all the cases. The top and bottom boundaries have radii $r_t=1$ and $r_b=0.55$, respectively.

The viscosity can be temperature-dependent and the non-dimensional viscosity, $\eta$, is given by

$\eta=exp[E(0.5-T)]$, (1)

where $E$ is the activation energy and $T$ is non-dimensional temperature. Rayleigh number is defined as

$Ra=\frac{\rho g \alpha \Delta T d^3}{\eta_{0.5} \kappa}$, (2)

where $d$ is the thickness of the mantle (i.e., the length scale) and $\eta_{0.5}$ is the viscosity at $T=0.5$. That is, Ra is what we sometimes call $Ra_{0.5}$. This also implies that time and velocity scalings are defined using the thickness of the mantle, which is somewhat different from those in CitcomS which uses planetary radius as the length scale. All the CitcomS results in this benchmark study are re-scaled to use the mantle thickness as length scaling, to be consistent with those used in other studies.

Boundary conditions are free-slip at the top and bottom boundaries and isothermal with non-dimensional temperatures of 0 and 1 at the top and bottom boundaries, respectively.

Initial conditions come in two types. While most cases use initial temperature that is given as a function of coordinates (i.e., perturbation at some given spherical harmonics superimposed on a conductive temperature profile), some other cases use temperature field from other calculations as initial temperature. For the first type of initial condition, the initial temperature is given by

$T(r, \theta, \phi) = r_b/(r_b-r_t)(1-r_t/r) + [\varepsilon_c cos(m \phi)+\varepsilon_s sin(m \phi)] p_{lm}(\theta, \phi) sin[\pi (r-r_b)/(r_t-r_b)]$, (3)

where $\varepsilon_c$ and $\varepsilon_s$ are the magnitudes of perturbation for cosine and sinine terms respectively, $l$ and $m$ are spherical harmonic degree and order, $r$, $\theta$, and $\phi$ are radial, co-latitude and longitude coordinates, and $p_{lm}$ is a normalized associated Legendre polynomial that is related to the associated Legendre polynomial $P_{lm}$ as:

$p_{lm} (\theta, \phi) = \sqrt{\frac{(2l+1)(l-m)!}{(1+\delta_{m0}) 2 \pi (l+m)!}} P_{lm}(\theta, \phi)$. (4)

The magnitude of the perturbations (i.e., $\varepsilon_c$ and $\varepsilon_s$) are 0.01, but for some cases, we only use the cosine term with $\varepsilon_s=0.0$. For some cases, perturbations can be given at two different sets of harmonics (e.g., $l=4$, $m=0$ and $l=4$, $m=4$ for all the cubic symmetry cases). Detailed initial conditions for each case will be discussed later.

We do not use random perturbations for benchmark studies here.

All the cases are computed to a steady state for global properties. Only steady state results are quantified. For some cases, we also give an output file that lists time-dependence of the results for better comparisons with other codes.

For every 20 timesteps, we compute Nusselt numbers for both the top and bottom boundaries, averaged temperature for the whole mantle or shell $<T>$, and averaged RMS velocity $<v_{rms}>$.

$<T> = \int_\Omega T d\Omega / \int_\Omega d\Omega$, (5)

$<v_{rms}> = [ \int_\Omega v^2 d\Omega / \int_\Omega d\Omega]^{1/2}$. (6)

When a case reaches a steady state, we then compute the time-averaged values and standard deviation of Nusselt numbers, averaged temperature, and averaged RMS velocity over a certain period of time.

For every 20 timesteps, we also compute the maximum and minimum radial velocity and temperature at the middle depth of the mantle (i.e., at $r=0.775$). Again, we compute their time-averaged values and standard deviations over a certain period of time after a steady state is reached.

These global and local properties are also used in Stemmer et al {2006} in their benchmark study.

Two sequences of cases are presented here, one associated with tetrahedral symmetry and the other with cubic symmetry. $Ra_{0.5}$ ranges from 7e3 to 1e5 and viscosity variations due to temperature-dependence ranges from 1 (i.e., isoviscous) to 1e6 in stagnant lid regime. These cases are similar to what Ratcliff et al {1996} and Stemmer et al. {2006} had presented. However, we added more cases with larger viscosity contrasts. For cases with $Ra_{0.5}=7e3$, we use 12x32x32x32 elements, and for cases with $Ra_{0.5}=1e5$, 12x48x48x48 elements are used. Redial resolution is refined near the top and bottom boundary layers.

In addition to those outputs in Tables, we also provide files that list the full time-dependence of our standard outputs. For some cases, thermal structure images are also provided.

For all the cases in this category, we use $Ra_{0.5}=7e3$, but these cases have different temperature dependent viscosity, $Dh$ of 1, 10, 20, 100, 1000, 1e4, 1e5, and 1e6 (Cases BM1A-BM1H in Table 1 for definitions of all cases. The initial condition is the same for these cases with $\varepsilon_c=\varepsilon_s=0.01$ in equation 3. The first three cases, BM1A, 1B, and 1C were also presented in Stemmer et al {2006}, Ratcliff et al. {1996}, and Yoshida and Kageyama {2004}; and Zhong et al. {2000} and Bercovici et al. {1989} computed case BM1A. Only the first five cases display tetrahedral symmetry. The last two cases, BM1G and 1H are in stagnant-lid regime, showing a large number of small plumes. Case BM1F with $Dh$ of 1e4 is a transitional case.

Table 2 shows the output results and comparisons with previous studies when available. Note that in addition to average values, Table 2 also gives the time duration over which the averages are computed and standard deviations. It is clear that for Cases BM1A, 1B and 1C, CitcomS' results compare well with Stemmer et al. {2006}, Ratcliff et al. {1996} and Yoshida and Kageyama {2004}.

Files for time-dependence of the outputs from $t=0$ to steady state and representative snapshots of residual temperature field for selective cases can be downloaded at Thermal Convection Benchmarks for CitcomS. The format of these output files is also explained there.

For cubic symmetry cases with $Ra_{0.5}=7e3$, we use the initial condition similar to Ratcliff et al. {1996}:

$T(r, \theta, \phi) = r_b/(r_b-r_t)(1-r_t/r) + \varepsilon_c [ p_{40}(\theta, \phi)+5/7 cos (4 \phi) p_{44} (\theta, \phi)] sin[\pi (r-r_b)/(r_t-r_b)]$, (7)

where $sin(m f)$ is excluded. We have computed cases with temperature dependent viscosity, $Dh$ of 1, 20, 30, 100, 1000, 1e4, 1e5 and 1e6 (Cases BM2A-BM2H in Table 1 for definitions of all cases). The first three cases, BM2A, 2B, and 2C, can again be compared with previous studies. The first six cases (BM2A-BM2F) display cubic symmetry with 6 upwelling plumes. Different from tetrahedral cases, now BM2F with $Dh$ of 1e4 also shows flow cubic symmetry, the same as cases with smaller viscosity contrast. The last two cases BM2G and BM2H are again in stagnant lid regime, with flow pattern breaking away from the cubic symmetry. Notice that the slight difference in Nu between BM1 and BM2 cases for the same Ra.

Table 3 shows the output results and comparisons with previous studies when available. Again, CitcomS compares well with all the previous studies for BM2A, BM2B and BM2C.

For cases with $Ra_{0.5}=1e5$, we have computed cases with temperature dependent viscosity, $Dh$ of 1, 10, 30, and 100, so far (cases BM3A-3D, in Table 1). We use the same initial condition as in equation 7 only for the case with uniform viscosity that leads to cubic symmetry flow. For other $Ra_{0.5}=1e5$ cases with variable viscosity, we use the final temperature field from either the uniform viscosity case (i.e., BM3A) or other cases as initial condition to help preserve the cubic symmetry flow. Table 1 lists the detailed initial conditions. We found that if we use equation 7 as initial conditions for these high Ra cases with variable viscosity, we can no long achieve cubic symmetry flow (at least after 25000 time steps for BM3B). More plumes develop from the CMB than expected from cubic symmetry and the plumes are stable for a long time.

Table 3 shows the results and comparisons with Ratcliff et al {1996} for cases BM3A, BM3B and BM3C for which Ratcliff et al {1996} also published RMS velocity and Nu results. While the results from CitcomS compares reasonably with Ratcliff et al {1996}, the deviations are significantly larger than those for small Ra cases presented earlier. We think that the difference is caused by the resolution. Ratcliff et al. {1996} used the same resolution (i.e., 32, 64, and 128 cells in radial, co-latitude and longitude directions, respectively, for cubic symmetry cases) for different Ra cases. We think that cases with $Ra_{0.5}=1e5$ should require higher resolution.

Again, files for time-dependence of the outputs from $t=0$ to steady state and representative snapshots of residual temperature field for selective cases can be downloaded at Thermal Convection Benchmarks for CitcomS.

Here we present results of Nussult numbers, RMS velocity, averaged temperature, and maximum and minimum flow velocity and temperature at the mid-mantle depth for 20 cases from CitcomS. For these cases, $Ra_{0.5}$ is either 7e3 or 1e5, and viscosity contrast varies from 1 to 1e6. The style of convection varies from mobile lid to stagnant lid regimes. For nine of the 20 cases, we could compare with previously published results {Stemmer et al., 2006; Ratcliff et al., 1996; Yoshida and Kageyama, 2004}. Comparisons show that CitcomS's results are generally consistent with these previous studies, although at high Ra, the relatively low resolution from previous studies may degrade the agreement. Time dependent results from CitcomS are compiled into a file for each case, and all of them can be downloaded for more detailed comparisons. Compared with benchmarks in Zhong et al. {2000}, here we added significantly more thermal convection cases.

We encourage our colleagues to send in their results for these cases. We think that more benchmark cases are needed to test time-dependent convection at even higher Ra, and internal heating convection. We welcome suggestions.

Acknowledgement: Mr. Nan Zhang, a graduate student at the University of Colorado, helped with calculations and analyses.