3D Spherical Mantle Convection Benchmarks
Presents initial benchmark results for CitcomS and comparisons with Stemmer et al {2006} and Ratcliff et al {1996} and Yoshida and Kageyama {2004} when the solutions are available from these codes. Comments and suggestions are requested.
1) Introduction
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 Ra0.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, …).
2) Definition of Benchmark Cases
2.1) General features
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 rt=1 and rb=0.55, respectively.
2.2) Viscosity, Rayleigh number and scalings.
The viscosity can be temperature-dependent and the non-dimensional viscosity, η, is given by
η=exp ( E(0.5−T) ), (1)
where E is the activation energy and T is non-dimensional temperature. Rayleigh number is defined as
Ra=(ρgαΔTd3)/(η0.5κ) , (2)
where d is the thickness of the mantle (i.e., the length scale) and η0.5 is the viscosity at T=0.5. That is, Ra is what we sometimes call Ra0.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.
2.3) Boundary conditions
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.
2.4) Initial conditions
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,θ,ϕ)=rb/(rb−rt)(1−rt/r)+( εccos(mϕ)+εssin(mϕ)]plm(θ,ϕ)sin(π(r−rb)/(rt−rb )]), (3)
where εc and εs are the magnitudes of perturbation for cosine and sine terms respectively, l and m are spherical harmonic degree and order, r, θ, and ϕ are radial, co-latitude and longitude coordinates, and plm is a normalized associated Legendre polynomial that is related to the associated Legendre polynomial Plm as:
plm(θ,ϕ)=√ ((2l+1)(l−m)! / (1+δm0)2π(l+m)! ) Plm(θ,ϕ). (4)
The magnitude of the perturbations (i.e., εc and εs) are 0.01, but for some cases, we only use the cosine term with ε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.
3) Quantifying Outputs of Benchmarks
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.
3.1) Global properties
For every 20 time steps, we compute Nusselt numbers for both the top and bottom boundaries, averaged temperature for the whole mantle or shell , and averaged RMS velocity < vrms >.
< T >=∫ΩTdΩ/∫ΩdΩ, (5)
< vrms >=( ∫Ωv2dΩ/∫ΩdΩ )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.
3.2) Local properties
For every 20 time steps, 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.
4) Results
Two sequences of cases are presented here, one associated with tetrahedral symmetry and the other with cubic symmetry. Ra0.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 Ra0.5=7e3, we use 12×32×32×32 elements, and for cases with Ra0.5=1e5, 12×48×48×48 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.
4.1) The cases with tetrahedral symmetry and its variations
For all the cases in this category, we use Ra0.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 εc=ε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.
4.2) The cases with cubic symmetry and its variations
For cubic symmetry cases with Ra0.5=7e3, we use the initial condition similar to Ratcliff et al. {1996}:
T(r,θ,ϕ)=rb/(rb−rt)(1−rt/r)+εc ( p40(θ,ϕ)+5/7cos(4ϕ)p44(θ,ϕ) ) sin(π(r−rb)/(rt−rb)), (7)
where sin(mf) 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 Ra0.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 Ra0.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 Ra0.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.
5) Conclusions
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, Ra0.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.