# Constraints on core–mantle boundary topography from models of thermal and thermochemical convection

Constraints on core–mantle boundary topography from models of thermal and thermochemical... Abstract Mantle flow induces dynamic topography at the core–mantle boundary (CMB), with distribution and amplitude that depend on details of the flow. To assess whether observations of CMB topography can give constraints on deep mantle structure, we determine CMB dynamic topography associated with different models of mantle convection, including thermochemical and purely thermal models. We investigate the influence of key controlling parameters, specifically the thermal viscosity ratio (ΔηT) and, for thermochemical models, the density contrast (ΔρC) and viscosity ratio (ΔηC) between primordial and regular materials. In purely thermal models, plume clusters induce positive topography with an amplitude that decreases with increasing ΔηT. In thermochemical models with moderate density contrasts, around 100–200 kg m−3, reservoirs of dense material induce depressions in CMB topography, surrounded by a ridge of positive topography. The average depression depth and ridge height increase with increasing ΔρC and ΔηC, but decrease with increasing ΔηT. We find that for purely thermal models or thermochemical models with ΔρC ∼ 90 kg m−3 and less, the long-wavelength (spherical harmonic degrees up to l = 4) dynamic topography and shear wave velocity anomalies predicted by thermochemical distributions anticorrelate. By contrast, for models with ΔρC ≥ 100 kg m−3 and ΔηC > 1, long-wavelength dynamic topography and shear wave velocity anomalies correlate well. This potentially provides a test to infer the nature, that is, either purely or mostly thermal (ΔρC ≤ 100 kg m−3 m−3) or strongly thermochemical (ΔρC ≥ 100 kg m−3), of the low shear wave velocity provinces observed by global tomographic images. The presence of post-perovskite, provided that its viscosity is similar to that of bridgmanite, does not alter these conclusions. Composition and structure of the mantle, Mantle processes, Structure of the Earth, Dynamics: convection currents, and mantle plumes 1 INTRODUCTION The seismic structure of the lowermost mantle is dominated by two large low shear wave velocity provinces (LLSVPs) below Africa and the Pacific, where shear wave speed decreases by a few per cent compared to its horizontal average. LLSVPs are observed by all global seismic tomography models, with some differences in the detailed structure and amplitude of the velocity drop (for a comparison of models published between 2005 and 2010, see Ritsema et al.2011). Their presence is further supported by a variety of seismic data and methods (Lekić et al.2012), suggesting that they are not simple artefacts of tomographic models. In particular, because they are seen by normal mode data (Ishii & Tromp 1999; Trampert et al.2004; Mosca et al.2012), they are unlikely the result of wave-propagation complexities. While the presence of LLSVPs is generally accepted, their nature is still a matter of debate (for recent discussions on this problem, see Davies et al.2015; Deschamps et al.2015; Garnero et al.2016). LLSVPs were first associated with superplumes or clusters of smaller plumes, that is, buoyant structures triggered by purely thermal (hot) anomalies. Several hints, however, including the anticorrelation between shear wave and bulk-sound velocity anomalies (e.g. Ishii & Tromp 1999; Masters et al.2000; Trampert et al.2004), and the de-correlation between density and seismic velocity anomalies (Trampert et al.2004), indicated that LLSVPs cannot be explained by purely thermal anomalies, but partially originate from compositional anomalies (Trampert et al.2004; Mosca et al.2012). Alternatively, it has been proposed that the anticorrelation between shear wave and bulk-sound speeds could be mostly explained by complex wave-propagation effects (Schuberth et al.2012), or by the presence of post-perovskite (pPv; Davies et al.2012). Thus, a clear diagnosis on the nature of LLSVPs, purely thermal or thermochemical, cannot be based on current seismic velocity data only, but requires additional observational constraints independent of seismic velocities. Two promising observables are long-period electromagnetic C-responses, which are related to the distribution of electrical conductivity in the deep mantle and are different for purely thermal and thermochemical structures (Deschamps & Khan 2016), and seismic attenuation, which is mostly sensitive to temperature and can be inferred from waveform data (e.g. Konishi et al.2017). Another potential constraint is the core–mantle boundary (CMB) topography. Using Very Long Baseline Interferometry (VLBI) data, Gwinn et al. (1986) showed that the CMB dynamical flattening (strictly speaking, the dynamical flattening of the outer core), that is, the spherical harmonic degree and order l = 2 and m = 0 of CMB topography, is larger than its theoretical value for a hydrostatic Earth model by 5 per cent. If more complexities are taken into account, the excess CMB dynamical flattening reduces to 3.8 per cent (Mathews et al.2002). Seismic normal modes data further indicate that the peak-to-peak amplitude of the spherical harmonic degree l = 2 in CMB topography should not exceed 5.0 km (Koelemeijer et al.2012). Analysis of seismic phases reflecting on the CMB underside (PKKP) and upperside (PcP) suggests that 95 per cent of the CMB topography is distributed between −4.0 and 4.0 km (Garcia & Souriau 2000). Several studies, using different CMB transmitted or reflected seismic phases, have also attempted to draw global maps, but generally disagree on both the pattern and amplitude of CMB topography (Morelli & Dziewonski 1987; Doornbos & Hilton 1989; Sze & van der Hilst 2003; Tanaka 2010). For spherical harmonic degrees up to l = 4, the peak-to-peak amplitude ranges from 3.5 km (Koper et al. 2003) to 12 km (Morelli & Dziewonski 1987). Insights on CMB topography may further be obtained from geodynamic models. In a fluid animated by convection, flow generates normal stress at its interfaces, which in turn induces dynamic topography (e.g. Ricard et al.1984; Hager & Richards 1989). Plumes stretch the base and push the surface of the fluid upwards (Fig. 1a), inducing positive topography on these boundaries. In contrast, downwellings pull the surface and push the base of the system downwards, leading to negative topography (Fig. 1b). Dynamic topography is sensitive to mantle viscosity and its lateral variations, and to the density jump at the boundaries. In the Earth's mantle, the detailed CMB dynamic topography thus strongly depends on the mode of mantle dynamics, in particular on lateral variations in density due to large-scale chemical heterogeneities. The 2-D-Cartesian simulations of Hansen & Yuen (1989) indicated that piles of dense material at the bottom of the system strongly depress the CMB, inducing negative topography (Figs 1c and d). Lassak et al. (2007) pointed out that in models combining large thermal viscosity contrast and piles with small chemical density contrast, the negative CMB topography induced by slabs is much larger than that caused by the pile. Lassak et al. (2010) further showed that imposing a large thermal viscosity contrast flattens out the topography induced by thermochemical piles, because these regions are hotter, and thus less viscous, than surroundings. Using models of convection in which density anomalies are scaled from seismic tomography, Yoshida (2008) showed that if LLSVPs are purely thermal structures, they should result in positive topography with amplitude around 5 km, while, on the contrary, if they are thermochemical piles, they would induce local depressions. These models, however, are limited to moderate chemical density contrasts, around 60 kg m−3 and lower, and assume that the viscosity is independent of the type of material, dense or regular. Figure 1. View largeDownload slide Influence of mantle dynamics on core–mantle boundary (CMB) topography. Thin black arrows denote the stress applied on the CMB. (a) Thermal plumes pull the CMB upwards, inducing positive topography. (b) Slabs stretch the CMB downwards, inducing negative topography. (c,d) Thermochemical piles depress the CMB, inducing negative topography with amplitude that depends on the chemical density contrast between the material of these piles and the regular mantle. Some elements of this sketch are taken from Garnero et al. (2016). Figure 1. View largeDownload slide Influence of mantle dynamics on core–mantle boundary (CMB) topography. Thin black arrows denote the stress applied on the CMB. (a) Thermal plumes pull the CMB upwards, inducing positive topography. (b) Slabs stretch the CMB downwards, inducing negative topography. (c,d) Thermochemical piles depress the CMB, inducing negative topography with amplitude that depends on the chemical density contrast between the material of these piles and the regular mantle. Some elements of this sketch are taken from Garnero et al. (2016). In this study, we performed models of convection in 3-D-spherical geometry, including models with chemical density contrast larger than 100 kg m−3, which may be more appropriate for LLSVPs, models with an intrinsic viscosity ratio between primordial and regular materials, whose effects on CMB dynamic topography have not yet been investigated, and models accounting for the pPv phase transition. Results from these models provide a criteria to discriminate between possible models of LLSVPs, purely thermal or slightly thermochemical on one hand, or strongly thermochemical, that is, with chemical density contrast around 100 kg m−3 or more, on the other hand. 2 MANTLE CONVECTION MODELLING We performed numerical simulations of purely thermal and thermochemical convection in 3-D-spherical geometry using StagYY (Tackley 1998, 2008), a software that solves the conservation equations of energy, momentum, and mass for a compressible, anelastic fluid with infinite Prandtl number. Thermochemical models distinguish two types of material, regular (pyrolitic) mantle and primordial material, and further solve the conservation equation of composition. Primordial material accounts for chemical heterogeneities that may be present at the bottom of the mantle as a result of early differentiation. An additional source of chemical heterogeneity, which we do not take into account in our models, is the recycling of subducted oceanic crust (MORB). Recent models of convection (Li et al.2014a) indicate that subducted MORB may be partially incorporated into pre-existing reservoirs of primordial material, and partially re-entrained to the surface by plumes, in agreement with the Basal Melange (BAM) scenario (Tackley 2012). Because subducted MORB is denser than pyrolitic mantle, their effect on CMB topography is likely similar to that of primordial material, and neglecting its contribution may not alter our conclusions. For details on the numerical methods employed to solve the conservation equations, we refer to Tackley (2008). The setup used in this paper is essentially that used in Deschamps et al. (2015) except for the viscosity law, which now includes yield stress, and for the explored range of control parameters. Details of this setup can be found in previous studies (Deschamps & Tackley 2009), but for clarity we discuss them again below. Spherical geometry is modelled by a set of Yin and Yang grids (Kageyama & Sato 2004) with resolution equivalent to a spherical grid with 256 points in latitude and 512 points in longitude at each depth. The vertical resolution is fixed to 128 or 64 nodes for models with and without pPv, respectively. Grid refinement is imposed at the top and at the bottom of the shell, allowing a good resolution of reservoirs of primordial material and thermal boundary layers. The ratio between the radius of the core and the total radius is set to its Earth value, that is, f = 0.55. The bottom and surface boundaries are free slip. The system is heated both from the bottom and from within, and the heating rate of primordial material is increased by a factor 10 compared to that of regular material. The total mantle heating rate is equivalent to a surface heat flux of 10 mW m−2, which is about half the value expected for Earth's mantle (Jaupart et al.2015), and leads to a ratio of internal to basal heating close to 0.3. Compressibility generates additional sinks and sources of heat, that are controlled by the dissipation number, Di. Because it depends on thermal expansion, Di varies with depth. In all our models the surface value of the dissipation number is fixed to DiS = 1.2, and the depth variations of thermal expansion (Fig. 2) imply that its volume average is equal to 0.43. Figure 2. View largeDownload slide Thermodynamical reference state as a function of depth. Profiles shown are for temperature (Tref), density (ρref), thermal expansion (αref) and Grüneisen parameter (γref). The discontinuity in density at non-dimensional depth of 0.228 corresponds to the limit between upper and lower mantle. Figure 2. View largeDownload slide Thermodynamical reference state as a function of depth. Profiles shown are for temperature (Tref), density (ρref), thermal expansion (αref) and Grüneisen parameter (γref). The discontinuity in density at non-dimensional depth of 0.228 corresponds to the limit between upper and lower mantle. The conservation equations involve the definition of a reference state, which here consists of a set of vertical profiles for density, temperature, and thermal expansion calculated using thermodynamic relationships for the Earth's mantle (Tackley 1998). Fig. 2 plots the variations of the reference state with depth. Note that thermal expansion decreases by a factor 5 from surface to bottom, and that the Grüneisen parameter, γref, varies in such a way that its product with reference density is constant with depth. Temperature is scaled with respect to the super-adiabatic temperature difference across the system, ΔTS, and density and thermal expansion are scaled with respect to their surface values. Table 1 lists the thermodynamic properties used to build the reference state, together with other scalings and properties. Table 1. Parameters and scalings of numerical experiments. Parameter  Symbol  Value  Units  Non-dimensional  Non-dimensional parameters          Reference Rayleigh number  RaS      3.0 × 108  Surface dissipation number  DiS      1.2  Volume average dissipative number  Di      0.43  Total internal heating  HC  10  mW m−2  4.8  Compositional heating ratio  ΔHC  10                Compositional parameters (thermochemical cases only)          Buoyancy ratio  Bz      0.15–0.70  Volume fraction of dense material (%)  Xprim      2.5–7.0            Physical & thermodynamical parameters          Acceleration of gravity  g  9.81  m s−2  1.0  Mantle thickness  D  2891  km  1.0  Reference adiabat  Tas  1600  K  0.64  Super-adiabatic temperature difference  ΔTS  2500  K  1.0  Surface density  ρS  3300  kg m−3  1.0  Surface thermal expansion  αS  5.0 × 10−5  K−1  1.0  Surface thermal diffusivity  κS  6.24 × 10−7  m2 s−1  1.0  Heat capacity  CP  1200  J kg−1 K−1  1.0  Surface conductivity  kS  3.0  W m−1 K−1  1.0  Surface Grüneisen parameter  γS  1.091      Density jump at z = 660 km  Δρ660  400  kg m−3  0.1212  Clapeyron slope at z = 660 km  Γ660  −2.5  MPa K−1  −0.0668  Post-perovskite density jump  ΔρpPv  62  kg m−3  0.0188  Clapeyron slope of post-perovskite  ΓpPv  13  MPa K−1  0.3474  CMB temperature  TCMB  3750  K  1.5  Density jump at CMB  ΔρCMB  5280  kg m−3  1.6            Viscosity law          Reference viscosity  η0  1.6 × 1021  Pa s  1.0  Viscosity ratio at z = 660 km  Δη660  30      Logarithmic thermal viscosity ratio  Ea  13.816–23.026      Logarithmic vertical viscosity ratio  Va  2.303      Compositional viscosity ratio  ΔηC  1–102      Surface yield stress  σ0  290  MPa  7.5 × 106  Yield stress gradient  $${\dot{\sigma }_z}$$  0.01  Pa/Pa  0.01  Parameter  Symbol  Value  Units  Non-dimensional  Non-dimensional parameters          Reference Rayleigh number  RaS      3.0 × 108  Surface dissipation number  DiS      1.2  Volume average dissipative number  Di      0.43  Total internal heating  HC  10  mW m−2  4.8  Compositional heating ratio  ΔHC  10                Compositional parameters (thermochemical cases only)          Buoyancy ratio  Bz      0.15–0.70  Volume fraction of dense material (%)  Xprim      2.5–7.0            Physical & thermodynamical parameters          Acceleration of gravity  g  9.81  m s−2  1.0  Mantle thickness  D  2891  km  1.0  Reference adiabat  Tas  1600  K  0.64  Super-adiabatic temperature difference  ΔTS  2500  K  1.0  Surface density  ρS  3300  kg m−3  1.0  Surface thermal expansion  αS  5.0 × 10−5  K−1  1.0  Surface thermal diffusivity  κS  6.24 × 10−7  m2 s−1  1.0  Heat capacity  CP  1200  J kg−1 K−1  1.0  Surface conductivity  kS  3.0  W m−1 K−1  1.0  Surface Grüneisen parameter  γS  1.091      Density jump at z = 660 km  Δρ660  400  kg m−3  0.1212  Clapeyron slope at z = 660 km  Γ660  −2.5  MPa K−1  −0.0668  Post-perovskite density jump  ΔρpPv  62  kg m−3  0.0188  Clapeyron slope of post-perovskite  ΓpPv  13  MPa K−1  0.3474  CMB temperature  TCMB  3750  K  1.5  Density jump at CMB  ΔρCMB  5280  kg m−3  1.6            Viscosity law          Reference viscosity  η0  1.6 × 1021  Pa s  1.0  Viscosity ratio at z = 660 km  Δη660  30      Logarithmic thermal viscosity ratio  Ea  13.816–23.026      Logarithmic vertical viscosity ratio  Va  2.303      Compositional viscosity ratio  ΔηC  1–102      Surface yield stress  σ0  290  MPa  7.5 × 106  Yield stress gradient  $${\dot{\sigma }_z}$$  0.01  Pa/Pa  0.01  View Large Because the fluid properties (density, viscosity, thermal diffusivity, and thermal expansion) are allowed to vary throughout the system, the definition of the Rayleigh number is non-unique. In our simulations, we prescribed a reference Rayleigh number Ra0, defined at surface values of the thermodynamic parameters and reference viscosity η0. In all simulations shown here, Ra0 is set to 3.0 × 108, leading to an effective Rayleigh number (i.e. the Rayleigh number at the volume average viscosity) between 106 and 2.0 × 106, depending on the case. In thermochemical cases, the compositional field is modelled with a collection of 100 million particle tracers (200 million for cases with pPv), leading to an average number of tracers per cell of 15, which is enough to properly model entrainment (Tackley & King 2003). Tracers are of two types, modelling the regular mantle and primordial material, respectively, and are advected following a fourth-order Runge–Kutta method. At each time step, the compositional field is inferred from the concentration C of particles of primordial material in each cell. This concentration varies between 0 (corresponding to a cell filled with regular material only) and 1 (for a cell filled with primordial material only). The primordial material is initially distributed in a basal layer, the thickness of which is controlled by fixing the volume fraction of dense material, Xprim. In most of our calculations, we fixed Xprim to 3.5 per cent, corresponding to an upper estimate of the volume of LLSVPs (Hernlund & Houser 2008). In addition, we conducted calculations for Xprim = 2.5 per cent and Xprim = 7.0 per cent. The primordial material is assumed to be denser than the regular (pyrolitic) mantle, and the density contrast between the two materials is controlled by the buoyancy ratio. In our calculations, the buoyancy ratio is defined with respect to a reference density that increases with depth following a thermodynamical model of Earth's mantle,   $${B_z} = \frac{{\Delta {\rho _c}\left( z \right)}}{{{\alpha _S}\rho \left( z \right)\Delta {T_S}}},$$ (1)where Δρc(z) is the density contrast between dense and regular material, αS the surface thermal expansion, ρ(z) the reference density at depth z, and ΔTS the super-adiabatic temperature jump. In our models, the reference density increases with depth but Bz remains constant, implying that the chemical density contrast increases with depth proportionally to the reference density. With this definition, the buoyancy ratio required to obtain a chemical stratification of the system is smaller than when the surface density is taken as reference density. For equivalent chemical density contrast at the bottom of the system, however, results of the simulations do not differ substantially. We set Bz to 0.23 in most of our models. At the bottom of the system, and taking αS = 5.0 × 10−5 K−1, ρbot = 4950 kg m−3 and ΔTS = 2500 K, this leads to a density contrast between dense and regular material around 142 kg m−3, that is, about 2.5 per cent of PREM density (Dziewonski & Anderson 1981) at these depths, in agreement with the peak-to-peak chemical density anomalies mapped by normal modes data (Trampert et al.2004; Mosca et al.2012). Because the buoyancy ratio is certainly the most sensitive parameter controlling the dynamics of the system, we performed additional calculations with Bz between 0.15 and 0.70. This latter value allows maintaining a layer of primordial material that covers nearly all the CMB. Viscosity is allowed to vary with depth, temperature, and composition. An additional viscosity ratio Δη660 = 30 is added at the 660 km phase transition, and to avoid the formation of stagnant lid at the top of the system, we impose a yield stress. The viscosity is then fully described by   $$\eta \ = \frac{1}{{\frac{1}{{{\eta _b}}} + \frac{1}{{{\eta _Y}}}}},$$ (2)where   $${\eta _Y} = \frac{{{\sigma _0} + {{\dot{\sigma }}_z}P}}{{2\dot{e}}}\$$ (3)is the yield viscosity, and   \begin{eqnarray} {\eta _{\rm{b}}}\ \left( {z,T,C} \right) &=& {\eta _0}\ \left[ {1 + 29H\left( {z - 660} \right)} \right]\nonumber\\ &&\times \exp \left[ {{V_{\rm{a}}}\frac{z}{D} + {E_{\rm{a}}}\frac{{\Delta {T_{\rm{S}}}}}{{\left( {T + {T_{{\rm{off}}}}} \right)}} + {K_{\rm{a}}}C} \right].\end{eqnarray} (4) The yield viscosity (eq. 3) is defined from the yield stress, $${\sigma _Y} = \ {\sigma _0} + {\dot{\sigma }_z}P$$, and the second invariant of the stress tensor, $$\dot{e}.$$ The yield stress is set to σ0 at the surface and increases with pressure following a gradient $${\dot{\sigma }_z}$$. In eq. (4), η0 is a reference viscosity, H the Heaviside step function, z the depth, D the mantle thickness, ΔTS the superadiabatic temperature difference across the system and Toff the temperature offset, which is added to the temperature to reduce the viscosity jump across the top thermal boundary layer. The reference viscosity η0 is defined for the surface value of the reference adiabat (i.e. Tas = 0.64ΔTS), and at regular composition (C = 0). The viscosity variations with temperature are controlled by Ea, modelling the activation energy. To quantify the thermally induced increase of viscosity, we define a potential thermal viscosity ratio as ΔηT = exp(Ea). However, due to the adiabatic increase of temperature and to the temperature offset, which we fixed to Toff = 0.88ΔTS, the effective top-to-bottom thermal viscosity contrast is smaller than ΔηT by about two orders of magnitude. Here, we fixed Ea to 20.723 (corresponding to ΔηT = 109) in most simulations, and did additional calculations with values of Ea in the range 13.816–23.032 (106 ≤ ΔηT ≤ 1010). The viscosity variations with depth are controlled by Va, modelling the activation volume. In all simulations performed here, we fixed the value of Va to 2.303, leading to an increase of viscosity from top to bottom by a factor 300. This includes the viscosity jump at 660 km, but excludes the decrease due to adiabatic increase of temperature and the thermally induced increase in thermal boundary layers. The viscosity variations with composition are controlled by the parameter Ka, and the viscosity ratio between primordial and regular material (or chemical viscosity ratio) is given by ΔηC = exp(Ka). The viscosity of bridgmanite is larger than that of ferro-periclase by up to 3 orders of magnitude (Yamazaki & Karato 2001). If LLSVPs are enriched in bridgmanite by up to about 10 per cent, as suggested by probabilistic tomography (Trampert et al.2004; Mosca et al.2012), they should be intrinsically more viscous than surrounding mantle. Recent numerical experiments (Li et al.2014b; Deschamps et al.2015) further indicate that intrinsically more viscous primordial material increases the topography of reservoirs of dense material and the sharpness of their edges. The value of the chemical viscosity ratio depends on poorly constrained parameters, including the rock deformation model and the distribution of phases (bridgmanite and ferro-periclase) within the aggregate, and is thus difficult to estimate. In most of the thermochemical models presented in this study, we impose primordial material to be more viscous than regular material with ΔηC = 30. Finally, it should be pointed out that our viscosity law does not account for viscosity variations with grain size, which may influence mantle dynamics and evolution. The detailed effects of grain size remain a matter of debate. For diffusion creep, viscosity is usually assumed to increase with grain size following a power law with exponent in the range 2–3 (e.g. Korenaga & Karato 2008). Large grain size further favours the transition from diffusion to dislocation creep regimes, and it has been suggested that this transition may occur for grain size as low as 0.1–1 mm (Li et al.1996). Due to various physical and chemical processes, grain size may strongly vary in space and time. For instance, because in early Earth grain growth occurred in a partially molten environment, initial grain size may have been rather large, while recrystallization, in particular during the ringwoodite to bridgmanite phase transition, leads to smaller grain size. Regions with larger grain size in primordial mantle had larger viscosity, and depending on the vigour of convection, they may have mixed efficiently or, on the contrary, poorly with regions of lower viscosity (Solomatov & Reese 2008). If thermochemical reservoirs in our models resulted from early mantle differentiation, their intrinsic viscosity may differ from that of the regular mantle by up to a few orders of magnitude. A simple way to account for this effect would then be to systematically vary the chemical density contrast, ΔηC, up to very large values, for example, 103 and more. This, however, neglects recrystallization effects that may affect the primordial material during its circulation in the mantle. Fig. 3, plotting viscosity profiles along the initial temperature profile, illustrates the radial viscosity variations in our models when thermal and chemical effects are taken into account. The effect of intrinsically more viscous (ΔηC = 30) dense material at the bottom of the mantle is denoted by the dashed curve. Note that for Ea = 20.723, which is the case of most of our models, the regular increase with depth (Va = 2.303) in the lower mantle is nearly compensated by the decrease due to adiabatic increase of temperature. Figure 3. View largeDownload slide Radial viscosity profiles along initial temperature profile for purely thermal models or thermochemical models with chemical viscosity contrast ΔηC = 1 (plain curves), and thermochemical models with ΔηC = 30 (dotted curve). Two values of the thermal viscosity contrast (controlled by the parameter Ea) are shown. Figure 3. View largeDownload slide Radial viscosity profiles along initial temperature profile for purely thermal models or thermochemical models with chemical viscosity contrast ΔηC = 1 (plain curves), and thermochemical models with ΔηC = 30 (dotted curve). Two values of the thermal viscosity contrast (controlled by the parameter Ea) are shown. The phase change at 660 km is modelled with a discontinuous phase transition controlled by defining a point on the phase boundary and a Clapeyron slope, Γ660. Here, we imposed z = 660 km and T = 1900 K as anchor point, and Γ660 = −2.5 MPa K−1. The density contrast at the phase transition is Δρ660 = 400 kg m−3 and is scaled with the surface density. In cases that include pPv phase, the pPv phase transition is modelled with a reference temperature of 2700 K at 2700 km depth. Lateral deviations in the transition depth are determined using the phase function approach of Christensen & Yuen (1985). The Clapeyron slope and the density contrast are set to ΓpPv = 13 MPa K−1 (Tateno et al.2009) and ΔρpPv = 62 kg m−3 (corresponding to a relative contrast of ∼1.0 per cent, Murakami et al.2004), respectively. Determining the stability field of pPv further requires fixing the temperature at the CMB, TCMB. Based on the numerical simulations of Li et al. (2015), which pointed out that values of TCMB around 3750 K describe best lower mantle structure and dynamics, we assumed TCMB = 3750 K. The dynamic topography at the CMB, hCMB, is calculated from the normal stress, σzz, induced by the flow on this boundary, and takes into account self-gravitational effects following Zhang & Christensen (1993). At each point of the CMB, dynamic topography is given by   $${h_{{\rm{CMB}}}} = \frac{{{\sigma _{zz}} + {{\rm{\Phi }}_{{\rm{CMB}}}}\Delta {\rho _{{\rm{CMB}}}}}}{{\Delta {\rho _{{\rm{CMB}}}}g}},$$ (5)where ΔρCMB is the density difference between the mantle and the outer core, ΦCMB the perturbation of gravitational potential at CMB, and g the acceleration of gravity. The normal stress is given by   $$\ {\sigma _{zz}} = \ 2\eta \left( {\frac{{\partial {v_r}}}{{\partial r}} - \frac{1}{3}\nabla \cdot {\boldsymbol{v}}} \right),$$ (6) where v is velocity and vr its vertical component. Note that the divergence of velocity in the right hand side of eq. (6) is due to compressibility. Gravitational potential perturbations are obtained by solving Poisson's equation for a distribution of density anomalies (with respect to a reference state) that includes thermal and chemical effects, and variations in CMB topography. For comparison, we also calculated dynamic topography neglecting self-gravitation, but only found small differences, mostly in amplitude, with the topography accounting for self-gravity. For a given set of parameters, the rms amplitude of dynamic topography is, on average, smaller by 5 per cent if self-gravitational effects are included. It is usually assumed that the mechanical boundary condition, free or free-slip, does not affect the boundary topography, because this topography is very small compared to the lateral lengthscale, and therefore slopes are also very small. This has been well tested for the surface boundary (Crameri et al.2012). For large slopes, such as at subduction zones trenches, differences may however be important, and a free surface may help single-sided subduction (e.g. Crameri & Tackley 2014). The initial temperature distribution is built by adding small random perturbations to an adiabatic radial profile with thin superadiabatic boundary layers at the top and bottom of the shell. Purely thermal experiments are run up to a dimensional time of 4.5 Gyr. Thermochemical experiments start with a long (up to 1.5–2.0 Gyr) transient phase during which the bottom layer of dense material heats up, and are run up to a dimensional time of 15 Gyr. Models do not reach stationary state, that is, the location of plumes and downwellings vary with time. However, after a period of time that depends on each case, average velocity and surface and bottom heat fluxes stabilize and oscillate around constant values, thus defining a dynamic quasi-equilibrium. All the snapshots presented and discussed in Section 3 are taken in the quasi-equilibrium stage. In thermochemical models, the detailed shape of thermochemical piles varies with time, but the selected snapshots are generally representative of the thermochemical structure obtained for each case. Snapshots are taken at 2.7 Gyr for purely thermal cases, and 11.2 or 15.0 Gyr for thermochemical cases. Note that the experiments discussed here are not designed to model the detailed evolution of the Earth's mantle, for which accurate initial conditions are not yet known. This would further require adding time-decreasing radioactive heating and CMB temperature in our models. Instead, the aim is to obtain various characteristic thermochemical structures, for example, plume clusters or thermochemical piles, to estimate the influence of these structures on CMB dynamic topography. Thus, time indication in our experiments should not be used to discriminate between possible models of mantle dynamics, or to interpret early mantle evolution. 3 RESULTS Using the setup described in Section 2, we performed a series of 21 numerical experiments with various values of the buoyancy ratio, volume fraction of dense material, and thermal and chemical viscosity ratios (Table 2). All models are in the mobile-lid regime, with horizontal surface velocities of around 1 cm yr−1. The surface heat flux is in the range 35–40 mW m−2, depending on the model. These are rather low values compared to heat flux measured at the surface of the Earth, the average of which is around 86 mW m−2 (Davies 2013). It should however be kept in mind that the observed surface heat flux is likely overestimating the heat flux at the top of the mantle because it includes a contribution related to crustal heat production, which is difficult to estimate, but likely around 8 TW (Jaupart et al.2015). In cratonic regions, for instance, this contribution may be large and the mantle heat flux could be as low as 15–20 mW m−2 (Jaupart & Mareschal 2007). In addition, because our models do not include plate tectonics, they do not distinguish between continental areas and oceanic regions, where observed heat flux is usually highest. Finally, our models do not include secular cooling of the mantle, which contributes up to half of Earth's surface heat loss (Jaupart et al.2015). The bottom heat flux is around 80–100 mW m−2 depending on the model, equivalent to a total power around 12–15 TW, that is, in good agreement with the expected core heat flow, in the range 5–17 TW (Jaupart et al.2015). Table 2. Run properties and corresponding output core–mantle boundary (CMB) topography. The prefix T and TC denote purely thermal and thermochemical cases, respectively. Bz is the buoyancy ratio, Xprim the volume fraction of dense material, ΔηC the intrinsic (chemical) viscosity ratio, and Ea the logarithmic thermal viscosity contrast (eq. 3). Other properties are listed in Table 1 and are the same for all runs. CMB topography statistics includes rms of global topography, and average topography 〈h〉 in primordial reservoirs or plume clusters. For primordial reservoirs, we distinguish the average depth of their interiors 〈h−〉, and the average elevation of their surrounding ridges, 〈h+〉. All topographies are given in km. Run  Bz  Xprim  ΔηC  Ea  pPv  Resolution  hCMB (km)                rms  〈h−〉  〈h+〉  T1  –  –  –  20.723  no  512 × 256 × 64  2.7  –  1.9  T2  –  –  –  13.816  no  512 × 256 × 64  2.7  –  2.9  T3  –  –  –  23.026  no  512 × 256 × 64  2.4  –  1.6  T1-pPv  –  –  –  20.723  yes  512 × 256 × 128  2.4  –  2.5  T2-pPv  –  –  –  13.816  yes  512 × 256 × 128  3.1  –  3.9  TC1  0.23  3.5  30  20.723  no  512 × 256 × 64  2.8  −1.5  1.5  TC2  0.23  3.5  30  13.816  no  512 × 256 × 64  2.6  −1.7  1.8  TC3  0.23  3.5  30  23.026  no  512 × 256 × 64  2.7  −1.3  1.4  TC4  0.15  3.5  30  20.723  no  512 × 256 × 64  3.1  −0.6  1.3  TC5  0.50  3.5  30  20.723  no  512 × 256 × 64  3.0  −2.4  2.0  TC6  0.70  3.5  30  20.723  no  512 × 256 × 64  3.4  −2.5  2.8  TC7  0.23  3.5  1  20.723  no  512 × 256 × 64  2.0  −0.7  1.0  TC8  0.23  3.5  10  20.723  no  512 × 256 × 64  2.4  −1.2  1.4  TC9  0.23  3.5  100  20.723  no  512 × 256 × 64  2.5  −1.7  1.5  TC10  0.23  2.5  30  20.723  no  512 × 256 × 64  2.4  −1.7  1.3  TC11  0.23  7.0  30  20.723  no  512 × 256 × 64  2.6  −1.3  1.4  TC1-pPv  0.23  3.5  30  20.723  yes  512 × 256 × 128  2.8  −1.9  1.5  TC2-pPv  0.23  3.5  30  13.816  yes  512 × 256 × 128  3.8  −2.2  2.4  TC4-pPv  0.15  3.5  30  20.723  yes  512 × 256 × 128  3.0  −0.9  1.3  TC6-pPv  0.70  3.5  30  20.723  yes  512 × 256 × 128  3.9  −2.5  2.7  TC7-pPv  0.23  3.5  1  20.723  yes  512 × 256 × 128  2.7  −0.9  1.2  Run  Bz  Xprim  ΔηC  Ea  pPv  Resolution  hCMB (km)                rms  〈h−〉  〈h+〉  T1  –  –  –  20.723  no  512 × 256 × 64  2.7  –  1.9  T2  –  –  –  13.816  no  512 × 256 × 64  2.7  –  2.9  T3  –  –  –  23.026  no  512 × 256 × 64  2.4  –  1.6  T1-pPv  –  –  –  20.723  yes  512 × 256 × 128  2.4  –  2.5  T2-pPv  –  –  –  13.816  yes  512 × 256 × 128  3.1  –  3.9  TC1  0.23  3.5  30  20.723  no  512 × 256 × 64  2.8  −1.5  1.5  TC2  0.23  3.5  30  13.816  no  512 × 256 × 64  2.6  −1.7  1.8  TC3  0.23  3.5  30  23.026  no  512 × 256 × 64  2.7  −1.3  1.4  TC4  0.15  3.5  30  20.723  no  512 × 256 × 64  3.1  −0.6  1.3  TC5  0.50  3.5  30  20.723  no  512 × 256 × 64  3.0  −2.4  2.0  TC6  0.70  3.5  30  20.723  no  512 × 256 × 64  3.4  −2.5  2.8  TC7  0.23  3.5  1  20.723  no  512 × 256 × 64  2.0  −0.7  1.0  TC8  0.23  3.5  10  20.723  no  512 × 256 × 64  2.4  −1.2  1.4  TC9  0.23  3.5  100  20.723  no  512 × 256 × 64  2.5  −1.7  1.5  TC10  0.23  2.5  30  20.723  no  512 × 256 × 64  2.4  −1.7  1.3  TC11  0.23  7.0  30  20.723  no  512 × 256 × 64  2.6  −1.3  1.4  TC1-pPv  0.23  3.5  30  20.723  yes  512 × 256 × 128  2.8  −1.9  1.5  TC2-pPv  0.23  3.5  30  13.816  yes  512 × 256 × 128  3.8  −2.2  2.4  TC4-pPv  0.15  3.5  30  20.723  yes  512 × 256 × 128  3.0  −0.9  1.3  TC6-pPv  0.70  3.5  30  20.723  yes  512 × 256 × 128  3.9  −2.5  2.7  TC7-pPv  0.23  3.5  1  20.723  yes  512 × 256 × 128  2.7  −0.9  1.2  View Large Table 2 further lists the rms amplitude of dynamic topography for each model. To better quantify the dynamic topography caused by plumes, downwellings, and reservoirs of dense material, we define several regions above the CMB. The plume, or plume cluster, region is defined as the area where the temperature is larger than Tplume = Tm + cplume(Tmax − Tm), where Tmax and Tm are the maximum and average temperatures above the CMB, respectively, and cplume a constant that we fixed to 0.5. The region where downwellings accumulate is defined as the region where the temperature is lower than Tslab = Tm − cslab(Tm − Tmin), where Tmin is the minimum temperature above the CMB, and cslab a constant that we fixed, again, to 0.5. In thermochemical cases, reservoirs of dense material are defined by the region where the concentration in dense particles, C, is larger than 0.5. Note that because the transition between the dense reservoir and its surroundings is very sharp, the size and extension of the reservoir do not vary with threshold values of C larger than 0.3. By contrast, fixing cplume = cslab = 0.5 is rather subjective, and changing these values leads to variations in the extents of plume and downwelling regions, inducing in turn variations in the estimated average topography within these regions. For instance, reducing cslab increases the area of downwelling regions by including regions with smoother topography, thus decreasing the average topography. Importantly, while the extent of plume and slab areas may influence the estimated average CMB topography in these regions, it does not affect the trend, that is, the fact that plumes induce elevations in the CMB, and slabs depressions. Figs 4 and 5 show snapshots of the temperature, composition and CMB dynamic topography for cases T1 and TC1, respectively, which we take as reference cases for purely thermal and thermochemical models. In case T1, plumes are generated at the CMB and are grouping in two main clusters (Figs 4a and b). Plumes are interconnected through hot ridges, and the typical angular distance between two plumes is 20°. Cold downwellings reach the CMB in the wide regions located between plume clusters. At their base, plumes generate upwards normal stress, and therefore positive dynamic topography (Figs 4c and d). In the plume clusters region, the average amplitude of dynamic topography is 1.9 km. By contrast, as they reach the CMB, downwellings induce negative dynamic topography, that is, depressions of the CMB. These depressions are strongly localized, but have large amplitude. In case T1, the slab stacking region covers less than 0.5 per cent of the total CMB area, but within this region the rms amplitude of dynamic topography is 17 km. A distribution histogram of the CMB topography (Fig. 6a) indicates that throughout most of the CMB (∼90 per cent of its total surface) dynamic topography ranges between −2.0 and 3.0 km. Downwellings slightly skew the distribution towards topography deeper than −5.0 km. Note that for convenience, the histogram is truncated to the interval [−10; 10] km, implying that topography caused by downwellings is not entirely represented. However, for topography deeper than −10 km the frequency is very small, less than 10−3 for a 1 km-stacking bin. Figure 4. View largeDownload slide Snapshot of a purely thermal case without post-perovskite (T1). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.72. (c,d) Core–mantle boundary dynamic topography (colour scale). The orange and blue dotted lines in plots (c) and (d) indicate the limits of the plumes and downwelling regions, respectively. Black dotted lines in plots (b) and (c) denote the path along which topographic profiles in Figs 15(a) and (b) are drawn. Figure 4. View largeDownload slide Snapshot of a purely thermal case without post-perovskite (T1). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.72. (c,d) Core–mantle boundary dynamic topography (colour scale). The orange and blue dotted lines in plots (c) and (d) indicate the limits of the plumes and downwelling regions, respectively. Black dotted lines in plots (b) and (c) denote the path along which topographic profiles in Figs 15(a) and (b) are drawn. Figure 5. View largeDownload slide Snapshot of a thermochemical case without post-perovskite (TC1). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.67. (c,d) Isosurface of the composition for C = 0.5. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in plots (e) and (f) indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. Black dotted lines in plots (e) and (f) denote the path along which the topographic profile in Fig. 15(c) is drawn. Figure 5. View largeDownload slide Snapshot of a thermochemical case without post-perovskite (TC1). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.67. (c,d) Isosurface of the composition for C = 0.5. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in plots (e) and (f) indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. Black dotted lines in plots (e) and (f) denote the path along which the topographic profile in Fig. 15(c) is drawn. Figure 6. View largeDownload slide Distribution histograms of dynamic topography at the CMB for cases T1 (a), TC1 (b), T1-pPv (c) and TC1-pPv (d) (see Table 2 for the properties of these cases). Frequency is normalized to the total area of the CMB. For convenience, histograms are truncated to the interval [−10; 10] km, implying that topography caused by downwellings is not entirely represented. Figure 6. View largeDownload slide Distribution histograms of dynamic topography at the CMB for cases T1 (a), TC1 (b), T1-pPv (c) and TC1-pPv (d) (see Table 2 for the properties of these cases). Frequency is normalized to the total area of the CMB. For convenience, histograms are truncated to the interval [−10; 10] km, implying that topography caused by downwellings is not entirely represented. In case TC1, the thermochemical structure at the bottom of the system is dominated by two large antipodal reservoirs of primordial material culminating about 500 km above the CMB (Figs 5c and d). Reservoirs are hotter than average, and plumes are generated at their edges and in their central areas (Figs 5a and b). Large downwellings reach the CMB in regions located in between the two thermochemical reservoirs. Each reservoir of dense material induces a large depression throughout most of its base, surrounded by a thin ridge of positive topography on its edges (Figs 5e and f). On average, the depth of the depression and the height of the ridge are equal to −1.5 km and 1.5 km, respectively (Table 2). As in case T1, downwellings induce deep local depressions with rms amplitude of 6.2 km. Interestingly, the CMB topography distribution is different from that obtained for purely thermal model T1 (Fig. 6b), and now follows a bimodal distribution, with a peak between −2.0 and −1.0 km, corresponding to the interior of the reservoir of primordial material, and a second, less pronounced, peak between 1.0 and 2.0 km, corresponding to its edges. Overall, most of the CMB area (>90 per cent) has a topography between −3.0 and 4.0 km. As in case T1, the distribution is slightly skewed by downwellings towards topography deeper than −5.0 km. Figs 7 and 8 show snapshots for cases T1-pPv and TC1-pPv, which are similar to T1 and TC1, respectively, except that they include the pPv phase at the bottom of the system. Adding the pPv does not substantially modify the convection patterns obtained for cases T1 and TC1, except that downwellings are globally less well developed. For the CMB temperature we imposed, pPv is stable outside the plume clusters regions for the purely thermal case (T1-pPv; Figs 7c and d) and outside the reservoir of primordial material for the thermochemical case (TC1-pPv; Figs 8e and f). Note that in these models, the viscosities of pPv and bridgmanite are similar. Weak pPv, with a viscosity smaller than that of bridgmanite by three orders of magnitude, may slightly modify the convection pattern, in particular for the thermochemical cases (e.g. Li et al.2015). The effects on dynamic topography are similar to those obtained without the pPv phase. In purely thermal models, plume clusters induce positive topography with an average height of ∼2.4 km (Figs 7e and f), while in thermochemical models, reservoirs of primordial material cause a depression of the CMB in their interiors and a ridge on their edges, with average depth and height equal to −1.9 km and 1.5 km, respectively (Figs 8g and h). Again, downwellings cause deep local troughs. Distributions of topography are slightly more dispersed, with most of the topography ranging between −4.0 km and 4.0 km for model T1-pPv, and between −4.0 and 5.0 km for model TC1-pPv, but are globally unchanged (Figs 6c and d). Figure 7. View largeDownload slide Snapshot of a purely thermal case with post-perovskite (T1-pPv). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.80. (c,d) Stability field of the post-perovskite phase. (e,f) Core–mantle boundary dynamic topography (colour scale). The orange and blue dotted lines in plots (e) and (f) indicate the limits of the plumes and downwelling regions, respectively. Figure 7. View largeDownload slide Snapshot of a purely thermal case with post-perovskite (T1-pPv). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.80. (c,d) Stability field of the post-perovskite phase. (e,f) Core–mantle boundary dynamic topography (colour scale). The orange and blue dotted lines in plots (e) and (f) indicate the limits of the plumes and downwelling regions, respectively. Figure 8. View largeDownload slide Snapshot of a thermochemical case with post-perovskite (TC1-pPv). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.67. (c,d) Isosurface of the composition for C = 0.5. (e,f) Stability field of the post-perovskite phase. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in plots (g) and (h) indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. Figure 8. View largeDownload slide Snapshot of a thermochemical case with post-perovskite (TC1-pPv). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.67. (c,d) Isosurface of the composition for C = 0.5. (e,f) Stability field of the post-perovskite phase. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in plots (g) and (h) indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. Increasing the volume fraction of dense material, Xprim, obviously increases the size of reservoirs of primordial material, in particular the fraction of CMB area they cover, but does not substantially modify their stability (Li et al.2014b). All other parameters being equal to those of case TC1, our calculations for Xprim = 2.5 per cent (case TC10) and Xprim = 7.0 per cent (case TC11) further indicate that dynamic topography is only slightly affected by a change in Xprim (Table 2). Other parameters of the model, in particular the buoyancy ratio and the thermal and compositional viscosity contrasts, induce more substantial changes on the flow pattern. We now further explore the consequences of these changes on the CMB dynamic topography. 3.1 Buoyancy ratio Fig. 9 shows snapshots of models with buoyancy ratio Bz = 0.15 (case TC4) and Bz = 0.7 (case TC6), corresponding to chemical density contrasts of 93 and 433 kg m−3, respectively. In case TC4, the chemical density contrast is too small to prevent strong entrainment of primordial material, leading to the formation of narrow, elevated thermochemical piles with sharp edges (Figs 9a–d). The dynamic topography induced by these piles is strongly attenuated compared to that of case TC1. The average depth of central depression reaches only 0.6 km, and the average height of surrounding ridge 1.3 km (Table 2). As observed by Lassak et al. (2007), the depressions induced by slabs are much larger, with, in case TC4, an average depth of 11.6 km. The distribution of topography is no longer bimodal but peaks between 1.0 and 2.0 km (Fig. 10a), similar to the distributions obtained for purely thermal cases. In case TC6, on the contrary, the large chemical density contrast strongly inhibits entrainment of primordial material, leading to a nearly stratified system in which the core is almost fully covered by a layer of primordial material (Figs 9e–h). Compared to case TC1, dynamic topography is more pronounced, with average depth of depression and height of ridge around −2.5 and 2.8 km, respectively (Table 2). The negative CMB topography caused by slabs, on the contrary, is attenuated, with an average depth of 4.2 km. In addition, topographic highs around the reservoir reach large values, 10 km and more, strongly skewing the distribution of topography towards high values (Fig. 10b). Figure 9. View largeDownload slide Snapshots of thermochemical cases TC4 (a–d) with buoyancy ratio B = 0.15, and TC6 (e–h), with buoyancy ratio B = 0.70. Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the composition for C = 0.65. (c,d) Core–mantle boundary dynamic topography (colour scale). (e,f) Isosurface of the composition for C = 0.5. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in topography plots indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. Black dotted lines in plots (c) and (d) denote the path along which the topographic profile in Fig. 15(d) is drawn. Figure 9. View largeDownload slide Snapshots of thermochemical cases TC4 (a–d) with buoyancy ratio B = 0.15, and TC6 (e–h), with buoyancy ratio B = 0.70. Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the composition for C = 0.65. (c,d) Core–mantle boundary dynamic topography (colour scale). (e,f) Isosurface of the composition for C = 0.5. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in topography plots indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. Black dotted lines in plots (c) and (d) denote the path along which the topographic profile in Fig. 15(d) is drawn. Figure 10. View largeDownload slide Distribution histograms of dynamic topography at the CMB for six thermochemical cases, (a) TC4, (b) TC6, (c) TC2, (d) TC3, (e) TC7 and (f) TC9 (see Table 2 for the properties of each case). Frequency is normalized to the total area of the CMB. For convenience, histograms are truncated to the interval [−10; 10] km, implying that topography caused by downwellings is not entirely represented. Figure 10. View largeDownload slide Distribution histograms of dynamic topography at the CMB for six thermochemical cases, (a) TC4, (b) TC6, (c) TC2, (d) TC3, (e) TC7 and (f) TC9 (see Table 2 for the properties of each case). Frequency is normalized to the total area of the CMB. For convenience, histograms are truncated to the interval [−10; 10] km, implying that topography caused by downwellings is not entirely represented. Overall, increasing the buoyancy ratio of primordial material increases the amplitude of dynamic topography generated by thermochemical piles and reservoirs. Additional cases including the pPv phase transition (cases TC4-pPv and TC6-pPv in Table 2) indicate that this trend is not altered when pPv is present. 3.2 Thermal viscosity contrast Thermal viscosity contrast has a strong influence on the shape and stability of reservoirs of primordial material (Deschamps & Tackley 2009; Li et al.2014b). These reservoirs are hotter, and therefore less viscous than their surroundings. This viscosity contrast increases with Ea, which in turn reduces the mechanical coupling of dense material and their entrainment by plumes. Figs 11(a)–(d) show a snapshot of thermochemical case TC2 with thermal viscosity contrast Ea = 13.816, corresponding to a maximum (or potential) bottom-to-top viscosity contrast, ΔηT, of 106. Compared to case TC1, for which Ea = 20.723 and ΔηT = 109, the shape and stability of the reservoirs of dense material are substantially altered. Reservoirs have smaller length-scale, larger topography, and sharper edges, as already observed in the simulations of Li et al. (2014b). Changes in thermal viscosity contrast strongly affect the amplitude of CMB dynamic topography beneath reservoirs. Viscosity in reservoirs decreases with increasing thermal viscosity contrast, leading to flatter topography (Table 2). The average depth of depression and height of surrounding ridge are around −1.7 and 1.8 km for case TC2 (Ea = 13.816), and −1.3 and 1.4 km for case TC3 (Ea = 23.026). We observe a similar effect for purely thermal models, that is, the topography induced by plumes decreases in amplitude with increasing thermal viscosity contrast (Table 2). In our models, the average elevation drops from 2.9 km for model T2 (Ea = 13.816) to 1.6 km for model T3 (Ea = 23.026). Interestingly, in thermochemical models, reducing the thermal viscosity contrast changes the distribution of topography. Contrary to cases TC1 (Fig. 6b) and TC3 (Fig. 10d), the distribution of topography for case TC2 (Fig. 10c) is not bimodal. In addition, it is slightly more dispersed, with most of the topography being in the range −5.0 and 5.0 km. Again, the addition of pPv does not modify this trend, as indicated by comparison between cases TC1-pPv and TC2-pPv, and T1-pPv and T2-pPv (Table 2). Figure 11. View largeDownload slide Snapshots of thermochemical cases TC2 (a–d), with Ea = 13.816, and TC7 (e–h), with ΔηC = 1. Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the composition for C = 0.65. (c,d) Core–mantle boundary dynamic topography (colour scale). (e,f) Isosurface of the composition for C = 0.5. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in topography plots indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. Figure 11. View largeDownload slide Snapshots of thermochemical cases TC2 (a–d), with Ea = 13.816, and TC7 (e–h), with ΔηC = 1. Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the composition for C = 0.65. (c,d) Core–mantle boundary dynamic topography (colour scale). (e,f) Isosurface of the composition for C = 0.5. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in topography plots indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. 3.3 Chemical viscosity contrast Increasing the chemical viscosity ratio between primordial and regular materials, ΔηC, has the opposite effect to that of increasing the thermal viscosity contrast. The shape and stability of reservoirs are only slightly modified. Reservoirs are flatter if primordial and regular material have similar intrinsic viscosity (ΔηC = 1; case TC7, Figs 11e and f) than if primordial material is more viscous (e.g. case TC1, Figs 5c and d). The amplitude of CMB dynamic topography beneath the reservoirs increases with increasing intrinsic viscosity. Peak-to-peak amplitude rises from 1.8 km for case TC7 (ΔηC = 1) to 3.2 km for case TC9 (ΔηC = 102) (Table 2). For ΔηC = 1, most of CMB topography ranges between −2.0 and 3.0 km with a peak between 0 and −1.0 km (Fig. 10e), while for ΔηC = 10 or more, the CMB topography is more dispersed and follows a bimodal distribution (Figs 6b and 10f). Finally, comparison between cases TC7 and TC7-pPv shows that the addition of pPv does not alter the influence of ΔηC. 4 COMPARISON WITH SEISMOLOGICAL CONSTRAINTS A key question is whether CMB dynamic topography, provided it can be mapped from seismic data, can discriminate between purely thermal and thermochemical models of Earth's mantle. CMB topography predicted by our models strongly differ both in amplitude and pattern depending on the properties of each model, suggesting that CMB topography may bring important insights on lowermost mantle structure. Available seismic models on CMB topography, however, are limited to long wavelengths, up to spherical harmonic degree 4, for which purely thermal and thermochemical models may have similar signatures (Lassak et al.2010). In this section, we examine more in details the long-wavelength structure of our CMB topography models. 4.1 Amplitude: comparison with observed CMB topography Table 3 and Fig. 12(a) quantitatively compare CMB topography predicted by our models with available seismological constraints for different spherical harmonic filters. Normal modes data further suggest that the amplitude of spherical harmonic degree l = 2 of topography should be smaller than 5 km (Koelemeijer et al.2012). Our models fit this observation well (Table 3), with the notable exception of thermochemical models with very strong density contrast (TC5, TC6 and TC6-pPv), around 300 kg m−3 and more. These later models induce a strong chemical layering (Figs 9e and f) that is not seen by seismic velocity data, and are therefore very unlikely. Note, also, that the peak-to-peak amplitude of CMB topography predicted by purely thermal model T1 is slightly larger than the 5 km limit suggested by seismic normal modes, and that models with small density contrast (TC4 and TC4-pPv), around 90 kg m−3 and less, are close to this limit. Several studies have attempted mapping long-wavelength (up to spherical harmonic degree l = 4) CMB topography (Morelli & Dziewonski 1987; Doornbos & Hilton 1989; Koper et al.2003; Sze & van der Hilst 2003; Tanaka 2010). The resulting maps differ strongly both in pattern and amplitude, a major difficulty being that CMB topography trades off with lowermost mantle seismic velocity anomalies, and upper outer core structure. Studies including all spherical harmonic degrees up to l = 4 have peak-to-peak amplitude between 3.5 km (Koper et al.2003) and 12 km (Morelli & Dziewonski 1987). In general, for this spherical harmonic filter, the dynamic topography predicted by purely thermal models or thermochemical models with small density contrast have peak-to-peak amplitude ranging from 6 km to 11 km, somewhat larger than that of thermochemical models with moderate density contrast (from ∼100 to 200 kg m−3), which ranges between 3.5 and 7.5 km, depending on input parameters (Table 3 and Fig. 12a). Using degrees l = 2 and l = 4 only, Tanaka (2010) found a peak-to-peak amplitude of 5 km, in fair agreement with the peak-to-peak topography predicted by thermochemical models with moderate density contrast (between 2.5 and 5.5 km), but slightly lower than that of purely thermal models, which are around 6 km and more, except for model T3 (Table 3 and Fig. 12a). Figure 12. View largeDownload slide (a) Peak-to-peak amplitude of dynamic topography predicted by our models of convection (colour code) for three spherical harmonic filters (l = 2, l = 2 and 4 and l = 0–4). For comparison, amplitudes observed by Morelli & Dziewonski (1987) (MD87), Doornbos & Hilton (1989) (DH89), Sze & van der Hilst (2003) (SV03), Koper et al. (2003) (KPF03), Tanaka (2010) (T10), and Koelemeijer et al. (2012) (KDT12) are also indicated. (b) Correlation coefficient between long-wavelength (l = 0–4) shear-velocity anomalies and CMB dynamic topography predicted by each model of convection. See also Table 3 for numerical values. Figure 12. View largeDownload slide (a) Peak-to-peak amplitude of dynamic topography predicted by our models of convection (colour code) for three spherical harmonic filters (l = 2, l = 2 and 4 and l = 0–4). For comparison, amplitudes observed by Morelli & Dziewonski (1987) (MD87), Doornbos & Hilton (1989) (DH89), Sze & van der Hilst (2003) (SV03), Koper et al. (2003) (KPF03), Tanaka (2010) (T10), and Koelemeijer et al. (2012) (KDT12) are also indicated. (b) Correlation coefficient between long-wavelength (l = 0–4) shear-velocity anomalies and CMB dynamic topography predicted by each model of convection. See also Table 3 for numerical values. Table 3. Columns 2–7 list the root mean square (rms) and peak-to-peak amplitude of long-wavelength dynamic topography for four different spherical harmonic filters. For comparison, the peak-to-peak amplitudes of CMB topography observed by Morelli & Dziewonski (1987) (MD87), Doornbos & Hilton (1989) (DH89), Sze & van der Hilst (2003) (SV03), Koper et al. (2003) (KPF03), Tanaka (2010) (T10) and Koelemeijer et al. (2012) (KDT12) are also indicated. All topographies are given in kilometres. Column 8 lists the fraction of CMB area with topography between −4.0 and 4.0 km, x4km. Observed value is from Garcia & Souriau (2000) (GS00). The last column lists the correlation between long-wavelength (l = 0–4) dynamic topography and shear velocity anomalies, χh-dlnVs, predicted by each model of convection.     Rms    Peak  to peak  amplitude  x4km  χh-dlnVs    l = 2  l = 2, 4  l = 0–4  l = 2  l = 2, 4  l = 0–4    l = 0–4  Run  T1  1.3  1.3  1.6  5.1  6.6  8.5  0.94  −0.86  T2  0.7  1.9  2.1  2.0  7.1  8.7  0.89  −0.94  T3  0.8  0.9  1.2  2.7  4.1  5.8  0.95  −0.54  T1-pPv  1.3  1.3  1.5  4.5  5.9  7.7  0.91  −0.83  T2-pPv  0.8  2.1  2.3  3.4  9.8  11.5  0.81  −0.98  TC1  0.7  0.8  1.0  2.3  3.5  4.4  0.88  0.71  TC2  0.6  0.8  1.0  2.3  3.1  4.5  0.89  −0.13  TC3  0.7  1.2  1.3  2.9  5.7  6.7  0.92  0.33  TC4  1.1  1.2  1.6  4.4  5.4  7.6  0.89  −0.68  TC5  1.3  1.6  2.5  5.2  7.5  9.8  0.82  0.93  TC6  2.5  2.3  2.2  7.3  9.1  8.2  0.75  0.88  TC7  0.2  0.6  0.7  0.9  2.7  3.5  0.95  0.27  TC8  0.4  0.5  0.9  1.2  2.9  4.0  0.93  0.58  TC9  0.5  0.7  1.0  1.7  3.0  4.0  0.90  0.82  TC10  0.4  0.7  0.9  1.8  3.4  5.1  0.93  0.58  TC11  1.0  1.2  1.5  3.4  5.8  6.8  0.92  0.65  TC1-pPv  1.1  1.3  1.6  4.6  5.2  7.1  0.87  0.54  TC2-pPv  0.4  0.8  1.5  1.5  3.9  6.8  0.72  −0.29  TC4-pPv  1.1  1.2  1.5  4.8  5.6  8.6  0.89  −0.50  TC6-pPv  2.3  2.4  2.9  9.1  11.5  14.9  0.70  0.89  TC7-pPv  0.6  0.9  1.2  2.2  4.3  6.0  0.91  0.05  Observed  KDT12        ≤5.0          T10          5.0        MD87            12.0      DH89            8.0      SV03            4.0      KPF03            3.5      GS00              0.95        Rms    Peak  to peak  amplitude  x4km  χh-dlnVs    l = 2  l = 2, 4  l = 0–4  l = 2  l = 2, 4  l = 0–4    l = 0–4  Run  T1  1.3  1.3  1.6  5.1  6.6  8.5  0.94  −0.86  T2  0.7  1.9  2.1  2.0  7.1  8.7  0.89  −0.94  T3  0.8  0.9  1.2  2.7  4.1  5.8  0.95  −0.54  T1-pPv  1.3  1.3  1.5  4.5  5.9  7.7  0.91  −0.83  T2-pPv  0.8  2.1  2.3  3.4  9.8  11.5  0.81  −0.98  TC1  0.7  0.8  1.0  2.3  3.5  4.4  0.88  0.71  TC2  0.6  0.8  1.0  2.3  3.1  4.5  0.89  −0.13  TC3  0.7  1.2  1.3  2.9  5.7  6.7  0.92  0.33  TC4  1.1  1.2  1.6  4.4  5.4  7.6  0.89  −0.68  TC5  1.3  1.6  2.5  5.2  7.5  9.8  0.82  0.93  TC6  2.5  2.3  2.2  7.3  9.1  8.2  0.75  0.88  TC7  0.2  0.6  0.7  0.9  2.7  3.5  0.95  0.27  TC8  0.4  0.5  0.9  1.2  2.9  4.0  0.93  0.58  TC9  0.5  0.7  1.0  1.7  3.0  4.0  0.90  0.82  TC10  0.4  0.7  0.9  1.8  3.4  5.1  0.93  0.58  TC11  1.0  1.2  1.5  3.4  5.8  6.8  0.92  0.65  TC1-pPv  1.1  1.3  1.6  4.6  5.2  7.1  0.87  0.54  TC2-pPv  0.4  0.8  1.5  1.5  3.9  6.8  0.72  −0.29  TC4-pPv  1.1  1.2  1.5  4.8  5.6  8.6  0.89  −0.50  TC6-pPv  2.3  2.4  2.9  9.1  11.5  14.9  0.70  0.89  TC7-pPv  0.6  0.9  1.2  2.2  4.3  6.0  0.91  0.05  Observed  KDT12        ≤5.0          T10          5.0        MD87            12.0      DH89            8.0      SV03            4.0      KPF03            3.5      GS00              0.95    View Large Constraints on peak-to-peak amplitude may further be inferred from seismic phases reflecting either on the core side or on the mantle side. Using PKKP and PcP phases, Garcia & Souriau (2000) showed that 95 per cent of CMB topography with wavelength larger than 300 km should be between −4.0 and 4.0 km. Our thermal and thermochemical models are in good agreement with this observation (Figs 6 and 10; Table 3). The fraction of CMB area with topography ranging between −4.0 and 4.0 km calculated from unfiltered models is higher than 85 per cent except, again, for of thermochemical models with very large density contrast (Table 3). Importantly, the observations of Garcia & Souriau (2000) do not rule out the presence of deep troughs associated with downwellings since, according to our models, such features would be very localized and would cover a small fraction of the CMB area. Overall, the peak-to-peak amplitude of long-wavelength dynamic topography is unlikely to provide a clear diagnosis on the nature, thermal or thermochemical, of the deep mantle. On average, thermochemical models with density contrast in the range 100–200 kg m−3 have peak-to-peak amplitude lower than other models, but peak-to-peak amplitude further depends on other parameters, in particular to the sensitivity of viscosity to temperature. 4.2 Pattern: correlation between long-wavelength topography and shear wave velocity Lassak et al. (2010) showed that CMB dynamic topography induced by purely thermal and thermochemical models of convection strongly differs at short wavelengths, but are not substantially different at long wavelengths, that is, that current seismic maps of dynamic topography, which are limited to spherical harmonic degree l = 4, would be unable to discriminate between purely thermal and thermochemical structures. They further pointed out that long-wavelength topography and shear wave velocity would be anticorrelated. This conclusion may however change, depending on details of the thermochemical model. Left and middle columns in Fig. 13 compares the long-wavelength (spherical harmonic filter l = 0–4) distributions of lowermost mantle shear wave velocity anomalies and CMB dynamic topography for models T1, TC1, TC4, and TC7 (Table 2). For each model, we calculated shear wave velocity anomalies from the thermal and compositional distributions they predict, using the approach and seismic sensitivities described in Deschamps et al. (2012). We then averaged out these anomalies according to the radial parametrization of Masters et al. (2000), in which the lowermost layer samples the depth range 2700–2891 km. Non-dimensional temperature is rescaled with the super-adiabatic temperature jump ΔTS = 2500 K, leading to rms temperature anomalies typically around −400 K and 550 K in plume and slab regions, respectively, and 450 K within thermochemical reservoirs. For thermochemical models, we assumed that the dense material is enriched in iron and bridgmanite by 3 per cent and 18 per cent compared to regular (pyrolitic) material. For purely thermal model T1, plumes cause slower than average seismic velocity and positive topography. As a result, long-wavelength seismic velocity anomalies and CMB dynamic topography are strongly anticorrelated (Figs 13a and b). In thermochemical models with small density contrast, typically ∼90 kg m−3 and less, as model TC4 (Bz = 0.15, corresponding to ΔρC = 93 kg m−3) and the models of Lassak et al. (2010), thermochemical piles induce shallow depressions (Figs 9a–d). However, when filtered to long-wavelengths, topography beneath piles is positive, and long-wavelengths shear wave velocity anomalies and CMB topography are, again, strongly anticorrelated (Figs 13g and h). CMB topography may thus be unable to discriminate between a purely thermal and a weakly (ΔρC < 100 kg m−3) thermochemical structure of the mantle. However, for larger density contrasts, in particular model TC1 (Bz = 0.23, corresponding to ΔρC = 142 kg m−3), reservoirs of primordial material induce slower than average shear wave velocity (because they are hotter than average and enriched in iron), and cause a depression in the CMB (Fig. 5). Shear wave velocity anomalies and CMB topography are then strongly correlated (Figs 13d and e). For model TC7, in which intrinsic viscosities of primordial and regular material are similar (ΔηC = 1), the depression caused by reservoirs of primordial material is less pronounced (Figs 11g–h), resulting in a poorer, but still significant correlation between the maps of shear wave velocity anomalies and CMB dynamic topography (Figs 13j and k). Our models are not expected to strictly reproduce observed seismic patterns above the CMB. It is however interesting to compare long-wavelength shear wave velocity anomalies predicted by models of convection (Fig. 13, left column) with observed tomography from model SB10L18 (Masters et al.2000) at z = 2790 km depth (Fig. 14a). For model T1, shear wave velocity anomalies averaged in the layer 2700–2891 km vary in the range –0.7 to 0.7 per cent with rms amplitude around 0.3 per cent. This is much smaller than SB10L18, which varies between –2.5 and 1.4 per cent, with rms amplitude of 0.9 per cent. Explaining amplitudes of shear wave velocity anomalies observed in SB10L18 would require much larger temperature anomalies, around ± 1200 K. In thermochemical models, by contrast, shear-velocity anomalies have rms amplitude comparable to that of SB10L18. For model TC1, for instance, they vary between –2.1 and 1.7 per cent with rms amplitude around 1.1 per cent. Shear-velocity anomalies patterns predicted by models TC1 and TC4 are dominated by spherical harmonic degrees 2 and 3, and their power spectra are overall in good agreement with that of SB10L18 (Fig. 14b). The pattern predicted by model TC7 is more complex, with a strong contribution of spherical harmonic degree 4 that is not seen in SB10L18. Figure 13. View largeDownload slide Shear wave velocity anomalies in the lowermost mantle (2700–2891 km) (left column), CMB dynamic topography (middle column) and CMB isostatic topography (right column) predicted by purely thermal model T1 (Fig. 4) and thermochemical models TC1 (Fig. 5), TC4 (Fig. 9a-d), and TC7 (Fig. 11e-h), filtered for spherical harmonic degrees up to l = 4. Model properties are listed in Table 2. Shear wave velocity anomalies are calculated from the thermal and compositional fields predicted by each model using the seismic sensitivities of Deschamps et al. (2012), and are averaged out in the depth range 2700–2891 km, corresponding to the lowermost layer of Masters et al. (2000) tomographic model. For thermochemical models, primordial material is assumed to be enriched in iron and bridgmanite by 3 per cent and 18 per cent, respectively, compared to regular (pyrolitic) material. Isostatic topography is calculated from density anomalies averaged in the lowermost 400 km of the mantle (2591–2891 km). These density anomalies are directly estimated from the thermal and compositional fields predicted by each model, without using seismic sensitivity. Relative density anomalies, not shown here, are fully correlated with isostatic topography, with scaling factor between isostatic topography (in km) and relative density anomalies (in per cent), Rh, indicated at the bottom right of each plot. Interval of contour levels are 0.5 per cent for dlnVS and 1 km for CMB topography. Root mean square (rms) of each distribution is indicated at the bottom left of each plot. Figure 13. View largeDownload slide Shear wave velocity anomalies in the lowermost mantle (2700–2891 km) (left column), CMB dynamic topography (middle column) and CMB isostatic topography (right column) predicted by purely thermal model T1 (Fig. 4) and thermochemical models TC1 (Fig. 5), TC4 (Fig. 9a-d), and TC7 (Fig. 11e-h), filtered for spherical harmonic degrees up to l = 4. Model properties are listed in Table 2. Shear wave velocity anomalies are calculated from the thermal and compositional fields predicted by each model using the seismic sensitivities of Deschamps et al. (2012), and are averaged out in the depth range 2700–2891 km, corresponding to the lowermost layer of Masters et al. (2000) tomographic model. For thermochemical models, primordial material is assumed to be enriched in iron and bridgmanite by 3 per cent and 18 per cent, respectively, compared to regular (pyrolitic) material. Isostatic topography is calculated from density anomalies averaged in the lowermost 400 km of the mantle (2591–2891 km). These density anomalies are directly estimated from the thermal and compositional fields predicted by each model, without using seismic sensitivity. Relative density anomalies, not shown here, are fully correlated with isostatic topography, with scaling factor between isostatic topography (in km) and relative density anomalies (in per cent), Rh, indicated at the bottom right of each plot. Interval of contour levels are 0.5 per cent for dlnVS and 1 km for CMB topography. Root mean square (rms) of each distribution is indicated at the bottom left of each plot. Figure 14. View largeDownload slide Long-wavelength shear wave velocity anomalies in the lowermost mantle. (a) Observed tomography (SB10L18, Masters et al. 2000) filtered for spherical harmonic degrees up to l = 4. Interval of contour levels is 0.5 per cent. (b) Power spectra of SB10L18 (histogram) and selected models shown in Fig. 13 (colour code) up to spherical harmonic degree l = 8. Figure 14. View largeDownload slide Long-wavelength shear wave velocity anomalies in the lowermost mantle. (a) Observed tomography (SB10L18, Masters et al. 2000) filtered for spherical harmonic degrees up to l = 4. Interval of contour levels is 0.5 per cent. (b) Power spectra of SB10L18 (histogram) and selected models shown in Fig. 13 (colour code) up to spherical harmonic degree l = 8. To quantitatively compare long-wavelength patterns of topography and shear wave velocity predicted by models of convection, we calculated the correlation χh-dlnVs between these distributions (Table 3 and Fig. 12b). Clearly, CMB topography and shear-velocity anomalies for purely thermal models and thermochemical models with small density contrast anticorrelate, with χh-dlnVs ≤ –0.5, while they correlate for most thermochemical models with moderate and strong density contrast, with χh-dlnVs ≥ 0.5. Interesting exceptions are cases with moderate thermal viscosity contrast (ΔηT = 106; TC2 and TC2-pPv), for which CMB topography and shear-velocity anomalies are de-correlated or slightly anticorrelated, and, again, models with no chemical viscosity contrast (ΔηC = 1; TC7), for which distributions of topography and shear-velocity are only slightly correlated. These results suggest that long-wavelength CMB dynamic topography potentially provides a test to discriminate between purely thermal models of lower mantle and thermochemical models in which LLSVPs consist of material substantially denser (ΔρC > 100 kg m−3, i.e. ∼1.8 per cent of PREM and more) and intrinsically more viscous than the surrounding mantle. Interestingly, Tanaka (2010) found depressions of ∼2 km deep beneath Africa and the Pacific, that is, fairly correlated with LLSVPs. If correct, and following our models, this would support the hypothesis that LLSVPs are thermochemical structures denser than the surrounding mantle, with density contrast typically in the range 100–200 kg m−3. This scenario is further supported by available estimates of chemical density anomalies (Trampert et al.2004; Mosca et al.2012), which have a peak-to-peak amplitude around 2.0–3.0 per cent. In addition, because primordial material may be enriched in bridgmanite (Trampert et al.2004; Mosca et al.2012), which is more viscous than ferro-periclase by up to 3 orders of magnitude (Yamazaki & Karato 2001), LLSVPs may be intrinsically more viscous than the surrounding mantle. Other models of CMB topography, however, do not observe clear correlation or anticorrelation between shear-velocity and CMB topography patterns. For instance, the topography observed by Doornbos & Hilton (1989) is dominated by spherical harmonic degree 1 with a depression beneath India and a topographic high beneath Western Europe and Africa, while the pattern obtained by Morelli & Dziewonski (1987) is more complex, with topographic highs beneath South Indian Ocean, Northern Atlantic and Northern Pacific. Soldati et al. (2012) performed inversion of P-wave data and reflected seismic phases (PcP and PKP) jointly for compressional velocity anomalies and CMB topography, using a geodynamically constrained method. More precisely, density anomalies were scaled with seismic velocity anomalies and combined with CMB topography kernels calculated using a spectral method (e.g. Forte & Peltier 1991). Resulting CMB topography maps have peak-to-peak amplitude around 8 km, and are fairly anticorrelated with compressional velocity maps above the CMB. In particular, LLSVPs are associated with topography highs with strong positive topography of up to about 5 km, beneath the central Pacific and southern tip of Africa. It has further been pointed out that CMB topography models inferred from geodynamically constrained tomography fit normal mode constraints better than models deduced directly from seismic tomography (Soldati et al.2013). It should be noted, however, that the approach used by Soldati et al. (2012) implicitly assumes that chemical anomalies are small, and that they are, at a given depth, fully correlated with seismic velocity anomalies. Furthermore it does not account for lateral variations in viscosity (e.g. related to variations in temperature), which, as illustrated by our calculations and those of Lassak et al. (2010), may substantially affect the amplitude of topography. 5 DISCUSSION Our calculations show strong differences between the CMB dynamic topographies induced by purely thermal and thermochemical models of convection. In purely thermal models plumes elevate the CMB, while hot thermochemical piles depress it. In both cases, increasing the thermal viscosity contrast decreases the amplitude of CMB topography. In thermochemical models, the buoyancy ratio, measuring the density contrast between primordial and regular material, and the thermal and compositional viscosity contrasts further affect the shape and stability of reservoirs of primordial material, and the amplitude of dynamic topography, as sketched in plots c and d of Fig. 1. Increasing the chemical density contrast and viscosity ratio enhances the depression beneath reservoirs. 5.1 Influence of post-perovskite Our calculations further indicate that the presence of pPv at the bottom of the mantle does not substantially influence the dynamic topography. All other parameters being the same, the presence of pPv slightly modifies the amplitude of topography, but the effects due to varying one specific parameter are unchanged. This conclusion assumes that the viscosity of pPv is similar or close to that of bridgmanite. If, in contrast, and as suggested by ab initio calculations (Ammann et al.2010), pPv is less viscous than bridgmanite by two to three orders of magnitude, dynamic topography may be strongly altered. Both the numerical models of Yoshida (2008) and the CMB maps of Soldati et al. (2012) report a strong attenuation of CMB topography when a low viscosity layer is imposed in the lowermost 200 km of the mantle completely surrounding the core. In our models, the pPv stability region does not cover the entire core surface, but forms discontinuous patches located outside plume clusters and, in the thermochemical cases, outside reservoirs of primordial material. In these regions, and if pPv is much less viscous than bridgmanite, CMB topography may thus flatten out compared to that found in our models. 5.2 Comparison with previous models Our models are in overall agreement with those of Lassak et al. (2010). In particular, the influences of buoyancy ratio and thermal viscosity contrast on CMB topography they observed are similar to what we see in our simulations. Note, however, that due to different choices of the super-adiabatic temperature jump and reference thermal expansion (Table 1), and of a different definition of the buoyancy ratio (eq. 1), the chemical density contrast at the bottom of the system for a given value of buoyancy ratio is much larger in our models. For instance, in the models of Lassak et al. (2010) with B = 0.6 it is around 60 kg m−3 (1 per cent of PREM), whereas in our models with Bz = 0.23 it reaches 142 kg m−3 (2.5 per cent of PREM). As a result, in model TC7 thermochemical reservoirs induce a 0.7 km deep depression in dynamic topography (Figs 8e–h), and a plateau with low elevation in Lassak et al. (2010) models. In model TC1 (Fig. 5), this trend is amplified by the fact that primordial material is intrinsically more viscous than the regular material by a factor 30, which deepens the depression beneath reservoirs of primordial material. To investigate the influence of LLSVPs on dynamic topography, Yoshida (2008) ran models of convection in which density anomalies are scaled from seismic tomography model SMEAN (Becker & Boschi 2002), and arrived at similar conclusions regarding the effect of the chemical density contrast. In purely thermal models, LLSVPs result in positive topography of ∼5 km. Models assuming piles denser by 1 per cent lead to a relatively flat positive topography beneath LLSVPs, while a larger density contrast (2 per cent) causes relatively flat topography with local depressions ∼2 km deep, not unlike our model TC7. Interestingly, this negative topography flattens out when the viscosity of the lowermost mantle is decreased by one or two orders of magnitude. Note, however, that Yoshida (2008) did not run models in which LLSVP viscosity is intrinsically different from that of surrounding material. 5.3 Predicted isostatic topography A recent analysis of Stoneley mode splitting data, which are sensitive to the lowermost mantle structure, indicates that LLSVPs might be lighter than surrounding mantle (Koelemeijer et al.2017), in contradiction with earlier normal mode studies (Ishii & Tromp 1999; Trampert et al.2004; Mosca et al.2012). Stoneley mode are further sensitive to CMB topography, and density estimates thus trade-off with this topography. Assuming that long-wavelength (>1000 km) CMB topography is due to isostatic compensation, Koelemeijer et al. (2017) found that, for degree 2, their data set is best fitted if LLSVPs are lighter than PREM by about 1.4 per cent and CMB is elevated by about 1 km. For degrees 4 and 6, best fits are, by contrast, obtained for LLSVPs denser than average mantle by about 2.5 per cent and a CMB depressed by 1–5 km. When all structural degrees are taken into account, light LLSVPs with CMB topography around 1 km or more are again preferred, but the overall fit is much reduced. In addition, when accounting for seismic normal modes used in previous studies, density models do not have preference for denser or lighter LLSVPs. Finally, it should be kept in mind that splitting data are based on the self-coupling approximation (Woodhouse & Girnius 1982), which might bias the inference of density structure (Al-Attar et al.2012; Yang & Tromp 2015). Lateral variations in temperature and composition predicted by our models can be used to estimate lateral density anomalies and reconstruct CMB isostatic topography. For purely thermal models, one expects that CMB dynamic and isostatic topographies have similar trends. On the contrary, because thermochemical piles and reservoirs are hotter than average, patterns of dynamic and isostatic topographies may differ substantially, in particular for low buoyancy ratio. Right column in Fig. 13 plots CMB isostatic topography estimated from density anomalies averaged in the lowermost 400 km of the mantle for models T1, TC1, TC4 and TC7. Density anomalies, not shown here, are, by definition, fully correlated with isostatic topography, with scaling factor between topography (in km) and relative density anomalies (in per cent) around −0.27 (the exact value of this scaling factor depends on the average density in the lowermost 400 km layer, which is very similar for all 4 models shown here). For the calculation of density anomalies, temperature anomalies are rescaled with the super adiabatic temperature jump, and thermal expansion and density are rescaled with their surface values (Table 1). Reference density and thermal expansion vary as indicated in Fig. 2, and density anomalies are calculated with respect to horizontal average (which is different and usually higher than the reference density). For comparison with the dynamic topography shown in the middle column of Fig. 13, we filtered density anomalies and isostatic topography for spherical harmonic degrees l = 0–4. For purely thermal model T1, isostatic topography has a pattern similar to that of dynamic topography, but its amplitude is smaller by about a factor 4. Within plume regions, the average filtered density anomaly and isostatic topography are −0.1 per cent and 0.4 km (Fig. 13c), respectively, and are thus consistent with Stoneley mode splitting measurements of Koelmeijer et al. (2017). In thermochemical models, thermal effects alter the effective density within reservoirs and piles of dense material. For small density contrast, as in model TC4 (Fig. 13i), thermal effects nearly compensate chemical effects within piles, leading to average density anomaly and isostatic topography of about 0.3 per cent and −1.0 km, respectively. For larger density contrast, as in models TC1 and TC7, density anomalies within piles are still dominated by chemical effects, resulting in relatively large isostatic topographies (Figs 13f and l). For instance, in model TC1, density anomalies within thermochemical reservoirs and isostatic topography beneath them are equal to 0.8 per cent and −3.2 km on average. Such reservoirs provide poorer fit to Koelemeijer et al. (2017) observations than light elevated plume clusters in purely thermal models, but are still in reasonable agreement with these data (see their fig. 3 for RLL = −0.6 and H = −5, corresponding to LLSVP total density anomaly of 0.8 per cent and isostatic topography of −2.5 km). Isostatic topography at CMB may be considered as a first order estimate of CMB topography, taking into account only the effects of hydrostatic pressure. Dynamic topography further accounts for stresses related to mantle flow, in particular slab pushing or plume pulling. This appears very clearly when comparing middle and right columns in Fig. 13. For instance, in case T1, hot, lighter plumes induce positive isostatic topography. They further pull the CMB upwards, thus increasing CMB topography. As a result, compared to isostatic topography (Fig. 13c), dynamic topography has a similar pattern, but more amplitude (Fig. 13b). For thermochemical cases, dynamical and isostatic effects are opposing each other. Dynamic and isostatic topographies have slightly different patterns, and depressions beneath thermochemical reservoirs are shallower in dynamic topography (compare, for instance, Figs 13e and f for case TC1). Note that for small chemical density contrasts, as in case TC4, dynamic effects lead to very shallow depressions (Fig. 9c), and when filtered for long-wavelengths, dynamic topography beneath thermochemical piles is positive (Fig. 13h). Because it accounts for dynamics effects related to mantle flow, we expect dynamic topography to provide a better description of CMB topography than isostatic topography does. Interestingly, in dynamic topography predicted by case TC1 and filtered with l = 0–4, the depth of the depressions beneath thermochemical reservoirs is on average −1.0 km (Fig. 13e), which still reasonably fits Stoneley mode observations. 5.4 CMB excess ellipticity An important constraint on dynamic topography is provided by the excess ellipticity of the CMB, that is, its spherical harmonic degree l = 2 and order m = 0 term (Y20). This excess can be deduced from Very Long Baseline Interferometry (VLBI) measurements of the retrograde annual nutation of the Earth's rotation axis. More precisely, the period of the free core nutation (FCN) depends on the CMB dynamic flattening, ef, and the resonance between FCN and retrograde annual nutation increases the amplitude of this annual nutation. The FCN period deduced from VLBI data, first measured by Gwinn et al. (1986), indicates that the CMB dynamical flattening is larger than the value expected for a hydrostatic Earth, with an excess dynamical flattening of 5 per cent. Assuming that this flattening is due to a homogeneous non-hydrostatic layer at the top of the outer core, it can be interpreted as an excess difference between the equatorial and polar radii of the CMB of 490 ± 110 m. Taking into account more complexities, including ocean tides, electromagnetic coupling between core and mantle, and anelastic effects in the mantle, Mathews et al. (2002) found a slightly smaller excess, of ∼390 m, corresponding to an excess CMB dynamical flattening of 3.8 per cent. This excess of CMB ellipticity is usually attributed to mantle dynamics (Gwinn et al.1986; Forte et al.1995; Lumb et al.1995). Our models do not define a polar axis and associated equatorial plane that would allow the calculation of excess CMB ellipticity. Nonetheless, one may point out that excess CMB ellipticity would be difficult to reconcile with large, deep depressions in CMB dynamic topography, as those predicted in model TC1 (Bz = 0.23 and ΔηC = 30), located along the equator. This problem may be solved in models with a slightly lower buoyancy ratio, around Bz = 0.2 (ΔρC = 124 kg m−3), which is still consistent with density anomalies observed in the lowermost mantle (Trampert et al.2004; Mosca et al.2012), and/or lower chemical density contrast (typically, ΔηC around 10 and less). In these models, the depressions associated with reservoirs of primordial material are much attenuated, as illustrated by model TC7, and would therefore have a limited influence on excess CMB ellipticity. Interestingly, such models still predict long-wavelength topographies that are distinct from those obtained for purely thermal convection (Fig. 13), with, for models with ΔηC around 10 a correlation with shear wave velocity anomalies that is still important (around 0.6 for model TC8). In addition, the link between dynamic topography and observed excess ellipticity may be affected by other factors. First, it should be kept in mind that the period of the FCN constrains the dynamic flattening, and that interpreting this flattening in terms of geometric flattening requires the prescription of a non-hydrostatic layer at the top of the outer core, with thickness and density as free parameters. Moreover, it has been pointed out that other components of CMB topography than its Y20 term may influence the period of the FCN (Wu & Wahr 1997), and therefore the VLBI measurements and the estimates of the excess CMB ellipticity inferred from them. Finally, gravitational and viscous couplings between specific mantle structures (e.g. LLSVPs) and the core may also affect the period of the FCN. For instance, a recent study (Chao 2017) showed that the 6 yr periodic variations of the length of the day are well explained by gravitational coupling between the inner core and the lowermost mantle. Interestingly, this coupling implies mass excess in regions overlapping LLSVPs (Ding & Chao, submitted), thus contradicting the Stoneley mode analysis of Koelemeijer et al. (2017), but supporting earlier normal mode studies (Ishii & Tromp 1999; Trampert et al.2004; Mosca et al.2012). The effect of this coupling on the period of the FCN, however, remains to be investigated. 5.5 Implications for core dynamics Another interesting question that remains to be investigated is the effect of dynamic topography on outer core flow. Calkins et al. (2012) run core flow models in cylindrical geometry with ridges (i.e. trenches on the mantle side) of Gaussian shape, and found that this topography results in stationary Rossby waves, which, in turn, increase the strength of the zonal flow and induce flow structures that remain fixed relative to the mantle. More recently, Tarduno et al. (2015) suggested that a step-like topography in CMB due to the African LLSVP could induce upwellings at the top of the core beneath the south of Africa, resulting in reversed magnetic flux at the CMB and persistent low magnetic field at the surface of the Earth. The gradient topography predicted by our models are rather modest, as illustrated in Fig. 15 (note that the scale of vertical axis is exaggerated by about a factor 115 compared to that of the horizontal axis). For instance, the positive topography induced by plume cluster in model T1 gently increases up to its culminating point (Fig. 15a). The depressions induced by thermochemical reservoirs in model TC1 are flat with marked slopes on their edges, but the gradient of these slopes are typically around 2–4 m km−1 (Fig. 15c). In models with small chemical density contrast, the depressions induced by piles are attenuated and much shallower than the negative topography beneath slabs (Fig. 15d). Slabs induce deeper and steeper depressions on more localised areas, but again, topographic gradients remain limited, up to 15 m km−1 (Fig. 15b). Thus, if CMB topography of real Earth is well described by the dynamic topography predicted by our models, it may have a limited influence on the outer core flow. Figure 15. View largeDownload slide Topographic profiles along different great circles for (a,b) model T1 (Fig. 4), (c) model TC1 (Fig. 5), and (d) model TC4 (Figs 9a–d). Note that the scale of the vertical axis is exaggerated by about a factor 115 compared to that of the horizontal axis, that is, profiles appear much steeper than in real data. Figure 15. View largeDownload slide Topographic profiles along different great circles for (a,b) model T1 (Fig. 4), (c) model TC1 (Fig. 5), and (d) model TC4 (Figs 9a–d). Note that the scale of the vertical axis is exaggerated by about a factor 115 compared to that of the horizontal axis, that is, profiles appear much steeper than in real data. 6 CONCLUSIONS We calculated the dynamic topography at CMB induced by different models of mantle convection, and investigated the influence of important input parameters of these models. In all models, cold downwellings, modelling slabs, induce very deep, but local depressions. In purely thermal models, plumes or plume clusters induce positive topography that decreases in amplitude as the thermal viscosity ratio increases. For moderate buoyancy ratios, corresponding to compositional density contrast (ΔρC) in the range 100–200 kg m−3, consistent with estimates of chemical density anomalies in the lowermost mantle (Trampert et al.2004; Mosca et al.2012), large, stable reservoirs of primordial material are generated. In their interiors, these reservoirs cause depressions in the CMB topography, surrounded by a ridge at their edges. The depth of the depressions and the height of the ridges increase with increasing ΔρC and chemical viscosity contrast (ΔηC). The presence of the pPv phase at the bottom of the mantle does not alter these conclusions, provided that the viscosity of pPv is similar to that of bridgmanite. If, in contrast, pPv viscosity is lower by 2 orders of magnitude or more (Amman et al.2010), the CMB topography may substantially flatten out, as suggested by models in which a weak (low viscosity) layer is imposed at the bottom of the mantle (Yoshida 2008; Soldati et al.2012). Additional calculations are, however, needed to determine the influence of discontinuous patches of weak pPv, as those obtained in models T1-pPv (Figs 7c and d) and TC1-pPv (Figs 8e and f). The peak-to-peak amplitude of CMB dynamic topography predicted by our models for different spherical harmonic filters is comparable to that estimated by available seismic constraints on CMB topography. Purely thermal models usually have larger amplitude than thermochemical models with similar properties, and for spherical harmonic filters including degree l = 2 and degrees l = 2 and l = 4, the amplitude predicted by purely thermal models is overall slightly larger than the maximum amplitude estimated from seismic data. Because it is sensitive to other model parameters, peak-to-peak amplitude alone is however unlikely to discriminate between purely thermal and thermochemical structures of lower mantle. On another hand, a comparison between long-wavelength patterns of shear wave velocity anomalies and topography may provide important clues. The positive topography induced by plumes, when filtered to long wavelengths, anticorrelates with shear wave velocity anomalies induced by the thermal field. By contrast, for thermochemical models, long-wavelength CMB topography correlates well with seismic velocity anomalies for models with ΔρC around 140 kg m−3 or more and intrinsically more viscous primordial material (ΔηC > 1). Shear wave velocity and CMB topography maps progressively decorrelate with decreasing ΔρC or decreasing ΔηC, and for ΔρC around 90 kg m−3 and less, CMB topography and shear wave velocity anticorrelates, as for purely thermal cases. Long-wavelength pattern of CMB topography, provided they can be measured accurately, therefore provides a potential constraint to discriminate between different possible thermochemical structures of lowermost mantle. Acknowledgements We are grateful to Nicolas Flament and to the editor (Gaël Choblet) for their constructive and careful reviews that helped improving a first version of this article. We are also grateful to Hagay Amit and Andrew Valentine for additional discussions on core–mantle coupling and density mapping from normal mode data. This research was funded by the Ministry of Science and Technology (MOST) of Taiwan grant 105-2116-M-001-017, Academia Sinica (Taipei) grant AS-102-CDA-M02, and Orchid grant 103PMRIA0100024. All data and models used in this study are available upon request to FD (frederic@earth.sinica.edu.tw). REFERENCES Al-Attar D., Woodhouse J.H., Deuss A., 2012. Calculation of normal mode spectra in laterally heterogeneous earth models using an iterative direct solution method, Geophys. J. Int. , 189, 1038– 1046. Google Scholar CrossRef Search ADS   Amman M.W., Brodholt J.P., Wookey J., Dobson D.P., 2010. First principles constraints on diffusion in lower-mantle minerals and a weak D″ layer, Nature , 465, 462– 465. Google Scholar CrossRef Search ADS PubMed  Becker T.W., Boschi L., 2002. A comparison of tomographic and geodynamic mantle models, Geochem. Geophys. Geosyst. , 3, doi:10.1029/2001GC000168. Calkins M.A., Noir J., Eldredge J.D., Aurnou J.M., 2012. The effects of boundary topography on convection in Earth's core, Geophys. J. Int. , 189, 799– 814. Google Scholar CrossRef Search ADS   Chao B.F., 2017. Dynamics of axial torsional libration under the mantle-inner core gravitational interaction, J. geophys. Res. , 121, doi:10.1002/2016JB013515. Christensen U.R., Yuen D.A., 1985. Layered convection induced by phase transitions, J. geophys. Res. , 90( 10), 291– 210, 300. Crameri F., Tackley P.J., 2014. Spontaneous development of arcuate single-sided subduction in global 3-D mantle convection models with a free surface, J. geophys. Res. , 119, 5921– 5942. Google Scholar CrossRef Search ADS   Crameri F.et al.  , 2012. comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the ‘sticky air’ method, Geophys. J. Int. , 189, 38– 54. Google Scholar CrossRef Search ADS   Davies D.R., Goes S., Davies J.H., Schuberth B.S.A., Bunge H.-P., Ritsema J., 2012. Reconciling dynamic and seismic models of Earth's lower mantle: The dominant role of thermal heterogeneity, Earth planet. Sci. Lett. , 353–354, 253– 269. Google Scholar CrossRef Search ADS   Davies D.R., Goes S., Lau H.C.P., 2015. Thermally dominated deep mantle LLSVPs: a review, in The Earth's Heterogeneous Mantle , pp. 441– 478, eds Khan A., Deschamps F., Springer. Google Scholar CrossRef Search ADS   Davies J.H., 2013. Global map of solid Earth heat flow, Geochem. Geophys. Geosyst. , 14, doi:10.1029/ggge.20271. Deschamps F., Tackley P.J., 2009. Searching for models of thermo-chemical convection that explain probabilistic tomography II - Influence of physical and compositional parameters, Phys. Earth planet. Inter. , 176, 1– 18. Google Scholar CrossRef Search ADS   Deschamps F., Cobden L., Tackley P.J., 2012. The primitive nature of large low shear-wave velocity provinces, Earth planet. Sci. Lett. , 349–350, 198– 208. Google Scholar CrossRef Search ADS   Deschamps F., Li Y., Tackley P.J., 2015. Large-scale thermo-chemical structure of the deep mantle: observations and models, in The Earth's Heterogeneous Mantle , pp. 479– 515, eds Khan A., Deschamps F., Springer. Google Scholar CrossRef Search ADS   Deschamps F., Khan A., 2016. Electrical conductivity as a constraint on lower mantle thermo-chemical structure, Earth planet. Sci. Lett. , 450, 108– 119. Google Scholar CrossRef Search ADS   Ding H., Chao B.-F., A 6-year westward rotary motion in the Earth: detection and possible MICG coupling mechanism, submitted to Earth Planet. Sci. Lett.  Doornbos D.J., Hilton T., 1989. Models of the core-mantle boundary and the travel times of internally reflected core phase, J. geophys. Res. , 94, 15 741–15 751. Google Scholar CrossRef Search ADS   Dziewonski A.M., Anderson D.L., 1981. Preliminary reference Earth model, Phys. Earth Planet. Inter. , 25, 297– 356. Google Scholar CrossRef Search ADS   Forte A.M., Peltier R., 1991. Viscous flow models of global geophysical observables 1. Forward problems, J. geophys. Res. , 96, 20 131–20 159. Google Scholar CrossRef Search ADS   Forte A.M., Mitrovica J.X., Woodward R.L., 1995. Seismic-geodynamic determination of the origin of excess ellipticity of the core–mantle boundary, Geophys. Res. Lett. , 22, 1013– 1016. Google Scholar CrossRef Search ADS   Garcia R., Souriau A., 2000. Amplitude of the core–mantle boundary topography estimated by stochastic analysis of core phases, Phys. Erth planet. Inter. , 117, 345– 359. Google Scholar CrossRef Search ADS   Garnero E.J., McNamara A., Shim S.-H., 2016. Continent-sized anomalous zones with low seismic velocity at the base of Earth's mantle, Nat. Geosci. , 9, 481– 489. Google Scholar CrossRef Search ADS   Gwinn C.R., Herring T.A., Shapiro I.I., 1986. Geodesy by radio interferometry: studies of the forced nutations of the Earth 2. Interpretation, J. geophys. Res. , 91, 4755– 4765. Google Scholar CrossRef Search ADS   Hager B.H., Richards M.A., 1989. Long wavelength variations in Earth's geoid: physical models and dynamical implications, Phil. Trans. R. Soc. A , 328, 309– 327. Google Scholar CrossRef Search ADS   Hansen U., Yuen D.A., 1989. Dynamical influences from thermo-chemical instabilities at the core–mantle boundary, Geophys. Res. Lett. , 16, 629– 632. Google Scholar CrossRef Search ADS   Hernlund J., Houser C., 2008. On the statistical distribution of seismic velocities in Earth's deep mantle, Earth planet. Sci. Lett. , 265, 423– 437. Google Scholar CrossRef Search ADS   Ishii M., Tromp J., 1999. Normal-mode and free-air gravity constraints on lateral variations in velocity and density of Earth's mantle, Science , 285, 1231– 1236. Google Scholar CrossRef Search ADS PubMed  Jaupart C., Mareschal J.-C., 2007. Heat flow and thermal structure of the lithosphere, in Treatise on Geophysics , 6, pp 218– 251, ed. Schubert, G., Elsevier. Jaupart C., Labrosse S., Lucazeau F., Mareschal J.-C., 2015. Temperature, heat, and energy in the mantle of the Earth, in Treatise on Geophysics , 2nd edn, 7, pp. 218– 251, ed. Schubert G., Elsevier. Kageyama A., Sato T., 2004. “Yin-Yang grid”: an overset grid in spherical geometry, Geochem. Geophys. Geosys. , 5, doi:10.1029/2004GC000734. Koelemeijer P.J., Deuss A., Trampert J., 2012. Normal mode sensitivity to Earth's D″ layer and topography on the core–mantle boundary: what we can and cannot see, Geophys. J. Int. , 190, 553– 568. Google Scholar CrossRef Search ADS   Koelemeijer P., Deuss A., Ritsema J., 2017. Density structure of Earth's lowermost mantle from Stoneley mode splitting observations, Nat. Commun. , 8, doi:10.1038/ncomms15241. Konishi K., Fuji N., Deschamps F., 2017. Elastic and anelastic structure of the lowermost mantle beneath the western Pacific from waveform inversion, Geophys. J. Int. , 208, 1290– 1304. Google Scholar CrossRef Search ADS   Korenaga J., Karato S.-I., 2008. A new analysis of experimental data on olivine rheology, J. geophys. Res. , 113, B02403, doi:10.1029/2007JB005100. Google Scholar CrossRef Search ADS   Koper K.D., Pyle M.L., Franks J.M., 2003. Constraints on aspherical core structure from PKiKP-PcP differential travel times, J. geophys. Res. , 108, doi:10.1029/2002JB001995. Lassak T.M., McNamara A.K., Zhong S., 2007. Influence of thermo-chemical piles on topography at Earth's core-mantle boundary, Earth planet. Sci. Lett. , 261, 443– 455. Google Scholar CrossRef Search ADS   Lassak T.M., McNamara A.K., Garnero E.J., Zhong S., 2010. Core–mantle boundary topography as a possible constraint on lower mantle chemistry and dynamics, Earth planet. Sci. Lett. , 289, 232– 241. Google Scholar CrossRef Search ADS   Lekić V., Cottaar S., Dziewonski A.M., Romanowicz B., 2012. Cluster analysis of global lower mantle tomography: A new class of structure and implications for chemical heterogeneity, Earth planet. Sci. Lett. , 357–358, 68– 77. Google Scholar CrossRef Search ADS   Li P., Karato S.-I., Wang Z., 1996. High-temperature creep in fine-grained polycrystalline CaTiO3, an analogue material of (Mg,Fe)SiO3 perovskite, Phys. Earth planet. Inter. , 95, 19– 36. Google Scholar CrossRef Search ADS   Li M., McNamara A.K., Garnero E.J., 2014a. Chemical complexity of hotspots caused by cycling oceanic crust through mantle reservoirs, Nat. Geosci. , 7, 366– 370. Google Scholar CrossRef Search ADS   Li Y., Deschamps F., Tackley P.J., 2014b. The stability and structure of primordial reservoirs in the lower mantle: insights from models of thermo-chemical convection in 3-D spherical geometry, Geophys. J. Int. , 199, 914– 930. Google Scholar CrossRef Search ADS   Li Y., Deschamps F., Tackley P.J., 2015. Effects of the post-perovskite phase transition properties on the stability and structure of primordial reservoirs in the lower mantle of the Earth, Earth planet. Sci. Lett. , 432, 1– 12. Google Scholar CrossRef Search ADS   Lumb L.I., Jarvis G.T., Aldridge K.D., DeLandro-Clarke W., 1995. The period of the free core nutation: towards a dynamical basis for an “extra-fattening” of the core-mantle boundary, Phys. Earth planet. Inter. , 90, 255– 271. Google Scholar CrossRef Search ADS   Masters G., Laske G., Bolton H., Dziewonski A.M., 2000. The relative behavior of shear velocity, bulk sound speed, and compressional velocity in the mantle: implication for thermal and chemical structure, in Earth's Deep Interior: Mineral Physics and Tomography from the Atomic to the Global Scale , pp. 63– 87, Geophysical Monograph Ser. 117, eds Karato,, S.-I. et al., American Geophysical Union. Google Scholar CrossRef Search ADS   Mathews P.M., Herring T.A., Buffett B., 2002. Modeling of nutation and precession: new nutation series for nonrigid Earth and insights into the Earth's interior, J. geophys. Res. , 107, doi:10.1029/2001JB000390. Morelli A., Dziewonski A.M., 1987. Topography of the core-mantle boundary and lateral homogeneity of the liquid core, Nature , 325, 678– 683. Google Scholar CrossRef Search ADS   Mosca I., Cobden L., Deuss A., Ritsema J., Trampert J., 2012. Seismic and mineralogical structures of the lower mantle from probabilistic tomography, J. geophys. Res. , 117, doi:10.1029/2011JB008851. Murakami M., Hirose K., Kawamura K., Sata N., Ohishi Y., 2004. Post-perovskite phase transition in MgSiO3, Science , 304, 855– 858. Google Scholar CrossRef Search ADS PubMed  Ricard Y., Fleitout L., Froideveaux C., 1984. Geoid heights and lithospheric stresses for a dynamic Earth, Ann. Geophys. , 2, 267– 286. Ritsema J., Deuss A., vanHeijst H.-J., Woodhouse J.H., 2011. S40RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements, Geophys. J. Int. , 184, 1223– 1236. Google Scholar CrossRef Search ADS   Schuberth B.S.A., Zaroli C., Nolet G., 2012. Synthetic seismograms for a synthetic Earth: long-period P- and S-wave traveltime variations can be explained by temperature alone, Geophys. J. Int. , 188, 1393– 1412. Google Scholar CrossRef Search ADS   Soldati G., Boschi L., Forte A.M., 2012. Tomography of core–mantle boundary and lowermost mantle coupled by geodynamics, Geophys. J. Int. , 189, 730– 746. Google Scholar CrossRef Search ADS   Soldati G., Koelemeijer P., Boschi L., Deuss A., 2013. Constraints on core–mantle boundary topography from normal mode splitting, Geochem. Geophys. Geosys. , 14, doi:10.1002/ggge.20115. Solomatov V.S., Reese C.C., 2008. Grain size variations in the Earth's mantle and the evolution of primordial chemical heterogeneities, J. geophys. Res. , 113, B07408, doi:10.1029/2007JB005319. Google Scholar CrossRef Search ADS   Sze E.K.M., van der Hilst R.D., 2003. Core mantle boundary topography from short period PcP, PKP, and PKKP data, Phys. Earth planet. Inter. , 135, 27– 46. Google Scholar CrossRef Search ADS   Tackley P.J., 1998. Three-dimensional simulations of mantle convection with a thermo-chemical CMB boundary layer: D″?, in The Core-Mantle Boundary Region , Geodynamical Ser., 28, pp. 231– 253, eds Gurnis M. et al., American Geophysical Union. Tackley P.J., 2008. Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid, Phys. Earth planet. Inter. , 171, 7– 18. Google Scholar CrossRef Search ADS   Tackley P.J., 2012. Dynamics and evolution of the deep mantle resulting from thermal, chemical, phase and melting effects, Earth-Sci. Rev. , 110, 1– 25. Google Scholar CrossRef Search ADS   Tackley P.J., King S.D., 2003. Testing the tracer ratio method for modeling active compositional fields in mantle convection simulations, Geochem. Geophys. Geosys. , 4, doi:10.1029/2001GC000214. Tanaka S., 2010. Constraints on the core‐mantle boundary topography from P4KP‐PcP differential travel times, J. geophys. Res. , 115, B04310, doi:10.1029/2009JB006563. Tarduno J.A., Watkeys M.K., Huffman T.N., Cottrell R.D., Blackman E.G., Wendt A., Scribner C.A., Wagner C.L., 2015. Antiquity of the South Atlantic Anomaly and evidence for top-down control on the geodynamo, Nat. Commun ., 6, doi:10.1038/ncomms8865. Tateno S., Hirose K., Sata N., Ohishi Y., 2009. Determination of post-perovskite phase transition boundary up to 4400 K and implications for thermal structure in D″ layer, Earth planet. Sci. Lett. , 277, 130– 136. Google Scholar CrossRef Search ADS   Trampert J., Deschamps F., Resovsky J.S., Yuen D.A., 2004. Probabilistic tomography maps significant chemical heterogeneities in the lower mantle, Science , 306, 853– 856. Google Scholar CrossRef Search ADS PubMed  Woodhouse J.H., Girnius T.P., 1982. Surface waves and free oscillations in a regionalized Earth model, Geophys. J. Int. , 68, 653– 673. Google Scholar CrossRef Search ADS   Wu X., Wahr J.M., 1997. Effects of non-hydrostatic core-mantle boundary topography and core dynamics on Earth rotation, Geophys. J. Int. , 128, 18– 42. Google Scholar CrossRef Search ADS   Yamazaki D., Karato S.-I., 2001. Some mineral physics constraints on the rheology and geothermal structure of Earth's lower mantle, Am. Mineral. , 86, 385– 391. Google Scholar CrossRef Search ADS   Yang H.-Y., Tromp J., 2015. Synthetic free-oscillation spectra: an appraisal of various mode-coupling methods, Geophys. J. Int. , 203, 1179– 1192. Google Scholar CrossRef Search ADS   Yoshida M., 2008. Core-mantle boundary topography estimated from numerical simulations of instantaneous mantle flow, Geochem. Geophys. Geosys. , 9, doi:10.1029/2008GC002008. Zhang S., Christensen U., 1993. Some effects of lateral viscosity varuiations on geoid and surface velocities induced by density anomalies in the mantle, Geophys. J. Int. , 114, 531– 547. Google Scholar CrossRef Search ADS   © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Geophysical Journal International Oxford University Press

# Constraints on core–mantle boundary topography from models of thermal and thermochemical convection

Geophysical Journal International, Volume 212 (1) – Jan 1, 2018
25 pages

/lp/ou_press/constraints-on-core-mantle-boundary-topography-from-models-of-thermal-SmcEah00Am
Publisher
Oxford University Press
ISSN
0956-540X
eISSN
1365-246X
D.O.I.
10.1093/gji/ggx402
Publisher site
See Article on Publisher Site

### Abstract

Abstract Mantle flow induces dynamic topography at the core–mantle boundary (CMB), with distribution and amplitude that depend on details of the flow. To assess whether observations of CMB topography can give constraints on deep mantle structure, we determine CMB dynamic topography associated with different models of mantle convection, including thermochemical and purely thermal models. We investigate the influence of key controlling parameters, specifically the thermal viscosity ratio (ΔηT) and, for thermochemical models, the density contrast (ΔρC) and viscosity ratio (ΔηC) between primordial and regular materials. In purely thermal models, plume clusters induce positive topography with an amplitude that decreases with increasing ΔηT. In thermochemical models with moderate density contrasts, around 100–200 kg m−3, reservoirs of dense material induce depressions in CMB topography, surrounded by a ridge of positive topography. The average depression depth and ridge height increase with increasing ΔρC and ΔηC, but decrease with increasing ΔηT. We find that for purely thermal models or thermochemical models with ΔρC ∼ 90 kg m−3 and less, the long-wavelength (spherical harmonic degrees up to l = 4) dynamic topography and shear wave velocity anomalies predicted by thermochemical distributions anticorrelate. By contrast, for models with ΔρC ≥ 100 kg m−3 and ΔηC > 1, long-wavelength dynamic topography and shear wave velocity anomalies correlate well. This potentially provides a test to infer the nature, that is, either purely or mostly thermal (ΔρC ≤ 100 kg m−3 m−3) or strongly thermochemical (ΔρC ≥ 100 kg m−3), of the low shear wave velocity provinces observed by global tomographic images. The presence of post-perovskite, provided that its viscosity is similar to that of bridgmanite, does not alter these conclusions. Composition and structure of the mantle, Mantle processes, Structure of the Earth, Dynamics: convection currents, and mantle plumes 1 INTRODUCTION The seismic structure of the lowermost mantle is dominated by two large low shear wave velocity provinces (LLSVPs) below Africa and the Pacific, where shear wave speed decreases by a few per cent compared to its horizontal average. LLSVPs are observed by all global seismic tomography models, with some differences in the detailed structure and amplitude of the velocity drop (for a comparison of models published between 2005 and 2010, see Ritsema et al.2011). Their presence is further supported by a variety of seismic data and methods (Lekić et al.2012), suggesting that they are not simple artefacts of tomographic models. In particular, because they are seen by normal mode data (Ishii & Tromp 1999; Trampert et al.2004; Mosca et al.2012), they are unlikely the result of wave-propagation complexities. While the presence of LLSVPs is generally accepted, their nature is still a matter of debate (for recent discussions on this problem, see Davies et al.2015; Deschamps et al.2015; Garnero et al.2016). LLSVPs were first associated with superplumes or clusters of smaller plumes, that is, buoyant structures triggered by purely thermal (hot) anomalies. Several hints, however, including the anticorrelation between shear wave and bulk-sound velocity anomalies (e.g. Ishii & Tromp 1999; Masters et al.2000; Trampert et al.2004), and the de-correlation between density and seismic velocity anomalies (Trampert et al.2004), indicated that LLSVPs cannot be explained by purely thermal anomalies, but partially originate from compositional anomalies (Trampert et al.2004; Mosca et al.2012). Alternatively, it has been proposed that the anticorrelation between shear wave and bulk-sound speeds could be mostly explained by complex wave-propagation effects (Schuberth et al.2012), or by the presence of post-perovskite (pPv; Davies et al.2012). Thus, a clear diagnosis on the nature of LLSVPs, purely thermal or thermochemical, cannot be based on current seismic velocity data only, but requires additional observational constraints independent of seismic velocities. Two promising observables are long-period electromagnetic C-responses, which are related to the distribution of electrical conductivity in the deep mantle and are different for purely thermal and thermochemical structures (Deschamps & Khan 2016), and seismic attenuation, which is mostly sensitive to temperature and can be inferred from waveform data (e.g. Konishi et al.2017). Another potential constraint is the core–mantle boundary (CMB) topography. Using Very Long Baseline Interferometry (VLBI) data, Gwinn et al. (1986) showed that the CMB dynamical flattening (strictly speaking, the dynamical flattening of the outer core), that is, the spherical harmonic degree and order l = 2 and m = 0 of CMB topography, is larger than its theoretical value for a hydrostatic Earth model by 5 per cent. If more complexities are taken into account, the excess CMB dynamical flattening reduces to 3.8 per cent (Mathews et al.2002). Seismic normal modes data further indicate that the peak-to-peak amplitude of the spherical harmonic degree l = 2 in CMB topography should not exceed 5.0 km (Koelemeijer et al.2012). Analysis of seismic phases reflecting on the CMB underside (PKKP) and upperside (PcP) suggests that 95 per cent of the CMB topography is distributed between −4.0 and 4.0 km (Garcia & Souriau 2000). Several studies, using different CMB transmitted or reflected seismic phases, have also attempted to draw global maps, but generally disagree on both the pattern and amplitude of CMB topography (Morelli & Dziewonski 1987; Doornbos & Hilton 1989; Sze & van der Hilst 2003; Tanaka 2010). For spherical harmonic degrees up to l = 4, the peak-to-peak amplitude ranges from 3.5 km (Koper et al. 2003) to 12 km (Morelli & Dziewonski 1987). Insights on CMB topography may further be obtained from geodynamic models. In a fluid animated by convection, flow generates normal stress at its interfaces, which in turn induces dynamic topography (e.g. Ricard et al.1984; Hager & Richards 1989). Plumes stretch the base and push the surface of the fluid upwards (Fig. 1a), inducing positive topography on these boundaries. In contrast, downwellings pull the surface and push the base of the system downwards, leading to negative topography (Fig. 1b). Dynamic topography is sensitive to mantle viscosity and its lateral variations, and to the density jump at the boundaries. In the Earth's mantle, the detailed CMB dynamic topography thus strongly depends on the mode of mantle dynamics, in particular on lateral variations in density due to large-scale chemical heterogeneities. The 2-D-Cartesian simulations of Hansen & Yuen (1989) indicated that piles of dense material at the bottom of the system strongly depress the CMB, inducing negative topography (Figs 1c and d). Lassak et al. (2007) pointed out that in models combining large thermal viscosity contrast and piles with small chemical density contrast, the negative CMB topography induced by slabs is much larger than that caused by the pile. Lassak et al. (2010) further showed that imposing a large thermal viscosity contrast flattens out the topography induced by thermochemical piles, because these regions are hotter, and thus less viscous, than surroundings. Using models of convection in which density anomalies are scaled from seismic tomography, Yoshida (2008) showed that if LLSVPs are purely thermal structures, they should result in positive topography with amplitude around 5 km, while, on the contrary, if they are thermochemical piles, they would induce local depressions. These models, however, are limited to moderate chemical density contrasts, around 60 kg m−3 and lower, and assume that the viscosity is independent of the type of material, dense or regular. Figure 1. View largeDownload slide Influence of mantle dynamics on core–mantle boundary (CMB) topography. Thin black arrows denote the stress applied on the CMB. (a) Thermal plumes pull the CMB upwards, inducing positive topography. (b) Slabs stretch the CMB downwards, inducing negative topography. (c,d) Thermochemical piles depress the CMB, inducing negative topography with amplitude that depends on the chemical density contrast between the material of these piles and the regular mantle. Some elements of this sketch are taken from Garnero et al. (2016). Figure 1. View largeDownload slide Influence of mantle dynamics on core–mantle boundary (CMB) topography. Thin black arrows denote the stress applied on the CMB. (a) Thermal plumes pull the CMB upwards, inducing positive topography. (b) Slabs stretch the CMB downwards, inducing negative topography. (c,d) Thermochemical piles depress the CMB, inducing negative topography with amplitude that depends on the chemical density contrast between the material of these piles and the regular mantle. Some elements of this sketch are taken from Garnero et al. (2016). In this study, we performed models of convection in 3-D-spherical geometry, including models with chemical density contrast larger than 100 kg m−3, which may be more appropriate for LLSVPs, models with an intrinsic viscosity ratio between primordial and regular materials, whose effects on CMB dynamic topography have not yet been investigated, and models accounting for the pPv phase transition. Results from these models provide a criteria to discriminate between possible models of LLSVPs, purely thermal or slightly thermochemical on one hand, or strongly thermochemical, that is, with chemical density contrast around 100 kg m−3 or more, on the other hand. 2 MANTLE CONVECTION MODELLING We performed numerical simulations of purely thermal and thermochemical convection in 3-D-spherical geometry using StagYY (Tackley 1998, 2008), a software that solves the conservation equations of energy, momentum, and mass for a compressible, anelastic fluid with infinite Prandtl number. Thermochemical models distinguish two types of material, regular (pyrolitic) mantle and primordial material, and further solve the conservation equation of composition. Primordial material accounts for chemical heterogeneities that may be present at the bottom of the mantle as a result of early differentiation. An additional source of chemical heterogeneity, which we do not take into account in our models, is the recycling of subducted oceanic crust (MORB). Recent models of convection (Li et al.2014a) indicate that subducted MORB may be partially incorporated into pre-existing reservoirs of primordial material, and partially re-entrained to the surface by plumes, in agreement with the Basal Melange (BAM) scenario (Tackley 2012). Because subducted MORB is denser than pyrolitic mantle, their effect on CMB topography is likely similar to that of primordial material, and neglecting its contribution may not alter our conclusions. For details on the numerical methods employed to solve the conservation equations, we refer to Tackley (2008). The setup used in this paper is essentially that used in Deschamps et al. (2015) except for the viscosity law, which now includes yield stress, and for the explored range of control parameters. Details of this setup can be found in previous studies (Deschamps & Tackley 2009), but for clarity we discuss them again below. Spherical geometry is modelled by a set of Yin and Yang grids (Kageyama & Sato 2004) with resolution equivalent to a spherical grid with 256 points in latitude and 512 points in longitude at each depth. The vertical resolution is fixed to 128 or 64 nodes for models with and without pPv, respectively. Grid refinement is imposed at the top and at the bottom of the shell, allowing a good resolution of reservoirs of primordial material and thermal boundary layers. The ratio between the radius of the core and the total radius is set to its Earth value, that is, f = 0.55. The bottom and surface boundaries are free slip. The system is heated both from the bottom and from within, and the heating rate of primordial material is increased by a factor 10 compared to that of regular material. The total mantle heating rate is equivalent to a surface heat flux of 10 mW m−2, which is about half the value expected for Earth's mantle (Jaupart et al.2015), and leads to a ratio of internal to basal heating close to 0.3. Compressibility generates additional sinks and sources of heat, that are controlled by the dissipation number, Di. Because it depends on thermal expansion, Di varies with depth. In all our models the surface value of the dissipation number is fixed to DiS = 1.2, and the depth variations of thermal expansion (Fig. 2) imply that its volume average is equal to 0.43. Figure 2. View largeDownload slide Thermodynamical reference state as a function of depth. Profiles shown are for temperature (Tref), density (ρref), thermal expansion (αref) and Grüneisen parameter (γref). The discontinuity in density at non-dimensional depth of 0.228 corresponds to the limit between upper and lower mantle. Figure 2. View largeDownload slide Thermodynamical reference state as a function of depth. Profiles shown are for temperature (Tref), density (ρref), thermal expansion (αref) and Grüneisen parameter (γref). The discontinuity in density at non-dimensional depth of 0.228 corresponds to the limit between upper and lower mantle. The conservation equations involve the definition of a reference state, which here consists of a set of vertical profiles for density, temperature, and thermal expansion calculated using thermodynamic relationships for the Earth's mantle (Tackley 1998). Fig. 2 plots the variations of the reference state with depth. Note that thermal expansion decreases by a factor 5 from surface to bottom, and that the Grüneisen parameter, γref, varies in such a way that its product with reference density is constant with depth. Temperature is scaled with respect to the super-adiabatic temperature difference across the system, ΔTS, and density and thermal expansion are scaled with respect to their surface values. Table 1 lists the thermodynamic properties used to build the reference state, together with other scalings and properties. Table 1. Parameters and scalings of numerical experiments. Parameter  Symbol  Value  Units  Non-dimensional  Non-dimensional parameters          Reference Rayleigh number  RaS      3.0 × 108  Surface dissipation number  DiS      1.2  Volume average dissipative number  Di      0.43  Total internal heating  HC  10  mW m−2  4.8  Compositional heating ratio  ΔHC  10                Compositional parameters (thermochemical cases only)          Buoyancy ratio  Bz      0.15–0.70  Volume fraction of dense material (%)  Xprim      2.5–7.0            Physical & thermodynamical parameters          Acceleration of gravity  g  9.81  m s−2  1.0  Mantle thickness  D  2891  km  1.0  Reference adiabat  Tas  1600  K  0.64  Super-adiabatic temperature difference  ΔTS  2500  K  1.0  Surface density  ρS  3300  kg m−3  1.0  Surface thermal expansion  αS  5.0 × 10−5  K−1  1.0  Surface thermal diffusivity  κS  6.24 × 10−7  m2 s−1  1.0  Heat capacity  CP  1200  J kg−1 K−1  1.0  Surface conductivity  kS  3.0  W m−1 K−1  1.0  Surface Grüneisen parameter  γS  1.091      Density jump at z = 660 km  Δρ660  400  kg m−3  0.1212  Clapeyron slope at z = 660 km  Γ660  −2.5  MPa K−1  −0.0668  Post-perovskite density jump  ΔρpPv  62  kg m−3  0.0188  Clapeyron slope of post-perovskite  ΓpPv  13  MPa K−1  0.3474  CMB temperature  TCMB  3750  K  1.5  Density jump at CMB  ΔρCMB  5280  kg m−3  1.6            Viscosity law          Reference viscosity  η0  1.6 × 1021  Pa s  1.0  Viscosity ratio at z = 660 km  Δη660  30      Logarithmic thermal viscosity ratio  Ea  13.816–23.026      Logarithmic vertical viscosity ratio  Va  2.303      Compositional viscosity ratio  ΔηC  1–102      Surface yield stress  σ0  290  MPa  7.5 × 106  Yield stress gradient  $${\dot{\sigma }_z}$$  0.01  Pa/Pa  0.01  Parameter  Symbol  Value  Units  Non-dimensional  Non-dimensional parameters          Reference Rayleigh number  RaS      3.0 × 108  Surface dissipation number  DiS      1.2  Volume average dissipative number  Di      0.43  Total internal heating  HC  10  mW m−2  4.8  Compositional heating ratio  ΔHC  10                Compositional parameters (thermochemical cases only)          Buoyancy ratio  Bz      0.15–0.70  Volume fraction of dense material (%)  Xprim      2.5–7.0            Physical & thermodynamical parameters          Acceleration of gravity  g  9.81  m s−2  1.0  Mantle thickness  D  2891  km  1.0  Reference adiabat  Tas  1600  K  0.64  Super-adiabatic temperature difference  ΔTS  2500  K  1.0  Surface density  ρS  3300  kg m−3  1.0  Surface thermal expansion  αS  5.0 × 10−5  K−1  1.0  Surface thermal diffusivity  κS  6.24 × 10−7  m2 s−1  1.0  Heat capacity  CP  1200  J kg−1 K−1  1.0  Surface conductivity  kS  3.0  W m−1 K−1  1.0  Surface Grüneisen parameter  γS  1.091      Density jump at z = 660 km  Δρ660  400  kg m−3  0.1212  Clapeyron slope at z = 660 km  Γ660  −2.5  MPa K−1  −0.0668  Post-perovskite density jump  ΔρpPv  62  kg m−3  0.0188  Clapeyron slope of post-perovskite  ΓpPv  13  MPa K−1  0.3474  CMB temperature  TCMB  3750  K  1.5  Density jump at CMB  ΔρCMB  5280  kg m−3  1.6            Viscosity law          Reference viscosity  η0  1.6 × 1021  Pa s  1.0  Viscosity ratio at z = 660 km  Δη660  30      Logarithmic thermal viscosity ratio  Ea  13.816–23.026      Logarithmic vertical viscosity ratio  Va  2.303      Compositional viscosity ratio  ΔηC  1–102      Surface yield stress  σ0  290  MPa  7.5 × 106  Yield stress gradient  $${\dot{\sigma }_z}$$  0.01  Pa/Pa  0.01  View Large Because the fluid properties (density, viscosity, thermal diffusivity, and thermal expansion) are allowed to vary throughout the system, the definition of the Rayleigh number is non-unique. In our simulations, we prescribed a reference Rayleigh number Ra0, defined at surface values of the thermodynamic parameters and reference viscosity η0. In all simulations shown here, Ra0 is set to 3.0 × 108, leading to an effective Rayleigh number (i.e. the Rayleigh number at the volume average viscosity) between 106 and 2.0 × 106, depending on the case. In thermochemical cases, the compositional field is modelled with a collection of 100 million particle tracers (200 million for cases with pPv), leading to an average number of tracers per cell of 15, which is enough to properly model entrainment (Tackley & King 2003). Tracers are of two types, modelling the regular mantle and primordial material, respectively, and are advected following a fourth-order Runge–Kutta method. At each time step, the compositional field is inferred from the concentration C of particles of primordial material in each cell. This concentration varies between 0 (corresponding to a cell filled with regular material only) and 1 (for a cell filled with primordial material only). The primordial material is initially distributed in a basal layer, the thickness of which is controlled by fixing the volume fraction of dense material, Xprim. In most of our calculations, we fixed Xprim to 3.5 per cent, corresponding to an upper estimate of the volume of LLSVPs (Hernlund & Houser 2008). In addition, we conducted calculations for Xprim = 2.5 per cent and Xprim = 7.0 per cent. The primordial material is assumed to be denser than the regular (pyrolitic) mantle, and the density contrast between the two materials is controlled by the buoyancy ratio. In our calculations, the buoyancy ratio is defined with respect to a reference density that increases with depth following a thermodynamical model of Earth's mantle,   $${B_z} = \frac{{\Delta {\rho _c}\left( z \right)}}{{{\alpha _S}\rho \left( z \right)\Delta {T_S}}},$$ (1)where Δρc(z) is the density contrast between dense and regular material, αS the surface thermal expansion, ρ(z) the reference density at depth z, and ΔTS the super-adiabatic temperature jump. In our models, the reference density increases with depth but Bz remains constant, implying that the chemical density contrast increases with depth proportionally to the reference density. With this definition, the buoyancy ratio required to obtain a chemical stratification of the system is smaller than when the surface density is taken as reference density. For equivalent chemical density contrast at the bottom of the system, however, results of the simulations do not differ substantially. We set Bz to 0.23 in most of our models. At the bottom of the system, and taking αS = 5.0 × 10−5 K−1, ρbot = 4950 kg m−3 and ΔTS = 2500 K, this leads to a density contrast between dense and regular material around 142 kg m−3, that is, about 2.5 per cent of PREM density (Dziewonski & Anderson 1981) at these depths, in agreement with the peak-to-peak chemical density anomalies mapped by normal modes data (Trampert et al.2004; Mosca et al.2012). Because the buoyancy ratio is certainly the most sensitive parameter controlling the dynamics of the system, we performed additional calculations with Bz between 0.15 and 0.70. This latter value allows maintaining a layer of primordial material that covers nearly all the CMB. Viscosity is allowed to vary with depth, temperature, and composition. An additional viscosity ratio Δη660 = 30 is added at the 660 km phase transition, and to avoid the formation of stagnant lid at the top of the system, we impose a yield stress. The viscosity is then fully described by   $$\eta \ = \frac{1}{{\frac{1}{{{\eta _b}}} + \frac{1}{{{\eta _Y}}}}},$$ (2)where   $${\eta _Y} = \frac{{{\sigma _0} + {{\dot{\sigma }}_z}P}}{{2\dot{e}}}\$$ (3)is the yield viscosity, and   \begin{eqnarray} {\eta _{\rm{b}}}\ \left( {z,T,C} \right) &=& {\eta _0}\ \left[ {1 + 29H\left( {z - 660} \right)} \right]\nonumber\\ &&\times \exp \left[ {{V_{\rm{a}}}\frac{z}{D} + {E_{\rm{a}}}\frac{{\Delta {T_{\rm{S}}}}}{{\left( {T + {T_{{\rm{off}}}}} \right)}} + {K_{\rm{a}}}C} \right].\end{eqnarray} (4) The yield viscosity (eq. 3) is defined from the yield stress, $${\sigma _Y} = \ {\sigma _0} + {\dot{\sigma }_z}P$$, and the second invariant of the stress tensor, $$\dot{e}.$$ The yield stress is set to σ0 at the surface and increases with pressure following a gradient $${\dot{\sigma }_z}$$. In eq. (4), η0 is a reference viscosity, H the Heaviside step function, z the depth, D the mantle thickness, ΔTS the superadiabatic temperature difference across the system and Toff the temperature offset, which is added to the temperature to reduce the viscosity jump across the top thermal boundary layer. The reference viscosity η0 is defined for the surface value of the reference adiabat (i.e. Tas = 0.64ΔTS), and at regular composition (C = 0). The viscosity variations with temperature are controlled by Ea, modelling the activation energy. To quantify the thermally induced increase of viscosity, we define a potential thermal viscosity ratio as ΔηT = exp(Ea). However, due to the adiabatic increase of temperature and to the temperature offset, which we fixed to Toff = 0.88ΔTS, the effective top-to-bottom thermal viscosity contrast is smaller than ΔηT by about two orders of magnitude. Here, we fixed Ea to 20.723 (corresponding to ΔηT = 109) in most simulations, and did additional calculations with values of Ea in the range 13.816–23.032 (106 ≤ ΔηT ≤ 1010). The viscosity variations with depth are controlled by Va, modelling the activation volume. In all simulations performed here, we fixed the value of Va to 2.303, leading to an increase of viscosity from top to bottom by a factor 300. This includes the viscosity jump at 660 km, but excludes the decrease due to adiabatic increase of temperature and the thermally induced increase in thermal boundary layers. The viscosity variations with composition are controlled by the parameter Ka, and the viscosity ratio between primordial and regular material (or chemical viscosity ratio) is given by ΔηC = exp(Ka). The viscosity of bridgmanite is larger than that of ferro-periclase by up to 3 orders of magnitude (Yamazaki & Karato 2001). If LLSVPs are enriched in bridgmanite by up to about 10 per cent, as suggested by probabilistic tomography (Trampert et al.2004; Mosca et al.2012), they should be intrinsically more viscous than surrounding mantle. Recent numerical experiments (Li et al.2014b; Deschamps et al.2015) further indicate that intrinsically more viscous primordial material increases the topography of reservoirs of dense material and the sharpness of their edges. The value of the chemical viscosity ratio depends on poorly constrained parameters, including the rock deformation model and the distribution of phases (bridgmanite and ferro-periclase) within the aggregate, and is thus difficult to estimate. In most of the thermochemical models presented in this study, we impose primordial material to be more viscous than regular material with ΔηC = 30. Finally, it should be pointed out that our viscosity law does not account for viscosity variations with grain size, which may influence mantle dynamics and evolution. The detailed effects of grain size remain a matter of debate. For diffusion creep, viscosity is usually assumed to increase with grain size following a power law with exponent in the range 2–3 (e.g. Korenaga & Karato 2008). Large grain size further favours the transition from diffusion to dislocation creep regimes, and it has been suggested that this transition may occur for grain size as low as 0.1–1 mm (Li et al.1996). Due to various physical and chemical processes, grain size may strongly vary in space and time. For instance, because in early Earth grain growth occurred in a partially molten environment, initial grain size may have been rather large, while recrystallization, in particular during the ringwoodite to bridgmanite phase transition, leads to smaller grain size. Regions with larger grain size in primordial mantle had larger viscosity, and depending on the vigour of convection, they may have mixed efficiently or, on the contrary, poorly with regions of lower viscosity (Solomatov & Reese 2008). If thermochemical reservoirs in our models resulted from early mantle differentiation, their intrinsic viscosity may differ from that of the regular mantle by up to a few orders of magnitude. A simple way to account for this effect would then be to systematically vary the chemical density contrast, ΔηC, up to very large values, for example, 103 and more. This, however, neglects recrystallization effects that may affect the primordial material during its circulation in the mantle. Fig. 3, plotting viscosity profiles along the initial temperature profile, illustrates the radial viscosity variations in our models when thermal and chemical effects are taken into account. The effect of intrinsically more viscous (ΔηC = 30) dense material at the bottom of the mantle is denoted by the dashed curve. Note that for Ea = 20.723, which is the case of most of our models, the regular increase with depth (Va = 2.303) in the lower mantle is nearly compensated by the decrease due to adiabatic increase of temperature. Figure 3. View largeDownload slide Radial viscosity profiles along initial temperature profile for purely thermal models or thermochemical models with chemical viscosity contrast ΔηC = 1 (plain curves), and thermochemical models with ΔηC = 30 (dotted curve). Two values of the thermal viscosity contrast (controlled by the parameter Ea) are shown. Figure 3. View largeDownload slide Radial viscosity profiles along initial temperature profile for purely thermal models or thermochemical models with chemical viscosity contrast ΔηC = 1 (plain curves), and thermochemical models with ΔηC = 30 (dotted curve). Two values of the thermal viscosity contrast (controlled by the parameter Ea) are shown. The phase change at 660 km is modelled with a discontinuous phase transition controlled by defining a point on the phase boundary and a Clapeyron slope, Γ660. Here, we imposed z = 660 km and T = 1900 K as anchor point, and Γ660 = −2.5 MPa K−1. The density contrast at the phase transition is Δρ660 = 400 kg m−3 and is scaled with the surface density. In cases that include pPv phase, the pPv phase transition is modelled with a reference temperature of 2700 K at 2700 km depth. Lateral deviations in the transition depth are determined using the phase function approach of Christensen & Yuen (1985). The Clapeyron slope and the density contrast are set to ΓpPv = 13 MPa K−1 (Tateno et al.2009) and ΔρpPv = 62 kg m−3 (corresponding to a relative contrast of ∼1.0 per cent, Murakami et al.2004), respectively. Determining the stability field of pPv further requires fixing the temperature at the CMB, TCMB. Based on the numerical simulations of Li et al. (2015), which pointed out that values of TCMB around 3750 K describe best lower mantle structure and dynamics, we assumed TCMB = 3750 K. The dynamic topography at the CMB, hCMB, is calculated from the normal stress, σzz, induced by the flow on this boundary, and takes into account self-gravitational effects following Zhang & Christensen (1993). At each point of the CMB, dynamic topography is given by   $${h_{{\rm{CMB}}}} = \frac{{{\sigma _{zz}} + {{\rm{\Phi }}_{{\rm{CMB}}}}\Delta {\rho _{{\rm{CMB}}}}}}{{\Delta {\rho _{{\rm{CMB}}}}g}},$$ (5)where ΔρCMB is the density difference between the mantle and the outer core, ΦCMB the perturbation of gravitational potential at CMB, and g the acceleration of gravity. The normal stress is given by   $$\ {\sigma _{zz}} = \ 2\eta \left( {\frac{{\partial {v_r}}}{{\partial r}} - \frac{1}{3}\nabla \cdot {\boldsymbol{v}}} \right),$$ (6) where v is velocity and vr its vertical component. Note that the divergence of velocity in the right hand side of eq. (6) is due to compressibility. Gravitational potential perturbations are obtained by solving Poisson's equation for a distribution of density anomalies (with respect to a reference state) that includes thermal and chemical effects, and variations in CMB topography. For comparison, we also calculated dynamic topography neglecting self-gravitation, but only found small differences, mostly in amplitude, with the topography accounting for self-gravity. For a given set of parameters, the rms amplitude of dynamic topography is, on average, smaller by 5 per cent if self-gravitational effects are included. It is usually assumed that the mechanical boundary condition, free or free-slip, does not affect the boundary topography, because this topography is very small compared to the lateral lengthscale, and therefore slopes are also very small. This has been well tested for the surface boundary (Crameri et al.2012). For large slopes, such as at subduction zones trenches, differences may however be important, and a free surface may help single-sided subduction (e.g. Crameri & Tackley 2014). The initial temperature distribution is built by adding small random perturbations to an adiabatic radial profile with thin superadiabatic boundary layers at the top and bottom of the shell. Purely thermal experiments are run up to a dimensional time of 4.5 Gyr. Thermochemical experiments start with a long (up to 1.5–2.0 Gyr) transient phase during which the bottom layer of dense material heats up, and are run up to a dimensional time of 15 Gyr. Models do not reach stationary state, that is, the location of plumes and downwellings vary with time. However, after a period of time that depends on each case, average velocity and surface and bottom heat fluxes stabilize and oscillate around constant values, thus defining a dynamic quasi-equilibrium. All the snapshots presented and discussed in Section 3 are taken in the quasi-equilibrium stage. In thermochemical models, the detailed shape of thermochemical piles varies with time, but the selected snapshots are generally representative of the thermochemical structure obtained for each case. Snapshots are taken at 2.7 Gyr for purely thermal cases, and 11.2 or 15.0 Gyr for thermochemical cases. Note that the experiments discussed here are not designed to model the detailed evolution of the Earth's mantle, for which accurate initial conditions are not yet known. This would further require adding time-decreasing radioactive heating and CMB temperature in our models. Instead, the aim is to obtain various characteristic thermochemical structures, for example, plume clusters or thermochemical piles, to estimate the influence of these structures on CMB dynamic topography. Thus, time indication in our experiments should not be used to discriminate between possible models of mantle dynamics, or to interpret early mantle evolution. 3 RESULTS Using the setup described in Section 2, we performed a series of 21 numerical experiments with various values of the buoyancy ratio, volume fraction of dense material, and thermal and chemical viscosity ratios (Table 2). All models are in the mobile-lid regime, with horizontal surface velocities of around 1 cm yr−1. The surface heat flux is in the range 35–40 mW m−2, depending on the model. These are rather low values compared to heat flux measured at the surface of the Earth, the average of which is around 86 mW m−2 (Davies 2013). It should however be kept in mind that the observed surface heat flux is likely overestimating the heat flux at the top of the mantle because it includes a contribution related to crustal heat production, which is difficult to estimate, but likely around 8 TW (Jaupart et al.2015). In cratonic regions, for instance, this contribution may be large and the mantle heat flux could be as low as 15–20 mW m−2 (Jaupart & Mareschal 2007). In addition, because our models do not include plate tectonics, they do not distinguish between continental areas and oceanic regions, where observed heat flux is usually highest. Finally, our models do not include secular cooling of the mantle, which contributes up to half of Earth's surface heat loss (Jaupart et al.2015). The bottom heat flux is around 80–100 mW m−2 depending on the model, equivalent to a total power around 12–15 TW, that is, in good agreement with the expected core heat flow, in the range 5–17 TW (Jaupart et al.2015). Table 2. Run properties and corresponding output core–mantle boundary (CMB) topography. The prefix T and TC denote purely thermal and thermochemical cases, respectively. Bz is the buoyancy ratio, Xprim the volume fraction of dense material, ΔηC the intrinsic (chemical) viscosity ratio, and Ea the logarithmic thermal viscosity contrast (eq. 3). Other properties are listed in Table 1 and are the same for all runs. CMB topography statistics includes rms of global topography, and average topography 〈h〉 in primordial reservoirs or plume clusters. For primordial reservoirs, we distinguish the average depth of their interiors 〈h−〉, and the average elevation of their surrounding ridges, 〈h+〉. All topographies are given in km. Run  Bz  Xprim  ΔηC  Ea  pPv  Resolution  hCMB (km)                rms  〈h−〉  〈h+〉  T1  –  –  –  20.723  no  512 × 256 × 64  2.7  –  1.9  T2  –  –  –  13.816  no  512 × 256 × 64  2.7  –  2.9  T3  –  –  –  23.026  no  512 × 256 × 64  2.4  –  1.6  T1-pPv  –  –  –  20.723  yes  512 × 256 × 128  2.4  –  2.5  T2-pPv  –  –  –  13.816  yes  512 × 256 × 128  3.1  –  3.9  TC1  0.23  3.5  30  20.723  no  512 × 256 × 64  2.8  −1.5  1.5  TC2  0.23  3.5  30  13.816  no  512 × 256 × 64  2.6  −1.7  1.8  TC3  0.23  3.5  30  23.026  no  512 × 256 × 64  2.7  −1.3  1.4  TC4  0.15  3.5  30  20.723  no  512 × 256 × 64  3.1  −0.6  1.3  TC5  0.50  3.5  30  20.723  no  512 × 256 × 64  3.0  −2.4  2.0  TC6  0.70  3.5  30  20.723  no  512 × 256 × 64  3.4  −2.5  2.8  TC7  0.23  3.5  1  20.723  no  512 × 256 × 64  2.0  −0.7  1.0  TC8  0.23  3.5  10  20.723  no  512 × 256 × 64  2.4  −1.2  1.4  TC9  0.23  3.5  100  20.723  no  512 × 256 × 64  2.5  −1.7  1.5  TC10  0.23  2.5  30  20.723  no  512 × 256 × 64  2.4  −1.7  1.3  TC11  0.23  7.0  30  20.723  no  512 × 256 × 64  2.6  −1.3  1.4  TC1-pPv  0.23  3.5  30  20.723  yes  512 × 256 × 128  2.8  −1.9  1.5  TC2-pPv  0.23  3.5  30  13.816  yes  512 × 256 × 128  3.8  −2.2  2.4  TC4-pPv  0.15  3.5  30  20.723  yes  512 × 256 × 128  3.0  −0.9  1.3  TC6-pPv  0.70  3.5  30  20.723  yes  512 × 256 × 128  3.9  −2.5  2.7  TC7-pPv  0.23  3.5  1  20.723  yes  512 × 256 × 128  2.7  −0.9  1.2  Run  Bz  Xprim  ΔηC  Ea  pPv  Resolution  hCMB (km)                rms  〈h−〉  〈h+〉  T1  –  –  –  20.723  no  512 × 256 × 64  2.7  –  1.9  T2  –  –  –  13.816  no  512 × 256 × 64  2.7  –  2.9  T3  –  –  –  23.026  no  512 × 256 × 64  2.4  –  1.6  T1-pPv  –  –  –  20.723  yes  512 × 256 × 128  2.4  –  2.5  T2-pPv  –  –  –  13.816  yes  512 × 256 × 128  3.1  –  3.9  TC1  0.23  3.5  30  20.723  no  512 × 256 × 64  2.8  −1.5  1.5  TC2  0.23  3.5  30  13.816  no  512 × 256 × 64  2.6  −1.7  1.8  TC3  0.23  3.5  30  23.026  no  512 × 256 × 64  2.7  −1.3  1.4  TC4  0.15  3.5  30  20.723  no  512 × 256 × 64  3.1  −0.6  1.3  TC5  0.50  3.5  30  20.723  no  512 × 256 × 64  3.0  −2.4  2.0  TC6  0.70  3.5  30  20.723  no  512 × 256 × 64  3.4  −2.5  2.8  TC7  0.23  3.5  1  20.723  no  512 × 256 × 64  2.0  −0.7  1.0  TC8  0.23  3.5  10  20.723  no  512 × 256 × 64  2.4  −1.2  1.4  TC9  0.23  3.5  100  20.723  no  512 × 256 × 64  2.5  −1.7  1.5  TC10  0.23  2.5  30  20.723  no  512 × 256 × 64  2.4  −1.7  1.3  TC11  0.23  7.0  30  20.723  no  512 × 256 × 64  2.6  −1.3  1.4  TC1-pPv  0.23  3.5  30  20.723  yes  512 × 256 × 128  2.8  −1.9  1.5  TC2-pPv  0.23  3.5  30  13.816  yes  512 × 256 × 128  3.8  −2.2  2.4  TC4-pPv  0.15  3.5  30  20.723  yes  512 × 256 × 128  3.0  −0.9  1.3  TC6-pPv  0.70  3.5  30  20.723  yes  512 × 256 × 128  3.9  −2.5  2.7  TC7-pPv  0.23  3.5  1  20.723  yes  512 × 256 × 128  2.7  −0.9  1.2  View Large Table 2 further lists the rms amplitude of dynamic topography for each model. To better quantify the dynamic topography caused by plumes, downwellings, and reservoirs of dense material, we define several regions above the CMB. The plume, or plume cluster, region is defined as the area where the temperature is larger than Tplume = Tm + cplume(Tmax − Tm), where Tmax and Tm are the maximum and average temperatures above the CMB, respectively, and cplume a constant that we fixed to 0.5. The region where downwellings accumulate is defined as the region where the temperature is lower than Tslab = Tm − cslab(Tm − Tmin), where Tmin is the minimum temperature above the CMB, and cslab a constant that we fixed, again, to 0.5. In thermochemical cases, reservoirs of dense material are defined by the region where the concentration in dense particles, C, is larger than 0.5. Note that because the transition between the dense reservoir and its surroundings is very sharp, the size and extension of the reservoir do not vary with threshold values of C larger than 0.3. By contrast, fixing cplume = cslab = 0.5 is rather subjective, and changing these values leads to variations in the extents of plume and downwelling regions, inducing in turn variations in the estimated average topography within these regions. For instance, reducing cslab increases the area of downwelling regions by including regions with smoother topography, thus decreasing the average topography. Importantly, while the extent of plume and slab areas may influence the estimated average CMB topography in these regions, it does not affect the trend, that is, the fact that plumes induce elevations in the CMB, and slabs depressions. Figs 4 and 5 show snapshots of the temperature, composition and CMB dynamic topography for cases T1 and TC1, respectively, which we take as reference cases for purely thermal and thermochemical models. In case T1, plumes are generated at the CMB and are grouping in two main clusters (Figs 4a and b). Plumes are interconnected through hot ridges, and the typical angular distance between two plumes is 20°. Cold downwellings reach the CMB in the wide regions located between plume clusters. At their base, plumes generate upwards normal stress, and therefore positive dynamic topography (Figs 4c and d). In the plume clusters region, the average amplitude of dynamic topography is 1.9 km. By contrast, as they reach the CMB, downwellings induce negative dynamic topography, that is, depressions of the CMB. These depressions are strongly localized, but have large amplitude. In case T1, the slab stacking region covers less than 0.5 per cent of the total CMB area, but within this region the rms amplitude of dynamic topography is 17 km. A distribution histogram of the CMB topography (Fig. 6a) indicates that throughout most of the CMB (∼90 per cent of its total surface) dynamic topography ranges between −2.0 and 3.0 km. Downwellings slightly skew the distribution towards topography deeper than −5.0 km. Note that for convenience, the histogram is truncated to the interval [−10; 10] km, implying that topography caused by downwellings is not entirely represented. However, for topography deeper than −10 km the frequency is very small, less than 10−3 for a 1 km-stacking bin. Figure 4. View largeDownload slide Snapshot of a purely thermal case without post-perovskite (T1). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.72. (c,d) Core–mantle boundary dynamic topography (colour scale). The orange and blue dotted lines in plots (c) and (d) indicate the limits of the plumes and downwelling regions, respectively. Black dotted lines in plots (b) and (c) denote the path along which topographic profiles in Figs 15(a) and (b) are drawn. Figure 4. View largeDownload slide Snapshot of a purely thermal case without post-perovskite (T1). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.72. (c,d) Core–mantle boundary dynamic topography (colour scale). The orange and blue dotted lines in plots (c) and (d) indicate the limits of the plumes and downwelling regions, respectively. Black dotted lines in plots (b) and (c) denote the path along which topographic profiles in Figs 15(a) and (b) are drawn. Figure 5. View largeDownload slide Snapshot of a thermochemical case without post-perovskite (TC1). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.67. (c,d) Isosurface of the composition for C = 0.5. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in plots (e) and (f) indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. Black dotted lines in plots (e) and (f) denote the path along which the topographic profile in Fig. 15(c) is drawn. Figure 5. View largeDownload slide Snapshot of a thermochemical case without post-perovskite (TC1). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.67. (c,d) Isosurface of the composition for C = 0.5. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in plots (e) and (f) indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. Black dotted lines in plots (e) and (f) denote the path along which the topographic profile in Fig. 15(c) is drawn. Figure 6. View largeDownload slide Distribution histograms of dynamic topography at the CMB for cases T1 (a), TC1 (b), T1-pPv (c) and TC1-pPv (d) (see Table 2 for the properties of these cases). Frequency is normalized to the total area of the CMB. For convenience, histograms are truncated to the interval [−10; 10] km, implying that topography caused by downwellings is not entirely represented. Figure 6. View largeDownload slide Distribution histograms of dynamic topography at the CMB for cases T1 (a), TC1 (b), T1-pPv (c) and TC1-pPv (d) (see Table 2 for the properties of these cases). Frequency is normalized to the total area of the CMB. For convenience, histograms are truncated to the interval [−10; 10] km, implying that topography caused by downwellings is not entirely represented. In case TC1, the thermochemical structure at the bottom of the system is dominated by two large antipodal reservoirs of primordial material culminating about 500 km above the CMB (Figs 5c and d). Reservoirs are hotter than average, and plumes are generated at their edges and in their central areas (Figs 5a and b). Large downwellings reach the CMB in regions located in between the two thermochemical reservoirs. Each reservoir of dense material induces a large depression throughout most of its base, surrounded by a thin ridge of positive topography on its edges (Figs 5e and f). On average, the depth of the depression and the height of the ridge are equal to −1.5 km and 1.5 km, respectively (Table 2). As in case T1, downwellings induce deep local depressions with rms amplitude of 6.2 km. Interestingly, the CMB topography distribution is different from that obtained for purely thermal model T1 (Fig. 6b), and now follows a bimodal distribution, with a peak between −2.0 and −1.0 km, corresponding to the interior of the reservoir of primordial material, and a second, less pronounced, peak between 1.0 and 2.0 km, corresponding to its edges. Overall, most of the CMB area (>90 per cent) has a topography between −3.0 and 4.0 km. As in case T1, the distribution is slightly skewed by downwellings towards topography deeper than −5.0 km. Figs 7 and 8 show snapshots for cases T1-pPv and TC1-pPv, which are similar to T1 and TC1, respectively, except that they include the pPv phase at the bottom of the system. Adding the pPv does not substantially modify the convection patterns obtained for cases T1 and TC1, except that downwellings are globally less well developed. For the CMB temperature we imposed, pPv is stable outside the plume clusters regions for the purely thermal case (T1-pPv; Figs 7c and d) and outside the reservoir of primordial material for the thermochemical case (TC1-pPv; Figs 8e and f). Note that in these models, the viscosities of pPv and bridgmanite are similar. Weak pPv, with a viscosity smaller than that of bridgmanite by three orders of magnitude, may slightly modify the convection pattern, in particular for the thermochemical cases (e.g. Li et al.2015). The effects on dynamic topography are similar to those obtained without the pPv phase. In purely thermal models, plume clusters induce positive topography with an average height of ∼2.4 km (Figs 7e and f), while in thermochemical models, reservoirs of primordial material cause a depression of the CMB in their interiors and a ridge on their edges, with average depth and height equal to −1.9 km and 1.5 km, respectively (Figs 8g and h). Again, downwellings cause deep local troughs. Distributions of topography are slightly more dispersed, with most of the topography ranging between −4.0 km and 4.0 km for model T1-pPv, and between −4.0 and 5.0 km for model TC1-pPv, but are globally unchanged (Figs 6c and d). Figure 7. View largeDownload slide Snapshot of a purely thermal case with post-perovskite (T1-pPv). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.80. (c,d) Stability field of the post-perovskite phase. (e,f) Core–mantle boundary dynamic topography (colour scale). The orange and blue dotted lines in plots (e) and (f) indicate the limits of the plumes and downwelling regions, respectively. Figure 7. View largeDownload slide Snapshot of a purely thermal case with post-perovskite (T1-pPv). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.80. (c,d) Stability field of the post-perovskite phase. (e,f) Core–mantle boundary dynamic topography (colour scale). The orange and blue dotted lines in plots (e) and (f) indicate the limits of the plumes and downwelling regions, respectively. Figure 8. View largeDownload slide Snapshot of a thermochemical case with post-perovskite (TC1-pPv). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.67. (c,d) Isosurface of the composition for C = 0.5. (e,f) Stability field of the post-perovskite phase. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in plots (g) and (h) indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. Figure 8. View largeDownload slide Snapshot of a thermochemical case with post-perovskite (TC1-pPv). Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the non-dimensional temperature for T = 0.67. (c,d) Isosurface of the composition for C = 0.5. (e,f) Stability field of the post-perovskite phase. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in plots (g) and (h) indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. Increasing the volume fraction of dense material, Xprim, obviously increases the size of reservoirs of primordial material, in particular the fraction of CMB area they cover, but does not substantially modify their stability (Li et al.2014b). All other parameters being equal to those of case TC1, our calculations for Xprim = 2.5 per cent (case TC10) and Xprim = 7.0 per cent (case TC11) further indicate that dynamic topography is only slightly affected by a change in Xprim (Table 2). Other parameters of the model, in particular the buoyancy ratio and the thermal and compositional viscosity contrasts, induce more substantial changes on the flow pattern. We now further explore the consequences of these changes on the CMB dynamic topography. 3.1 Buoyancy ratio Fig. 9 shows snapshots of models with buoyancy ratio Bz = 0.15 (case TC4) and Bz = 0.7 (case TC6), corresponding to chemical density contrasts of 93 and 433 kg m−3, respectively. In case TC4, the chemical density contrast is too small to prevent strong entrainment of primordial material, leading to the formation of narrow, elevated thermochemical piles with sharp edges (Figs 9a–d). The dynamic topography induced by these piles is strongly attenuated compared to that of case TC1. The average depth of central depression reaches only 0.6 km, and the average height of surrounding ridge 1.3 km (Table 2). As observed by Lassak et al. (2007), the depressions induced by slabs are much larger, with, in case TC4, an average depth of 11.6 km. The distribution of topography is no longer bimodal but peaks between 1.0 and 2.0 km (Fig. 10a), similar to the distributions obtained for purely thermal cases. In case TC6, on the contrary, the large chemical density contrast strongly inhibits entrainment of primordial material, leading to a nearly stratified system in which the core is almost fully covered by a layer of primordial material (Figs 9e–h). Compared to case TC1, dynamic topography is more pronounced, with average depth of depression and height of ridge around −2.5 and 2.8 km, respectively (Table 2). The negative CMB topography caused by slabs, on the contrary, is attenuated, with an average depth of 4.2 km. In addition, topographic highs around the reservoir reach large values, 10 km and more, strongly skewing the distribution of topography towards high values (Fig. 10b). Figure 9. View largeDownload slide Snapshots of thermochemical cases TC4 (a–d) with buoyancy ratio B = 0.15, and TC6 (e–h), with buoyancy ratio B = 0.70. Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the composition for C = 0.65. (c,d) Core–mantle boundary dynamic topography (colour scale). (e,f) Isosurface of the composition for C = 0.5. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in topography plots indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. Black dotted lines in plots (c) and (d) denote the path along which the topographic profile in Fig. 15(d) is drawn. Figure 9. View largeDownload slide Snapshots of thermochemical cases TC4 (a–d) with buoyancy ratio B = 0.15, and TC6 (e–h), with buoyancy ratio B = 0.70. Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the composition for C = 0.65. (c,d) Core–mantle boundary dynamic topography (colour scale). (e,f) Isosurface of the composition for C = 0.5. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in topography plots indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. Black dotted lines in plots (c) and (d) denote the path along which the topographic profile in Fig. 15(d) is drawn. Figure 10. View largeDownload slide Distribution histograms of dynamic topography at the CMB for six thermochemical cases, (a) TC4, (b) TC6, (c) TC2, (d) TC3, (e) TC7 and (f) TC9 (see Table 2 for the properties of each case). Frequency is normalized to the total area of the CMB. For convenience, histograms are truncated to the interval [−10; 10] km, implying that topography caused by downwellings is not entirely represented. Figure 10. View largeDownload slide Distribution histograms of dynamic topography at the CMB for six thermochemical cases, (a) TC4, (b) TC6, (c) TC2, (d) TC3, (e) TC7 and (f) TC9 (see Table 2 for the properties of each case). Frequency is normalized to the total area of the CMB. For convenience, histograms are truncated to the interval [−10; 10] km, implying that topography caused by downwellings is not entirely represented. Overall, increasing the buoyancy ratio of primordial material increases the amplitude of dynamic topography generated by thermochemical piles and reservoirs. Additional cases including the pPv phase transition (cases TC4-pPv and TC6-pPv in Table 2) indicate that this trend is not altered when pPv is present. 3.2 Thermal viscosity contrast Thermal viscosity contrast has a strong influence on the shape and stability of reservoirs of primordial material (Deschamps & Tackley 2009; Li et al.2014b). These reservoirs are hotter, and therefore less viscous than their surroundings. This viscosity contrast increases with Ea, which in turn reduces the mechanical coupling of dense material and their entrainment by plumes. Figs 11(a)–(d) show a snapshot of thermochemical case TC2 with thermal viscosity contrast Ea = 13.816, corresponding to a maximum (or potential) bottom-to-top viscosity contrast, ΔηT, of 106. Compared to case TC1, for which Ea = 20.723 and ΔηT = 109, the shape and stability of the reservoirs of dense material are substantially altered. Reservoirs have smaller length-scale, larger topography, and sharper edges, as already observed in the simulations of Li et al. (2014b). Changes in thermal viscosity contrast strongly affect the amplitude of CMB dynamic topography beneath reservoirs. Viscosity in reservoirs decreases with increasing thermal viscosity contrast, leading to flatter topography (Table 2). The average depth of depression and height of surrounding ridge are around −1.7 and 1.8 km for case TC2 (Ea = 13.816), and −1.3 and 1.4 km for case TC3 (Ea = 23.026). We observe a similar effect for purely thermal models, that is, the topography induced by plumes decreases in amplitude with increasing thermal viscosity contrast (Table 2). In our models, the average elevation drops from 2.9 km for model T2 (Ea = 13.816) to 1.6 km for model T3 (Ea = 23.026). Interestingly, in thermochemical models, reducing the thermal viscosity contrast changes the distribution of topography. Contrary to cases TC1 (Fig. 6b) and TC3 (Fig. 10d), the distribution of topography for case TC2 (Fig. 10c) is not bimodal. In addition, it is slightly more dispersed, with most of the topography being in the range −5.0 and 5.0 km. Again, the addition of pPv does not modify this trend, as indicated by comparison between cases TC1-pPv and TC2-pPv, and T1-pPv and T2-pPv (Table 2). Figure 11. View largeDownload slide Snapshots of thermochemical cases TC2 (a–d), with Ea = 13.816, and TC7 (e–h), with ΔηC = 1. Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the composition for C = 0.65. (c,d) Core–mantle boundary dynamic topography (colour scale). (e,f) Isosurface of the composition for C = 0.5. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in topography plots indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. Figure 11. View largeDownload slide Snapshots of thermochemical cases TC2 (a–d), with Ea = 13.816, and TC7 (e–h), with ΔηC = 1. Left and right columns show two opposite views of the same snapshot. (a,b) Isosurface of the composition for C = 0.65. (c,d) Core–mantle boundary dynamic topography (colour scale). (e,f) Isosurface of the composition for C = 0.5. (g,h) Core–mantle boundary dynamic topography (colour scale). The white and blue dotted lines in topography plots indicate the limits of the reservoirs of dense material and of downwelling regions, respectively. 3.3 Chemical viscosity contrast Increasing the chemical viscosity ratio between primordial and regular materials, ΔηC, has the opposite effect to that of increasing the thermal viscosity contrast. The shape and stability of reservoirs are only slightly modified. Reservoirs are flatter if primordial and regular material have similar intrinsic viscosity (ΔηC = 1; case TC7, Figs 11e and f) than if primordial material is more viscous (e.g. case TC1, Figs 5c and d). The amplitude of CMB dynamic topography beneath the reservoirs increases with increasing intrinsic viscosity. Peak-to-peak amplitude rises from 1.8 km for case TC7 (ΔηC = 1) to 3.2 km for case TC9 (ΔηC = 102) (Table 2). For ΔηC = 1, most of CMB topography ranges between −2.0 and 3.0 km with a peak between 0 and −1.0 km (Fig. 10e), while for ΔηC = 10 or more, the CMB topography is more dispersed and follows a bimodal distribution (Figs 6b and 10f). Finally, comparison between cases TC7 and TC7-pPv shows that the addition of pPv does not alter the influence of ΔηC. 4 COMPARISON WITH SEISMOLOGICAL CONSTRAINTS A key question is whether CMB dynamic topography, provided it can be mapped from seismic data, can discriminate between purely thermal and thermochemical models of Earth's mantle. CMB topography predicted by our models strongly differ both in amplitude and pattern depending on the properties of each model, suggesting that CMB topography may bring important insights on lowermost mantle structure. Available seismic models on CMB topography, however, are limited to long wavelengths, up to spherical harmonic degree 4, for which purely thermal and thermochemical models may have similar signatures (Lassak et al.2010). In this section, we examine more in details the long-wavelength structure of our CMB topography models. 4.1 Amplitude: comparison with observed CMB topography Table 3 and Fig. 12(a) quantitatively compare CMB topography predicted by our models with available seismological constraints for different spherical harmonic filters. Normal modes data further suggest that the amplitude of spherical harmonic degree l = 2 of topography should be smaller than 5 km (Koelemeijer et al.2012). Our models fit this observation well (Table 3), with the notable exception of thermochemical models with very strong density contrast (TC5, TC6 and TC6-pPv), around 300 kg m−3 and more. These later models induce a strong chemical layering (Figs 9e and f) that is not seen by seismic velocity data, and are therefore very unlikely. Note, also, that the peak-to-peak amplitude of CMB topography predicted by purely thermal model T1 is slightly larger than the 5 km limit suggested by seismic normal modes, and that models with small density contrast (TC4 and TC4-pPv), around 90 kg m−3 and less, are close to this limit. Several studies have attempted mapping long-wavelength (up to spherical harmonic degree l = 4) CMB topography (Morelli & Dziewonski 1987; Doornbos & Hilton 1989; Koper et al.2003; Sze & van der Hilst 2003; Tanaka 2010). The resulting maps differ strongly both in pattern and amplitude, a major difficulty being that CMB topography trades off with lowermost mantle seismic velocity anomalies, and upper outer core structure. Studies including all spherical harmonic degrees up to l = 4 have peak-to-peak amplitude between 3.5 km (Koper et al.2003) and 12 km (Morelli & Dziewonski 1987). In general, for this spherical harmonic filter, the dynamic topography predicted by purely thermal models or thermochemical models with small density contrast have peak-to-peak amplitude ranging from 6 km to 11 km, somewhat larger than that of thermochemical models with moderate density contrast (from ∼100 to 200 kg m−3), which ranges between 3.5 and 7.5 km, depending on input parameters (Table 3 and Fig. 12a). Using degrees l = 2 and l = 4 only, Tanaka (2010) found a peak-to-peak amplitude of 5 km, in fair agreement with the peak-to-peak topography predicted by thermochemical models with moderate density contrast (between 2.5 and 5.5 km), but slightly lower than that of purely thermal models, which are around 6 km and more, except for model T3 (Table 3 and Fig. 12a). Figure 12. View largeDownload slide (a) Peak-to-peak amplitude of dynamic topography predicted by our models of convection (colour code) for three spherical harmonic filters (l = 2, l = 2 and 4 and l = 0–4). For comparison, amplitudes observed by Morelli & Dziewonski (1987) (MD87), Doornbos & Hilton (1989) (DH89), Sze & van der Hilst (2003) (SV03), Koper et al. (2003) (KPF03), Tanaka (2010) (T10), and Koelemeijer et al. (2012) (KDT12) are also indicated. (b) Correlation coefficient between long-wavelength (l = 0–4) shear-velocity anomalies and CMB dynamic topography predicted by each model of convection. See also Table 3 for numerical values. Figure 12. View largeDownload slide (a) Peak-to-peak amplitude of dynamic topography predicted by our models of convection (colour code) for three spherical harmonic filters (l = 2, l = 2 and 4 and l = 0–4). For comparison, amplitudes observed by Morelli & Dziewonski (1987) (MD87), Doornbos & Hilton (1989) (DH89), Sze & van der Hilst (2003) (SV03), Koper et al. (2003) (KPF03), Tanaka (2010) (T10), and Koelemeijer et al. (2012) (KDT12) are also indicated. (b) Correlation coefficient between long-wavelength (l = 0–4) shear-velocity anomalies and CMB dynamic topography predicted by each model of convection. See also Table 3 for numerical values. Table 3. Columns 2–7 list the root mean square (rms) and peak-to-peak amplitude of long-wavelength dynamic topography for four different spherical harmonic filters. For comparison, the peak-to-peak amplitudes of CMB topography observed by Morelli & Dziewonski (1987) (MD87), Doornbos & Hilton (1989) (DH89), Sze & van der Hilst (2003) (SV03), Koper et al. (2003) (KPF03), Tanaka (2010) (T10) and Koelemeijer et al. (2012) (KDT12) are also indicated. All topographies are given in kilometres. Column 8 lists the fraction of CMB area with topography between −4.0 and 4.0 km, x4km. Observed value is from Garcia & Souriau (2000) (GS00). The last column lists the correlation between long-wavelength (l = 0–4) dynamic topography and shear velocity anomalies, χh-dlnVs, predicted by each model of convection.     Rms    Peak  to peak  amplitude  x4km  χh-dlnVs    l = 2  l = 2, 4  l = 0–4  l = 2  l = 2, 4  l = 0–4    l = 0–4  Run  T1  1.3  1.3  1.6  5.1  6.6  8.5  0.94  −0.86  T2  0.7  1.9  2.1  2.0  7.1  8.7  0.89  −0.94  T3  0.8  0.9  1.2  2.7  4.1  5.8  0.95  −0.54  T1-pPv  1.3  1.3  1.5  4.5  5.9  7.7  0.91  −0.83  T2-pPv  0.8  2.1  2.3  3.4  9.8  11.5  0.81  −0.98  TC1  0.7  0.8  1.0  2.3  3.5  4.4  0.88  0.71  TC2  0.6  0.8  1.0  2.3  3.1  4.5  0.89  −0.13  TC3  0.7  1.2  1.3  2.9  5.7  6.7  0.92  0.33  TC4  1.1  1.2  1.6  4.4  5.4  7.6  0.89  −0.68  TC5  1.3  1.6  2.5  5.2  7.5  9.8  0.82  0.93  TC6  2.5  2.3  2.2  7.3  9.1  8.2  0.75  0.88  TC7  0.2  0.6  0.7  0.9  2.7  3.5  0.95  0.27  TC8  0.4  0.5  0.9  1.2  2.9  4.0  0.93  0.58  TC9  0.5  0.7  1.0  1.7  3.0  4.0  0.90  0.82  TC10  0.4  0.7  0.9  1.8  3.4  5.1  0.93  0.58  TC11  1.0  1.2  1.5  3.4  5.8  6.8  0.92  0.65  TC1-pPv  1.1  1.3  1.6  4.6  5.2  7.1  0.87  0.54  TC2-pPv  0.4  0.8  1.5  1.5  3.9  6.8  0.72  −0.29  TC4-pPv  1.1  1.2  1.5  4.8  5.6  8.6  0.89  −0.50  TC6-pPv  2.3  2.4  2.9  9.1  11.5  14.9  0.70  0.89  TC7-pPv  0.6  0.9  1.2  2.2  4.3  6.0  0.91  0.05  Observed  KDT12        ≤5.0          T10          5.0        MD87            12.0      DH89            8.0      SV03            4.0      KPF03            3.5      GS00              0.95        Rms    Peak  to peak  amplitude  x4km  χh-dlnVs    l = 2  l = 2, 4  l = 0–4  l = 2  l = 2, 4  l = 0–4    l = 0–4  Run  T1  1.3  1.3  1.6  5.1  6.6  8.5  0.94  −0.86  T2  0.7  1.9  2.1  2.0  7.1  8.7  0.89  −0.94  T3  0.8  0.9  1.2  2.7  4.1  5.8  0.95  −0.54  T1-pPv  1.3  1.3  1.5  4.5  5.9  7.7  0.91  −0.83  T2-pPv  0.8  2.1  2.3  3.4  9.8  11.5  0.81  −0.98  TC1  0.7  0.8  1.0  2.3  3.5  4.4  0.88  0.71  TC2  0.6  0.8  1.0  2.3  3.1  4.5  0.89  −0.13  TC3  0.7  1.2  1.3  2.9  5.7  6.7  0.92  0.33  TC4  1.1  1.2  1.6  4.4  5.4  7.6  0.89  −0.68  TC5  1.3  1.6  2.5  5.2  7.5  9.8  0.82  0.93  TC6  2.5  2.3  2.2  7.3  9.1  8.2  0.75  0.88  TC7  0.2  0.6  0.7  0.9  2.7  3.5  0.95  0.27  TC8  0.4  0.5  0.9  1.2  2.9  4.0  0.93  0.58  TC9  0.5  0.7  1.0  1.7  3.0  4.0  0.90  0.82  TC10  0.4  0.7  0.9  1.8  3.4  5.1  0.93  0.58  TC11  1.0  1.2  1.5  3.4  5.8  6.8  0.92  0.65  TC1-pPv  1.1  1.3  1.6  4.6  5.2  7.1  0.87  0.54  TC2-pPv  0.4  0.8  1.5  1.5  3.9  6.8  0.72  −0.29  TC4-pPv  1.1  1.2  1.5  4.8  5.6  8.6  0.89  −0.50  TC6-pPv  2.3  2.4  2.9  9.1  11.5  14.9  0.70  0.89  TC7-pPv  0.6  0.9  1.2  2.2  4.3  6.0  0.91  0.05  Observed  KDT12        ≤5.0          T10          5.0        MD87            12.0      DH89            8.0      SV03            4.0      KPF03            3.5      GS00              0.95    View Large Constraints on peak-to-peak amplitude may further be inferred from seismic phases reflecting either on the core side or on the mantle side. Using PKKP and PcP phases, Garcia & Souriau (2000) showed that 95 per cent of CMB topography with wavelength larger than 300 km should be between −4.0 and 4.0 km. Our thermal and thermochemical models are in good agreement with this observation (Figs 6 and 10; Table 3). The fraction of CMB area with topography ranging between −4.0 and 4.0 km calculated from unfiltered models is higher than 85 per cent except, again, for of thermochemical models with very large density contrast (Table 3). Importantly, the observations of Garcia & Souriau (2000) do not rule out the presence of deep troughs associated with downwellings since, according to our models, such features would be very localized and would cover a small fraction of the CMB area. Overall, the peak-to-peak amplitude of long-wavelength dynamic topography is unlikely to provide a clear diagnosis on the nature, thermal or thermochemical, of the deep mantle. On average, thermochemical models with density contrast in the range 100–200 kg m−3 have peak-to-peak amplitude lower than other models, but peak-to-peak amplitude further depends on other parameters, in particular to the sensitivity of viscosity to temperature. 4.2 Pattern: correlation between long-wavelength topography and shear wave velocity Lassak et al. (2010) showed that CMB dynamic topography induced by purely thermal and thermochemical models of convection strongly differs at short wavelengths, but are not substantially different at long wavelengths, that is, that current seismic maps of dynamic topography, which are limited to spherical harmonic degree l = 4, would be unable to discriminate between purely thermal and thermochemical structures. They further pointed out that long-wavelength topography and shear wave velocity would be anticorrelated. This conclusion may however change, depending on details of the thermochemical model. Left and middle columns in Fig. 13 compares the long-wavelength (spherical harmonic filter l = 0–4) distributions of lowermost mantle shear wave velocity anomalies and CMB dynamic topography for models T1, TC1, TC4, and TC7 (Table 2). For each model, we calculated shear wave velocity anomalies from the thermal and compositional distributions they predict, using the approach and seismic sensitivities described in Deschamps et al. (2012). We then averaged out these anomalies according to the radial parametrization of Masters et al. (2000), in which the lowermost layer samples the depth range 2700–2891 km. Non-dimensional temperature is rescaled with the super-adiabatic temperature jump ΔTS = 2500 K, leading to rms temperature anomalies typically around −400 K and 550 K in plume and slab regions, respectively, and 450 K within thermochemical reservoirs. For thermochemical models, we assumed that the dense material is enriched in iron and bridgmanite by 3 per cent and 18 per cent compared to regular (pyrolitic) material. For purely thermal model T1, plumes cause slower than average seismic velocity and positive topography. As a result, long-wavelength seismic velocity anomalies and CMB dynamic topography are strongly anticorrelated (Figs 13a and b). In thermochemical models with small density contrast, typically ∼90 kg m−3 and less, as model TC4 (Bz = 0.15, corresponding to ΔρC = 93 kg m−3) and the models of Lassak et al. (2010), thermochemical piles induce shallow depressions (Figs 9a–d). However, when filtered to long-wavelengths, topography beneath piles is positive, and long-wavelengths shear wave velocity anomalies and CMB topography are, again, strongly anticorrelated (Figs 13g and h). CMB topography may thus be unable to discriminate between a purely thermal and a weakly (ΔρC < 100 kg m−3) thermochemical structure of the mantle. However, for larger density contrasts, in particular model TC1 (Bz = 0.23, corresponding to ΔρC = 142 kg m−3), reservoirs of primordial material induce slower than average shear wave velocity (because they are hotter than average and enriched in iron), and cause a depression in the CMB (Fig. 5). Shear wave velocity anomalies and CMB topography are then strongly correlated (Figs 13d and e). For model TC7, in which intrinsic viscosities of primordial and regular material are similar (ΔηC = 1), the depression caused by reservoirs of primordial material is less pronounced (Figs 11g–h), resulting in a poorer, but still significant correlation between the maps of shear wave velocity anomalies and CMB dynamic topography (Figs 13j and k). Our models are not expected to strictly reproduce observed seismic patterns above the CMB. It is however interesting to compare long-wavelength shear wave velocity anomalies predicted by models of convection (Fig. 13, left column) with observed tomography from model SB10L18 (Masters et al.2000) at z = 2790 km depth (Fig. 14a). For model T1, shear wave velocity anomalies averaged in the layer 2700–2891 km vary in the range –0.7 to 0.7 per cent with rms amplitude around 0.3 per cent. This is much smaller than SB10L18, which varies between –2.5 and 1.4 per cent, with rms amplitude of 0.9 per cent. Explaining amplitudes of shear wave velocity anomalies observed in SB10L18 would require much larger temperature anomalies, around ± 1200 K. In thermochemical models, by contrast, shear-velocity anomalies have rms amplitude comparable to that of SB10L18. For model TC1, for instance, they vary between –2.1 and 1.7 per cent with rms amplitude around 1.1 per cent. Shear-velocity anomalies patterns predicted by models TC1 and TC4 are dominated by spherical harmonic degrees 2 and 3, and their power spectra are overall in good agreement with that of SB10L18 (Fig. 14b). The pattern predicted by model TC7 is more complex, with a strong contribution of spherical harmonic degree 4 that is not seen in SB10L18. Figure 13. View largeDownload slide Shear wave velocity anomalies in the lowermost mantle (2700–2891 km) (left column), CMB dynamic topography (middle column) and CMB isostatic topography (right column) predicted by purely thermal model T1 (Fig. 4) and thermochemical models TC1 (Fig. 5), TC4 (Fig. 9a-d), and TC7 (Fig. 11e-h), filtered for spherical harmonic degrees up to l = 4. Model properties are listed in Table 2. Shear wave velocity anomalies are calculated from the thermal and compositional fields predicted by each model using the seismic sensitivities of Deschamps et al. (2012), and are averaged out in the depth range 2700–2891 km, corresponding to the lowermost layer of Masters et al. (2000) tomographic model. For thermochemical models, primordial material is assumed to be enriched in iron and bridgmanite by 3 per cent and 18 per cent, respectively, compared to regular (pyrolitic) material. Isostatic topography is calculated from density anomalies averaged in the lowermost 400 km of the mantle (2591–2891 km). These density anomalies are directly estimated from the thermal and compositional fields predicted by each model, without using seismic sensitivity. Relative density anomalies, not shown here, are fully correlated with isostatic topography, with scaling factor between isostatic topography (in km) and relative density anomalies (in per cent), Rh, indicated at the bottom right of each plot. Interval of contour levels are 0.5 per cent for dlnVS and 1 km for CMB topography. Root mean square (rms) of each distribution is indicated at the bottom left of each plot. Figure 13. View largeDownload slide Shear wave velocity anomalies in the lowermost mantle (2700–2891 km) (left column), CMB dynamic topography (middle column) and CMB isostatic topography (right column) predicted by purely thermal model T1 (Fig. 4) and thermochemical models TC1 (Fig. 5), TC4 (Fig. 9a-d), and TC7 (Fig. 11e-h), filtered for spherical harmonic degrees up to l = 4. Model properties are listed in Table 2. Shear wave velocity anomalies are calculated from the thermal and compositional fields predicted by each model using the seismic sensitivities of Deschamps et al. (2012), and are averaged out in the depth range 2700–2891 km, corresponding to the lowermost layer of Masters et al. (2000) tomographic model. For thermochemical models, primordial material is assumed to be enriched in iron and bridgmanite by 3 per cent and 18 per cent, respectively, compared to regular (pyrolitic) material. Isostatic topography is calculated from density anomalies averaged in the lowermost 400 km of the mantle (2591–2891 km). These density anomalies are directly estimated from the thermal and compositional fields predicted by each model, without using seismic sensitivity. Relative density anomalies, not shown here, are fully correlated with isostatic topography, with scaling factor between isostatic topography (in km) and relative density anomalies (in per cent), Rh, indicated at the bottom right of each plot. Interval of contour levels are 0.5 per cent for dlnVS and 1 km for CMB topography. Root mean square (rms) of each distribution is indicated at the bottom left of each plot. Figure 14. View largeDownload slide Long-wavelength shear wave velocity anomalies in the lowermost mantle. (a) Observed tomography (SB10L18, Masters et al. 2000) filtered for spherical harmonic degrees up to l = 4. Interval of contour levels is 0.5 per cent. (b) Power spectra of SB10L18 (histogram) and selected models shown in Fig. 13 (colour code) up to spherical harmonic degree l = 8. Figure 14. View largeDownload slide Long-wavelength shear wave velocity anomalies in the lowermost mantle. (a) Observed tomography (SB10L18, Masters et al. 2000) filtered for spherical harmonic degrees up to l = 4. Interval of contour levels is 0.5 per cent. (b) Power spectra of SB10L18 (histogram) and selected models shown in Fig. 13 (colour code) up to spherical harmonic degree l = 8. To quantitatively compare long-wavelength patterns of topography and shear wave velocity predicted by models of convection, we calculated the correlation χh-dlnVs between these distributions (Table 3 and Fig. 12b). Clearly, CMB topography and shear-velocity anomalies for purely thermal models and thermochemical models with small density contrast anticorrelate, with χh-dlnVs ≤ –0.5, while they correlate for most thermochemical models with moderate and strong density contrast, with χh-dlnVs ≥ 0.5. Interesting exceptions are cases with moderate thermal viscosity contrast (ΔηT = 106; TC2 and TC2-pPv), for which CMB topography and shear-velocity anomalies are de-correlated or slightly anticorrelated, and, again, models with no chemical viscosity contrast (ΔηC = 1; TC7), for which distributions of topography and shear-velocity are only slightly correlated. These results suggest that long-wavelength CMB dynamic topography potentially provides a test to discriminate between purely thermal models of lower mantle and thermochemical models in which LLSVPs consist of material substantially denser (ΔρC > 100 kg m−3, i.e. ∼1.8 per cent of PREM and more) and intrinsically more viscous than the surrounding mantle. Interestingly, Tanaka (2010) found depressions of ∼2 km deep beneath Africa and the Pacific, that is, fairly correlated with LLSVPs. If correct, and following our models, this would support the hypothesis that LLSVPs are thermochemical structures denser than the surrounding mantle, with density contrast typically in the range 100–200 kg m−3. This scenario is further supported by available estimates of chemical density anomalies (Trampert et al.2004; Mosca et al.2012), which have a peak-to-peak amplitude around 2.0–3.0 per cent. In addition, because primordial material may be enriched in bridgmanite (Trampert et al.2004; Mosca et al.2012), which is more viscous than ferro-periclase by up to 3 orders of magnitude (Yamazaki & Karato 2001), LLSVPs may be intrinsically more viscous than the surrounding mantle. Other models of CMB topography, however, do not observe clear correlation or anticorrelation between shear-velocity and CMB topography patterns. For instance, the topography observed by Doornbos & Hilton (1989) is dominated by spherical harmonic degree 1 with a depression beneath India and a topographic high beneath Western Europe and Africa, while the pattern obtained by Morelli & Dziewonski (1987) is more complex, with topographic highs beneath South Indian Ocean, Northern Atlantic and Northern Pacific. Soldati et al. (2012) performed inversion of P-wave data and reflected seismic phases (PcP and PKP) jointly for compressional velocity anomalies and CMB topography, using a geodynamically constrained method. More precisely, density anomalies were scaled with seismic velocity anomalies and combined with CMB topography kernels calculated using a spectral method (e.g. Forte & Peltier 1991). Resulting CMB topography maps have peak-to-peak amplitude around 8 km, and are fairly anticorrelated with compressional velocity maps above the CMB. In particular, LLSVPs are associated with topography highs with strong positive topography of up to about 5 km, beneath the central Pacific and southern tip of Africa. It has further been pointed out that CMB topography models inferred from geodynamically constrained tomography fit normal mode constraints better than models deduced directly from seismic tomography (Soldati et al.2013). It should be noted, however, that the approach used by Soldati et al. (2012) implicitly assumes that chemical anomalies are small, and that they are, at a given depth, fully correlated with seismic velocity anomalies. Furthermore it does not account for lateral variations in viscosity (e.g. related to variations in temperature), which, as illustrated by our calculations and those of Lassak et al. (2010), may substantially affect the amplitude of topography. 5 DISCUSSION Our calculations show strong differences between the CMB dynamic topographies induced by purely thermal and thermochemical models of convection. In purely thermal models plumes elevate the CMB, while hot thermochemical piles depress it. In both cases, increasing the thermal viscosity contrast decreases the amplitude of CMB topography. In thermochemical models, the buoyancy ratio, measuring the density contrast between primordial and regular material, and the thermal and compositional viscosity contrasts further affect the shape and stability of reservoirs of primordial material, and the amplitude of dynamic topography, as sketched in plots c and d of Fig. 1. Increasing the chemical density contrast and viscosity ratio enhances the depression beneath reservoirs. 5.1 Influence of post-perovskite Our calculations further indicate that the presence of pPv at the bottom of the mantle does not substantially influence the dynamic topography. All other parameters being the same, the presence of pPv slightly modifies the amplitude of topography, but the effects due to varying one specific parameter are unchanged. This conclusion assumes that the viscosity of pPv is similar or close to that of bridgmanite. If, in contrast, and as suggested by ab initio calculations (Ammann et al.2010), pPv is less viscous than bridgmanite by two to three orders of magnitude, dynamic topography may be strongly altered. Both the numerical models of Yoshida (2008) and the CMB maps of Soldati et al. (2012) report a strong attenuation of CMB topography when a low viscosity layer is imposed in the lowermost 200 km of the mantle completely surrounding the core. In our models, the pPv stability region does not cover the entire core surface, but forms discontinuous patches located outside plume clusters and, in the thermochemical cases, outside reservoirs of primordial material. In these regions, and if pPv is much less viscous than bridgmanite, CMB topography may thus flatten out compared to that found in our models. 5.2 Comparison with previous models Our models are in overall agreement with those of Lassak et al. (2010). In particular, the influences of buoyancy ratio and thermal viscosity contrast on CMB topography they observed are similar to what we see in our simulations. Note, however, that due to different choices of the super-adiabatic temperature jump and reference thermal expansion (Table 1), and of a different definition of the buoyancy ratio (eq. 1), the chemical density contrast at the bottom of the system for a given value of buoyancy ratio is much larger in our models. For instance, in the models of Lassak et al. (2010) with B = 0.6 it is around 60 kg m−3 (1 per cent of PREM), whereas in our models with Bz = 0.23 it reaches 142 kg m−3 (2.5 per cent of PREM). As a result, in model TC7 thermochemical reservoirs induce a 0.7 km deep depression in dynamic topography (Figs 8e–h), and a plateau with low elevation in Lassak et al. (2010) models. In model TC1 (Fig. 5), this trend is amplified by the fact that primordial material is intrinsically more viscous than the regular material by a factor 30, which deepens the depression beneath reservoirs of primordial material. To investigate the influence of LLSVPs on dynamic topography, Yoshida (2008) ran models of convection in which density anomalies are scaled from seismic tomography model SMEAN (Becker & Boschi 2002), and arrived at similar conclusions regarding the effect of the chemical density contrast. In purely thermal models, LLSVPs result in positive topography of ∼5 km. Models assuming piles denser by 1 per cent lead to a relatively flat positive topography beneath LLSVPs, while a larger density contrast (2 per cent) causes relatively flat topography with local depressions ∼2 km deep, not unlike our model TC7. Interestingly, this negative topography flattens out when the viscosity of the lowermost mantle is decreased by one or two orders of magnitude. Note, however, that Yoshida (2008) did not run models in which LLSVP viscosity is intrinsically different from that of surrounding material. 5.3 Predicted isostatic topography A recent analysis of Stoneley mode splitting data, which are sensitive to the lowermost mantle structure, indicates that LLSVPs might be lighter than surrounding mantle (Koelemeijer et al.2017), in contradiction with earlier normal mode studies (Ishii & Tromp 1999; Trampert et al.2004; Mosca et al.2012). Stoneley mode are further sensitive to CMB topography, and density estimates thus trade-off with this topography. Assuming that long-wavelength (>1000 km) CMB topography is due to isostatic compensation, Koelemeijer et al. (2017) found that, for degree 2, their data set is best fitted if LLSVPs are lighter than PREM by about 1.4 per cent and CMB is elevated by about 1 km. For degrees 4 and 6, best fits are, by contrast, obtained for LLSVPs denser than average mantle by about 2.5 per cent and a CMB depressed by 1–5 km. When all structural degrees are taken into account, light LLSVPs with CMB topography around 1 km or more are again preferred, but the overall fit is much reduced. In addition, when accounting for seismic normal modes used in previous studies, density models do not have preference for denser or lighter LLSVPs. Finally, it should be kept in mind that splitting data are based on the self-coupling approximation (Woodhouse & Girnius 1982), which might bias the inference of density structure (Al-Attar et al.2012; Yang & Tromp 2015). Lateral variations in temperature and composition predicted by our models can be used to estimate lateral density anomalies and reconstruct CMB isostatic topography. For purely thermal models, one expects that CMB dynamic and isostatic topographies have similar trends. On the contrary, because thermochemical piles and reservoirs are hotter than average, patterns of dynamic and isostatic topographies may differ substantially, in particular for low buoyancy ratio. Right column in Fig. 13 plots CMB isostatic topography estimated from density anomalies averaged in the lowermost 400 km of the mantle for models T1, TC1, TC4 and TC7. Density anomalies, not shown here, are, by definition, fully correlated with isostatic topography, with scaling factor between topography (in km) and relative density anomalies (in per cent) around −0.27 (the exact value of this scaling factor depends on the average density in the lowermost 400 km layer, which is very similar for all 4 models shown here). For the calculation of density anomalies, temperature anomalies are rescaled with the super adiabatic temperature jump, and thermal expansion and density are rescaled with their surface values (Table 1). Reference density and thermal expansion vary as indicated in Fig. 2, and density anomalies are calculated with respect to horizontal average (which is different and usually higher than the reference density). For comparison with the dynamic topography shown in the middle column of Fig. 13, we filtered density anomalies and isostatic topography for spherical harmonic degrees l = 0–4. For purely thermal model T1, isostatic topography has a pattern similar to that of dynamic topography, but its amplitude is smaller by about a factor 4. Within plume regions, the average filtered density anomaly and isostatic topography are −0.1 per cent and 0.4 km (Fig. 13c), respectively, and are thus consistent with Stoneley mode splitting measurements of Koelmeijer et al. (2017). In thermochemical models, thermal effects alter the effective density within reservoirs and piles of dense material. For small density contrast, as in model TC4 (Fig. 13i), thermal effects nearly compensate chemical effects within piles, leading to average density anomaly and isostatic topography of about 0.3 per cent and −1.0 km, respectively. For larger density contrast, as in models TC1 and TC7, density anomalies within piles are still dominated by chemical effects, resulting in relatively large isostatic topographies (Figs 13f and l). For instance, in model TC1, density anomalies within thermochemical reservoirs and isostatic topography beneath them are equal to 0.8 per cent and −3.2 km on average. Such reservoirs provide poorer fit to Koelemeijer et al. (2017) observations than light elevated plume clusters in purely thermal models, but are still in reasonable agreement with these data (see their fig. 3 for RLL = −0.6 and H = −5, corresponding to LLSVP total density anomaly of 0.8 per cent and isostatic topography of −2.5 km). Isostatic topography at CMB may be considered as a first order estimate of CMB topography, taking into account only the effects of hydrostatic pressure. Dynamic topography further accounts for stresses related to mantle flow, in particular slab pushing or plume pulling. This appears very clearly when comparing middle and right columns in Fig. 13. For instance, in case T1, hot, lighter plumes induce positive isostatic topography. They further pull the CMB upwards, thus increasing CMB topography. As a result, compared to isostatic topography (Fig. 13c), dynamic topography has a similar pattern, but more amplitude (Fig. 13b). For thermochemical cases, dynamical and isostatic effects are opposing each other. Dynamic and isostatic topographies have slightly different patterns, and depressions beneath thermochemical reservoirs are shallower in dynamic topography (compare, for instance, Figs 13e and f for case TC1). Note that for small chemical density contrasts, as in case TC4, dynamic effects lead to very shallow depressions (Fig. 9c), and when filtered for long-wavelengths, dynamic topography beneath thermochemical piles is positive (Fig. 13h). Because it accounts for dynamics effects related to mantle flow, we expect dynamic topography to provide a better description of CMB topography than isostatic topography does. Interestingly, in dynamic topography predicted by case TC1 and filtered with l = 0–4, the depth of the depressions beneath thermochemical reservoirs is on average −1.0 km (Fig. 13e), which still reasonably fits Stoneley mode observations. 5.4 CMB excess ellipticity An important constraint on dynamic topography is provided by the excess ellipticity of the CMB, that is, its spherical harmonic degree l = 2 and order m = 0 term (Y20). This excess can be deduced from Very Long Baseline Interferometry (VLBI) measurements of the retrograde annual nutation of the Earth's rotation axis. More precisely, the period of the free core nutation (FCN) depends on the CMB dynamic flattening, ef, and the resonance between FCN and retrograde annual nutation increases the amplitude of this annual nutation. The FCN period deduced from VLBI data, first measured by Gwinn et al. (1986), indicates that the CMB dynamical flattening is larger than the value expected for a hydrostatic Earth, with an excess dynamical flattening of 5 per cent. Assuming that this flattening is due to a homogeneous non-hydrostatic layer at the top of the outer core, it can be interpreted as an excess difference between the equatorial and polar radii of the CMB of 490 ± 110 m. Taking into account more complexities, including ocean tides, electromagnetic coupling between core and mantle, and anelastic effects in the mantle, Mathews et al. (2002) found a slightly smaller excess, of ∼390 m, corresponding to an excess CMB dynamical flattening of 3.8 per cent. This excess of CMB ellipticity is usually attributed to mantle dynamics (Gwinn et al.1986; Forte et al.1995; Lumb et al.1995). Our models do not define a polar axis and associated equatorial plane that would allow the calculation of excess CMB ellipticity. Nonetheless, one may point out that excess CMB ellipticity would be difficult to reconcile with large, deep depressions in CMB dynamic topography, as those predicted in model TC1 (Bz = 0.23 and ΔηC = 30), located along the equator. This problem may be solved in models with a slightly lower buoyancy ratio, around Bz = 0.2 (ΔρC = 124 kg m−3), which is still consistent with density anomalies observed in the lowermost mantle (Trampert et al.2004; Mosca et al.2012), and/or lower chemical density contrast (typically, ΔηC around 10 and less). In these models, the depressions associated with reservoirs of primordial material are much attenuated, as illustrated by model TC7, and would therefore have a limited influence on excess CMB ellipticity. Interestingly, such models still predict long-wavelength topographies that are distinct from those obtained for purely thermal convection (Fig. 13), with, for models with ΔηC around 10 a correlation with shear wave velocity anomalies that is still important (around 0.6 for model TC8). In addition, the link between dynamic topography and observed excess ellipticity may be affected by other factors. First, it should be kept in mind that the period of the FCN constrains the dynamic flattening, and that interpreting this flattening in terms of geometric flattening requires the prescription of a non-hydrostatic layer at the top of the outer core, with thickness and density as free parameters. Moreover, it has been pointed out that other components of CMB topography than its Y20 term may influence the period of the FCN (Wu & Wahr 1997), and therefore the VLBI measurements and the estimates of the excess CMB ellipticity inferred from them. Finally, gravitational and viscous couplings between specific mantle structures (e.g. LLSVPs) and the core may also affect the period of the FCN. For instance, a recent study (Chao 2017) showed that the 6 yr periodic variations of the length of the day are well explained by gravitational coupling between the inner core and the lowermost mantle. Interestingly, this coupling implies mass excess in regions overlapping LLSVPs (Ding & Chao, submitted), thus contradicting the Stoneley mode analysis of Koelemeijer et al. (2017), but supporting earlier normal mode studies (Ishii & Tromp 1999; Trampert et al.2004; Mosca et al.2012). The effect of this coupling on the period of the FCN, however, remains to be investigated. 5.5 Implications for core dynamics Another interesting question that remains to be investigated is the effect of dynamic topography on outer core flow. Calkins et al. (2012) run core flow models in cylindrical geometry with ridges (i.e. trenches on the mantle side) of Gaussian shape, and found that this topography results in stationary Rossby waves, which, in turn, increase the strength of the zonal flow and induce flow structures that remain fixed relative to the mantle. More recently, Tarduno et al. (2015) suggested that a step-like topography in CMB due to the African LLSVP could induce upwellings at the top of the core beneath the south of Africa, resulting in reversed magnetic flux at the CMB and persistent low magnetic field at the surface of the Earth. The gradient topography predicted by our models are rather modest, as illustrated in Fig. 15 (note that the scale of vertical axis is exaggerated by about a factor 115 compared to that of the horizontal axis). For instance, the positive topography induced by plume cluster in model T1 gently increases up to its culminating point (Fig. 15a). The depressions induced by thermochemical reservoirs in model TC1 are flat with marked slopes on their edges, but the gradient of these slopes are typically around 2–4 m km−1 (Fig. 15c). In models with small chemical density contrast, the depressions induced by piles are attenuated and much shallower than the negative topography beneath slabs (Fig. 15d). Slabs induce deeper and steeper depressions on more localised areas, but again, topographic gradients remain limited, up to 15 m km−1 (Fig. 15b). Thus, if CMB topography of real Earth is well described by the dynamic topography predicted by our models, it may have a limited influence on the outer core flow. Figure 15. View largeDownload slide Topographic profiles along different great circles for (a,b) model T1 (Fig. 4), (c) model TC1 (Fig. 5), and (d) model TC4 (Figs 9a–d). Note that the scale of the vertical axis is exaggerated by about a factor 115 compared to that of the horizontal axis, that is, profiles appear much steeper than in real data. Figure 15. View largeDownload slide Topographic profiles along different great circles for (a,b) model T1 (Fig. 4), (c) model TC1 (Fig. 5), and (d) model TC4 (Figs 9a–d). Note that the scale of the vertical axis is exaggerated by about a factor 115 compared to that of the horizontal axis, that is, profiles appear much steeper than in real data. 6 CONCLUSIONS We calculated the dynamic topography at CMB induced by different models of mantle convection, and investigated the influence of important input parameters of these models. In all models, cold downwellings, modelling slabs, induce very deep, but local depressions. In purely thermal models, plumes or plume clusters induce positive topography that decreases in amplitude as the thermal viscosity ratio increases. For moderate buoyancy ratios, corresponding to compositional density contrast (ΔρC) in the range 100–200 kg m−3, consistent with estimates of chemical density anomalies in the lowermost mantle (Trampert et al.2004; Mosca et al.2012), large, stable reservoirs of primordial material are generated. In their interiors, these reservoirs cause depressions in the CMB topography, surrounded by a ridge at their edges. The depth of the depressions and the height of the ridges increase with increasing ΔρC and chemical viscosity contrast (ΔηC). The presence of the pPv phase at the bottom of the mantle does not alter these conclusions, provided that the viscosity of pPv is similar to that of bridgmanite. If, in contrast, pPv viscosity is lower by 2 orders of magnitude or more (Amman et al.2010), the CMB topography may substantially flatten out, as suggested by models in which a weak (low viscosity) layer is imposed at the bottom of the mantle (Yoshida 2008; Soldati et al.2012). Additional calculations are, however, needed to determine the influence of discontinuous patches of weak pPv, as those obtained in models T1-pPv (Figs 7c and d) and TC1-pPv (Figs 8e and f). The peak-to-peak amplitude of CMB dynamic topography predicted by our models for different spherical harmonic filters is comparable to that estimated by available seismic constraints on CMB topography. Purely thermal models usually have larger amplitude than thermochemical models with similar properties, and for spherical harmonic filters including degree l = 2 and degrees l = 2 and l = 4, the amplitude predicted by purely thermal models is overall slightly larger than the maximum amplitude estimated from seismic data. Because it is sensitive to other model parameters, peak-to-peak amplitude alone is however unlikely to discriminate between purely thermal and thermochemical structures of lower mantle. On another hand, a comparison between long-wavelength patterns of shear wave velocity anomalies and topography may provide important clues. The positive topography induced by plumes, when filtered to long wavelengths, anticorrelates with shear wave velocity anomalies induced by the thermal field. By contrast, for thermochemical models, long-wavelength CMB topography correlates well with seismic velocity anomalies for models with ΔρC around 140 kg m−3 or more and intrinsically more viscous primordial material (ΔηC > 1). Shear wave velocity and CMB topography maps progressively decorrelate with decreasing ΔρC or decreasing ΔηC, and for ΔρC around 90 kg m−3 and less, CMB topography and shear wave velocity anticorrelates, as for purely thermal cases. Long-wavelength pattern of CMB topography, provided they can be measured accurately, therefore provides a potential constraint to discriminate between different possible thermochemical structures of lowermost mantle. Acknowledgements We are grateful to Nicolas Flament and to the editor (Gaël Choblet) for their constructive and careful reviews that helped improving a first version of this article. We are also grateful to Hagay Amit and Andrew Valentine for additional discussions on core–mantle coupling and density mapping from normal mode data. This research was funded by the Ministry of Science and Technology (MOST) of Taiwan grant 105-2116-M-001-017, Academia Sinica (Taipei) grant AS-102-CDA-M02, and Orchid grant 103PMRIA0100024. All data and models used in this study are available upon request to FD (frederic@earth.sinica.edu.tw). REFERENCES Al-Attar D., Woodhouse J.H., Deuss A., 2012. Calculation of normal mode spectra in laterally heterogeneous earth models using an iterative direct solution method, Geophys. J. Int. , 189, 1038– 1046. Google Scholar CrossRef Search ADS   Amman M.W., Brodholt J.P., Wookey J., Dobson D.P., 2010. First principles constraints on diffusion in lower-mantle minerals and a weak D″ layer, Nature , 465, 462– 465. Google Scholar CrossRef Search ADS PubMed  Becker T.W., Boschi L., 2002. A comparison of tomographic and geodynamic mantle models, Geochem. Geophys. Geosyst. , 3, doi:10.1029/2001GC000168. Calkins M.A., Noir J., Eldredge J.D., Aurnou J.M., 2012. The effects of boundary topography on convection in Earth's core, Geophys. J. Int. , 189, 799– 814. Google Scholar CrossRef Search ADS   Chao B.F., 2017. Dynamics of axial torsional libration under the mantle-inner core gravitational interaction, J. geophys. Res. , 121, doi:10.1002/2016JB013515. Christensen U.R., Yuen D.A., 1985. Layered convection induced by phase transitions, J. geophys. Res. , 90( 10), 291– 210, 300. Crameri F., Tackley P.J., 2014. Spontaneous development of arcuate single-sided subduction in global 3-D mantle convection models with a free surface, J. geophys. Res. , 119, 5921– 5942. Google Scholar CrossRef Search ADS   Crameri F.et al.  , 2012. comparison of numerical surface topography calculations in geodynamic modelling: an evaluation of the ‘sticky air’ method, Geophys. J. Int. , 189, 38– 54. Google Scholar CrossRef Search ADS   Davies D.R., Goes S., Davies J.H., Schuberth B.S.A., Bunge H.-P., Ritsema J., 2012. Reconciling dynamic and seismic models of Earth's lower mantle: The dominant role of thermal heterogeneity, Earth planet. Sci. Lett. , 353–354, 253– 269. Google Scholar CrossRef Search ADS   Davies D.R., Goes S., Lau H.C.P., 2015. Thermally dominated deep mantle LLSVPs: a review, in The Earth's Heterogeneous Mantle , pp. 441– 478, eds Khan A., Deschamps F., Springer. Google Scholar CrossRef Search ADS   Davies J.H., 2013. Global map of solid Earth heat flow, Geochem. Geophys. Geosyst. , 14, doi:10.1029/ggge.20271. Deschamps F., Tackley P.J., 2009. Searching for models of thermo-chemical convection that explain probabilistic tomography II - Influence of physical and compositional parameters, Phys. Earth planet. Inter. , 176, 1– 18. Google Scholar CrossRef Search ADS   Deschamps F., Cobden L., Tackley P.J., 2012. The primitive nature of large low shear-wave velocity provinces, Earth planet. Sci. Lett. , 349–350, 198– 208. Google Scholar CrossRef Search ADS   Deschamps F., Li Y., Tackley P.J., 2015. Large-scale thermo-chemical structure of the deep mantle: observations and models, in The Earth's Heterogeneous Mantle , pp. 479– 515, eds Khan A., Deschamps F., Springer. Google Scholar CrossRef Search ADS   Deschamps F., Khan A., 2016. Electrical conductivity as a constraint on lower mantle thermo-chemical structure, Earth planet. Sci. Lett. , 450, 108– 119. Google Scholar CrossRef Search ADS   Ding H., Chao B.-F., A 6-year westward rotary motion in the Earth: detection and possible MICG coupling mechanism, submitted to Earth Planet. Sci. Lett.  Doornbos D.J., Hilton T., 1989. Models of the core-mantle boundary and the travel times of internally reflected core phase, J. geophys. Res. , 94, 15 741–15 751. Google Scholar CrossRef Search ADS   Dziewonski A.M., Anderson D.L., 1981. Preliminary reference Earth model, Phys. Earth Planet. Inter. , 25, 297– 356. Google Scholar CrossRef Search ADS   Forte A.M., Peltier R., 1991. Viscous flow models of global geophysical observables 1. Forward problems, J. geophys. Res. , 96, 20 131–20 159. Google Scholar CrossRef Search ADS   Forte A.M., Mitrovica J.X., Woodward R.L., 1995. Seismic-geodynamic determination of the origin of excess ellipticity of the core–mantle boundary, Geophys. Res. Lett. , 22, 1013– 1016. Google Scholar CrossRef Search ADS   Garcia R., Souriau A., 2000. Amplitude of the core–mantle boundary topography estimated by stochastic analysis of core phases, Phys. Erth planet. Inter. , 117, 345– 359. Google Scholar CrossRef Search ADS   Garnero E.J., McNamara A., Shim S.-H., 2016. Continent-sized anomalous zones with low seismic velocity at the base of Earth's mantle, Nat. Geosci. , 9, 481– 489. Google Scholar CrossRef Search ADS   Gwinn C.R., Herring T.A., Shapiro I.I., 1986. Geodesy by radio interferometry: studies of the forced nutations of the Earth 2. Interpretation, J. geophys. Res. , 91, 4755– 4765. Google Scholar CrossRef Search ADS   Hager B.H., Richards M.A., 1989. Long wavelength variations in Earth's geoid: physical models and dynamical implications, Phil. Trans. R. Soc. A , 328, 309– 327. Google Scholar CrossRef Search ADS   Hansen U., Yuen D.A., 1989. Dynamical influences from thermo-chemical instabilities at the core–mantle boundary, Geophys. Res. Lett. , 16, 629– 632. Google Scholar CrossRef Search ADS   Hernlund J., Houser C., 2008. On the statistical distribution of seismic velocities in Earth's deep mantle, Earth planet. Sci. Lett. , 265, 423– 437. Google Scholar CrossRef Search ADS   Ishii M., Tromp J., 1999. Normal-mode and free-air gravity constraints on lateral variations in velocity and density of Earth's mantle, Science , 285, 1231– 1236. Google Scholar CrossRef Search ADS PubMed  Jaupart C., Mareschal J.-C., 2007. Heat flow and thermal structure of the lithosphere, in Treatise on Geophysics , 6, pp 218– 251, ed. Schubert, G., Elsevier. Jaupart C., Labrosse S., Lucazeau F., Mareschal J.-C., 2015. Temperature, heat, and energy in the mantle of the Earth, in Treatise on Geophysics , 2nd edn, 7, pp. 218– 251, ed. Schubert G., Elsevier. Kageyama A., Sato T., 2004. “Yin-Yang grid”: an overset grid in spherical geometry, Geochem. Geophys. Geosys. , 5, doi:10.1029/2004GC000734. Koelemeijer P.J., Deuss A., Trampert J., 2012. Normal mode sensitivity to Earth's D″ layer and topography on the core–mantle boundary: what we can and cannot see, Geophys. J. Int. , 190, 553– 568. Google Scholar CrossRef Search ADS   Koelemeijer P., Deuss A., Ritsema J., 2017. Density structure of Earth's lowermost mantle from Stoneley mode splitting observations, Nat. Commun. , 8, doi:10.1038/ncomms15241. Konishi K., Fuji N., Deschamps F., 2017. Elastic and anelastic structure of the lowermost mantle beneath the western Pacific from waveform inversion, Geophys. J. Int. , 208, 1290– 1304. Google Scholar CrossRef Search ADS   Korenaga J., Karato S.-I., 2008. A new analysis of experimental data on olivine rheology, J. geophys. Res. , 113, B02403, doi:10.1029/2007JB005100. Google Scholar CrossRef Search ADS   Koper K.D., Pyle M.L., Franks J.M., 2003. Constraints on aspherical core structure from PKiKP-PcP differential travel times, J. geophys. Res. , 108, doi:10.1029/2002JB001995. Lassak T.M., McNamara A.K., Zhong S., 2007. Influence of thermo-chemical piles on topography at Earth's core-mantle boundary, Earth planet. Sci. Lett. , 261, 443– 455. Google Scholar CrossRef Search ADS   Lassak T.M., McNamara A.K., Garnero E.J., Zhong S., 2010. Core–mantle boundary topography as a possible constraint on lower mantle chemistry and dynamics, Earth planet. Sci. Lett. , 289, 232– 241. Google Scholar CrossRef Search ADS   Lekić V., Cottaar S., Dziewonski A.M., Romanowicz B., 2012. Cluster analysis of global lower mantle tomography: A new class of structure and implications for chemical heterogeneity, Earth planet. Sci. Lett. , 357–358, 68– 77. Google Scholar CrossRef Search ADS   Li P., Karato S.-I., Wang Z., 1996. High-temperature creep in fine-grained polycrystalline CaTiO3, an analogue material of (Mg,Fe)SiO3 perovskite, Phys. Earth planet. Inter. , 95, 19– 36. Google Scholar CrossRef Search ADS   Li M., McNamara A.K., Garnero E.J., 2014a. Chemical complexity of hotspots caused by cycling oceanic crust through mantle reservoirs, Nat. Geosci. , 7, 366– 370. Google Scholar CrossRef Search ADS   Li Y., Deschamps F., Tackley P.J., 2014b. The stability and structure of primordial reservoirs in the lower mantle: insights from models of thermo-chemical convection in 3-D spherical geometry, Geophys. J. Int. , 199, 914– 930. Google Scholar CrossRef Search ADS   Li Y., Deschamps F., Tackley P.J., 2015. Effects of the post-perovskite phase transition properties on the stability and structure of primordial reservoirs in the lower mantle of the Earth, Earth planet. Sci. Lett. , 432, 1– 12. Google Scholar CrossRef Search ADS   Lumb L.I., Jarvis G.T., Aldridge K.D., DeLandro-Clarke W., 1995. The period of the free core nutation: towards a dynamical basis for an “extra-fattening” of the core-mantle boundary, Phys. Earth planet. Inter. , 90, 255– 271. Google Scholar CrossRef Search ADS   Masters G., Laske G., Bolton H., Dziewonski A.M., 2000. The relative behavior of shear velocity, bulk sound speed, and compressional velocity in the mantle: implication for thermal and chemical structure, in Earth's Deep Interior: Mineral Physics and Tomography from the Atomic to the Global Scale , pp. 63– 87, Geophysical Monograph Ser. 117, eds Karato,, S.-I. et al., American Geophysical Union. Google Scholar CrossRef Search ADS   Mathews P.M., Herring T.A., Buffett B., 2002. Modeling of nutation and precession: new nutation series for nonrigid Earth and insights into the Earth's interior, J. geophys. Res. , 107, doi:10.1029/2001JB000390. Morelli A., Dziewonski A.M., 1987. Topography of the core-mantle boundary and lateral homogeneity of the liquid core, Nature , 325, 678– 683. Google Scholar CrossRef Search ADS   Mosca I., Cobden L., Deuss A., Ritsema J., Trampert J., 2012. Seismic and mineralogical structures of the lower mantle from probabilistic tomography, J. geophys. Res. , 117, doi:10.1029/2011JB008851. Murakami M., Hirose K., Kawamura K., Sata N., Ohishi Y., 2004. Post-perovskite phase transition in MgSiO3, Science , 304, 855– 858. Google Scholar CrossRef Search ADS PubMed  Ricard Y., Fleitout L., Froideveaux C., 1984. Geoid heights and lithospheric stresses for a dynamic Earth, Ann. Geophys. , 2, 267– 286. Ritsema J., Deuss A., vanHeijst H.-J., Woodhouse J.H., 2011. S40RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements, Geophys. J. Int. , 184, 1223– 1236. Google Scholar CrossRef Search ADS   Schuberth B.S.A., Zaroli C., Nolet G., 2012. Synthetic seismograms for a synthetic Earth: long-period P- and S-wave traveltime variations can be explained by temperature alone, Geophys. J. Int. , 188, 1393– 1412. Google Scholar CrossRef Search ADS   Soldati G., Boschi L., Forte A.M., 2012. Tomography of core–mantle boundary and lowermost mantle coupled by geodynamics, Geophys. J. Int. , 189, 730– 746. Google Scholar CrossRef Search ADS   Soldati G., Koelemeijer P., Boschi L., Deuss A., 2013. Constraints on core–mantle boundary topography from normal mode splitting, Geochem. Geophys. Geosys. , 14, doi:10.1002/ggge.20115. Solomatov V.S., Reese C.C., 2008. Grain size variations in the Earth's mantle and the evolution of primordial chemical heterogeneities, J. geophys. Res. , 113, B07408, doi:10.1029/2007JB005319. Google Scholar CrossRef Search ADS   Sze E.K.M., van der Hilst R.D., 2003. Core mantle boundary topography from short period PcP, PKP, and PKKP data, Phys. Earth planet. Inter. , 135, 27– 46. Google Scholar CrossRef Search ADS   Tackley P.J., 1998. Three-dimensional simulations of mantle convection with a thermo-chemical CMB boundary layer: D″?, in The Core-Mantle Boundary Region , Geodynamical Ser., 28, pp. 231– 253, eds Gurnis M. et al., American Geophysical Union. Tackley P.J., 2008. Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid, Phys. Earth planet. Inter. , 171, 7– 18. Google Scholar CrossRef Search ADS   Tackley P.J., 2012. Dynamics and evolution of the deep mantle resulting from thermal, chemical, phase and melting effects, Earth-Sci. Rev. , 110, 1– 25. Google Scholar CrossRef Search ADS   Tackley P.J., King S.D., 2003. Testing the tracer ratio method for modeling active compositional fields in mantle convection simulations, Geochem. Geophys. Geosys. , 4, doi:10.1029/2001GC000214. Tanaka S., 2010. Constraints on the core‐mantle boundary topography from P4KP‐PcP differential travel times, J. geophys. Res. , 115, B04310, doi:10.1029/2009JB006563. Tarduno J.A., Watkeys M.K., Huffman T.N., Cottrell R.D., Blackman E.G., Wendt A., Scribner C.A., Wagner C.L., 2015. Antiquity of the South Atlantic Anomaly and evidence for top-down control on the geodynamo, Nat. Commun ., 6, doi:10.1038/ncomms8865. Tateno S., Hirose K., Sata N., Ohishi Y., 2009. Determination of post-perovskite phase transition boundary up to 4400 K and implications for thermal structure in D″ layer, Earth planet. Sci. Lett. , 277, 130– 136. Google Scholar CrossRef Search ADS   Trampert J., Deschamps F., Resovsky J.S., Yuen D.A., 2004. Probabilistic tomography maps significant chemical heterogeneities in the lower mantle, Science , 306, 853– 856. Google Scholar CrossRef Search ADS PubMed  Woodhouse J.H., Girnius T.P., 1982. Surface waves and free oscillations in a regionalized Earth model, Geophys. J. Int. , 68, 653– 673. Google Scholar CrossRef Search ADS   Wu X., Wahr J.M., 1997. Effects of non-hydrostatic core-mantle boundary topography and core dynamics on Earth rotation, Geophys. J. Int. , 128, 18– 42. Google Scholar CrossRef Search ADS   Yamazaki D., Karato S.-I., 2001. Some mineral physics constraints on the rheology and geothermal structure of Earth's lower mantle, Am. Mineral. , 86, 385– 391. Google Scholar CrossRef Search ADS   Yang H.-Y., Tromp J., 2015. Synthetic free-oscillation spectra: an appraisal of various mode-coupling methods, Geophys. J. Int. , 203, 1179– 1192. Google Scholar CrossRef Search ADS   Yoshida M., 2008. Core-mantle boundary topography estimated from numerical simulations of instantaneous mantle flow, Geochem. Geophys. Geosys. , 9, doi:10.1029/2008GC002008. Zhang S., Christensen U., 1993. Some effects of lateral viscosity varuiations on geoid and surface velocities induced by density anomalies in the mantle, Geophys. J. Int. , 114, 531– 547. Google Scholar CrossRef Search ADS   © The Authors 2017. Published by Oxford University Press on behalf of The Royal Astronomical Society.

### Journal

Geophysical Journal InternationalOxford University Press

Published: Jan 1, 2018

## You’re reading a free preview. Subscribe to read the entire article.

### DeepDyve is your personal research library

It’s your single place to instantly
that matters to you.

over 18 million articles from more than
15,000 peer-reviewed journals.

All for just $49/month ### Explore the DeepDyve Library ### Search Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly ### Organize Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place. ### Access Get unlimited, online access to over 18 million full-text articles from more than 15,000 scientific journals. ### Your journals are on DeepDyve Read from thousands of the leading scholarly journals from SpringerNature, Elsevier, Wiley-Blackwell, Oxford University Press and more. All the latest content is available, no embargo periods. DeepDyve ### Freelancer DeepDyve ### Pro Price FREE$49/month
\$360/year

Save searches from
PubMed

Create lists to

Export lists, citations