Influence of the Ringwoodite-perovskite Transition on Mantle Convection in Spherical Geometry as a Function of Clapeyron Slope and Rayleigh Number

We investigate the influence on mantle convec-tion of the negative Clapeyron slope ringwoodite to per-ovskite and ferro-periclase mantle phase transition, which is correlated with the seismic discontinuity at 660 km depth. In particular, we focus on understanding the influence of the magnitude of the Clapeyron slope (as measured by the Phase Buoyancy parameter, P) and the vigour of convection (as measured by the Rayleigh number, Ra) on mantle con-vection. We have undertaken 76 simulations of isoviscous mantle convection in spherical geometry, varying Ra and P. Three domains of behaviour were found: layered convection for high Ra and more negative P , whole mantle convection for low Ra and less negative P , and transitional behaviour in an intervening domain. The boundary between the lay-ered and transitional domain was fit by a curve P = αRa β where α = −1.05, and β = −0.1, and the fit for the boundary between the transitional and whole mantle convection domain was α = −4.8, and β = −0.25. These two curves converge at Ra ≈2.5×10 4 (well below Earth mantle vigour) and P ≈ −0.38. Extrapolating to high Ra, which is likely earlier in Earth history, this work suggests a large transitional domain. It is therefore likely that convection in the Archean would have been influenced by this phase change, with Earth being at least in the transitional domain, if not the layered domain.


Introduction
Mantle convection has had a dominant control on Earth's surface evolution.It has been known for many years that mineral phase changes with negative Clapeyron slope (γ ), where γ = dP/dT, are capable of layering mantle convection (Olson and Yuen, 1982;Christensen, 1995).Layering occurs when the lateral temperature variations of convection produce laterally varying vertical deflections of the boundary away from its equilibrium position.The nature and magnitude of the deflections depend on the value of γ .If there is an appropriate restoring density contrast between the phases, the deflection is accompanied by buoyancy forces that act against the convective thermal buoyancy.If the Clapeyron slope and density difference are sufficient, the buoyancy forces resulting from the phase boundary deflections can overcome the local convective thermal buoyancy, resulting in layered convection (Fig. 1).
The mineral phase change in the olivine system from ringwoodite to ferro-periclase and Mg-perovskite is known to have a negative Clapeyron slope (i.e. is an endothermic reaction) and a restoring density increase.This reaction is widely believed to correspond to the seismic discontinuity at 660 km depth, separating the upper and lower mantle (Bernal, 1936;Ringwood, 1969;Shim et al., 2001).Mineralogical work produces estimates of the Clapeyron slope between −0.5 MPa K −1 and −3.5 MPa K −1 , with recent values for a dry mantle heading towards the less negative value (Katsura et al., 2003;Fei et al., 2004;Litasov et al., 2005;Hirose, 2002;Ito and Takahashi, 1989;Irifune et al., 1998).Recent measurements for a "wet" mantle have the more negative values (Ohtani and Litasov, 2006) The fact that the material below the nominal depth of the phase change can still be the lighter Ringwoodite phase is shown by the letter A on the phase diagram to the left.As a result, there is a downward deflection of the boundary, which leads to the lighter ringwoodite phase producing a buoyancy force and potentially layering.Part (B) on the right hand side illustrates that the same mechanism also can work in reverse, potentially leading to layering with hot upwellings.based on estimates of the deflection of the boundary and independent estimates of the thermal perturbations would suggest a Clapeyron slope around −2.0 to −3.5 MPa K −1 (Lebedev et al., 2002;Fukao et al., 2009).We note that Liu (1994) has pointed out that the slope of the Clapeyron curve for this reaction might not be negative at higher temperatures; in that case in these regions the convection would not be layered.Hirose (2002) demonstrated experimentally that this transition might occur at around 2070 K.This temperaturedependence to the phase transition leads to the possibility that upwellings and downwellings might be affected differently.We note that present-day estimates of mantle temperature at 660 km depth might be around 1880 K (Katsura et al., 2004); therefore, the negative Clapeyron slope could have played a dynamic role, at least for average and cold regions of the mantle, from early in Earth history.
There is also an uncertainty regarding the density contrast across this phase change, probably lying between 7.0 % and 9.3 % (Dziewonski and Anderson, 1981;Weidner and Wang, 1998;Billen, 2008).The mean absolute density, in contrast, is fairly well constrained from global radial seismological models at just over 4.15 Kg m −3 (Dziewonski and Anderson, 1981).The coefficient of thermal expansion at around 660 km depth is probably slightly less than 2 × 10 −5 K −1 for base of upper mantle and slightly more than 2 × 10 −5 K −1 for uppermost lower mantle (Steinberger and Calderwood, 2006).
Considering the factors outlined above, we can reasonably expect this layering process to affect Earth's mantle dynamics.Early simulations suggest that under some conditions this layering breaks down in an episodic manner such that the evolution can be very time-dependent (Davies, 1995).This layering effect has been shown to be sensitive to the vigour of convection, with a system more likely to layer at higher vigour (Christensen and Yuen, 1985).The motivation of this work is to understand the role that such a phase change might have had on Earth's thermal and geological evolution by investigating how the vigour of convection controls the behaviour in spherical models of mantle convection.
Most work to date has focussed on 2-D, 2-D axisymmetric, and 3-D Cartesian geometries (Christensen and Yuen, 1985;Machetel and Weber, 1991;Liu et al., 1991;Peltier and Solheim, 1992;Zhao et al., 1992;Weinstein, 1993;Steinbach et al., 1993;Solheim and Peltier, 1994;Machetel et al., 1995;Peltier, 1996;Marquart et al., 2001).While there has been notable work including this phase change in spherical geometry (Tackley et al., 1993(Tackley et al., , 1994;;Machetel et al., 1995;Bunge et al., 1997), there has only been limited work in spherical geometry to characterise the influence of the value of the slope of the phase change and the vigour of convection on the behaviour.Elements of our work in spherical geometry have been presented previously (Wolstencroft and Davies, 2008a,b).Yanagisawa et al. (2010) undertook a similar study simultaneously where they recently concluded that 3 domains of behaviour can be identified.Our methodology is similar: we characterise the behaviour of a large number of simulations at varying vigour of convection and strength of phase transition to define a regime diagram.We constrain the boundaries between the various behaviours and parameterise them.This allows extrapolation of the boundaries to very high vigour, which will be relevant to early Earth history but which are currently beyond computational means to simulate.After characterising the behaviour rigorously, we conclude by speculating on its possible consequences for Earth evolution.

Methods
Appropriate assumptions were made to focus on the objective of characterising how the mantle behaves as a function of both the vigour of convection and the magnitude of the negative Clapeyron slope (γ ) of the phase change reaction (i.e. the slope dP/dT of the phase change boundary in Pressure-Temperature space).For example, in our numerical simulations we assumed the mantle is a viscous, incompressible, isochemical fluid, with constant material properties with infinite Prandtl number subject to the Boussinesq approximation (Ricard, 2007).The resulting equations modelled are: the equation for conservation of mass the equation for the conservation of linear momentum and the equation for the conservation of energy where u: velocity, η: dynamic viscosity, p: hydrodynamic pressure, ρ: density, ρ: lateral variation in density away from the reference state, g: acceleration due to gravity, r: unit vector in radial direction, C p : specific heat at constant pressure, T : temperature, k: thermal conductivity, and H : radiogenic heat production per unit mass.
Non-dimensionalisation of these equations shows that the behaviour is controlled by only 3 non-dimensional parameters: the basal heated Rayleigh number Ra, the internal heating Rayleigh number Ra H , and the Phase Buoyancy parameter P (e.g.Bunge et al., 1997).
The basal-heated Rayleigh number Ra, is given by where α: coefficient of thermal expansion, T : superadiabatic temperature drop across the shell, D: mantle thickness, and κ: thermal diffusivity = k ρC p .The internal-heated Rayleigh number Ra H is given by Since only η and γ were varied, the ratio of Ra/Ra H = T k/ρH D 2 = 0.054 was always constant in this work.We therefore only need one of the two Rayleigh numbers to describe our experiments; we use Ra.This use of a fixed Ra/Ra H ratio is different to the approach adopted by Yanagisawa et al. (2010), where the Ra H was reduced as Ra increased.The impact of the different treatment of Ra H is that, in this study, the proportional contribution of the basal heat flux to total heat flux increases with Ra.In the work of Yanagisawa et al. (2010), this proportion remains approximately constant with Ra.
The Phase Buoyancy Parameter P is given by where γ : Clapeyron slope of the phase change, δρ: density change across the phase change, and ρ: mean density of the two phases.
Using these equations we investigated the form of layered mantle convection for varying Ra and P .The simulations were undertaken using a benchmarked version of TERRA (Baumgardner, 1985;Bunge et al., 1997;Davies and Davies, 2009), which is able to model very high Ra in spherical geometry (Wolstencroft et al., 2009).The details of how TERRA solves the dimensional form of the equations with shows the constant thermal gradient of a purely conductive regime.
The dashed line in Part (B) shows the whole mantle convection thermal structure, illustrating the large temperature gradients at its boundaries.The solid line represents the absolute radial mass flux, which is zero at the two boundaries.For whole mantle convecting it peaks in mid-mantle.(C) illustrates a layered system.In this case, the solid line shows that the absolute radial mass flux will be a minimum at the phase change, which also leads to an additional thermal boundary layer, illustrated by a step in the dashed line, across the phase change.(D) shows a cross-section and a radial surface through the thermal anomaly field.Blue is colder than the radial average and red hotter.The visualisation demonstrates a transitional case.There is a large degree of layering, with the 660 km phase change clearly visible by temperature contrasts across the interface, but there is also some passage of material across.
pressure, temperature and velocity as the free variables, are presented in Baumgardner (1985); Bunge and Baumgardner (1995), and Yang and Baumgardner (2000).TERRA uses dimensional variables and the basic parameters of the simulations are listed in Table 1.Convection was isochemical, with free-slip and isothermal boundary conditions and a component of internal heating.Spatial resolution was carefully selected to resolve the features of each simulation.This resulted in the highest vigour run requiring a spatial resolution of around 15 km near the surface and around 22 km midmantle.The very lowest vigour runs only required a resolution of around 88 km.Using a resolution appropriate to the simulation helps make this large parameter survey possible.A total of 76 simulations were undertaken, with γ ranging from −2 to −30 MPa K −1 (the equivalent phase buoyancy parameter P for this range of γ ranges from −0.0554 to −0.831) and Ra from 5.19 × 10 3 to 8.49 × 10 7 (Table 2).We note that for a sub-set of the simulations, a positive Clapeyron slope phase change was included at 410 km depth equal to 1.5 MPa K −1 .Since the "660" Clapeyron slope has very large negative values, whether the small slope 410 km phase change is included or set to zero usually makes little difference (see Table 2).In fact at low Ra, where this might make a difference, many simulations were undertaken with and without; no difference was found.The phase change is implemented at a depth of 660 km depth using the sheet mass anomaly method of Tackley et al. (1993) implemented in TERRA by Bunge et al. (1997).This method applies an appropriate body force at the discontinuity depth, dependent on the local temperature and phase change parameters.Tackley et al. (1993) shows that this is an excellent method especially when the radial resolution cannot resolve the likely boundary deflections, as in our simulations.The sheet mass anomaly method effectively also assumes an infinitesimally thin phase loop.The sharp discontinuity observed seismically at short periods suggests that this phase change does have a thin loop (Shearer and Masters, 1992).Mineralogy experiments also suggest a thin loop (Ito and Takahashi, 1989;Wood, 1990).Earlier simulations with numerical methods including the width of the loop suggest that the large scale dynamics are largely insensitive to its width within reasonable parameter values (Christensen and Yuen, 1985;Tackley, 1995).In order to be consistent with the assumption of incompressiblity, effects related to latent heat have been disregarded.While latent heat is known to have an influence at low Ra (Schubert and Turcotte, 1971), with increasing Ra, boundary deflection is known to dominate (Christensen, 1995).The work of Ita and King (1994) and Christensen and Yuen (1985) showing that Boussinesq, extended Boussinesq and compressible simulations give similar results reinforces that latent heat effects are of minor importance at high Ra.Therefore, since we are interested in how the behaviours might extrapolate to high Ra (where influence of latent heat would be expected to be irrelevant), this simplification is sensible.We assume a density contrast across the phase change δρ of δρ/ρ ≈ 9 %.The constancy of δρ implies a simple linear relationship between our Clapeyron slope (γ ) and the Phase Buoyancy Parameter (P ).We should remind the reader that the individual values of parameters, including δρ, γ , are ultimately not important.Only the resulting values of the controlling non-dimensional parameters, Ra, Ra H and P , are important.
Each simulation was initiated by a radially uniform, laterally small scale, random structure.The results were independent of this initial condition since they were run for an extended length of time, until a thermal quasi-steady state was reached.This was achieved as follows.The surface and basal heat flux over time were plotted for each case.As the cases contain no evolution in heat sources, these curves show a "flattening" when balance is achieved between heat input and heat output.Cases were run until this balance was achieved.The transitional class of cases do not generally show flattening but tend to oscillate about a particular SHF value instead.Cases were classified as layered, whole mantle or transitional using the following process: 1.The total absolute radial mass flux across the mantle and the average temperature profile were plotted as shown in Fig. 2.
2. If the mass flux at 660 km depth was approx.90% or greater of the maximum mass flux of the model, the case was classified whole mantle convection.
3. If the mass flux at 660 km depth was approx.10% or lower than the maximum mass flux, the case was classified as layered.
4. If the case fell between these classifications and demonstrated periodic variability in surface heat flux, it was classified as transitional.The temperature profile was used to verify this state since a layered system develops a characteristic internal thermal boundary layer due to the lack of advective heat transport across the phase change (Fig. 2).Whole mantle convection shows no such internal thermal boundary layer.
5. Additionally, the cases were visualised in 3-D using isosurfaces, cross sections and specific depth surfaces (e.g.Fig. 2).These plots allow judgments to be made as to what proportion of up or down welling features penetrate the phase boundary and what the nature of the coupling is across the boundary, e.g.thermal or viscous.This final step was important for correctly classifying boundary cases where summary data values were inconclusive.

Results
The results show that at low Ra and high P (i.e.low absolute value) whole mantle convection is preferred, while at high Ra and more negative P , layered convection is preferred (Table 2, Fig. 3).An intermediate domain of transitional behaviour is found between the two end-member behaviours.
In previous work, this transitional region has been termed partially layered; we avoid this term as it does not accurately represent the wide range of behaviours displayed.for the boundary between the layered regime and the transitional regime, and for the boundary between the transitional regime and the whole mantle convection regime.We note that the line fits to the data are not perfect.There is a suggestion that the form of the curves, i.e. simple power law, might not be the true form of the relationship.Without guidance of an alternative relationship, we argue it is best to keep with a simple relationship, which in this case has a long history of usage in the field (Christensen and Yuen, 1985).We now go on to compare these results with earlier work and discuss their possible implications.

Domains of convective mode
As described in the Introduction, there has been a long history of investigating the effect of Ra and γ on the behaviour of mantle convection with mineral phase changes.As mentioned above, our work shows three domains of behaviour: a layered convection domain at high Ra and more negative P , a whole mantle convection domain at low Ra and less negative P , and a transitional domain at intermediate values of Ra and P .
To repeat, this work focuses just on the buoyancy effect resulting from the deflection of the phase change boundary.By undertaking incompressible convection simulations, which to be self-consistent ignore the effects of the latent heat of the reaction, we are removing its influence.
We do observe episodic behaviour in the transitional domain.It is beyond the scope of this project to undertake sufficient simulations to characterise accurately where this behaviour occurs on the regime diagram, but we do note that it does not extend across the whole transitional domain and seems to occur for parameter values closer to the layered domain.It would be interesting to constrain the boundary of this behaviour.
The range of Ra and particularly P used in this study are well beyond the likely values for Earth.Our reasons for doing this were to be able to better constrain the boundaries of the various convective domains.However, these values may have application when considering other materials (such as water ice) in other bodies such as ice planets.

Fits to domain boundaries
The earliest reference, to our knowledge, for a curve fit to the boundary between layered and whole mantle convection was from Christensen and Yuen (1985) 2-D Cartesian work, where the form P = αRa β was assumed and the values found to be: α = −4.4 and β = −0.2.The work presented here is similar to that of Yanagisawa et al. (2010).They also used the code TERRA and very similar input parameters to ours.The key differences between the two sets of simulations were www.solid-earth.net/2/315/2011/Solid Earth, 2, 315-326, 2011   that this study used a fixed Ra/Ra H and a 150 K greater temperature change across the mantle.In some cases there was also a difference in the Clapeyron slope assumed for the phase change at 410 km depth: 0 and 1.5 MPa K −1 in our work against 0 and 2.5 MPa K −1 in the work of Yanagisawa et al. (2010).Yanagisawa et al. (2010), like us, also define a transitional regime and their fits to the two boundaries have α = −4.7 and β = −0.2 for the boundary between layered and transitional, and α = −12 and β = −0.33 for the boundary between whole mantle convection and transitional.
Figure 4 shows these previous domain boundaries with our curves.We note that these curves agree well at medium-high Ra where most of the previous data points lie, but that there are differences at lower Ra and in the width of the transitional domain.
Therefore, while we would not expect our and Yanagisawa et al. ( 2010) results to be identical (interpretation of the boundary cases might also differ), they should be sufficiently similar that it is worthwhile to consider the curves together (Fig. 4).We note that our curves also satisfy most of their points.This gives further confidence in our curves, especially as we extrapolate to high Ra, where the Yanagisawa et al. ( 2010) study has many data-points.
We note that the two curves in our study converge at low Rayleigh number -we estimate this point to be around Ra ≈ 2.5 × 10 4 and P ≈ −0.38.The curves of Yanagisawa et al. (2010) do not converge at low Ra; in fact, they have only a limited number of simulations in this range.If it is correct that the two curves should converge at low Ra, then it tightens considerably the fit to the points and therefore the extrapolation of this relationship to high Ra.At very high Ra, the two curves will approach the Ra axis but never cross it.This is as one would expect since a mineral reaction with a non-negative Clapeyron slope cannot layer the flow with the resulting boundary deflection.

Limitations of modelling
It is interesting to speculate what this research implies for mantle convective behaviour earlier in Earth history.Before doing so, we would like to emphasise that the simulations were intentionally simple to allow understanding and complete investigation of the parameter space.As a result, there are limitations when applying these outcomes to Earth.The models are without depth or temperature-dependent viscosity, do not have plates or continents, and thus can not display behaviours such as slab-rollback (Goes et al., 2008;Yanagisawa et al., 2010).These simplifications could be expected to affect the detail of the conclusions here; hopefully, the broad trends in this space will remain valid as more sophisticated models are investigated in the future.More sophisticated discussions should also consider other mineral reactions, both in the olivine and garnet systems, the effect of latent heat (Christensen, 1998) and volume change (Krien and Fleitout, 2010) of the phase changes.We note that the olivine component of the mantle probably makes up around 60 % of the mantle; therefore, when applying this work to the mantle the effective value of P should be reduced by a similar proportion.Since there is some uncertainty regarding the proportion of olivine in the transition zone (see for example Anderson and Bass (1986) who argue for a piclogitic transition zone), our approach of not making this correction in the simulations makes it simpler for others to use this work.

Rayleigh number
Applying a single number like a Rayleigh number is clearly fraught for the actual Earth, since many properties are spatially variable and not constant.Accepting this, to advance the discussion we estimate that present day mantle Ra might be ≈ 10 7 , assuming, e.g.mean viscosity η ≈ 5 × 10 21 Pa s, κ ≈10 −6 m 2 s −1 , α ≈2 × 10 −5 K −1 , super-adiabatic temperature drop of ≈2550 K (Steinberger and Calderwood, 2006).Earlier in Earth history, a hotter mantle is to be expected due to the dissipation of gravitational energy from the formation era and higher radioactivity.Limited observations support this assumption that the mantle was hotter earlier in Earth history (Nisbet et al., 1995;Green, 1975).A hotter mantle would translate to a higher Rayleigh number.To illustrate the potential implications, we will assume that the only changing parameter in the Ra as a function of time is viscosity, with it being lower at higher temperatures.Mantle rheology is a field with large uncertainties but an activation energy of 500 KJ mol −1 (Korenaga and Karato, 2008) would suggest that mantle viscosity might decrease by approximately an order of magnitude for every 100 degrees increase in temperature.
Magmatic products also suggest that the mantle was hotter in the Archean, with estimates varying from 100 K to 300-500 K hotter at 3 Ga depending upon interpretations such as how wet the komatiite source region was and what is representative of average mantle (Lee et al., 2009;Grove and Parman, 2004;Nisbet et al., 1995;Abbott et al., 1994).Assuming a 200 K hotter mantle at ≈ 3 Ga, we might expect Earth to have a viscosity 2 orders of magnitude lower and therefore Ra 2 orders of magnitude greater at Ra ≈ 10 9 , but clearly with significant uncertainty.2010) with the curve from Christensen and Yuen (1985) overlaid on the domain boundary lines of this study.While we note a reasonable agreement between the points and the lines from this study, all the points do not satisfy our lines.This could be the result of the slightly different parameters and the measure of subjectivity in defining the behaviour category of a simulation.Key differences between our curves and those from previous work are the wider transitional region at high Ra and the meeting point of our curves.

Transitional Layered Convection
Whole-Layered

Present-day Earth
Using the estimates above for the Rayleigh number and the Phase Buoyancy parameter, we mark Earth's current position on the domain diagram (Fig. 5) with a large Earth whose size suggests an approximate sense of the uncertainty of the parameters.The approximate Ra and P values would suggest that we are today either in the transitional regime or just in the whole mantle convection regime and very unlikely to be in the layered regime.Seismological evidence seems to strongly support this, with observations of subducting slabs descending from the upper to the lower mantle (Grand et al., 1997;van der Hilst et al., 1997;Creager and Jordan, 1984) and also observations of stagnant slabs, which might reflect some resistance at this boundary (Fukao et al., 1992)

Earth evolution
If we speculate as to where Earth was on this regime diagram 3 Ga ago, we might expect it to be at a Ra approximately 2 orders of magnitude greater and similar P .This would suggest, in a simple interpretation, that Earth would have most likely occupied the deeper regions of the transitional domain in the past.It is possible that in very early Earth history, the mantle operated in a layered convection mode.Therefore, as we look at our results, we would predict that over Earth history the mantle evolved from a layered/transitional regime to a dominantly whole mantle convection regime.In such a scenario it is possibile that Earth might have had episodic mantle convection in its earlier history as previously suggested (Condie, 2001(Condie, , 1998;;Parman, 2007;Pearson et al., 2007;Frimmel, 2008;Ernst and Buchan, 2002).We also note that the evidence presented for episodicity is controversial; zircon peaks might reflect preservation (Hawkesworth et al., 2009), while it has been argued that the isotopic peaks may not be statistically robust (Rudge, 2008).Phase-change induced mantle avalanches could initiate superplumes/superevents (Condie, 1998).Such events have been suggested to affect not just magmatic outputs, but also core-generated magnetic fields (Larson, 1991).Clearly there is the potential for multiple observations to be affected.Future work to better constrain the input parameters, Earth history and undertake more realistic simulations are encouraged by the work to date.

Conclusions
We discover 3 domains of behaviour for a spherical geometry convecting mantle with a negative Clapeyron slope phase change simulating the ringwoodite to ferro-periclase and Mg-Perovskite transition at 660 km depth: a whole mantle convection domain, a layered convection domain, and a transitional domain.The boundaries separating the domains converge at the low near-critical Rayleigh number, while the transitional domain (which includes episodic behaviour) is very broad at realistic Clapeyron slope.By extrapolating power law fits of these well constrained domain boundaries to high Rayleigh number (convective vigour), we suggest that it is likely that the transitional domain and possibly also the layered domain will be of interest during early Earth history and therefore for understanding Earth evolution.This work encourages more realistic simulations to be undertaken in the future as more computational resources become available.

Fig. 1 .
Fig. 1.Illustration of the layering mechanism modelled.The left hand side figure shows a schematic version of the phase diagram of the olivine component -Ringwoodite is stable in the upper left region (labelled with R) while Perovskite and Ferropericlase is stable in the region (labelled with Pv) below the phase boundary line.The right hand side: Part (A) illustrates a cold downwelling impinging on the phase change.The fact that the material below the nominal depth of the phase change can still be the lighter Ringwoodite phase is shown by the letter A on the phase diagram to the left.As a result, there is a downward deflection of the boundary, which leads to the lighter ringwoodite phase producing a buoyancy force and potentially layering.Part (B) on the right hand side illustrates that the same mechanism also can work in reverse, potentially leading to layering with hot upwellings.

Fig. 2 .
Fig. 2. Criteria for defining whether convection is layered.(A)shows the constant thermal gradient of a purely conductive regime.The dashed line in Part (B) shows the whole mantle convection thermal structure, illustrating the large temperature gradients at its boundaries.The solid line represents the absolute radial mass flux, which is zero at the two boundaries.For whole mantle convecting it peaks in mid-mantle.(C) illustrates a layered system.In this case, the solid line shows that the absolute radial mass flux will be a minimum at the phase change, which also leads to an additional thermal boundary layer, illustrated by a step in the dashed line, across the phase change.(D) shows a cross-section and a radial surface through the thermal anomaly field.Blue is colder than the radial average and red hotter.The visualisation demonstrates a transitional case.There is a large degree of layering, with the 660 km phase change clearly visible by temperature contrasts across the interface, but there is also some passage of material across.

MFig. 3 .
Fig. 3. Layering status of modelled cases as a function of convective vigour (Ra) and layering strength of the phase change (Buoyancy Parameter, P ).Key defines the symbols used to indicate the nature of the convection.The solid line is a power-law fit attempting to divide the whole and transitional cases, while the dashed line divides the layered and transitional cases.The vertical dash-dot line represents the whole-layered boundary, which only occurs at very low Ra and very negative Clapeyron slope values.The equivalent Clapeyron slope for our parameters is shown on the right hand side vertical axis.

Fig. 4 .
Fig. 4. Data points and curves from the study of Yanagisawa et al. (2010) with the curve fromChristensen and Yuen (1985) overlaid on the domain boundary lines of this study.While we note a reasonable agreement between the points and the lines from this study, all the points do not satisfy our lines.This could be the result of the slightly different parameters and the measure of subjectivity in defining the behaviour category of a simulation.Key differences between our curves and those from previous work are the wider transitional region at high Ra and the meeting point of our curves.

Fig. 5 .
Fig. 5.The figure of Earth represents an estimate of where the present-day Earth might fit on a Rayleigh number versus Buoyancy Parameter plot.The lines mark the boundaries between the 3 regimes found in this work.The arrow shows an approximate potential mantle evolution route through the parameter space over Earth history.
To better Solid Earth, 2, 315-326, 2011 www.solid-earth.net/2/315/2011/constrain the boundary, a large number of simulations were run at two values of P (−0.221 and −0.332), spanning small ranges of Ra.We have fit curves of the form P = αRa β to these boundaries.The best fit curves defining our boundaries are P = −1.05Ra −0.1 (7)