ŒDIPUS: a new tool to study the dynamics of planetary interiors (original) (raw)

Journal Article

,

Search for other works by this author on:

,

2Department of Geophysics, Faculty of Mathematics and Physics, Charles University, Prague, Czech Republic

Search for other works by this author on:

,

Search for other works by this author on:

Search for other works by this author on:

Received:

29 November 2006

Accepted:

23 February 2007

• PDF
• Split View
• • Article contents
  • Figures & tables
  • Video
  • Audio
  • Supplementary Data
• Cite

Cite

G. Choblet, O. Čadek, F. Couturier, C. Dumoulin, ŒDIPUS: a new tool to study the dynamics of planetary interiors, Geophysical Journal International, Volume 170, Issue 1, July 2007, Pages 9–30, https://doi.org/10.1111/j.1365-246X.2007.03419.x
Close

Navbar Search Filter Mobile Enter search term Search

Summary

We present a new numerical method to describe the internal dynamics of planetary mantles through the coupling of a dynamic model with the prediction of geoid and surface topography. Our tool is based on the simulation of thermal convection with variable viscosity in a spherical shell with a finite-volume formulation. The grid mesh is based on the ‘cubed sphere’ technique that divides the shell into six identical blocks. An investigation of various numerical advection schemes is proposed: we opted for a high-resolution, flux-limiter method. Benchmarks of thermal convection are then presented on steady-state tetrahedral and cubic solutions and time-dependent cases with a good agreement with the few recent programs developed to solve this problem.

A dimensionless framework is proposed for the calculation of geoid and topography introducing two dimensionless numbers: such a formulation provides a good basis for the systematic study of the geoid and surface dynamic topography associated to the convection calculations. The evaluation of geoid and surface dynamic topography from the gridded data is performed in the spectral domain. The flow solver is then tested extensively against a precise spectral program, producing response functions for geoid as well as bottom and surface topographies. For a grid mesh of a reasonable size (6 × 64 × 64 × 64) a very good agreement (to within ∼1 per cent) is found up to spherical harmonic degree 15.

1 Introduction

Several of today-s questions in the study of planetary interiors imply the use of models simulating thermal convection in a spherical geometry. This is either implicitly required by the internal dynamics (effect of curvature on thermal instabilities and heat flux budget, correct scaling of mass flux when addressing compositional gradients, modelling plate motion on a sphere) or due to a better comparison of these with geophysical data often expressed as a decomposition into spherical harmonics (e.g. in the case of data from space exploration). In addition, in order to address many problems a numerical model of thermal convection should handle large viscosity gradients: (i) the effective viscosity of silicates or ices is very sensitive to thermodynamical parameters, inducing very specific dynamical regimes such as the asymptotic conductive lid regime (e.g. Davaille & Jaupart 1993); (ii) in order to mimic the plate behaviour, more complex rheologies have been proposed that imply even larger viscosity gradients (e.g. Tackley 2000a,b; see also Bercovici 2003 for a review).

The first programs simulating thermal convection in a spherical shell were based on a spectral approach (Glatzmaier 1988; Zhang 1995) since this allows an elegant treatment in the case of a viscosity field that is constant or radius-dependent. However, these methods lose part of their interest in the case of large lateral viscosity variations. A first attempt to solve variable viscosity convection with a grid-based method is proposed by Hsui et al. (1995) after a finite element discretization developed in the isoviscous case (program Terra) by Baumgardner (1985). Another technique based on spherical coordinates is proposed by Ratcliff et al. (1996). To the authors' knowledge, the first robust program solving Stokes flow and thermal convection in a spherical shell for a fluid with strongly variable viscosity (CITCOMs) is then presented by Zhong et al. (2000). Based on the finite element program CITCOM (Moresi & Gurnis 1996), this version uses brick elements (eight nodes) and the sphere is divided in 12 blocks of approximately equal size. This mesh is well suited for parallel computing and the FEM formulation ensures the robustness of the program. Though it seems necessary to have several different tools to solve this problem (especially when spherical shell laboratory experiments are not a simple set-up to build), new numerical solutions have not been proposed until very recently. In the framework of finite volumes or finite differences methods, the necessity to propose complex non-conform meshes to solve the pole problem for the spherical geometry induced a relatively heavy computational treatment. Alternatively, the more precise and more ‘elegant’ spectral formulation for such a geometry (through the use of spherical harmonic functions) loses a significant part of its interest when azimuthal viscosity variations are introduced. Three finite volume formulations have been proposed recently (Kageyama & Sato 2004; Choblet 2005; Stemmer et al. 2006) all resulting mainly from a new geometrical treatment. The first one, based on the ‘Yin Yang’ grid uses the classical spherical coordinate system and the poles are excluded. The two latter ones use a decomposition of the sphere into six blocks obtained through the projection of a cube, the ‘cubed sphere’ decomposition (while among the five platonic solids, Terra used the icosahedrons) and either formulate the equations directly in the natural (but non-orthogonal) coordinate system (Choblet 2005) or uses the Cartesian coordinates of a ‘smoothed cubed sphere’ mesh where the non-orthogonal aspect is reduced iteratively (Harder & Hansen 2005; Stemmer et al. 2006).

Observational constraints on the style of thermal convection in planetary bodies are still rare. Two potentially important sources of information about the internal dynamics of a planet are the geoid and topography. In a viscous mantle, the stresses arising from mantle flow deform all density interfaces, including the surface. The pattern of such a ‘dynamic’ topography and its magnitude bear evidence about the viscosity structure of the planet, the mass flux through the major internal interfaces and the characteristic wavelengths of mantle convection (e.g. Schubert et al. 2001). Important information about mantle convection can also be contained in the observed gravitational field which is a superposition of the gravitational signal due to density heterogeneities inside the body and the signals of density interfaces [surface, core-mantle boundary (CMB), etc.] which may be deflected by mantle flow (Ricard et al. 1984; Hager & Clayton 1989; King 1995). Prediction of the long-wavelength component of the geoid requires the integration of mass anomalies over the whole body and can only be done in a spherical geometry. This integration can significantly be simplified if the density anomalies and dynamic topographies are expressed in terms of series of spherical harmonic functions (Burša & Pec 1993). That is why the momentum equation has often been solved in a spectral domain which allows the equations governing the flow to be directly coupled with the Laplace-Poisson equation for the gravitational potential (Ricard et al. 1984; Hager & Clayton 1989). Unfortunately, as already mentioned, the spectral methods are not suitable in the case of models with large lateral viscosity variations (Cadek & Fleitout 2003). In this paper, we therefore, combine the solution of the Navier-Stokes equation in a finite volume formulation with a spectral method used for evaluating the dynamic topography and the geoid.

Origin, Evolution and Dynamics of the Interiors of Planets Using Simulation (ŒDIPUS) may seem a rather grandiloquent acronym for the tool we present in this paper (may the authors be forgiven for this!): as indicated in Section 2 we so far focus on the purely thermal convection aspect (Section 2.1) with temperature-dependent rheology (Section 2.2) of planetary interior dynamics. Our final objective is however to include the compositional aspect in the program (Section 2.3). The numerical methods are described in Section 3: a new version of the ‘cubed sphere’ mesh is presented (Section 3.1), the flow solver, already described by Choblet (2005), is shortly evocated (Section 3.2) and the advection scheme, problematic in this non-orthogonal coordinate system (see the remarks by Hernlund & Tackley 2003) is discussed (Section 3.3). For the sake of clarity, we describe the treatment of the geoid and the dynamic topography in an independent Section 4 since the formalism (non-canonical dimensionless problem) and the numerical method (based on a decomposition in the spectral domain) are specific. The last Section 5 compiles tests we performed to ensure the validity of the method, concerning the flow solver, thermal convection and our computation of geoid and dynamic topography.

2 Physical Model

2.1 Thermal convection

We solve the classical conservation equations describing thermal convection of a variable viscosity, infinite Prandtl number fluid in the Boussinesq approximation (aside from the viscosity variations) here expressed in their dimensionless form:

1

2

3

Ra = (αρ0_g_Δ_T d_3)/(κμ0) is the surface Rayleigh number and is the dimensionless heating rate (in the whole paper, • denotes the dimensionless variable while stands for the variable with a physical dimension; see Table 1 for the notations). Note that dimensionless length variables (mostly r) are scaled using d (and not the radius of the outer layer): dimensionless radii thus vary from r b = f/(1 −f) to r t = 1/(1 −f), where is the ratio between the radii of the inner () and outer (, radius of the planet) shells.

Table 1.

Variables and symbols.

2.2 Rheology

The temperature dependence of viscosity is modelled in this study in terms of a simple exponential function:

4

where γ, the viscosity parameter, is the fourth dimensionless parameter introduced in the study of variable viscosity convection in a spherical shell (besides Ra, h and f). Relationship (4) is often used in mantle-convection studies so that a comparison with earlier results is easier. Such an approximation of the Arrhenius law () considered to describe the effective viscosity of rocks and ices due to solid-state thermally activated creep is termed the Frank-Kamenetskii approximation (Frank-Kamenetskii 1969). In the case of the stagnant lid regime, most of the viscosity variations occur in the cold conductive layer at the surface. Thus, if the two functions are scaled properly, heat transfer and flow are not much affected by the strong difference in the cold conductive region (Reese et al. 1999). A more important point in the context of the present study concerns stress in the conductive lid (where the two viscosity laws differ significantly) since it is used to compute geoid and dynamic topography (see Section 4.1). A recent study by Solomatov (2004) indicates that the stress distribution in the lid is similar for the two laws except within a thin layer beneath the surface.

As mentioned earlier, more realistic aspects can certainly be proposed when dealing with planetary dynamics in terms of rheology: the stress or grain size dependence, the influence of water (among other volatiles) or melt fraction, and the necessity to incorporate more exotic rheologies to produce plate-like behaviour. However, these complexities are mostly conceptual ones: the aim of the present study being to describe the numerical tool, all the above complexities may possibly be described by our numerical method as long as it handles large enough viscosity gradients (with the notable exception of anisotropic viscosity).

2.3 (No) compositional dynamics

It is now well established that density and viscosity gradients of compositional origin in the Earth-s mantle are likely to play a significant role (i.e. as important as the thermal ones) in its dynamics, at least in two crucial regions: the lithosphere (both continental, e.g. van Gerven et al. 2004 and oceanic, e.g. Choblet & Parmentier, 2001) and the deepest part of the lower mantle (e.g. Ishii & Tromp 2004; Deschamps & Trampert 2004). It is also believed that most of the lunar interior dynamics subsequent to the freezing of a primordial magma ocean are mostly shaped by compositional stratification and overturn (e.g. Hess & Parmentier 2005). Moreover, some elements such as radioactive elements and volatiles in silicates and impurities in ices, even if they do not affect density nor viscosity due to very small concentrations, play an important role in the study of the evolution of planets. A recent good example is the enigma of the methane content of Titan-s atmosphere, maybe of internal origin (Tobie et al. 2006). All these effects have direct consequences on the geoid and topography.

For all these reasons, it seems important that dynamic models include the possibility of the transport of chemical species through the pure advection of a scalar field with the flow (this would also enable the description of the transport of physical quantities such as grain size). Such a program simulating thermochemical convection in a spherical shell is developed for the first time by McNamara (2004). Though the advection scheme proposed in this tool is built to minimize numerical diffusion (see Section 3.3), even this small diffusion might not be acceptable when the real diffusion is negligible. In the framework of our program, we believe that the most appropriate method may be the tracers method described by Tackley & King (2003) rather than more complex ‘surface reconstruction methods’. Indeed, though the grid mesh is not strictly uniform and though the coordinates system is not orthogonal, both the algorithm proposed by Tackley & King (2003) in a cartesian context and the properties of the method are conserved.

3 Numerical Method

In this section, we mostly refer to the finite volume method used to solve the thermal convection problem. The mesh (Section 3.1) and the techniques used in the flow solver (Section 3.2) are already described in an earlier paper (Choblet 2005). New aspects developed here are mostly a new orientation of the coordinate system (cf.Section 3.1), a benchmark of advection schemes in the cubed sphere context (Section 3.3) and the communication between blocks (cf.Appendix B).

3.1 Mesh

The cubed sphere method proposed by Ronchi et al. (1996) is first designed to avoid coordinate singularities at the poles. Other interesting aspects of the method are: (i) the fact that the physical grid defined on the spherical surface is reasonably close to a uniform grid, (ii) the decomposition of the spherical shell into six identical regions, with the same metric and (iii) an efficient parallel implementation. The mapping of the spherical surface is obtained by the projection of the sides of a circumscribed cube onto this surface. Fig. 1 illustrates this technique and introduces the cubed sphere coordinate system (r, ξ, η). Note that it is not an orthogonal system. In addition, it appears that angular coordinates ξ and η are such that a symmetry between them is expected in all equations. We explicitly write the equations in this coordinate system in the program, thus avoiding the computation of the metric tensor elements during the numerical computation. As proposed by Ronchi et al. (1996), the following auxiliary variables are introduced in order to simplify the expressions of differential operators:

5

This surface is then mapped onto a rectangular grid, where the equations are solved using a standard finite volumes technique designed for regular meshes. The gridding technique can be illustrated by two sets of great circles spaced with a uniform angular increment (Δξ and Δη) with η and η varying between −π/4 and π/4. The mesh is then extended radially between the two spherical surfaces bounding the shell (r b) and (r t), by concentric spheres with identical angular discretization, spaced by a radius increment Δ_r_ (see Fig. 2). Note that if the mesh associated to one of the blocks displays the same number of cells in the three directions, the azimuthal extent is usually larger than the radial dimension, assuming f is not too small: (π/2(1 −f) times smaller on surface cells and π_f_/2(1 −f) times on the cells adjacent to the inner boundary). This is often the configuration used in this study's calculations since radial gradients play an important role in thermal convection.

Figure 1.

The sphere decomposed into six blocks and the curvilinear coordinate system on block 0. The main difference with spherical coordinates is that curves of constant values of the second coordinate η are great circles (or ‘horizontal meridians’) and not parallels.

Figure 2.

The grid mesh as six ‘cubic’ blocks. The azimuthal dimension of cells varies linerarly with the radius.

The curved numerical domain is thus divided into non-uniform cells C ijk(i, j, k corresponding to the indexes in the _r_-, ξ- and η-directions, respectively). The staggered grid mesh approach is adopted and the choice is made to use contravariant coordinates providing a better evaluation of fluxes for this configuration (see Choblet 2005, for further details). For technical reasons, we modified the coordinates orientation on blocks 2, 3 and 4. If there are N i grid cells in a given direction, for example, direction ξ, there are N _i_+ 1 discrete values of the contravariant coordinate of the velocity field _V_ξ between ξ = −π/4 (V_ξ_j =0) and ). Thus, V_ξ_j =0 on block 1 will be computed from the value on block 0 (and also value _V_η (ξ = π/4) on block 0, due to non-orthogonality). In order to write strictly the same algorithm on each block, it is thus necessary that boundaries between blocks involve coordinates values of both -π/4 and π/4 (while, with the classical orientation proposed by Ronchi et al. 1996, coordinate η = π/4 on block 1 coincides with coordinate ξ = π/4 on block 4…). The new orientation satisfies this aspect. Relationships between this new coordinate system and the Cartesian and spherical basis are indicated in Table 2.

Table 2.

The modified ‘cubed sphere’ coordinate system.

3.2 Flow solver

In the cubed sphere coordinate system, eqs (1) and (2) provide, respectively, the following projections:

6

and, along r,

7

along ξ,

8

along η,

9

where the components t•• of the deviatoric part of the stress tensor (π = −pI+τ) are

10

11

12

13

14

15

Note that non-orthogonality between ξ and η leads to longer than usual expressions for the classical differential operators. The expressions of the discrete equations corresponding to [eqs (6)–(15)](#m6 m15) can be found in Choblet (2005).

The flow solver is decoupled from the treatment of the energy equation: the velocity field computed at time step n is used for the advection of heat in the computation of temperature at time _n_+ 1. The pressure and velocity fields are then updated simultaneously using the new temperature field to compute buoyancy and viscosity. We use a FAS-multigrid algorithm with V-cycles, a Gauss-Seidel smoother, linear transfer operators and coarse grid operators on all grids derived from the fine grid discretization (see Choblet 2005, for further details).

3.3 Advection-diffusion scheme

The equation describing conservation of energy is in the dimensionless form:

16

A classical drawback of grid based methods used to solve the advection problem is that discrete schemes usually provide solutions that satisfy better a different partial differential equation than the simple advective transport eq. (16): this modified equation includes diffusive (for first order schemes) or dispersive (second order schemes) terms. Moreover, in the context of non-orthogonal coordinate systems, the numerical diffusion is anisotropic: for the cubed-sphere coordinates, it is larger along ξ and η directions (e.g. Hernlund & Tackley 2003). In order to solve these problems, we have tested three advection methods in the ‘cubed sphere’ context described in Appendix A: (i) the classical first order ‘upwind’ method produces large amounts of numerical diffusion and it is only used as a reference scheme, (ii) the ‘MPDATA’ method, used in several programs simulating thermal convection in a Cartesian context to study mantle dynamics (e.g. Tackley 1994) and (iii) the ‘van Leer’ method used more commonly in the context of compressible fluids (van Leer 1979).

The test we propose is based on a solid rotation velocity field around the _x_-axis, that is, proportional to

17

We initially impose a sphere of radius 0.1 with an anomalous temperature 1 within a 0 temperature field in block 0 of the cubed sphere. It is then advected with the rotation velocity field so that if diffusion were negligible, the field should be invariant. We measure numerical diffusion through the volume of the region bounded by isotherm (T = 0.2) (Fig. 3). Fig. 4 displays a comparison of the different schemes using this criterion: increasing the number of correction steps (iord) with ‘antidiffusive’ velocities in the MPDATA scheme provides better results. Even with iord = 2, the volume of the hot sphere remains smaller on a mesh including 323 grid cells than the reference upwind scheme with a twice larger resolution (643 grid cells). Tests of the van Leer methods, whatever slope is used, return even better results: the volume obtained with slope 2 (geometrical average, see Appendix A) after a similar time is twice smaller than the one obtained for iord = 3. Slope 1 (MC limiter, see Appendix A) leads to a volume that increases less (40 per cent smaller than the increase obtained for slope 2). Even though the treatment of the energy equation corresponds only to a small fraction of the global costs associated to the simulation of thermal convection with our method, the CPU cost of the various advection methods is displayed in Table 3. The van Leer methods are both twice more time-consuming than the simple upwind scheme. Increasing the number of corrective steps in MPDATA increases the computation time (iord = 4 already costs twice as much as the van Leer methods). For all these reasons, we will use the quite simple and quite efficient van Leer method (with MC limiter slope 1) in the following computations.

Figure 3.

Structure of the modified ‘cubed sphere’ coordinate system on other blocks: the original orientation of the curvilinear coordinates is modified due to the array structure of velocity fields in a staggered grid mesh approach (see text for further details).

Figure 4.

Comparison of the different advection schemes: normalized volume delimited by the 0.2 isotherm as a function of time. Regular lines refer to a 323 grid mesh for block 0, bold lines refer to a 643 grid mesh. The lines without symbols are associated to the reference upwind scheme. Squares refer to the MPDATA scheme with increasing values of the number of ‘antidiffusive’ steps associated to the shade: iord = 2 (black), iord = 3 (grey), iord = 4 (white). Circles refer to the van-Leer method with either slope 1 (white circles) or slope 2 (black circles), see Appendix A for further details on the slopes. Time 0.314 approximately corresponds to ten rotation cycles for the velocity field we chose.

Table 3.

CPU cost of the advection schemes.

Further remarks on the numerical diffusion process are illustrated on Fig. 5. While the reference upwind scheme provides a strong flattening of the volume compared to the initial spherical shape, this feature is attenuated in the van Leer case. The main reason for this flattening is that the velocity field we employ for the test has a zero V r component. It should thus be noted that this flattening does not reflect anisotropy in the numerical diffusion process. Anisotropy can be observed on the ‘azimuthal’ projections (right in Fig. 5) where a significant distortion exists for the upwind scheme while the projection of the volume remains circular in the van Leer case (coordinates ξ and η playing a role that should be strictly equivalent in this test).

Figure 5.

Volume of the T = 0.2 isotherm after one full rotation on a 643 grid mesh. Upper panel: results for the upwind difference scheme. Lower panel: results for the van Leer method with slope 2.

4 Evaluation Of The Geoid And Dynamic Topography

In this section, we give basic formulas relating the stresses due to thermal convection in planetary interiors with surface observables, namely the dynamic topography and the geoid. We will follow the approach described by Ricard et al. (1984) and Richards & Hager (1984). For the sake of clarity, we first present these equations using physical quantities (denoted with tildes) and then we rewrite them in dimensionless forms, consistent with descriptions considered in previous sections.

4.1 Equations

The stresses arising from the viscous flow in the mantle deform the surface of the planet and the bottom boundary (sometimes referred to in this text as the CMB, as is the case for the Earth and most terrestrial planets). The deviation of the deformed surface from a spherical shape,

18

where is the radius of the surface at a point with spherical coordinates θ and ϕ, and is the mean radius of the planet, can be obtained from the equilibrium of radial forces evaluated on a sphere of radius . In a first approximation, we can write

19

where is the radial component of the traction force at a spherical boundary of radius ,

20

is its average value,

21

and g is the gravitational acceleration which is supposed here to be constant throughout the shell. The topography induced by mantle flow is usually referred to as dynamic topography. In evaluating the forces due to topographic masses (the right-hand side of eq. 19), we only considered the lithostatic pressure and neglected the effect of self-gravitation. Eq. (19) would have a more complex form if the viscous mantle is overlaid by an elastic or viscoelastic layer, as in Zhong (2002), or a membrane (e.g. Ribe 1992). A similar equation of equilibrium for the surface can be written for the CMB. If we neglect the variations of pressure in the liquid core, associated with gravitational effects of density anomalies in the mantle and the surface and CMB topographic masses, we can express the topography of the CMB as follows:

22

where Δρ_c_ is the density jump at the CMB.

The anomalies of gravitational potential induced by thermal convection in the shell are a superposition of the gravitational signals from topographies and , and the gravitational signal induced by density anomalies in the shell. Using the integral expression for gravitational potential and eqs (19) and (22), we obtain

23

where G is the gravitational constant, are the density anomalies in the mantle, and is the distance between the points and . If potential is known we can determine the height of the geoid from the Bruns theorem (Lambeck 1988),

24

We will show in Section 4.2 that the evaluation of integrals in eq. (23) can easily be done in spectral domain.

Comparison of the predicted geoid and dynamic topography with the observed data requires the predicted quantities to be expressed in physical units. With the characteristic scales defined above, the radial traction and density fields used in the computation of geoid and topography are obtained from results of the thermal convection model with the following relationship: . In this context, the dimensionless dynamic topography is

25

where dimensionless number E = αΔ_T_, measuring thermal expansion, is introduced. The dimensionless version of eq. (23) is

26

The dimensionless height of the geoid is then

27

where G*= G_ρ0_d/g is a gravitational dimensionless parameter. The scales introduced here are determined by the scaling of the internal dynamics problem ([eqs 1–3](#m1 m3)). Therefore, their values are not really significant for the geoid and the dynamic topography: for example, length scale d (thickness of the convective shell) is significantly larger than the characteristic values of or

4.2 Spectral evaluation

The computation of the surface dynamic topography, eq. (19), requires the radical component of the traction σ_rr_ at the upper boundary to be evaluated from a dimensionless version of eq. (20). This quantity is directly computed on the grid described in Section 3.1. For many applications, however, the surface topography needs to be expressed in terms of a spherical harmonic series. In this paper, for the sake of simplicity of notation, we use the complex spherical harmonic basis {Y_ℓ_m} which is orthonormal over the unit sphere:

28

where the asterisk denotes complex conjugation and δ_ij_ is the Kronecker delta. The amplitude t s of the dynamic topography at a point {θ, ϕ} can then be approximated as follows,

29

where ℓmax is the maximum degree considered and (t s)ℓ_m_ are the spherical harmonic coefficients of degree ℓ and order m,

30

(for more details, see Jones 1985 and Varshalovich et al. 1988). If t s is given at discrete points {θ_i_, ϕ_i_}, i = 1, 2, …, N g, the integral in eq. (30) can be replaced by a sum and the coefficients (t s)ℓ_m_ can be evaluated using the following simple formula:

31

where Δ_s_ i is the surface area corresponding to the _i_-th point, normalized so that . Eq. (25) thus reduces to

32

Note that the summation in eq. (29) is done from degree 1 to ℓmax. The degree-zero term corresponds to the mean value of t s over the sphere and is by definition equal to zero. The choice of the cut-off degree ℓmax in eq. (29) depends on the density of gridpoints. Since only the spherical harmonics whose half-wavelengths are greater than the size of the maximum grid step can be retrieved from the grid data, ℓmax must satisfy the condition

33

where λ_i_ is the grid step.

The calculation of the geoid, eq. (27), requires the evaluation of the volume and surface integrals in eq. (26) for gravitational potential. This integration is usually also carried out in spectral domain. First, we express all the quantities appearing in eq. (26) in terms of spherical harmonic expansions:

34

where X stands for U, T and evaluated on both interfaces (r = r t and r = r b). The terms of degree zero are not included in summation as we assume that all quantities considered represent deviations from average values. Substituting eq. (34) into (26), using the addition theorem (e.g. Burša & Pĕč 1993), and invoking the orthonormality property (28), we obtain,

35

where U_ℓ_m are the spherical harmonic coefficients of the external (r ≥ r t) gravitational potential. The spherical harmonic coefficients of the geoid, eq. (24), are then simply

36

5 Tests And Benchmark

5.1 Flow solver and response functions

The flow solver was tested successfully against analytical solutions for both an isoviscous fluid and a fluid with radius-dependent viscosity (Choblet 2005). These tests are based on the differential rotations of the bounding spheres.

To further test the flow solver, we have computed response functions, that is, the topography and geoid produced by a normalized density anomaly located at a varying depth (for details, see, e.g. Ricard et al. 1984; Hager & Clayton 1989) and compared our results with those obtained with a spectral method described by Čadek (1989). We thus solve only eqs (1), (2), (32), (35) and (36): eq. (3) describing the conservation of energy is not considered. Boundary conditions are no-slip on top and free-slip on the bottom. We use the following values of the dimensionless parameters:

37

The dimensionless temperature field is

38

Three viscosity functions have been considered. The first case is isoviscous, and two radius-dependent viscosities are tested (Fig. 6): profile 1 mimics the radial variations observed in a stagnant lid regime with a range of approximately four orders of magnitude while profile 2 is a simple step-like viscosity variation. Note that we do not introduce azimuthal gradients in the tested viscosity field because the spectral method we compare our results with, would lose its interest in this case.

Figure 6.

Viscosity profiles used for the ‘response function’ benchmark of the flow solver: the solid line curve is test profile 1 - the dotted step-like curve is test profile 2.

We use both the grid-based flow solver described above (cf.Section 3.2, the grid mesh includes 6 × 64 × 64 × 64 cells corresponding to the discretization used in classical thermal convection calculations) and a solver based on a spectral method in order to produce the velocity and pressure fields that are then used to compute dimensionless values of the dynamic topography and geoid (eqs 32 and 36). For each couple (ℓ, m), dynamic topographies at the surface, t s, and at the CMB, t c, as well as the geoid N are evaluated as a function of the mean depth of the density anomaly. We then compute the spherical harmonic coefficients (t s)ℓ_m_, (t c)ℓ_m_ and N_ℓ_m using the integrations indicated in Section 4.2. Besides computing the spherical harmonic coefficients at the same degree as the load, we have also ensured that the coefficients of geoid and topographies at other degrees and orders effectively equal zero. We have also tested that the response function at degree ℓ are identical for any m ≤ℓ.

Figs 7–9 display some results of these tests obtained for degrees 2, 5, 10 and 15. Concerning the absolute values of the response functions, the isoviscous case displays the classical result that increasing degree ℓ results in a sensibility to shallower (for surface dynamic topography and geoid) or deeper (for bottom topography) structures and also in a smaller amplitude. It is also interesting to note that the effect of the strong and continuous viscosity increase, mimicking a stagnant lid but with only radial viscosity variations (profile 1, flat response function for low degrees between r = 1.8 and 2 for topography, cf.Fig. 8, especially for low degrees) is also observed for profile 2 (Fig. 9) where the viscosity jump across the ‘lid’ is much smaller and discontinuous (30). Finally, the results indicate how the response functions reflect the radial derivative of the viscosity profile (profile 1 is smooth, profile 2 is discontinuous).

Figure 7.

Response functions obtained by the spectral method (solid line) and by the ŒDIPUS solver (open circles) and the associated error. On the left panel, dimensionless geoid, on the right panel, dimensionless dynamic topographies (blue: ‘CMB’ topography, red: surface topography). The dashed line with the same colour indicates the ‘error’ (i.e. the absolute difference between two values at a given radius). Note that the scale for this difference is 100 times larger than the one corresponding to the related variable in order to represent values of about a percent of the maximum value. Four spherical harmonic degrees are presented: l = 2, 5, 10, 15. In this first case, the material is isoviscous.

Figure 8.

Same than Fig. 7 for curved viscosity profile 1 (cf.Fig. 6).

Figure 9.

Same than Fig. 7 for step-like viscosity profile 2 (cf.Fig. 6).

Concerning the differences between the two methods (‘error’ variable in Figs 7–9), the spectral method is considered to be the reference due to its better precision, ‘error’ thus mostly refer to the grid-based flow solver. In the isoviscous case (cf.Fig. 7), both results are in very good agreement. The maximum discrepancy is always observed for the largest degree (here ℓ = 15). For dynamic topographies, it is located near the two interfaces and reaches a value of 1 per cent. These are responsible for the surface error on geoid (also 1 per cent of the peak value located at _r_∼ 1.9). These errors are likely to be caused by the discrete treatment of boundary conditions. In the case of profile 1 (Fig. 8), the shape of the error profiles is affected by the viscosity profile. This is also obvious in the case of step-like profile 2 (Fig. 9): the two discontinuities in the viscosity profile induce discrepancies between the two results above and below. This leads to 1–1.5 per cent errors even for lower degrees (e.g. for ℓ = 5 between r = 1.4 and 1.6). Finally, we consider that the very good agreement between these results proves that ŒDIPUS provides reliable estimations of geoid and topography for planetary applications where degrees ℓ≤ 15 will be considered. It should be noted that for the present size of the grid mesh, the results deteriorate strongly when degrees ℓ≥ 20 are studied. Hence, a refinement of the grid mesh (at least in the azimuthal direction, 6 × 64 × 128 × 128 cells) is needed to reduce the discrepancy to acceptable values (less than 5 per cent for geoid). Eq. (33) indicates that, in the case where each block of the cubed sphere is associated to a _N_2 mesh (with N = _N_ξ = _N_η), degrees up to 2_N_− 1 could be considered (ℓmax = 127 if N = 64, and only N = 11 cells would be needed for ℓmax = 20…). This relationship clearly overestimates the maximum reliable degree for a given ‘grid size’ in the angular directions. Instead, our tests indicate that approximately eight times more grid cells are needed. An empirical rule would thus be that each block of the cubed sphere is divided in N = 4 ×ℓmax cells in each angular direction if spherical harmonics up to degree ℓmax are to be considered with confidence. A more rigorous approach, however, is that each time higher degrees of topography and geoid are computed, a short benchmark procedure on response functions up to this degree be made for the mesh.

5.2 Thermal convection

Benchmarks for thermal convection have already been produced in the earlier publication (Choblet 2005) including tests on the onset of convection for spherical harmonics perturbations and tests on the value of the Nusselt number for stationary convection at relatively low Rayleigh numbers. However, since we now use a different numerical method to simulate advection that is proven to be more accurate (e.g. it produces less numerical diffusion, cf.Section 3.3), we present four results referenced in Table 4. Table 5 displays a new benchmark based on recent compilations (Yoshida & Kageyama 2004; Stemmer et al. 2006): two cases (numbered 1 and 2) have been isolated among the list, one corresponding to a tetrahedral solution, one to a cubic one, these two steady platforms being identified by linear stability analysis and obtained by several numerical models for relatively low values of the Rayleigh number.

Table 4.

List of the presented calculations with the values of the dimensionless parameters: ratio between the two radii bounding the spherical shell (f), surface Rayleigh number (Ra), exponential viscosity parameter (γ), dimensionless volumetric heating (h), mechanical boundary condition on the two spherical interfaces. The comment in the last column refers to the shape of the initial temperature pertubation for the two first cases.

Table 5.

Benchmark tests for an isoviscous fluid. For the two runs, f = 0.55 (close to the Earth's mantle curvature). Case 1: Ra = 7 × 103, tetrahedral solution - Case 2: Ra = 3.5 × 103, cubic solution. The code in the first column refers to the models used for comparison: B89 (Bercovici et al.1989), H96 (Harder & Christensen 1996), R96 (Ratcliff et al. 1996), Z00 (Zhong et al. 2000), Y04 (Yoshida & Kageyama 2004), S06 (Stemmer et al. 2006), C06 is the present study. The number of cells in the radial direction (N r) as well as the approximate mesh size used to produce the results appear in column 2. Surface (_Nu_top) and bottom (_Nu_bot) Nusselt numbers are indicated as well as the rms of the velocity field (_V_rms).

We used a 6 × 32 × 32 × 32 cells grid mesh to perform both calculations. This mesh provides a satisfactory resolution for these cases with small Rayleigh numbers. Mechanical boundary conditions are stress-free on both spherical surfaces bounding the shell. The results show a good agreement with solutions obtained earlier with various numerical methods. A comforting observation is also the small difference between bottom and surfaces heat fluxes that indicates that the scheme is indeed globally conservative. The techniques used in the various programs compared here seems to have a small effect on global scalar values such as the Nusselt number. Part of the scattering might be attributable to the variations of the mesh resolution (a precise analysis is proposed by Stemmer et al. 2006); the specific treatment of numerical diffusion proposed in this study (cf.Section 3.3) is not observable. Global averages are probably not a sufficient criterion for a very precise benchmark. In addition, the use of higher order accuracy schemes in time would also produce variations among the various programs if a time-dependent case were studied (our explicit method is only order 1 in time). A precise benchmark in the future should also include such a test.

We further propose two time-dependent results for thermal convection. One (Case 3) is obtained for an isoviscous fluid with purely basal heating, the other (Case 4) corresponds to a stagnant-lid regime and includes approximately 20 per cent volumetric heating. Both calculations were performed with no-slip radial boundaries and both have reached a statistical steady state on a 6 × 64 × 64 × 64 cells grid mesh. A typical duration for such a run is about a 100 hr on 48 (rather slow, 500 MHz) SGI processors R14000. Time-averaged profiles are presented in Fig. 10 (with therefore, a relatively smooth aspect). A classical result is the asymmetry between the cold and hot boundary layers due to their different curvature. If both boundary layers had a similar thickness, the temperature difference across the hot one would be approximately 1/(1 +_f_2) and _f_2/(1 +_f_2) across the cold one. The temperature of the well-mixed interior is thus linked to the ratio between the areas of both surfaces (in this case, this temperature is indeed close to 1/5, cf.Fig. 10 and Table 6). Both the introduction of temperature-dependent viscosity and the addition of volumetric heating naturally increases this value (cf.Fig. 10). The minimum and maximum temperature profiles also indicate that in both cases, the average is closer to the minimum value than to the maximum. In Case 3, the whole layer participates to convective motion: the difference between the minimum temperature and the average value at a given depth thus simply reflects the temperature difference across the cold boundary layer; the difference between the maximum and the average reflects the difference within the hot boundary layer. Since, for geometry reasons indicated above, the boundary layer temperature difference is (_f_−2) larger across the hot layer than across the cold one, the average radial profile is closer to the minimum profile. In Case 4, though the temperature difference across the cold boundary layer is larger than across the hot one, the effective cold sublayer that participates to convection is only a fraction of it (see Table 6). Assuming a simple force balance where the downwelling velocity is mostly due to the buoyancy of the unstable part of the cold boundary layer and the upwelling velocity to the buoyancy of the hot boundary layer, we expect that in both cases, the downwelling velocity is smaller than the upwelling velocity: this is what is observed. The asymmetry is even more pronounced when temperature-dependence of viscosity is included. This reflects the fact that, in a spherical shell, viscosity variations may be significant in the hot boundary layer while they are smaller than one order of magnitude in the Cartesian case. In addition, in the stagnant lid case, radial velocities vanish at the bottom of the lid.

Figure 10.

Time-averaged radial profiles corresponding to Cases 3 (top) and 4 (bottom) introduced above. Instantaneous data have been averaged for time periods corresponding to several convective overturns. Left: Radially averaged temperature as well as minimum and maximum temperatures at a given radius. Centre: Upwelling (circles) and downwelling (diamonds) velocities. Right: Radial advective q adv and diffusive q di heat fluxes and total power (circles) through a given radius-radial advection vanishes on the (r = r b) and (r = r t) boundaries while radial diffusion is small in the convective region. Dimensionless power P(r) (see text) should be constant throughout the layer for a steady-state flow. This value is plotted with symbols on the right-hand panel.

Table 6.

Global scalar values characteristic of Cases 3 and 4. These values are time-averages corresponding to the profiles displayed in Fig. 10. Nusselt numbers, boundary layers (and stagnant lid) thicknesses, temperature differences across the boundary layers (and the stagnant lid). Boundary layer thicknesses refer to the radius where radial advection and radial diffusion are equal. The stagnant lid thickness is defined with the radius of the intersection between the tangent to the advective profile (cf.Fig. 10) through its inflexion point and the zero-flux axis (cf.Davaille & Jaupart 1993).

Concerning radial heat fluxes, though the classical feature of two radii close to the bottom and top of the domain delimiting zones where radial diffusion exceeds radial advection (hot and cold boundary layers) from a zone where the contrary is true (well-mixed convecting core), it is interesting to note that the conservation of the power from one radius to the other implies that the sum of the advective and convective part of the radial flux diminishes as the spherical area increases (i.e. as the square of the radius) when no volumetric heating is taken into account (explaining the parabolic shape of the advective profile, Fig. 10). When a constant volumetric heating rate is introduced, the global power is proportional to P(r) = _r_2[q adv(r) +q di(r)]- (_r_3-r_3_b)h/3. The conservation of this quantity throughout the spherical shell gives a significant amount of information on the state of the flow in terms of distance from steady-state and possible resolution problems. In both Cases presented in Fig. 10, a statistical steady state is clearly reached [P(r) is nearly independent of _r_]. However, a departure from the nominal value can be observed at depths corresponding to the unstable boundary layers: this is more obvious for Case 4, especially for the hot bottom boundary layer (the deviation is larger than 20 per cent). This feature could be related to time period for the averaging that is too small. Here, it indicates that the convective regime would imply a more accurate grid mesh. It is interesting to note that though the temperature difference across the hot boundary layer is smaller for Case 4 than for Case 3, the resolution limit does not seem to be reached for Case 3. This is probably due to the fact that the gradient is about 20 per cent larger in Case 4 (see the bottom Nusselt number values in Table 6).

5.3 Coupled computation: internal dynamics, topography and geoid

Dynamic topography and geoid have been computed for the four thermal convection results presented above (Table 4). We have used the following values for the dimensionless parameters: E = 5 × 10−2 and G*= 2 × 10−2. It should be noted that this choice, while it is based on plausible orders of magnitude for a generic planetary mantle (cf.Table 7), remains nevertheless rather arbitrary. Note that both the amplitude of the presented geoid and the relative influence of (i) topographies at the interfaces and (ii) density anomalies throughout the spherical shell, are affected by this choice. Finally, it should be noted that since the dimensionless value of geoid N implies the use of dimensionless parameter G*, an alternative variable may be the dimensionless potential at the surface of the shell U(r t) that does not require the introduction of further dimensionless quantities. Indeed, N is simply proportional to U(r t). However, if dynamic topographies are to be computed from numerical simulations of thermal convection, E needs to be introduced.

Table 7.

Acceptable range for dimensionless parameters E and G* (as well as internal Rayleigh number Ra i and dimensionless radius of curvature f of the shell as an indication). The two end-members values for E (from Mercury, or a small silicate mantle, to Venus or a large silicate mantle) are not associated to the same objects than for G* (from Europa, or a thin icy crust, to Venus, or a large silicate mantle). Note that most of the crude estimates within this table can be debated (especially the values of d, Δ_T_ and μ0). Gravitational constant G = 6.67 × 10−11 m3 s−2 kg−1.

Figs 11–14 display the results of the computation for Cases 1–4, respectively. According to Section 4.2, the data needed to produce these maps are the radial stress σ_rr_ on the interfaces and the 3-D temperature field T. These gridded data sets (with variable areas) are first transposed into the spherical harmonics basis (coefficients are computed to a maximum degree ℓmax of 30). Dimensionless topography and geoid are then computed according to the procedure described in 4.2 and finally plotted on a latitude-longitude grid. Thus, degrees shorter than ℓmax = 30 do not appear on topography and geoid fields in Figs 11–14.

Figure 11.

Surface dynamic topography (T) and geoid (N) obtained for stationary temperature and flow fields corresponding to Case 1.

Figure 12.

Surface dynamic topography (T) and geoid (N) obtained for stationary temperature and flow fields corresponding to Case 2.

Figure 13.

Surface dynamic topography (T) and geoid (N) obtained for instananeous temperature and flow fields corresponding to Case 3.

Figure 14.

Surface dynamic topography (T) and geoid (N) obtained for instananeous temperature and flow fields corresponding to Case 4.

A first result is that both the tetrahedral and cubic Cases 1 and 2 (Figs 11 and 12, respectively) lead to topography and geoid variations that still exhibit, among others, the strong spherical harmonic coefficients inherited from the initial perturbation: (ℓ = 3, m = 2) for tetrahedral Case 1, (ℓ = 4, m = 0) and (ℓ = 4, m = 4) for cubic Case 2. The amplitudes are rather large (on the order of a few percents of the planet's radius for surface dynamic topography) and indicate that these cases are probably not relevant for a comparison in the planetary context.

Time-dependent Cases 3 and 4 lead to shorter wavelengths variations of both surface dynamic topography and geoid. Again, these features directly reflect the wavelength observed for thermal convection in the spherical shell: the power spectrum associated to the results in Case 3 peaks for degree (ℓ = 8) while the maximum is obtained for higher degrees in Case 4 (ℓ = 12 and larger, the uncertainty associated to degrees higher than 15 in the spectral evaluation of geoid and topography being non-negligible). This result is mostly induced by value of the internal Rayleigh number Ra i (i.e. obtained for a viscosity value characteristic of the convective sublayer): Ra i = Ra = 1.4 × 106 for isoviscous Case 3 while Ra _i_∼ 8 × 107 for Case 4. The major effect of the viscous conductive lid in Case 4, compared to the isoviscous Case 3, is the increase of the amplitude of the topography while the geoid has the same order of magnitude in both cases. It is important to note that, as demonstrated by Solomatov (2004), an Arrhenius-type viscosity law would provide a larger deviatoric stress field near the surface of the stagnant lid. Its influence on dynamic topography is probably only weak since topography is mainly maintained by variations in pressure, which mostly reflects the density distribution. The effect on geoid might be worth examinating.

6 Conclusion

The variations of geoid and surface dynamic topography contain important information about internal dynamic processes occurring within silicate mantles of terrestrial planets or surface icy layers of the large satellites of giant planets. Though, ultimately, a precise estimate of the internal structure will also require ground-based geophysical experiments (see, e.g. Verhoeven et al. 2005), observations such as the geoid and topography can be obtained from experiments in an orbiter mission.

In this study, we propose a dimensionless framework for this coupled problem: internal dynamics and the resulting geoid and dynamic topography. The characteristic scales are inherited from infinite Prandtl number thermal convection, with the four dimensionless parameters Ra, f, γ and h. Two additional dimensionless parameters are introduced in order to produce dimensionless values for geoid and topography: E and G* (Table 7 indicates plausible ranges for this two terms, much smaller than the one associated to the Rayleigh number). As well as the velocity scale based on thermal diffusion κ/d is not significant for velocities associated to natural thermal convection, the length scale based on the shell thickness is much larger than realistic values for N, t s. Finally, a good alternative to dimensionless geoid may be the surface value of the dimensionless gravitational potential, U(r t).

In order to model internal dynamics, we propose an efficient numerical method based on the ‘cubed-sphere’ grid mesh where the spherical shell is decomposed on six non-conformal blocks. We present a high-resolution scheme (with a slope-limiter non-linear method) to treat the advection term in the energy equation as well as a benchmark of various schemes. The flow solver is already described for one of the blocks in (Choblet 2005). It is tested extensively against a very precise spectral method for flows with a radial-dependence of viscosity. The response functions for geoid and dynamic topographies indicate a very good agreement up to degree ℓ = 15 for a 6 × 64 × 64 × 64 grid mesh in the flow solver (a correct description of higher degrees requires finer meshes). We believe that such an agreement is needed to apply thermal convection results to the computation of geoid and topography for planetary applications. Since various new methods for thermal convection in a spherical shell now exist, a precise benchmark should include such a calculation of geoid and topography.

However, the physical model proposed here is still oversimplified. For example, the adjunction of an increase of viscosity with depth (due to pressure dependence and/or to a viscous stratification within the mantle) has been shown to reduce significantly the characteristic degree of spherical convection (Bunge 1997; Roberts 2006; Yoshida & Kageyama 2006). Similarly, compositional effects on density and viscosity would lead to different structures and thus different results for geoid and topography. More, the rheological description of the stagnant lid is clearly the weakest point of many numerical codes used to model the flow in planetary interiors and ŒDIPUS is not an exception. The effect of a more realistic Arrhenius type viscosity law has already been discussed. The simple introduction of an elasticity component (as in Zhong 2002, for example) would also modify significantly the results: in our viscous models, a surface dynamic topography is obtained even for a very viscous lid while an elastic behaviour would largely reduce the viscous component. Indeed, the observed geoid and topography are superpositions of various contributions among which the dynamic one may be the least important (e.g. Le Stunff & Ricard 1995). A large portion of the topographic load may be maintained by isostasy and/or elastic flexure, depending on the wavelength of the topographic features, the thickness of the lithosphere and the radius of the body (e.g. Turcotte et al. 1981; Pauer et al. 2006). A careful analysis of the data is thus necessary for a correct interpretation of the geoid and topography in terms of viscous flow models.

For these reasons, a direct comparison of the topography and geoid obtained from a convection simulation with planetary observables is likely to be misleading in many cases. An intermediate approach may include the computation of dimensionless geoid and topography for a large range of model parameters and for various physical models of the internal dynamics. Next steps can then be the a posteriori correction of the predicted geoids and topographies for lithospheric effects (e.g. by including an elastic layer on the top or considering a membrane correction for a stiff uppermost layer) and their comparison with observed data, at least partly reduced to isostatic effects or at wavelengths where the dynamic contribution is presumably dominant.

Acknowledgments

We thank S. Zhong and an anonymous reviewer for their comments. The work benefited from Charles University grant GAUK 280/2006/B-GEO/MFF as well as the research project MSM 0021620800 of the Ministry of Education in the Czech Republic and from the Programme National de Planétologie as well as ANR ETHER in France.

References

,

1985

.

Three-dimensional treatment of convective flow in the earth s mantle

,

J. Stat. Phys.

,

39

,

501

–

511

.

,

2003

.

The generation of plate tectonics from mantle convection

,

Earth planet. Sci. Lett.

,

205

,

107

–

121

.

,

1989

.

Three-dimensional spherical models of convection in the Earth-s mantle

,

Science

,

244

,

950

–

955

.

,

1997

.

A sensitivity study of three-dimensional spherical mantle convection at 108 Rayleigh number: Effects of depth-dependent viscosity, heating mode, and an endothermic phase change

,

J. Geophys. Res.

,

102

,

11991

–

12008

.

,

1993

.

Gravity Field and Dynamics of the Earth

,

Springer-Verlag, Berlin

.

,

1989

.

Spherical tensor approach to the solution of the mantle stress problem

,

Studia Geoph. Geod.

,

33

,

117

–

197

.

,

2003

.

Effect of lateral viscosity variations in the top 300 km on the geoid and dynamic topography

,

Geophys. J. Int.

,

152

,

566

–

580

.

,

2005

.

Modelling thermal convection with large viscosity gradients in one block of the cubed sphere

,

J. Comput. Phys.

,

205

,

269

–

291

.

,

2001

.

Mantle upwelling and melting beneath slow spreading centers: effects of realistic rheology and melt productivity

,

Earth Planet. Sci. Lett.

,

184

,

589

–

604

.

,

1993

.

Transient high-rayleigh-number thermal convection with large viscosity variations

,

J. Fluid Mech.

,

253

,

141

–

166

.

,

2004

.

Towards a lower mantle reference temperature and composition

,

Earth planet. Sci. Lett.

,

222

,

161

–

175

.

,

1969

.

Diffusion and Heat Transfer in Chemical Kinetics.

Plenum Press

,

New York.

,

1988

.

Numerical simulations of mantle convection: time-dependent, three-dimensional, compressible, spherical shell.

,

Geophys. Astrophys. Fluid Dyn.

,

43

,

223

–

264

.

,

1959

.

A finite difference scheme for numerical computation of discontinuous solutions of equations of fluid dynamics

,

Mat. Sb.

,

47

,

271

–

306

.

,

1989

.

Constraints on the structure of mantle convection using seismic observations, flow models, and the geoid

, in

Mantle Convection: Plate Tectonics and Global Dynamics

,

Vol. 78

, pp.

675

—

763

, ed.,

Gordon and Breach

,

New York.

,

1996

.

A one-plume model of martian mantle convection

,

Nature

,

380

,

507

–

509

.

,

2005

.

A finite-volume solution method for thermal convection and dynamo problems in spherical shells

,

Geophys. J. Int.

,

161

,

522

–

532

.

,

2003

.

Three-dimensional spherical shell convection at infinite prandtl number using the ‘cubed sphere’ method.

, in

Computational Fluid and Solid Mechanics

, p.

931

, ed.,

Elsevier

.

,

2005

.

A model for the thermal and chemical evolution of the Moon-s interior: implications for the onset of mare volcanism

,

Earth Planet. Sci. Lett.

,

134

,

501

–

514

.

,

1995

.

A preliminary study of the effects of some flow parameters in the generation of poloidal and toroidal energies within a 3d spherical thermal convection system with variable viscosity

,

Pure Appl. Geophys.

,

145

,

487

–

503

.

,

2004

.

Constraining large-scale mantle heterogeneity using mantle and inner-core sensitive normal modes

,

Phys. Earth Planet. Inter.

,

146

,

113

–

124

.

,

1985

.

Spherical Harmonics and Tensors for Classical Field Theory

,

Research Studies Press Ltd.

,

Letchworth, England.

,

2004

.

Yin-yang grid: an overset grid in spherical geometry

,

Geochem. Geophys. Geosyst.

,

5

,

Q09005

,

.

,

1995

.

Models of mantle viscosity

, in

Mineral Physics & Crystallography, A Handbook of Physical Constants

, pp.

227

–

236

, ed.,

AGU

,

Washington, DC.

,

1988

.

Physical Geodesy

,

Claredon Press

,

Oxford.

,

1995

.

Topography and geoid due to lithospheric mass anomalies

,

Geophys. J. Int.

,

122

,

982

–

990

.

,

2004

.

Thermochemical structures within a spherical mantle: Superplumes or piles?

,

J. Geophys. Res.

,

109

,

7402

–+.

,

1996

.

Constraints on the lateral srength of slabs from three-dimensionnal dynamic flow models

,

Earth Planet. Sci. Lett.

,

138

,

15

–

28

.

,

2006

.

Modeling the dynamic component of the geoid and topography of Venus

,

J. Geophys. Res.

,

111

,

11012

–+.

,

1996

.

Steady tetrahedral and cubic patterns of spherical shell convection with temperature-dependent viscosity

,

J. Geophys. Res.

,

101

,

25 473

–

25 484

.

,

1999

.

Non-Newtonian stagnant lid convection and magmatic resurfacing on Venus

,

Icarus

,

139

,

67

–

80

.

,

1984

.

Geoid heights and lithospheric stresses for a dynamic earth

,

Ann. Geophys.

,

2

,

267

–

285

.

,

1984

.

Geoid anomalies in a dynamic Earth

,

J. Geophys. Res.

,

89

,

5987

–

6002

.

,

2006

.

Degree-1 convection in the Martian mantle and the origin of the hemispheric dichotomy

,

J. Geophys. Res.

,

111

6013

–+. Stemmer, K., Harder, H., & Hansen, U., 2006. A new method to simulate convection with strongly temperature- and pressure-dependent viscosity in a spherical shell: Applications to the Earth-s mantle, Phys. Earth Planet. Inter., 157, 223—249.

,

1996

.

The ‘cubed sphere’: a new method for the solution of partial differential equations in spherical geometry

,

J. Comput. Phys.

,

124

,

93

–

114

.

,

2001

.

Mantle Convection in the Earth and Planets.

,

Cambridge University Press

,

Cambridge.

,

1984

.

A fully multidimensional positive definite advection transport algorithm with small implicit diffusion

,

J. Comput. Phys.

,

54

,

325

–

362

.

,

2004

.

Initiation of subduction by small-scale convection

,

J. Geophys. Res.

,

109

,

1412

–+.

,

2006

.

A new method to simulate convection with strongly temperature- and pressure-dependent viscosity in a spherical shell: Applications to the Earth's mantle

,

Phys. Earth Planet. Inter.

,

157

,

223

–

249

.

,

1994

.

Three-dimensional models of mantle convection: influence of phase transitions and temperature-dependent viscosity

,

PhD thesis

,

California Institute of Technology

,

Pasadena.

,

2000

.

Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations 1. Pseudo-plastic yielding

,

Geochem. Geophys. Geosyst.

,

1

,

.

,

2000

.

Self-consistent generation of tectonic plates in time-dependent, three-dimensional mantle convection simulations 2. strain-weakening and asthenosphere

,

Geochem. Geophys. Geosyst.

,

1

,

.

,

2003

.

Testing the tracer ratio method for modeling active compositional fields in mantle convection simulations

,

Geochem. Geophys. Geosyst.

,

4

,

.

,

2006

.

Episodic outgassing as the origin of atmospheric methane on Titan

,

Nature

,

440

,

61

–

64

.

,

1981

.

Role of membrane stresses in the support of planetary topography

,

J. Geophys. Res.

,

86

,

3951

–

3959

.

,

2004

.

Geophysical evidence for chemical variations in the Australian Continental Mantle

,

Geophys. Res. Lett.

,

31

,

17 607

.

,

1979

.

Towards the ultimate conservative difference scheme V. A second order sequel to Godunov-s method

,

J. Comput. Phys.

,

32

,

101

.

,

1988

.

Quantum Theory of Angular Momentum

,

World Scientific

.

et al. ,

2005

.

Interior structure of terrestrial planets: modeling Mars' mantle and its electromagnetic, geodetic, and seismic properties

,

J. Geophys. Res.

,

110

,

4009

.

,

2004

.

Application of the yin-yang grid to a thermal convection of a boussinesq fluid with infinite prandtl number in a three-dimensional spherical shell

,

Geophys. Res. Lett.

,

31

,

L12 609

.

,

2006

.

Low-degree mantle convection with strongly temperature- and depth-dependent viscosity in a three-dimensional spherical shell

,

J. Geophys. Res.

,

111

,

3412

.

,

1995

.

The influences of lower mantle viscosity stratification on 3D spherical-shell mantle convection

,

Earth Planet. Sci. Lett.

,

132

,

157

–

166

.

,

2002

.

Effects of lithosphere on the long-wavelength gravity anomalies and their implications for the formation of the Tharsis rise on Mars

,

J. Geophys. Res.

,

107

,

8

–

1

.

,

2000

.

Role of temperature-dependent viscosity and surface plates in spherical shell models of mantle convection

,

J. Geophys. Res.

,

105

,

11 063

–

11 082

.

Appendix

Appendix A: Description of the Advection Schemes

The motivation to introduce high resolution methods arise from the problem stated in ‘Godunov's theorem’: There are no monotone, linear schemes for the linear advection equation of second or higher order of accuracy. When only linear schemes are considered, one is forced to choose between first order schemes that produce numerical diffusion or second order schemes that are not monotonic (i.e. produce spurious, non-physical oscillations). High resolution methods thus introduce a non-linear treatment of the advective flux term.

Neglecting for the moment the diffusive terms of right-hand side, the advective part of eq. (16) can be written as follows

(A1)

where is the Jacobian of the coordinate transformation from the Cartesian basis to the natural basis associated to the cubed sphere coordinates and contravariant velocity components refer to this natural basis (while V_• refer to the physical basis, see Choblet (2005) for further details). Integrating eq. (A1) over the volume of cell C and during time step Δ_t, one obtains the conservation form

(A2)

with

(A3)

and

(A4)

where index I refer to the direction (r, ξ or η) so that S I is the surface of the cell C corresponding to I = constant (unless explicitly written, the coordinates are r, ξ, η, t).

Let us now consider the discrete problem. The staggered grid mesh approach is used so that velocities are localized on cell boundaries (corresponding to r_−_i, r+i, ξ−j, ξ+j, η−k, η+k) and scalar quantities in the centre (corresponding to r i, ξ_j_, η_k_) (see Choblet 2005, for further details). The Godunov approach (Godunov 1959) then consists in the introduction of a function constant in the volume of cell C but varying through time: γ_n_ ijk (t). We thus obtain for the discrete flux on the cell boundaries

(A5)

(here, and in the following lines, unless explicitly written, the indexes are i, j, k, n). The first order accurate upwind scheme corresponds then to the choice of a constant value for γ corresponding either to the ‘left’ or the ‘right’ value T, depending on the sign of , for example,

(A6)

with Δ_l_ξ, Δ_l_η the dimensions of the cell. A better (second order) precision can be obtained by adopting a (tri)linear function for γ, introducing slopes σ_r_, σξ, ση, constant in the time interval [t, t+Δ_t_]

(A7)

Note that the case σ• = 0 provides again the Godunov-upwind method from a different point of view. Classical choices for the slope including a centred (e.g. σ_r_ = (T i+1−T i_−1)/2Δ_r, Fromm), upwind (e.g. σ_r_ = (T i_−_T i_−1)/Δ_r, Beam-Warming) or downwind (e.g. σ_r_ = (T i+1−T i)/Δ_r_, Lax-Wendroff) however may induce oscillations in regions of large T gradients: in this context, a way to avoid these oscillations is to limit the slope's values (slope- or flux-limiter methods) thus inducing the non-linear aspect of the method. Van Leer's solution (van Leer 1979) is then to consider the monotonized central-difference limiter (or MC-limiter) defined as follows, for example in the _r_-direction:

(A8)

with the minmod function:

(A9)

This slope is termed slope 1 in our comparative study. In smooth regions, this reduces to the centred (Fromm) slope and gives a reasonably sharp resolution in large gradients regions by comparing it with twice the one-sided slopes. This is the reason why the use of non-linear limiters does not only make sense in the context of compressible flows where it was first introduced but also for incompressible flows where large gradients occur. Another solution is to compute the slope as a geometrical average:

(A10)

The Multidimensional Positive Definite Advection Transport Algorithm (MPDATA) scheme corrects the classical upwind procedure by reapplying it with an antidiffusive velocity field. The latter is derived from the truncation error analysis of the upwind scheme assuming a uniform flow. It is clear that this procedure can then be reapplied so that a good precision can be reached, iteratively (for a detailed description, see Smolarkiewicz 1984). This method has been used successfully to model mantle dynamics in the Cartesian context. However, the corrective steps are quite time-consuming in the cubed-sphere coordinate system. Since a rigorous truncation error analysis in a multidimensional complex flow context is not possible, we do not study the precise effect of both high-resolution methods introduced here (van Leer slope limiter, MPDATA) on the physics of the flow in terms of oscillations. The single effect on numerical diffusion is treated for a simple solid rotation flow.

Appendix B: Communication Between Blocks

Within each of the six blocks, in the framework of parallel computation, communication between processors is identical to a classical domain decomposition in a Cartesian parallelepipedic volume (cf.Choblet 2005). The main object here is thus to describe communication between the six blocks. As can be seen on Fig. 3, the ‘new’ orientation of the coordinate system induces only two configurations: either the coordinates ‘match’ on two adjacent blocks (one of the coordinate is the same modulo π/2, e.g. ξ for blocks 0 and 1 or η for blocks 1 and 4) or they are inversed (ξ becomes η, e.g. from 1 to 2, or vice versa, e.g. from 3 to 0).

Three kinds of variables may then be considered, due to the discretization:

• (i)
  scalar variables and velocity component V r, though they are ‘located’ at different depths in the staggered grid mesh approach, are however similar in terms of angular coordinates ξ and η. The buffer planes outside the block, at −π/4 -Δ/2 and π/4 +Δ/2, are filled using the values from inside the adjacent blocks (e.g. located, respectively at −π/4 +Δ/2 and π/4 −Δ/2 in the case of a simple configuration), a simple 1-D interpolation being necessary since only one of the coordinates is different between the two blocks (e.g. η0 and η1 are different).
• (ii)
  when dealing with azimuthal vector components, a 2-D transformation is needed due to non-orthogonality: if we consider the example of blocks 0 and 1, corresponding to a simple configuration,
  
  (B1)
  (these relations are similar in the case of an inverse configuration): in the case of normal velocity component (here _V_ξ), the discrete buffer planes located at -π/4 −Δ and π/4 are filled so that only the discrete values between -π/4 and π/4 −Δ are computed; again a simple 1-D interpolation is performed, with an additional scaling factor corresponding to coefficient in eq. (B1).
• (iii)
  in the case of a velocity component that is tangent to the boundary (_V_η for the boundary between 0 and 1), as in the scalar case, the buffer planes outside the block, at -π/4 −Δ/2 and π/4 +Δ/2, are filled using the values of the corresponding tangent velocity from within the adjacent blocks (at -π/4 +Δ/2 and π/4 −Δ/2) and of the normal velocity at -π/4 and π/4 −Δ.

It should be noted that in spite of non-orthogonality, these interpolation and transformation procedures are very simple compared to the one induced by most multiblocks algorithms.

© 2007 The Author Journal compilation © 2007 RAS

Citations

Views

Altmetric

Metrics

Total Views 682

469 Pageviews

213 PDF Downloads

Since 12/1/2016

Citations

42 Web of Science

×

Email alerts

Citing articles via

• Latest

• Most Read

• Most Cited

More from Oxford Academic